常微分方程的数值求解：步长、精度与稳定性

画面里一条光滑的轨迹，背后可能只是一次次有限步长的推进。步长影响每一步的误差，也决定衰减会不会被算成增长、振子的能量会不会慢慢偏离。双摆与 Rössler 吸引子介绍了各自的模型和 RK4 实现；轨道实验使用 velocity Verlet。判断这些积分方法是否合适，可以先从有精确解的衰减方程和简谐振子入手。

一步究竟在近似什么

初值问题给出变化率和起点：

$$y'(t)=f(t,y(t)),\qquad y(t_0)=y_0.$$

状态 $y$ 可以是一个数，也可以是位置、速度等组成的向量。取步长 $h>0$，在 $t_n=t_0+nh$ 计算近似值 $y_n\approx y(t_n)$。精确的一步满足

$$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y(t))\,dt.$$

困难在于积分里的 $y(t)$ 尚未知。数值方法用已知状态和若干次斜率计算来近似这段增量。绘图时再插值，能让曲线更平滑，却不会提高已经算出的状态精度。这里的误差主要来自时间离散化；浮点舍入是另一层问题，见浮点数与稳定计算。

显式 Euler：从一步误差到整体误差

Driscoll 与 Braun 的《Fundamentals of Numerical Computation》§6.2给出显式 Euler 方法：用起点的斜率走完这一步。

$$y_{n+1}=y_n+h f(t_n,y_n).$$

取一个能直接核对的例子：

$$y'=-2y,\qquad y(0)=1,\qquad y(t)=e^{-2t}.$$

令 $h=0.1$，则 $y_{n+1}=0.8y_n$，所以 $y_n=0.8^n$。

步数 $n$	时间 $t_n$	Euler 近似 $y_n$	精确解 $e^{-2t_n}$
0	0.0	1.000000000	1.000000000
1	0.1	0.800000000	0.818730753
2	0.2	0.640000000	0.670320046
10	1.0	0.107374182	0.135335283

局部一步误差假定这一步从精确状态出发，只看这一步丢掉多少。对光滑的解，Taylor 展开给出

$$y(t_n+h)-\bigl[y(t_n)+hf(t_n,y(t_n))\bigr] =\frac{h^2}{2}y''(t_n)+O(h^3).$$

因此 Euler 的局部一步误差为 $O(h^2)$。在起点，$y''(0)=4$，首项是 $2h^2$；实际绝对误差在 $h=0.1$ 时为 $0.018730753$，在 $h=0.05$ 时为 $0.004837418$，接近缩小到四分之一。教材把上述误差除以 $h$ 后称为局部截断误差，因此书中的这个量是 $O(h)$；比较误差阶时，要先看定义是否除过步长。

整体误差比较实际迭代的 $y_n$ 与精确的 $y(t_n)$，包括此前各步误差的传播。上述教材的收敛定理要求单位局部误差有一致界，且一步增量函数对状态的导数有一致界。在固定的有限时间区间、右端函数与解足够光滑且满足这些条件时，Euler 的整体误差为 $O(h)$。到达同一终点需要约 $1/h$ 步，且旧误差会影响后续斜率，不能只拿最后一步的局部误差代表全程误差。

在 $T=1$，$h=0.1$ 与 $0.05$ 的整体绝对误差分别约为 $0.027961101$ 和 $0.013758629$，比值约为 $2.032$。步长减半后，一阶方法的整体误差通常接近减半。

RK4：多算斜率换取精度

教材 §6.4 的经典四阶 Runge–Kutta 方法每步计算四次斜率。这里的 $k_i$ 表示斜率，尚未乘上 $h$：

$$\begin{aligned} k_1&=f(t_n,y_n),\\ k_2&=f(t_n+h/2,y_n+hk_1/2),\\ k_3&=f(t_n+h/2,y_n+hk_2/2),\\ k_4&=f(t_n+h,y_n+hk_3),\\ y_{n+1}&=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4). \end{aligned}$$

在固定有限区间上，若右端函数与解足够光滑，并满足前述误差传播条件，局部一步误差为 $O(h^5)$，整体误差为 $O(h^4)$。若已进入渐近收敛范围，且舍入误差尚未占主导，步长减半后整体误差应接近原来的 $1/16$。

对 $y'=-2y$，代入四次斜率可得到一步的乘数

$$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},\qquad z=-2h.$$

它与 $e^z$ 的前五项一致。于是 $T=1$ 的近似值可以独立用 $R(-2h)^{1/h}$ 核对。下方程序得到 $h=0.1,0.05,0.025$ 时的终点绝对误差分别为 $4.265\times10^{-6}$、$2.452\times10^{-7}$、$1.470\times10^{-8}$，相邻比值为 $17.396$ 和 $16.682$，正在接近 16。高阶方法仍要满足自己的稳定性条件。

稳定性：衰减为什么会被算成增长

教材 §11.3用测试方程 $y'=\lambda y$ 定义绝对稳定性：固定 $h$，让步数趋于无穷，数值解是否保持有界。设 $z=h\lambda$，显式 Euler 给出

$$y_{n+1}=(1+z)y_n,\qquad y_n=(1+z)^n y_0.$$

所以稳定区域为 $|1+z|\le1$，即复平面中以 $-1$ 为圆心、半径为 1 的闭圆盘。要让衰减模态确实趋于零，需要严格满足 $|1+z|<1$。对实数 $\lambda<0$，这变成

$$0<h<\frac{2}{|\lambda|}.$$

对 $\lambda=-2$，界限是 $h=1$。在界限上，乘数为 $-1$，非零解会来回变号而不衰减；超过界限，振幅开始增长。取 $h=1.1$，从 $y_0=1$ 得到

$$1,\;-1.2,\;1.44,\;-1.728,\;2.0736.$$

真实解始终为正且衰减。即使仍在稳定区域内，当 $0.5<h<1$ 时 Euler 也会变号，说明稳定并不保证定量精确或保持正值。一般的复数 $\lambda$ 必须检查 $|1+h\lambda|$，不能只用 $|\lambda|$ 套实数界限。

RK4 在这个测试方程上的稳定区域为 $|R(z)|\le1$，也有边界。这与整体误差阶描述的是两件事：$h\to0$ 时在固定终点收敛，不保证固定 $h$ 后可以无限延长模拟。

刚性与隐式 Euler

教材 §11.4从多个时间尺度解释刚性：快速衰减的模态会限制显式方法的步长，而读者关心的慢变化要观察很久。例如

$$x'=-x,\qquad w'=-1000w,\qquad x(0)=w(0)=1.$$

两个衰减时间尺度分别为 1 和 $0.001$。快速分量已经很小时，慢分量仍在变化；显式 Euler 为了让快速模态衰减仍需 $h<0.002$。舍入和离散化误差也可能重新激发这个模态。刚性取决于方程、方法和所需精度，不能只凭曲线看上去陡不陡判断。

隐式 Euler 把斜率取在未知的下一状态：

$$y_{n+1}=y_n+h f(t_{n+1},y_{n+1}).$$

对测试方程，移项得到

$$y_{n+1}=\frac{y_n}{1-h\lambda}.$$

由乘数直接得到稳定区域 $|1-z|\ge1$，包含整个左半平面，因此它是 A 稳定的。对实数 $\lambda<0$，任意 $h>0$ 都让该模态衰减。比如快速分量取 $h=0.1$ 时，乘数为 $1/101$。这允许跨过快速瞬态，但若需要看清瞬态，仍须缩小步长。隐式 Euler 的整体误差仍是一阶；稳定性不会自动提高精度。

教材 §6.7介绍了隐式方法的实现：通常要先求出方程的根，才能确定下一状态。对一般非线性方程，隐式 Euler 每步要解 $G(u)=0$，其中

$$G(u)=u-y_n-hf(t_{n+1},u).$$

Newton 迭代的一次修正满足

$$\bigl[I-hJ_f(t_{n+1},u^{(k)})\bigr]\Delta u=-G(u^{(k)}), \qquad u^{(k+1)}=u^{(k)}+\Delta u.$$

$J_f$ 是 $f$ 对状态的 Jacobian 矩阵，$I$ 是单位矩阵。计算或近似 Jacobian、求解线性方程组、反复迭代都要花时间。还需要检查迭代是否收敛，让求解误差小于时间离散化的误差预算；一次失败的隐式求解不能当作合格的一步。

振子的长期能量：显式、隐式与辛方法

取单位质量、单位角频率的简谐振子：

$$q'=p,\qquad p'=-q,\qquad (q(0),p(0))=(1,0).$$

精确解为 $q(t)=\cos t$、$p(t)=-\sin t$，能量 $E=(q^2+p^2)/2$ 恒为 $1/2$。显式 Euler 同时用旧的 $q,p$ 更新：

$$q_{n+1}=q_n+hp_n,\qquad p_{n+1}=p_n-hq_n.$$

将两个平方相加，交叉项抵消，得到

$$E_{n+1}=(1+h^2)E_n.$$

任意固定正步长都会不断增加能量。$h=0.1$、$T=100$ 共 1000 步时，$E=\tfrac12(1.01)^{1000}\approx10479.577819$。隐式 Euler 恰好相反，能量每步除以 $1+h^2$，同一终点约为 $0.000024$：轨迹有界，但凭空加入了阻尼。

Hairer 的《Geometric Numerical Integration》第二讲给出辛 Euler：对这种可分离的 Hamilton 系统，先更新动量，再用新动量更新位置。

$$p_{n+1}=p_n-hq_n,\qquad q_{n+1}=q_n+hp_{n+1}.$$

“辛”指保留 Hamilton 流的几何结构；在这个二维例子中，更新矩阵的行列式为 1，保留相空间面积。它仍是一阶方法，也不逐步保持原能量。对这个振子，代入更新式可以验证，它精确保留的是修正能量

$$H_h=\frac12(q^2+p^2-hqp).$$

固定 $0<h<2$ 时这个二次型正定。利用 $|qp|\le E$ 和初始 $H_h=1/2$，可得

$$\frac{1/2}{1+h/2}\le E_n\le\frac{1/2}{1-h/2}.$$

因此 $h=0.1$ 时能量始终在约 $0.476190$ 与 $0.526316$ 之间振荡；这里的界限来自这个线性振子，不能直接套到所有 Hamilton 系统。恒定的 $h$ 也是这个守恒式的条件，逐步改变 $h$ 后不能沿用它。能量有界仍可能积累相位误差。

同一讲义的定理 2给出的二阶 velocity Verlet 采用半步动量、整步位置、再半步动量。对 $q''=a(q)$，公式是

$$\begin{aligned} p_{n+1/2}&=p_n+\frac h2a(q_n),\\ q_{n+1}&=q_n+hp_{n+1/2},\\ p_{n+1}&=p_{n+1/2}+\frac h2a(q_{n+1}). \end{aligned}$$

它也是适用于可分离 Hamilton 系统的辛方法；轨道实验用它推进中心引力模型。下面的算例用 $a(q)=-q$，可以直接比较长期能量变化与终点状态误差。

一个可运行的比较

将下列代码保存为 ode_compare.py，运行 python3 ode_compare.py。程序只用 Python 标准库。Euler 与 RK4 函数支持任意长度的状态序列；f(t, y) 必须返回相同数量的导数分量，否则抛出 ValueError。两个隐式函数是针对这两道线性例题直接解出的更新式，辛 Euler 和 Verlet 函数则专用于这个振子。

from math import cos, exp, hypot, sin


def slope(f, t, y):
    values = tuple(f(t, y))
    if len(values) != len(y):
        raise ValueError("f(t, y) must have the same length as y")
    return values


def euler(f, t, y, h):
    return tuple(a + h * b for a, b in zip(y, slope(f, t, y)))


def rk4(f, t, y, h):
    def shift(k, scale):
        return tuple(a + scale * b for a, b in zip(y, k))
    k1 = slope(f, t, y)
    k2 = slope(f, t + h / 2, shift(k1, h / 2))
    k3 = slope(f, t + h / 2, shift(k2, h / 2))
    k4 = slope(f, t + h, shift(k3, h))
    return tuple(a + h * (b + 2 * c + 2 * d + e) / 6
                 for a, b, c, d, e in zip(y, k1, k2, k3, k4))


def decay(t, y):
    return (-2 * y[0],)


def decay_implicit(f, t, y, h):
    return (y[0] / (1 + 2 * h),)


def oscillator(t, y):
    q, p = y
    return (p, -q)


def oscillator_implicit(f, t, y, h):
    q, p = y
    return ((q + h * p) / (1 + h * h),
            (p - h * q) / (1 + h * h))


def symplectic_euler(f, t, y, h):
    q, p = y
    p = p - h * q
    return (q + h * p, p)


def verlet(f, t, y, h):
    q, p = y
    half_p = p - h * q / 2
    q = q + h * half_p
    return (q, half_p - h * q / 2)


print("local Euler: h, one-step error")
for h in (0.1, 0.05):
    print(f"{h:.3f} {abs(1 - 2 * h - exp(-2 * h)):.9f}")

print("decay at T=1: method, h, y, absolute error")
for name, step in (("Euler", euler), ("implicit", decay_implicit), ("RK4", rk4)):
    errors = []
    for n in (10, 20, 40):
        h, y = 1 / n, (1.0,)
        for j in range(n):
            y = step(decay, j * h, y, h)
        error = abs(y[0] - exp(-2))
        errors.append(error)
        print(f"{name:8s} {h:.3f} {y[0]:.9f} {error:.3e}")
    print("error ratios: " + " ".join(f"{a / b:.3f}"
          for a, b in zip(errors, errors[1:])))

print("Euler with h=1.1: y0 through y4")
y, values = (1.0,), [1.0]
for j in range(4):
    y = euler(decay, j * 1.1, y, 1.1)
    values.append(y[0])
print(" ".join(f"{value:.4f}" for value in values))

print("oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error")
for name, step in (("Euler", euler), ("implicit", oscillator_implicit),
                   ("RK4", rk4), ("symplectic", symplectic_euler), ("Verlet", verlet)):
    y, energies = (1.0, 0.0), [0.5]
    for j in range(1000):
        y = step(oscillator, j * 0.1, y, 0.1)
        energies.append((y[0] ** 2 + y[1] ** 2) / 2)
    error = hypot(y[0] - cos(100), y[1] + sin(100))
    print(f"{name:10s} {min(energies):.6f} {max(energies):.6f} "
          f"{energies[-1]:.6f} {error:.3e}")

输出：

local Euler: h, one-step error
0.100 0.018730753
0.050 0.004837418
decay at T=1: method, h, y, absolute error
Euler    0.100 0.107374182 2.796e-02
Euler    0.050 0.121576655 1.376e-02
Euler    0.025 0.128512157 6.823e-03
error ratios: 2.032 2.016
implicit 0.100 0.161505583 2.617e-02
implicit 0.050 0.148643628 1.331e-02
implicit 0.025 0.142045682 6.710e-03
error ratios: 1.966 1.983
RK4      0.100 0.135339548 4.265e-06
RK4      0.050 0.135335528 2.452e-07
RK4      0.025 0.135335298 1.470e-08
error ratios: 17.396 16.682
Euler with h=1.1: y0 through y4
1.0000 -1.2000 1.4400 -1.7280 2.0736
oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error
Euler      0.500000 10479.577819 10479.577819 1.438e+02
implicit   0.000024 0.500000 0.000024 9.935e-01
RK4        0.499993 0.500000 0.499993 8.332e-05
symplectic 0.476190 0.526316 0.521321 5.665e-02
Verlet     0.498750 0.500000 0.499724 4.222e-02

能量最小值和最大值包括初值及全部网格点；state_error 是终点与 $(\cos100,-\sin100)$ 之差的欧氏范数。这里 RK4 的有限时间误差很小，辛方法的能量变化有界，但状态仍有偏差。选择方法时，应分别看终点精度、相位和长期守恒量，不能仅按阶数排序。

自适应步长与解的检查

教材 §6.5介绍嵌入式误差估计：同一步计算两个不同阶数的近似，用差值估计较低阶公式的局部误差。若估计误差超过目标，就拒绝这一步并用更小的 $h$ 重试；足够小时才接受。若估计量 $\eta$ 随步长约按 $h^{r+1}$ 变化，可用

$$h_{\mathrm{new}}=s h\left(\frac{\mathrm{tol}}{\eta}\right)^{1/(r+1)},\qquad 0<s<1,$$

作为调节依据，另限制步长增长和缩小的幅度。$r$ 是被估计公式的阶数，$s$ 为安全因子。指数要匹配误差估计量的阶数；这里给的是局部误差控制思路，不是所有求解器通用的更新公式。

SciPy 的 solve_ivp 文档说明：RK45 用嵌入式 5(4) 阶公式，以四阶公式估计误差、用五阶公式推进。它与上面的经典 RK4 不同。局部误差的尺度为 atol + rtol * abs(y)；接近零时需要合适的 atol，不同尺度的分量可分别设置它。非刚性问题可选 RK45，刚性问题可选 Radau 或 BDF。这些是局部误差控制设置，不能直接当成终点误差的保证。

检查一个解时，把比较目标固定下来：

1. 同一终点，减半步长。 对阶数为 $r$ 的方法，在渐近范围内计算 $y_h,y_{h/2},y_{h/4}$。差值比 $\|y_h-y_{h/2}\|/\|y_{h/2}-y_{h/4}\|$ 应接近 $2^r$；$y_{h/2}$ 的误差可估为 $\|y_{h/2}-y_h\|/(2^r-1)$；最细网格 $y_{h/4}$ 的误差则用 $\|y_{h/4}-y_{h/2}\|/(2^r-1)$ 估计。这些估计需要稳定、足够光滑、主导误差项不为零，且舍入尚未主导。比较整条轨迹时，使用共同的时间点。
2. 检查模型应有的守恒量。 对这个振子检查能量，也检查相位或状态误差；对中心引力模型可检查能量与角动量。守恒量接近不代表轨迹精确，耗散模型也不应强求能量守恒。
3. 对照独立参考。 有精确解时先用它；没有时，选适合方程的参考求解器，收紧容差，并再次收紧来确认参考值稳定。用 solve_ivp 对照固定步长程序时，t_eval 只指定保存结果的时间点，不指定内部步长。两次计算要使用相同方程、参数、初值和比较时间。

混沌系统还会放大微小状态差异。对双摆和 Rössler 吸引子，先在有限的短时间内检查步长收敛；长期再比较与模型相符的守恒量或统计特征。一条长期轨迹变得不同，本身无法区分初值敏感性与积分器误差。

我的阅读 · 打开过的笔记会记在当前浏览器里。