常微分方程模型及其数值解.ppt

常微分方程数值解,杨瑞琰,0、导言,在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的导数之间的一种平衡关系。而正确地构建这种平衡关系，需要对实际问题的深入浅出的刻画，根据物理的和非物理的原理、定律或定理，作出合理的假设和简化并将它升华成数学问题。,另一个是方程的求解和结果分析。对一些常系数的或特殊函数形式的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法数值解法。本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。,欧拉方法和龙格库塔方法,一阶常微分方程初值问题的一般形式为,y=(x,y),axb,(4),y(a)=,其中(x,y)是已知函数,为给定的初值.,如果函数(x,y)在区域axb,-y0为Lipschitz常数,则初值问题(4)有唯一解.,所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.,a=x0x1x2xnxN=b,其中剖分节点xn=a+nh,n=0,1,N,h称为剖分步长.数值解法就是求精确解y(x)在剖分节点xn上的近似值yny(xn),n=1,2,n.,假设初值问题(4)的解y=y(x)唯一存在且足够光滑.对求解区域a,b做剖分,我们采用数值积分方法来建立差分公式.,1构造数值解法的基本思想,在区间xn,xn+1上对方程(4)做积分,则有,对右边的积分应用左矩形公式，则有,梯形公式,o,x,y,a,b,左矩形公式,y=(x),右矩形公式,中矩形公式,对右边的积分应用左矩形公式，则有,因此，建立节点处近似值yn满足的差分公式,称之为Euler公式.,称为梯形公式.,若对(6)式右边的积分应用梯形求积公式，则可导出差分公式,利用Euler方法求初值问题,解此时的Euler公式为,称为Euler中点公式或称双步Euler公式.,若在区间xn-1,xn+1上对方程(4)做积分,则有,对右边的积分应用中矩形求积公式，则得差分公式,例3,的数值解.此问题的精确解是y(x)=x/(1+x2).,分别取步长h=0.2,0.1,0.05,计算结果如下,functionoutx,outy=MyEuler(fun,x0,xt,y0,PointNum)%MyEuler用前向差分的欧拉方法解微分方程%fun表示f（x，y）%x0,xt表示自变量的初值和终值%y0表示函数在x0处的值,其可以为向量形式%PointNum表示自变量在x0,xt上取的点数ifnargin5|PointNum=0%如果函数仅输入4个参数值，则PointNum默认值为100PointNum=100;endifnargin4%y0默认值为0y0=0;end,MyEuler.m,h=(xt-x0)/PointNum;%计算步长hx=x0+0:PointNum*h;%自变量数组y(1,:)=y0(:);%将输入存为行向量，输入为列向量形式fork=1:PointNumf=feval(fun,x(k),y(k,:);%计算f(x,y)在每个迭代点的值f=f(:);y(k+1,:)=y(k,:)+h*f;%对于所取的点x迭代计算y值endouty=y;outx=x;%plot(x,y)%画出方程解的函数图,functionf=myfun011(x,y)f=1/(1+x*x)-2*y*y;,myfun011.m,functionex01()x1,y1=MyEuler(myfun011,0,2,0,10);%欧拉法所得的解h1=2/10;%计算步长x11,y11=MyEuler(myfun011,0,2,0,20);%欧拉法所得的解h2=2/20;%计算步长x12,y12=MyEuler(myfun011,0,2,0,100);%欧拉法所得的解h2=2/100;%计算步长plot(x1,y1,+b,x11,y11,og,x12,y12,-r)legend(h=0.2的欧拉法解,h=0.1的欧拉解,h=0.05的欧拉解),ex01.m,Euler中点公式则不然,计算yn+1时需用到前两步的值yn,yn-1,称其为两步方法,两步以上的方法统称为多步法.,在Euler公式和梯形公式中,为求得yn+1,只需用到前一步的值yn,这种差分方法称为单步法,这是一种自开始方法.,隐式公式中,每次计算yn+1都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.,在Euler公式和Euler中点公式中,需要计算的yn+1已被显式表示出来,称这类差分公式为显式公式,而梯形公式中,需要计算的yn+1隐含在等式两侧,称其为隐式公式.,从数值积分的角度来看,梯形公式,计算数值解的精度要比Euler公式好,但它属于隐式公式,不便于计算.,实际上,常将Euler公式与梯形公式结合使用:,2改进的Euler方法,由迭代法收敛的角度看，当,(是给定的精度要求)时,取,就可以保证迭代公式收敛,而当h很小时,收敛是很快的.,而且,只要,通常采用只迭代一次的算法:,称之为改进的Euler方法.这是一种单步显式方法.,改进的Euler方法也可以写成,y=y-2x/y,0x1,的数值解,取步长h=0.1.精确解为y(x)=(1+2x)1/2.,例4求初值问题,y(0)=1,解(1)利用Euler方法,计算结果如下:,(2)利用改进Euler方法,functionXout,Yout=MyEulerPro(fun,x0,xt,y0,PointNumber)%MyEulerPro用改进的欧拉法解微分方程ifnargin5|PointNumber=0%如果函数仅输入4个参数值，则PointNumer默认值为100PointNumer=100;endifnargin4%y0默认值为0y0=0;end,MyEulerPro.m,h=(xt-x0)/PointNumber;%计算所取的两离散点之间的距离x=x0+0:PointNumber*h;%表示出离散的自变量xy(1,:)=y0(:);fori=1:PointNumber%迭代计算过程f1=h*feval(fun,x(i),y(i,:);f1=f1(:);f2=h*feval(fun,x(i+1),y(i,:)+f1);f2=f2(:);y(i+1,:)=y(i,:)+1/2*(f1+f2);endXout=x;Yout=y;,functionf=myfun012(x,y)f=y-2*x/y;,myfun012.m,functionex02()%ex02比较改进欧拉法，简单欧拉方法以及微分方程符号解x3,y3=MyEulerPro(myfun012,0,1,1,10);x,y1=MyEuler(myfun012,0,1,1,10);%欧拉法所得的解fork=1:11t(k)=sqrt(1+2*x(k);endplot(x,y1,-b,x3,y3,og,x,t,:*r)legend(简单欧拉法解,改进欧拉解,解析解),ex02.m,在节点xn+1的误差y(xn+1)-yn+1,不仅与yn+1这一步计算有关,而且与前n步计算值yn,yn-1,y1都有关.,为了简化误差的分析,着重研究进行一步计算时产生的误差.即假设yn=y(xn),求误差y(xn+1)-yn+1,这时的误差称为局部截断误差,它可以反映出差分公式的精度.,3差分公式的误差分析,如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.,研究差分公式阶的重要手段是Taylor展开式,一元函数和二元函数的Taylor展开式为:,另外,在yn=y(xn)的条件下,考虑到y(x)=(x,y(x),则有,y(xn)=(xn,y(xn)=(xn,yn)=n,y(xn)=(xn,y(xn)=x(xn,yn)+y(xn,yn)(xn,yn),yn+1=yn+h(xn,yn),对Euler方法，有,=yn+(xn,yn)h+O(h2),从而有:y(xn+1)-yn+1=O(h2),所以Euler方法是一阶方法.,再看改进Euler方法,因为,可得,所以,改进的Euler方法是二阶方法.,而,从而有:y(xn+1)-yn+1=O(h3),4Taylor展开方法,设y(x)是初值问题（4)的精确解,利用Taylor展开式可得,称之为p阶Taylor展开方法.,因此，可建立节点处近似值yn满足的差分公式,其中,所以,此差分公式是p阶方法.,由于Taylor展开方法涉及很多复合函数(x,y(x)的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而,Taylor展开方法给出了一种构造单步显式高阶方法的途径.,Euler方法可写为,可见,公式的局部截断误差为:y(xn+1)-yn+1=O(hp+1).,5Runge-Kutta方法,构造差分公式,改进的Euler方法可写为,其中i,i,ij为待定参数.,若此公式的局部截断误差为,由于,yn+1=yn+h1n+h2(n+hxn+hnyn)+O(h3),O(h3),称此公式为p阶Runge-kutta方法,简称p阶R-K方法.,对于p=2的情形,应有,=yn+h(1+2)n+h22(xn+nyn)+O(h3),所以,只要令,1+2=1,2=1/2,2=1/2(8),一般地,参数由(8)确定的一族差分公式(7)统称为二阶R-K方法.,称之为中点公式,或可写为,若取=1,则得1=2=1/2,=1,此时公式(7)就是改进的Euler公式;,若取1=0,则得2=1,=1/2,公式(7)为,高阶R-K公式可类似推导.,下面列出常用的三阶、四阶R-K公式.,四阶标准R-K公式,三阶R-K公式,解四阶标准R-K公式为,例3用四阶标准R-K方法求初值问题,y=y-2x/y,0x1,y(0)=1,的数值解,取步长h=0.2.,计算结果如下:,functionx,y=MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta方法解微分方程形为y(t)=f(x,y(x)%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式%x范围为x0,xt,初值为y0,PointNum为离散点数,varargin为可选输入项可传适当参数给函数f(x,y)ifnargin4|PointNum=0PointNum=100;endifnargin3y0=0;end,MyRunge_Kutta.m,y(1,:)=y0(:);%初值存为行向量形式h=(xt-x0)/(PointNum-1);%计算步长x=x0+0:PointNum*h;%得x向量值fork=1:PointNum%迭代计算f1=h*feval(fun,x(k),y(k,:),varargin:);f1=f1(:);%得公式中k1f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin:);f2=f2(:);%得公式中k2f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin:);f3=f3(:);%得公式中k3f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin:);f4=f4(:);%得公式中k4y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;%得y(n+1)end,functionex03()%用三种不同方法解微分方程clear,clfx0=0;xt=1;Num=100;h=(xt-x0)/(Num-1);x=x0+0:Num*h;a=1;yt=sqrt(1+2*x);%真值解y0=1;%设定函数初值PointNum=4;%设定取点数x1,ye=MyEuler(myfun012,x0,xt,y0,PointNum);x2,yh=MyEulerPro(myfun012,x0,xt,y0,PointNum);x3,yr=MyRunge_Kutta(myfun012,x0,xt,y0,PointNum);plot(x,yt,k,x1,ye,b:,x2,yh,g:,x3,yr,r:)legend(真值,简单欧拉法解,改进欧拉法解,Runge-Kutta法解)holdonplot(x1,ye,bo,x2,yh,b+,x3,yr,r*),ex03.m,也可以构造隐式R-K方法,其一般形式为,称之为p级隐式R-K方法,同显式R-K方法一样确定参数.如,是二级二阶隐式R-K方法,也就是梯形公式.但是p级隐式R-K方法的阶可以大于p,例如,一级隐式中点公式为,或写为,它是二阶方法.,RK方法的MATLAB实现,对于微分方程（组）的初值问题龙格库塔方法可用如下MATLAB命令实现其计算：t，x=ode23(f，ts，x0，options)t，x=ode45(f，ts，x0，options)其中ode23用的是3级2阶龙格库塔公式，ode45用的是以Runge-Kutta-Fehberg命名的5级4阶公式。,命令的输入f是待解方程写成的函数m文件：functiondx=f(t,x)dx=f1;f2;fn;若ts=t0,t1,t2,tf，则输出在指定时刻t0,t1,t2,tf的函数值；等分点时用ts=t0:k:tf，输出在t0,tf内等分点处的函数值。x0为函数初值（n维向量）。options可用于设定误差限（options缺省时设定相对误差10-3，绝对误差10-6），命令为：options=odeset(reltol,rt,abstol,at)其中rt，at分别为设定的相对误差和绝对误差。命令的输出t为指定的ts，x为相应的函数值（n维向量）。注意，计算步长是根据误差限自动调整的，并不是输入中指定的ts的分点。,%用ode23得到微分方程解fun=inline(-3*y2+2*x.2+3*x,x,y);%用inline构造函数f(x,y)x,y=ode23(fun,0,1,1);%可得到x,y输出向量值ode23(fun,0,1,1)%可得到输出得函数图fun=inline(-3*y2+2*x.2+3*x,x,y);%用inline构造函数f(x,y)ode45(fun,0,1,1)%可得到输出得函数图x,y=ode45(fun,0,1,1);,functiondx=myfun02(x,y)dx=zeros(2,1);dx(1)=y(1)*(1-0.1*y(2);dx(2)=y(2)*(-0.5+0.02*y(1);clear,clcx,y=ode45(myfun02,015,252);plot(x,y(:,1),-,x,y(:,2),*)%画出y(1),y(2)的函数图legend(y1函数图,y2函数图),用自编函数MyEulerPro（）MyRun