常微分方程的线性多步法

p 1
r
p 1
( p 1) ( k ) p k ] y p 1 ( x n ) O( h p 2 )
k 1
第八章常微分方程数值解法
2 p 使 y xn ，h，h ， ，h 的系数为零，得到关于 k 和 k 的线性方程组

k 0 r k 1
y xn1 y xn
xn1
xn
f x，y x dx。 （8.4.1）
如推导Newton-Cotes求积公式一样，用等距节点的插值多项式来替代被积函 数，再对插值多项式积分，这样就得到一系列求积公式。 例如，用梯形方法计算积分项

x n 1
xn
h f x，y x dx f xn，y xn f xn1，y xn1 ， 2
也可以从得到的显式和隐式 Adama公式，有局部截断误差的定义来求出方
法的局部截断误差。表 8-8中列出了它们的局部截断误差的主项，有表 8-8 可以看出，Adams隐式方法的局部截断误差要小。
表 8-8
r 0 1
5 3 h y xn 12
2
3 4 4 h y xn 8
3
251 5 5 h y xn 720
（8.4.5）
现利用Taylor展开定理，确定线性多步公式（8.4.4）中的待定参数 k， k ， 使她达到
p
阶精度，即 Tn1 O h
p 1
。
对（8.4.5）式的右端各项在
xn 点处作Taylor展开有
第八章常微分方程数值解法
kh j y j x kh p1 y p1 x O h p 2 ， y xn k n n p 1! j! j 0 p kh j y j x kh p1 y p1 x O h p1 。 y x n k n n ! p! j 1 j 1
第八章常微分方程数值解法
yn 1 yn h rj f n j，
j 0
r
（8.4.2）
其中
r 1 1 x n1 tk rj l j x dx dt，j 0， ，r。 1 xn 0 h k 0，k j k j
由此可得（8.4.2）中的系数，其具体数值见表8-6。公式（8.4.2）是一个r+1 步的显式公式，称为Adams显式公式。r=0时，即为Euler公式。
可见，三秒末跳伞员的末速度约有 21 ft / sec 。 若将模型修改为 p=1.1，取 h=0.2，则有计算结果：
v0 0，v1 5.3216，v2 8.8911，v3 11 .2565， v4 12 .8630，v5 13 .9411，v6 146674 ，v7 15 .1552， v8 15 .4830，v9 15 .7030，v10 15 .8508，v11 15 .9500， v12 16 .0165，v13 16 .0612，v14 16 .0912，v15 16 .1113。
公式有
r r
x x 一般地，设由 r 1 个数据点 xn，f n ， n1，f n1 ， ， n r，f n r
Pr x f n j l j x ，l j x
j 0
k 0，k j

x xn k 。 xn j xn k
k 1


第八章常微分方程数值解法 下面我们构造几个著名的四阶线性多步公式，考虑下列形式的公式
yn1 0 yn 1 yn1 2 yn 2 3 yn 3
5 4
h 1 f n1 0 f n 1 f n1 2 f n 2 3 f n 3 ，
p 为 F kv ，其中 1 p 2 ，比例系数 k 依赖于物体的大小、形状，空气
的密度和粘度。跳伞员下落的速度可描述为下列模型：
第八章常微分方程数值解法
dv p m k v mg，v0 0， dt
负号表示下降。显然，当 1< p <2 时，适合于数值方法求解。 设 k / m =1.5，g=32，先用中点法提供开始值，再用下列两步而阶方法
第八章常微分方程数值解法 代入（8.4.1）式有
y xn1 y xn
h f xn，y xn f xn1，y xn1 。 2
据此即可导出公式（8.1.4）。
f ， 构造插值多项式 Pr x ，这里， k f xk，yk xk x0 kh 。 运用插值
第八章常微分方程数值解法 可见三秒末跳伞员的末速度减慢了。计算结果如下图所示
0
-5
-10
-15
-20
-25
0
0.5
1
1.5
2
2.5
3
+ 表示 p=1时的解，* 表示 p=1.1时的解
第八章常微分方程数值解法
节点。这样的插值多项式 Pr x 在求积区间 x n，x n1 上逼近 f x，y x 是一
多项式 Pr x ，则有
Pr x f n j 1l j x ，l j x
j 0
r
k 0，k j

r
x x n k 1 。 xn j 1 x n k 1
于是我们有如下的计算公式
第八章常微分方程数值解法
yn 1 yn h rj f n j 1，
基于数值积分可以构造出一系列求解常微分方程的计算公式，下面介绍
基于 Taylor 展开的待定系数法，它可灵活地构造出线性多步法。对固定的系
数，可以选取待定系数使线性多步法的阶尽可能高。还可以根据需要，确定 显
式还是隐式。 设构造如下具有 p 阶精度的线性多步公式 yn1 0 yn 1 yn1 r yn r h 1 f n1 0 f n r f n r 。 (8.4.4) 当 1 0 时，则（8.4.4）为显式多步式。当 1 0 时，（8.4.4）为隐式 多 步式。它们的局部截断误差为
j 0
r
（8.4.3）
其中
r 0 1 xn1 tk rj l j x dx dt，j 0，，r， 1， xn 1 h k 0，k j k j
其具体数值见表 8-7。公式（8.4.3）是隐式公式，称为 Adams隐式公式。 r=0，1时分别为隐式 Euler公式和梯形公式。
1 2 Adams显式公式 2 h y x n
1 2 1 3 1 h y x n h4 y 4 xn Adams隐式公式 h y xn 2 12 24

19 5 5 h y xn 720
第八章常微分方程数值解法
8.4.2
基于Taylor展开的方法
第八章常微分方程数值解法
8.4
线性多步法
8.4.1 基于数值积分的方法
8.4.2
基于Taylor展开的方法
8.4.2
基于Taylor展开的方法
第八章常微分方程数值解法
8.4
线性多步法
常微分方程初值问题（8.1.1）的数值解法中，除了Runge-
Kutta型公式等单步法之外，还有另一种类型的解法，即某一步的
s 1
r s yn h r 0 f x n1，yn1 rj f n j 1 ，s 0， ， 1， j 1


第八章常微分方程数值解法 迭代收敛条件为 h r 0 L 1，其中 L为f关于y 的Lipschitz常数 利用插值多项式的余项，可以求出 Adams方法的局部截断误差。当然
（8.4.7）
Tn1
3 3 h5 [1 1 k 5 k k ] y 5 x n O h6 。 5! k 1 k 1

（8.4.8）
由于 r=3，p=4 ，由（8.4.6）得到5个方程，而（8.4.7）中有9个为知量， 因此，（8.4.7）中有4个自由度。 若取 1 0，1 2 3 0 ，由（8.4.6）式得到其他5个待定参数的 方程组，解之得
r
k
1，
j
k
k
j k
k 1
r
j 1
（8.4.6）
k 1，j 1， ，p。 2，
而且得到线性多步法的局部截断误差
Tn1
r h p 1 [1 k p 1! k 1 r p
p 1
k
p 1 k k ] y p 1 x n O h p 2 。
第八章常微分方程数值解法
表 8-6
j 0 1 3 23 55 -1 -16 -59 -2774 5 37 2616 -9 -1274 251 1 2 3 4
0 j
21 j 12 2 j 24 3 j
720 4 j 1901
应用实例： 考虑跳伞员的下落速度。 自由落体运动可用牛顿第二定律描述：F=ma。实验表明，空气阻力模型
p
将它们代入（8.4.5）式整理后得
Tn1
r h! (1 k ) y xn [1 ( k ) k k 0 j 1 j! k 1
r
p
j
j ( k ) k ] y j xn
k 1 r
r
j 1
h [1 ( k ) k p 1! k 1
第八章常微分方程数值解法
Tn1
r r y x n1 k y xn k h k f xn k，y xn k ， k 1 k 0
利用原微分方程，有
Tn1
r r y x Fra Baidu bibliotek1 k y xn k h k y x n k 。 k 1 k 0
第八章常微分方程数值解法
表8-7
j 0 1 1 5 9 251 1 8 19 646 -1 -5 -264 1 106 -19 1 2 3 4
0 j
2 1 j 12 2 j 24 3 j 720 4 j
对于隐式公式（8.4.3），需要用迭代求解。确定 yn 1 的迭代公式为
yn1
在上述Adams显式公式的推导中，选用了xn，xn1， ，xn r 作为插值
个外推结果。为了改善逼近效果，我们变外推为内推，即改用xn1，xn， ，xn r 1
x x 为插值节点，用数据点 xn1，f n1 ， n，f n ， ， n r 1，f n r 1 构造插值
公式不仅与前一步解的值有关，而且与前若干步解的值有关，利 用前面多步的信息预测下一步的值，这就是多步法的基本思想，
可以期望获得较高的精度。构造多步法有多种途径，下面先讨论
基于数值积分的方法。
第八章常微分方程数值解法
8.4.1 基于数值积分的方法
将（8.1.1）中的方程在区间 xn，xn 1 上积分，可以得到
vi 1 vi
h fi fi 1 ，i 1，， 2 2
求其他需要计算的值。当 p=1 时，取 h=0.2 有
第八章常微分方程数值解法
v0 0，v1 5.4400，v2 9.3920，v3 12 .3816， v4 14 .6187，v5 16 .2975，v6 17 .5564，v7 18 .5007， v8 19 .2088，v9 19 .7400，v10 20 .1383，v11 20 .4371， v12 20 .6611，v13 20 .8292，v14 20 .9552，v15 21 .0497。