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

《常微分方程模型及其数值解.ppt》由会员分享，可在线阅读，更多相关《常微分方程模型及其数值解.ppt（53页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程数值解常微分方程数值解杨瑞琰0、导言 在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的导数之间的一种平衡关系。而正确地构建这种平衡关系，需要对实际问题的深入浅出的刻画，根据物理的和非物理的原理、定律或定理，作出合理的假设和简化并将它升华成数学问题。另一个是方程的求解和

2、结果分析。对一些常系数的或特殊函数形式的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法数值解法。本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。欧拉方法和龙格欧拉方法和龙格库塔方法库塔方法一阶常微分方程初值问题的一般形式为 y=(x,y),axb(4)y(a)=其中(x,y)是已知函数,为给定的初值.如果函数(x,y)在区域axb,-y0为LipschitzLipschitz常

3、数常数,则初值问题(4)有唯一解.所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.a=x0 x1x2xnxN=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 1 构造数值解法的基本思想构造数值解法的基本思想 在区间xn,xn+1上对方程(4)做积分,则有对右边的积分应用左矩形公式，则有梯形公式oxyab左矩形公式y=(x)右矩形公式中矩形公式对

4、右边的积分应用左矩形公式，则有因此，建立节点处近似值yn满足的差分公式称之为EulerEuler公式公式.称为梯形公式梯形公式.若对(6)式右边的积分应用梯形求积公式，则可导出差分公式 利用Euler方法求初值问题 解解 此时的Euler公式为称为EulerEuler中点公式中点公式或称双步双步EulerEuler公式公式.若在区间xn-1,xn+1上对方程(4)做积分,则有对右边的积分应用中矩形求积公式，则得差分公式例例3的数值解.此问题的精确解是y(x)=x/(1+x2).分别取步长h=0.2,0.1,0.05,计算结果如下function outx,outy=MyEuler(fun,x0

5、,xt,y0,PointNum)%MyEuler 用前向差分的欧拉方法解微分方程%fun 表示f（x，y）%x0,xt表示自变量的初值和终值%y0表示函数在x0处的值,其可以为向量形式%PointNum表示自变量在x0,xt上取的点数if nargin5|PointNum=0%如果函数仅输入4个参数值，则PointNum默认值为100 PointNum=100;endif nargin4%y0默认值为0 y0=0;endh=(xt-x0)/PointNum;%计算步长hx=x0+0:PointNum*h;%自变量数组y(1,:)=y0(:);%将输入存为行向量，输入为列向量形式for k=1:

6、PointNum f=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)%画出方程解的函数图function f=myfun011(x,y)f=1/(1+x*x)-2*y*y;myfun011.mfunction ex01()x1,y1=MyEuler(myfun011,0,2,0,10);%欧拉法所得的解h1=2/10;%计算步长x11,y11=MyEuler(myfun011,0,2,0,20);%欧拉法所得的解h2=2/20

7、;%计算步长x12,y12=MyEuler(myfun011,0,2,0,100);%欧拉法所得的解h2=2/100;%计算步长plot(x1,y1,+b,x11,y11,og,x12,y12,-r)的欧拉法解的欧拉解的欧拉解)hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.0

8、00.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.

9、00481-0.00227Euler中点公式则不然,计算yn+1时需用到前两步的值yn,yn-1,称其为两步方法两步方法,两步以上的方法统称为多步法多步法.在Euler公式和梯形公式中,为求得yn+1,只需用到前一步的值yn,这种差分方法称为单步法单步法,这是一种自开始方法.隐式公式中,每次计算yn+1都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.在Euler公式和Euler中点公式中,需要计算的yn+1已被显式表示出来,称这类差分公式为显式公式显式公式,而梯形公式中,需要计算的yn+1隐含在等式两侧,称其为隐式公式隐式公式.从数值积分的角度来看,梯形公式计算数值解的精度要比E

10、uler公式好,但它属于隐式公式,不便于计算.实际上,常将Euler公式与梯形公式结合使用:2 改进的改进的Euler方法方法 由迭代法收敛的角度看，当 (是给定的精度要求)时,取 就可以保证迭代公式收敛,而当h很小时,收敛是很快的.而且,只要 通常采用只迭代一次的算法:称之为改进的改进的Euler方法方法.这是一种单步显式方法.改进的Euler方法也可以写成 y=y-2x/y ,0 x1的数值解,取步长h=0.1.精确解为y(x)=(1+2x)1/2.例例4 求初值问题 y(0)=1 解解 (1)利用Euler方法计算结果如下:(2)利用改进Euler方法function Xout,Yout

11、=MyEulerPro(fun,x0,xt,y0,PointNumber)%MyEulerPro 用改进的欧拉法解微分方程if nargin5|PointNumber=0%如果函数仅输入4个参数值，则PointNumer默认值为100 PointNumer=100;endif nargin4%y0默认值为0 y0=0;endh=(xt-x0)/PointNumber;%计算所取的两离散点之间的距离x=x0+0:PointNumber*h;%表示出离散的自变量xy(1,:)=y0(:);for i=1:PointNumber%迭代计算过程 f1=h*feval(fun,x(i),y(i,:);f

12、1=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;function f=myfun012(x,y)f=y-2*x/y;myfun012.mfunction ex02()%ex02 比较改进欧拉法，简单欧拉方法以及微分方程符号解x3,y3=MyEulerPro(myfun012,0,1,1,10);x,y1=MyEuler(myfun012,0,1,1,10);%欧拉法所得的解for k=1:11t(k)=sqrt(1+2*x(k);endplot(x,y1

13、,-b,x3,y3,og,x,t,:*r)legend(简单欧拉法解,改进欧拉解,解析解)ex02.mnxnEuler方法yn改进Euler法yn精确解y(xn)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.26

14、49911.3416411.4142141.4832401.5491931.6124521.6733201.732051 在节点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阶方法阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.研究差

15、分公式阶的重要手段是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)4 Taylor展开方法展开方法

16、设y(x)是初值问题（4)的精确解,利用Taylor展开式可得称之为p阶Taylor展开方法.因此，可建立节点处近似值yn满足的差分公式其中所以,此差分公式是p阶方法.由于Taylor展开方法涉及很多复合函数(x,y(x)的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而,Taylor展开方法给出了一种构造单步显式高阶方法的途径.Euler方法可写为 可见,公式的局部截断误差为:y(xn+1)-yn+1=O(hp+1).5 Runge-Kutta方法方法 构造差分公式 改进的Euler方法可写为其中i,i,ij为待定参数.若此公式的局部截断误差为由于 yn+1=yn+

17、h1n+h2(n+hxn+hn yn)+O(h3)O(h3),称此公式为p p阶阶Runge-kuttaRunge-kutta方法方法,简称p p阶阶R-KR-K方法方法.对于p=2的情形,应有 =yn+h(1+2)n+h22(xn+n yn)+O(h3)所以,只要令 1+2=1,2=1/2,2=1/2 (8)一般地,参数由(8)确定的一族差分公式(7)统称为二二阶阶R-KR-K方法方法.称之为中点公式中点公式,或可写为若取=1,则得1=2=1/2,=1,此时公式(7)就是改进的Euler公式;若取1=0,则得2=1,=1/2,公式(7)为 高阶R-K公式可类似推导.下面列出常用的三阶、四阶R

18、-K公式.四阶标准四阶标准R-KR-K公式公式 三阶三阶R-KR-K公式公式 解解 四阶标准R-K公式为例例3 用四阶标准R-K方法求初值问题 y=y-2x/y ,0 x1 y(0)=1的数值解,取步长h=0.2.计算结果如下:function x,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,

19、y)if nargin 4|PointNum=0 PointNum=100;endif nargin 3 y0=0;endMyRunge_Kutta.my(1,:)=y0(:);%初值存为行向量形式h=(xt-x0)/(PointNum-1);%计算步长x=x0+0:PointNum*h;%得x向量值for k=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(

20、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)endfunction ex03()%用三种不同方法解微分方程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,y

21、e=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法解)hold onplot(x1,ye,bo,x2,yh,b+,x3,yr,r*)nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001

22、.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321 也可以构造隐式R-K方法,其一般形式为称之为p p级隐式级隐式R-KR-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

23、级2阶龙格库塔公式，ode45用的是以Runge-Kutta-Fehberg命名的 5级4阶公式。命令的输入f是待解方程写成的函数m文件：function dx=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分别为设定的相对误差和绝对误差。命令的输

24、出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);function dx=myfun02(x

25、,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,0 15,25 2);plot(x,y(:,1),-,x,y(:,2),*)%画出y(1),y(2)的函数图legend(y1函数图,y2函数图)用自编函数MyEulerPro（）MyRunge_Kutta（）求多阶微分方程clear,clcx0=1;%设定函数变量起始值xt=2;%设定函数变量终止值y0=1 10 30;%y初值N=20;x,yE=MyEulerPro(myfun03,x0,xt,y0,N)

26、;x,yR=MyRunge_Kutta(myfun03,x0,xt,y0,N);plot(x,yE(:,1),:r*,x,yR(:,1),-go)legend(改进欧拉法解,R-K法解)function dy=myfun03(x,y)dy=zeros(3,1)%初始化变量dydy(1)=y(2);%dy(1)表示y的一阶导数,其等于y的第二列值dy(2)=y(3);%dy(2)表示y的二阶导数dy(3)=2*y(3)/x3+3*y(2)/x3+3*exp(2*x)/x3%dy(3)表示y的三阶导数用ode23 ode45 ode113解多阶微分方程clear,clcx23,y23=ode23(

27、myfun03,1,2,1 10 30);x45,y45=ode45(myfun03,1,2,1 10 30);x113,y113=ode113(myfun03,1,2,1 10 30);figure(1)%第一幅图plot(x23,y23(:,1),*r,x45,y45(:,1),ob,x113,y113(:,1),+g)%作出各种函数所得结果legend(ode23解,ode45解,ode113解)title(ODE函数求解结果)figure(2)plot(x45,y45)%以ode45为例作出函数以及其各阶导数图legend(y,y一阶导数,y两阶导数)title(y,y一阶导数,y二阶导数函数图)