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

《第八章 常微分方程的初值问题1.ppt》由会员分享，可在线阅读，更多相关《第八章 常微分方程的初值问题1.ppt（53页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第八章 常微分方程的初值问题常微分方程（常微分方程（ODE）：只含有未知函数的导数的方程。：只含有未知函数的导数的方程。ODE的阶数：的阶数：最高阶导数的阶数最高阶导数的阶数ODE问题问题初值问题：初值问题：边值问题：边值问题：已知初始点处的条件已知初始点处的条件已知初始点和末端点处的条件已知初始点和末端点处的条件本章只讨论初值问题本章只讨论初值问题 一阶一阶ODE问题的形式问题的形式 1、如果、如果 f 与与 y 无关，则计算变为第无关，则计算变为第5章（数值积分）中章（数值积分）中讨论的任一种直接积分方法应用讨论的任一种直接积分方法应用 初始条件初始条件2、如果、如果 f 是是 y 的函数

2、的函数，积分过程将不同于前者。，积分过程将不同于前者。若若 f 是是 y 的线性函数，如：的线性函数，如：f=ay+b 其中其中a,b是常数或是是常数或是 t 的函数，的函数，此时原方程称为线性此时原方程称为线性ODE若若 f 不是线性函数，方程就称为非线性不是线性函数，方程就称为非线性ODE。一、求一、求ODEODE的解析解的解析解dsolvedsolve输出变量列表输出变量列表=dsolve(eq1,eq2,.,eqn,cond1,cond2,.,condn,v1,v2,vn)其中其中 eq1、eq2、.、eqn为微分方程，为微分方程，cond1、cond2、.、condn为初值条件，为初

3、值条件，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(

4、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)例：

5、例：注注4 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)或

6、或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

7、(_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的一种简的一种简便方法。尽管精度不高，但由于简单，特别适用于便

8、方法。尽管精度不高，但由于简单，特别适用于快速编程求解。它有三种格式：快速编程求解。它有三种格式：二、用二、用Euler Euler 方法求数值解方法求数值解（a）向前向前Euler 法法（b）改进的改进的Euler 法法（c）向后向后Euler法法介绍这些方法是为了了解初值问题求解的基本思想。介绍这些方法是为了了解初值问题求解的基本思想。一阶一阶ODE 对对 两边从两边从x0 到到 x 积分得：积分得：(积分方程）积分方程）1 1、向前、向前EulerEuler法法推导推导1：设节点为设节点为向前向前Euler法：法：用向前差分公式代替导数：用向前差分公式代替导数：此处，此处，y(xn)表示

9、表示 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为显

10、式一阶微分方程右端的函数为显式一阶微分方程右端的函数%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下面再用向前欧拉法数值求解，为此先编写向前欧拉法下面再用向前欧拉法数值求解，为此先编写向前欧拉法的程序的程序最后通过图形比较用向前欧拉得到的数值解和解析解的最后通过图形比较用向前欧拉得到的数值解

11、和解析解的误差误差 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开始计算，每步

12、都要求解一个关于开始计算，每步都要求解一个关于yn+1的方程的方程（一般是一个非线性方程），可用如下的迭代法计算：（一般是一个非线性方程），可用如下的迭代法计算：利用此算法，可得：利用此算法，可得：利用利用 （为允许误差）来控制是否继续为允许误差）来控制是否继续迭代迭代迭代法太麻烦，迭代法太麻烦，实际上，当实际上，当h取得很小时，只让上式中取得很小时，只让上式中的第二式迭代一次就可以，即的第二式迭代一次就可以，即改进的改进的Euler法（也叫欧拉预估法（也叫欧拉预估校正法）校正法）预估算式预估算式校正算式校正算式改进的改进的Euler法法=向前欧拉法向前欧拉法+梯形法梯形法向后向后Euler法

13、依赖于向后差分近似，其格式为：法依赖于向后差分近似，其格式为：3 3、向后、向后EulerEuler法法 精度与向前欧拉法相同。如果精度与向前欧拉法相同。如果 f 为非线性函数，则与为非线性函数，则与改进的改进的Euler算法一样，在每一步中需要采用迭代法。算法一样，在每一步中需要采用迭代法。该方法有两个优点：该方法有两个优点：（a）绝对稳定；绝对稳定；（b）如果解为正，则可保证数值计算结果也为正。如果解为正，则可保证数值计算结果也为正。三、龙格库塔（三、龙格库塔（Runge-kutta）方法方法 Euler法的一个重要缺陷是精度太低。为了获得法的一个重要缺陷是精度太低。为了获得高精度，就要减

14、小高精度，就要减小 h，这不仅会增加计算时间，也这不仅会增加计算时间，也会产生舍入误差。会产生舍入误差。1 1、二阶、二阶Runge-kutta方法方法或或其实就是欧拉预估其实就是欧拉预估校正方法（或改进的欧拉法）校正方法（或改进的欧拉法）例例 用改进的欧拉法（即二阶用改进的欧拉法（即二阶龙格龙格-库塔法）数值求解并库塔法）数值求解并与向前欧拉法、解析解画与向前欧拉法、解析解画图比较图比较前面已求出解析解：前面已求出解析解：function x,y=Euler_mend(fun,x0,y0,xmax,h)%fun为显式一阶微分方程右端的函数为显式一阶微分方程右端的函数%x0,y0为初始条件为初

15、始条件,满足满足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)

16、;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 3、三阶龙格库塔方法、三阶龙格库塔方法常见常见的三阶龙格库塔方法的格式为：的三阶龙格库塔方法的格式为：二阶龙格库塔方法截断误差为二阶龙格库塔方法截断误差为精度不高，希望通过增加计算精度不高，希望通过增加计算f的次数提高截断误差的的次数提高截断误差的阶数，为此引入阶数，为此引入4 4、四阶龙格库塔方法、四

17、阶龙格库塔方法常见常见的四阶龙格的四阶龙格-库塔方法有两种，库塔方法有两种，一种基于一种基于1/3辛普森法，格式：辛普森法，格式：另一种基于另一种基于3/8辛普森法，格式：辛普森法，格式：例例 用二阶龙格用二阶龙格-库塔法和四阶库塔法和四阶龙格龙格-库塔法数值求解并与、库塔法数值求解并与、解析解画图比较解析解画图比较前面已求出解析解：前面已求出解析解：先来编写四阶龙格先来编写四阶龙格-库塔法（基于库塔法（基于1/31/3辛普森法）的程序：辛普森法）的程序：function x,y=RK4(fun,x0,y0,xmax,h)%fun为显式一阶微分方程右端的函数为显式一阶微分方程右端的函数%x0,

18、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);

19、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求解：求解：从图形上看，似乎二阶龙格从图形上看，似乎二阶龙格库塔法与四阶龙格库塔法与四阶龙格-库塔法库塔法同样接近解析解，故再从数值结果比较看看哪种误差小同样接近解析解，故再从数值结果比较看看哪种误差小er

20、r=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.000

21、0 0.0281 0.0000四、二阶ODE问题 二阶二阶ODE的一般形式为：的一般形式为：其中其中 是常数或是是常数或是 的函数，后两个方程为初的函数，后两个方程为初始条件。始条件。如果如果 与与u无关，则该方程称为线性无关，则该方程称为线性ODE;如果三者之中有一个是如果三者之中有一个是u或或 的函数，称为非线性的。的函数，称为非线性的。下面着重研究用向前下面着重研究用向前Euler法求解二阶法求解二阶ODE，及，及MATLAB程序。程序。所以二阶所以二阶ODE分解为两个一阶分解为两个一阶ODE：设：设：定义定义则则上述两个一阶上述两个一阶ODE写为：写为：其其向前向前Euler法的格式为

22、：法的格式为：例例 求解二阶求解二阶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)再

23、用向前再用向前Euler法求解法求解五、五、matlab相关命令数值求解相关命令数值求解ODEX,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用于控制求解过

24、程，用于控制求解过程，p1，p2为方程中的参数为方程中的参数X,Y=求解函数求解函数(fun,x0,xf,y0,option，p1,p2,.)其中：其中：（2）x0,xf为求解区间，为求解区间，y0 为初值条件；为初值条件；（3）option（可省略）为控制选项（用于设置误差限，（可省略）为控制选项（用于设置误差限，求解方程最大允许步长等等），求解方程最大允许步长等等），（4）p1,p2为微分方程中的附加参数为微分方程中的附加参数（5）Matlab在数值求解时自动对求解区间进行分割，在数值求解时自动对求解区间进行分割，X(向量向量)中返回的是分割点的值中返回的是分割点的值(自变量自变量)，Y(

25、向量向量)中返回的是解函数在这些分割点上的函数值。中返回的是解函数在这些分割点上的函数值。（6）求解函数可以是求解函数可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb求解函数求解函数 ODE类型类型特点特点说明说明ode45非刚性非刚性单步法；单步法；4，5 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为(x)3大部分场合的首选方法大部分场合的首选方法ode23非刚性非刚性单步法；单步法；2，3 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为(x)3使用于精度较低的情形使用于精度较低的情形ode113非刚性非刚性多步法；多步

26、法；Adams算法；高低精算法；高低精度均可到度均可到 10-310-6计算时间比计算时间比 ode45 短短ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15s刚性刚性多步法；多步法；Gears 反向数值微反向数值微分；精度中等分；精度中等若若 ode45 失效时，可失效时，可尝试使用尝试使用ode23s刚性刚性单步法；单步法；2 阶阶Rosebrock 算算法；低精度法；低精度当精度较低时，计算时当精度较低时，计算时间比间比 ode15s 短短ode23tb刚性刚性梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，计算时间比间比ode1

27、5s短短没有一种算法可以有效地解决所有的没有一种算法可以有效地解决所有的 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

28、+2*x2+2*x,x,y);求初值问题求初值问题 的数值解，求解范围的数值解，求解范围为为 0,0.5例例 先定义需要求解的显式微分方程为一个函数先定义需要求解的显式微分方程为一个函数在命令窗口用在命令窗口用inline定义定义最后在命令窗口调用函数求解最后在命令窗口调用函数求解x,y=ode23(fun2,0,0.5,1);如果需求解的问题是如果需求解的问题是高阶高阶常微分方程，则需常微分方程，则需将其化为一将其化为一阶常微分方程组阶常微分方程组，此时需用，此时需用函数文件来定义该常微分方函数文件来定义该常微分方程组。程组。令令 ，则原方程可化为，则原方程可化为 求解求解 Ver der

29、Pol 初值问题初值问题例例：前面演示过，该方程没有解析解，该方程不是一阶前面演示过，该方程没有解析解，该方程不是一阶显式微分方程，故需要进行变换显式微分方程，故需要进行变换先用函数文件定义一阶显式微分方程组先用函数文件定义一阶显式微分方程组function y=vdp_eq1(t,x,flag,mu)y=x(2);-mu*(x(1)2-1)*x(2)-x(1);l 再编写脚本文件再编写脚本文件 vdpl.m，在命令窗口直接运行该文件。，在命令窗口直接运行该文件。x0=-0.2;-0.7;tf=20;mu=1;t1,y1=ode45(vdp_eq1,0,tf,x0,mu);mu=2;t2,y2

30、=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)

31、,y2(:,2),:)两种结果一样两种结果一样法法2 2 用用inlineinline定义方程（注意调用格式不同）定义方程（注意调用格式不同）(a)不同不同u下的时间响应曲线下的时间响应曲线(b)相平面曲线相平面曲线 下面，改变下面，改变u u的值，并设仿真终止时间为的值，并设仿真终止时间为30003000，看看计，看看计算结果如何算结果如何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,t

32、2,y2,:)发现，经过很长时间的等待，得出一个错误信息发现，经过很长时间的等待，得出一个错误信息.事实上，由于变步长所采用的步长过小，而要求仿真的事实上，由于变步长所采用的步长过小，而要求仿真的终止时间比较大，导致输出的终止时间比较大，导致输出的y矩阵过大，超出了计算机矩阵过大，超出了计算机存储空间的容限，所以这个问题不适合用存储空间的容限，所以这个问题不适合用ode45来求解来求解,应该采用应该采用ode15s来求解来求解 在许多领域，常会遇到一类特殊的常微分方程，其中在许多领域，常会遇到一类特殊的常微分方程，其中一些解变化缓慢，另一些解变化快，且相差比较悬殊，这一些解变化缓慢，另一些解变

33、化快，且相差比较悬殊，这些方程称为些方程称为刚性方程（刚性方程（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.(a)x1(t)曲线曲线(b)(变化比较平变化比较平滑）滑）(b)x2(t)曲线曲线(c)（在某些点上变化较（在某些点上变化较快）快）其实，许多传统的刚性问题可以采用普通函数来求解，不必刻其实，许多传统的刚性问题可以采用普通函数来求解，不必刻意选择刚性问题的解法，当然，有的却必须用刚性问题的解法意选择刚性问题的解法，当然，有的却必须用刚性问题的解法