第十五章_常微分方程解法.doc

1、-179-第十五章 常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 ，于是对于用微分方程解决实际问题来说，数值解法就2xyd是一个十分重要的手段。1 常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是(1)0)(,yabxaxfd在下面的讨论中，我们总假定函数 连续，且关于 满足李普希兹(Lipschit

2、z)条),(fy件，即存在常数 ，使得L|),(| yLxf 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。所谓数值解法，就是求问题（1）的解 在若干点)(xbaN210处的近似值 的方法， 称为问题（1）的数值解，),2(Nny ),ny称为由 到 的步长。今后如无特别说明，我们总取步长为常量 。nxh1nx h建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：（i）用差商近似导数若用向前差商 代替 代入（1）中的微分方程，则得hynn)(1)(nx),0,) fxy化简得 )(,)(1nnnxyf如果用 的近似值 代入上式右端，所得结果作为 的近似值，记为 ，)xyy

3、 )(1nxy1ny则有),0(),(1 hfnn（2）这样，问题（1）的近似解可通过求解下述问题)( ),1(),01ayyxfnn -180-（3）得到，按式（3）由初值 可逐次算出 。式（3）是个离散化的问题，称为0y,21y差分方程初值问题。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。（ii）用数值积分方法将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得(4) 1 ),10(),()(1nxn ndxyfyx 右边的积分用矩形公式或梯形公式计算。（iii ）Taylor 多项式近似将函数 在 处展开，取一次 Taylor 多项式近似，则得)

4、(n )(,)()(1 nnnxyhfxyhxy再将 的近似值 代入上式右端，所得结果作为 的近似值 ，得到离nxy 11ny散化的计算公式),(1nnf以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的 Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。2 欧拉（Euler）方法2.1 Euler 方法Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出 的近似值 。这组公式求问题（1）)(nxy),2(ny的数值解称为向前 Euler 公式。如果在微分方程离散化时，用向后差商代替

5、导数，即，则得计算公式hxynnn)()( 11)( ),10(),011aynyxhfn（5）用这组公式求问题（1）的数值解称为向后 Euler 公式。向后 Euler 法与 Euler 法形式上相似，但实际计算时却复杂得多。向前 Euler 公式是显式的，可直接求解。向后 Euler 公式的右端含有 ，因此是隐式公式，一般要1ny用迭代法求解，迭代公式通常为 ),20(),(1)1(0 kxhfyknnkn（6）2.2 Euler 方法的误差估计对于向前 Euler 公式（3）我们看到，当 时公式右端的 都是近似的，,ny-181-所以用它计算的 会有累积误差，分析累积误差比较复杂，这里先

6、讨论比较简单的1ny所谓局部截断误差。假定用（3）式时右端的 没有误差，即 ，那么由此算出ny)(nxy,()(1nnxhf（7）局部截断误差指的是，按（7）式计算由 到 这一步的计算值 与精确值11ny之差 。为了估计它，由 Taylor 展开得到的精确值 是)(1nxy1)(nyx )(x)(2)()( 3hOxyhxynn（8）（7） 、 （8）两式相减（注意到 ）得,f)()()( 2321xyxynn （9）即局部截断误差是 阶的，而数值算法的精度定义为：2h若一种算法的局部截断误差为 ，则称该算法具有 阶精度。)(1phOp显然 越大，方法的精度越高。式（9）说明，向前 Euler

7、 方法是一阶方法，因此p它的精度不高。3 改进的 Euler 方法3.1 梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分，即 )(,)(,2)(, 11 nnx xyfxyfhdxyfn并用 代替 ，则得计算公式,ny1n),(),(11 nn ff这就是求解初值问题（1）的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 (10) ),210( ),(),(1)()0(1kyxfyxfhyknnnn由于函数 关于 满足 Lipschitz 条件，容易看出),xf | )1

8、()()(1)( knknk yhLy-182-其中 为 Lipschitz 常数。因此，当 时，迭代收敛。但这样做计算量较大。L120hL如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法改进 Euler 法。3.2 改进 Euler 法按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将 Euler 公式与梯形公式结合使用：先用 Euler 公式求 的一个初步近似值 ，称为预测值，1ny1ny然后用梯形公式校正求得近似值 ，即1n 校 正预 测 ),(),(211 nnn yxfyxfhy（11）式（11）称为由 Euler 公式和梯形

9、公式得到的预测校正系统，也叫改进 Euler 法。为便于编制程序上机，式（11）常改写成)(21,)(qpnpnqnpyyhxf（12）改进 Euler 法是二阶方法。4 龙格库塔（RungeKutta）方法回到 Euler 方法的基本思想用差商代替导数上来。实际上，按照微分中值定理应有10),()(1 hxyhxynnn注意到方程 就有),f)(,)(1fnnnn（13）不妨记 ，称为区间 上的平均斜率。可见给出一,hxyfK,1x种斜率 ， （13）式就对应地导出一种算法。向前 Euler 公式简单地取 为 ，精度自然很低。改进的 Euler 公式可理),(nyfK解为 取 ， 的平均值，

10、其中 ，这种处),(nxf 1n ),(1nnyxhfy理提高了精度。如上分析启示我们，在区间 内多取几个点，将它们的斜率加权平均作为,nx，就有可能构造出精度更高的计算公式。这就是龙格 库塔方法的基本思想。K4.1 首先不妨在区间 内仍取 2 个点，仿照（13）式用以下形式试一下,1n-183- 1,0),()12121hkyxfkynn（14）其中 为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们,21分析局部截断误差 ，因为 ，所以（14）可以化为1)(ny)(nx )()(, ,(,)(21121 2hOxyfhkfxkfhnnnnn（15）其中 在点 作了 Taylor

11、 展开。 （15）式又可表为2k)(,nyx )()( 32211 hfxyyxnn 注意到 )()()( 321 Ohxynnnn中 ， ，可见为使误差 ，只须令f yf )1hyx， ， 2121（16）待定系数满足（16）的（15）式称为 2 阶龙格库塔公式。由于（16）式有 4 个未知数而只有 3 个方程，所以解不唯一。不难发现，若令 ， ，即2111为改进的 Euler 公式。可以证明，在 内只取 2 点的龙格库塔公式精度最高,1nx为 2 阶。4.2 4 阶龙格库塔公式要进一步提高精度，必须取更多的点，如取 4 点构造如下形式的公式： ),(),)( )36251434 32231

12、21 32hkhkyxfkf kkhynnnn （17）其中待定系数 共 13 个，经过与推导 2 阶龙格库塔公式类似、但更复杂的ii,-184-计算，得到使局部误差 的 11 个方程。取既满足这些方程、又)()(51hOyxn较简单的一组 ，可得ii,),(2),()2(6343121 431hkyxfkhkyxfkynnn（18）这就是常用的 4 阶龙格库塔方法（简称 RK 方法） 。5 线性多步法以上所介绍的各种数值解法都是单步法，这是因为它们在计算 时，都只用到1ny前一步的值 ，单步法的一般形式是ny),10(),(1 Nnhyxnn （19）其中 称为增量函数，例如 Euler 方

13、法的增量函数为 ，改进 Euler 法),(hx ,(yxf的增量函数为 ),(),(2, hfyxfyfy如何通过较多地利用前面的已知信息，如 ，来构造高精度的算法rnn1计算 ，这就是多步法的基本思想。经常使用的是线性多步法。1ny让我们试验一下 ，即利用 计算 的情况。1r1,ny从用数值积分方法离散化方程的（4）式1)(,)(1nxn dxfyx出发，记 ， ，式中被积函数 用二节点ff, 1nfy)(,xyf， 的插值公式得到（因 ，所以是外插。),(1nn )()(1 )(, 111nnnnfxfxhfxyf（20）此式在区间 上积分可得,1nx-185-123)(,1 nx fh

14、dxyfn于是得到)(11nnf（21）注意到插值公式（20）的误差项含因子 ，在区间 上积分后)(1nnxx,1nx得出 ，故公式（21）的局部截断误差为 ，精度比向前 Euler 公式提高 1 阶。3h3hO若取 可以用类似的方法推导公式，如对于 有,2r r)9759(43211 nnn fffhy（22）其局部截断误差为 。)5O如果将上面代替被积函数 用的插值公式由外插改为内插，可进一步减)(,xyf小误差。内插法用的是 ，取 时得到的是梯形公式，取 时11rnn 3r可得)59(242111 nnnn fffhy（23）与（22）式相比，虽然其局部截断误差仍为 ，但因它的各项系数（

15、绝对值）大(hO为减小，误差还是小了。当然， （23）式右端的 未知，需要如同向后 Euler 公式一1nf样，用迭代或校正的办法处理。6 一阶微分方程组与高阶微分方程的数值解法6.1 一阶微分方程组的数值解法设有一阶微分方程组的初值问题021)(),ii myayxf ),21(mi（24）若记 ， ， ，则初值Tmy,(21T),(0210 Tmff),(21问题（24）可写成如下向量形式0)(,yaxf（25）如果向量函数 在区域 ：,(yxfDmRb,连续，且关于 满足 Lipschitz 条件，即存在 ，使得对 ，0L,bax，都有mRy21,-186-2121),(),(yLyxf

16、f 那么问题（25）在 上存在唯一解 。ba)(x问题（25）与（1）形式上完全相同，故对初值问题（1）所建立的各种数值解法可全部用于求解问题（25） 。6.2 高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有 阶常微分方程初值问题m)1(0)1()1(00)( , mmm yayay bxxf（26）引入新变量 ，问题（26）就化为一阶微分方程初值问)(21, m题 )1(01 211 )1(0232 1 ),( mmmyayxfyyya （27）然后用 6.1 中的数值方法求解问题（27） ，就可以得到问题（26）的数值解。最后需要指出的是，在化

17、学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组)(xAydx（28）其中 为 阶方阵。若矩阵 的特征值 满足关系Rym,),21(mi),21(0emii |Re|n|ax1 iiii 则称方程组（28）为刚性方程组或 Stiff 方程组，称数 |/|1imiimis为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长 不作任何限制。h7 Matlab 解法7.1 Matlab 数值解7.1.1 非刚性常微分方程的解

18、法Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113 ，其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的-187-首选方法，ode23 采用二三阶 RK 方法，ode113 采用的是多步法，效率一般比 ode45高。Matlab 的工具箱中没有 Euler 方法的功能函数。（I）对简单的一阶方程的初值问题0)(,yxf改进的 Euler 公式为)(21,)(qpnpnqnpyyhxf我们自己编写改进的 Euler 方法函数 eulerpro.m 如下：function x,y=eulerpro(fun,x0,xfinal,y0,

19、n);if nargin5,n=50;endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end例1 用改进的Euler方法求解， ，xy2)5.0(1(y解 编写函数文件 doty.m 如下：function f=doty(x,y);f=-2*y+2*x2+2*x;在Matlab命令窗口输入：x,y=eulerpro(doty,0,0.5,1,10)即可求得数值解。（II）

20、ode23，ode45，ode113的使用Matlab的函数形式是t,y=solver(F,tspan，y0)这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程右端的函数。tspan=t0,tfinal是求解区间，y0是初值。),(yxf例2 用RK方法求解， ，x2)5.0(1(y解 同样地编写函数文件 doty.m 如下：function f=doty(x,y);f=-2*y+2*x2+2*x;在Matlab命令窗口输入：x,y=ode45(doty,0,0.5,1)-188-即可求得数值解。7.1.2 刚性常微分方程的解法Matlab的工具箱提

21、供了几个解刚性常微分方程的功能函数，如ode15s， ode23s，ode23t ，ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。7.1.3 高阶微分方程 解法),(1()( nnytfy例 3 考虑初值问题10(0 y解 （i）如果设 ，那么,3211)0( 31232yy初值问题可以写成 的形式，其中)(,YtFY。;321yY（ii）把一阶方程组写成接受两个参数 和 ，返回一个列向量的 M 文件 F.m：tfunction dy=F(t,y);dy=y(2);y(3);3*y(3)+y(2)*y(1);注意：尽管不一定用到参数 和 ，M 文件必须接受此两参数。这里向量 必

22、须是ty dy列向量。（iii）用 Matlab 解决此问题的函数形式为T,Y=solver(F,tspan,y0)这里 solver 为 ode45、ode23、ode113，输入参数 F 是用 M 文件定义的常微分方程组，tspan=t0 tfinal是求解区间，y0 是初值列向量。在 Matlab 命令窗口输入T,Y=ode45(F,0 1,0;1;-1)就得到上述常微分方程的数值解。这里 Y 和时刻 T 是一一对应的，Y(:,1)是初值问题的解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。例 4 求 van der Pol 方程0)1(2yy的数值解，这里 是一参数。解 （i）

23、化成常微分方程组。设 ，则有,21y221)(yy（ii）书写 M 文件（对于 ）vdp1.m:function dy=vdp1(t,y);dy=y(2);(1-y(1)2)*y(2)-y(1);（iii）调用 Matlab 函数。对于初值 ，解为0)(,2)0(yT,Y=ode45(vdp1,0 20,2;0);（iv ）观察结果。利用图形输出解的结果：plot(T,Y(:,1),-,T,Y(:,2),-)title(Solution of van der Pol Equation,mu=1);xlabel(time t);-189-ylabel(solution y);legend(y1,

24、y2);0 0123 1 例 5 van der Pol 方程， （刚性）10解 (i)书写 M 文件 vdp1000.m:function dy=vdp1000(t,y);dy=y(2);1000*(1-y(1)2)*y(2)-y(1);（ii）观察结果t,y=ode15s(vdp1000,0 3000,2;0);plot(t,y(:,1),o)title(Solution of van der Pol Equation,mu=1000);xlabel(time t);ylabel(solution y(:,1);7.2 常微分方程的解析解在 Matlab 中，符号运算工具箱提供了功能强大的

25、求解常微分方程的符号运算命令dsolve。常微分方程在 Matlab 中按如下规定重新表达：符号 D 表示对变量的求导。 Dy 表示对变量 y 求一阶导数，当需要求变量的 n 阶导数时，用 Dn 表示，D4y 表示对变量 y 求 4 阶导数。由此，常微分方程 在 Matlab 中，将写成D2y+2*Dy=y。y27.2.1 求解常微分方程的通解无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：dsolve(diff_equation)dsolve( diff_equation，var)式中 diff_equation 为待解的常微分方程，第 1 种格式将以变量 t 为自变量进行求解，

26、第 2 种格式则需定义自变量 var。例 6 试解常微分方程0)2(yxy解 编写程序如下：syms x ydiff_equ=x2+y+(x-2*y)*Dy=0;dsolve(diff_equ,x)-190-7.2.2 求解常微分方程的初边值问题求解带有初边值条件的常微分方程的使用格式为：dsolve(diff_equation，condition1，condition2，var)其中 condition1，condition2 ， 即为微分方程的初边值条件。例 7 试求微分方程，xy 4)2(,7)1(,8)(y的解。解 编写程序如下：y=dsolve(D3y-D2y=x,y(1)=8,Dy

27、(1)=7,D2y(2)=4,x)7.2.3 求解常微分方程组求解常微分方程组的命令格式为：dsolve(diff_equ1，diff_equ2，var)dsolve(diff_equ1，diff_equ2，condition1，condition2，var)第 1 种格式用于求解方程组的通解，第 2 种格式可以加上初边值条件，用于具体求解。例 8 试求常微分方程组： xfgcosin3的通解和在初边值条件为 的解。1)5(,3)(,0)2(gf解 编写程序如下：clc,clearequ1=D2f+3*g=sin(x);equ2=Dg+Df=cos(x);general_f,general_g

28、=dsolve(equ1,equ2,x)f,g=dsolve(equ1,equ2,Df(2)=0,f(3)=3,g(5)=1,x)7.2.4 求解线性常微分方程组（i）一阶齐次线性微分方程组， ，AX nx1nnaA 11这里的表示对 求导数。 是它的基解矩阵。 ， 的解为tteAX 0)(t。0)()ettA例 9 试解初值问题，XX2013 12)0(解 编写程序如下：syms ta=2,1,3;0,2,-1;0,0,2;x0=1;2;1;-191-x=expm(a*t)*x0（ii）非齐次线性方程组由参数变易法可求得初值问题，)(tfAX0Xt的解为.tsAt dfeet00 )()()

29、(例 10 试解初值问题， 。teXXt2cos123 10)(X解 编写程序如下：clc,clearsyms t sa=1,0,0;2,1,-2;3,2,1;ft=0;0;exp(t)*cos(2*t);x0=0;1;1;x=expm(a*t)*x0+int(expm(a*(t-s)*subs(ft,s),s,0,t);x=simple(x)习 题 十 五1. 用欧拉方法和龙格库塔方法求微分方程数值解，画出解的图形，对结果进行分析比较。（i） （Bessel 方程，令2,2,0)(22 yynxyx，精确解 。)1nsi（ii ） ，幂级数解)(,1)(,co yxy8642!59!1x2. 一只小船渡过宽为 的河流，目标是起点 正对着的另一岸 点。已知河水dAB流速 与船在静水中的速度 之比为 。1v2vk（i）建立小船航线的方程，求其解析解。（ii ）设 m， m/s， m/s，用数值解法求渡河所需时间、012v任意时刻小船的位置及航行曲线，作图，并与解析解比较。