常微分方程的初值问题

常微分方程的初值问题第一页，共五十三页，2022年，8月28日一阶ODE问题的形式1、如果f与y无关，则计算变为第5章（数值积分）中讨论的任一种直接积分方法应用

初始条件2、如果f是y的函数，积分过程将不同于前者。若f是y

的线性函数，如：f=ay+b

其中a,b是常数或是

t的函数，此时原方程称为线性ODE若f不是线性函数，方程就称为非线性ODE。第二页，共五十三页，2022年，8月28日一、求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'''第三页，共五十三页，2022年，8月28日例：求微分方程的通解，并验证。>>

y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')

结果为y=(1/2*x^2+C1)*exp(-x^2)>>

symsx;diff(y)+2*x*y-x*exp(-x^2)

结果为ans=0注2：如果省略初值条件，则表示求通解；注3：如果省略自变量，则默认自变量为t

例：

dsolve('Dy=2*x')%dy/dt=2x结果为ans=2*x*t+C1第四页，共五十三页，2022年，8月28日例：求微分方程在初值条件下的特解，并画出解函数的图形。>>

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)例：第五页，共五十三页，2022年，8月28日注4：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。例：求微分方程组在初值条件下的特解[x,y]=dsolve('Dx+5*x+y=exp(t)',...'Dy-x-3*y=0','x(0)=1','y(0)=0','t')第六页，共五十三页，2022年，8月28日[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的输出个数只能为一个或与方程个数相等。第七页，共五十三页，2022年，8月28日例：用dsolve求解著名的VanderPol方程

>>symsmu;>>y=dsolve('D2y+mu*(y^2-1)*Dy+y')结果：Warning,cannotfindanexplicitsolutiony=&where(_a,{[diff(_b(_a),_a)*_b(_a)+mu*_b(_a)*_a^2-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：若找不到解析解，则返回其积分形式。只有很少一部分微分方程（组）能求出解析解。

大部分微分方程（组）只能利用数值方法求数值解。第八页，共五十三页，2022年，8月28日

Euler（欧拉）方法是求解一阶ODE的一种简便方法。尽管精度不高，但由于简单，特别适用于快速编程求解。它有三种格式：二、用Euler方法求数值解（a）向前Euler法（b）改进的Euler法（c）向后Euler法介绍这些方法是为了了解初值问题求解的基本思想。第九页，共五十三页，2022年，8月28日一阶ODE对两边从x0到x

积分得：

(积分方程）第十页，共五十三页，2022年，8月28日1、向前Euler法推导1：设节点为向前Euler法：用向前差分公式代替导数：此处，y(xn)表示xn处的理论解，yn表示y(xn)的近似解第十一页，共五十三页，2022年，8月28日一阶ODE对两边从xn

到xn+1

积分得：推导2：yn近似代替y(xn)用矩形代替右边的积分向前Euler法第十二页，共五十三页，2022年，8月28日例求出解析解和数值解并画图比较先用dsolve求解析解y=dsolve('Dy=y+2*x/(y^2)','y(0)=1','x')结果为y=1/3*(-18-54*x+45*exp(3*x))^(1/3)解析解：第十三页，共五十三页，2022年，8月28日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;forn=1:N-1x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));end下面再用向前欧拉法数值求解，为此先编写向前欧拉法的程序第十四页，共五十三页，2022年，8月28日最后通过图形比较用向前欧拉得到的数值解和解析解的误差

fun=inline('y+2*x/y^2','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])第十五页，共五十三页，2022年，8月28日向前欧拉方法截断误差为第十六页，共五十三页，2022年，8月28日对两边从xn

到xn+1

积分得：2、改进的Euler法yn近似代替y(xn)用梯形代替右边的积分梯形法第十七页，共五十三页，2022年，8月28日从n=0开始计算，每步都要求解一个关于yn+1的方程（一般是一个非线性方程），可用如下的迭代法计算：利用此算法，可得：利用（为允许误差）来控制是否继续迭代第十八页，共五十三页，2022年，8月28日迭代法太麻烦，实际上，当h取得很小时，只让上式中的第二式迭代一次就可以，即改进的Euler法（也叫欧拉预估—校正法）预估算式校正算式改进的Euler法=向前欧拉法+梯形法第十九页，共五十三页，2022年，8月28日向后Euler法依赖于向后差分近似，其格式为：

3、向后Euler法精度与向前欧拉法相同。如果f为非线性函数，则与改进的Euler算法一样，在每一步中需要采用迭代法。该方法有两个优点：

（a）绝对稳定；

（b）如果解为正，则可保证数值计算结果也为正。第二十页，共五十三页，2022年，8月28日三、龙格－库塔（Runge-kutta）方法

Euler法的一个重要缺陷是精度太低。为了获得高精度，就要减小h，这不仅会增加计算时间，也会产生舍入误差。1、二阶Runge-kutta方法第二十一页，共五十三页，2022年，8月28日或其实就是欧拉预估—校正方法（或改进的欧拉法）第二十二页，共五十三页，2022年，8月28日例用改进的欧拉法（即二阶龙格-库塔法）数值求解并与向前欧拉法、解析解画图比较前面已求出解析解：第二十三页，共五十三页，2022年，8月28日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;forn=1:N-1x(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先编写改进的欧拉法的程序：第二十四页，共五十三页，2022年，8月28日再分别调用Euler_bf.m和Euler_mend.m求解：

fun=inline('y+2*x/y^2','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])第二十五页，共五十三页，2022年，8月28日第二十六页，共五十三页，2022年，8月28日3、三阶龙格－库塔方法常见的三阶龙格－库塔方法的格式为：二阶龙格－库塔方法截断误差为精度不高，希望通过增加计算f的次数提高截断误差的阶数，为此引入第二十七页，共五十三页，2022年，8月28日4、四阶龙格－库塔方法常见的四阶龙格-库塔方法有两种，一种基于1/3辛普森法，格式：第二十八页，共五十三页，2022年，8月28日另一种基于3/8辛普森法，格式：第二十九页，共五十三页，2022年，8月28日例用二阶龙格-库塔法和四阶龙格-库塔法数值求解并与、解析解画图比较前面已求出解析解：先来编写四阶龙格-库塔法（基于1/3辛普森法）的程序：第三十页，共五十三页，2022年，8月28日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;forn=1:N-1x(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第三十一页，共五十三页，2022年，8月28日fun=inline('y+2*x/y^2','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求解：第三十二页，共五十三页，2022年，8月28日从图形上看，似乎二阶龙格—库塔法与四阶龙格-库塔法同样接近解析解，故再从数值结果比较看看哪种误差小第三十三页，共五十三页，2022年，8月28日err=[abs(y'-y1'),abs(y'-y2')]结果为err=000.00090.00000.00160.00000.00220.00000.00260.00000.00310.00000.00350.00000.00400.00000.00470.00000.00540.00000.00620.00000.00730.00000.00840.00000.00980.00000.01150.00000.01330.00000.01550.00000.01810.00000.02100.00000.02430.00000.02810.0000第三十四页，共五十三页，2022年，8月28日四、二阶ODE问题二阶ODE的一般形式为：其中是常数或是的函数，后两个方程为初始条件。如果与u无关，则该方程称为线性ODE;如果三者之中有一个是u或的函数，称为非线性的。下面着重研究用向前Euler法求解二阶ODE，及MATLAB程序。第三十五页，共五十三页，2022年，8月28日所以二阶ODE分解为两个一阶ODE：设：第三十六页，共五十三页，2022年，8月28日定义则上述两个一阶ODE写为：其向前Euler法的格式为：第三十七页，共五十三页，2022年，8月28日例求解二阶ODE解：设令则原方程的向量形式为：向前Euler法的格式为：第三十八页，共五十三页，2022年，8月28日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法求解第三十九页，共五十三页，2022年，8月28日第四十页，共五十三页，2022年，8月28日五、matlab相关命令数值求解ODE[X,Y]

=求解函数(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为方程中的参数第四十一页，共五十三页，2022年，8月28日[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第四十二页，共五十三页，2022年，8月28日求解函数ODE类型特点 说明ode45非刚性单步法；4，5阶R-K方法；累计截断误差为(△x)3大部分场合的首选方法ode23非刚性单步法；2，3阶R-K方法；累计截断误差为(△x)3使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到10-3～10-6

计算时间比ode45

短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear’s反向数值微分；精度中等若ode45

失效时，可尝试使用ode23s刚性单步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解函数，对于不同的ODE，可以调用不同的求解函数。第四十三页，共五十三页，2022年，8月28日求初值问题的数值解，求解范围为[0,0.5]例先定义需要求解的显式微分方程为一个函数functionyd=fun1(x,y)yd=-2*y+2*x^2+2*x最后在命令窗口调用函数求解[x,y]=ode23(‘fun1’,[0,0.5],1);第四十四页，共五十三页，2022年，8月28日fun2=inline('-2*y+2*x^2+2*x','x','y');求初值问题的数值解，求解范围为[0,0.5]例先定义需要求解的显式微分方程为一个函数在命令窗口用inline定义最后在命令窗口调用函数求解[x,y]=ode23(fun2,[0,0.5],1);第四十五页，共五十三页，2022年，8月28日如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。令

，则原方程可化为求解VerderPol初值问题例：前面演示过，该方程没有解析解，该方程不是一阶显式微分方程，故需要进行变换第四十六页，共五十三页，2022年，8月28日先用函数文件定义一阶显式微分方程组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：第四十七页，共五十三页，2022年，8月28日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