数值分析常微分方程数值解

dydxf(x,y)，axb (9-1) y(a)y0的数值解法具有典型性。

常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。用解析方法只能求出线性常系数等特殊类型的方程的解。对非线性方程来说，解析方法一般是为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。因此研究微分方程的数值解法是非常必要的。

只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。由常微分方程的理论知，如果(9-1)中的f(x,y)满足条件

（1）f(x,y)在区域D{(x,y) axb，y}上连续； （2）f(x,y)在上关于满足Lipschitz条件，即存在常数，使得

f(x,y)f(x,y)Lyy

则初值问题(9-1)在区间[a,b]上存在惟一的连续解yy(x)。在下面的讨论中，我们总假定方程满足以上两个条件。

所谓数值解法，就是求问题(9-1)的解yy(x)在若干点

ax0x1x2xNb

处的近似值yn(n1,2,,N)的方法。yn(n1,2,,N)称为问题(9-1)的数值解，hxn1xn称为由到xn1的步长。今后如无特别说明，我们总假定步长为常量。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数

在问题(9-1)中，若用向前差商

y(xn1)y(xn)h代替y(xn)，则得

y(xn1)y(xn)hf(xn,yn(xn)) (n0,1,,N1)

y(xn)用其近似值代替，所得结果作为y(xn1)的近似值，记为yn1，则有 yn1ynhf(xn,yn) (n0,1,,N1)

这样，问题(9-1)的近似解可通过求解下述问题

yn1ynhf(xn,yn) (n0,1,,N1)y0y(x0)(9-2)

得到，按式(9-2)由初值经过步迭代，可逐次算出y1,y2,yN。此方程称为差分方程。

需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 (2) 用数值积分法

将问题(9-1)中的微分方程在区间[xn,xn1]上两边积分，可得

y(xn1n1)y(xn)xxf(x,y(x))dx (n0,1,,N1)(9-3)

n用yn1,分别代替y(xn1),y(xn)，若对右端积分采用取左端点的矩形公式，即

xn1xf(x,y(x))dxhf(xnn,yn)

同样可得出显式公式(9-2)。

类似地，对右端积分采用其它数值积分方法，又可得到不同的计算公式。 (3) 用Taylor多项式近似。把y(xn1)在点处Taylor展开，取一次多项式近似，则得

y(x)hy(xh2n1)y(xnn)2!y()2 y(x)hf(x(xhnn,yn))2!y() [xn,xn1]设h1，略去余项，并以代替y(xn)，便得 yn1ynhf(xn,yn)

以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。

上面我们给出了求解初值问题(9-1)的一种最简单的数值公式(9-2)。虽然它的精度比较低，实践中很少采用，但它的导出过程能较清楚地说明构造数值解公式的基本思想，且几何意义明确，因此它在理论上仍占有一定的地位。

1简单的数值方法和基本概念

1.1 Euler法与向后Euler法

一、Euler法

Euler方法就是用差分方程初值问题

yn1ynhf(xn,yn) (n0,1,,N1)y0y(a)(9-4)

的解来近似微分方程初值问题(9-1)的解，即由公式(9-4)依次算出y(xn)的近似值yn(n1,2,)。

从几何上看，微分方程yf(x,y)在xoy平面上确定了一个向量场：点(x,y)处的方向斜率为f(x,y)。问题(9-1)的解yy(x)代表一条过点(x0,y0)的曲线，称为积分曲线，且此曲线上每点的切向都与向量场在这点的方向一致。从点P0(x0,y0)出发，以f(x0,y0)为斜率作一直线段，与直线xx1交于点P1(x1,y1)，显然有

y1y0hf(x0,y0)，再从出发，以f(x1,y1)为斜率作直线段推进到xx2上一点P2(x2,y2)，其余类推，这样得

到解曲线的一条近似曲线，它就是折线P0PP12。因此Euler方法又称为Euler折线法。

二、向后Euler法

在微分方程离散化时，用向后差商代替导数，即y(x(xn1)y(xn)n1)yh，则得到如下差分方程

yn1ynhf(xn1,yn1) (n0,1,,N1)y0y(x(9-5)

0)用这组公式求问题(9-1)的数值解称为向后Euler法。

向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。Euler法计算yn1的公式中不含有yn1，这样

的公式称为显式公式；向后Euler法计算yn1的公式中含有yn1，称为隐式公式。显式公式与隐式公式各有特点。显式公式的优点是使用方便，计算简单，效率高。其缺点是计算精度低，稳定性差；隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。

上面隐式公式中，在求解yn1时，yn,xn1为已知，yn1是方程yn1ynhf(xn1,yn1)的根。一般说来，这是一个非线性方程，因此我们通过构造简单迭代法来求解。迭代格式为

y(0)n1ynhf(xn,yn)y(k1)x(k) n1ynhf(n1,yn1) (k0,1,2,)由于f(x,y)满足Lipschitz条件，所以

y(k1)n1y(k)(k)n1hf(xn1,yn1)f(xn1,yn1)hLyn1yn1

由此可知，只要hL1，迭代法就收敛到解yn1。

1.2 梯形公式

利用数值积分方法将微分方程离散化时，若用梯形公式计算式(9-3)中右端积分，即

xn1xf(x,y(x))dxh[f(xn,y(xn))f(xn1,y(xn1))] n2并用yn,yn1代替y(xn),y(xn1)，则得计算公式

yyhn1n2[f(xn,yn)f(xn1,yn1)](9-6)

这就是求解初值问题(9-1)的梯形公式。

梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为

y(0)n1ynhf(xn,yn)(k1)hyn1yn2f(x(k) n,yn)f(xn1,yn1) (k0,1,)由于函数f(x,y)关于满足Lipschitz条件，所以

y(k1)1yn1h2f(x(k)hL(k)nn1,yn1)f(xn1,yn1)2yn1yn1 其中为Lipschitz常数。因此，当0hL21时，迭代法就收敛到解yn1。

1.3 局部截断误差与方法的精度

为了刻画近似解的准确程度，引入局部截断误差与方法精度的概念。

定义 假设在某一步的近似解是准确的，即yny(xn)（这个假设称为局部化假设）。在此前提下，用某公式推算所得yn1，我们称

Rn1y(xn1)yn1

为该公式（即该方法）的局部截断误差。

定义9.2 如果某种方法的局部截断误差是

Rp1n1y(xn1)yn1O(h)

则称该方法是阶方法，或具有阶精度。显然越大，方法的精度越高。

1）Euler法的截断误差

假设问题的解y(x)充分光滑，且前步计算结果是精确的，即

yiy(xi)，y(xi)f(xi,y(xi))(in)

于是Euler法的截断误差是

Rn1y(xn1)yn1y(xn1)ynhf(xn,yn)y(xhy(x

n1)y(xn)n)h22y(xn)O(h3)(9-7) 这里h2y(x22n)称为局部截断误差主项。显然Rn1O(h)。 2）向后Euler法的截断误差。计算公式是

yn1ynhf(xn1,yn1)

将f(x,y)对用微分中值定理，有

f(xn1,yn1)f(xn1,y(xn1))fy(xn1,)(yn1y(xn1))（在yn1与y(xn1)之间）将f(xn1,y(xn1))在处Taylor展开

f(xn1,y(xn1))y(xn1)y(xn)hy(xn)O(h2)

于是

yn1y(xn)hy(xn)h2y(xn)hfy(xn1,)(yn1y(xn1))O(h3)

将方程的解作Taylor展开

y(x)hy(xh2n1)y(xnn)2y(xn)O(h3)

因此

y(xn1)ynh212y(xhn)2fy(xn1,)(yn1y(xn1))O(h3)

故

Rn1y(xn1)y11hf,)h2n12y(xn)O(h3)y(xn1h2

1hfy(xn1,)32y(xn)O(h)h22y(xn)O(h3)(9-8)

3）梯形法的计算公式是

yhn1yn2[f(xn,yn)f(xn1,yn1)]

将f(x,y)对用微分中值定理，有

f(xn1,yn1)f(xn1,y(xn1))fy(xn1,)(yn1y(xn1))（在yn1与y(xn1)之间）将f(xn1,y(xn1))在处Taylor展开

(xy(xh2fn1,n1))y(xn1)y(xn)hy(xn)2y(xn)O(h3)

于是

2yy(xh2)h34y(xhn1n)hy(xn)y(xnn)2fy(xn1,)(yn1y(xn1))O(h4)

将方程的解作Taylor展开

h2h3y(xn1)y(xn)hy(xn)2y(xn)6y(xn)O(h4)

因此

3y(xhn1)yny(xh1n)2fy(xn1,)(yn1y(xn1))O(h412) 故

Rn1y(xn1)y11hfh3n1y(x4n)O(h)y(xn1,)12[1hfh3y(xn1,)]12y(xn)O(h4)(9-9)

h3y(xn)O(h412)所以梯形法是二阶方法。

1.4 改进的Euler法

我们看到，虽然梯形方法提高了精度，但其算法复杂，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数f(x,y)的值，而迭代又要反复进行若干次，计算量很大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法。

具体地说，我们先用欧拉公式求得一个初步的近似值yn1，称之为预测值，预测值yn1的精度可能很差，再用梯形公式将它校正一次得yn1，称为校正值。这样的预测校正系统通常称为改进的欧拉公式。即

预测：yn1ynhf(xn,yn) 校正：yn1ynh2f(xn,yn)f(xn1,yn1) 为了编制程序上机，将上式改写成

ypynhf(xn,yn)yqynhf(xnh,yp)(9-10) yn112(ypyq)算法

(1) 输入a,b,f(x,y)，整数，初值

(2) 置hbaN,n0,xa,yy0 输出(x,y)

(3) 计算

ypyhf(x,y)yqyhf(xh,yp)

12(ypyq)y,xhx 输出(x,y)

(4) 若nN1，置n1n，转3；否则停机。 改进欧拉法的截断误差。计算公式是

yhn1yn2[f(xn,yn)f(xn1,ynhf(xn,yn))]

在(xn,y(xn))处作Taylor展开式，注意到yny(xn)，有

yy(xh2xhn1n)f(n,y(xn))2[f(xn,y(xn))hfx(xn,y(xn))hfy(xn,y(xn))f(xn,y(xn))O(h2)]

y(xh2n)hy(xn)y(x32n)O(h)于是

R3n1y(xn1)yn1O(h)(9-11)

所以改进Euler法是二阶方法。

2 龙格-库塔（Runge-Kutta）法

2.1 Runge-Kutta法的基本思想及一般形式

设初值问题(9-1)的解yy(x)C1[a,b]，按微分中值定理，必存在[xn,xn1]，使

y(xn1)y(xn)hy()ynhf(,y())(9-12) y*nhK其中叫做y(x)在[xn,xn1]上的平均斜率。只要对平均斜率提供一种算 法，公式(9-12)就给出一种数值解公式。例如，用K1f(xn,yn)代替，就得到前进Euler公式；用K2f(xn1,yn1)替代，就得到后退Euler公式；如果用K1,K2的算术平均值替代，则可得到2阶精度的梯形公式，如图9-3。可以设想，如果在[xn,xn1]上能多预测几个点的yEn1 斜率值，用它们的加权平均值替代，就有望得到具有较高精度的数值yyTBn1 n1 解公式，这就是Runge-Kutta法的基本思想。

Runge-Kutta公式的一般形式是

r图9-3

yn1ynhciKii1K1f(x(i2,,r)n,yn)(9-13)

i1Kif(xnih,ynhijKj)j1

其中是yy(x)在xnih (0i1)点的斜率预测值。ci,i,ij均为常数。选取这些常数的原则是使公式(9-13)具有尽可能高的精度。公式(9-13)叫做级显式Runge-Kutta公式，简称RK方法。

2.2 二阶RK法的推导

r2的RK方法的近似公式为 yn1ynh(c1K1c2K2)K1f(xn,yn)(9-14) K2f(xn2h,ynh21K1)这里c1,c2,2,21均为待定常数，我们希望适当选取这些系数，使公式阶数尽量高。先展开，按照二元函数Taylor

级数：f(as,bt)n1k0k!sxtf(a,b，得

ky)Kf(x12n,yn)2hfx(xn,yn)21hfy(xn,yn)f(xn,yn)2![22h2fxx(xn,yn)(9-15)

222221h2fxy(xn,yn)f(xn,yn)21hfyy(xn,yn)f2(xn,yn)]为了叙述方便，把f(xn,yn)及其偏导数中的xn,yn省略不写，将式(9-15)代入(9-14)的第一式得

y2n1ynh(c1c2)fhc2(2fx21fyf)h3c2(22

22fxx2221fxyf21fyyf2)再展开y(xn1)，注意到公式

yf(x,y)yfxfyf yf2xx2fxyffyyffy(fxfyf)得

23y(xn1)y(xn)hy(xhn)2!y(xhn)3!y(xn)

hfh2y2(fh3nxfyf)6(fxx2fxyffyyf2fyfxfyf)于是，局部截断误差为

Rn1y(xn1)yn1h(1ch212c11c2)f22fx2c221fyf

h31612c21112222fxx3c2221fxyf62c221fyyffy(fxfyf)(9-16)

要使式(9-16)的局部截断误差为O(h3)，则应要求f,fx,fyf的系数为零，于是有

c1c21c220.5(9-17) c2210.5方程有4个未知数，3个方程，所以有无穷多组解，它的每组解代入式(9-14)得到的近似公式，局部截断误差均为

O(h3)，故这些方法统称为二阶方法。例如，取c1c20.5，得2211，近似公式为

yn1yn0.5h(K1K2)K1f(xn,yn)(9-18) K2f(xnh,ynhK1)这就是改进Euler公式。如果取c10,c21，有22112，近似公式为 yn1ynhK2K1f(xn,yn)(9-19) K2f(xn0.5h,yn0.5hK1)这也是常用的二阶公式，称为中点公式。

注：对于一般函数f(x,y)，由于fy(fxfyf)0，所以不论参数如何选取，只能有Rn1Oh3。这说明(9-14)至多是二阶的方法。

类似地，对r3和r4的情形，通过更复杂的计算，可以导出三阶和四阶RK公式，其中常用的三阶和四阶RK公式为

yn1yhn(K14K2K3)6K1f(xn,yn)(9-20) K2f(xhhn2,yn2K1)K3f(xnh,ynhK12hK2)和

yn1ynh(K12K22K3K4)6K1f(xn,yn)Khh2f(xn,ynK221)(9-21) K3f(xhhn2,yn2K2)K4f(xnh,ynhK3)式(9-21)称为经典的四阶RK公式，通常说四阶RK方法就是指用式(9-21)求解。

算法

(1) 输入a,b,f(x,y),N,y0 (2) 置hbaN,n1,x0a (3) 计算

xx0hK1f(x0,y0)Kf(xhh202,y02K1)Khh

3f(x02,y02K2)K4f(x0h,y0hK3)yyh06(K12K22K3K4)输出(x,y)

(4) 若nN，则置n1n,xx0,yy0转3；否则停机。

例4 用Euler法（h0.025），改进Euler法（h0.05）和经典R-K方法（h0.1）计算初值问题

yy1 (0x0.5)y(0)0. 其结果如下表9-1

表9-1

Euler方法 改进Euler法 经典R-K法 (h0.025) (h0.05) (h0.1) y(xn) 0 0 0 0 0 注：1）通过比较RK四阶经典公式和前面用Euler法、改进Euler法（它是一种二阶RK方法）的计算结果，我们发现，四阶RK方法的精度高得多。就计算量来说，虽然Euler法、改进Euler法每步只需计算一个或二个函数值，而四阶RK方法每步需计算四个函数值，但由于放大了步长，计算量几乎相同。可见，选择有效的算法是非常重要的。

2）RK方法的导出基于Taylor展开，故它要求问题的解具有较高的光滑度。当解充分光滑时，四阶RK方法确实优于改进Euler法。如果解的光滑性较差，则用四阶RK方法求得的数值解，其精度可能反而比改进Euler法的差。因此，在实际计算时，需根据问题的具体情况选择合适的算法。

3）当r1,2,3,4时，RK公式的最高阶数恰是，当r4时，RK公式的最高阶数不是，如r5时，RK公式的最高阶数仍为4，r6时，RK公式的最高阶数仍为5。由此看来，四阶RK公式是最经济的。

2.3 变步长的RK方法

若问题(9-1)的解函数y(x)变化是不均匀的，用等步长求数值解，则可能产生有些点处精度过高，有些点处精度过低的情况，为保证一定的精度，必须取较小的步长。这样做既增加了计算量，也导致误差的严重积累，因而实际计算时，往往采用事后估计误差、自动调整步长的RK方法。

设用阶方法计算y(h)n1，记从出发，以步长计算一步所得y(xn1)的近似值yn1，以步长

h2计算两步得y(xn1)的近似值y(h)n21。由于阶公式的局部截断误差为O(hp1)，且y(p1)在小区间[xn,xn1]上变化不大，故有

y(x(h)chp1n1)yn1(9-22) p1y(xn1)y(hn2)12ch2(9-23)

将式(9-24)乘以减(9-23)得

(h(2p1)y(xn1)2pyn2)1y(h)n10

从而有

hy(xn1)yn2p(h)(2)12p1yn1yn1(9-24) 这样，可以事后误差估计(9-24)中的大小来选择合适的步长。

变步长的RK方法与数值积分做法相同，可以从的值来选择合适的步长，从而得到合乎精度要求的yn1。具体作法是：设要求精度是，如果，就反复加倍步长进行计算，直到为止。这时上一次的步长就是合适的步长，上一次计算所得的结果，就是合乎精度要求的yn1；如果，则反复折半步长，直到为止。这时，最后一次步长就是合适的步长，而这最后一次的计算结果就是满足精度要求的yn1。这种计算过程中自动选择

步长的方法，叫变步长方法。

以上所介绍的各种数值解法都是单步法，这时因为它们在计算yn1时，都只用到前一步的值，单步法的一般形式是

yn1ynh(xn,yn,h) (n0,1,,N1) (9-25)

其中(x,y,h)为增量函数，例如Euler方法的增量函数为f(x,y)，改进Euler法的增量函数为

(x,y,h)12f(x,y)f(xh,yhf(x,y)) 3 相容性、收敛性与稳定性

3.1 相容性与收敛性

上面所介绍的方法都是用离散化的手段，将微分方程初值问题化为差分方程初值问题求解的。这些转化是否合理？即当h时，差分方程是否能无限逼近微分方程，差分方程的解是否能无限逼近微分方程初值问题的准确解

y(xn)，这就是相容性与收敛性问题。

用单步法(9-25)求解初值问题(9-1)，即用差分方程初值问题

yn1ynh(xn,yn,h)y(x0)y (9-26) 0的解作为问题(9-1)的近似解，如果近似是合理的，则应有

y(xh)y(xh)h(x,y(x),h)0 (h0) (9-27)

（体会：要使得差分方程能无限逼近微分方程，我们将微分方程的精确解代人差分方程中，当h0时，应该满足）

其中y(x)为问题(9-1)的精确解。因为

limy(xh)y(xh)h0hy(x)f(x,y)

故由(9-27)得

limh0(x,y,h)f(x,y)

如果增量函数(x,y(x),h)关于连续，则有

(x,y,0)f(x,y) (9-28)

如果单步法的增量函数(x,y,h)满足条件(9-28)，则称单步法(9-25)与初值问题(9-1)相容。通常称(9-28)为单步

法的相容条件。

满足相容条件(9-28)是可以用单步法求解初值问题(9-1)的必要条件。

容易验证Euler法和改进Euler法均满足相容性条件。事实上，对Euler法，增量函数为

(x,y,h)f(x,y)

自然满足条件(9-28)，改进Euler法的增量函数为

(x,y,h)12f(x,y)f(xh,yhf(x,y))

因为f(x,y)连续，从而有

(x,y,0)12f(x,y)f(x,y)f(x,y)

所以Euler法和改进Euler法均与初值问题(9-1)相容。一般地，如果单步法有阶精度（p1），则其截断误差为

y(xh)y(x)h(x,y(x),h)O(hp1)

上式两端同除以，得

y(xh)y(xh)h(x,y,h)O(hp)

令h0，如果(x,y(x),h)连续，则有 y(x)(x,y,0)0

所以p1的单步法均与问题(9-1)相容。由此即得各阶RK方法与初值问题(9-1)相容。

定义9.4 一种数值方法称为是收敛的，如果对于任意初值及任意固定的x(a,b]，都有

limh0yny(x) (xanh)

其中y(x)为初值问题(9-1)的精确解。

如果我们取消局部化假定，使用某单步法公式，从出发，一步一步地推算到xn1处的近似值yn1。若不计各步的舍人误差，而每步都有局部截断误差，这些局部截断的积累就是整体截断误差。

定义9.5我们称

en1y(xn1)yn1

为某数值方法的整体误差。其中y(x)为初值问题(9-1)的精确解，yn1为不计舍人误差时用某数值方法从开始，逐步得到的在xn1处的近似值（不考虑舍人误差的情况下，截断误差的积累）。

定理设单步法(9-25)具有阶精度，其增量函数(x,y,h)关于满足Lipschitz条件，问题(9-1)的初值是精确的，即

y(x0)y0，则单步法的整体截断误差为 en1y(xn1)yn1O(hp)

证明 由已知，(x,y,h)关于满足Lipschitz条件，故存在L0，使得对任意的y1,y2及x[a,b]，

0hh0，都有

(x,y1,h)(x,y2,h)Ly1y2

记yn1y(xn)h(xn,y(xn),h)，因为单步法具有阶精度，故存在M0，使得

Rn1y(xn1)yp1n1Mh

从而有

en1y(xn1)yn1y(xn1)yn1yn1yn1Mhp1y(xn)h(xn,y(xn),h)ynh(xn,yn,h) Mhp1y(xn)ynh(xn,y(xn),h)(xn,yn,h)Mhp1(1hL)en反复递推得

eMhp1(1hL)Mhp1n1(1hL)en11(1hL)(1hL)nMhp1(1hL)n1e0

(1hL)n11Mhp1(1hL)n1hLe0因为y(x0)y0，即e00，又(n1)hba，于是 baba(1hL)n1(1hL)hehln(1hL)eL(ba)

所以

eMn1LhpeL(ba)1O(hp) 推论1 设单步法具有(p1)阶精度，增量函数(x,y,h)在区域： axb, y, 0hh0 上连续，且关于满足Lipschitz条件，则单步法是收敛的。

当f(x,y)在区域D:axb,y上连续，且关于满足Lipschitz条件时，改进Euler方法，各阶RK方法的增量函数(x,y,h)在区域上连续，且关于满足Lipschitz条件，因而它们都是收敛的。

关于单步法收敛的一般结果是：

定理设增量函数(x,y,h)在区域上连续，且关于满足Lipschitz条件，则单步法收敛的充分必要条件是相容性条件(9-28)。

3.2 稳定性

稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长h0时，方法的整体截断误差是否趋于零的问题。而稳定性则是讨论舍人误差的积累能否对计算结果有严重影响的问题。

若一种数值方法在节点值上有一个大小为的扰动，于以后各节点ym(mn)上产生的偏差均不超过，则称该方法是稳定的。

我们以Euler方法为例进行讨论。假设由于舍人误差，实际得到的不是而是ynynn，其中是误差。由此再计算一步，得到

yn1ynhf(xn,yn)

把它与不考虑舍人误差的Euler计算公式相减，并记n1yn1yn1，就有

n1nhf(xn,yn)f(xn,yn)1hfy(xn,)n

其中ffyy。如果满足条件

1hfy(xn,)1， (9-29)

则从到yn1的计算，误差是不增的，可以认为计算是稳定的。如果条件(9-29)不满足，则每步误差将增大。当

fy0时，显然条件(9-29)不可能满足，我们认为问题本身具有先天的不稳定性。当fy0时，为了满足收敛性要

求(9-29)，有时要很小，这样步数就相当多，这时误差的累积可能是十分严重的，出现数值不稳定现象。

现在我们考虑fy0的初值问题

yy,y(x0)a (9-30)

其解是yaexp(xx0)。今设在xx0的初值有误差，实际求解的问题是

yy,y(x0)aa (9-31)

它的准确解是y(aa)exp(xx0)。因此如果用很精确的方法，求出(9-31)的解，则它和(9-30)的解在处误差为，而在处的误差则是(a)exp(xx0)。它随着的增大而增大，与所选的数值方法无关，是问题本身固有的特性。

再看一个fy0的例子，初值问题为

y100y,y(x0)a

其准确解是yaexp100(xx0)。

用欧拉法求解有

yn1(1100h)yn (9-32)

若取h0.025，则欧拉公式的具体形式为

yn11.5yn

同样讨论有初始误差，则在xxnn的误差是1.5a。数值不稳定。

式(9-32)中取h0.005，则欧拉公式的具体形式为yn10.5yn。对初始误差，在xxnn的误差是0.5a。数值稳定。

以上讨论表明稳定性不但与方法有关，也与步长的大小有关，当然也与方程中的f(x,y)有关。为了只考察数值方法本身，一般只检验数值方法用于求解模型方程的稳定性，模型方程为

yy (9-33)

其中为复数。这个方程的分析比较简单，一般的方程可以通过局部线性化为这种形式，例如在(x,y)的邻域

yf(x,y)f(x,y)fx(x,y)(xx)fy(x,y)(yy)

略去高阶项，再作变量替换就得到uu的形式。事实上，方程可简写为yb(axy)c，作变换

uaxyacb即可得ubu。对于个方程的方程组，可以线性化为YAY，其中为mm的Jacobi矩阵fy。若有个不同的特征值T1,2,,m，则可对角化为QAQdiag(1,,m)，再作变换Y=QU，得到i个非耦合的方程UiiUi，其中可以复数，所以一般讨论(9-33)中的为复数。

对于方程(9-33)，若Re0，类似以上分析，认为方程是不稳定的。所以我们考虑Re0的情形，这时不同的数值方法可能是数值稳定或不稳定的。当一个单步法用于试验方程yy，从计算一步得

yn1E(h)yn (9-34)

其中E(h)依赖于所选的方法。因为通过点(xn,yn)试验方程的解曲线（它满足yf(x,y),y(xn)yn）为

ynexp(xxn)，而一个阶单步法的局部截断误差在y(xn)yn时有

Tp1n1y(xn1)yn1O(h)，

所以有

ynexp(h)E(h)ynO(hp1) (9-35)

这样可以看出E(h)是的一个近似值。

由(9-34)可以看到，若计算中有误差，则计算yn1时将产生误差E(h)，所以有下面定义。

定义9.6 如果(9-34)式中，E(h)1，则称单步法(9-25)是绝对稳定的。在复平面上复变量满足E(h)1的区域，称为方法(9-25)的绝对稳定区域，它与实轴的交称为绝对稳定区间。

在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致。事实上，E(h)1时也可以认为误差不增长。

(1) Euler方法的稳定性

Euler方法用于模型方程(9-33)，得yn1(1h)yn，所以有E(h)1h。所以绝对稳定条件是

1h1，它的绝对稳定区域是复平面上以(1,0)为中心的单位圆。而为实数时，绝对稳定区间是(2,0)。

(2) 梯形公式的稳定性

1h对模型方程，梯形公式的具体表达式yhn1yn(ynyn1)，即yn221y1hn，所以梯形公式的绝对2稳定条件为1h21h2。化简得Re(h)0，因此梯形公式的绝对稳定区域为平面的左平面。特别地，当为负实数时，对任意的h0，梯形公式都是稳定的。

Im(h) Im(h) Re(h)

Re(h)

图9-4 图9-5

(3) RK方法的稳定性

与前面的讨论相仿，将RK方法用于模型方程(9-33)，可得二、三、四阶RK方法的绝对稳定区域分别为

1h1(h)221

1h12(h)216(h)31

1h112(h)26(h)3124(h)41

当为实数时，三、四阶显式RK方法的绝对稳定区域分别为2.51h0，2.78h0。 表9-2 常用几个公式的绝对稳定区间

方法 阶数 区间 Euler公式（显式） 1 (2,0) 后退Euler公式（隐式） 1 (,0)(2,) 梯形公式（隐式） 2 (,0) 二级R-K公式（显式） 2 (2,0) 二级R-K公式（隐式） 2 (,0) 三级R-K公式（显式） 3 (2.51,0) 经典R-K公式（显式） 4 (2.78,0) 例5 分别取h0.1和0.2，用经典R-K方法计算

y20y (0x1), y(0)1

解 本例20，分别为和，前者在绝对稳定区间内，后者则不在，这里列出计算误差

表9-3

yn(h0.1) y(xn)yn yn(h0.2) y(xn)yn 1 1 5 25 125 4.57247104 1.52415104 625 5.08052105 1.88167106 3125 从表9-3可以看出，当h0.1时，虽然计算结果的相对误差很大，但计算过程是稳定的。而当h0.2时，计算过程不稳定。这是因为h0.1时，h2属于稳定区间[2.78,0]，当h0.2时，h4不属于稳定区间。

对于一般方程yf(x,y)，在考虑稳定性对步长的时，常用

fy代替（因为前面变换中bfy），只要

步长的选取使得在所用方法的绝对稳定区域内即可。例如若用Euler方法求解初值问题

y20xy1x2 0x2 y(0)1因为fy20x1x20且max0x210，要使1h1，只要

4 线性多步法

前面所讲的各种数值解法，都是单步法。因为它们在计算yn1时，都仅仅用到前面的一个信息；如果在计算值

yn1时，能够比较充分地利用前面的已知信息，如yn,yn1,,ynr，那么就可望使所得到的yn1更加精确。这就是

多步法的基本思想。

多步法中最常用的是线性多步法，其一般公式是

ryn1iynihrifni (9-36)

i0i1其中k,k均为常数，ynky(xnk),fkf(xk,yk) (kn1,,nr)。

当10时，上式称为显式，否则称为隐式。

线性多步法与Runge-Kutta法都是高精度方法，不同的是线性多步法是利用前面已求得的y(x)在各节点处的近似值及导数近似值的线性组合来逼近y(xn1)的，而Runge-Kutta法则是用y(x)在[xn,xn1]内若干点处的一阶导数预测值的线性组合，来逼近y(x)在区间[xn,xn1]内的平均斜率，从而求得y(xn1)的近似值。构造线性多步法的公式有多种途径，最常用的是Taylor展开方法和数值积分的方法.

定义设y(x)是初值问题(9-1)的精确解，线性多步法(9-36)在xn1上的局部截断误差为

rrTn1y(xn1)kynkhk 0kkynk1若Tp1n1O(h)，则称方法(9-36)是阶的，p1则称方法(9-36)与问题(9-1)的第一个方程相容的。

4.1 线性多步法公式的推导

利用Taylor展开导出线性多步公式(9-36)的基本方法是：将线性多步公式(9-36)在处进行Taylor展开，然后与

y(xn1)在处的Taylor展开式相比较，要求它们前面的项重合，由此确定参数k,k。设初值问题(9-1)的解y(x)充

分光滑，记y(j)(j)ny(xn) (j0,1,)，把它在处展成Taylor级数则有

p(x1y(x)yn)(j)(xxn)(p1)nyn (9-37)

j1j!(p1)!yn及

py(x)(xxj1n)(j)(xxpn)(p1)(j1)!ynp!yn (9-38)

j1用xnkx0(nk)h替代式(9-37)、(9-38)中的，则得

pyy(x(kh)j(j)(kh)p1(p1)nknk)ynj!yn(p1)!yn

j1及

pyy(x(kh)j1y(j)(kh)p(p1)nknk)nj1(j1)!p!yn

把这两个式子代入式(9-36)，得

rpjp1yy(kh)(j)(kh)(p1)n1knk0ynynj1j!(p1)!rp(kh)j1phk(j)(kh)(p1)k11)!ynj1(jp!yn(9-39)

rphjrjr(j)kyn0j1j!k(k)jk0kk(k)j11ynkhp1(p1)!rp1rp(pk(k)(p1)k1k1)k(k)1yn为了使式(9-36)具有阶精度，只须使式(9-39)的前p1项与y(xn1)在处的Taylor展开式

py(xhj(j)hp1(p1)n1)ynynyn(9-40)

j1j!(p1)!的前p1项系数对应相等即可。对比关于的同次项系数，得到确定k,k的方程组：

rk1k0 (j1,,p) (9-41) rr(k)jkj(k)j1k1k1k1可见，只要式(9-36)的系数k,k满足式(9-41)，方法(9-36)就具有阶精度。式(9-40)减去(9-39)，可得局部截断误差为rTn1hp1r(p1)!1(k)p1k(p1)(k)p(p1)O(hp2) (9-42) k1kk1yn4.2常用的线性多步法公式

式(9-41)共有p1个方程，2r3个待定常数：0,1,,r,1,0,1,r，只要2r3p1，即

rp21，就可以由式(9-41)，定出具有阶精度的公式(9-36)的系数k,k。 下面介绍几个常用的线性多步法公式。 一、阿达姆斯（Adams）公式 取p4,r3，有

3k1k0 (j1,2,3,4)(9-43) 33(k)jj(k)j1kk0kk11此时方程有9个未知数，5个方程，有四个自由变量。特别取12310，可解得

555937911,024,124,224,324

相应的线性多步法公式为

yhn1yn24(55fn59fn137fn29fn3)(9-44) 因为10，式(9-44)称为Adams显式公式，它是四阶公式，局部截断误差为

5Th33n15!1(k)5)4(5)k5(kkk1k1ynO(h6)251 (9-45)

h5y(5)nO(h6720)如果令12330，由式(9-43)可得解

9195101,124,024,124,224

相应的线性多步法公式为

yhn1yn24(9fn119fn5fn1fn2)(9-46) 因为10，式(9-46)称为Adams隐式公式，它是四阶公式，局部截断误差为

T19h5y(5)6n1720nO(h) (9-47)

二、米尔恩（Milne）公式

对方程组(9-42)，令01210，解出

31,84803,13,23,30

相应的线性多步法公式

yn1yn34h3(2fnfn12fn2) (9-48) 称为Milne公式，其局部截断误差为

T14n145h5y(5)nO(h6) (9-49) Milne公式是四阶四步显式公式。

三、海明（Hamming）公式

对方程组(9-41)，取r2，并令120，可得Hamming公式

y13n18(9ynyn2)8h(fn12fnfn1) (9-50)

其局部截断误差为

Tn1140h5y(5)nO(h6) (9-51) Hamming公式是四阶三步隐式公式。

由于线性多步法公式(9-36)，只有在知道了前面的yn,yn1,,ynr之后，才能使用。所以开头的y1,y2,,yr的

值必须先用同阶的单步法（Taylor级数法或Runge-Kutta法）求出，然后才能用线性多步法公式（例如，若使用四阶

线性多步法公式，必须用同阶单步法，算出y1,y2,y3）。

例6 分别用四阶Adams显式和隐式公式求初值问题

yxy 0x1y(0)0 的数值解，取h0.1。

解 根据题意，xnnh0.1n,fnxnyn，由四阶Adams显式公式(9-43)有

yhn1yn24(55fn59fn137fn29fn3)

124(18.5yn5.9yn13.7yn20.9yn30.24n0.12)(n3,,9)

由四阶Adams隐式公式(9-44)有

yn1yhn24(9fn119fn5fn1fn2)1

24(0.9yn122.1yn0.5yn10.1yn20.24n0.12)(n2,,9)

由上式可解出

y1n124.9(22.1yn0.5yn10.1yn20.24n0.12)(n2,,9)

利用精确解yexx1求出起步值（一般可用四阶RK公式计算起步值）后，按上面的公式计算，结果见下表

表9.4

Adams显式法 Adams隐式法 y(xn)yn y(xn)yn 0 0 2.1107 2.87106 3.8107 4.82106 5.2107 6.77106 6.3107 8.09106 7.1107 9.19106 7.7107 9.95106 8.1107 1.052105 8.4107 从表9.4可以看出，四阶Adams隐式法比显式法的精度高，比较这两种方法的局部截断误差公式(9-45)与(9-47)，可以说明这一现象。一般地，同阶的隐式法比显式法精确，而且数值稳定性也好。但在隐式公式中，通常很难用迭代法求解，这样又增加了计算量。因此实际计算时，很少单独使用显式公式或隐式公式，而是将它们联合使用：先用显式公式求出y(xn1)的预测值，记作，再用隐式公式对预测值进行校正，求出y(xn1)的近似值。

4.3 预测-校正系统

一般地，采用同阶的显式公式与隐式公式配对使用，即把由显式求出的yn1（记yn1）作为y(xn1)的预测值，然后再代入隐式公式进行校正，求出更接近y(xn1)的值yn1。这样就构成了预测-校正系统。采用的预测校正系统有两种：

Adams显式-隐式公式

yn1ynh(55fn59fn137fn29fn3)24 (9-52) yn1ynh24(9f(xn1,yn1)19fn5fn1fn2)Milne-Hamming预测校正系统

yn1yn34h(2fnfn132fn2)yn118(9ynyn2)3 (9-53) 8h(f(xn1,yn1)2fnfn1)说明：(1) 以上两种预测校正公式均为四阶公式，其起步值通常用四阶RK公式计算。

(2) 有时为提高精度，校正公式可迭代进行多次，但迭代次数一般超过3次。为减少一次迭代所产生的误差，常常用局部截断误差进一步修正预测值与校正值，得到更精确的预测-校正公式。下面分别对Adams预测-校正公式与Milne-Hamming公式进行讨论。

由Adams公式的局部截断误差公式

y(x251n1)yn1720h5y(5)nO(h6) (9-54) y(x195n1)yn1720hy(5)nO(h6) (9-55)

两式相减得

y270n1yn1720h5y(5)nO(h6) 从而有

h5y(5)720n270(yn1yn1) 将上式代人(9-54)，(9-55)，得

y(x251n1)yn1720(yn1yn1) (9-56) y(xn1)yn119720(yn1yn1) (9-57)

用式(9-56)，(9-57)修正公式(9-52)，就得到多环节的Adams预测-校正公式

pn1ynh(55fn59fn137fn29fn3) 预测24mn1pn1251(cnpn) 改进270cn1yh (9-58)

n(9f(xn1,yn1)19fn5fn1fn2) 校正24y19n1cn1270(cn1pn1) 改进第二式本该m251n1pn1270(cn1pn1)，但因计算的公式在后，只好用cnpn近似替代。 完全类似地，可以导出多环节的Milne-Hamming预测-校正公式（也称为Hamming方法）

pn1yn34h(2fnfn12fn2) 预测3mn1pn1112(cp) 改进121nncn118(9ynyn2)3 (9-59)

8hf(xn1,mn1)2fnfn1 校正y9n1cn1121(cn1pn1) 改进上述预测-校正公式的优点是每算一步只需计算两个函数值，计算量小于四阶RK方法，而且在计算过程中已大致估计出误差。它们不足之处在于必须利用别的方法进行启动。

(1) 输入a,b,f(x,y),N,y0 (2) 置hbaN,x0a,n1 (3) 计算

fn1f(xn1,yn1)K1hfn1Khfxh,yK12n12n12KxhK3hfn12,y2n12

K4hfxn1h,yn1K3yy1nn16(K12K22K3K4)xnanh输出(xn,yn)

(4) 若n3，置n1n，返回3；否则，置n1n,0p0,0c0，转6。 (5) 计算

f3f(x3,y3)xx3hpy403(2f3f22f0)mp112121(c0p0) c18(9y33y1)8hf(x,m)2f3f2yc9121(cp)输出(x,y)

(6)若nN，置n1n,xj1xj,yj1yj,fj1fj (j01,2),xx3,yy3,pp0,

cc0，转6；否则停机。

5 一阶微分方程组与高价微分方程的数值解法

6.1 一阶微分方程组的数值解法 设有一阶微分方程组的初值问题

yifi(x,y1,y2,,ym)y (i1,,m)(9-60) i(a)yi0若记y(yT,yT1,y2,,ym)0(y10,y20,,ym0),f(fT1,f2,,fm)，则初值问题(9-60)可写成如下向量形式

yf(x,y)y(a)y(9-61) 0如果向量函数f(x,y)在区域：axb,yRm连续，且关于满足Lipschitz条件，即存在L0，使得对x[a,b],y1,y2Rm，都有

f(x,y1)f(x,y2)Ly1y2

那么问题(9-61)在[a,b]是存在唯一解y=y(x)。

问题(9-61)与问题(9-1)形式上完全相同，故对初值问题(9-61)所建立的各种数值解法可全部用于求解问题(9-61)，只需将换成向量，f(x,y)换成f(x,y)即可。

例如，对问题(9-61)，Euler方法的计算公式为

yn1ynhf(xn,yn)(9-62)

其中yTn(y1n,y2n,,ymn)，其分量形式为

yin1yinhfi(xn,y1n,,ymn) (i1,,m)

改进Euler法的计算公式为

yn1ynh(K1K2)2K1f(xn,y(9-63) n)K2f(xnh,ynhK1)经典四阶RK公式为

yn1yhn(K12K22K3K4)6K1f(xn,yn)Kf(xh,yhK2n2n21)(9-) K3f(xnh,yhnK2)22K4f(xnh,ynhK3)其分量形式为

yin1yinh(Ki12Ki22Ki3Ki4)6Ki1fi(xn,y1n,y2n,,ymn)Kf(xh,yhK,,yhKi2in21n211mn2m1)(9-65) Ki3fhhhi(xn,y1nK12,,ymnKm2)222Ki4fi(xnh,y1nhK13,,ymnhKm3)一般地，用于问题(9-61)的单步法公式为

yn1ynhΦ(xn,yn,h)

其中Φ(xn,yn,h)(1,2,,m)T是m2元的维向量函数。为便于理解，下面用含有两个方程的方程组

（m2）具体说明。

设有初值问题

yf(x,y,z),y(a)y0zg(x,y,z),z(a)z(9-66) 0它的经典四阶RK公式为

yn1ynh(K112K122K13K14)6(9-67) zn1zhn6(K212K222K23K24)K11f(xn,yn,zn)K21g(xn,yn,zn)Khhh12f(xn2,yn2K11,zn2K21)K22g(xnh,ynhK11,zhK21)22n2K13f(xhh(9-68) n,yhnK12,zK)22n222Kg(xhhh23n2,yn2K12,zn2K22)K14f(xnh,ynhK13,znhK23)K24g(xnh,ynhK13,znhK23)计算过程为从节点处的近似解yn,zn出发，按式(9-68)顺序计算K11,K21,K12,K22,K13,K23,K14,K24，将所得结果代人式(9-67)，即得xn1处的近似解。

6.2高价微分方程的数值解法

高价微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有阶常微分方程初值问题

y(m)f(x,y,y,,y(m1)) axby(a)y(1)(m1)0,y(a)y0,,y(a)y(m1)(9-69) 0引入新变量y1y,y2y,,ymy(m1)，问题(9-69)就化为一阶微分方程组初值问题

y1y2 y1(a)y0y2yy(1)3 2(a)y0 (9-70) ym1y ym1(a)y(m2)m0ymf(x,y(m1)1,y2,,ym) ym(a)y0例7 用四阶RK方法求解初值问题

y5e2xsinx2y2y 0x1y(0)2,y(0)3(9-71) 解 引入新变量zy，问题(9-71)化为如下形式 yz, y(0)2z5e2xsinx2y2z, z(0)3 按式(9-69)，(9-70)求解。因为x00,y02,z03，由式(9-70)得

K11z03K215e2x0sinx02y02z0462K12zh02K2130.13.1K225ehsinhhh22y02K112z02K211.6238224K13zh02K223.0811911K235e0.1sin0.05220.123.1230.121.62382241.5762046K14z0hK223.1576205K245e0.2sin0.1220.13.0811911230.11.57620461.03186将上述结果代人式(9-69)得

y(xn)yn 0 -2 -3 0 1.85106 106 106 105 105 105 105

0.7673908 105 105 105