第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）统称为三阶龙格-库
塔格式。