常微分方程数值解法教材

1、常微分方程数值解法科学技术中的许多问题，在数学中往往归结为微分方程的求解问题。例如天文学中研究星体运动,空间技术中研究物体飞行等，都需要求解常微分方程初值问题。最简单的常微分方程初值问题是包=f (x, y), a 乞 x 乞 b,dx(*)y(xg) =y.定理如果f(x,y)在区域D二(x,y) a乞x乞b, y :内连续，且关于y满足Lipschitz条件，那么初值问题 (*)的解存在且惟一。因为数值解是求微分方程解y(x)的近似值，所以总假定微分方程的解存在惟一，即初值问题(*)中的f(x,y)满足定理的条件。除特殊情形外，初值问题(*) 一般求不出解析解，即使有的能求出解析解，其函数

2、表示式也比较复杂，计算量比较大，况且实际问题往往只要求在某一时刻解的函数值，因此，有必要研究初值问题(*)的数值解法。所谓初值问题(*)的数值解法，就是寻求解y(x)在区间a,b上的一系列点% ： x2 ：: X3 ：:xn =： 上的近似值,y2,，yn,.记h严x x-(i =1,2, |)表示相邻两个节点的间距，称为步长。求微分方程数值解的主要问题：(1) 如何将微分方程 史=f (x, y)离散化，并建立求其数值解的递推公式；dx(2) 递推公式的局部截断误差、数值数yn与精确解y(Xn)的误差估计；(3) 递推公式的稳定性与收敛性.、Euler方法Euler法是最早的解决一阶常微分方

3、程初值问题的一种数值方法，虽然它不够精确，很少被采 用，但是它在某种程度上反映了数值方法的基本思想和特征。考虑初值问题dy = f(x, y),(1.1).y(Xo) = y。，(1.2)为了求得它在等距离散点X, ：X2 ”：以Xn ”： 上的数值解，首先将(1.1)离散化。设h =Xj -Xj丄(i =1,2,1 H)，将式(1.1)离散化的办法有 Taylor展开法、数值微分法及数值积分法。如果在点xn将y(xn J = y(xn h)作Taylor展开，得 h2“ m(1.3)y(Xn 1) = y(Xn) hy (Xn) 2! y( n)， n(Xn, Xn 1)那么当h充分小时，略

4、去误差项一 y ( n)，用yn近似替代y(xn)、yn .1近似替代y(Xn d)，2h2并注意到y(Xn)二f(Xn，y(Xn)，便得(1.4)yo =y(x),YnYn hf (Xn,yn),n =0,1,|, N -1b a其中Xn =X0nh, h.用(1.4)求解(1.1)的方法称为Euler方法。N如果利用差商近似替代微商，那么可得y(Xn 1) - y(Xn)一一 ：T(Xn) = f(Xn,y(Xn).(1.5)h在(1.5)中若用yn近似替代y(xn)、yn d近似替代y(xn.j，同样得到递推公式(1.4). 如果在xn,xn上对y = f (x, y(x)积分，得Xn

5、+y(Xn+) y(Xn)= k f (x,y(x)dx.(1.6)那么对上式右端积分用左矩形求积公式，并用yn近似替代y(xn)、yn .1近似替代y(xn4)，也可得到递推公式(1.4).我们知道，在xOy平面上，微分方程 3 =f(x, y)的解称为积分曲线，积分曲线上一点dx(x, y)的切线斜率等于函数f(x, y)的值。如果在 D中每一点，都画上一条以f (x, y)在这点的值为斜率并指向x增加方向的有向线段，即在 D上作出了一个由方程 史二f(x,y)确定的 dx方向场，那么方程的解f = y(x).从几何上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向相

6、一致。从初始点P0(x0,y0)出发，过这点的积分曲线为y = y(x)，斜率为 f(X0,y).设在 x = x。附近y(x)可用过F0点的切线近似表示,切线方程为 y = y0 f (x0,y0)(x0).当 x = & 时,y(x1)的近似值为y0 f (冷，丫0)(为- x0)，并记为，这就是得到似公式m f(%,%)&-).当 X=X2时，y(X2) 的近似值为+ f (x1, y1)(x2 x1)，并记为y2.于是就得到当x=x2时计 算y(x2)的近似公式讨2=5 f (X1,yJ(X2 -xj.重复上面方法，一般可得当X二xn V的计算y(xn,)的近似公式yn 勺二 yn f

7、 (Xn,%)(Xn .i Xn).如果h詁-紀(i =1,2j|)，则上面公式就是(i.4).将P,Pi,，Pn连续起来，就得到一条折线，所以Euler方法又称为折线法。由公式(1.4)看出，已知y0便可算出 力.已知y1，便可算出y2，如此继续下去，这种只 用前一步的值yk便可计算出yk1的递推公式称为 单步法。若在(1.6)中，其右边的积分由数值积分的右矩形公式近似，并用yn替代y(xn)， yn d替代y(xn,)，则可得到yy(xo),、(1.7)Jn+ =yhf(Xn+, yn + ),并称(1.7)为后退Euler公式。Euler公式(1.4)是关于yn彳的一个直接计算公式，是显

8、式的；而公式(1.7)右端是含有ynr的一个函数方程，因此是隐式的，也是单步法。45图9-1若在公式(1.6)中,其右边积分用数值梯形积分公式近似，并用yn替代y(xn)， ynd替代f (xn 1)，则可得到梯形方法公式hYn Yn-f (xn.Yn) f &n 1, Yn 计几(1.8)2梯形方法同后退Euler方法一样都是隐式单步法.对于隐式方法，通常采用迭代法。对后退Euler方法，先Euler方法计算y“十，并将它作为初值丫补？，即= yn + hf (x“, y“),再把它代入(1.7)的右端，便得到后退Euler方法的迭代公式为(1.9)ynr = yn +hf(Xn,yn),%

9、 +hf(XnHt,yn*),k=O,1；同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为yn01 *n hf (Xn,yn),丫匕1yn fgyn) f(Xn1,ynk1),k = 0,1,、误差分析首先讨论Euler方法的误差不妨假定在计算yn 时用到前面一步的值是准确值 y(xn),即 yn =y(Xn).利用算法由y(Xn)计算出的y(Xn 1)的近似值为yn.1，则丫(人.1)-1称为局部 截断误差，也就是计算一步所产生的误差，记为Tn1.将 y(xn +)在xn处作Taylor展开Xn ：:Xn1 -y(Xn iH y(Xn h)二 y(Xn) hf(X. , y(Xn

10、)由Euler方法(1.4)得n 1 二 yn hf (人小)=y(Xn ) hf (人，y(Xn ).上面两式相减得Tn 1 二 丫区 i) n 1 二与 丫 ( )(1.11)若y(X)在解的存在区间a,b上充分光滑，且记 M2=maxy(x),则Euler方法局部截X哥a,b J断误差为h2Tn1 EM2T2=O(h2)(1.12)由于计算yn1时(n=o除外)，一般用到的不是y(Xn)，而是近似值yn,因为每一步计算 除局部截断误差外，还应考虑前一步不准确而引起的误差。我们称这种误差为总体截断误差或方法误差。体截断误差有一个积累过程， 它与计算步数有关。 把第n步的总体截断误差记为 e

11、n二y(Xn) - yn.对任意的第n - 1步有en+ = y(Xn+)ynM|y(XnG n+ M 一 丫4由11)知丫(召卅)一估=人十|，故只须估计出yn1 - yn 1, = y(Xn) hf (x., y(Xn) - y. -hf (x., y.)兰 |y(Xn) yn| +h f (Xn, 丫区)- f(X., Yn).因为f(x, y)关于y满足Lip条件，所以战十yn 兰 en +hL 編=(1 + hL) en这里L是Lip常数。综上所述lej 兰 Tn十(1 + hL)q(1.13)(1.佝给出了第n 1步总体截断误差与第 n步总体截断误差之间的关系，它对一对n都成立。特

12、别当 n=0 时,ej 兰Tj+(1+hL)e0，由(1.2)知：eo=y(x0) y0= 0，因此en ZTn (hL)en二|Tn +(1+hL)|Tn+(1 + hL)|en二 Tn (1 hL)Tnd (1 hL)2 en, _|l| 0时，en 0，也就是说当h充分小时近似值yn能和y(xj充分接近，即数值解是收敛的。对后退Euler方法，类似上面的讨论可知，其局部截断误差为y(Xn 1) - yn 厂hj23y (Xn) O(h ),即局部截断误差关于 h是二阶的，整体截断误差关于h是一阶的.梯形方法公式是将(1.6)中右边积分用梯形求积公式得到的，而梯形方法公式的截断误h 3差为

13、-yn), Xn ：n ”： X. 1，因此梯形方法公式(1.8)的局部截断误差关于h是三阶的,12总体截断误差关于h是二阶的。可以证明：对后退 Euler方法的迭代公式(1.9),如果h充分小，使得0 ： hL :：: 1，那么迭代 过程是收敛的；对梯形方法的迭代公式(1.10)，如果h充分小，使得0：匹1，那么迭代过2程也是收敛的(这里L是Lip常数)。三、改进的Euler方法梯形方法的迭代公式(1.10)比Euler方法精度高，但其计算较复杂，在应用公式(1.10)进行计 算时，每迭代一次，都要重新计算函数 f(x, y)的值,且还要判断何时可以终止或转下一步计算.为了控制计算量和简化计

14、算法，通常只迭代一次就转入下一步计算.具体地说，我们先用Euler公 式求得一个初步的近似值 1，称之为预测值，然后用公式(1.10)作一次迭代得yn-1,即将0n1校 正一次.这样建立的预测一校正方法称为 改进的Euler方法：预测：ynyn hf(Xn，yn),校正：yn 1 二 yn h f(Xn，yn) f (Xn 1，Vn 1).(1.15)这个计算公式也可以表示为yp 二 ynhf (Xn,yn), yc = yn hf (Xn 1, yp),1% 卅=(yp +yJ2取步长h =0.1，分别用Euler方法及改进的Euler方法求解初值问题0乞 x1.y(0) -1.解 这个初值

15、问题的准确解为y(x) =1：(2ex -x-1).根据题设知f (x, y)二-y(1 xy).(1) Euler方法的计算式为yn 1 二 yn -0.1 yn(1 Xn%), 由 y0 = y(0)二1,得y =1 - 0.1 1 (1 0 1) =0.9,y2 =0.9-0.1 0.9 (1 0.1 0.9) =0.8019,这样继续计算下去，其结果列于表 9.1.改进的Euler方法的计算式为yp = yn -0.1xyn(1 + Xnyn),yc= yn -.1 yp(1 召 卩),1yn 十=3(yp +yc),2由 yo = y(0) - 1，得yp =1 0.1叫1 汉(1+

16、0) = 0.9,ych0.1 0.9 (1 0.1 0.9) =0.9019,1 丄如=? (0.9 +0.9019) =0.90095fyp =0.90095 0.仆0.90095 汇(1+0.仆0.90095) =0.80274,yc =0.90095 -0.1 0.80274(1 0.2 0.80274) =0.80779,1丄y2 =(0.80274 +0.80779) =0.80526L 2这样继续计算下去，其结果列于表Euler方法改进的Euler方法准确值Xnynyny(Xn)0.10.90000000.90095000.90062350.20.80190000.8052632

17、0.80463110.30.70884910.71532790.71442980.40.62289020.63256510.63145290.50.54508150.55761530.55634600.60.47571770.49055100.48918000.70.41456750.43106810.42964450.80.36108010.37863970.37720450.90.31454180.33262780.33121291.00.27418330.29235930.2909884从表9.1可以看出，Euler方法的计算结果只有 2位有效数字，而改进的Euler方法确有3位有效数字

18、,这表明改进的Euler方法的精度比Euler方法高. 2 Runge-Kutta 方法上一节介绍的Euler方法、梯形方法及改进的Euler方法都是单步法，即计算yn+i时，只2用到yn Euler方法及后退的Euler方法的局部截断误差是 0(h)，总体截断误差是 0(h)，梯 形方法的总体截断误差是 0(h2). 一个方法的总体截断误差若为 0(hp)，则我们称它为p阶 方法。体截断误差和局部截断误差的关系是 ：总体截断误差=0(hJ 局部截断误差一般地，方法的总体截断误差阶越高，该方法的精度也越高。在上一节我们利用 y(Xn i)在点Xo的Taylor展开的前两项导出 Euler公式很

19、自然地，我们想到：若将y(xn1)在点X。的Taylor展开多取几项，则有希望获得高阶方法。但是直接利用Taylor展开取多项的办法，需要计算f (x, y)的高阶导数，计算量较大。因此在建立Runge-Kutta 方法时不采用求高阶导数的办法，而是用计算不同点上的函数值，然后对这些值作线性组合,构造近似公式，把近似公式和解的Taylor展开相比较，使它们在前面的若干项相同，从而使近似公式达到一定的阶。二阶 Runge-Kutta 方法我们已经知道Euler方法每一步只计算一次f(x,y)的值，总体截断误差为 0(h)，它的计算式可以写成yn1 二 yn hk(y。已知)& = f (Xn，y

20、n)而改进的Euler方法，每迭代一步要计算两个函数的值，其总体误差为0(h2)，它的计算公式可改写成hh= yn +：ki +：k2,22* ki = f (Xn, Yn),k2 = f (Xn +h, yn +hki).以计算两次函数值为例，设一般计算式为yn 1=ynh(Gk1 C2k2),* kj = f (Xn, yn),(2.1)匕=f (Xn + ph,yn +qhki).适当选择参数g,C2,p,q的值，使得在y(Xn)二yn的假设下，截断误差y(XnHi)yH阶尽可能高。为此，将yn yn h(Ciki C2k2)在点(Xn,yn)作Taylor展开，因为ki = f (Xn

21、, y(Xn) =y (Xn),k2的展开式为k2 二 f (Xn,y(Xn) ph(Xn,y(Xn) qhki(Xn, y(Xn)cxcy12!ccph = +qkh excy jf (Xn, y(Xn) HI3nf (Xn,y(Xn) ph f (Xn, y (x.) qhf (x., y( x.) ( x., y(Xn) - 0(h ),所以yn 勺二 y(Xn) (GC2)hf (Xn，y.)C2ph2 fx(Xn，yn) 9 fy(Xn，yn) f (Xn，yn) O(h)P再将y(Xn 1)在x = Xn作Taylor展开h23y(Xn 1)= y(Xn) hy (Xn) 2!【f

22、x(Xn,yn)f y ( Xn , yn ) f (Xn , Yn) O(h ).若要求局部截断误差达到O(h3)，则通过比较上面两式知，参数c1,c2, p,q必须满足(2.2)(2.3) 1 1c1 c2 = 1, c2 p, c2q =2 21这是四个未知量，三个方程式的方程组，因此解不惟一。当我们选取&=, p = q=1时,2(2.2)和(2,.3)的前三项完全一致。R-K方法)的一种迭代方法将参数值代入(2.1)便得到二阶Runge-Kutta方法(简记为hh,yn =y-k-k2,22也=f (Xn, yn),k2 = f (Xn +h,yn +hkj,它的局部截断误差为 0(

23、h2)，(2.4)实际上就是改进的Euler方法(1.15).1若取q=0, c2=1, p=q ，则得到二阶R-K方法的又一种迭代公式2，广yn4f = yn +hk2,* K = f (Xn,yn),k2 = f(Xn +；, yn +;kj,L.22或写成h1ynynhf (Xn,nhf (Xn, yn).22(2.4)(2.5)公式(2.4)、(2.5)都是显式单步二阶公式.可以证明了不论这四个参数如何选择，都不能使局部截断误差达到0(h4).这说明在计算两次函数值的情况下，局部截断误差的阶最高为0(h3)，要再提高阶就必须增加计算函数值的次数。、三阶及四阶 Runge-Kutta方法

24、为了提高方法的精度，考虑每步计算三次函数f (x, y)值。根据两次计算函数值的做法，很自然地取yn 1的形式为yn = yn +h(Ciki +劭2 +叭)，kl = f (Xn,yn),k f (Xn 叭 yn b2ihkj,k f (Xn a2h, yn bs/k bszhk?).适当选择参数c、G、5Ora2、ah、di、g、5,使截断误差y(xn)yn申的阶尽可能高。由于在推导公式时只考虑局部截断误差，故设yn二y(xn)。类似前面二阶 r-k公式的推导方法，将yn “在(Xn,yn)处作Taylor展开，然后再将 y(xn 在x = xn处作Taylor展开(这 里不详细写出展开式

25、)，只要两个展开式的前四项相同，便有y(xn .J - yn彳=0(h4)，而要两个展开式的前四项相同，参数必须满足G 0203 hk3 = f (x-h, yn +：k2), 2k f (Xn - h, yn - hk3).公式(2.8)就是一种常用的 四阶Runge-Kutta公式，也称为 标准的(或经典的)四阶R-K方法。三阶、四阶Runge-Kutta公式都是单步显式公式。02 0312 2a C2a 2 C3(2.6)这是8个未知数6个方程的方程组，解不是惟一的，可以得到很多公式。当我们取方程组的1 1141一组解为a1,a2=1,b21,b= -1,b32= 2, s , c2,

26、c3时，就得到一2 2 6 6 6种常用的三阶Runge-Kutta公式hyn+ = % +二(匕 +4k2 +k3),6(2.7)匕=f区皿),k2 = f (Xn +?,yn +尹)， 皿=f (Xn +h, yn hk1 +2hk2).从上面导出公式(2.7)的过程知，它的局部截断误差为o(h4).如果每步计算四次函数f (x,y)的值，完全类似的，可以导出局部截断误差为0(h5)的四阶Runge-Kutta公式。详细推导过程这里略去，只给出结果1Yn+=y-(k2k2kk4),6k1 = f (xn, Yn), 丄h丄h(2.8)伙2 = f (Xn +；,yn2 2例2试用Euler

27、方法、改进的Euler方法及四阶经典 R-K方法在不同步长下计算初值问题単=y(1 xy),0 沁叮，dxy(0)在0.2、0.4、0.8、1.0处的近似值，并比较它们的数值结果.解 对上述三种方法，每执行一步所需计算f (x, y) - -y(1 xy)的次数分别为1、2、4。为了公正起见，上述三种方法的步长之此应为1:2:4。因此，在用 Euler方法、改进的 Euler方法及四阶经典 R-K方法计算0。2、0。4、0。8、1。0处的近似值时，它们的步长应分别取 为0。05、0。1、0。2，以使三种方法的计算量大致相等。Euler方法的计算格式为yn 1 二 yn -0.05 yn(1 X

28、nyn).改进的Eluer方法的计算格式为yp =yn 0.1xyn(1 +Xnyn), yc =yn -0.1xyp(1+Xn*yp),1yn卅二尹卩f四阶经典R-K方法的计算格式为=yn *6* (k1 *2k2 *2k3 * k4 ),k1 = yn(1 +xnyn),0 20- kJ,2号 *2),20 2低2=-Wn +-Xk1)1+(Xnk3 = 一Wn +52)口 +(Xn2k4 二-(Yn0.2 k3)1(Xn0.2)( Yn 0.2 & )初始值均为y0二y(0) =1，将计算结果列于表 9.2.表9.2Euler方法(步长 h=0.05)改进的Euler方法(步长 h=0.

29、1)四阶经典R-K方法(步长 h=0.2)准确解Xnynynyny(Xn)0.20.80318660.80526320.80463630.80463110.40.62717770.63256510.63146530.63145290.60.48255860.49055100.48919790.48918000.80.36930360.37863970.37722490.37720451.00.28274820.29235930.29100860.2909884从表9.2可以看出，在计算量大致相等的情况下，Euler方法计算的结果只有2位有效数字，改进的Euler方法计算的结果有 3位有效数字，

30、而四阶经典R-K方法计算的结果却有5位有效数字，这与理论分析是一致的。例1和例2的计算结果说明，在解决实际问题时，选择恰当的算法是非常必要的。需要指出的是Runge-Kutta方法的基于Taylor展开法，因而要求解具有足够的光滑性。 如 果解的光滑性差，使用四阶 Runge-Kutta方法求得数值解的精度，可能不如改进的Euler方法精度高。因此，在实际计算时，要根据具体问题的特性，选择合适的算法。三、变步长的 Runge-Kutta方法上面导出的Runge-Kutta方法都是定步长的，单从每一步来看，步长h越小，局部截断误差也越少，但随着步长的减小，在一定范围内要进行的步数就会增加，而步数

31、增加不仅增加 计算量，还有可能相起舍入误差的积累过大。由于Runge-Kutta方法是单步法，每一步计算步长都是独立的，所以步长的选择具有较大的灵活性。因此根据实际问题的具体情况合理选择 每一步的步长是非常有意义的。正如第五章建立变长的Simps on求积公式那样，我们来建立变步长的Runge-Kutta公式。以四阶标准Runge-Kutta公式为例进行说明， 从基点x0出发，先选一个步长h，利用公式 (2.8)求出的近似值记为 y；h),由于公式的局部截断误差为O(h5),所以有(k). 2y(xjy1ch ,(2.9)其中c为常数，然后将步长h进行折半，即取步长为-,由基点x0出发，到仝0

32、 X1算一步，由2 2x 亠 x(h)01的X1再算一步，将求得的近似值记为 y12 ,因此有2y(N)-ypTc.Y j.(2.10)12丿由(2.9)及(2.10)得y(xj-yf 丄y(G-y1h) 16.一般地，从xn出发，照上面的做法也可得到y(Xn-y俘、丄y(Xn+)-y黑 16，或写成y(Xnd 丄丫存)_ ynh1.(2.11)15(2.11)是事后估计式，记心=皿丫)-yf1(2,12)对给定的步长 h,若二 (；是预先指定的精度)，说明步长h太大，应折半进行计算，直至 A 呂为止，这时取yf2作为近似值；如果 A s为止，这时 再将步长折半一次，就把这次所得结果作为近似值

33、。变步长的Runge-Kutta方法的计算步骤如下：设误差上限为；，误差最小下限为，M . 1，步长最大值为ho，从Xn出发进行计算。M 步长为h .(1) 用步长h和Runge-Kutta公式计算丫补阳，用步长号计算两步得丫畀)，并计算也；(2) 若厶；，说明步长过大，应将 h折半，返回(1)重新计算；(3) 若，说明步长过小，在下一步将h放大，但不超过h0.M这种通过加倍和减半处理步长的Runge-Kutta方法就称为变步长 Runge-Kutta方法。从表面上看，为了选择步长，每一步的计算量增加了，但从总体考虑还是合算的。 3线性多步法上一节介绍的Runge-Kutta方法是单步法，计算

34、 yn1时，只用到y.,而已知信息、 yn,等没有被直接利用。可以设想如果充分利用已知信息ynv, yn,，来计算yn1，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。构 造Runge-Kutta公式是利用Taylor展开的方法，本节利用数值积分公式来构造线性多步法公式。 为清晰简明起见，只介绍简单的线性多步法，至于一般情形可完全类似推导出来。考虑如下形式的求解公式k 4k丫.厂*4 2 ：ifn_i,(3.1)i =0=4其中 fi 三 f(Xj,%),: i (i 二 0,1|l|,k-1), (j 二-1,0,ll),k-1)为与 n无关的常数。由于按

35、公式(3.计算yn 1时，需要知道yn, yn4,IH, yn* 1这k个值，所以称为 k步法，又yn4(i =0,1/ ,k-1)及fn4(i = T,o,1，k-1)都是线性的，所以称(3.1)为线性多步(这里 是k步)法。当- 0，公式中不含fn 1，这时公式就显式的；当 4 0，称为隐式的。我们知道常微分初值问题(1.1)与积分Xn_1 y(Xn_tt) = y(Xn) + f(t,y(t)dt(3.2)Xn是等价的。对d = f (x,y(x)作多项式插值，利用插值型求积公式，便可得到相应的线性多 dx步法公式。下面我们讨论最简单的多步法 Admas方法、Adams显式与隐式公式为了

36、导出求微分方程(1.1)的数值解法，我们将(3.2)右边积分的被积函数用插值多项式来近 似替代。从插值角度来看，插值多项式次数越高越精确，但不能过高(因为高次插值会出现 龙格现象)，这里，我们取三次插值多项。插值节点除X,，Xn d外，通常还要在XnXnv内再由于已知取两点作为插值节点，但取区间Xn , Xn 1 内的点时，函数值又不知道。 yn，ynd,yn,，所以我们取插值节点Xn,XnJ,Xn,Xn .1，作三次插值多项式 的=(t j()(tW、fg,yg)(xn+ _xn)(xj _x n _l)(xn 卅一人 _2 )0 Xn1)(t-Xn)(t 心 f(xn,y(xn)(XXn1

37、)(XXnJ)(X-Xn)n n(一气严-胃厂心)、f (ye)(Xn- Xn 1)(Xn-Xn)(Xn-Xn)+(t Xn 卅)(t Xn_1)(t Xn)(Xn-Xn1)(Xn-Xn)(Xn/-Xn)以心 X21 3(t -Xn)(t -人住-XnRf (Xn 1 , y(Xn 1) 6h13 (t -Xn 1)(t -Xn4)(t -Xn/)f (Xn, yg)(3.3)2h13(t -人(t -Xj(t -XnJf (Xn丫区)2h1 3(t Xn J(t Xn)(t Xn J f (X,y(Xn) 6h其中人=备一备0二n1,n,n 1).插值余项为F(4)g)(3.4)(3.5)y

38、(Xn)、b(t)(t -Xn i)(t -Xn)(t -Xn4)(t -Xn),4!其中 F(x)二 f(X, y(X),Xn 才：Xn 1，故有f (t,y(t) = P3(t) +3(t).将(3.5)代入(3.2)右端的积分，并略去r3(t)，令t =xnuh，再将y(xnj)、y(xn)、y(Xn 1)分别用近似值ynN、ynd、yn、yn 1表示，得1Jyn 1 in hf(Xn 1, Yn 1)0U(U 1)(U 2)dU1 1-2hf(Xn，Vn) 0(U -1)(U1)(U 2)du(3.6)112hf (Xnj,ynj) 0(U -1)U(U 2)dU11-hf(Xn,yn

39、) 0(U -1)U(U 1)dU1 =ynh9f (Xn 1, yn 1) 19f(Xn,yn) - 5f (Xn，yn) f(Xn/,yn2).(外推)公式。24(3.6)称为四阶隐式Adams若插值节点取为Xn、Xn、XnN、Xn：，则三次插值多项式为1P3(t)3(t -Xn)(t -Xn(t -Xn)f (Xn,y(Xn)6h13(t -Xn)(t -Xn/t-Xnf (Xn，y(Xn)2h16h3(t Xn)(t Xn)(t Xn)f (Xn/,y(Xn/)(t -Xn)(t -Xn)(t -Xnf (Xn；,y(Xn),插值余项为MjFSgxStxSXn,宀 Xn,于是有f(t,

40、y(t)二 p3(t) r3(t).(3.7)将(3.7)代入(3.2)右端，略去余项，积分，并用yi近似替代y(Xj) (i二n-3, n-2, n-1, n)，便得1ynynh55 f (Xn,yn)-59f(Xyn4)37 f (x.2y./ -9f(x.，y.；)(3.8)24这就是四阶显式 Adams (内插)公式。由余项表示式(3.4)得，四阶隐式 Adams公式的截断误差为竹+ 1(4)起T十=(tXn 十)(tXn)(tXnJ(tXnQdt.冷24令t =xn uh，由第二积分中值定理得Tn 1F )( )(t Xn i)(t Xn )(t Xn_|)(t Xn _2 )dU2

41、4(3.9) 19 .5(5)”卄=一莎 h y ()凡匸：:：Xn 1.这表明四阶隐式 Adams方法的局部截断误差为0(h5)。同理可得四阶显式 Adams方法的误差为251Tn卑=7h5yL),Xn Xn.(3.10)即四阶显式Adams方法的局部误差也为 O(h5)。、初始值的计算在Adams公式(3.6)或(3.8)中需要多个初始值才能进行计算，例如用公式(3.6)计算yn 时就需要知道yn、yn二、ynN，而微分方程(1.1)只提供了一个初始值，还有两个初值需要用其 他方法来计算。因此，一般来说线性多步法不是自动开始的方法。因为微分方程的初始值在 定解时起着重要作用，所以补充的初始

42、值精度要高一些。对于四阶Adams公式而言(隐式的或显式的)，局部截断误差为 O(h5)，因此用于确定补充初始值的其他方法的局部截断误差至5少应不低于O(h )，否则由于初始值不准确，会影响到最后的结果。对四阶Adams公式而言，最常用的补充初始值的方法是用四阶Runge-Kutta方法，最好用较小步长，如来计算Adams2公式中所缺少的那个初值。还可以用Taylor展开法，如将y(x)在x = 处作Taylor展开1 2y(x) =y(x) y (x)(x-X0)弓 y (x)(x-X0),取它的前若干项，使余项满足精度要求。四阶显式Adams公式，在补充缺少的初值后就可计算了。每步只算一个

43、函数值，局部截 断误差可达到 O(h5)，这是单步法难以达到的。而隐式Adams公式，每步需要解一个非线性方程。与显式 Adams公式相比，在计算上增加了许多困难。但是隐式方法的精度及稳定性比 显式方法好得多(同阶的Adams显式与隐式相比)。三、预测一校正方法在实际应用中，为了保留隐式Adams方法的优点，一般将显式公式与隐式公式联合使用, 用显式公式提供预测值，用隐式公式加以校正，这样，即不用求解非线性方程，又能达到较 高的精度。这种方法称为预测一校正方法。而一般Adams预测一校正过程 是h(3.11)y =yn +(55yn +59yn 二十 37yn 9yn)， 24:yn J =

44、yn +9 f(Xn+,yni)+19y：5yn+ y；/,I24其中 y 二 f(Xj,yj, i 二n, n 1, n 2, n 3.(3.11)表明一般Adams预测一校正过程就是用显式Adams公式给出较好的近似值，再用隐式Adams公式进行迭代可以证明，在一定条件下(k)yn 1 r yn 1(kr，),于是，可在迭代几次后，取ynk1作为yn 1的近似值。在实际应用中，我们希望用隐式公式迭代时，迭代的次数不要太多，如果不进行迭代最好。为此，仿照改进的 Euler方法，我们也可以构造下列Adams预测一校正公式：h预测：yn 1 二 yn 刃(55yn -59yn3跖/ -刖心),Y

45、n 1 二 f (Xn 1, Yn 1)；h(3.12)校正：yn 1 二 yn 24(9yn 1 T9yn -5yn-y,yn 1=f(Xn.1,yn1).上述预测一校正公式是四步法，它在计算yn1时不但要用到前一步的信息yn、yn,而且要用到更前面三步的信息yn、yn,、ynJ3，因此它不能自行启动。在实际计算时，可借助某种单步法，例如可用四阶Runge-Kutta方法提供启动值y1、y2、y3。般来说，Adams预测一校正公式(3.12)虽然计算量小，但其精度比一般 Adams预测一校正过程(3.11)要低。为了改善简单Adams预测一校正公式(3.12)的精度，我们通过对显式及隐式Ad

46、ams方法进行误差分析，得到下列修正的Adams预测一校正公式。用四阶显式Adams公式作预测，第n步得到的结果为 yn51，用四阶隐式Adams公式计算第n步得到结果为yd。由(3.9)及(3.10)得Tnh = y(Xndl) yhTn1 = y(Xn1)-yn1251 h5y(5)( 1), Xn A ：1 ： Xn.720195 .h yOXn/ 耳2 Xn4t.720(3.13)将两式相减，得y5)(養5ysc p 195yn 1 - yn 1h720假定解y(X)的五阶导数y5)(x)在求解区间内连续，且变化不大，则在-Xn,Xn-1 内必存在一点使19 , 52515 (5)27

47、0,5 (5) , 7290 hy(2)莎hy f 莎hy (),于是h5y(5)()720270cp (yn 1 f yn 1 ).这样就得到p 251 cp 、y(Xn1) - yn1(ynlLyn l),(3.14)270c19 cp 、y(Xn 1) - yn 1 ： - 270 ( yn 1 - yn 1).(3.14)是(3.13)两个公式的修正公式，由此可给出修正的Adams预测一校正公式预测：Ynp+ =Yn + (55fn 59fn+37 f：/ 9件)，24修正：m：1 =Yn270(y Yn )，求f:m： 1 二 f(Xn1，mn1),(3.15)校正：chyn 1 = yn(9mn 1 19fn -5fnfn/),24修正：YnYn 1(Yn 1 Yn1)，270求导：Yn 1 =f(X 1，Yn1).修正的Adams预测一校正公式(3.14)也是四步法，它在计算yn 时要用到前面四步的信息Yn、yn、Yn、Yn N、丫2和丫：-，因此它在开始计算之前必须先给出启动值、丫2、y3和Ys -y3p。同Adams预测一校正公式(3.12) 样y1、y2、y3可用其他四阶单步法(如四 阶经典Runge-Kutta方法)提供，而令 y； - yf =0。由(3.15)得c19cp