常微分方程
目录
1 常微分方程初值问题 3
2 插值型数值微商近似导数 3
2.1 一阶插商近似一阶导数-Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 高阶格式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3.1 Euler 法的误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3 Runge-Kutta 方法 6
3.1 思想 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 显式 Runge-Kutta 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2.1 一级 Runge-Kutta 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2.2 二级 Runge-Kutta 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3 隐式 Runge-Kutta 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4 单步法的性质 11
4.1 截断误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2 相容性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.3 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.4 稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
5 线性多步法 13
1
5.1 一阶常微分方程多步法的一般形式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5.2 基于数值积分的线性多步法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5.3 预估-校正法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
6 多步法的性质 14
6.1 误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.2 差分方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.2.1 差分方程与特征方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.2.2 特征方程的根 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.3 差分方法的相容性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
6.4 差分方法的收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.5 差分方法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
7 常微分方程组的数值解法 17
1 常微分方程初值问题
常微分方程初值问题为
𝑢
′
= 𝑓 (𝑥, 𝑢), 𝑥 ∈ [𝑎, 𝑏]
𝑢(𝑎) =
˜
𝑢
期望解是唯一的, 稳定的, 初值的微小改变不会导致解的巨大改变
若 𝑓 (𝑡, 𝑢) 是实函数, 并且在矩形区域 (𝑥, 𝑢) ∈ [𝑎, 𝑏] × (−∞, +∞) 上连续, 并且 𝑓 (𝑥, 𝑢) 对 𝑢 满足 Lipschitz
条件:
∃𝐿 > 0, ∀𝑢
1
, 𝑢
2
∈ ℝ, | 𝑓 (𝑥, 𝑢
1
) − 𝑓 (𝑥, 𝑢
2
)| ≤ 𝐿|𝑢
1
− 𝑢
2
|
2 插值型数值微商近似导数
2.1 一阶插商近似一阶导数-Euler 法
前差
𝑢
′
(𝑥)
𝑥=𝑥
0
=
𝑢(𝑥
0
+ ℎ) − 𝑢(𝑥
0
)
ℎ
+ 𝑂(ℎ) ≈
𝑣
𝑛+1
− 𝑣
𝑛
ℎ
得到向前 Euler 法
𝑣
𝑗+1
= 𝑣
𝑗
+ ℎ 𝑓 (𝑥
𝑗
, 𝑣
𝑗
)
𝑣
0
=
˜
𝑢
这种可以直接代入求值的方法称为显式方法
后差
𝑢
′
(𝑥)
𝑥=𝑥
0
=
𝑢(𝑥
0
) − 𝑢(𝑥
0
− ℎ)
ℎ
+ 𝑂(ℎ) ≈
𝑣
𝑛
− 𝑣
𝑛−1
ℎ
得到向后 Euler 法
𝑣
𝑗+1
= 𝑣
𝑗
+ ℎ 𝑓 (𝑥
𝑗+1
, 𝑣
𝑗+1
)
𝑣
0
=
˜
𝑢
这并不能直接代入求值, 需要求解方程, 称为隐式方法. 隐式方法的优点是稳定性好
中心差
𝑢
′
(𝑥)
𝑥=𝑥
0
=
𝑢(𝑥
0
+ ℎ) − 𝑢(𝑥
0
− ℎ)
2ℎ
+ 𝑂(ℎ
2
) ≈
𝑣
𝑛+1
− 𝑣
𝑛−1
2ℎ
得到中心差法
𝑣
𝑗+1
= 𝑣
𝑗−1
+ 2ℎ 𝑓 (𝑥
𝑗
, 𝑣
𝑗
)
𝑣
0
=
˜
𝑢, 𝑣
1
=
˜
𝑢 + ℎ 𝑓 (𝑥
0
,
˜
𝑢)
它是显示多步法
2.2 高阶格式
可以通过插值型数值微商近似导数得到高阶格式, 如
𝑣
𝑛+1
= 𝑎
0
𝑣
𝑛−2
+ 𝑎
1
𝑣
𝑛−1
+ 𝑎
2
𝑣
𝑛
+ 𝑏 𝑓 (𝑥
𝑛
, 𝑣
𝑛
)
其中
𝑣
𝑛
≈
𝑢
(
𝑥
𝑛
)
取
𝑥
0
≡ 𝑥
𝑛−2
, 𝑥
1
≡ 𝑥
𝑛−1
, 𝑥
2
≡ 𝑥
𝑛
, 𝑥
3
≡ 𝑥
𝑛+1
利用拉格朗日插值多项式
𝑢(𝑥) = 𝑢
𝑛−2
𝐿
0
(𝑥) + 𝑢
𝑛−1
𝐿
1
(𝑥) + 𝑢
𝑛
𝐿
2
(𝑥) + 𝑢
𝑛+1
𝐿
3
(𝑥)
对其求导得
𝑢
′
(𝑥) = 𝑢
𝑛−2
𝐿
′
0
(𝑥) + 𝑢
𝑛−1
𝐿
′
1
(𝑥) + 𝑢
𝑛
𝐿
′
2
(𝑥) + 𝑢
𝑛+1
𝐿
′
3
(𝑥)
即得导数近似
2.3 误差
衡量误差需要考察三个方面:
1. 局部截断误差 (近似程度)
2. 数值解是否收敛 (步长趋于 0 时, 数值解是否趋于解析解)
3. 数值解是否稳定 (计算机舍入误差是否会随计算步骤增加而无限放大)
2.3.1 Euler 法的误差
考察局部截断误差, 即假定前面的计算是准确的, 局部截断误差为用一阶差商近似导数的误差
𝑇
𝑗
≡
𝑢(𝑥 + ℎ) − 𝑢(𝑥)
ℎ
− 𝑓 (𝑥, 𝑢)
𝑥=𝑥
𝑗
−
d𝑢(𝑥)
d𝑥
− 𝑓 (𝑥, 𝑢)
𝑥=𝑥
𝑗
由于一般取
d𝑢
d𝑥
− 𝑓 (𝑥, 𝑢) = 0, 所以
𝑇
𝑗
≡
𝑢(𝑥 + ℎ) − 𝑢(𝑥)
ℎ
− 𝑓 (𝑥, 𝑢)
𝑥=𝑥
𝑗
前差
𝑇
𝑛
≡
𝑢(𝑥
𝑛+1
) − 𝑢(𝑥
𝑛
)
ℎ
𝑥=𝑥
𝑛
− 𝑢
′
(𝑥
𝑛
) =
ℎ
2
𝑢
′′
(𝜉
𝑛
), 𝜉
𝑛
∈ (𝑥
𝑛
, 𝑥
𝑛+1
)
后差
𝑇
𝑛
≡
𝑢(𝑥
𝑛
) − 𝑢(𝑥
𝑛−1
)
ℎ
𝑥=𝑥
𝑛
− 𝑢
′
(𝑥
𝑛
) = −
ℎ
2
𝑢
′′
(𝜉
𝑛
), 𝜉
𝑛
∈ (𝑥
𝑛−1
, 𝑥
𝑛
)
中心差
𝑇
𝑛
≡
𝑢(𝑥
𝑛+1
) − 𝑢(𝑥
𝑛−1
)
2ℎ
𝑥=𝑥
𝑛
− 𝑢
′
(𝑥
𝑛
) =
ℎ
2
6
𝑢
′′′
(𝜉
𝑛
), 𝜉
𝑛
∈ (𝑥
𝑛−1
, 𝑥
𝑛+1
)
若截断误差满足
𝑇 = 𝑂 (ℎ
𝑝
)
则称该方法具有 𝑝 阶精度
考虑整体误差与收敛性
𝑒
𝑛
≡ 𝑢(𝑥
𝑛
) − 𝑣
𝑛
对于前差 Euler 法
𝑇
𝑗
=
𝑢
𝑗+1
− 𝑢
𝑗
ℎ
− 𝑓 (𝑥
𝑗
, 𝑢
𝑗
) =
1
2
ℎ𝑢
′′
(𝜉
𝑗
)
那么
𝑢
𝑗+1
= 𝑢
𝑗
+ ℎ 𝑓 (𝑥
𝑗
, 𝑢
𝑗
) + ℎ𝑇
𝑗
因而
𝑒
𝑗+1
= 𝑒
𝑗
+ ℎ[ 𝑓 (𝑥
𝑗
, 𝑢
𝑗
) − 𝑓 (𝑥
𝑗
, 𝑣
𝑗
)] + ℎ𝑇
𝑗
取绝对值得到不等式
𝑒
𝑗+1
≤ ℎ
𝑓 (𝑥
𝑗
, 𝑢
𝑗
) − 𝑓 (𝑥
𝑗
, 𝑣
𝑗
)
+ ℎ
𝑇
𝑗
有 Lipschitz 条件, 所以
𝑓 (𝑥
𝑗
, 𝑢
𝑗
) − 𝑓 (𝑥
𝑗
, 𝑣
𝑗
)
≤ 𝐿
𝑢
𝑗
− 𝑣
𝑗
= 𝐿
𝑒
𝑗
所以
𝑒
𝑗+1
≤ (1 + 𝐿ℎ)
𝑒
𝑗
+ ℎ
𝑇
𝑗
认为 𝑢 的性质足够好, 那么取最大值
𝑇 = max
𝑗
𝑇
𝑗
得到
𝑒
𝑗+1
≤ (1 + 𝐿ℎ)
𝑒
𝑗
+
𝑇
𝐿
递推得到
|
𝑒
𝑛
|
≤ (1 + 𝐿ℎ)
𝑛+1
|
𝑒
0
|
+
𝑇
𝐿
若初值准确, 则 𝑒
0
= 0. 若 ℎ → 0, 则得到 𝑇 → 0, 所以 𝑒
𝑛
→ 0, 即数值解收敛于解析解
考虑稳定性, 对于两个初值
˜
𝑢, 𝑧
0
, 由数值方法得到的数值解分别为 𝑣
𝑛
, 𝑧
𝑛
, 则希望
|
𝑣
𝑛
− 𝑧
𝑛
|
≤ 𝐶
|
˜
𝑢 − 𝑧
0
|
设
˜
𝑢 初值的计算结果为 𝑣
𝑗+1
,𝑧
0
初值的计算结果为 𝑧
𝑗+1
, 希望估计其差值
𝑣
𝑗+1
− 𝑧
𝑗+1
𝑣
𝑗+1
− 𝑧
𝑗+1
=
𝑣
𝑗
+ ℎ 𝑓 (𝑥
𝑗
, 𝑣
𝑗
) − 𝑧
𝑗
− ℎ 𝑓 (𝑥
𝑗
, 𝑧
𝑗
)
≤
𝑣
𝑗
− 𝑧
𝑗
+ ℎ
𝑓 (𝑥
𝑗
, 𝑣
𝑗
) − 𝑓 (𝑥
𝑗
, 𝑧
𝑗
)
由于满足 Lipschitz 条件, 所以
𝑓 (𝑥
𝑗
, 𝑣
𝑗
) − 𝑓 (𝑥
𝑗
, 𝑧
𝑗
)
≤ 𝐿
𝑣
𝑗
− 𝑧
𝑗
所以
|
𝑣
𝑛+1
− 𝑧
𝑛+1
|
≤ (1 + 𝐿ℎ)
|
𝑣
𝑛
− 𝑧
𝑛
|
≤ · · · ≤ (1 + ℎ𝐿)
𝑛+1
|
˜
𝑢 − 𝑧
0
|
而
(1 + 𝐿ℎ)
𝑛+1
≤ 𝑒
(𝑛+1) 𝐿ℎ
≤ 𝑒
(𝑏−𝑎)𝐿
这是一个与步长无关的常数, 所以数值解是稳定的
3 Runge-Kutta 方法
3.1 思想
对于常微分方程初值问题
𝑢
′
= 𝑓 (𝑥, 𝑢), 𝑥 ∈ [𝑎, 𝑏]
𝑢(𝑎) =
˜
𝑢
由中值定理, 对于函数 𝑢(𝑥) 应有
𝑢(𝑥 + ℎ) = 𝑢(𝑥) + ℎ𝑢
′
(𝜉), 𝜉 ∈ (𝑥, 𝑥 + ℎ)
希望找到一个方法能够准确地估计 𝑢
′
(𝜉), 即 (𝑥, 𝑥 + ℎ) 的平均斜率. 由于 𝑢
′
(𝑥) = 𝑓 (𝑥, 𝑢), 所以
𝑢(𝑥 + ℎ) = 𝑢(𝑥) + ℎ 𝑓 (𝜉, 𝑢(𝜉)), 𝜉 ∈ (𝑥, 𝑥 + ℎ)
可以将 (𝑥, 𝑥 + ℎ) 划分为 𝑅 个小区间, 即取
𝑥 + 𝑎
1
ℎ, 𝑥 + 𝑎
2
ℎ, · · · , 𝑥 + 𝑎
𝑅
ℎ
设可以求得这点的斜率 𝐾, 即
𝑥 + 𝑎
1
ℎ 𝑥 + 𝑎
2
ℎ · · · 𝑥 + 𝑎
𝑅
ℎ
𝐾
1
𝐾
2
· · · 𝐾
𝑅
则认为平均斜率为它们的加权平均
𝐾
∗
=
𝑅
Õ
𝑟=1
𝑐
𝑟
𝐾
𝑟
也就是说
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ
𝑅
Õ
𝑟=1
𝑐
𝑟
𝐾
𝑟
理论上 𝐾
𝑟
应该为
𝐾
∗
𝑟
= 𝑓 (𝑥
𝑛
+ 𝑎
𝑟
ℎ, 𝑢(𝑥
𝑛
+ 𝑎
𝑟
ℎ))
但是 𝑢(𝑥
𝑛
+ 𝑎
𝑟
ℎ) 是未知的, 考虑简单地估计
𝐾
1
= 𝑓 (𝑥
𝑛
, 𝑢
𝑛
)
𝐾
2
= 𝑓 (𝑥
𝑛
+ 𝑎
2
ℎ, 𝑢
𝑛
+ 𝑏
21
ℎ𝐾
1
)
𝐾
3
= 𝑓 (𝑥
𝑛
+ 𝑎
3
ℎ, 𝑢
𝑛
+ 𝑏
31
ℎ𝐾
1
+ 𝑏
32
ℎ𝐾
2
)
· · ·
𝐾
𝑅
= 𝑓 (𝑥
𝑛
+ 𝑎
𝑅
ℎ, 𝑢
𝑛
+ 𝑏
𝑅1
ℎ𝐾
1
+ 𝑏
𝑅2
ℎ𝐾
2
+ · · · + 𝑏
𝑅,𝑅−1
ℎ𝐾
𝑅−1
)
由于可以显式地求解每一个 𝐾
𝑟
, 所以称为显式 Runge-Kutta 方法. 与之对应的是隐式 Runge-Kutta
方法, 它是显示方法的推广,𝐾
𝑟
的计算需要求解方程
𝐾
1
= 𝑓 (𝑥
𝑛
, 𝑢
𝑛
+ 𝑏
11
ℎ𝐾
1
+ 𝑏
12
ℎ𝐾
2
+ · · · + 𝑏
1,𝑅−1
ℎ𝐾
𝑅−1
)
𝐾
2
= 𝑓 (𝑥
𝑛
+ 𝑎
2
ℎ, 𝑢
𝑛
+ 𝑏
21
ℎ𝐾
1
+ 𝑏
22
ℎ𝐾
2
+ · · · + 𝑏
2,𝑅−1
ℎ𝐾
𝑅−1
)
· · ·
𝐾
𝑅
= 𝑓 (𝑥
𝑛
+ 𝑎
𝑅
ℎ, 𝑢
𝑛
+ 𝑏
𝑅1
ℎ𝐾
1
+ 𝑏
𝑅2
ℎ𝐾
2
+ · · · + 𝑏
𝑅,𝑅−1
ℎ𝐾
𝑅−1
)
即
𝑢(𝑥
𝑛
+ 𝑎
𝑟
ℎ) ≈ 𝑢
𝑛
+ ℎ
𝑅
Õ
𝑠=1
𝑏
𝑟 𝑠
𝐾
𝑠
那么
𝐾
𝑟
= 𝑓 (𝑥
𝑛
+ 𝑎
𝑟
ℎ, 𝑢
𝑛
+ ℎ
𝑅
Õ
𝑠=1
𝑏
𝑟 𝑠
𝐾
𝑠
), 𝑟 = 1, 2, · · · , 𝑅
3.2 显式 Runge-Kutta 方法
若有 𝑎
1
= 0, 并且 𝑟 < 𝑠 时有 𝑏
𝑟 𝑠
= 0, 则称为 𝑅 级显式 Runge-Kutta 方法. 一般的 𝑅 级显式 Runge-Kutta
法可以写为” 列阵形式”
0
𝑎
2
𝑏
21
𝑎
3
𝑏
31
𝑏
32
.
.
.
.
.
.
.
.
.
.
.
.
𝑎
𝑅
𝑏
𝑅1
𝑏
𝑅2
· · · 𝑏
𝑅,𝑅−1
𝑐
1
𝑐
2
· · · 𝑐
𝑅−1
𝑐
𝑅
其中的系数 𝑐
𝑟
, 𝑏
𝑟 𝑠
, 𝑎
𝑟
需要通过截断误差的要求给出, 系数不一定唯一. 对 𝑢 做展开即可得到误差
𝑇
𝑛
=
𝑢
𝑛+1
− 𝑢
𝑛
ℎ
−
𝑅
Õ
𝑟=1
𝑐
𝑟
𝐾
𝑟
!
𝑥=𝑥
𝑛
其中
𝑢(𝑥
𝑛+1
) − 𝑢(𝑥
𝑛
)
ℎ
=
𝑢
′
+
ℎ
2
𝑢
′′
+ · · · +
ℎ
𝑝−1
𝑝!
𝑢
( 𝑝 )
𝑥=𝑥
𝑛
,𝑢=𝑢
𝑛
+ 𝑂(ℎ
𝑝
)
利用微分方程可以求得导数
𝑢
′
= 𝑓 (𝑥, 𝑢)
𝑢
′′
=
d
d𝑥
𝑓 (𝑥, 𝑢) = 𝑓
𝑥
+ 𝑓 𝑓
𝑢
𝑢
′′′
=
d
d𝑥
( 𝑓
𝑥
+ 𝑓 𝑓
𝑢
) = 𝑓
𝑥𝑥
+ 2 𝑓 𝑓
𝑥
+ 𝑓
2
𝑓
𝑥𝑢
+ 𝑓 𝑓
𝑢𝑢
· · ·
后一项 𝐾
𝑟
的求和为
𝑅
Õ
𝑟=1
𝑐
𝑟
𝐾
𝑟
=
𝑅
Õ
𝑟=1
𝑐
𝑟
𝑓 (𝑥, 𝑢(𝑥))
𝑥=𝑥
𝑛
可以由此给出截断误差
3.2.1
一级
Runge-Kutta
方法
𝑅 = 1, 要求 𝑇 = 𝑂 (ℎ)
3.2.2 二级 Runge-Kutta 方法
𝑅 = 2, 要求 𝑇 = 𝑂 (ℎ
2
), 有
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ(𝑐
1
𝐾
1
+ 𝑐
2
𝐾
2
)
其中
𝐾
1
= 𝑓 (𝑥
𝑛
, 𝑢
𝑛
), 𝐾
2
= 𝑓 (𝑥
𝑛
+ 𝑎
2
ℎ, 𝑢
𝑛
+ 𝑏
21
ℎ𝐾
1
)
展开即可得方程
𝑐
1
+ 𝑐
2
= 1
𝑐
2
𝑎
2
=
1
2
𝑐
2
𝑏
21
=
1
2
3.3 隐式 Runge-Kutta 方法
考虑
𝑢(𝑥
𝑛
+ ℎ) = 𝑢(𝑥
𝑛
) +
∫
𝑥
𝑛
+ℎ
𝑥
𝑛
𝑢
′
(𝑥)d𝑥
由于 𝑢
′
(𝑥) = 𝑓 (𝑥, 𝑢), 所以
𝑢(𝑥
𝑛
+ ℎ) = 𝑢(𝑥
𝑛
) +
∫
𝑥
𝑛
+ℎ
𝑥
𝑛
𝑓 (𝑥, 𝑢(𝑥))d𝑥
考虑估计右侧积分, 最准确的方式是使用 Gauss-Legendre 求积公式, 将其归一化到 [−1, 1] 上
∫
𝑥
𝑛
+ℎ
𝑥
𝑛
𝑓 (𝑥, 𝑢(𝑥))d𝑥 =
ℎ
2
∫
1
−1
𝑓
𝑥
𝑛
+
ℎ
2
(1 + 𝜉), 𝑢
𝑥
𝑛
+
ℎ
2
(1 + 𝜉)
d𝜉
就可以对其使用 Gauss-Legendre 求积公式
∫
1
−1
𝑓
𝑥
𝑛
+
ℎ
2
(1 + 𝜉), 𝑢
𝑥
𝑛
+
ℎ
2
(1 + 𝜉)
d𝜉 ≈
𝑅
Õ
𝑖=1
𝑐
𝑖
𝑓
𝑥
𝑛
+
ℎ
2
(1 + 𝜉
𝑖
), 𝑢
𝑥
𝑛
+
ℎ
2
(1 + 𝜉
𝑖
)
其中 𝜉
𝑖
为 Gauss-Legendre 求积公式的节点, 为了简单起见, 不妨记
𝑎
𝑖
=
1
2
(1 + 𝜉
𝑖
)
那么公式化为了
∫
𝑥
𝑛
+ℎ
𝑥
𝑛
𝑓 (𝑥, 𝑢(𝑥))d𝑥 ≈
𝑅
Õ
𝑖=1
𝑐
𝑖
𝑓
[
𝑥
𝑛
+ 𝑎
𝑖
ℎ, 𝑢
(
𝑥
𝑛
+ 𝑎
𝑖
ℎ
)]
𝑐
𝑖
为 Gauss-Legendre 求积公式的系数, 它可以如下求出
𝑐
𝑖
=
∫
1
−1
𝑅
Ö
𝑗=1, 𝑗≠𝑖
𝜉 − 𝜉
𝑗
𝜉
𝑖
− 𝜉
𝑗
d𝜉
不妨再记
𝐾
𝑖
= 𝑓
[
𝑥
𝑛
+ 𝑎
𝑖
ℎ, 𝑢
(
𝑥
𝑛
+ 𝑎
𝑖
ℎ
)]
那么最初的递推式就变为了
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ
𝑅
Õ
𝑖=1
𝑐
𝑖
𝐾
𝑖
为了得到 𝐾
𝑖
, 注意到 𝐾
𝑖
的定义
𝐾
𝑖
= 𝑓
[
𝑥
𝑛
+ 𝑎
𝑖
ℎ, 𝑢
(
𝑥
𝑛
+ 𝑎
𝑖
ℎ
)]
需要求出 𝑢(𝑥
𝑛
+ 𝑎
𝑖
ℎ), 但是 𝑢(𝑥) 是未知的, 但是有积分式
𝑢(𝑥
𝑛
+ 𝑎
𝑖
ℎ) = 𝑢(𝑥
𝑛
) +
∫
𝑥
𝑛
+𝑎
𝑖
ℎ
𝑥
𝑛
𝑓 (𝑥, 𝑢(𝑥))d𝑥
只需要求解后面的积分即可, 注意到 𝐾
𝑖
是 𝑓 (𝑥, 𝑢(𝑥)) 在 𝑥 = 𝑥
𝑛
+ 𝑎
𝑖
ℎ 处的值, 因而可以构造 𝑓 (𝑥, 𝑢(𝑥)) 的
插值多项式
𝑓 (𝑥, 𝑢(𝑥)) ≈
𝑅
Õ
𝑖=1
𝐾
𝑖
𝐿
𝑖
(𝑥)
用它代替原积分式中的 𝑢
′
(𝑥), 即
𝑢(𝑥
𝑛
+ 𝑎
𝑖
ℎ) ≈ 𝑢(𝑥
𝑛
) +
∫
𝑥
𝑛
+𝑎
𝑖
ℎ
𝑥
𝑛
𝑅
Õ
𝑗=1
𝐾
𝑗
𝐿
𝑗
(𝑥)d𝑥 = 𝑢(𝑥
𝑛
) +
𝑅
Õ
𝑗=1
𝐾
𝑗
∫
𝑥
𝑛
+𝑎
𝑖
ℎ
𝑥
𝑛
𝐿
𝑗
(𝑥)d𝑥
为了将其化为更一般的形式, 考虑令
𝑥 = 𝑥
𝑛
+ 𝜉ℎ
那么
∫
𝑥
𝑛
+𝑎
𝑖
ℎ
𝑥
𝑛
𝐿
𝑗
(𝑥)d𝑥 =
∫
𝑎
𝑖
0
𝐿
𝑗
(𝑥
𝑛
+ 𝜉ℎ)ℎd𝜉
而
𝐿
𝑗
(𝑥
𝑛
+ 𝜉ℎ) =
𝑅
Ö
𝑘=1,𝑘≠ 𝑗
(𝑥
𝑛
+ 𝜉ℎ) − (𝑥
𝑛
+ 𝑎
𝑘
ℎ)
(𝑥
𝑛
+ 𝑎
𝑗
ℎ) − (𝑥
𝑛
+ 𝑎
𝑘
ℎ)
=
𝑅
Ö
𝑘=1,𝑘≠ 𝑗
𝜉 − 𝑎
𝑘
𝑎
𝑗
− 𝑎
𝑘
因而
∫
𝑥
𝑛
+𝑎
𝑖
ℎ
𝑥
𝑛
𝐿
𝑗
(𝑥)d𝑥 = ℎ
∫
𝑎
𝑖
0
𝑙
𝑗
(𝜉)d𝜉
其中
𝑙
𝑗
(𝜉) =
𝑅
Ö
𝑘=1,𝑘≠ 𝑗
𝜉 − 𝑎
𝑘
𝑎
𝑗
− 𝑎
𝑘
记
𝑏
𝑖 𝑗
=
∫
𝑎
𝑖
0
𝑙
𝑗
(𝜉)d𝜉
那么原积分式就可以化为
𝑢(𝑥
𝑛
+ 𝑎
𝑖
ℎ) ≈ 𝑢(𝑥
𝑛
) + ℎ
𝑅
Õ
𝑗=1
𝑏
𝑖 𝑗
𝐾
𝑗
代回 𝐾
𝑖
的定义, 就可以得到
𝐾
𝑖
≈ 𝑓
"
𝑥
𝑛
+ 𝑎
𝑖
ℎ, 𝑢
(
𝑥
𝑛
)
+ ℎ
𝑅
Õ
𝑗=1
𝑏
𝑖 𝑗
𝐾
𝑗
#
𝑎
𝑖
=
1
2
(1 + 𝜉
𝑖
), 𝑏
𝑖 𝑗
=
∫
𝑎
𝑖
0
𝑙
𝑗
(𝜉)d𝜉, 𝑐
𝑖
=
∫
1
−1
𝑙
𝑖
(𝜉)d𝜉
其中 𝜉
𝑖
为 Gauss-Legendre 求积公式的节点,𝑙
𝑗
(𝑥) 为 Lagrange 插值多项式的基函数
4 单步法的性质
一般一阶常微分方程
𝑢
′
= 𝑓 (𝑥, 𝑢)
的显式单步格式为
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ𝜑(𝑥
𝑛
, 𝑣
𝑛
, ℎ)
4.1 截断误差
截断误差是假定前面的计算都精确的前提下, 单纯由数值近似导致的误差. 一般显示单步格式在 𝑥
𝑛
处的
截断误差定义为
𝑇
𝑛
≡
𝑢(𝑥
𝑛+1
) − 𝑢(𝑥
𝑛
)
ℎ
− 𝜑(𝑥
𝑛
, 𝑢(𝑥
𝑛
), ℎ)
−
[
𝑢
′
− 𝑓 (𝑥, 𝑢)
]
𝑥=𝑥
𝑛
,𝑢=𝑢( 𝑥
𝑛
)
对于给定的数值方法, 若 𝑝 是使得其截断误差为
𝑇
𝑛
= 𝑂 (ℎ
𝑝
)
的最大整数, 则称该方法具有 𝑝 阶精度
4.2 相容性
若局部截断误差 𝑇
𝑛
满足
𝑇
𝑛
= 𝑂 (ℎ
𝑝
), 𝑝 > 0
则称该单步法与
𝑢
′
=
𝑓
(
𝑥, 𝑢
)
是相容的
.
也可以用极限的方式定义
,
若有
lim
ℎ→0
𝜑(𝑥, 𝑢, ℎ) = 𝑓 (𝑥, 𝑢)
则称该单步法与 𝑢
′
= 𝑓 (𝑥, 𝑢) 是相容的
4.3 收敛性
若对于 𝑥 = 𝑥
𝑛
时有
lim
ℎ→0
𝑣
𝑛
= 𝑢(𝑥
𝑛
)
则称该单步法
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ𝜑(𝑥
𝑛
, 𝑣
𝑛
, ℎ)
是收敛的, 整体误差为
𝑒
𝑛
≡ 𝑢(𝑥
𝑛
) − 𝑣
𝑛
若 𝜑(𝑥, 𝑢, ℎ) 关于 𝑢 满足 Lipschitz 条件
∃𝐿 > 0, ∀𝑢
1
, 𝑢
2
∈ ℝ, |𝜑(𝑥, 𝑢
1
, ℎ) − 𝜑(𝑥, 𝑢
2
, ℎ)| ≤ 𝐿|𝑢
1
− 𝑢
2
|
并且局部截断误差有最大值 𝑇 , 则有整体误差
|
𝑒
𝑛
|
≤ (1 + 𝐿ℎ)
𝑛+1
|
𝑒
0
|
+
𝑇
𝐿
若有 lim
ℎ→0
𝑇 = 0, 则单步法是收敛的, 并且是相容的
对于 𝑚 级 Runge-Kutta 方法
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ
𝑚
Õ
𝑟=1
𝑐
𝑟
𝐾
𝑟
若 𝑓 (𝑥, 𝑢) 关于 𝑢 满足 Lipschitz 条件
∃𝐿 > 0, ∀𝑢
1
, 𝑢
2
∈ ℝ, | 𝑓 (𝑥, 𝑢
1
) − 𝑓 (𝑥, 𝑢
2
)| ≤ 𝐿|𝑢
1
− 𝑢
2
|
则有整体误差
|
𝑒
𝑛
|
≤ (1 + 𝐿ℎ)
𝑛+1
|
𝑒
0
|
+
𝑇
𝐿
1
其中 𝑇 为最大局部截断误差,𝐿
1
定义为
𝐿
1
≡ 2𝐿
𝑚
Õ
𝑟=1
|
𝑐
𝑟
|
4.4 稳定性
若 𝜑(𝑥, 𝑢, ℎ) 关于 𝑢 满足 Lipschitz 条件
∃𝐿 > 0, ∀𝑢
1
, 𝑢
2
∈ ℝ, |𝜑(𝑥, 𝑢
1
, ℎ) − 𝜑(𝑥, 𝑢
2
, ℎ)| ≤ 𝐿|𝑢
1
− 𝑢
2
|
则相容的单步法 𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ𝜑(𝑥
𝑛
, 𝑣
𝑛
, ℎ) 是稳定的
即若存在常数 𝐶 > 0, ℎ
0
> 0 使得对于任意的初值
˜
𝑢, 𝑧
0
, 由数值方法得到的数值解为 𝑣
𝑛
, 𝑧
𝑛
满足
|
𝑣
𝑛
− 𝑧
𝑛
|
≤ 𝐶
|
˜
𝑢 − 𝑧
0
|
考虑以下的方程
𝑢
′
= 𝜆𝑢
𝑢(𝑎) =
˜
𝑢
其中 𝜆 是复数, 并且实部小于零 ℜ(𝜆) < 0, 将数值方法应用于该方程, 若初值发生微小变化, 数值解的变
化也是微小的, 则称该数值方法是绝对稳定的, 也称为 A 稳定的
5 线性多步法
5.1 一阶常微分方程多步法的一般形式
对于一阶常微分方程的初值问题
𝑢
′
= 𝑓 (𝑥, 𝑢), 𝑥 ∈ [𝑎, 𝑏]
𝑢(𝑎) =
˜
𝑢
多步法有一般形式
𝑣
𝑛+1
=
𝑘
Õ
𝑗=1
𝛼
𝑗
𝑣
𝑛+1− 𝑗
+ ℎ𝜙(𝑥
𝑛+1
, 𝑥
𝑛
, · · · , 𝑥
𝑛+1−𝑘
, 𝑣
′
𝑛+1
, 𝑣
′
𝑛
, · · · , 𝑣
′
𝑛+1−𝑘
, ℎ)
线性多步法的 𝜙 是线性函数, 那么
𝑣
𝑛+1
=
𝑘
Õ
𝑗=1
𝛼
𝑗
𝑣
𝑛+1− 𝑗
+ ℎ
𝑘
Õ
𝑗=0
𝛽
𝑗
𝑣
′
𝑛+1− 𝑗
若 𝛽
0
≠ 0, 则为隐式格式
5.2 基于数值积分的线性多步法
将微分方程在 [𝑥
𝑛− 𝑝
, 𝑥
𝑛+1
] 上积分, 得到
𝑢(𝑥
𝑛+1
) = 𝑢(𝑥
𝑛− 𝑝
) +
∫
𝑥
𝑛+1
𝑥
𝑛− 𝑝
𝑓 (𝑥, 𝑢(𝑥))d𝑥
用数值积分公式近似积分, 或是用插值函数近似被积函数, 可以得到线性多步法, 称为差分方法
以 𝑥
𝑛
, 𝑥
𝑛−1
, · · · , 𝑥
𝑛−𝑞
为节点, 构造数值积分公式, 为显示 Adams 方法
以 𝑥
𝑛+1
, 𝑥
𝑛
, · · · , 𝑥
𝑛−𝑞+1
为节点, 构造 𝑓 (𝑥, 𝑢(𝑥)) 的插值多项式, 为隐式 Adams 方法
5.3 预估-校正法
先使用显示方法得到一个预估值
𝑣
𝑛+1
= 𝑣
𝑛
+ ℎ 𝑓 (𝑥
𝑛
, 𝑣
𝑛
)
再使用隐式方法得到一个校正值
𝑣
𝑛+1
= 𝑣
𝑛
+
ℎ
2
[
𝑓 (𝑥
𝑛
, 𝑣
𝑛
) + 𝑓 (𝑥
𝑛+1
, 𝑣
𝑛+1
)
]
6 多步法的性质
6.1 误差
多步法的局部截断误差与单步法相同, 假定前面的计算是准确的, 线性多步法的局部截断误差为
𝑇
𝑛
=
1
ℎ
𝑢
𝑛+1
−
𝐽
Õ
𝑗=1
𝛼
𝑗
𝑢
𝑛+1− 𝑗
− ℎ
𝐽
Õ
𝑗=0
𝛽
𝑗
𝑓
𝑛+1− 𝑗
!
− (𝑢
′
− 𝑓 (𝑥, 𝑢))
𝑥=𝑥
𝑛
=
1
ℎ
𝑢
𝑛+1
−
𝐽
Õ
𝑗=1
𝛼
𝑗
𝑢
𝑛+1− 𝑗
!
−
𝐽
Õ
𝑗=0
𝛽
𝑗
𝑓
𝑛+1− 𝑗
= 𝑂 (ℎ
𝑝
)
相容性的定义也同单步法, 同样有推论
lim
ℎ→0
𝑇
𝑛
= 0 ⇒ 相容
6.2 差分方程
6.2.1 差分方程与特征方程
常系数差分方程的一般形式
𝑎
0
𝑦
𝑗
+ 𝑎
1
𝑦
𝑗+1
+ · · · + 𝑎
𝑛
𝑦
𝑗+𝑛
= 𝑐
𝑗+𝑛
当 𝑐
𝑗+𝑛
= 0 时, 称为齐次差分方程; 当 𝑎
0
与 𝑎
𝑛
都不为零时称为 𝑛 阶差分方程
6.2.2 特征方程的根
对于齐次差分方程, 若给定了初始条件
𝑦
0
, 𝑦
1
, · · · , 𝑦
𝑛−1
则其解唯一确定. 可以有特征方程
𝑎
0
+ 𝑎
1
𝑥 + 𝑎
2
𝑥
2
+ · · · + 𝑎
𝑛
𝑥
𝑛
= 0
若 𝑥 是特征方程的根, 则
𝑦
𝑗
= 𝛼𝑥
𝑗
为齐次差分方程的解
若特征方程有 𝑛 个互异的根 𝑥
1
, 𝑥
2
, · · · , 𝑥
𝑛
, 则则差分方程的通解为
𝑦
𝑗
=
𝑛
Õ
𝑘=1
𝛼
𝑘
𝑥
𝑗
𝑘
其中 𝛼
𝑘
由初始条件 𝑦
0
, 𝑦
1
, · · · , 𝑦
𝑛−1
唯一确定
若特征方程有重根及对应的重数为
𝑥
1
(𝑛
1
), 𝑥
2
(𝑛
2
), · · · , 𝑥
𝑚
(𝑛
𝑚
)
重数之和自然为 𝑛, 则差分方程的通解为
𝑦
𝑗
=
𝑛
1
Õ
𝑖=1
𝛼
(𝑖)
1
𝑗
𝑖−1
𝑥
𝑗
1
+
𝑛
2
Õ
𝑖=1
𝛼
(𝑖)
2
𝑗
𝑖−1
𝑥
𝑗
2
+ · · · +
𝑛
𝑚
Õ
𝑖=1
𝛼
(𝑖)
𝑚
𝑗
𝑖−1
𝑥
𝑗
𝑚
其中 𝛼
(𝑖)
𝑘
由初始条件 𝑦
0
, 𝑦
1
, · · · , 𝑦
𝑛−1
唯一确定
6.3 差分方法的相容性
一阶常微分方程的初值问题
𝑢
′
= 𝑓 (𝑥, 𝑢), 𝑥 ∈ [𝑎, 𝑏]
𝑢(𝑎) =
˜
𝑢
𝐾 多步法的一般形式为
𝛼
0
𝑢
𝑛
+ · · · + 𝛼
𝐾
𝑢
𝑛+𝐾
= ℎ
[
𝛽
0
𝑓 (𝑥
𝑛
, 𝑢
𝑛
) + · · · + 𝛽
𝐾
𝑓 (𝑥
𝑛+𝐾
, 𝑢
𝑛+𝐾
)
]
其中 𝛼
𝐾
≠ 0, 𝛼
0
与 𝛽
0
不全为零. 若其局部截断误差满足
𝑇
𝑛
= 𝑂 (ℎ
𝑝
)
因而相容的多步法至少是一阶精度的. 对于线性 𝐾 步法
𝐾
Õ
𝑗=0
𝛼
𝑗
𝑢
𝑛+ 𝑗
= ℎ
𝐾
Õ
𝑗=0
𝛽
𝑗
𝑓 (𝑥
𝑛+ 𝑗
, 𝑢
𝑛+ 𝑗
)
考虑两个特征多项式
𝜌(𝑥) =
𝐾
Õ
𝑗
=
0
𝛼
𝑗
𝑥
𝑗
, 𝜎(𝑥) =
𝐾
Õ
𝑗
=
0
𝛽
𝑗
𝑥
𝑗
它的相容性条件为
𝜌(1) = 0, 𝜌
′
(1) = 𝜎(1) = 1
6.4 差分方法的收敛性
𝐾 步法需要 𝐾 个初始值
𝑥
𝑗
= 𝑎 + 𝑗 ℎ, 𝑗 = 0, 1, · · · , 𝐾 − 1
ℎ → 0 时初始值
˜
𝑢
𝑗
趋于
˜
𝑢. 若此时线性 𝐾 步法的解 𝑢
𝑛
收敛于微分方程的解 𝑢(𝑥
𝑛
), 则称该方法是收敛的
线性 𝐾 步法收敛的充要条件是其特征方程
𝑎
0
+ 𝑎
1
𝑥 + 𝑎
2
𝑥
2
+ · · · + 𝑎
𝑛
𝑥
𝑛
= 0
的所有根的模都小于等于 1, 并且模为 1 的根是单根
6.5 差分方法的稳定性
与常微分方程相容的线性 𝐾 步法其稳定性与收敛性等价
将数值方法应用于模型常微分方程
𝑢
′
= 𝜆𝑢
𝑢(𝑎) =
˜
𝑢
其中 𝜆 是实部小于零的复数. 若初值发生微小变化, 数值解的变化也是微小的, 则称该数值方法是 A 稳
定的
线性 𝐾 步法绝对稳定的充要条件是
存在左半复平面上的一个区域, 使得 𝜆ℎ 取自该区域时, 齐次差分方程
𝐾
Õ
𝑗=0
(𝛼
𝑗
− 𝜆ℎ)𝑦
𝑛+ 𝑗
= 0
其特征方程
𝐾
Õ
𝑗=0
(𝛼
𝑗
− 𝜆ℎ)𝑥
𝑗
= 0
根的模小于 1
7 常微分方程组的数值解法
对于常微分方程组
𝑢
′
1
= 𝑓
1
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
𝑢
′
2
= 𝑓
2
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
· · ·
𝑢
′
𝑚
= 𝑓
𝑚
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
𝑢
1
(𝑎) = 𝜂
1
𝑢
2
(𝑎) = 𝜂
2
· · ·
𝑢
𝑚
(𝑎) = 𝜂
𝑚
可以使用向量的形式. 令
𝑈 =
𝑢
1
𝑢
2
.
.
.
𝑢
𝑚
, 𝐹(𝑥, 𝑈) =
𝑓
1
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
𝑓
2
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
.
.
.
𝑓
𝑚
(𝑥, 𝑢
1
, 𝑢
2
, · · · , 𝑢
𝑚
)
, 𝜂 =
𝜂
1
𝜂
2
.
.
.
𝜂
𝑚
那么微分方程组可以写为
𝑈
′
= 𝐹 (𝑥, 𝑈)
𝑈(𝑎) = 𝜂
可以使用上述数值方法求解
对于高次的常微分方程可以转化为一阶常微分方程组, 如
𝑢
′′
= 𝑓 (𝑥, 𝑢, 𝑢
′
), 𝑥 ∈ [𝑎, 𝑏]
𝑢
′
(𝑎) = 𝑧
0
𝑢(𝑎) = 𝑦
0
令 𝑧 = 𝑢
′
, 则
𝑢
′
= 𝑧
𝑧
′
= 𝑓 (𝑥, 𝑢, 𝑧)
𝑢(𝑎) = 𝑦
0
𝑧(𝑎) = 𝑧
0
