第八章 常微分方程的初值问题1.ppt

单击此处编辑母版标题样式,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第八章 常微分方程的初值问题,常微分方程（,ODE,）：只含有未知函数的导数的方程。,ODE,的阶数：最高阶导数的阶数,ODE,问题,初值问题：,边值问题：,已知初始点处的条件,已知初始点和末端点处的条件,本章只讨论初值问题,一阶,ODE,问题的形式,1,、如果,f,与,y,无关，则计算变为第,5,章（数值积分）中讨论的任一种直接积分方法应用,初始条件,2,、如果,f,是,y,的函数，积分过程将不同于前者。,若,f,是,y,的线性函数，如：,f=,ay+b,其中,a,b,是常数或是,t,的函数，,此时原方程称为线性,ODE,若,f,不是线性函数，方程就称为非线性,ODE,。,一、求,ODE,的解析解,dsolve,输出变量列表,=dsolve,(eq1,eq2,.,eqn,cond1,cond2,.,condn,v1,v2,vn,),其中,eq1,、,eq2,、,.,、,eqn,为微分方程，,cond1,、,cond2,、,.,、,condn,为初值条件，,v1,v2,vn,为自变量。,注,1,：,微分方程中用,D,表示对,自变量,的导数，如：,Dy,y,；,D2y y,；,D3y y,例：,求微分方程 的通解，并验证。,y=dsolve(Dy+2*x*y=x*exp(-x2),x),结果为,y=(1/2*x2+C1)*exp(-x2),syms,x;diff(y)+2*x*y-x*exp(-x2),结果为,ans,=0,注,2,：,如果省略初值条件，则表示求通解；,注,3,：,如果省略自变量，则默认自变量为,t,例：,dsolve(Dy,=2*x)%,dy/dt,=2x,结果为,ans,=,2*x*t+C1,例：,求微分方程 在初值条件 下的特解，并画出解函数的图形。,y=,dsolve(x,*,Dy+y-exp(x,)=0,y(1)=2*exp(1),x),结果为,y=,(exp(x)+exp(1)/x,ezplot(y,)%,此时的,y,已经是符号变量，故不用,ezplot(y,),dsolve(D3y=-y,y(0)=1,Dy(0)=0,D2y(0)=0),结果为,ans=,1/3*exp(-t)+2/3*exp(1/2*t)*cos(1/2*3(1/2)*t),例：,注,4,：,解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构,（,structure,）,类型的数据。,例：,求微分方程组,在初值条件,下的特解,x,y=dsolve(Dx+5*,x+y,=,exp(t,),.,Dy-x-3*y=0,x(0)=1,y(0)=0,t),x,y,=dsolve(Dx+5*,x+y,=exp(t),Dy-x-3*y=0,.,x(0)=1,y(0)=0,t),或,r=dsolve(Dx+5*,x+y,=,exp(t,),Dy-x-3*y=0,.,x(0)=1,y(0)=0,t),这里返回的,r,是一个 结构类型 的数据,r.x%,查看解函数,x(t),r.y%,查看解函数,y(t),dsolve,的输出个数只能为一个 或 与方程个数相等。,例：,用,dsolve,求解著名的,Van,der,Pol,方程,syms mu;,y=dsolve(D2y+mu*(y2-1)*Dy+y),结果：,Warning,cannot find an explicit solution,y=,&,where(_a,diff(_b(_a),_a,)*_,b(_a)+mu,*_,b(_a,)*_a2-mu*_,b(_a)+_a,=0,_,b(_a,)=,diff(y(t),t,),_a=,y(t,),y(t,)=_a,t=Int(1/_b(_a),_a)+C1),注,5,：,若找不到解析解，则返回其积分形式。,只有很少一部分微分方程（组）能求出解析解。大部分微分方程（组）只能利用,数值方法,求数值解。,Euler,（,欧拉）方法是求解一阶,ODE,的一种简,便方法。尽管精度不高，但由于简单，特别适用于快速编程求解。它有三种格式：,二、用,Euler,方法求数值解,（,a,）,向前,Euler,法,（,b,）,改进的,Euler,法,（,c,）,向后,Euler,法,介绍这些方法是为了了解初值问题求解的基本思想。,一阶,ODE,对,两边从,x0,到,x,积分得：,(,积分方程）,1,、向前,Euler,法,推导,1,：,设节点为,向前,Euler,法：,用向前差分公式代替导数：,此处，,y,(,x,n,),表示,x,n,处的理论解，,y,n,表示,y,(,x,n,),的近似解,一阶,ODE,对,两边从,x,n,到,x,n+1,积分得：,推导,2,：,y,n,近似代替,y,(,x,n,),用矩形代替右边的积分,向前,Euler,法,例,求出,解析解和数值解并画图比较,先用,dsolve,求解析解,y=dsolve(Dy=y+2*x/(y2),y(0)=1,x),结果为,y=1/3*(-18-54*x+45*exp(3*x)(1/3),解析解：,function x,y=Euler_bf(fun,x0,y0,xmax,h),%fun,为显式一阶微分方程右端的函数,%x0,y0,为初始条件,满足,y(x0)=y0,%xmax,为,x,的取值最大值,%h,为将,x,等分后的步长,N=(xmax-x0)/h+1;%N,为总的节点数,x(1)=x0;y(1)=y0;,for n=1:N-1,x(n+1)=x(n)+h;,y(n+1)=y(n)+h*feval(fun,x(n),y(n);,end,下面再用向前欧拉法数值求解，为此先编写向前欧拉法的程序,最后通过图形比较用向前欧拉得到的数值解和解析解的误差,fun=inline(y+2*x/y2,x,y);,x1,y1=Euler_bf(fun,0,1,2,0.1);,x2,y2=Euler_bf(fun,0,1,2,0.05);,x=0:0.01:2;,y=1/3*(-18-54*x+45*exp(3*x).(1/3);,plot(x1,y1,*,x2,y2,x,x,y),axis(0,2,0,9),向前欧拉方法截断误差为,对,两边从,x,n,到,x,n+1,积分得：,2,、改进的,Euler,法,y,n,近似代替,y,(,x,n,),用梯形代替右边的积分,梯形法,从,n=0,开始计算，每步都要求解一个关于,y,n+1,的方程（一般是一个非线性方程），可用如下的迭代法计算：,利用此算法，可得：,利用,（为允许误差）来控制是否继续迭代,迭代法太麻烦，,实际上，当,h,取得很小时，只让上式中的第二式迭代一次就可以，即,改进的,Euler,法（也叫欧拉预估,校正法）,预估算式,校正算式,改进的,Euler,法,=,向前欧拉法,+,梯形法,向后,Euler,法依赖于向后差分近似，其格式为：,3,、向后,Euler,法,精度与向前欧拉法相同。如果,f,为非线性函数，则与改进的,Euler,算法一样，在每一步中需要采用迭代法。该方法有两个优点：（,a,）,绝对稳定；（,b,）,如果解为正，则可保证数值计算结果也为正。,三、龙格库塔（,Runge-kutta,）,方法,Euler,法的一个重要缺陷是精度太低。为了获得,高精度，就要减小,h,，,这不仅会增加计算时间，也,会产生舍入误差。,1,、二阶,Runge-kutta,方法,或,其实就是欧拉预估,校正方法（或改进的欧拉法）,例,用改进的欧拉法（即二阶龙格,-,库塔法）数值求解并与向前欧拉法、解析解画图比较,前面已求出解析解：,function x,y=Euler_mend(fun,x0,y0,xmax,h),%fun,为显式一阶微分方程右端的函数,%x0,y0,为初始条件,满足,y(x0)=y0,%xmax,为,x,的取值最大值,%h,为将,x,等分后的步长,N=(xmax-x0)/h+1;%N,为总的节点数,x(1)=x0;y(1)=y0;,for n=1:N-1,x(n+1)=x(n)+h;,k1=h*feval(fun,x(n),y(n);,k2=h*feval(fun,x(n+1),y(n)+k1);,y(n+1)=y(n)+1/2*(k1+k2);,end,先编写改进的欧拉法的程序：,再分别调用,Euler_bf.m,和,Euler_mend.m,求解：,fun=inline(y+2*x/y2,x,y);,x1,y1=Euler_bf(fun,0,1,2,0.1);,x2,y2=Euler_mend(fun,0,1,2,0.1);,x=0:0.01:2;,y=1/3*(-18-54*x+45*exp(3*x).(1/3);,plot(x1,y1,*,x2,y2,s,x,y),axis(0,2,0,9),3,、三阶龙格库塔方法,常见,的三阶龙格库塔方法的格式为：,二阶龙格库塔方法截断误差为,精度不高，希望通过增加计算,f,的次数提高截断误差的阶数，为此引入,4,、四阶龙格库塔方法,常见,的四阶龙格,-,库塔方法有两种，,一种基于,1/3,辛普森法，格式：,另一种基于,3/8,辛普森法，格式：,例,用二阶龙格,-,库塔法和四阶龙格,-,库塔法数值求解并与、解析解画图比较,前面已求出解析解：,先来编写四阶龙格,-,库塔法（基于,1/3,辛普森法）的程序：,function x,y=RK4(fun,x0,y0,xmax,h),%fun,为显式一阶微分方程右端的函数,%x0,y0,为初始条件,满足,y(x0)=y0,%xmax,为,x,的取值最大值,%h,为将,x,等分后的步长,N=(xmax-x0)/h+1;%N,为总的节点数,x(1)=x0;y(1)=y0;,for n=1:N-1,x(n+1)=x(n)+h;,k1=h*feval(fun,x(n),y(n);,k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2);,k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2);,k4=h*feval(fun,x(n+1),y(n)+k3);,y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);,end,fun=inline(y+2*x/y2,x,y);,x1,y1=Euler_mend(fun,0,1,2,0.1);,x2,y2=RK4(fun,0,1,2,0.1);,x=0:0.1:2;,y=1/3*(-18-54*x+45*exp(3*x).(1/3);,plot(x1,y1,*,x2,y2,s,x,y),axis(0,2,0,9),再分别调用,Euler_mend.m,和,RK4.m,求解：,从图形上看，似乎二阶龙格,库塔法与四阶龙格,-,库塔法同样接近解析解，故再从数值结果比较看看哪种误差小,err=abs(y-y1),abs(y-y2),结果为,err=0 0,0.0009 0.0000,0.0016 0.0000,0.0022 0.0000,0.0026 0.0000,0.0031 0.0000,0.0035 0.0000,0.0040 0.0000,0.0047 0.0000,0.0054 0.0000,0.0062 0.0000,0.0073 0.0000,0.0084 0.0000,0.0098 0.0000,0.0115 0.0000,0.0133 0.0000,0.0155 0.0000,0.0181 0.0000,0.0210 0.0000,0.0243 0.0000,0.0281 0.0000,四、二阶,ODE,问题,二阶,ODE,的一般形式为：,其中 是常数或是 的函数，后两个方程为初,始条件。,如果 与,u,无关，则该方程称为线性,ODE;,如果三者之中有一个是,u,或 的函数，称为非线性的。,下面着重研究用向前,Euler,法求解二阶,ODE,，及,MATLAB,程序。,所以二阶,ODE,分解为两个一阶,ODE,：,设：,定义,则,上述两个一阶,ODE,写为：,其,向前,Euler,法的格式为：,例,求解二阶,ODE,解：设,令,则原方程的向量形式为：,向前,Euler,法的格式为：,h=0.05;t_max=5;n=1;,y(:,1)=0;1;t(1)=0;,while t(n)t_max,y(:,n+1)=y(:,n)+h*,f_def(y(:,n),t,);,t(n+1)=t(n)+h;,n=n+1;,end,axis(0 5-1 1);plot(t,y(1,:),t,y(2,:),:),function f=f_def(y,t),f=y(2);(-5*abs(y(2)*y(2)-20*y(1);,先用函数文件定义,f(u,v,t,),再用向前,Euler,法求解,五、,matlab,相关命令数值求解,ODE,X,Y,=,求解函数,(fun,x0,xf,y0,option,，,p1,p2,.),其中：,（,1,）,fun,是用,inline,或函数文件定义的显式常微分方程的,函数名,函数文件格式：,或,inline,格式：,function yd=,函数名,(x,y,flag,p1,p2,),yd=.,函数名,=inline(,显式方程,x,y,flag,p1,p2,.),x,为自变量，,y,为因变量，,yd,为,y,的导数,flag,用于控制求解过程，,p1,，,p2,为方程中的参数,X,Y,=,求解函数,(fun,x0,xf,y0,option,，,p1,p2,.),其中：,（,2,）,x0,xf,为求解区间，,y0,为初值条件；,（,3,）,option,（可省略）为控制选项（用于设置误差限，求解方程最大允许步长等等），,（,4,）,p1,p2,为微分方程中的附加参数,（,5,）,Matlab,在数值求解时自动对求解区间进行分割，,X(,向量,),中返回的是分割点的值,(,自变量,),，,Y(,向量,),中返回的是解函数在这些分割点上的函数值。,（,6,）,求解函数可以是,ode45,、,ode23,、,ode113,、,ode15s,、,ode23s,、,ode23t,、,ode23tb,求解函数,ODE,类型,特点,说明,ode45,非刚性,单步法；,4,，,5,阶,R-K,方法；累计截断误差为,(x),3,大部分场合的首选方法,ode23,非刚性,单步法；,2,，,3,阶,R-K,方法；累计截断误差为,(x),3,使用于精度较低的情形,ode113,非刚性,多步法；,Adams,算法；高低精度均可到,10,-3,10,-6,计算时间比,ode45,短,ode23t,适度刚性,采用梯形算法,适度刚性情形,ode15s,刚性,多步法；,Gears,反向数值微分；精度中等,若,ode45,失效时，可尝试使用,ode23s,刚性,单步法；,2,阶,Rosebrock,算法；低精度,当精度较低时，计算时间比,ode15s,短,ode23tb,刚性,梯形算法；低精度,当精度较低时，计算时间比,ode15s,短,没有一种算法可以有效地解决所有的,ODE,问题，因此,MATLAB,提供了多种,ODE,求解函数，对于不同的,ODE,，可以调用不同的,求解函数,。,求初值问题,的数值解，求解范围为,0,0.5,例,先定义需要求解的显式微分方程为一个函数,function yd=fun1(x,y),yd=-2*y+2*x2+2*x,最后在命令窗口调用函数求解,x,y=ode23(fun1,0,0.5,1);,fun2=inline(-2*y+2*x2+2*x,x,y);,求初值问题,的数值解，求解范围为,0,0.5,例,先定义需要求解的显式微分方程为一个函数,在命令窗口用,inline,定义,最后在命令窗口调用函数求解,x,y=ode23(fun2,0,0.5,1);,如果需求解的问题是,高阶,常微分方程，则需,将其化为一阶常微分方程组,，此时需用,函数文件来定义该常微分方程组。,令,，则原方程可化为,求解,Ver,der,Pol,初值问题,例：,前面演示过，该方程没有解析解，该方程不是一阶显式微分方程，故需要进行变换,先用函数文件定义一阶显式微分方程组,function y=vdp_eq1(t,x,flag,mu),y=x(2);,-mu*(x(1)2-1)*x(2)-x(1);,再编写脚本文件,vdpl.m,，在命令窗口直接运行该文件。,x0=-0.2;-0.7;,tf,=20;,mu=1;t1,y1=ode45(,vdp_eq1,0,tf,x0,mu);,mu=2;t2,y2=ode45(,vdp_eq1,0,tf,x0,mu);,plot(t1,y1,t2,y2,:),figure;,plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:),法,1,：,vdp_eq2=inline(x(2);-mu*(x(1)2-1)*x(2)-x(1),t,x,flag,mu,);,x0=-0.2;-0.7;,tf,=20;,mu=1;t1,y1=ode45(,vdp_eq2,0,tf,x0,mu);,mu=2;t2,y2=ode45(,vdp_eq2,0,tf,x0,mu);,plot(t1,y1,t2,y2,:),figure;,plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:),两种结果一样,法,2,用,inline,定义方程（注意调用格式不同）,(a),不同,u,下的时间响应曲线,(b),相平面曲线,下面，改变,u,的值，并设仿真终止时间为,3000,，看看计算结果如何,vdp_eq2=inline(x(2);-mu*(x(1)2-1)*x(2)-x(1),t,x,flag,mu,);,x0=-0.2;-0.7;,tf,=3000;,mu=1000;t,y=ode45(vdp_eq2,0,tf,x0,mu);,plot(t1,y1,t2,y2,:),发现，经过很长时间的等待，得出一个错误信息,.,事实上，由于变步长所采用的步长过小，而要求仿真的终止时间比较大，导致输出的,y,矩阵过大，超出了计算机存储空间的容限，所以这个问题不适合用,ode45,来求解,应该采用,ode15s,来求解,在许多领域，常会遇到一类特殊的常微分方程，其中一些解变化缓慢，另一些解变化快，且相差比较悬殊，这些方程称为,刚性方程（,Stiff,方程）,刚性方程一般不适合用,ode45,求解，应采用,ode15s,求解,vdp_eq2=inline(x(2);-mu*(x(1)2-1)*x(2)-x(1),t,x,flag,mu,);,h_opt,=odeset(RelTol,1e-6);,x0=-0.2;-0.7;,tf,=3000;,tic,mu,=1000;t,y=,ode15s,(vdp_eq2,0,tf,x0,h_opt,mu);,toc,plot(t,y(:,1);,figure;,plot(t,y(:,2),结果：,Elapsed time is 4.320568 seconds.,x1(t),曲线,(,变化比较平滑）,x2(t),曲线,（在某些点上变化较快）,其实，许多传统的刚性问题可以采用普通函数来求解，不必刻,意选择刚性问题的解法，当然，有的却必须用刚性问题的解法,