Runge-Kutta 方法简介
Euler 方法及其改进形式都是在区间 [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ] y n + 1 y_{n+1} y n + 1
RK 方法的核心思想是在 [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ] k i k_i k i h h h
二阶 Runge-Kutta 方法
数学基础
二阶 RK 方法的一般计算格式为：
{ k 1 = h f ( x n , y n ) k 2 = h f ( x n + α h , y n + β k 1 ) y n + 1 = y n + ω 1 k 1 + ω 2 k 2 \begin{cases}
k_1 = h f(x_n, y_n) \\
k_2 = h f(x_n + \alpha h, y_n + \beta k_1) \\
y_{n+1} = y_n + \omega_1 k_1 + \omega_2 k_2
\end{cases} ⎩ ⎨ ⎧ k 1 = h f ( x n , y n ) k 2 = h f ( x n + α h , y n + β k 1 ) y n + 1 = y n + ω 1 k 1 + ω 2 k 2 其中 ω 1 , ω 2 , α , β \omega_1, \omega_2, \alpha, \beta ω 1 , ω 2 , α , β
系数推导 ：
类似 改进 Euler 中的推导过程，
记 k ~ 1 = h f ( x n , y ( x n ) ) , k ~ 2 = h f ( x n + α h , y ( x n ) + β k ~ 1 ) \tilde{k}_1 = hf(x_n, y(x_n)),\;\tilde{k}_2 = hf(x_n+\alpha h, y(x_n)+\beta \tilde{k}_1) k ~ 1 = h f ( x n , y ( x n )) , k ~ 2 = h f ( x n + α h , y ( x n ) + β k ~ 1 )
y ~ n + 1 = y ( x n ) + ω 1 k ~ 1 + ω 2 k ~ 2 = y ( x n ) + ω 1 k ~ 1 + ω 2 h ( f + f x α h + f y β k ~ 1 ) = y ( x n ) + h ( ω 1 + ω 2 ) f + h 2 ω 2 ( α f x + β f y f ) + O ( h 3 ) \begin{aligned}
\tilde{y}_{n+1}
&= y(x_n) + \omega_1 \tilde{k}_1 + \omega_2 \tilde{k}_2 \\
&= y(x_n) + \omega_1 \tilde{k}_1 + \omega_2 h (f + f_x \alpha h + f_y \beta \tilde{k}_1) \\
&= y(x_n) + h (\omega_1+\omega_2) f + h^2 \omega_2 (\alpha f_x + \beta f_y f) + O(h^3)
\end{aligned} y ~ n + 1 = y ( x n ) + ω 1 k ~ 1 + ω 2 k ~ 2 = y ( x n ) + ω 1 k ~ 1 + ω 2 h ( f + f x α h + f y β k ~ 1 ) = y ( x n ) + h ( ω 1 + ω 2 ) f + h 2 ω 2 ( α f x + β f y f ) + O ( h 3 ) 又因为：
y ( x n + 1 ) = y ( x n ) + h f + h 2 f x + f y f 2 + O ( h 3 ) y(x_{n+1}) = y(x_n) + h f + h^2 \frac{f_x + f_y f}{2} + O(h^3) y ( x n + 1 ) = y ( x n ) + h f + h 2 2 f x + f y f + O ( h 3 ) 希望达到 2 阶精度，匹配系数得：
{ ω 1 + ω 2 = 1 α ω 2 = 1 2 β ω 2 = 1 2 \begin{cases}
\omega_1 + \omega_2 = 1 \\
\alpha \omega_2 = \frac{1}{2} \\
\beta \omega_2 = \frac{1}{2}
\end{cases} ⎩ ⎨ ⎧ ω 1 + ω 2 = 1 α ω 2 = 2 1 β ω 2 = 2 1 方程个数小于未知数个数，解不唯一。
注 ：一般会取 α = β \alpha = \beta α = β k 2 k_2 k 2 x x x y y y
计算公式
取法一 （改进 Euler 法）：
ω 1 = 1 2 , ω 2 = 1 2 , α = 1 , β = 1 \omega_1 = \frac{1}{2}, \quad \omega_2 = \frac{1}{2}, \quad \alpha = 1, \quad \beta = 1 ω 1 = 2 1 , ω 2 = 2 1 , α = 1 , β = 1 { k 1 = h f ( x n , y n ) k 2 = h f ( x n + h , y n + k 1 ) y n + 1 = y n + 1 2 ( k 1 + k 2 ) \begin{cases}
k_1 = h f(x_n, y_n) \\
k_2 = h f(x_n + h, y_n + k_1) \\
y_{n+1} = y_n + \frac{1}{2}(k_1 + k_2)
\end{cases} ⎩ ⎨ ⎧ k 1 = h f ( x n , y n ) k 2 = h f ( x n + h , y n + k 1 ) y n + 1 = y n + 2 1 ( k 1 + k 2 ) 取法二 （中点方法）：
ω 1 = 0 , ω 2 = 1 , α = 1 2 , β = 1 2 \omega_1 = 0, \quad \omega_2 = 1, \quad \alpha = \frac{1}{2}, \quad \beta = \frac{1}{2} ω 1 = 0 , ω 2 = 1 , α = 2 1 , β = 2 1 { k 1 = h f ( x n , y n ) k 2 = h f ( x n + h 2 , y n + k 1 2 ) y n + 1 = y n + k 2 \begin{cases}
k_1 = h f(x_n, y_n) \\
k_2 = h f\left(x_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right) \\
y_{n+1} = y_n + k_2
\end{cases} ⎩ ⎨ ⎧ k 1 = h f ( x n , y n ) k 2 = h f ( x n + 2 h , y n + 2 k 1 ) y n + 1 = y n + k 2 Butcher 表 ：
0 α β ω 1 ω 2 \begin{array}{c|cc}
0 & & \\
\alpha & \beta & \\
\hline
& \omega_1 & \omega_2
\end{array} 0 α β ω 1 ω 2 改进 Euler：
0 1 1 1 / 2 1 / 2 \begin{array}{c|cc}
0 & & \\
1 & 1 & \\
\hline
& 1/2 & 1/2
\end{array} 0 1 1 1/2 1/2 中点方法：
0 1 / 2 1 / 2 0 1 \begin{array}{c|cc}
0 & & \\
1/2 & 1/2 & \\
\hline
& 0 & 1
\end{array} 0 1/2 1/2 0 1 算法实现
中点方法 ：
def rk2_midpoint (f, a, b, y0, N): """ 二阶 Runge-Kutta 中点方法 """ h = (b - a) / (N - 1 ) hh = h * 0.5 x = a y = [ 0.0 ] * N y[ 0 ] = y0 for i in range ( 1 , N): k1 = h * f(x, y[i - 1 ]) k2 = h * f(x + hh, y[i - 1 ] + k1 * 0.5 ) y[i] = y[i - 1 ] + k2 x += h return y 稳定性分析
对模型方程 y ′ = λ y ( λ ∈ C ) y' = \lambda y\ (\lambda \in \mathbb{C}) y ′ = λ y ( λ ∈ C )
y n + 1 = ( 1 + λ h + 1 2 λ 2 h 2 ) y n y_{n+1} = (1 + \lambda h + \frac{1}{2}\lambda^2 h^2) y_n y n + 1 = ( 1 + λh + 2 1 λ 2 h 2 ) y n 故绝对稳定条件 为 ∣ 1 + λ h + 1 2 λ 2 h 2 ∣ < 1 |1 + \lambda h + \frac{1}{2}\lambda^2 h^2| < 1 ∣1 + λh + 2 1 λ 2 h 2 ∣ < 1
当 λ ∈ R \lambda \in \mathbb{R} λ ∈ R λ < 0 \lambda < 0 λ < 0
h < − 2 λ h < -\frac{2}{\lambda} h < − λ 2 四阶 Runge-Kutta 方法
数学基础
四阶 RK 的一般形式含 4 个斜率 k 1 , k 2 , k 3 , k 4 k_1, k_2, k_3, k_4 k 1 , k 2 , k 3 , k 4
{ k 1 = h f ( x n , y n ) k 2 = h f ( x n + α 1 h , y n + β 1 k 1 ) k 3 = h f ( x n + α 2 h , y n + β 2 k 1 + γ 1 k 2 ) k 4 = h f ( x n + α 3 h , y n + β 3 k 1 + γ 2 k 2 + ζ 1 k 3 ) y n + 1 = y n + ω 1 k 1 + ω 2 k 2 + ω 3 k 3 + ω 4 k 4 \begin{cases}
k_1 = h f(x_n, y_n) \\
k_2 = h f(x_n + \alpha_1 h, y_n + \beta_1 k_1) \\
k_3 = h f(x_n + \alpha_2 h, y_n + \beta_2 k_1 + \gamma_1 k_2) \\
k_4 = h f(x_n + \alpha_3 h, y_n + \beta_3 k_1 + \gamma_2 k_2 + \zeta_1 k_3) \\
y_{n+1} = y_n + \omega_1 k_1 + \omega_2 k_2 + \omega_3 k_3 + \omega_4 k_4
\end{cases} ⎩ ⎨ ⎧ k 1 = h f ( x n , y n ) k 2 = h f ( x n + α 1 h , y n + β 1 k 1 ) k 3 = h f ( x n + α 2 h , y n + β 2 k 1 + γ 1 k 2 ) k 4 = h f ( x n + α 3 h , y n + β 3 k 1 + γ 2 k 2 + ζ 1 k 3 ) y n + 1 = y n + ω 1 k 1 + ω 2 k 2 + ω 3 k 3 + ω 4 k 4 希望局部截断误差为 O ( h 5 ) O(h^5) O ( h 5 )
计算公式
取满足方程组的一组经典解，得到标准四阶 RK 方法 ：
{ k 1 = h f ( x n , y n ) k 2 = h f ( x n + h 2 , y n + k 1 2 ) k 3 = h f ( x n + h 2 , y n + k 2 2 ) k 4 = h f ( x n + h , y n + k 3 ) y n + 1 = y n + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) \begin{cases}
k_1 = h f(x_n, y_n) \\
k_2 = h f\left(x_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right) \\
k_3 = h f\left(x_n + \frac{h}{2}, y_n + \frac{k_2}{2}\right) \\
k_4 = h f(x_n + h, y_n + k_3) \\
y_{n+1} = y_n + \frac{1}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)
\end{cases} ⎩ ⎨ ⎧ k 1 = h f ( x n , y n ) k 2 = h f ( x n + 2 h , y n + 2 k 1 ) k 3 = h f ( x n + 2 h , y n + 2 k 2 ) k 4 = h f ( x n + h , y n + k 3 ) y n + 1 = y n + 6 1 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) Butcher 表 ：
0 1 / 2 1 / 2 1 / 2 0 1 / 2 1 0 0 1 1 / 6 1 / 3 1 / 3 1 / 6 \begin{array}{c|cccc}
0 & & & & \\
1/2 & 1/2 & & & \\
1/2 & 0 & 1/2 & & \\
1 & 0 & 0 & 1 & \\
\hline
& 1/6 & 1/3 & 1/3 & 1/6
\end{array} 0 1/2 1/2 1 1/2 0 0 1/6 1/2 0 1/3 1 1/3 1/6 注 ：
虽然 k 2 , k 3 k_2, k_3 k 2 , k 3 k 3 k_3 k 3 k 2 k_2 k 2 k 1 k_1 k 1
算法实现
def rk4 (f, a, b, y0, N): """ 经典四阶 Runge-Kutta 方法 """ h = (b - a) / (N - 1 ) hh = h * 0.5 x = a y = [ 0.0 ] * N y[ 0 ] = y0 for i in range ( 1 , N): k1 = h * f(x, y[i - 1 ]) k2 = h * f(x + hh, y[i - 1 ] + k1 * 0.5 ) k3 = h * f(x + hh, y[i - 1 ] + k2 * 0.5 ) k4 = h * f(x + h, y[i - 1 ] + k3) y[i] = y[i - 1 ] + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0 x += h return y 稳定性分析
对模型方程 y ′ = λ y y' = \lambda y y ′ = λ y
R ( z ) = 1 + z + z 2 2 + z 3 6 + z 4 24 , z = λ h R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}, \quad z = \lambda h R ( z ) = 1 + z + 2 z 2 + 6 z 3 + 24 z 4 , z = λh 故绝对稳定条件为：∣ R ( z ) ∣ ≤ 1 |R(z)| \leq 1 ∣ R ( z ) ∣ ≤ 1
当 λ ∈ R \lambda \in \mathbb{R} λ ∈ R λ < 0 \lambda < 0 λ < 0
h < − 2.78 λ h < -\frac{2.78}{\lambda} h < − λ 2.78 注 ：四阶 RK 的绝对稳定区间比二阶 RK 略大，但二者均为条件稳定，对刚性方程（∣ λ ∣ ≫ 0 |\lambda| \gg 0 ∣ λ ∣ ≫ 0 h h h
线性多步法简介
一般形式
若 y n + 1 y_{n+1} y n + 1 y n y_n y n y n , y n − 1 , … , y n − k y_n, y_{n-1}, \dots, y_{n-k} y n , y n − 1 , … , y n − k 多步法 。其一般线性形式为
y n + 1 = ∑ i = 0 k α i y n − i + h ∑ i = − 1 k β i f ( x n − i , y n − i ) y_{n+1} = \sum_{i=0}^{k} \alpha_i y_{n-i} + h \sum_{i=-1}^{k} \beta_i f(x_{n-i}, y_{n-i}) y n + 1 = i = 0 ∑ k α i y n − i + h i = − 1 ∑ k β i f ( x n − i , y n − i )
若 β − 1 = 0 \beta_{-1} = 0 β − 1 = 0 f ( x n + 1 , y n + 1 ) f(x_{n+1}, y_{n+1}) f ( x n + 1 , y n + 1 )
若 β − 1 ≠ 0 \beta_{-1} \neq 0 β − 1 = 0
构造原理
基于恒等式：
y ( x n + 1 ) − y ( x n − k ) = ∫ x n − k x n + 1 f ( ξ , y ( ξ ) ) d ξ y(x_{n+1}) - y(x_{n-k}) = \int_{x_{n-k}}^{x_{n+1}} f(\xi, y(\xi))\,\mathrm{d}\xi y ( x n + 1 ) − y ( x n − k ) = ∫ x n − k x n + 1 f ( ξ , y ( ξ )) d ξ 令 F ( x ) ≡ f ( x , y ( x ) ) F(x) \equiv f(x, y(x)) F ( x ) ≡ f ( x , y ( x )) ψ ( x ) \psi(x) ψ ( x ) F ( x ) F(x) F ( x )
注 ：当 k = 0 k=0 k = 0
隐式方法的局部截断误差简化
对隐式多步法（β − 1 ≠ 0 \beta_{-1} \neq 0 β − 1 = 0 y ~ n + 1 \tilde{y}_{n+1} y ~ n + 1
y ~ n + 1 = ∑ i = 0 k α i y ( x n − i ) + h ∑ i = 0 k β i y ′ ( x n − i ) + h β − 1 f ( x n + 1 , y ~ n + 1 ) \tilde{y}_{n+1} = \sum_{i=0}^{k}\alpha_i y(x_{n-i}) + h\sum_{i=0}^{k}\beta_i y'(x_{n-i}) + h\beta_{-1}f(x_{n+1},\tilde{y}_{n+1}) y ~ n + 1 = i = 0 ∑ k α i y ( x n − i ) + h i = 0 ∑ k β i y ′ ( x n − i ) + h β − 1 f ( x n + 1 , y ~ n + 1 ) 若将隐式项中的 y ~ n + 1 \tilde{y}_{n+1} y ~ n + 1 y ( x n + 1 ) y(x_{n+1}) y ( x n + 1 )
T n + 1 = y ( x n + 1 ) − [ ∑ i = 0 k α i y ( x n − i ) + h ∑ i = − 1 k β i y ′ ( x n − i ) ] T_{n+1} = y(x_{n+1}) - \left[\sum_{i=0}^{k}\alpha_i y(x_{n-i}) + h\sum_{i=-1}^{k}\beta_i y'(x_{n-i})\right] T n + 1 = y ( x n + 1 ) − [ i = 0 ∑ k α i y ( x n − i ) + h i = − 1 ∑ k β i y ′ ( x n − i ) ] 由微分中值定理：
f ( x n + 1 , y ~ n + 1 ) − f ( x n + 1 , y ( x n + 1 ) ) = − f y ⋅ l n + 1 f(x_{n+1},\tilde{y}_{n+1}) - f(x_{n+1},y(x_{n+1})) = -f_y \cdot l_{n+1} f ( x n + 1 , y ~ n + 1 ) − f ( x n + 1 , y ( x n + 1 )) = − f y ⋅ l n + 1 故在 p 阶精度时，隐式修正项为 h β − 1 f y ⋅ l n + 1 = O ( h p + 2 ) h\beta_{-1}f_y \cdot l_{n+1} = O(h^{p+2}) h β − 1 f y ⋅ l n + 1 = O ( h p + 2 )
l n + 1 = T n + 1 + O ( h p + 2 ) l_{n+1} = T_{n+1} + O(h^{p+2}) l n + 1 = T n + 1 + O ( h p + 2 ) 也即隐式项的修正不影响精度结果。实际推导局部截断误差时，可直接将隐式项中的 y n + 1 y_{n+1} y n + 1 y ( x n + 1 ) y(x_{n+1}) y ( x n + 1 )
出发值计算
线性多步法需 k k k y 0 , y 1 , … , y k − 1 y_0, y_1, \dots, y_{k-1} y 0 , y 1 , … , y k − 1 单步法 （如四阶 RK）以小步长计算出发值，且单步法阶数应不低于多步法阶数，以保证整体精度不被污染。
Adams 外推公式
思想 ：
插值节点取积分区间左侧的 x n , x n − 1 , … , x n − k x_n, x_{n-1}, \dots, x_{n-k} x n , x n − 1 , … , x n − k [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ] 外推 积分，得到显式 格式。
二阶 Adams-Bashforth（两点线性插值）
对 F ( x ) F(x) F ( x ) x n , x n − 1 x_n, x_{n-1} x n , x n − 1
ψ 1 ( x ) = x − x n − 1 h F ( x n ) − x − x n h F ( x n − 1 ) \psi_1(x) = \frac{x - x_{n-1}}{h} F(x_n) - \frac{x - x_n}{h} F(x_{n-1}) ψ 1 ( x ) = h x − x n − 1 F ( x n ) − h x − x n F ( x n − 1 ) 余项：
R 1 ( x ) = 1 2 ( x − x n ) ( x − x n − 1 ) F ′ ′ ( η ) , η ∈ ( x n − 1 , x n ) R_1(x) = \frac{1}{2}(x - x_n)(x - x_{n-1}) F''(\eta), \quad \eta \in (x_{n-1}, x_n) R 1 ( x ) = 2 1 ( x − x n ) ( x − x n − 1 ) F ′′ ( η ) , η ∈ ( x n − 1 , x n ) 在 [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ] h h h
∫ x n x n + 1 ψ 1 ( ξ ) d ξ = h 2 [ 3 F ( x n ) − F ( x n − 1 ) ] \int_{x_n}^{x_{n+1}} \psi_1(\xi)\,\mathrm{d}\xi = \frac{h}{2}\left[3F(x_n) - F(x_{n-1})\right] ∫ x n x n + 1 ψ 1 ( ξ ) d ξ = 2 h [ 3 F ( x n ) − F ( x n − 1 ) ] 得到计算格式 ：
y n + 1 = y n + h 2 [ 3 f ( x n , y n ) − f ( x n − 1 , y n − 1 ) ] y_{n+1} = y_n + \frac{h}{2}\left[3f(x_n, y_n) - f(x_{n-1}, y_{n-1})\right] y n + 1 = y n + 2 h [ 3 f ( x n , y n ) − f ( x n − 1 , y n − 1 ) ] 局部截断误差 ：
l n + 1 = ∫ x n x n + 1 R 1 ( ξ ) d ξ = 5 12 h 3 F ′ ′ ( η ) l_{n+1} = \int_{x_n}^{x_{n+1}} R_1(\xi)\,\mathrm{d}\xi = \frac{5}{12} h^3 F''(\eta) l n + 1 = ∫ x n x n + 1 R 1 ( ξ ) d ξ = 12 5 h 3 F ′′ ( η ) 四阶 Adams-Bashforth（四点三次插值）
对 F ( x ) F(x) F ( x ) x n , x n − 1 , x n − 2 , x n − 3 x_n, x_{n-1}, x_{n-2}, x_{n-3} x n , x n − 1 , x n − 2 , x n − 3 ψ 3 ( x ) \psi_3(x) ψ 3 ( x ) [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ]
计算格式 ：
y n + 1 = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) y_{n+1} = y_n + \frac{h}{24}\left(55f_n - 59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) y n + 1 = y n + 24 h ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) 其中 f i ≡ f ( x i , y i ) f_i \equiv f(x_i, y_i) f i ≡ f ( x i , y i )
局部截断误差 ：
l n + 1 = 251 720 h 5 F ( 4 ) ( η ) + O ( h 6 ) , η ∈ ( x n − 3 , x n ) l_{n+1} = \frac{251}{720} h^5 F^{(4)}(\eta) + O(h^6), \quad \eta \in (x_{n-3}, x_n) l n + 1 = 720 251 h 5 F ( 4 ) ( η ) + O ( h 6 ) , η ∈ ( x n − 3 , x n ) 从而四阶 Adams-Bashforth 方法具有 4 阶精度。
Adams 内插公式
思想 ：
插值节点包含积分区间右端点，取为 x n + 1 , x n , … , x n − k x_{n+1}, x_n, \dots, x_{n-k} x n + 1 , x n , … , x n − k [ x n , x n + 1 ] [x_n, x_{n+1}] [ x n , x n + 1 ] 内插 积分，得到隐式 格式。
二阶 Adams-Moulton（梯形公式）
对 F ( x ) F(x) F ( x ) x n , x n + 1 x_n, x_{n+1} x n , x n + 1 计算格式 为：
y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] y_{n+1} = y_n + \frac{h}{2}\left[f(x_n, y_n) + f(x_{n+1}, y_{n+1})\right] y n + 1 = y n + 2 h [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] 局部截断误差 ：
l n + 1 = − 1 12 h 3 F ′ ′ ( η ) + O ( h 4 ) , η ∈ ( x n , x n + 1 ) l_{n+1} = -\frac{1}{12} h^3 F''(\eta) + O(h^4), \quad \eta \in (x_{n}, x_{n+1}) l n + 1 = − 12 1 h 3 F ′′ ( η ) + O ( h 4 ) , η ∈ ( x n , x n + 1 ) 四阶 Adams-Moulton（四点三次插值）
对 F ( x ) F(x) F ( x ) x n + 1 , x n , x n − 1 , x n − 2 x_{n+1}, x_n, x_{n-1}, x_{n-2} x n + 1 , x n , x n − 1 , x n − 2
计算格式 ：
y n + 1 = y n + h 24 ( 9 f n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) y_{n+1} = y_n + \frac{h}{24}\left(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}\right) y n + 1 = y n + 24 h ( 9 f n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) 局部截断误差 ：
l n + 1 = − 19 720 h 5 F ( 4 ) ( η ) + O ( h 6 ) , η ∈ ( x n − 2 , x n + 1 ) l_{n+1} = -\frac{19}{720} h^5 F^{(4)}(\eta) + O(h^6), \quad \eta \in (x_{n-2}, x_{n+1}) l n + 1 = − 720 19 h 5 F ( 4 ) ( η ) + O ( h 6 ) , η ∈ ( x n − 2 , x n + 1 ) 结论 ：四阶 Adams-Moulton 为隐式方法，局部截断误差系数小于同阶 Adams-Bashforth，精度更优且稳定性更好。
预估-校正法
基本思想
预估（Predictor, P） ：用显式方法（如 Adams-Bashforth）计算 y n + 1 y_{n+1} y n + 1 p n + 1 p_{n+1} p n + 1 求值（Evaluation, E） ：计算 f ( x n + 1 , p n + 1 ) f(x_{n+1}, p_{n+1}) f ( x n + 1 , p n + 1 ) 校正（Corrector, C） ：用隐式方法（如 Adams-Moulton）迭代求解 y n + 1 y_{n+1} y n + 1
改进 Euler 的 PEC 形式
P : p n + 1 = y n + h f ( x n , y n ) E : f ~ n + 1 = f ( x n + 1 , p n + 1 ) C : y n + 1 = y n + h 2 ( f n + f ~ n + 1 ) \begin{aligned}
\mathrm{P}: &\quad p_{n+1} = y_n + h f(x_n, y_n) \\
\mathrm{E}: &\quad \tilde{f}_{n+1} = f(x_{n+1}, p_{n+1}) \\
\mathrm{C}: &\quad y_{n+1} = y_n + \frac{h}{2}\left(f_n + \tilde{f}_{n+1}\right)
\end{aligned} P : E : C : p n + 1 = y n + h f ( x n , y n ) f ~ n + 1 = f ( x n + 1 , p n + 1 ) y n + 1 = y n + 2 h ( f n + f ~ n + 1 ) 注 ：此即改进 Euler 法的实现形式，仅进行一次校正（迭代次数为 1）。
Adams 预估-校正法
PEC 格式 （四阶）：
P : p n + 1 = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) E : f ~ n + 1 = f ( x n + 1 , p n + 1 ) C : y n + 1 = y n + h 24 ( 9 f ~ n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) \begin{aligned}
\mathrm{P}: &\quad p_{n+1} = y_n + \frac{h}{24}\left(55f_n - 59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\
\mathrm{E}: &\quad \tilde{f}_{n+1} = f(x_{n+1}, p_{n+1}) \\
\mathrm{C}: &\quad y_{n+1} = y_n + \frac{h}{24}\left(9\tilde{f}_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}\right)
\end{aligned} P : E : C : p n + 1 = y n + 24 h ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) f ~ n + 1 = f ( x n + 1 , p n + 1 ) y n + 1 = y n + 24 h ( 9 f ~ n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) P(EC)2 ^2 2 （四阶，两次校正）：
P : p n + 1 = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) E : f ~ n + 1 ( 1 ) = f ( x n + 1 , p n + 1 ) C : y n + 1 ( 1 ) = y n + h 24 ( 9 f ~ n + 1 ( 1 ) + 19 f n − 5 f n − 1 + f n − 2 ) E : f ~ n + 1 ( 2 ) = f ( x n + 1 , y n + 1 ( 1 ) ) C : y n + 1 ( 2 ) = y n + h 24 ( 9 f ~ n + 1 ( 2 ) + 19 f n − 5 f n − 1 + f n − 2 ) \begin{aligned}
\mathrm{P}: &\quad p_{n+1} = y_n + \frac{h}{24}\left(55f_n - 59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\
\mathrm{E}: &\quad \tilde{f}_{n+1}^{(1)} = f(x_{n+1}, p_{n+1}) \\
\mathrm{C}: &\quad y_{n+1}^{(1)} = y_n + \frac{h}{24}\left(9\tilde{f}_{n+1}^{(1)} + 19f_n - 5f_{n-1} + f_{n-2}\right) \\
\mathrm{E}: &\quad \tilde{f}_{n+1}^{(2)} = f(x_{n+1}, y_{n+1}^{(1)}) \\
\mathrm{C}: &\quad y_{n+1}^{(2)} = y_n + \frac{h}{24}\left(9\tilde{f}_{n+1}^{(2)} + 19f_n - 5f_{n-1} + f_{n-2}\right)
\end{aligned} P : E : C : E : C : p n + 1 = y n + 24 h ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) f ~ n + 1 ( 1 ) = f ( x n + 1 , p n + 1 ) y n + 1 ( 1 ) = y n + 24 h ( 9 f ~ n + 1 ( 1 ) + 19 f n − 5 f n − 1 + f n − 2 ) f ~ n + 1 ( 2 ) = f ( x n + 1 , y n + 1 ( 1 ) ) y n + 1 ( 2 ) = y n + 24 h ( 9 f ~ n + 1 ( 2 ) + 19 f n − 5 f n − 1 + f n − 2 ) 注 ：校正公式的截断误差优于预估公式。实际计算中通常限制校正迭代次数为 1–2 次，以平衡计算开销与精度收益。
误差修正与 PMECME 方法
利用 Adams-Bashforth（预估）与 Adams-Moulton（校正）的同阶截断误差进行事后误差估计与修正 ：
设 p n + 1 p_{n+1} p n + 1 c n + 1 c_{n+1} c n + 1
y ( x n + 1 ) − p n + 1 = 251 720 h 5 y ( 5 ) ( ξ n ) , x n − 3 < ξ n < x n + 1 y ( x n + 1 ) − c n + 1 = − 19 720 h 5 y ( 5 ) ( η n ) , x n − 2 < η n < x n + 1 \begin{aligned}
y(x_{n+1}) - p_{n+1} &= \frac{251}{720} h^5 y^{(5)}(\xi_n), \quad x_{n-3} < \xi_n < x_{n+1} \\
y(x_{n+1}) - c_{n+1} &= -\frac{19}{720} h^5 y^{(5)}(\eta_n), \quad x_{n-2} < \eta_n < x_{n+1}
\end{aligned} y ( x n + 1 ) − p n + 1 y ( x n + 1 ) − c n + 1 = 720 251 h 5 y ( 5 ) ( ξ n ) , x n − 3 < ξ n < x n + 1 = − 720 19 h 5 y ( 5 ) ( η n ) , x n − 2 < η n < x n + 1 两式相减，由中值定理得：
c n + 1 − p n + 1 = 270 720 h 5 y ( 5 ) ( ζ n ) c_{n+1} - p_{n+1} = \frac{270}{720} h^5 y^{(5)}(\zeta_n) c n + 1 − p n + 1 = 720 270 h 5 y ( 5 ) ( ζ n ) 从而得到事后误差估计 ：
y ( x n + 1 ) − p n + 1 ≈ 251 270 ( c n + 1 − p n + 1 ) y ( x n + 1 ) − c n + 1 ≈ − 19 270 ( c n + 1 − p n + 1 ) \begin{aligned}
y(x_{n+1}) - p_{n+1} &\approx \frac{251}{270}(c_{n+1} - p_{n+1}) \\
y(x_{n+1}) - c_{n+1} &\approx -\frac{19}{270}(c_{n+1} - p_{n+1})
\end{aligned} y ( x n + 1 ) − p n + 1 y ( x n + 1 ) − c n + 1 ≈ 270 251 ( c n + 1 − p n + 1 ) ≈ − 270 19 ( c n + 1 − p n + 1 ) PMECME 修正方案 ：
P : p n + 1 = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) M : m n + 1 = p n + 1 + 251 270 ( c n − p n ) E : f n + 1 = f ( x n + 1 , m n + 1 ) C : c n + 1 = y n + h 24 ( 9 f n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) M : y n + 1 = c n + 1 − 19 270 ( c n + 1 − p n + 1 ) E : f n + 1 = f ( x n + 1 , y n + 1 ) \begin{aligned}
\mathrm{P}: &\quad p_{n+1} = y_n + \frac{h}{24}\left(55f_n - 59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\
\mathrm{M}: &\quad m_{n+1} = p_{n+1} + \frac{251}{270}(c_n - p_n) \\
\mathrm{E}: &\quad f_{n+1} = f(x_{n+1}, m_{n+1}) \\
\mathrm{C}: &\quad c_{n+1} = y_n + \frac{h}{24}\left(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}\right) \\
\mathrm{M}: &\quad y_{n+1} = c_{n+1} - \frac{19}{270}(c_{n+1} - p_{n+1}) \\
\mathrm{E}: &\quad f_{n+1} = f(x_{n+1}, y_{n+1})
\end{aligned} P : M : E : C : M : E : p n + 1 = y n + 24 h ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) m n + 1 = p n + 1 + 270 251 ( c n − p n ) f n + 1 = f ( x n + 1 , m n + 1 ) c n + 1 = y n + 24 h ( 9 f n + 1 + 19 f n − 5 f n − 1 + f n − 2 ) y n + 1 = c n + 1 − 270 19 ( c n + 1 − p n + 1 ) f n + 1 = f ( x n + 1 , y n + 1 ) 其中两步 M \mathrm{M} M
注 ：PMECME 可以在不显著增加计算量的前提下将有效精度提升约 1–2 阶。
常微分方程组与高阶微分方程
常微分方程组
问题形式 ：
{ d y d x = f ( x , y , z ) d z d x = g ( x , y , z ) y ( x 0 ) = y 0 , z ( x 0 ) = z 0 \begin{cases}
\displaystyle\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y, z) \\[6pt]
\displaystyle\frac{\mathrm{d}z}{\mathrm{d}x} = g(x, y, z) \\[6pt]
y(x_0)=y_0,\ z(x_0)=z_0
\end{cases} ⎩ ⎨ ⎧ d x d y = f ( x , y , z ) d x d z = g ( x , y , z ) y ( x 0 ) = y 0 , z ( x 0 ) = z 0 数值处理 ：将各未知量视为向量 Y = [ y , z ] T \mathbf{Y} = [y, z]^T Y = [ y , z ] T F = [ f , g ] T \mathbf{F} = [f, g]^T F = [ f , g ] T
注 ：方程组求解时，各分量必须同步推进，即先计算所有分量的斜率，再统一更新，避免使用已更新的分量去计算其他分量（这会改变方法的数学性质）。
高阶微分方程
降阶法 ：m m m m − 1 m-1 m − 1
例 ：二阶方程
d 2 y d x 2 = g ( x , y , y ′ ) , y ( x 0 ) = y 0 , y ′ ( x 0 ) = y 0 ′ \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} = g(x, y, y')
,\ y(x_0)=y_0,\ y'(x_0)=y'_0 d x 2 d 2 y = g ( x , y , y ′ ) , y ( x 0 ) = y 0 , y ′ ( x 0 ) = y 0 ′ 令 z = d y d x z = \dfrac{\mathrm{d}y}{\mathrm{d}x} z = d x d y
{ d y d x = z d z d x = g ( x , y , z ) y ( x 0 ) = y 0 , z ( x 0 ) = y 0 ′ ≡ z 0 \begin{cases}
\displaystyle\frac{\mathrm{d}y}{\mathrm{d}x} = z \\[6pt]
\displaystyle\frac{\mathrm{d}z}{\mathrm{d}x} = g(x, y, z) \\[6pt]
y(x_0)=y_0,\ z(x_0)=y'_0 \equiv z_0
\end{cases} ⎩ ⎨ ⎧ d x d y = z d x d z = g ( x , y , z ) y ( x 0 ) = y 0 , z ( x 0 ) = y 0 ′ ≡ z 0 注 ：高阶 RK 方法（如 RK4）对向量值 F \mathbf{F} F k i k_i k i
复杂度与稳定性总览
方法 类型 局部误差 全局误差 每步 f f f 稳定区间（实 λ < 0 \lambda<0 λ < 0 存储需求 Euler 显式单步 O ( h 2 ) O(h^2) O ( h 2 ) O ( h ) O(h) O ( h ) 1 h < − 2 / λ h < -2/\lambda h < − 2/ λ O ( 1 ) O(1) O ( 1 ) 改进 Euler 显式单步 O ( h 3 ) O(h^3) O ( h 3 ) O ( h 2 ) O(h^2) O ( h 2 ) 2 h < − 2 / λ h < -2/\lambda h < − 2/ λ O ( 1 ) O(1) O ( 1 ) RK2（中点） 显式单步 O ( h 3 ) O(h^3) O ( h 3 ) O ( h 2 ) O(h^2) O ( h 2 ) 2 h < − 2 / λ h < -2/\lambda h < − 2/ λ O ( 1 ) O(1) O ( 1 ) RK4 显式单步 O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 ) 4 h < − 2.78 / λ h < -2.78/\lambda h < − 2.78/ λ O ( 1 ) O(1) O ( 1 ) Adams-Bashforth 4 显式多步 O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 ) 1 较小 O ( 4 ) O(4) O ( 4 ) Adams-Moulton 4 隐式多步 O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 ) 迭代 较大 O ( 4 ) O(4) O ( 4 ) PECE (Adams) 预估-校正 O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 ) 2 中等 O ( 4 ) O(4) O ( 4 ) PMECME 修正预估-校正 O ( h 6 ) O(h^6) O ( h 6 ) O ( h 5 ) O(h^5) O ( h 5 ) 2 中等 O ( 4 ) O(4) O ( 4 )
注 ：
Adams-Moulton 的绝对稳定域显著大于同阶显式方法，适合刚性方程；预估-校正法在保持显式方法计算效率的同时，通过一次校正获得接近隐式的稳定性与精度。