第十五章常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象， 并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解， 而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的 方程如d . y2 x2，于是对于用微分方程解决实际问题来说，数值解法就是一个十dx分重要的手段。1常微分方程的离散化dy下面主要讨论一阶常微分方程的初值问题，其一般形式是=f (x, y) a x b(1)y(a) = y在下面的讨论中，我们总假定函数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 )用差商近似导数若用向前差商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)是个离散化的问题，称为差 分方程初值问题。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。(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)以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的 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 = y(a)用这组公式求问题(1)的数值解称为向后 Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有yn ，因此是隐式公式，一般要用迭代法求解，迭代公式通常为(6)ynI =yn+hf(Xn,yn)n= yn +hf(Xn4t, yn*)(k =0,1,2，)2.2 Euler方法的误差估计对于向前Euler公式(3)我们看到，当n =1,2/时公式右端的 y都是近似的， 所以用它计算的yn !会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的 所谓局部截断误差。假定用(3)式时右端的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阶的，而数值算法的精度定义为：若一种算法的局部截断误差为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)的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为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)求解时，每步可以只迭代一次，由此导出一种新的方法一改进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)1yn 卅=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，精度自然很低。改进的 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/ /:为待定系数，看看如何确定它们使(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(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 +十妇 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 , ynkik3当2hk2+寸k41这就是常用的4阶龙格一库塔方法(简称f(Xn h,yn hk3)RK方法)。5线性多步法 以上所介绍的各种数值解法都是单步法， 前一步的值yn ,单步法的一般形式是 丫计=yn +h 申(Xn,