第十五章常微分方程的解法

1、第十五章常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象， 并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解， 而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的 方程如d . y2 x2，于是对于用微分方程解决实际问题来说，数值解法就是一个十dx分重要的手段。§1常微分方程的离散化dy下面主要讨论一阶常微分方程的初值问题，其一般形式是=f (x, y) a <x <b(1)y(a) = y°在下面

2、的讨论中，我们总假定函数f (x, y)连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L ,使得I f(x,y) f(x,y)伍 L |y y |这样，由常微分方程理论知，初值问题的解必定存在唯一。所谓数值解法，就是求问题(1)的解y(x)在若干点a 二 x° ： x ：:. x? ；：- Xn 二 b处的近似值yn( n=1,2, ,N)的方法，y“( n=1,2, ,N)称为问题(1)的数值解, hn二Xn 1 -Xn称为由X.到人1的步长。今后如无特别说明，我们总取步长为常量h。建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：(i )用差商近似导数

3、若用向前差商y(Xn J代替y'(Xn)代入(1)中的微分方程，则得hy(xn 1) y(Xn)一：f (Xn, y(Xn) (n -0,1,)h化简得y(Xn 1)： y(Xn) hf (Xn, y(Xn)如果用 y(xn)的近似值yn代入上式右端，所得结果作为y(xn+)的近似值，记为yn出，则有yn 十二 yn +hf (Xn, yn) (n =0,1，)(2)这样，问题(1)的近似解可通过求解下述问题:丫时=yn +hf(Xnn)(门=0,1，)丿(3 )卜0 = y(a)得到，按式(3)由初值 y可逐次算出y1, y2 /。式(3)是个离散化的问题，称为差 分方程初值问题。需

4、要说明的是，用不同的差商近似导数，将得到不同的计算公式。(ii )用数值积分方法将问题(1)的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端 积分，得Xn +y(Xn 半)y(Xn) = f (X, y(x)dx (n =0,1,)xn右边的积分用矩形公式或梯形公式计算。(iii) Taylor多项式近似将函数y(x)在Xn处展开，取一次 Taylor多项式近似，则得y(Xn 1)： y(Xn) hy'(Xn) *区)hf 区”区)再将y(Xn)的近似值yn代入上式右端，所得结果作为y(Xn.1)的近似值yn.1，得到离散化的计算公式yn 1 二 yn hf (Xn,yn)

5、以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的 Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误£。§欧拉(Euler)方法2.1 Euler 方法Euler方法就是用差分方程初值问题 (3)的解来近似微分方程初值问题 (1)的解, 即由公式(3)依次算出y(xn)的近似值yn(n 1,2/ )。这组公式求问题(1)的数值 解称为向前Euler公式。如果在微分方程离散化时，用向后差商代替导数，即y'(xny(Xn “ 心， 则得计算公式(5)Yn+ = % +hf (Xn*ynG (n =0,1，) M =

6、 y(a)用这组公式求问题(1)的数值解称为向后 Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有yn ，因此是隐式公式，一般要用迭代法求解，迭代公式通常为(6)yn°I =yn+hf(Xn,yn)n= yn +hf(Xn4t, yn*)(k =0,1,2，)2.2 Euler方法的误差估计对于向前Euler公式(3)我们看到，当n =1,2/时公式右端的 y都是近似的， 所以用它计算的yn !会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的 所谓局部截断误差。假定用(3)式

7、时右端的yn没有误差，即yn = y(Xn)，那么由此算出yn 1 = y(Xn) hf (Xn,y(Xn)( 7)局部截断误差指的是，按(7)式计算由Xn到Xn 1这一步的计算值yn 1与精确值y(Xn J (8)之差y(Xn 1)yn 1。为了估计它，由Taylor展开得到的精确值 y(Xn 1)是h23y(Xn .1)= y(Xn) hy'(Xn)y''(Xn) O(h )2(7)、(8)两式相减(注意到 y'=f(x, y)得y(Xn 1)- yn 1(9)h232y''(Xn) O(h ) ： O(h ) 2即局部截断误差是h2阶的，而

8、数值算法的精度定义为：若一种算法的局部截断误差为o(h p")，则称该算法具有p阶精度。显然p越大，方法的精度越高。式(9)说明，向前 Euler方法是一阶方法，因此 它的精度不高。§3 改进的Euler方法3.1 梯形公式(4)中之右端积分,利用数值积分方法将微分方程离散化时，若用梯形公式计算式 即Xn 1Xn 1-hXnf (x,y(x)dx f (Xn, y(Xn) f (Xn 1, y(Xn 1)Xn2并用yn, yn 1代替y(Xn), y(Xn 1)，则得计算公式hYn 1 二 yn 2【f(Xn,yn) f (Xn 1 , Yn 1)这就是求解初值问题(1)的

9、梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为yn? =yn +hf (Xn, yn)(10)ynf yn £f(Xn,yn) f(Xn 1, ynk1)2(k =0,1,2,)L由于函数f (x, y)关于y满足Lipschitz条件，容易看出 I (k 1) 化) hL (k) (k 书| yn 1 -yn 1 | I yn 1 yn 1 |2其中L为Lipschitz常数。因此，当0 ：仏-1时，迭代收敛。但这样做计算量较大。2如果实际计算时精度要求不太高，用公式(10)求解时，每步可以只

10、迭代一次，由此导出一种新的方法一改进Euler法。3.2 改进Euler法按式(5)计算问题(1)的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用： 先用Euler公式求yn 1的一个初步近似值 yn 1，称为预测值，然 后用梯形公式校正求得近似值 yn j，即(11)yn d -yn hf(Xn,yn) 彳hyn 1 =yn 2f (Xn，Yn )f X 1, Vn 1)校正式(11)称为由Euler公式和梯形公式得到的预测一校正系统，也叫改进 为便于编制程序上机，式(11)常改写成yp = yn +hf (Xn, yn)Tq = yn +hf (Xn +h,p)1y

11、n 卅=2(yp +yq)改进Euler法是二阶方法。预测Euler 法。(12)§4 龙格一库塔(RungeKutta)方法回到Euler方法的基本思想一用差商代替导数一上来。实际上，按照微分中值定理 应有y(xn J -y(xn)y'(xn 巾)，0 ： v : 1h 注意到方程y' = f (x, y)就有y(Xnd1)=y(Xn) +hf (Xn +Th,y(Xn +旳)(13)不妨记K = f(Xn7h,y(Xndh),称为区间Xn,Xn“上的平均斜率。可见给出一种斜率K , (13)式就对应地导出一种算法。向前Euler公式简单地取f (xn, yn)为K

12、，精度自然很低。改进的 Euler公式可理 解为K取f (Xn, yn) , f (Xn 1, % 1)的平均值，其中 1二n ' hf(Xn”n)，这种处 理提高了精度。如上分析启示我们，在区间Xn,Xn内多取几个点，将它们的斜率加权平均作为 K，就有可能构造出精度更高的计算公式。这就是龙格一库塔方法的基本思想。4.1首先不妨在区间Xn,Xn 1内仍取2个点，仿照(13)式用以下形式试一下yn 1 =yn hCK丞2)& 二 f (Xn, yn)( 14)k2 = f (Xn +ah, yn + PhkJ 0 ", P其中, -2/ /:为待定系数，看看如何确定它们

13、使(14)式的精度尽量高。为此我们分析局部截断误差y(Xn 1)- yn，1，因为yn =y(Xn)，所以(14 )可以化为(15)yn =y(Xn)+h(人匕 +ek2)匕=f (Xn,y(Xn) = y'(Xn)“2 = f(Xn +ah,y(Xn) + Bhk1)=f (Xn,y(Xn) +ethfx(Xn,y(Xn)+Phk1fy(Xn, y(Xn) +0(h2)其中k2在点(Xn, y(Xn)作了 Taylor展开。(15)式又可表为203yn 1 = y(Xn) ( 1 ，2)hy'(Xn)2： h ( fxffy) O(h )cc注意到h23y(Xn1) = y(

14、Xn) hy'(Xn) y''(Xn) O(h3)2y，可见为使误差y(Xn 1 ) - yn 1 =O(h3)，只须令1 p中 y = f , y'' = fx ff- = 1， '2 ， = 12 a待定系数满足(16)的(15)式称为2阶龙格一库塔公式。由于(1而只有3个方程，所以解不唯一。不难发现，若令、=，2二丄2(16)16)式有4个未知数(15)进的Euler公式。可以证明，在Xn，Xn 1内只取2点的龙格一库塔公式精度最高为2阶。4.2 4阶龙格一库塔公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：yn 十=yn

15、 +十妇 k2 + 嘉k3 + 阪)人=f (Xn,yn)*2 = f (Xn 中h, yn + 怖匕)(17)k3 = f (Xn +讣 yn + fhk1 +p3hk2)札=f (Xn +sh, yn + 04hk1 + 05hk2 + %hk3)其中待定系数'i/'i, 共13个，经过与推导2阶龙格一库塔公式类似、但更复杂的计 算，得到使局部误差y(Xn 1)- yn O(h5)的11个方程。取既满足这些方程、 又较简 单的一组'i/i, -i，可得(18)hy 12k2 2k3 kq)6=f (Xnn)h=f(Xn g,ynr h-f (xn 2 , ynkik

16、3当2hk2+寸k41这就是常用的4阶龙格一库塔方法(简称f(Xn h,yn hk3)RK方法)。§5线性多步法 以上所介绍的各种数值解法都是单步法， 前一步的值yn ,单步法的一般形式是 丫计=yn +h 申(Xn, yn,h)这是因为它们在计算yn 1时，都只用到(n =0,1/ ,N -1)(19)其中(x, y,h)称为增量函数，例如 Euler方法的增量函数为 f (x, y)，改进Euler法的 增量函数为1(x,y,h)f(x,y) f(x h, y hf(x, y)2如何通过较多地利用前面的已知信息，如yn, yn,，yn_r，来构造高精度的算法计算yn 1，这就是多

17、步法的基本思想。经常使用的是线性多步法。让我们试验一下r -1，即利用yn,ynv计算yn .1的情况。从用数值积分方法离散化方程的(4)式xny(Xn 1) -y(Xn)二.f (x, y(x)dxXn出发'记 f (Xn, yn)二 fn， f (X(Xn j,yn , yn d , ,yn_r ，n V, yn J )二 ffn 4)，(Xn, fn)的插值公式得到(因 X 一 X.)，所以是外插。X诂n -4XnJ - Xn式中被积函数f (x, y(x)用二节点x _Xn/f (x,y(x) = fn f Xn -Xnj1(X - Xn j) fn h此式在区间Xn, Xn上

18、积分可得Xn 1x f(x, y(x)dxXn于是得到(20)一(X - Xn) fn 4.hyn 1 _ yn (3fnifnV)2注意到插值公式(20)的误差项含因子(X -Xn 4)(X-Xn )，在区间Xn , Xn上积分后 得出h3，故公式(21)的局部截断误差为 O(h3)，精度比向前Euler公式提高1阶。(21)若取r =2,3,可以用类似的方法推导公式，如对于r = 3有hyn 1 二 yn 24(55fn -59fn-37 f -9f) 其局部截断误差为 0(h5)。f (x, y(x)用的插值公式由外插改为内插，可进yn 1, yn,，yn_r 1，取=1时得到的是梯形公

19、式，取=3时步减如果将上面代替被积函数小误差。内插法用的是可得(9 fn 1 19 fn - 5 f 24与(22)式相比，虽然其局部截断误差仍为 为减小，误差还是小了。当然， 样，用迭代或校正的办法处理。yn 1 = ynfn/)(23)50(h)，但因它的各项系数(绝对值)大(23)式右端的fn d未知，需要如同向后 Euler公式一§6 一阶微分方程组与高阶微分方程的数值解法6.1 一阶微分方程组的数值解法 设有一阶微分方程组的初值问题Vi = fi(x,浙2,ym)y(a) = yi0(i =12,m)(24)f =：( fi, f2，fm)T，则初值若记 y =(y1, y

20、2, ym)T，y。=(%。，y?。,ymo)T，问题(24)可写成如下向量形式y'=f(x, y)(a) = yo如果向量函数 f(x, y)在区域D :a _ x _ b, y Rm连续，且关于y满足Lipschitz条件，即存在L 0，使得对xa,b，y1, y Rm， 都有(25)I f(x, yj - f (x,y2)乞 L y1 -y2那么问题(25)在a,b上存在唯一解y = y(x)。问题(25)与(1)形式上完全相同，故对初值问题(1)所建立的各种数值解法可 全部用于求解问题(25)。6.2高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组

21、初值问题。 设有m阶常微分方程初值问题；y(m)=f (x,y,y',，y/ , / X(1)(mJ) / (,(a) =y°,y (a) = y° , y (a) = y°引入新变量y1 = y, y2二y',ym二y(m ，问题(26)就化为一阶微分方程初值问题(2) a 乞xb(mV).(m 4),y (a)二 yo(26)(m -4)Vi=y2yi(a) = y°(1)y； = y3y2(a) = y°=(27)ym= ymym(a)=y°y'm = f (x, yi, ym)ym(a) = y0mJ)

22、然后用6.1中的数值方法求解问题(27),就可以得到问题(26)的数值解。最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组dv Ay-T(x)( 28)dx其中 y Rm,A为m阶方阵。若矩阵 A的特征值i(i =1,2- ,m)满足关系Re i : 0 (i =1,2,m)max| Re min | Re i |1 岂卯1 im则称方程组(28)为刚性方程组或 Stiff方程组，称数s = max| Re i |/min | Re i |恒哲| | 1生g |i |为刚性比。对刚性方程组，用前面所介绍的方法求

23、解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长h不作任何限制。§ Matlab 解法7.1 Matlab数值解7.1.1非刚性常微分方程的解法Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45, ode23,ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比 ode45高。Matlab的工具箱中没有 Euler方法的功能函数。(I)对简单的一阶方程的初值问题；y = f(x,y)MX。)=

24、 y。改进的Euler公式为yp =yn +hf(Xn,yn)Ym +hf (Xn +h, Yp)1yn 十=2(Yp +yq)我们自己编写改进的 Euler方法函数eulerpro.m如下：functionx,y=eulerpro (fun, x0,xfi nal,y 0,n);if nargin<5,n=50;e ndh=(xfi nal-x0)/n;x(1)=xO;y(1)=yO;for i=1: nx(i+1)=x(i)+h;y1= y(i)+h*feval(fu n, x(i),y(i); y2=y(i)+h*feval(fu n,x(i+1),y1); y(i+1)=(y1+

25、y2)/2;end例1用改进的Euler方法求解y，- -2y 2x2 2x , (0 込 x 込0.5) , y(0) = 1 解编写函数文件doty.m如下：fun ctio n f=doty(x,y); f=-2*y+2*xA2+2*x;在Matlab命令窗口输入：x,y=eulerpro('doty',0,0.5,1,10)即可求得数值解。(II)ode23，ode45, ode113的使用Matlab的函数形式是t,y=solver('F',tspan， y0)这里solver为ode45，ode23, ode113，输入参数F是用M文件定义的微分方程

26、 y' = f(x, y)右端的函数。tspan=t0,tfinal是求解区间，y0是初值。例2用R防法求解2y' =2y 2x 2x，(0 乞 x 0.5)，y(0) =1解 同样地编写函数文件 doty.m如下：fun ctio n f=doty(x,y); f=-2*y+2*xA2+2*x;在Matlab命令窗口输入：x,y=ode45('doty',0,0.5,1)即可求得数值解。7.1.2刚性常微分方程的解法Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s,ode23t，ode23tb，这些函数的使用同上述非刚性微

27、分方程的功能函数。7.1.3 高阶微分方程 y(n f(t, y,y'/ ,y(n4)解法例3考虑初值问题ym_3y“_y'y=0丫(0)=0 W(0)=1 厂(0)=-1y1(0) =0y2(0) =1y3(0) 一 1解(i)如果设y1 y, y2y',y"，那么y2 二 y3.y'3 =3y3 y2 y1Y1 y2初值问题可以写成 Y'=F(t,Y),Y(0) =Y0的形式，其中y 二y1；y2；y3】。(ii)把一阶方程组写成接受两个参数t和y，返回一个列向量的M文件F.m：fun ctio ndy=F(t,y);dy=y(2);y(3

28、);3*y(3)+y(2)*y(1);注意：尽管不一定用到参数 t和y , M文件必须接受此两参数。这里向量dy必须是列向量。(iii )用Matlab解决此问题的函数形式为T,Y=solver('F',tspa n, 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)是解的导数

29、，Y(:,3)是解的二阶导数。例4 求van der Pol 方程 y''_P(1 一 y2)y'+y =0的数值解，这里.0是一参数。解(i )化成常微分方程组。设y1 = y, y2 = y'，则有y' y2丿 2 畀2 =卩(1yy2 - y1(ii)书写 M文件(对于"-1)vdp1.m:function dy=vdp1(t,y); dy=y(2);(1-y(1)A2)*y(2)-y(1);(iii )调用Matlab函数。对于初值y(0) = 2, y'(0) =0，解为T,Y=ode45('vdp1' ,0

30、20,2;0);(iv )观察结果。利用图形输出解的结果：plot(T,Y(:,1),'-',T,Y(:,2),'-')title('Soluti on of van der Pol Equati on,m u=1');xlabel( 'time t' );ylabel( 'solutio n y');legend('y1', 'y2');-2 - '：Soluti on of van der Pol Equati on ,mu=1Wy2y npuosVVV3 IIIIIIi

31、ll02468101214161820time t例5 van der Pol 方程，=1000 (刚性) 解(i) 书写M文件vdp1000.m:fun ction dy=vdp1000(t,y);dy=y(2);1000*(1-y(1)A2)*y(2)-y(1);(ii )观察结果t,y=ode15s( 'vdplOOO' ,0 3000,2;0);plot(t,y(:,1),'o')title('Solutio n of van der Pol Equatio n,mu=1000');xlabel( 'time t' );y

32、label( 'solution y(:,1)');7.2常微分方程的解析解在Matlab中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令 dsolve。常微分方程在Matlab中按如下规定重新表达：符号D表示对变量的求导。Dy表示对变量y求一阶导数，当需要求变量的n阶导数时，用Dn表示，D4y表示对变量y求4阶导数。由此，常微分方程y'' 2yy在Matlab中，将写成'D2y+2*Dy=y'。7.2.1求解常微分方程的通解无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：dsolve('diff_equati

33、on') dsolve(' diff_equation'， 'var')式中diff_equation为待解的常微分方程，第1种格式将以变量t为自变量进行求解， 第2种格式则需定义自变量 var。例6试解常微分方程x2 y (x -2y)y' = 0解编写程序如下：syms x ydiff_equ='xA2+y+(x-2*y)*Dy=0'dsolve(diff_equ,'x')7.2.2 -求解常微分方程的初边值问题求解带有初边值条件的常微分方程的使用格式为：dsolve('diff_equation

34、9; , 'condition1 , condition2 ,''var') 其中condition1 , condition2 ,即为微分方程的初边值条件。例7试求微分方程W=x , y(1)=8, y'(1)=7,y'' (2) =4 的解。解编写程序如下：y=dsolve( 'D3y-D2y=x', 'y(1)=8,Dy(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试求常微分方程组：”f''+3g