第3章常微分方程2

即为改进的欧拉公式（2.12）。
如取p=1/2 ，则 λ1 =0, λ2 =1 ，（3.3）就成为
（3.5）
8
（3.5）称为变形的欧拉格式。 由于（3.5）中的
是欧拉格式预测出来的区间[xi,xi+1]的中点xi+1/2的近
似解, K2 就近似地等于此中点的斜率
，
因此（3.5）就相当于用中点xi+1/2的斜率作为（3.1）
K 可看作是y = y (x) 在区间[xi,xi+1] 上的平均
斜率。这样，欧拉格式(2.2)相当于取(xi , yi) 点上斜
率
作为平均斜率K 的近似值，当然
2
是十分粗糙的，因此精度必然很低。而改进的欧拉 格式（2.12）可改写成
与（3.1）比较知，它相当于把(xi , yi) 和 (xi+1 , yi+1) 两个点上的斜率K1和K2 的算术平均值作为（3.1） 中的平均斜率K的近似值。其中K2 是通过已知信息
从节点 出发，先以h为步长跨一步到节点 ，求出一个近似值 。如计算公式是P阶的，则
22
(3.10)
当h值不大时，式中的系数C可近似地看作常数。
然后，将步长减半，即以 为步长，从节点 出发，
跨两步到节点 ，再求得一个近似值 。其中
每跨一步的截断误差
，故有
(3.11)
23
（3.11）2p （3.10）得：
(3.9)
经典龙格—库塔方法的程序框图见图7-6。
15
输入 a,b,n,a：
图7-6 stop
16
例2 试分别用欧拉方法
，改进的欧拉
方法
及经典R—K方法
，
求解下列初值问题，并比较三种方法所得结
果的精度：
解 三种方法的具体算式分别如下： 欧拉格式：
17
改进的欧拉公式： 经典R—K格式：
18
表7-2
，则与
y y f k
h
i1
i
kj i j1
j0
局部截断误差为
B 对于不同k值下的 和 kj
值可算出，分别列于
k
表7-4中。
40
从表7-4可知： 当k=0时
y y h f
i1
i
i1
即为隐式欧拉格式，其局部截断误差为：
R h 1
i1
2
2 y(
)
i1
当k=1时
y y h(f f )
i1
20
（2）经典R—K方法的局部截断误差为
为大致相同数量级下的常数，故有
注意： R—K方法的导出利用了泰勒展开，因此要求所 求的解有较好的光滑性，如果解的光滑性差，则采用经 典 R—K 法所得的数值解，其精度有可能反而不及改进 的欧拉法。因此在实际计算中，应根据问题的具体情况 来选择合适的算法。
21
还要用到更前面三步的信息
。它是一
种四步法，无法自行启动，需用其他四阶单步法
（如四阶R-K方法）先从 算出
作为其初值，
然后按（4.11）进行迭代。
阿达姆斯预测-校正方法的程序框图见图7-7。
例4 试用阿达姆斯预测-校正方法求解下列初值问
上的斜率K2 ，K3与点(xi , yi) 上的斜率K1 加权平均，
作为平均斜率K的近似值，即
10
其中K1和K2仍如（3.3）。
利用区间[xi,xi+q] 内的两个斜率K1和K2 ，加权 平均作为其平均斜率K，来预测 y(xi+q ) ：
从而得到
由此构成的计算格式为
11
(3.6)
类似于二阶龙格-库塔格式的导出过程，运用泰 勒展开的方法，可找出格式（3.6）的局部截断误差 为 ，从而具有三阶精度所必须满足的条件为：
中则是从精确解（已知
）算出。
（4.9）及 (4.10) 算出的结果分别列于表7-5中。
从表7-5可见，隐式的精度比同阶显式的要高。
47
表7-5
四步显式方法
三步隐式方法
精确解 y(xi )
48
4.3 阿达姆斯预测-校正方法 阿达姆斯隐式方法在计算过程中，一般要解超越 方程。例如下列定解问题
其阿达姆斯三步隐式方法的算式为：
的绝对值
B k
比显式的 Bk
要小;
（2）显式的计算工作量比隐式的小；
（3）隐式的稳定范围比显式的大。
例3 试分别用阿达姆斯四步显式和三步显式方法求 解下列初值问题，并比较两者结果的精度。
45
解 取h=0.1。两种方法的具体算式如下： 四步显式：
三步隐式：
（4.9）
46
即
（4. 10）
预测点（初值）可用同阶的R-K法计算。本例题
1 2 5 12
3 8
251 720
95 288
10987 60480
38
当k=3时
（4.7）
其中k=3时的（4.7）称为阿达姆斯四步显示方法，因
为它要用到前面四个节点上的 f 值，是一种常用的多
步算法，其精度为四阶。
（4.2）中的 f 也可用牛顿内插多项式来逼近。设
由
39
共k+1个数据来构造一个牛顿内插多项式 上面类似推导可得：
i2
i1
i
41
即为梯形格式，其局部截断误差为
R h 1 3 y( )
i1
12
i1
当k=2时
y y h (5 f 8 f f )
i1 i 12
i1
i
i1
R h y 1 i1 24
4
(4)
() i1
42
当k=3时
y y h (9 f 19 f 5 f f )
i1
i 24
i1
（3.4）
成立，格式（3.3）的局部截断误差就等于O(h3) ，从
而能具有二阶精度。
7
（3.4）中有三个待定系数: λ1, λ2 和 p ，但却只有两
个方程式，因此还有一个自由度。凡满足条件（3.4） 的一族格式（3.3）统称为二阶龙格-库塔格式。
当 p=1, λ1 = λ2 =1/2时，二阶R—K格式（3.3）
作量和获得较高的精度，可采用如下的计算公式：
（4.1）
28
其中 和 为常数，
。
（4.1）是用前面若干节点处的函数值与导数值的
线性组合来计算
的近似值 ，所以通称为
线性多步法。
当 时为单步法； 时为多步法（k步法）。 当 时为显式； 时为隐式。如中点欧拉格式 （2.8）就属于多步法（二步法）。
（4.1）的系数可利用泰勒展开后的待定系数法来 确定。
中的平均斜率K的近似值，故格式（3.5）也称为中
点格式。粗看起来，yi+1=yi+hK2中只含有一个斜率
值K2 ，但实际上K2 是通过K1 才能算出来的，因此，
式
9
中还隐含着K1 。这样，每完成一步仍需计算函数 f 值
两次，其计算工作量仍与改进的欧拉格式一样。 3.2 高阶龙格-库塔公式
要提高精度，可在区间[xi,xi+1]内取二个节点：
逼近（4.2）中的 f ，以导出一种线性多步法——阿达
姆斯（Adams）方法。
设由
共 个数据来构造一个 k 阶牛
顿向后插值多项式
，则由第二章可知，k 阶牛
30
顿向后插值多项式为 其余项为 将其代入（4.2），且利用
，可得：
31
当 上式中
有界，且
（4.3）
时，可以估计出
此即为截断误差项。于是可得
（4.4）
时为止，并以上一次步长的计算结果作
为
；
26
（2） 如果
，则反复减半步长进行计算，直到
时为止，并取其最后一次步长的计算结
果作为 。
这样做时，为了选择步长，每一步都要反复判别
增加了工作量，但在方程的解 变化剧烈的情况下， 总的计算工作量可以得到减少，结果还是合算的。
27
§4 亚当姆斯方法
4.1 线性多步方法
3.3 步长的自动选择
用数值法求接微分方程的过程中，选取适当的 步长至关重要。如步长太大则达不到精度要求；步 长太小则步数太多，不但会增加计算工作量，而且 可能导致舍入误差的严重积累。尤其是当微分方程 的解 变化较剧烈时，步长的合理取法是在变化 剧烈处步长应取得小些；在变化平缓处步长应取的 大些，也就是采用自动选择步长的变步长方法，即 根据精度要求先估计出下一步长的合理大小，然后 ，按此进行计算。这里介绍李查逊（Richardson） 外推法。
32
当k=0时： 即为显式欧拉格式，其局部断截误差为 当k=1时：
33
……..（4.5）
34
式中
将（4.5）代入（4.4）后，可得： 式中
（4.6）
35
代入（4.6）后得 其中
36
局部截断误差为
对应于不同的k值下的 和 值可算出，分 别列于表7-3中。
从表7-3可知： 当k=2时
37
表7-3
很难化成 显式表达式，只能用迭代的办法，这 就会增加计算工作量。
49
因此，在实际计算中，往往仿照改进的欧拉格式 的构造方法，把显式和隐式两种亚当姆斯格式结 合起来，构成预测一校正系统。
以四阶阿达姆斯方法为例，先由显示方法算 出近似值，作为隐式方法的预测值，然后再作校 正
预测
校正
（4.11）
50
用上式计算 时既要用到它前一步的信息 和 ，
12
(3.7)
其中共有七个待定系数：
，
但只有五个方程式，因此还有两个自由度。凡满足
条件（3.7）的一族格式（3.6）统称为三阶龙格-库
塔格式。
13
当待定系数取为
时的三阶龙格-库塔格式称为库塔格式，其具体 形式为：
(3.8)
14
继续推广这种处理过程，可得四阶龙格-库塔格 式。最常用的一种经典龙格-库塔格式的具体形式如 下：
j
k 012
01
1
3 2
1 2
2
23
16
12

12
5 12
3
55
59
24
24
37 24
4
1901 2774
720
720
2616 720
5
4777 7923 9982 1440 1440 1440
3
9 24
1274 720
7298 1440
4
251 720 2877 1440
5
476 1440
前面的几种步进方法，在计算yi+1时大多只用到前 一个节点上的近似值 yi ，而没有用到前几步计算所得出
的信息，故称单步法。实际上，经过多次单步法计算后，
已得出一系列近似值y0 , y1 , … , yi 和 y'0 , y'1 , … , y'i 等。为了充分利用这些信息来计算 yi+1 ，以减少计算工
yi来近似地预测的。
3
由此可以设想，取区间[xi,xi+1] 内的某一个节点
上的斜率K2 与(xi , yi) 点上的斜率K1作线性组合（即
加权平均），作为平均斜率K的近似值，即
为了要得到(xi +p, yi+p) 点上的斜率 K2 ， 需先预测
yi+p = yi +ph K1
4
根据预测值 yi+p 再来算出K2 ：
646 720 1487 1440
1 12
5 24
261 720
798 1440
1 24
106 720 482 1440
19 720
173 1440
27 1440
1 2
1 12
1 24
19 720
3 160
863 60480
44
阿达姆斯方法显式与隐式的比较如下：
（1）同一阶数下，隐式的局部截断误差的系数
4.2 显式和隐式阿达姆斯格式
线性多步法（4.1）中的某些特殊公式，也可由
29
数值积分法来导出。
将微分方程（1.1）的两边从 到 改写成等价的积分方程：
积分，即可
(4.2)
由于上式中的 f 内包含了未知函数 ，因此，
等式右边的积分不能直接算出 ，但可利用数值积分公 式导出一系列计算方法。现用牛顿向后插值多项式来
i
i1
i2
……（4.8）
R h y 19
i1
720
5
(5)
() i1
其中k=3时的（4.8）称为三步四阶阿达姆斯隐式算法， 它只用到前面三个节点上的 f 值，但其精度却为四阶。
43
表7-4
j
k0
01
11 2
25 12
39 24
4 251 720
5 475 1440
12345
1 2
8 12
19 24
由此构成的计算格式为
(3.3)
上式含有三个待定参数: λ1, λ2 和 p。适当选定其值
可使算法的局部截断误差为O(h3) ，即有二阶精度。
5
假定yi = y(xi ) ，分别将K1 和K2 作泰勒展开：
代入（3.3），得
将它与y(x ) 在xi+1 处的二阶泰勒展开式
6
进行比较系数后可知，只要
精确解
19
三种格式的 计算结果分别列于表7—2中。 上例中对三种方法采用了不同的步长h值，是为
了使它们所耗的计算工作量大致相同，以便于比较。 从表7—2可见，经典R—K方法的精度较改进的欧拉 方法又有了很大的提高。关于这一结论也可以从理 论上大致分析出来：
（1）欧拉方法的局部截断误差为
计算四步后截断误差为
即 由上式可见，如取
24
作为近似值，则 可能比 和 的精度要高。
当
时，可取
这种修正的方法与龙贝格数值积分法的思想是一致的。
（3.11）除以（3.10）得
25
由此可得误差事后估计式：
于是，可从步长减半前后的两次计算结果的偏差
来判别所选步长是否适当。当要求的数值精度为 时：
（1） 如果 ，则反复加倍步长进行计算，直到
§3 龙格—库塔方法
要进一步提高求解的精度，可用一种高精度的单 步法—龙格—库塔（Runge—Kutta）方法，简称 R—K法。它采用了间接使用泰勒级数法的技术。
3.1 龙格—库塔公式的导出
对于一阶常微分方程（1.1）的解 y = y (x) ，利
用微分中值定理得
即
1
也即 式中
(3.1) (3.2)