常微分方程的初值问题.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 solutiony=&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(xn)表示 xn 处的理论解，yn表示y(xn)的近似解,一阶ODE,对,两边从xn 到 xn+1 积分得：,推导2：,yn 近似代替y(xn),用矩形代替右边的积分,向前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),向前欧拉方法截断误差为,对,两边从xn 到 xn+1 积分得：,2、改进的Euler法,yn 近似代替y(xn),用梯形代替右边的积分,梯形法,从n=0开始计算，每步都要求解一个关于yn+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;endaxis(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 问题，因此MATLAB 提供了多种ODE求解函数，对于不同的ODE，可以调用不同的求解函数。,先定义需要求解的显式微分方程为一个函数,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);,先定义需要求解的显式微分方程为一个函数,在命令窗口用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);tocplot(t,y(:,1);figure;plot(t,y(:,2),结果：Elapsed time is 4.320568 seconds.,x1(t)曲线(变化比较平滑）,x2(t)曲线（在某些点上变化较快）,其实，许多传统的刚性问题可以采用普通函数来求解，不必刻意选择刚性问题的解法，当然，有的却必须用刚性问题的解法,