Transcript 第九章

第9章
常微分方程初值问题数值解法
9.1
9.2
9.3
9.4
9.5
9.6
1
9.7
引言
简单的数值方法
龙格-库塔方法
单步法的收敛性与稳定性
线性多步法
线性多步法的收敛性与稳定性
一阶方程组与刚性方程
9.1
引
言
考虑一阶常微分方程的初值问题
y  f ( x, y ),
x  [ x0 , b],
y ( x0 )  y0 .
（1.1）
（1.2）
如果存在实数 L  0 ，使得
f ( x, y1 )  f ( x, y2 )  L y1  y2 . y1 , y2  R,
（1.3）
则称 f 关于 y 满足利普希茨(Lipschitz）条件， L 称为 f 的
利普希茨常数（简称Lips.常数）.
2
定理1
设 f 在区域 D  {( x, y) a  x  b, y  R}上连
续，关于 y 满足利普希茨条件，则对任意 x0  [a, b], y0  R，
常微分方程(1.1),(1.2)式当 x  [a, b] 时存在唯一的连续
可微解 y (x) .
关于解对扰动的敏感性，有以下结论.
定理2
设 f 在区域 D（如定理1所定义）上连续，且
关于 y 满足利普希茨条件，设初值问题
y  f ( x, y),
y( x0 )  s
的解为 y ( x, s )，则
y ( x, s1 )  y ( x, s2 )  e
L x  x0
s1  s2 .
3
这个定理表明解对初值的敏感性，它与右端函数 f 有
关，当 f 的Lips.常数 L 比较小时，解对初值和右端函数相
对不敏感，可视为好条件.若 L 较大则可认为坏条件，即为
病态问题.
如果右端函数可导，由中值定理有
f ( x,  )
y ( x, y1 )  y ( x, y2 ) 
y1  y2 , 在y1 , y2之间.
y
若假定 f ( x, y ) 在域 D内有界，设 f ( x, y )  L，则
y
y( x, y1 )  y( x, y2 )  L y1  y2 .
y
4
它表明 f 满足利普希茨条件，且 L的大小反映了右端函
数 f 关于 y变化的快慢，刻画了初值问题(1.1)和(1.2)式是否
是好条件.
求解常微分方程的解析方法只能用来求解一些特殊类
型的方程，实际中归结出来的微分方程主要靠数值解法.
5
所谓数值解法，就是寻求解 y (x) 在一系列离散节点
x1  x2    xn  xn1  
上的近似值 y1 , y2 ,, yn , yn 1 ,.
相邻两个节点的间距 hn  xn 1  xn 称为步长.
如不特别说明，总是假定 hi  h(i  1,2,)为定数，
这时节点为 xn  x0  nh(i  0,1,2,) .
初值问题(1.1),(1.2)的数值解法的基本特点是采取
“步进式”.
即求解过程顺着节点排列的次序一步一步地向前推进.6
描述这类算法，只要给出用已知信息 yn , yn 1 , yn 2 , 
计算 yn 1 的递推公式.
首先对方程 y  f ( x, y )离散化，建立求数值解的递推
公式.
一类是计算 yn 1 时只用到前一点的值 y n，称为单步法.
另一类是用到 yn 1 前面 k 点的值 yn , yn 1 ,, yn k 1 ，
称为 k 步法.
其次，要研究公式的局部截断误差和阶，数值解 y n 与
精确解 y ( xn ) 的误差估计及收敛性，还有递推公式的计算
稳定性等问题.
7
简单的数值方法
9.2
9.2.1
欧拉法与后退欧拉法
在 xy平面上，微分方程 y  f ( x, y ) 的解 y  y (x) 称
作它的积分曲线.
积分曲线上一点 ( x, y ) 的切线斜率等于函数 f ( x, y )的
值.
如果按函数 f ( x, y )在 xy平面上建立一个方向场，那么，
积分曲线上每一点的切线方向均与方向场在该点的方向相
一致.
8
基于上述几何解释，从初始点 P0 ( x0 , y0 )出发，
先依方向场在该点的方向推进到 x  x1上一点 P1 ，然后再
从 P1 依方向场的方向推进到 x  x2 上一点 P ，循此前进
2
做出一条折线 P0 P1 P2 （图9-1).
图9-1
9
一般地，设已做出该折线的顶点 Pn ，过 Pn ( xn , yn )依
方向场的方向再推进到 Pn 1 ( xn 1 , yn 1 )，显然两个顶点
Pn , Pn 1 的坐标有关系
yn 1  yn
 f ( xn , yn ),
xn 1  xn
即
yn 1  yn  hf ( xn , yn ).
（2.1）
这就是著名的欧拉（Euler）方法.
欧拉法实际上是对常微分方程中的导数用均差近似，即
y ( xn 1 )  y ( xn )
 y( xn )  f ( xn , y ( xn ))
h
直接得到的.
10
若初值 y 0 已知，则依公式（2.1）可逐步算出
y1  y0  hf ( x0 , y0 ),
y2  y1  hf ( x1 , y1 ),

11
例1
求解初值问题
2x


y  y 
y


 y (0)  1.
(0  x  1),
表9  1
解
（2.2）
计算结果对比
欧拉公式的具体形式为
x
y
n
xn
n
yn
0.1
1.1000 2 xn 0.6
1.5090
0.2
1.1918 yn
0.7
1.5803
0 .3
1.2774
0.8
1.6498
0.4
1.3582
0.9
1.7178
yn 1  yn  h( yn 
取步长 h  0.1，计算结果见表9-1.
).
0.5 1.4351
1.7848
初值问题（2.2）的解为
y  1 21x.0，按这个解析式
子算出的准确值 y ( xn )同近似值 y n一起列在表9-1中，两者
相比较可以看出欧拉方法的精度很差.
12
表9  1
计算结果对比
xn
yn
y ( xn )
xn
yn
y ( xn )
0.1
1.1000
1.0954
0.6
1.5090
1.4832
0.2
1.1918
1.1832
0.7
1.5803
1.5492
0.3
1.2774
1.2649
0.8
1.6498
1.6125
0.4
1.3582
1.3416
0.9
1.7178
1.6733
0.5
1.4351
1.4142
1.0
1.7848
1.7321
还可以通过几何直观来考察欧拉方法的精度.
假设 yn  y ( xn ) ，即顶点 Pn 落在积分曲线 y  y (x)上，
那么，按欧拉方法做出的折线
过点 Pn 的切线（图9-2）.
Pn便是
Pn 1
y  y (x)
13
图9-2
从图形上看，这样定出的顶点 Pn 1 显著地偏离了原来
的积分曲线，可见欧拉方法是相当粗糙的.
为了分析计算公式的精度，通常可用泰勒展开将 y ( xn1 )
在 xn 处展开，则有
14
y ( xn1 )  y ( xn  h)
h2
 y ( xn )  y( xn )h 
y( n )
2
 n ( xn , xn1 ).
在 yn  y ( xn ) 的前提下， f ( xn , yn )  f ( xn , y( xn ))  y( xn )
于是可得欧拉法（2.1）的误差
h2
h2
y ( xn1 )  yn1 
y( n )
y( xn ),
2
2
（2.3）
称为此方法的局部截断误差.
yn 1  yn  hf ( xn , yn ). （2.1）
如果对方程 y  f ( x, y )从 xn 到 xn 1 积分，得
y ( xn 1 )  y ( xn )  
xn1
xn
f (t , y (t )) dt.
（2.4）
15
右端积分用左矩形公式 h f ( xn , y ( xn )) 近似.
再以 y n 代替 y ( xn ), y n1 代替 y ( xn 1 )也得到（2.1），
局部截断误差也是（2.3）.
如果在（2.4）中右端积分用右矩形公式 h f ( xn1 , y( xn1 ))
近似，则得另一个公式
yn 1  yn  hf ( xn 1 , yn 1 ),
（2.5）
称为后退的欧拉法.
它也可以通过利用均差近似导数 y( xn 1 )，即
y ( xn 1 )  y ( xn )
 y( xn 1 )  f ( xn 1 , y ( xn 1 ))
xn 1  xn
直接得到.
16
欧拉公式是关于 yn 1 的一个直接的计算公式，这类公
式称作是显式的；
后退欧拉公式的右端含有未知的 yn 1，它是关于 yn 1
的一个函数方程，这类公式称作是隐式的.
17
隐式方程通常用迭代法求解，而迭代过程的实质是逐
步显示化.
设用欧拉公式
yn( 01)  yn  h f ( xn , yn )
给出迭代初值 y n( 0)1 ，用它代入（2.5）式的右端，使之转化
为显式，直接计算得
yn(1)1  yn  h f ( xn1 , yn( 01) ),
然后再用 y n(1)1 代入（2.5）式，又有
yn( 21)  yn  h f ( xn1 , yn(1)1 ).
18
如此反复进行，得
yn( k11)  yn  h f ( xn1 , yn( k1) ),
(k  0,1,).
（2.6）
由于 f ( x, y ) 对 y 满足利普希茨条件(1.3).
由(2.6)减(2.5)得
yn( k11)  yn1  h f ( xn1 , yn( k1) )  f ( xn1 , yn1 )
 hL yn( k1)  yn1 .
由此可知，只要 hL  1 迭代法(2.6)就收敛到解 yn 1 .
从积分公式可以看到后退欧拉方法的公式误差与欧拉
法是相似的.
19
9.2.2
梯形方法
若用梯形求积公式近似等式（2.4）右端的积分

xn1
xn
f (t , y (t )) dt
xny
并分别用 yn , yn1代替 y( xn ),
则可得到比欧拉法
1 (x
n1 ),
（2.4）
y ( xn 1 )  y ( xn )  
f (t , y (t )) dt.
xn
精度高的计算公式
yn 1
h
 yn  [ f ( xn , yn )  f ( xn 1 , yn 1 )],
2
（2.7）
称为梯形方法.
梯形方法是隐式单步法，可用迭代法求解.
20
同后退的欧拉方法一样，仍用欧拉方法提供迭代初值，
则梯形法的迭代公式为
 yn( 01)  yn  h f ( xn , yn );


 yn( k11)  yn  h [ f ( xn , yn )  f ( xn1 , yn( k1) )]

2
（2.8）


(k  0,1,2,).
为了分析迭代过程的收敛性，将（2.7）与（2.8）式相
减，得
h
(k )
h
y

y

[
f
(
x
,
y
)

f
(
x
,
y
（2.7）
n 1
n 1 )],
yn 1n 1 yn  [ f (2xn , yn )n 1 f (n x1n 1 , yn 1 )],
2
( k 1)
n 1
21
于是有
hL
yn 1  yn( k1) ,
2
式中 L 为 f ( x, y )关于 y 的利普希茨常数.
yn 1  yn( k11) 
如果选取 h 充分小，使得
hL
 1,
2
则当 k  时有 yn( k1)  yn 1，
这说明迭代过程(2.8)是收敛的.
22
9.2.3
改进欧拉公式
梯形方法虽然提高了精度，但其算法复杂.
在应用迭代公式（2.8）进行实际计算时，每迭代一
次，都要重新计算函数 f ( x, y )的值.
而迭代又要反复进行若干次，计算量很大，而且往往
难以预测.
为了控制计算量，通常只迭代一两次就转入下一步的
计算，这就简化了算法.
具体地，先用欧拉公式求得一个初步的近似值 yn 1 ,
称之为预测值，
23
预测值 yn 1的精度可能很差，再用梯形公式（2.7）将
它校正一次，即按（2.8）式迭代一次得 y n1，这个结果称
校正值.
这样建立的预测-校正系统通常称为改进的欧拉公式：
预测
yn1  yn  h f ( xn , yn ),
h
（2.7）
h
y

y

[
f
(
x
,
y
)

f
(
x
,
y
)],
n

1
n
n
n
n

1
n

1
校正
yn 1 
2 yn  [ f ( xn , yn )  f ( xn 1 , yn 1 )]. （2.9）
h 2
( k 1)
（2.8）
yn1  yn  [ f ( xn , yn )  f ( xn1 , yn( k1) )]
2
也可以表为下列平均化形式
yc  yn  h f ( xn1 , y p ),
y p  yn  h f ( xn , yn ),
1
yn1  ( y p  yc ).
2
24
例2
用改进的欧拉方法求解初值问题（2.2）.
2x


y  y 
y


 y (0)  1.
解
(0  x  1),
（2.2）
这里 f ( x, y )  ( y  2 x ), 改进的欧拉公式为
y

2 xn
y

y

h
(
y

),
n
n
 p
yn


 y  y  h( y  2 xn1 ),
c
n
p

yp

1

y

( y p  yc ).
 n1
2
25
仍取 h  0.1，计算结果见表9-2.
表92 计算结果对比
xn
yn
y ( xn )
xn
yn
y ( xn )
1.0954
0.6
1.4860
1.4832
0.2
1.1841 1.1832
0.7
1.5525
1.5492
0.3
1.2662
1.2649
0.8
1.6153
1.6125
0.4
1.3434
1.3416
0.9
1.6782
1.6733
0.5
1.4164
1.4142
1.0
1.7379
1.7321
0.1 1.0959
同例1中欧拉法的计算结果比较，改进欧拉法明显改
善了精度.
26
9.2.4
单步法的局部截断误差与阶
初值问题(1.1)，(1.2)的单步法可用一般形式表示为
yn1  yn  h( xn , yn , yn1 ,（1.1）
h),
y

f
(
x
,
y
),


 y0 .
其中多元函数  与 f ( xy (, xy0))有关，
（1.2）
（2.10）
当  含有 yn 1 时，方法是隐式的，若不含 yn 1 则为
显式方法，所以显式单步法可表示为
（2.1）
2yn 1  yn  2hf ( xn , yn ).
h
h
（2.11）
 yynn1h (yxn(,yn n)
, h), y( xn ), （2.3）
y ( xynn11) 
2
2
 ( x, y, h) 称为增量函数，例如对欧拉法（2.1）有
 ( x, y, h)  f ( x, y ). 它的局部截断误差已由(2.3)给出. 27
对一般显式单步法则可如下定义.
定义1 设 y  y (x) 是初值问题(1.1)，(1.2)的准确解，
称
Tn 1
  f ( x, y ), （1.1） （2.12）
y

 y( xn 1 )  y ( xn ) h ( xn , y( xn ), h)
（1.2）
 y ( x0 )  y0 .
为显式单步法（2.11）的局部截断误差.
Tn 1 之所以称为局部的，是假设在 xn 前各步没有误差.
当 yn  yy( xn 
) 时，计算一步，则有
yn  h ( xn , yn , h),
n 1
（2.11）
y ( xn1 )  yn1  y ( xn1 ) [ yn  h ( xn , yn ,h)]
 y( xn1 )  y( xn )  h ( xn , y( xn ), h) Tn281.
局部截断误差可理解为用方法（2.11）计算一步的误差，即
在前一步精确的情况下用公式（2.11）计算产生的公式误差.
根据定义，欧拉法的局部截断误差
（2.11）
, yn(,xh),
Tn1 yny(1 xn1y)nyh(
xn()xn hf
,
y
(
x
n
n ))
 y( xn  h)  y( xn )  hy( xn )
h2

y( xn )  O(h 3 ),
2
即为（2.3）的结果.
h2
2
这里
y( xhn2) 称为局部截断误差主项.
显然
h
y ( xn1 ) 2yn1 
y( n )
y( xn ), （2.3）
2 Tn 1  O(2h 2 )
29
定义2
设 y (x) 是初值问题(1.1)，(1.2)的准确解，
若存在最大整数 p使显式单步法(2.11)的局部截断误差满足
 y  f ( x, y ), （1.1） p 1
Tn 1  y ( x  h)  y( x)  h ( x, y, h)  O(h ),
（1.2）
 y ( x0 )  y0 .
（2.13）
则称方法（2.11）具有
yn 1  yn ph阶精度.
 ( xn , yn , h),
（2.11）
若将（2.13）展开式写成
（2.10）
yn1  yn  h ( xn , yn , yn1p,1h),
Tn 1   ( xn , y ( xn )) h
 O(h p  2 ),
则  ( xn , y ( xn )) h p 1 称为局部截断误差主项.
以上定义对隐式单步法（2.10）也是适用的.
30
对后退欧拉法（2.5）其局部截断误差为
Tn1  y( xn1 )  y( xn )  hf ( xn1 , y( xn1 ))
h2
 hy( xn ) 
y( xn )  O ( h 3 )
2
 h[ y( xn )  hy( xn )  O ( h 2 )]
2
h
（2.5）
yn 
hf
1 
n O
1 ,(y
1 ),
 yn y
( xn()x
hn3).
2
2
h
这里 p  1，是1阶方法，局部截断误差主项为 
y( xn ) .
2
31
对梯形法（2.7）有
h
Tn1  y ( xn1 )  y ( xn )  [ y( xn )  y( xn1 )]
2
h2
h3
 hy( xh
y( xn ) 
y( xn )
n )
（2.7）
yn 1  yn  [ f (2xn , yn )  f 3
(!xn 1 , yn 1 )],
2
h
h2
 [ y( xn )  y( xn )  hy( xn ) 
y( xn )]  O(h 4 )
2
2
h3

y( xn )  O(h 4 ).
12
h3
所以梯形方法是二阶的，其局部误差主项为 
y( xn )
12
32
龙格-库塔方法
9.3
33
9.3.1
显式龙格-库塔法的一般形式
上节给出了显式单步法的表达式
yn 1  yn  h ( xn , yn , h),
其局部截断误差为
Tn 1  y ( x  h)  y ( x)  h ( x, y, h)  O(h p 1 ),
对欧拉法 Tn 1  O(h 2 )，即方法为 p  1 阶，
若用改进欧拉法，它可表为
yn 1
h
 yn  [ f ( xn , yn )  f ( xn  h, yn  hf ( xn , yn ))].
2
34
（3.1）
此时增量函数
1
 ( xn , yn , h)  [ f ( xn , yn )  f ( xn  h, yn  hf ( xn , yn ))].
2
（3.2）
与欧拉法的  ( xn , yn , h)  f ( xn , yn )相比，增加了计算一个
右函数 f 的值，可望 p  2 .
若要使得到的公式阶数 p更大， 就必须包含更多的
f 值.
从方程 y  f ( x, y ) 等价的积分形式（2.4），即
y ( xn 1 )  y ( xn )  
xn1
xn
f ( x, y ( x)) dx,
（3.3）
35
若要使公式阶数提高，就必须使右端积分的数值求积公式
精度提高，必然要增加求积节点.
为此可将（3.3）的右端用求积公式表示为

xn1
xn
r
f ( x, y ( x)) dx  h ci f ( xn  i h, y ( xn  i h)).
i 1
点数 r 越多，精度越高，上式右端相当于增量函数  ( x, y, h)，
x
y ( xn 1 )  y ( xn )  
f ( x, y ( x)) dx, （3.3）
x
为得到便于计算的显式方法，可类似于改进欧拉法，将公
n 1
n
式表示为
yn 1  yn  h ( xn , yn , h),
其中
（3.4）
36
r
 ( xn , yn , h)  ci K i ,
i 1
K1  f ( xn , yn ),
i1
K i  f ( xn  i h, yn  h ij K j )
i  2,, r , （3.5）
j 1
这里 ci , i , yij均为常数.
)n 1  yn  hf ( xn , yn ).
n1  yn  h ( xn , yn , hy
1
(3.4)和(3.5)称为
 ( xn , yn , h)  [ f ( xnr, 级显式龙格-库塔（Runge-Kutta）法，
yn )  f ( xn  h, yn  hf ( xn , yn ))].
2
简称R-K方法.
当 r 1,  ( xn , yn ,h)  f ( xn , yn ) 时，就是欧拉法，
（3.4）
y

y

h

(
x
,
y
,
h
),
n

1
n
n
n
此时方法的阶为 p  1.
37
当 r  2时，改进欧拉法(3.1)，(3.2)也是其中的一种，
下面将证明阶 p  2 .
要使公式(3.4),(3.5)具有更高的阶 p,就要增加点数 r .
下面就 r  2推导R-K方法.
yn 1  yn  h ( xn , yn , h),
（3.4）
r
 ( xn , yn , h)  ci K i ,
i 1
K1  f ( xn , yn ),
i  2,, r , （3.5）
i1
K i  f ( xn  i h, yn  h ij K j )
j 1
38
9.3.2
二阶显式R-K方法
对 r  2 的R-K方法，计算公式如下
 yn 1  yn  h(c1 K1  c2 K 2 ),

 K1  f ( xn , yn ),
 K  f ( x   h, y   hK ).
n
2
n
21
1
 2
（3.6）
这里 c1 , c2 , 2 , 21均为待定常数.
希望适当选取这些系数，使公式阶数 p 尽量高.
根据局部截断误差定义，（3.6）的局部截断误差为
Tn1  y ( xn1 )  y ( xn )
（3.7）
 h[c1 f ( xn , yn )  c2 f ( xn  2 h, yn  21hf n )],
39
这里 yn  y( xn ), f n  f ( xn , yn ) .
为得到 Tn 1的阶 p，要将上式各项在 ( xn , yn ) 处做泰
勒展开，由于 f ( x, y )是二元函数，故要用到二元泰勒展开，
各项展开式为
h2
h3
y ( xn 1 )  yn  hyn 
yn 
yn  O(h 4 ),
2
3!
其中
yn  f ( xn , yn )  f n ,
yn 
d
f ( xn , y ( xn ))  f x( xn , yn )  f y ( xn , yn ) f n ,
dx
 ( xn , yn )
yn f xx ( xn , yn )  2 f n f xy ( xn , yn )  f n2 f yy
 f y ( xn , yn )[ f x( xn , yn )  f n f y ( xn , yn )];
40
（3.8）
f ( xn  2 h, yn   21h f n )
 f n  f x( xn , yn )2 h  f y ( xn , yn ) 21h f n  O(h 2 ).
将以上结果代入局部截断误差公式则有
h2
Tn1  h f n  [ f x( xn , yn )  f y ( xn , yn ) f n ]
2
 h[c1 f n  c2 ( f n  2 f x( xn , yn ) h
3



f
(
x
,
y
)
f
h
]

O
(
h
)
y

y

h
(
c
K

c
K
),
 n 1
n 21 y
1n 1 n
2n
2
（3.6）

1
( xn , yn ),
 K1 f(1
 c1  c2 ) f n h  (  c2 2 ) f x( xn , yn ) h 2
 K  f ( x   h, y  2 hK ).
n
2
n
21
1
 2
1
 (  c2  21 ) f y ( xn , yn ) f n h 2  O ( h 3 ).
2
要使公式（3.6）具有 p  2 阶，必须使
41
即
1
1

c


0
,
 c2  21  0.
1 c1  c2  0,
2 2
2
2
1
1
c2 2  , c2  21  , c1  c2 1.
2
2
（3.9）
非线性方程组（3.9）的解是不唯一的.
h 0，则得
yn 1令 cy2n 
a 
[ f ( xn , yn )  f ( xn  h, yn  hf ( xn , yn ))].
2
1
c1  1  a, 2   21 
.
2a
这样得到的公式称为二阶R-K方法，如取 a  1 / 2，则
c1  c2  1 / 2,
2  21  1.
这就是改进欧拉法（3.1）.
42
若取 a  1，则 c1  0, c2  1, 2  21  1 / 2, 得计算公式 .
yn1  yn  hK 2 ,
K1  f ( xn , yn ),
（3.10）
h
h
K 2  f ( xn  , yn  K1 ).
2
2
 yn 1  yn  h(c1 K1  c2 K 2 ),
称为中点公式，相当于数值积分的中矩形公式.
（3.6）

 K1  f ( xn , yn ),
（3.10）也可表示为
 K  f ( x   h, y   hK ).
n
2
n
21
1
 2
h h h h
yn1y y

h
f
(
x

f ( xfn ,(yxn )).
n y  hf n( x , y n, y 
n1
n
n2
n
n , y n )).
2
2
2
4
r  2的R-K公式(3.6)的局部误差不可能提高到 O(h 43).
把 K 2多展开一项，从（3.8）的 yn
f y f x  f y f 的项是不能通过选择参数消掉的.
实际上要使 h 3 的项为零，需增加3个方程，要确定4个
参数 c1 , c2 , 2 , 21，这是不可能的.
故 r  2的显式R-K方法的阶只能是 p  2 ，而不能得
到三阶公式.
44
9.3.3
三阶与四阶显式R-K方法
要得到三阶显式R-K方法，必须 r  3.
此时（3.4）,（3.5）的公式表示为
 yn 1  yn  h(c1 K1  c2 K 2  c3 K 3 ),
 K  f ( x , y ),
 1
（3.11）
n
n
 yn 1  yn  h ( xn , yn , h),
（3.4）
K

f
(
x


h
,
y


hK
),
n
2
n
21
1
 2
r

K

f
(
x


h
 2 ( x , y ,nh)  3 , cynK,31hK1  32 hK 2 ).
n
n

i
i
其中 c1 , c2 , c3及 2 , 21i,13 , 31, 32 均为待定参数.
K1  f ( xn , yn ),
i  2,, r , （3.5）
公式（3.11）的局部截断误差为
i1
K i  f ( xn  i h, yn  h ij K j )
Tn 1  y( xn 1 )  y( xn )  h[c1jK
1 1  c2 K 2  c3 K 3 ].
45
只要将 K 2 , K 3 按二元函数泰勒展开，使 Tn 1  O(h 4 )，
可得待定参数满足方程

c1  c2  c3  1,

2   21 ,

     ,
32
31
 3

1
c2 2  c33  ,
2

1

2
2
c2 2  c33  3 ,

1
c32  32  .
6

（3.12）
46
这是8个未知数6个方程的方程组，解也不是唯一的.
所以这是一簇公式.
满足条件（3.12）的公式（3.11）统称为三阶R-K公式.
一个常见的公式为
h

 yn 1  yn  6 ( K1  4 K 2  K 3 ),

 K1  f ( xn , yn ),

h
h
 K 2  f ( xn  , yn  K1 ),
2
2


 K 3  f ( xn  h, yn  hK1  2hK 2 ).
此公式称为库塔三阶方法.
47
继续上述过程，经过较复杂的数学演算，可以导出各
种四阶龙格-库塔公式，下列经典公式是其中常用的一个：
h

 yn 1  yn  6 ( K1  2 K 2  2 K 3  K 4 ),

 K1  f ( xn , yn ),

h
h

 K 2  f ( xn  , yn  K1 ),
（3.13）
2
2

 K  f ( x  h , y  h K ),
n
n
2
 3
2
2

 K 4  f ( xn  h, yn  hK 3 ).

四阶龙格-库塔方法的每一步需要计算四次函数值 f ，
可以证明其截断误差为 O(h 5 ).
48
例3
设取步长 h  0.2，从 x  0到 x  1用四阶龙格-
库塔方法求解初值问题
2x


y  y 
y


 y (0)  1.
解
(0  x  1),
这里 f ( x, y )  y  2 x ，经典的四阶龙格-库塔公式为
y
h
yn1  yn  ( K1  2 K 2  2 K 3  K 4 ),
6
2 xn
K1  y n 
,
yn
49
2 xn  h
h
K 2  yn  K1 
,
h
2
yn  K1
2
表9  3
计算结果
xn
yn
y ( xn )
2 xn  h
h
K 3  yn  K 2 
,
h
2
yn  K 2
2
0.2
1.1832
1.1832
0.4
1.3417
1.3416
0.6
1.4833
1.4832
2( xn  h)
K 4  yn  hK 3 
.
yn  hK 3
0.8
1.6125
1.6125
1 .0
1.7321
1.7321
表9-3列出了计算结果，同时了出了相应的精确解.
比较例3和例2的计算结果，显然以龙格-库塔方法的精
度为高.
50
四阶龙格-库塔方法每一步要4次计算函数 f ,其计算量
比每一步只要2次计算函数 f 的二阶龙格-库塔方法-----改进的欧拉方法大一倍.
但由于这里放大了步长 (h  0.2) ，所以表9-3和表9-2
所耗费的计算量几乎相同.
龙格-库塔方法的推导基于泰勒展开方法，因而它要求
所求的解具有较好的光滑性质.
反之，如果解的光滑性差，那么，使用四阶龙格-库塔
方法求得的数值解，其精度可能反而不如改进的欧拉方法.
51
9.3.4
变步长的龙格-库塔方法
单从每一步看，步长越小，截断误差就越小，但随着
步长的缩小，在一定求解范围内所要完成的步数就增加了.
步数的增加不但引起计算量的增大，而且可能导致舍
入误差的严重积累.
因此同积分的数值计算一样，微分方程的数值解法也
有个选择步长的问题.
在选择步长时，需要考虑两个问题：
1° 怎样衡量和检验计算结果的精度？
52
2° 如何依据所获得的精度处理步长？
考察经典的四阶龙格-库塔公式
h

 yn 1  yn  6 ( K1  2 K 2  2 K 3  K 4 ),

 K1  f ( xn , yn ),

h
h

K

f
(
x

,
y

K1 ),
 2
n
n
（3.13）
2
2

 K  f ( x  h , y  h K ),
n
n
2
 3
2
2

 K 4  f ( xn  h, yn  hK 3 ).

从节点 xn 出发，先以 h为步长求出一个近似值 y n( h)1，53
由于公式的局部截断误差为 O(h 5 )，故有
（3.14）
y ( xn 1 )  yn( h)1  ch 5 ,
然后将步长折半，即取 h 为步长从 xn 跨两步到 xn 1 ，
2
5
h
再求得一个近似值 y ( 2 ) ，每跨一步的截断误差是 c h  ，
n 1
2
因此有
y ( xn 1 )  y
h
)
2
n 1
(
5
h
 2c  ,
2
（3.15）
比较（3.14）式和（3.15）式我们看到，步长折半后，
1
误差大约减少到
，
16
54
即有
h
)
2
n 1
(h)
n 1
y ( xn 1 )  y
y ( xn 1 )  y
(
1

.
16
由此易得下列事后估计式
y ( xn 1 )  y
h
)
2
n 1
(
h
( )
1

[ yn 21  yn( h)1 ].
15
这样，可以通过检查步长，折半前后两次计算结果的偏差
 y
h
)
2
n 1
(
 yn( h)1
来判定所选的步长是否合适.
具体地说，将区分以下两种情况处理：
55
1. 对于给定的精度  ，如果    ，反复将步长折半
进行计算，直至    为止.
这时取最终得到的 y
h
)
2
n 1
(
作为结果；
2. 如果    ，反复将步长加倍，直到    为止，
这时再将步长折半一次，就得到所要的结果.
这种通过加倍或折半处理步长的方法称为变步长方法.
表面上看，为了选择步长，每一步的计算量增加了，
但总体考虑往往是合算的.
56
单步法的收敛性与稳定性
9.4
9.4.1
收敛性与相容性
数值解法的基本思想是通过某种离散化手段将微分方
程转化为差分方程，如单步法(2.11)，即
yn 1  yn  h ( xn , yn , h),
（4.1）
它在 xn 处的解为 y n ，而初值问题（1.1）,（1.2）在 xn
处的精确解为 y ( xn ) ，
 y  f ( x, y ),
记 en  y( xn )  yn 称为整体截断误差.

 y ( x0 )  y0 .
（1.1）
（1.2）
57
收敛性就是讨论当 x  xn固定且 h  xn  x0  0 时
n
en  0 的问题.
定义3
若一种数值方法对于固定的 xn  x0  nh ,
当 h  0 时有 yn  y ( xn ) ，其中 y (x)是(1.1),(1.2)的准
确解，则称该方法是收敛的.
y  f ( x, y ), （1.1）

显然数值方法收敛是指 en  y( xn )  yn  0.
（1.2）
 y ( x0 )  y0 .
对单步法（4.1）有下述收敛性定理：
yn 1  yn  h ( xn , yn , h), （4.1）
58
定理3
假设单步法（4.1）具有 p 阶精度，且增量函数
 ( x, y, h) 关于 y 满足利普希茨条件
 ( x, y, h)   ( x, y, h)  L y  y ,
yn 1  yn  h ( xn , yn , h), （4.1）
（4.2）
又设初值 y 0 是准确的，即 y0  y ( x0 )，则其整体截断误差
y ( xn )  yn  O(h p ).
证明
（4.3）
设以 yn 1 表示取 yn  y ( xn ) 用公式（4.1）求得
的结果，即
（4.4）
yn 1  y ( xn )  h ( xn , y ( xn ), h), （4.1）
yn 1  yn  h ( xn , yn , h),
则 y ( xn 1 )  yn 1为局部截断误差，
59
由于所给方法具有 p阶精度，按定义2，存在定数 C，使
y( xn1 )  yn 1  Ch p 1.
又由式（4.4）与（4.1），得
yn1  yn1  y( xn )  yn  h  ( xn , y( xn ), h)  ( xn , yn ,h) .
yn 1  yn  h ( xn , yn , h), （4.1）
利用假设条件（4.2），有
yn 1  y ( xn )  h ( xn , y ( xn ), h), （4.4）
yn1  yn1  (1  hL ) y( xn )  yn ,
从而有
（4.2）
, hy)  
(y
x, y,hy)  
Lyy
 y), y
y
( x(nx,1 )y

(
x
n1
n1
n1
n1
n1
 (1 hL ) y( xn )  yn  Ch p1.
60
即对整体截断误差 en  y( xn )  yn 成立下列递推关系式
en1  (1  hL ) en  Ch p 1 ,
（4.5）
反复递推，可得
Ch p
en  (1  hL ) e0 
[(1  hL ) n  1].
L
n
（4.6）
再注意到当 xn  x0  nh  T时
(1  hL ) n  (e
hL
)n  e
TL
,
最终得下列估计式
en  e0 e
TL
Ch p TL

(e
 1).
L
（4.7）
61
由此可以断定，如果初值是准确的，即 e0  0，则（4.3）
式成立.
依据这一定理，判断单步法（4.1）的收敛性，归结为
（4.3）
y ( xn )  yn  O(h p ).
验证增量函数  能否满足利普希茨条件（4.2）.
对于欧拉方法，由于其增量函数  就是 f ( x, y )，故当
f ( x, y ) 关于 y 满足利普希茨条件时它是收敛的.
yn 1  yn  h ( xn , yn , h), （4.1）
再考察改进的欧拉方法，其增量函数
 ( x, y, h)   ( x, y, h)  L y  y , （4.2）

1
 ( xn , yn , h)  [ f ( xn , yn )  f ( xn  h, yn  hf ( xn , yn ))].
2
给出，这时有
62
1
 ( x, y, h)  ( x, y , h)  [ f ( x, y )  f ( x, y )
2
 f ( x  h, y  h f ( x, y))  f ( x  h, y  h f ( x, y )) ]
假设 f ( x, y )关于 y 满足利普希茨条件，记利普希茨常数为 L ,
则由上式推得
 ( x, y, h)   ( x, y , h)  L(1 
h
L) y  y .
2
设 h  h0 (h0为定数），上式表明 关于 y的利普希茨常数
h0
L  L(1 
L) ,
2
因此改进的欧拉方法也是收敛的.
63
类似地，也可验证其他龙格-库塔方法的收敛性.
定理3表明 p  1 时单步法收敛，并且当 y (x)是初值
问题（1.1），（1.2）的解，（4.1）具有 p 阶精度时，
有展开式
Tn1  y( x  h)  y ( x)  h ( x, y( x), h)
 y  f ( x, y ),y( x（1.1）
) 2
 y( x)h 
h 
y
(
x
)

y
.
（1.2）
0
0

2
( x,（4.1）
y ( x),0)h ]
yn 1h[yn (x,hy((xx),n0, )yn 
, hx),
 h[ y( x)  ( x, y( x),0)]  O(h 2 ).
所以 p  1的充要条件是 y( x)   ( x, y ( x),0)  0，
64
而 y( x)  f ( x, y ( x)) ，于是可给出如下定义：
定义4
若单步法（4.1）的增量函数 满足
 ( x, y,0)  f ( x, y ),
则称单步法（4.1）与初值问题（1.1），（1.2）相容.
相容性是指数值方法逼近微分方程(1.1)，即微分方程
(1.1)离散化得到的数值方法当 h  0 时可得到 y( x)  f ( x, y).
65
定理4
p 阶方法(4.1)与初值问题(1.1)，(1.2)相容
的充分必要条件是 p  1.
由定理3可知单步法（4.1）收敛的充分必要条件是(4.1)
是相容的.
以上讨论表明 p阶方法（4.1）当 p  1 时与(1.1),(1.2)
相容，反之相容方法至少是1阶的.
66
9.4.2
定义5
绝对稳定性与绝对稳定域
若一种数值方法在节点值 y n 上大小为  的扰动，
于以后各节点值 ym (m  n)上产生的偏差均不超过  ，则
称该方法是稳定的.
以欧拉法为例考察计算稳定性.
例4
考察初值问题
 y  100 y,

 y (0)  1.
其准确解 y ( x)  e 100x 是一个按指数曲线衰减得很快的函数，
如图9-3所示.
67
用欧拉法解方程 y  100 y得
yn1  (1  100h) yn .
若取 h  0.025，则欧拉公式的具体
形式为
yn 1  1.5 yn ,
计算结果列于表9-4的第2列.
可以看到，欧拉方法的解 y n
表9  4 计算结果对比
图9-3
节点
欧拉方法
0.025
 1.5
（图9-3中用×号标出）在准确值 y0( .x050
的上下波动，计算
2.25
n)
过程明显地不稳定.
0.075
 3.375
0.100
5.0625
68
但若取 h  0.005, yn 1  0.5 yn 则计算过程稳定.
再考察后退的欧拉方法，取 h  0.025时计算公式为
1
yn 1 
yn .
3.5
计算结果列于表9-4的第3列（图9-3中标以·号），这时计算
过程是稳定的.
表9 - 4 计算结果对比
节点
欧拉方法
后退欧拉方法
0.025
 1.5
0.2857
0.050
2.25
0.0816
0.075
 3.375
0.0233
0.100
5.0625
0.0067
69
这表明稳定性不但与方法有关，也与步长 h 的大小有关，
当然也与方程中的
有关.
f (x
, y)
为了只考察数值方法本身，通常只检验将数值方法用于
解模型方程的稳定性，模型方程为
y   y ,
（4.8）
其中  为复数.
对一般方程可以通过局部线性化化为这种形式.
例如在 ( x , y )的邻域，可展开为
y  f ( x, y )
 f ( x , y )  f x( x , y )( x  x )  f y ( x , y )( y  y ) ,
70
略去高阶项，再做变换即可得到 u  u的形式.
对于 m个方程的方程组，也可线性化为 y  Ay ，
这里 A为 m m的雅可比矩阵 (
f i
).
y j
若 A有 m个特征值 1 , 2 ,, m，则 i 还可能是复数，
所以，为了使模型方程结果能推广到方程组，方程 y  y
中  应为复数.
为保证微分方程本身的稳定性，还应假定 Re(  )  0 .
先研究欧拉方法的稳定性.
模型方程 y  y 的欧拉公式为
71
yn 1  (1  h ) yn .
（4.9）
设在节点值 y n 上有一扰动值  n ，它的传播使节点值 yn 1
产生大小为  n 1 的扰动值，
假设用 yn*  yn   n 按欧拉公式得出 yn*1  yn 1   n 1
的计算过程不再有新的误差，则扰动值满足
 n 1  (1  h ) n .
可见扰动值满足原来的差分方程（4.9）.
如果差分方程的解是不增长的，即有
yn 1  yn ,
则它就是稳定的.
72
显然，为要保证差分方程（4.9）的解不增长，只要选
取 h充分小， 使 yn 1  yn ,
即
（4.10）
1  h  1.
yn 1  (1  h ) yn . （4.9）
在   h 的复平面上，这是以
yn 1  yn  h ( xn , yn , h), （4.1）
(1,0)为圆心，1为半径的单位
圆内部（图9-4）.
y   y ,
（4.8）
图9-4
这个圆域称为欧拉法的绝对稳定域，一般情形可由下
面定义.
定义6
单步法（4.1）用于解模型方程（4.8），若得
到的解 yn 1  E (h ) yn，满足 E (h )  1，则称方法
（4.1）是绝对稳定的.
73
在   h 的平面上，使 E (h )  1的变量围成的
区域，称为绝对稳定域，它与实轴的交称为绝对稳定区间.
对欧拉法，E (h ) 1 h , 其绝对稳定域由
1  h  1.
给出，绝对稳定区间为  2  h  0.
在例5中   100,  2  100h  0，即 0  h  0.02
为绝对稳定区间.
例4中取 h  0.025故它是不稳定的，当取 h  0.005时
它是稳定的.
74
用二阶R-K方法解模型方程 y  y 可得到
yn 1
故
(h ) 2
 [1  h 
] yn ,
2
(h ) 2
E (h )  1  h 
.
2
2
(
h

)
绝对稳定域由 1  h 
 1 得到.
2
绝对稳定区间为  2  h  0，即 0  h  2 / .
类似可得三阶及四阶的R-K方法的 E (h )分别为
(h ) 2 (h ) 3
E (h )  1  h 

,
2!
3!
75
(h ) 2 (h ) 3 (h ) 4
E (h )  1  h 


.
2!
3!
4!
由 E (h )  1 可得到相应的绝对稳定域.
当  为实数时则得绝对稳定区间.分别为
三阶显式R-K方法： 2.51  h  0
即 0  h  2.51/ .
四阶显式R-K方法： 2.78  h  0
即 0  h  2.78 / .
图9-5
图9-5给出了R-K方法 p  1 到 p  4的绝对稳定域.
从以上讨论可知显式的R-K方法的绝对稳定域均为有限
域，都对步长 h 有限制.
如果 h 不在所给的绝对稳定区间内，方法就不稳定.
76
例5
y  20 y (0  x 1), y (0) 1, 分别取 h  0.1，
及 h  0.2,用经典的四阶R-K方法（3.13）计算.
解
本例   20, h 分别为  2及  4 .
前者在绝对稳定区间内，后者则不在. 经典的四阶R-K方法为
h

 yn 1  yn  6 ( K1  2 K 2  2 K 3  K 4 ),

 K1  f ( xn , yn ),

h
h

K

f
(
x

,
y

K1 ),
 2
n
n
（3.13）
2
2

 K  f ( x  h , y  h K ),
n
n
2
 3
2
2

77
K

f
(
x

h
,
y

hK
).
 4
n
n
3

这里 x0  0, y0 1,
用四阶R-K方法计算其误差见下表：
表9 - 5 计算结果
xn
0.2
0.4
0.6
h  0.1 0.93 10 1 0.12 10 1 0.14 10  2
h  0.2
4.98
25.0
125.0
0.8
1.0
0.15 10 3 0.17 10  4
625.0
3125.0
以上结果看到，如果步长 h 不满足绝对稳定条件，误
差增长很快.
78
对隐式单步法，可以同样讨论方法的绝对稳定性，
例如对后退欧拉法，用它解模型方程可得
yn 1 
1
yn ,
1  h
故
1
E ( h ) 
.
1  h
由 E ( h ) 
1
可得绝对稳定域为 1  h  1.
1  h
它是以 (1,0) 为圆心，1为半径的单位圆外部，故绝对
稳定区间为    h  0 .
79
当   0 时，则 0  h   ，即对任何步长均为稳定的.
对隐式梯形法，解模型方程 y  y 得
yn 1
h
1
2 y ,

h n
1
2
故
1  h / 2
E ( h ) 
.
1  h / 2
对 Re(  )  0有 E (h )  1  h / 2  1，
1  h / 2
故绝对稳定域为   h的左半平面，
80
绝对稳定区间为    h  0 ，即 0  h  时梯形法均
是稳定的.
隐式欧拉法与梯形方法的绝对稳定域均为{h Re( h )  0},
在具体计算中步长 h 的选择只需考虑计算精度及迭代收敛
性要求而不必考虑稳定性，具有这种特点的方法特别重要.
定义7
如果数值方法的绝对稳定域包含了{h Re( h )  0},
那么称此方法是A-稳定的.
由定义知A-稳定方法对步长 h 没有限制.
81
9.5
线性多步法
在逐步推进的求解过程中，计算 yn 1之前事实上已经
求出了一系列的近似值 y0 , y1 , , yn，
如果充分利用前面多步的信息来预测 yn 1，则可以期
望会获得较高的精度.
这就是构造所谓线性多步法的基本思想.
构造多步法的主要途径是基于数值积分方法和基于泰
勒展开方法，前者可直接由方程 y  f ( x, y ) 两端积分后
利用插值求积公式得到.
本节主要介绍基于泰勒展开的构造方法.
82
9.5.1
线性多步法的一般公式
如果计算 yn  k 时，除了使用 yn k 1 的值，还用到
yni (i  0,1,2,,k  2) 的值，则称此方法为线性多步法.
一般的线性多步法公式可表示为
k 1
k
i 0
i 0
yn  k    i yn i  h  i f n i ,
（5.1）
其中 yn i 为 y ( xn i )的近似， f ni  f ( xni , yni ), xni  x0  ih,
 i ,  i 为常数. 若  0及  0不全为零，则称为线性 k 步法.
计算时需先给出前面 k个近似值 y0 , y1 ,, yk 1 ,
再由（5.1）逐次求出 yk , yk 1 ,  .
83
k 1
k
i 0
i 0
yn  k    i yn i  h  i f n i ,
（5.1）
如果  k  0，称（5.1）为显式 k 步法，
这时 yn  k 可直接由（5.1）算出；
如果  k  0 ，则（5.1）称为隐式 k 步法，
求解时与梯形法（2.7）相同，要用迭代法方可算出 yn  k .
yn 1
h
 yn  [ f ( xn , yn )  f ( xn 1 , yn 1 )],
2
（2.7）
84
k 1
k
i 0
i 0
yn  k    i yn i  h  i f n i ,
（5.1）
（5.1）中系数  i及 i可根据方法的局部截断误差及阶
确定，其定义为：
定义8
设 y (x) 是初值问题(1.1)，(1.2)的准确解，
线性多步法（5.1）在 xn  k 上的局部截断误差为
  f ( x, y ), （1.1）
y

Tnk  L[ y ( xn ); h]

（5.2）
y
(
x
)

y
.
k 1
k

1
（1.2）
0
0

 y ( xnk )   i y ( xni )  h  i y( xni ).
i0
i0
若 Tn  k  O(h p 1 )，则称方法（5.1）是 p 阶的，
p  1 则称方法（5.1）与方程（1.1）是相容的.
85
由定义8，对 Tn  k 在 xn 处做泰勒展开，由于
(ih ) 2
y ( xn  ih )  y ( xn )  ih y( xn ) 
y( xn )
2!
(ih ) 3

y( xn )  ,
3!
(ih ) 2
y( xn  ih )  y( xn )  ihy( xn ) 
y( xn )  .
2!
代入（5.2）得
k 1
k 1
i0
i 0
Tnk  y ( xnk )   i y ( xni )  h  i y( xni ).（5.2）
Tn  k  c0 y ( xn )  c1hy( xn )  c2 h 2 y( xn )
   c p h p y ( p ) ( xn )   ,
（5.3）86
c0  1  ( 0  1     k 1 ),
其中
c1  k  [1  2 2    (k  1) k 1 ]  (  0  1     k ),
1 q
[k  (1  2 q  2    (k  1) q  k 1 ]
q!
1

[ 1  2 q 1  2    k q 1 k ]
(q  1)!
q  2,3, .
cp 
（5.4）
若在公式（5.1）中选择系数  i及  i ，使它满足
c0  c1    c p  0,
k 1
c p 1  0.
k
yn  k    i yn i  h  i f n i ,
p 阶的.
由定义可知此时所构造的多步法是
i 0
i 0
（5.1）
87
且
k 1
Tn k
k
y)n(xi ) h

1 
( pi1
n2 i ,

ipf
 ycnpk1hp 
y

O
(
h
).
n
i 0
i 0
（5.1）
（5.5）
称右端第一项为局部截断误差主项，c p 1 称为误差常数.
根据相容性定义， p  1, 即 c0  c1  0，
由（5.4）得
  0 1  k 1 1,
c0 
1 k(1 0  1k     k 1 ),
（5.6）
i




k
.



i
i

c1  k  [1 i21 2  i
0 (k  1) k 1 ]  (  0  1     k ),
故方法（5.1）与微分方程（1.1）相容的充分必要条件是
（5.6）成立.
88
当 k  1时，若 1  0，则由（5.6）可求得
 0 1,  0 1.
此时公式（5.1）为
  0 1  k 1 1,

ynk11  yn kh f n ,
i i    i  k .

 
i0
即为欧拉法. i1
（5.6）
从（5.4）可求得 c2  1 / 2  0，故方法为1阶精度，
且局部截断误差为
1 q
q
q
2 (  2  3  ( k  1) 
c pT   [1
k h
1x )  O2
k 1 ]


y
(
(
h
),
n 1 q!
n
2
1

[ 1  2 q 1  2    k q 1 k ]
这和第2节给出的定义及结果是一致的.
(q  1)!
89
对 k  1,
若 1  0 ，方法为隐式公式.
为了确定系数  0 ,  0 , 1，可由 c0  c1  c2  0 解得
0  1
1,  0q  1  1 / q2.
c p  [k  (1  2  2    (k  1) q  k 1 ]
于是得到公式
q!
1 h
q 1
q 1

[


2




k
k ]
1  f
2
yn 1 (qyn
(
f
),
n 1
1)!2
即为梯形法.
由（5.4）可求得 c3  1 / 12 , 故 p  2 ，所以梯形法
是二阶方法，其局部截断误差主项是  h 3 y( xn ) / 12.
这与第9.2节中的讨论也是一致的.
90
对 k  2 的多步法公式都可利用（5.4）确定系数  i ,  i ，
并由（5.5）给出局部截断误差.
Tn k  c p 1h p 1 y ( p 1) ( xn )  O(h p  2 ).
（5.5）
91
9.5.2
阿当姆斯显式与隐式公式
考虑形如
k
yn  k  yn  k 1  h  i f n i
（5.7）
i 0
的 k 步法，称为阿当姆斯（Adams)方法.
 k  0 为显式方法，  k  0为隐式方法，通常称为阿
当姆斯显式与隐式公式，也称Adams-Bashforth公式与Adams
-Monlton公式.
这类公式可直接由对方程 y  f ( x, y )两端从 xn k 1 到
xn  k 积分求得.
92
也可以利用(5.4)由 c1    c p  0 推出，对比
k
yn  k  yn  k 1  h  i f n i
i 0
1 q
[k  (1  2 q  2    (k  1) q  k 1 ]
qk !1
k
yn  k   
1 i yn i  h
q 
1 i f n i ,
q 1
 i 0
[ 1  2




k
k ]
i 0
2
(q  1)!
可知此时系数  0 1  k 2  0,  k 1 1 .
与
cp 
显然 c0  0 成立，下面只需确定系数  0 , 1 ,,  k ，
可令 c1    ck 1  0 ，求得  0 , 1 ,,  k .
若  k  0 ，则令 c1    ck  0 来求得  0 , 1 ,,  k 1 .
93
以 k  3 为例，由 c1  c2  c3  c4  0 ，根据







 0  1   2   3 1,
2( 1  2 2  3 3 )  5,
3( 1  4 2  9 3 ) 19,
4( 1  8 2  27  3 )  65.
若  3  0 ，则由前三个方程解得
0 
5
16
23
, 1   ,  2 
,
12
12
12
得到 k  3的阿当姆斯显式公式是
h
yn  3  yn  2 
(23 f n  2  16 f n 1  5 f n ).
12
（5.8）
94
由（5.4）求得 c4  3/8  0 ，所以（5.8）是三阶方法，
局部截断误差是
3
Tn 3  h 4 y ( 4 ) ( xn )  O( h 5 ).
8
h
（5.8）
yn 若
(23 f n  2  16 f n 1  5 f n ).
3  y n  2 ，则可解得
 3  0 12
1
5
19
3
0 
, 1  
, 2 
, 3  .
24
24
24
8
于是得 k  3的阿当姆斯隐式公式为
yn 3  yn  2 
h
(9 f n 3  19 f n  2  5 f n 1  f n ),
24
它是四阶方法，局部截断误差是
（5.9）
95
Tn 3
19 5 ( 5)

h y ( xn )  O(h 6 ).
720
（5.10）
类似的方法可求得阿当姆斯显式方法和隐式方法的公式，
表9-5及表9-6分别列出了
k  1,2时的阿当姆斯显式公式
,3,4
与阿当姆斯隐式公式，其中
为步数，
为方法的阶，
p
k
c p 1为误差常数.
96
表9 - 6 阿当姆斯显式公式
k
p
公式
1
1
yn 1  yn  h f n
2
2
yn  2  yn 1 
3
3
yn 3
4
4
yn  4
h
(3 f n 1  f n )
2
h
 yn  2 
(23 f n  2  16 f n 1  5 f n )
12
h
 yn 3 
(55 f n  3  59 f n  2  37 f n 1  9 f n )
24
c p 1
1
2
5
12
3
8
251
720
97
表9  7 阿当姆斯隐式公式
k
1
公式
p
2
yn 1
2 3
yn  2
3
yn 3
4
4 5
yn  4
h
 yn  ( f n 1  f n )
2
h
 yn 1  (5 f n  2  8 f n 1  f n )
12
h
 yn  2  (9 f n 3  19 f n  2  5 f n 1  f n )
24
h
 yn 3 
(251 f n  4  646 f n 3  264 f n  2
720
 106 f n 1  19 f n )
c p 1
1

12
1

24
19

720
3

160
98
例6
用四阶阿当姆斯显式和隐式方法解初值问题
y   y  x  1, y (0)  1.
取步长 h  0.1.
解
本题 f n   yn  xn 1, xn  nh  0.1n .
从四阶阿当姆斯显式公式得到
yn4  yn3 
h
(55 f n3  59 f n2  37 f n1  9 f n )
24
1
 (18.5 yn3  5.9 yn2  3.7 yn1  0.9 yn  0.24n  3.24).
24
从四阶阿当姆斯隐式公式得到
99
h
yn3  yn2  (9 f n3 19 f n2  5 f n1  f n )
24
1
 (0.9 yn3  22.1yn2  0.5 yn1  0.1yn  0.24n  3).
24
由此可直接解出 yn 3 而不用迭代，得到
yn  3
1

(22.1yn  2  0.5 yn 1  0.1yn  0.24n  3).
24.9
计算结果见表9-8，其中显式方法中的初值 y0 , y1 , y2 , y3 及
隐式方法中的 y0 , y1 , y2 均用准确解 y ( x)  e  x  x 得到，
对一般方程，可用四阶R-K方法计算初始近似.
(表略，见教材. )
100
从以上例子看到同阶的阿当姆斯方法，隐式方法要比
显式方法误差小.
这可以从两种方法的局部截断误差主项 c p1h p1 y ( p1) ( xn )
的系数大小得到解释.
这里 c p 1分别为 251 / 720及  19 / 720.
101
c米尔尼方法与辛普森方法
0  1  ( 0  1     k 1 ),
1 q
c p  [k  (1  2 q  2  (k 1) q  k 1 ]
考虑另一个 k q
! 4的显式公式
9.5.3
q1
yn  4  yn h(13 f[n
 2q21f


]0 f n ),
3 
n  2
1 fkn 1 

1
2
k
(q 1)!
其中  0 , 1 ,  2 ,  3 为待定常数.
可根据使公式的阶尽可能高这一条件来确定其数值.
由（5.4）可知 c0  0，再令 c1  c2  c3  c4  0 得到







 0  1   2   3  4,
2( 1  2 2  3 3 ) 16,
3( 1  4 2  9 3 )  64,
4( 1 8 2  27 3 )  256.
102
解此方程组得
8
4
8
 3  ,  2   , 1  ,  0  0.
3
3
3
于是得到四步显式公式
yn  4  yn 
4h
(2 f n 3  f n  2  2 f n 1 ),
3
（5.11）
称为米尔尼(Milne）方法.
由于 c  14 / 45，故方法为4阶，其局部截断误差为
5
Tn  4
14 5 ( 5)

h y ( xn )  O(h 6 ).
45
（5.12）
米尔尼方法也可以通过方程 y  f ( x, y )两端积分
103
y ( xn  4 )  y ( xn )  
xn 4
xn
f ( x, y ( x)) dx
得到. 若将方程 y  f ( x, y ) 从 xn 到 xn  2积分，可得
y ( xn  2 )  y ( xn )  
xn 2
xn
f ( x, y ( x)) dx.
右端积分通过辛普森求积公式就有
yn  2
h
 yn  ( f n  4 f n 1  f n  2 ).
3
（5.13）
称为辛普森方法.
它是隐式二步四阶方法，其局部截断误差为
Tn  2
h 5 ( 5)

y ( xn )  O(h 6 ).
90
（5.14）
104
9.5.4
汉明方法
辛普森公式是二步方法中阶数最高的，但它的稳定性差，
为了改善稳定性，考察另一类三步法公式
yn3   0 yn  1 yn1   2 yn 2  h( 1 f n1   2 f n 2  3 f n3 ),
其中系数  0 , 1 ,  2及 1 ,  2 ,  3为常数.
若希望导出的公式是四阶的，则系数中至少有一个自由
参数.
若取 1  1，则可得到辛普森公式.
若取 1  0，仍利用泰勒展开，由（5.4），令
c0  c1  c2  c3  c4  0,
105
则可得到









 0  2 1,
2 2  1   2   3  3,
4 2  2( 1  2  2  3 3 )  9,
8 2  3( 1  4  2  9  3 )  27,
16 2  4( 1  8 2  27  3 )  81.
解此方程组得
1
9
3
6
3
 0   ,  2  , 1   ,  2  ,  3  .
8
8
8
8
8
于是有
yn 3 
1
3h
(9 yn  2  yn ) 
( f n 3  2 f n  2  f n 1 ), （5.15）
106
8
8
称为汉明(Hamming)方法.
由于 c5  1 / 40，故方法是四阶的，且局部截断误差
Tn 3
h 5 ( 5)

y ( xn )  O(h 6 ).
40
（5.16）
107
9.5.5
预测-校正方法
对于隐式的线性多步法，计算时要进行迭代，计算量
较大.
为了避免进行迭代，通常采用显式公式给出 yn  k 的一
(0)
个初始近似，记为 y n  k，称为预测(predictor).
h
y接着计算
f ( xn , yn )  f ( xn 1 , yn 1 )]. （2.13）
f n  k[的值(evaluation)，再用隐式公式计算
n 1  y n 
2
yn  k ，称为校正（corrector）.
在（2.13）中用欧拉法做预测，再用梯形法校正，得
到改进欧拉法，它就是一个二阶预测-校正方法.
108
一般情况下，预测公式与校正公式都取同阶的显式方
法与隐式方法相匹配.
例如用四阶的阿当姆斯显式方法做预测，再用四阶阿
当姆斯隐式公式做校正，得到以下格式：
p
y
预测P： n  4  yn 3 
h
(55 f n 3  59 f n  2  37 f n 1  9 f n ),
24
求值E： f np 4  f ( xn  4 , ynp 4 ),
h
校正C： yn  4  yn 3  (9 f np 4  19 f n 3  5 f n  2  f n 1 ),
24
求值E： f n  4  f ( xn  4 , yn  4 ).
此公式称为阿当姆斯四阶预测-校正格式（PECE).
109
依据四阶阿当姆斯公式的截断误差，对于PECE的预测
步P有
y ( xn  4 )  ynp 4 
251 5 ( 5)
h y ( xn ),
720
对校正步C有
y ( xn  4 )  y n  4
19 5 ( 5)

h y ( xn ),
720
两式相减得
5
h y
( 5)
720
( xn )  
( ynp 4  yn  4 ),
270
于是有下列事后误差估计
110
251 p
y ( xn  4 )  y

( yn  4  yn  4 ),
270
19
y ( xn  4 )  y n  4 
( ynp 4  yn  4 ).
270
p
n4
容易看出
251

y y 
( yn4  ynp4 ), 

270

19
y n 4  y n 4 
( yn4  ynp4 ).


270
p
y
比 n  4 , yn  4 更好.
pm
n 4
p
n 4
（5.17）
但在 ynpm
的表达式中 yn  4 是未知的，因此计算时用
4
上一步代替，从而构造一种修正预测-校正格式(PMECME):
111
h
P： y  yn 3  (55 f n 3  59 f n  2  37 f n 1  9 f n ),
24
251 c
pm
p
M： yn  4  yn  4 
( yn 3  ynp3 ),
270
p
n4
E：
pm
f npm

f
(
x
,
y
4
n4
n  4 ),
C： y
c
n4
h
 yn 3  (9 f npm
 4  19 f n  3  5 f n  2  f n 1 ),
24
M： yn  4  ync 4 
E：
19
( ync 4  ynp 4 ),
270
f n  4  f ( xn  4 , yn  4 ).
在PMECME格式中已将（5.17）的 yn  4 及 yn  4分别
改为 ync 4及 yn  4.
112
利用米尔尼公式（5.11）和汉明公式（5.15）相匹配，
并利用截断误差（5.12），（5.16）改进计算结果，可类
似地建立四阶修正米尔尼-汉明预测-校正格式PMECME）：
4
p
y

y

h(2 f n 3  f n  2  2 f n 1 ),
P： n  4
n
3
4h
（5.11）
yn  4  yn 
(pm2 f n 3 
f

2
f
),
n

2
n

1
112
c
p
 ynp54 
(
y

y
),
M：3yn  4 14
( 5)
6 n  3（5.12）
n 3
Tn  4 
h 3yh 121
( xn )  O(h ).
1
yn 3  (9 yn  2 pmy45
( fpm
n) 
n  3  2 f n  2  f n 1 ), （5.15）
5
f ( xn8 (45,)yn  4 ),
8 E： f n  4  h
Tn 3  
y ( xn )  O(h 6 ). （5.16）
40
1
3
c
pm
y

(
9
y

y
)

h
(
f
C： n  4
n 3
n 1
n  4  2 f n  3  f n  2 ),
8
8
9
c
( ync 4  ynp 4 ),
M： yn  4  yn  4 
121
E： f n  4  f ( xn  4 , yn  4 ).
113
9.5.6
构造多步法公式的注记和例
构造多步法公式有基于数值积分和泰勒展开两种途径，
k 1
k
（5.1）
y


y

h

f
,


n

k
i
n

i
i
n

i
数值积分方法只能用于可转化为等价的积分方程的微分方程，
i 0
i 0
它是有局限性的.
而用泰勒展开则可构造任意多步法公式，其做法是根
据多步法公式的形式，直接在 xn 处做泰勒展开.
确定多步法（5.1）的系数  i 及  i (i  0,1,, k )时不
必套用系数公式（5.4），因为多步法公式不一定如（5.1）
的形式，而且套用公式容易记错.
114
例7
解初值问题
y  f ( x, y ),
y ( x0 )  y0 .
用显式二步法
yn1   0 yn  1 yn1  h(  0 f n  1 f n1 ),
其中
f n  f ( xn , yn ), f n1  f ( xn1 , yn1 ).
试确定参数  0 , 1 ,  0 , 1 使方法阶数尽可能高,并求局部
截断误差.
解
根据局部截断误差定义，用泰勒展开确定参数满
足的方程.
115
由于
Tn1  y ( xn  h)  0 y ( xn ) 1 y ( xn  h)
 h [  0 y( xn )  1 y( xn  h)]
h3
h2
y ( xn )
y ( xn ) 
 y ( xn )  hy ( xn ) 
3!
2
h 4 ( 4)
y ( xn )  O ( h ( 5 ) )  0 y ( xn )

4!
h3
h2
y ( xn )
y ( xn ) 
1[ y ( xn )  hy ( xn ) 
3!
2
h 4 ( 4)
y ( xn )  O ( h ( 5 ) )]

4!
  0 h y ( xn )  1h[ y ( xn )  hy ( xn )
h3 ( 4)
h2
y ( xn )  O ( h ( 4 ) )]
y ( xn ) 

3!
2
116
 (1   0  1 ) y ( xn )  (1  1   0  1 ) hy( xn )
1 1
1 1
1
2


 (  1  1 ) h y ( xn )  (  1  1 ) h 3 y( xn )
2 2
6 6
2
1
1
1
(

1  1 )h 4 y ( 4 ) ( xn )  O (h 5 )
24 24
6
为求参数  0 , 1 ,  0 , 1 使方法阶数尽量高，可令
1 0 1  0,
1 1
 1  1  0,
2 2
11   0  1  0,
1 1
1
 1  1  0.
6 6
2
即得方程组
117









 0 1 1,
1   0  1 1,
1  21 1,
1  31 1.
解得  0  4, 1  5,  0  4, 1  2，此时公式为三阶，
而且
Tn 1 
1 4 ( 4)
h y ( xn )  O ( h 5 )
6
即为所求局部截断误差.
而所得二步法为
yn1  4 yn  5 yn1  2h(2 f n  f n1 ).
118
证明存在  的一个值，使线性多步法
例8
yn 1   ( yn  yn 1 )  yn  2
1
 (3   )h( f n  f n 1 )
2
是四阶的.
证明 只要证明局部截断误差 Tn 1  O(h 5 )，则方法
为四阶.
仍用泰勒展开，由于
Tn 1  y ( xn  h)   [ y ( xn )  y ( xn  h)]  y ( xn  2h)
1
 (3   )h [ y( xn )  y( xn  h)]
2
119
h2
h3
h 4 ( 4)
 y ( xn )  hy( xn ) 
y ( xn ) 
y ( xn ) 
y ( xn )
2
3!
4!
2
3
h
h
 O ( h ( 5) )   [(  h) y ( xn ) 
y ( xn ) 
y ( xn )
2
3!
h 4 ( 4)

y ( xn )  O ( h ( 5) )]  [ y ( xn )  2hy( xn )
4!
( 2h) 2
( 2h) 3
( 2h) 4 ( 4 )

y( xn ) 
y ( xn ) 
y ( xn )  O ( h ( 5) )]
2
3!
4!
h
h2
 (3   )[ y( xn )  y ( xn )  hy ( xn ) 
y( xn )
2
2
h3 ( 4)

y ( xn )  O ( h ( 4 ) )]
3!
120
1 1
 [1  2  (3  )]hy( xn ) [    2
2 2
1
1 1
4
 (3  )]h 2 y( xn ) [   
2
6 6
3
1
1
1
2 1
 (3  )]h 3 y( xn ) [

   (3  )]
4
24 24
3 12
 h 4 y ( 4 ) ( xn )  O ( h 5 )
3 1
1
 (   )h 3 y( xn ) 
(9  )h 4 y ( 4 ) ( xn )  O(h 5 ).
4 12
24
当   9时，Tn 1  O(h 5 )，故方法是四阶的.
121
9.6
9.6.1
线性多步法的收敛性与稳定性
相容性及收敛性
线性多步法(5.1)式的相容性：在定义8中给出的局部
截断误差(5.2)中 Tn  k  O(h p 1 )，若 p  1称 k 步法(5.1)
与微分方程(1.1)相容，它等价于
lim
h 0
1
Tn  k  0.
h
（6.1）
对多步法(5.1)可引入多项式
k 1
 ( )      j j ,
k
j 0
（6.2）
122
和
k
 ( )    j j ,
（6.3）
j 0
分别称为线性多步法(5.1)的第一特征多项式和第二特征多
项式.
可以看出，如果(5.1)式给定，则  ( )和 ( )也完全确
定.反之也成立.
123
定理5
线性多步法(5.1)式与微分方程(1.1)相容的充
分必要条件是
 (1)  0,
 (1)   (1).
（6.4）
关于多步法(5.1)的收敛性：由于用多步法求数值解需
要 k个初值，而微分方程(1.1)只给出一个初值 y( x0 )  y0，
因此还要给出 k  1个初值才能用多步法(5.1)进行求解，即
k 1
k

 yn  k    j yn  j  h  j f n  j
j 0
j 0

 y   (h), i  0,1, , k  1,
i
 i
（6.5）
其中 y0 由微分方程初值给定，y1 , y2 ,, yk 1 可由相应单步法
给出.
124
设由(6.5)式在 x  xn 处得到的数值解为 y n ，这里
xn  x0  nh  [a, b] 为固定点，h 
定义9
b  a ，于是有下面定义.
n
设初值问题(1.1),(1.2)有精确解 y (x) .如果初
始条件 yi  i (h)满足条件
lim i (h)  y0 ,
h 0
i  0,1,, k  1
的线性 k 步法(6.5)在 x  xn 处的解 y n 有
lim
h 0
x  x0  nh
yn  y ( x),
则称线性 k 步法是收敛的.
定理6
设线性多步法(6.5)是收敛的，则它是相容的.
此定理的逆定理是不成立的.
125
例9
用线性二步法
 yn  2  3 yn 1  2 yn  h( f n 1  2 f n ),

 y0  0 (h), y1  i (h),
（6.6）
解初值问题 y  2 x, y (0)  0.
解
此初值问题精确解 y ( x)  x 2，而由(6.6)式知
 ( )   2  3  2,  ( )    2,
故有  (1)  0,
 (1)   (1)  1，故方法(6.6)是相容的，但
方法(6.6)的解并不收敛，在方法(6.6)中取初值
y0  0,
y1  h,
此时方法(6.6)为二阶差分方程
（6.7）
126
yn 2  3 yn1  2 yn  h(2 xn1  4 xn ),
y0  0, y1  h, （6.8）
其特征方程为
 ( )   2  3  2  0,
解得其根为 1  1及  2  2 .于是可求得(6.8)的解为
yn  (2n  1)h  n(n  1)h 2 ),
x  nh,
2n  1
n 1 2
lim yn  lim (
x
x )  ,
h 0
n 
n
n
n 
故方法不收敛.
多步法是否收敛与  ( )的根有关.
127
定义10
如果多步法（5.1）式的第一特征多项式  ( )
的根都在单位圆内或圆上，且在单位圆上的根为单根，则称
线性多步法（5.1）满足根条件.
定理7
线性多步法是相容的，则线性多步法（6.5）收
敛的充分必要条件是线性多步法（5.1）满足根条件.
128
9.6.2
稳定性与绝对稳定性
稳定性主要研究初始条件扰动与差分方程右端项扰动对
数值解的影响，假设多步法（6.5）有扰动 { n n  0,1,, N}，
则经过扰动后的解为 {zn n  0,1,, N }, N  b  a，它满足方
程
h
k 1
k

 z n  k    j z n  j  h (   j f ( xn  j , z n  j )   n  j )
j 0
j 0

 z   (h)   , i  0,1, , k  1.
i
i
 i
（6.9）
129
定义11
对初值问题(1.1)，(1.2)，由方法（6.5）得到的
差分方程解{ yn }0N ，由于有扰动{ n }0N ，使得方程（6.9）的解
为{z n }0N ，若存在常数 C 及 h0，使对所有 h  (0, h0 ) ，当  n  
0  n  N 时，有
zn  yn  C ,
则称多步法（5.1）是稳定的或称为零稳定的.
研究零稳定性就是研究 h  0 时差分方程（6.5）解 { yn }
的稳定性.它表明当初始扰动或右端扰动不大时，解的误差
也不大.
130
对多步法（5.1），当 h  0 时对应差分方程的特征方
程为  ( )  0 .故有
定理8
线性多步法（5.1）是稳定的充分必要条件是
它满足根条件.
关于绝对稳定性只要将多步法（5.1）用于解模型方程
（4.8），得到线性差分方程
k 1
k
j 0
j 0
yn  k    j yn  j  h   j yn  j .
（6.10）
131
利用线性多步法的第一、第二特征多项式  ( ),  ( ) ，令
 ( ,  )   ( )   ( ),   h.
（6.11）
此式称为线性多步法的稳定多项式，它是关于 的 k 次多项
式.
如果它的所有零点  r  r ( )(r  1,2,, k ) 满足  r  1，
则(6.10)式的解{ yn }当 n  时，有 yn  0.
132
定义12 对于给定的   h，如果稳定多项式(6.11)的零
点满足  r  1, r  1,2,, k ，则称线性多步法(5.1)关于此  值
是绝对稳定的.若在   h 的复平面的某个区域 R 中所有
值线性多步法(5.1)都是绝对稳定的，而在区域 R 外，方法是
不稳定的，则称 R为多步法(5.1) 的绝对稳定域. R与实轴的
交集称为线性多步法(5.1) 的绝对稳定区间.
 为实数时，可以只讨论绝对稳定区间.
由于线性多步法的绝对稳定域较为复杂，通常采用根轨
迹法.
133
阿当姆斯显式方法和隐式方法的绝对稳定域图形分别为
图9-6和图9-7，绝对稳定区间见表9-9.
图9-6
表9 - 9
图9-7
阿当姆斯公式绝对稳定区间
显式方法
隐式方法
k  p  1,  2  h  0，
k  1，p  2,    h  0，
k  p  2,  1  h  0，
k  2，p  3,  6.0  h  0，
k  p  3,  0.55  h  0， k  3，p  4,  3.0  h  0，
k  p  4,  0.30  h  0， k  4，p  5,  1.8  h  0，
134
例10 讨论辛普森方法
yn  2  yn 
的稳定性
解
h
( f n  4 f n 1  f n  2 )
3
辛普森方法的第第二特征多项式为
1 2
 ( )    1,  ( )  (  4  1).
3
2
 ( )  0的根分别为-1及1，它满足根条件，故方法是零稳
定的.但它的稳定多项式为
1
3
 ( ,  )   2  1   ( 2  4  1).
求绝对稳定域 R 的边界轨迹 R .若   R，则可令   ei ，
在  平面域 R 的边界轨迹 R 为
135
 (ei )
ei 2  1
3(e i  e -i )
3i sin 
   ( ) 

 i

.
i
-i
1
 (e )
2  cos
(ei 2  4ei  1) e  4  e
3
可看出  ( )在虚轴上，且对全部   [0,2 ], 3 sin   [ 3, 3 ],
2  cos
从而可知 R 为虚轴上从  3i 到 3i 的线段，故辛普森公式
的绝对稳定域为空集.即步长 h  0 ，此方法都不是绝对稳
定的，故它不能用于求解.
136
9.7
一阶方程组与刚性方程
9.7.1
一阶方程组
前面研究了单个方程 y  f 的数值解法，只要把 y
和 f 理解为向量，那么，所提供的各种计算公式即可应用
到一阶方程组的情形.
考察一阶方程组
yi  f i ( x, y1 , y2 ,, y N )
(i  1,2,, N )
的初值问题，初始条件给为
yi ( x0 )  yi0
(i  1,2, , N ).
137
若采用向量的记号，记
y  ( y1 , y2 ,, y N )T ,
y0  ( y10 , y20 ,, y N0 )T ,
f  ( f1 , f 2 ,, f N )T .
则上述方程组的初值问题可表示为
y  f ( x, y ), 

y ( x0 )  y0 . 
求解这一初值问题的四阶龙格-库塔公式为
h
yn 1  yn  (k1  2k 2  2k3  k 4 ),
6
式中
（7.1）
138
k1  f ( xn , yn ),
h
h
k 2  f ( xn  , yn  k1 ),
2
2
h
h
k3  f ( xn  , yn  k 2 ),
2
2
k4  f ( xn  h, yn  hk3 ).
考察两个方程的特殊情形：









y  f ( x, y, z ),
z   g ( x, y, z ),
y ( x0 )  y0 ,
z ( x0 )  z0 .
139
这时四阶龙格-库塔公式具有形式





h
yn1  yn  ( K1  2 K 2  2 K 3  K 4 ),
6
h
z n1  z n  ( L1  2 L2  2 L3  L4 ).
6
（7.2）
其中
 K1  f ( xn , yn , z n ),


 K  f ( x  h , y  h K , z  h L ),
2
n
n
1
n
1

2
2
2


h
h
h
K

f
(
x

,
y

K
,
z

L2 ),

3
n
n
2
n

2
2
2
（7.3）
140











K 4  f ( xn  h, yn  hK 3 , z n  hL3 ),
L1  g ( xn , yn , zn ),
h
h
h
L2  g ( xn  , yn  K1 , z n  L1 ),
2
2
2
h
h
h
L3  g ( xn  , yn  K 2 , z n  L2 ),
2
2
2
L4  g ( xn  h, yn  hK3 , zn  hL3 ),
（7.3）
这是一步法，利用节点 xn上的值 yn , z n，由（7.3）
式顺序计算 K1 , L1 , K 2 , L2 , K 3 , L3 , K 4 , L4，然后代入（7.2）
式即可求得节点 xn 1上的 yn 1 , zn 1 .
141
9.7.2
化高阶方程为一阶方程组
高阶微分方程（或方程组）的初值问题，原则上总可
以归结为一阶方程组来求解.
例如，考察下列 m阶微分方程
y ( m)  f ( x, y, y,, y ( m1) ),
（7.4）
初始条件为
y ( x0 )  y0 , y( x0 )  y0 ,, y ( m1) ( x0 )  y0( m1) . （7.5）
只要引进新的变量
y1  y, y2  y, , ym  y ( m 1) ,
142
即可将 m阶方程（7.4）化为如下的一阶方程组：
y1  y2 ,
( m)
( m 1)


y

f
(
x
,
y
,
y
,

,
y
),
y2  y3 ,




 1  ym ,
ym


 1  ym ,
ym

  f ( x, y1 , y2 ,, ym ).
ym

（7.4）
（7.6）
初始条件（7.5）则相应地化为
y1 ( x0 )  y0 , y2 ( x0 )  y0 ,, ym ( x0 )  y0( m1) . （7.7）
初值问题（7.4）,（7.5）和（7.6）,（7.7）是彼此等价的.
y ( x0 )  y0 , y( x0 )  y0 ,, y ( m1) ( x0 )  y0( m1) . （7.5）
143
特别地，对于下列二阶方程的初值问题：
 y  f ( x, y, y),

 y ( x0 )  y0 ,
 y( x )  y .
0
0

引进新的变量 z  y，即可化为下列一阶方程组的初值问题：
 y  z,
 z   f ( x, y , z ),


 y ( x0 )  y0 ,

 z ( x0 )  y0 .
144
针对这个问题应用四阶龙格-库塔公式（7.2），有
h

yn1  yn  ( K1  2 K 2  2 K 3  K 4 ), 

6

h
z n1  z n  ( L1  2 L2  2 L3  L4 ). 


6
由（7.3）式
K1  zn , L1  f ( xn , yn , zn ),
h
h
h
h
K 2  z n  L1 , L2  f ( xn  , yn  K1 , z n  L1 ),
2
2
2
2
h
h
h
h
K 3  z n  L2 , L3  f ( xn  , yn  K 2 , z n  L2 ),
2
2
2
2
K 4  zn  hL3 , L4  f ( xn  h, yn  hK3 , zn  hL3 ).
145
如果消去 K1 , K 2 , K 3 , K 4 ，则上述格式可表示为

h2
yn1  yn  hz n 
( L1  L2  L3 ), 

6

h
z n1  z n  ( L1  2 L2  2 L3  L4 ). 


6
这里
L1  f ( xn , yn , zn ),
h
h
h
L2  f ( xn  , yn  z n , z n  L1 ),
2
2
2
h
h
h2
h
L3  f ( xn  , yn  z n 
L1 , z n  L2 ),
2
2
4
2
h2
L4  f ( xn  h, yn  hzn 
L2 , z n  hL3 ).
2
146
9.7.3
刚性方程组
在求解方程组（7.1）时，经常出现解的分量数量级差
别很大的情形，这给数值求解带来很大困难，这种问题称
为刚性(stiff)问题.
考察以下例子.
给定系统







u   10000.25u  999.75v  0.5,
v  999.75u 10000.25v  0.5,
u (0) 1,
v(0)  1.
（7.8）
147
它可用解析方法求出准确解，方程右端系数矩阵
  10000.25
A  
 999.75


 10000.25 
999.75
的特征值为
1  0.5, 2  2000,
方程的准确解为
u (t )  e 0.5t  e 2000t  1,

 0.5t
 2000t
v
(
t
)


e

e
 1.

当 t   时，u (t )  1, v(t )  1称为稳态解，
u, v 中均含有快变分量 e 2000t 及慢变分量 e 0.5t .
148
对应于 2 的快速衰减的分量在 t  0.005秒时已衰减
到 e 10  0.
1
1
称 2  
为时间常数.

 0.0005
2 2000
当 t  10 2 时快变分量即可被忽略，而对应于 1 的
1
1
慢变分量，它的时间常数  1  

2，
1 0.5
它要计算到 t  10 1  20时，才能衰减到 e 10  0 ，
也就是说解 u, v必须计算到 t  20才能达到稳态解.
它表明方程的解分量变化速度相差很大，是一个刚性
方程组.
149
若用四阶龙格-库塔法, 步长选取要满足 h  2.78 / ，
即 h  2.78 / 2  0.00139
，才能使计算稳定.
而要计算到稳态解至少需要算到 t  20，则需计算
14388步.
这种用小步长计算长区间的现象是刚性方程数值求解
出现的困难，它是系统本身病态性质引起的.
150
对一般的线性系统
dy
 Ay (t )  g (t ),
dt
其中
（7.9）
y  ( y1 ,, y N )T R N , g  ( g1 ,, g N )T R N , A  R N  N .
若 A的特征值  j   j  i j ( j  1,, N , i 
 1)
相应的特征向量为  j ( j  1, , N )，则方程组（7.9）的
解为
N
y (t )   c j e j  j  (t ),
 t
（7.10）
j 1
其中 c j 为任意常数，可由初始条件 y (a)  y 0确定， (t )
为特解.
151
假定  j 的实部  j  Re(  j )  0 ,则当 t  时，
y (t )   (t ),  (t ) 为稳态解.
dy
 Ay (t )  g (t ),
dt
（7.9）
定义13 若线性系统（7.9）中 A的特征值  j 满足条件
Re(  j )  0 ( j  1, , N ) 且
s  max Re(  j ) / min Re(  j )  1,
1 j  N
1 j  N
则称系统（7.9）为刚性方程，称 s为刚性比.
152
刚性比 s  1时， A为病态矩阵，故刚性方程也称
病态方程.
通常 s  10就认为是刚性的. s越大病态越严重.
方程（7.8）的刚性比 s  4000，故它是刚性的.
对一般非线性方程组（7.1），可类似定义13，将 f 在点
(t , y (t )) 处线性展开，记 J (t )  f  R N  N.
y
y  f (的特征值为
x, y ), 
假定
,
（7.1）
J (t )
j (t ) j  1, , N

y ( x0 )  y0 . 
于是由定义13可知，当  j (t ) 满足条件
Re(  j (t ))  0 ( j  1, , N )
153
且
s (t )  max Re(  j (t )) / min Re(  j (t ))  1,
1 j  N
1 j  N
称系统（7.1）是刚性的， s (t ) 称为方程（7.1）的局部刚性比.
求刚性方程数值解时，若用步长受限制的方法就将出
现小步长计算大区间的问题，因此最好使用对步长
h不加
y  f ( x, y ), 
（7.1）

限制的方法，如前面已介绍的欧拉后退法及梯形法，即
y ( x0 )  y0 . 
A-稳定的方法，
所谓A-稳定就是指数值方法的绝对稳定域包含了  h
平面的左半平面.
154
这种方法当然对步长 h没有限制，但A-稳定方法要求
太苛刻，Dahlquist已证明所有显式方法都不是A-稳定的，
而隐式的A-稳定多步法阶数最高为2，且以梯形法误差常
数为最小.
这就表明本章所介绍的方法中能用于解刚性方程的方
法很少.
通常求解刚性方程的高阶线性多步法是吉尔（Gear)
方法，还有隐式龙格-库塔法（见文献[21]），这些方法
都有现成的数学软件可供使用.
155
156
157
158
159
160
161
162