第八章-常微分方程的初值问题VIP

第八章常微分方程的初值问题,常微分方程（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表示对自变量的导数，如：,Dyy；D2yy；D3yy,例：求微分方程的通解，并验证。,y=dsolve(Dy+2*x*y=x*exp(-x2),x)结果为y=(1/2*x2+C1)*exp(-x2),symsx;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)/xezplot(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求解著名的VanderPol方程,symsmu;y=dsolve(D2y+mu*(y2-1)*Dy+y)结果：Warning,cannotfindanexplicitsolutiony=M=70;g=9.8;x=0;y=0;h=0.1;n=0;%n为迭代次数xr(1)=t;yr(1)=v;,whilex20n=n+1;y=y+h*(-C/M*y2+g);x=x+h;yr(n+1)=v;xr(n+1)=x;endplot(xr,yr);axis(0,20,0,60),尽管向前Euler法非常简单，但有两种误差效应：第一种是截断误差，第二种误差源于计算的不稳定性。当方程的时间系数为负值，h又不是很小时，第二种误差就会增加。带有负时间系数的常见方程为其精确解为该问题的向前Euler格式为：,如果，数值解为正的趋于零；如果，随n的增加，解的符号交替改变；如果，则在每一步符号变化后，解的数量级增加，最终导致结果发散。,对,两边从xn到xn+1积分得：,2、改进的Euler法,yn近似代替y(xn),用梯形代替右边的积分,梯形法,从n=0开始计算，每步都要求解一个关于yn+1的方程（一般是一个非线性方程），可用如下的迭代法计算：,利用此算法，可得：,利用,（为允许误差）来控制是否继续迭代,迭代法太麻烦，实际上，当h取得很小时，只让上式中的第二式迭代一次就可以，即,改进的Euler法（也叫欧拉预估校正法）,预估算式,校正算式,改进的Euler法=向前欧拉法+梯形法,例2初值问题：,求解的值。,向前Euler法：,yf(1)=10;t(1)=0;h=0.1;n=1;fprintf(tyfn);fprintf(%3.1f%6.4fn,t(1),yf(1);whilet(n)1yf(n+1)=yf(n)+h*(-yf(n)1.5+1);t(n+1)=t(n)+h;n=n+1;fprintf(%3.1f%6.4fn,t(n),yf(n);end,迭代5步结果：,tyf0.010.00000.16.93770.25.21040.34.12100.43.38440.52.8618,此题无解析解,ym(1)=10;t(1)=0;h=0.1;n=1;fprintf(tymn);fprintf(%3.1f%6.4fn,t(1),ym(1);whilet(n)1y0(n+1)=ym(n)+h*(-ym(n)1.5+1);ym(n+1)=ym(n)+h/2*(-y0(n+1)1.5-ym(n)1.5+2);t(n+1)=t(n)+h;n=n+1;fprintf(%3.1f%6.4fn,t(n),ym(n);end,改进Euler法：,迭代5步结果：,两种方法的比较图：,plot(t,yf,-,t,ym,-);axis(0,1.2,0,10),tym0.010.00000.17.60520.25.99250.34.86160.44.04210.53.4320,向后Euler法依赖于向后差分近似，其格式为：,3、向后Euler法,精度与向前欧拉法相同。如果f为非线性函数，则与改进的Euler算法一样，在每一步中需要采用迭代法。该方法有两个优点：（a）绝对稳定；（b）如果解为正，则可保证数值计算结果也为正。,4、二阶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;whilet(n)t_maxy(:,n+1)=y(:,n)+h*f_def(y(:,n),t);t(n+1)=t(n)+h;n=n+1;endaxis(05-11);plot(t,y(1,:),t,y(2,:),:),functionf=f_def(y,t)f=y(2);(-5*abs(y(2)*y(2)-20*y(1);,先用函数文件定义f(u,v,t),再用向前Euler法求解,5、高阶ODE问题,设方程:,如果是常数或t的函数，则为线性的；,如果至少有一个是函数，则为非线性的；,初始条件：,上述微分方程可转化为四个一阶ODE：,向量形式为：,向前Euler法的格式为：,数值方法同样适用于积分方程。,如,设,对v微分：,则原方程变为：,向量形式为：,三、龙格库塔（Runge-kutta）方法,Euler法的一个重要缺陷是精度太低。为了获得高精度，就要减小h，这不仅会增加计算时间，也会产生舍入误差。,龙格库塔方法采用一种内迭代法来提高精度。,考虑：,积分：,将上式右端应用数值积分方法可得龙格库塔方法。,1、二阶Runge-kutta方法,将梯形公式用于上式积分,其中,由于右端未知，可以用向前Euler法所得的估计值来代替，由此得到的方法称为二阶龙格库塔方法。,得：,常见的二阶龙格库塔方法的形式为：,或,该方法其实就是欧拉预估校正方法,例4用二阶龙格库塔方法求解,其中,y(1)=0;t(1)=0;h=1e-4;n=1;whilet0.02k1=h*f_x10(y(n);k2=h*f_x10(y(n)+k1);y(n+1)=y(n)+1/2*(k1+k2);t(n+1)=t(n)+h;n=n+1;endplot(t,y)axis(0,0.02,0,0.6),functionf=f_x10(y)f=-400*y+200;,二阶龙格库塔方法截断误差为,精度不高，希望通过增加计算f的次数提高截断误差的阶数，为此引入,3、三阶龙格库塔方法,三阶龙格库塔方法的格式为：,4、四阶龙格库塔方法,常见四阶方法有两种，一种基于1/3辛普森法，格式：,另一种基于3/8辛普森法，格式：,例用四阶龙格库塔方法求解,其中,functionf=f_x10(y)f=-400*y+200;,下面用基于1/3辛普森法的龙格库塔方法求解,y(1)=0;t(1)=0;h=1e-4;n=1;whilet(n)0.02k1=h*f_x10(y(n);k2=h*f_x10(y(n)+k1/2);k3=h*f_x10(y(n)+k2/2);k4=h*f_x10(y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);t(n+1)=t(n)+h;n=n+1;endplot(t,y);axis(0,0.02,0,0.6);,四、matlab相关命令数值求解ODE,T,Y=solver(fun,x0,xf,y0,option，p1,p2,.),其中：（1）fun是用inline或函数文件定义的显式常微分方程的函数名,函数文件格式：,或inline格式：,functionyd=函数名(x,y,flag,p1,p2,)yd=.,函数名=inline(显式方程,x,y,flag,p1,p2,.),x为自变量，y为因变量，yd为y的导数,flag用于控制求解过程，p1，p2为方程中的参数,T,Y=solver(fun,x0,xf,y0,option，p1,p2,.),其中：（2）x0,xf为求解区间，y0为初值条件；（3）option（可省略）为控制选项（用于设置误差限，求解方程最大允许步长等等），（4）p1,p2为微分方程中的附加参数（5）Matlab在数值求解时自动对求解区间进行分割，T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。（6）solver为ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。,fun2=inline(-2*y+2*x2+2*x,x,y);,先定义需要求解的显式微分方程为一个函数,functionyd=fun1(x,y)yd=-2*y+2*x2+2*x,或在命令窗口用inline定义,最后在命令窗口调用函数求解,x,y=ode23(fun1,0,0.5,1);,x,y=ode23(fun2,0,0.5,1);,如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。,令，则原方程可化为,求解VerderPol初值问题,例：,前面演示过，该方程没有解析解，该方程不是一阶显式微分方程，故需要进行变换,先用函数文件定义一阶显式微分方程组,functiony=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,