常微分方程（十一）：数值方法

科学与工程中几乎所有有意思的微分方程都拒绝给出解析解：非线性向量场、变系数、上万个耦合状态变量——纸笔早在问题本身屈服之前就已经放弃。数值积分是穿过这道墙的方式。本章构建、评估、对比那一小套基本能解决你会遇到的所有 ODE 的算法，并给出判断积分器是否在欺骗你的诊断手段。

本章要点

• 欧拉方法背后的几何思想，及为什么经典 Runge-Kutta 几乎"白送"地把丢失的精度还给你
• 收敛阶$\mathcal{O}(h^p)$的清晰推导，以及它在 log-log 图上的真实含义
• 嵌入式 Runge-Kutta 实现的自适应步长（scipy.integrate.solve_ivp 的内核）
• 刚性问题：定义、识别、显式方法为何灾难性失败、隐式方法（BDF、Radau）如何拯救
• 复平面上的线性稳定区，以及它如何指导步长选择
• 多步法、Hamilton 流的辛积分器简介，外加生产环境的 checklist

前置知识：第 1-4 章 的 ODE 基础；Taylor 展开。

1. 为什么需要数值方法？

之前几章的解析方法——分离变量、积分因子、Laplace 变换、特征展开——很强大但脆弱：只对狭窄类别的方程有效，一旦真实问题脱离"对符号友好"的轨道立刻失效。例如$\frac{dy}{dx} = \sin(xy).$无封闭解。Navier-Stokes 方程、三体问题、任何含数种以上化学物种的反应网络、Lorenz 系统——全部击败符号方法。我们只能退而求其次：在网格点$x_n = x_0 + nh$上算离散近似$y_n \approx y(x_n)$。

接下来要回答三个问题：

• 更新规则$y_{n+1} = \Phi(y_n, h, f, \ldots)$该怎么选？
• 误差如何随$h$标度？
• 当方程"刚性"时（小$h$不是为了精度而是为了稳定）方法表现如何？

我们会一一回答。

2. 前向欧拉：几何思想

给定$\dot{y} = f(x, y)$、$y(x_0) = y_0$，最简单的更新就是用切线代替曲线：$y_{n+1} = y_n + h\,f(x_n, y_n).$这就是前向欧拉，也是所有显式单步格式的原型。直觉：每一步读出当前位置的斜率，沿那个方向走小距离$h$，再读斜率。像在大雾中走山坡，永远只信脚下的坡度。

精度阶

将真解 Taylor 展开：$y(x_n + h) = y(x_n) + h\,y'(x_n) + \tfrac{h^2}{2}y''(x_n) + \mathcal{O}(h^3).$与欧拉步$y_{n+1} = y_n + h f(x_n, y_n)$相减得局部截断误差$\mathcal{O}(h^2)$/步。在固定区间$[x_0, X]$上走$N=(X-x_0)/h$步，全局误差为$\mathcal{O}(h)$。步长减半，误差减半——双倍代价换微薄回报。

线性稳定性

把欧拉应用于测试方程$\dot{y} = \lambda y$（$\mathrm{Re}\,\lambda<0$）：$y_{n+1} = (1 + h\lambda)\,y_n.$放大因子$R(z) := 1 + z$（$z = h\lambda$）。稳定性要求$|R(z)| \le 1$。对实负$\lambda$有$h < 2/|\lambda|$；对纯虚$\lambda$（振荡问题），没有任何$h>0$稳定。欧拉方法至多是条件稳定，对纯振荡问题完全无用。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

import numpy as np

def euler(f, x0, y0, x_end, h):
    """前向欧拉。f(x, y) -> dy/dx。"""
    xs, ys = [x0], [y0]
    x, y = x0, y0
    while x < x_end - 1e-12:
        h_eff = min(h, x_end - x)
        y = y + h_eff * f(x, y)
        x = x + h_eff
        xs.append(x); ys.append(y)
    return np.array(xs), np.array(ys)

3. Heun、中点法、Runge-Kutta 家族

要修复欧拉的一阶精度，关键是每步多算几次$f$再聪明地组合。显式$s$级 Runge-Kutta 的通用结构：$k_i = f\!\left(x_n + c_i h,\; y_n + h\textstyle\sum_{jButcher 表。

Heun（改进欧拉），二阶

预测-校正两级方法：

1. 预测：$\tilde y = y_n + h\,f(x_n, y_n)$
2. 校正：$y_{n+1} = y_n + \tfrac{h}{2}\bigl[f(x_n, y_n) + f(x_n + h, \tilde y)\bigr]$

局部误差$\mathcal{O}(h^3)$，全局误差$\mathcal{O}(h^2)$。$h$减半，误差减为四分之一。已经好很多。

经典 RK4，四阶

四级工作马：

权重$(1,2,2,1)/6$与 Simpson 积分公式同构（不是巧合）。全局误差$\mathcal{O}(h^4)$，每步 4 次$f$求值——科学计算中最划算的交易之一。对光滑右端的非刚性问题，这个算法已经做了一个多世纪的工程默认值。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

def rk4(f, x0, y0, x_end, h):
    xs, ys = [x0], [y0]
    x, y = x0, y0
    while x < x_end - 1e-12:
        h_eff = min(h, x_end - x)
        k1 = f(x, y)
        k2 = f(x + h_eff/2, y + h_eff*k1/2)
        k3 = f(x + h_eff/2, y + h_eff*k2/2)
        k4 = f(x + h_eff, y + h_eff*k3)
        y = y + h_eff*(k1 + 2*k2 + 2*k3 + k4)/6
        x = x + h_eff
        xs.append(x); ys.append(y)
    return np.array(xs), np.array(ys)

并排对比：欧拉 vs RK4

这不是小差距。从一阶到四阶意味着固定$h$下的误差预算被平方两次。要让欧拉在$h=0.1$时达到 RK4 的精度，需要$h\approx 10^{-5}$——十万倍的代价。

4. 收敛阶可视化

最干净的检验方法阶的方式：在一系列步长下计算终点全局误差，画 log-log 图，斜率即为阶数。

两个实用教训：

• 阶数是硬上限：对$p$阶方法，无论怎么调都好不过$\mathcal{O}(h^p)$。要更快收敛只能换更高阶或外推。
• 舍入下界真实存在：光滑问题上 RK4 在$h=10^{-3}$已经摸到机器精度。再降$h$只是让步数变多、舍入累积加剧——反而会变差。

5. 自适应步长控制

真实解并非均匀光滑——有平台、有尖峰、有缓慢长尾。固定$h$要么在平台浪费，要么在尖峰危险。修复方式是自适应步长：在线选择$h$使局部误差贴近用户给定容差。

标准机制是嵌入式 Runge-Kutta：两个阶为$p$与$p+1$的方法共享中间级求值，差为局部误差$E_n$的估计。若$E_n>\text{tol}$则拒绝该步、缩小$h$重试；否则接受并选下一步$h_\text{new} = 0.9\,h\,\bigl(\text{tol}/E_n\bigr)^{1/(p+1)}.$0.9 是安全因子，指数来自$p+1$阶渐近。最流行的嵌入对是 Dormand-Prince RK4(5)——solve_ivp 内部的 “RK45”。每步 6 次$f$求值，同时给出 4 阶与 5 阶估计。

6. 刚性问题：显式方法的失败模式

某些方程在缓慢演化之上叠加少量快动力学。快部分衰减后你很想走大步——但显式方法不允许，它们要求$h \lesssim 1/|\lambda_{\max}|$永远成立，仅为了数值稳定。这类问题称为刚性。

教科书例子：高非线性的 van der Pol 振子$\ddot{y} - \mu(1-y^2)\dot{y} + y = 0.$$\mu=1000$时慢尺度$\sim \mu$，快尺度$\sim 1/\mu$，比值 $10^6$。用 RK45 积分这个就是灾难。

出路是隐式方法。最简单的是后向欧拉：$y_{n+1} = y_n + h\,f(x_{n+1}, y_{n+1}),$在新点求斜率。代价：每步要解关于$y_{n+1}$的代数方程，通常用 Newton 法。回报：后向欧拉的稳定区是整个左半平面（A-稳定）。无论问题多刚性，$h$都只受精度约束。

稳定区可视化

刚性问题的特征值沿负实轴拉得很远，显式区是窄走廊：要么把$h$缩到能装下（巨大代价），要么用 A-稳定方法（无$h$约束）。

工业级隐式家族：

• BDF（后向差分公式，1-5 阶）。LSODE、CVODE、SciPy 'BDF' 的默认刚性求解器。
• Radau IIA（3、5、9 阶）。强稳定性的隐式 RK；SciPy 'Radau'。
• 隐式中点 / 梯形规则。A-稳定，二阶，对 Hamilton 系统能量守恒。

实用准则：若 solve_ivp(..., method='RK45') 永远跑不完或拒绝容差，换 'BDF' 或 'Radau'。若刚/非刚区域交替，'LSODA' 自动检测。

7. 还应该知道的几种方法

多步法：Adams-Bashforth 和 Adams-Moulton

不在每步内部多算$f$，而是复用之前几步的值：$y_{n+1} = y_n + h \sum_{j=0}^{k-1} \beta_j f(x_{n-j}, y_{n-j}).$Adams-Bashforth（显式）与 Adams-Moulton（隐式）是教科书家族。预测-校正对将 AB 预测与 AM 校正组合。优点：启动后每步只需 1 次$f$求值；缺点：在不连续点附近脆弱，且需要单独的启动过程。

Hamilton 流的辛积分器

对能量守恒系统（行星轨道、分子动力学、格点规范理论），常规方法引起能量漂移：守恒量随舍入累积慢慢增减。辛积分器保持相空间的辛 2-形式；它们也无法严格守恒能量，但能把能量误差限制在真值附近的有界振荡内，即便积分上百万周期。

最简例子是$\ddot q = -\nabla V(q)$的 Stormer-Verlet（leapfrog）：$p_{n+1/2} = p_n - \tfrac{h}{2}\nabla V(q_n),\quad q_{n+1} = q_n + h\,p_{n+1/2},\quad p_{n+1} = p_{n+1/2} - \tfrac{h}{2}\nabla V(q_{n+1}).$二阶精度、时间可逆、辛。N 体天体物理的标准工具。

8. SciPy 实战

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

from scipy.integrate import solve_ivp
import numpy as np

# 通用非刚性问题：默认 RK45 即可。
sol = solve_ivp(lambda t, y: -2*y, [0, 5], [1.0],
                t_eval=np.linspace(0, 5, 100),
                rtol=1e-8, atol=1e-10)

# 刚性问题：换 BDF 或 Radau。
def vdp(t, y, mu=1000):
    return [y[1], mu*(1 - y[0]**2)*y[1] - y[0]]
sol = solve_ivp(vdp, [0, 3000], [2.0, 0.0],
                method='BDF',
                rtol=1e-6, atol=1e-9)

# 不知道刚不刚？让 LSODA 自己决定。
sol = solve_ivp(rhs, [0, T], y0, method='LSODA')

solve_ivp 最常用选项：

• method：'RK45'（默认）、'RK23'、'DOP853'（八阶高精度）、'Radau'、'BDF'、'LSODA'。
• rtol、atol：相对/绝对容差（按分量）。默认值偏松（$10^{-3}, 10^{-6}$）；正经计算应收紧到$10^{-8}, 10^{-10}$或更低。
• dense_output=True：返回插值器 sol.sol(t)，可在任意$t$查询，不限于积分器自选的步。
• events：传入函数，求解器自动定位其零点（碰撞、反弹、阈值穿越）。
• jac：为隐式方法提供解析 Jacobi 矩阵，相对差分 Jacobi 大幅提速。

实用可靠性 checklist

1. 用两套不同容差各跑一次。若 rtol 缩小 100 倍后答案明显改变，说明松的那次是错的。
2. 画步长历史。求解器死磕 min_step 必有问题。nfev（求值次数）爆表通常意味着用了显式方法去碰刚性问题。
3. 必须检查 success。solve_ivp 失败时默认沉默返回，必须自己 assert。
4. 长积分要用稠密输出 + 粗网格复核。漂移在单步检查中不可见，但在$10^6$周期积累后会很明显。

9. 方法选择速查表

一句话建议：先用 solve_ivp(method='RK45')；若慢或不稳定换 'BDF'；若关心守恒量再看辛方法。

习题

1. 阶数验证。实现 euler、heun、rk4。在$\dot y=-2y$, $y(0)=1$, $[0,3]$上取$h\in\{0.5, 0.25, 0.125, 0.0625\}$，画 log-log 误差并验证斜率为 1、2、4。
2. 稳定边界实测。把欧拉应用到$\dot y=\lambda y$（$\lambda=-10$），找出欧拉失稳的临界步长$h_*$并验证它与$2/|\lambda|$一致。
3. 刚性系统两种解法。用 'RK45' 与 'BDF' 各解一次 Robertson 化学动力学 $$\dot y_1 = -0.04 y_1 + 10^4 y_2 y_3,\quad \dot y_2 = 0.04 y_1 - 10^4 y_2 y_3 - 3\times 10^7 y_2^2,\quad \dot y_3 = 3\times 10^7 y_2^2,$$ 取$y(0)=(1,0,0)$, 区间$[0,10^{11}]$。比较运行时间与步数。
4. 自适应 vs 固定。用 solve_ivp 在 rtol=1e-6 下积分$\dot y=-y+100\,e^{-(t-1)^2/0.005}$, $[0,3]$。再用足以达到同等精度的固定步长 RK4 重做。汇报步数。
5. 辛 vs RK4 在 Kepler 轨道上。同步长用 leapfrog 与 RK4 各积分$10^4$周期的圆 Kepler 轨道，画出能量随时间的曲线。

参考文献

• Hairer, E., Norsett, S. P. & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems（第 2 版）. Springer.
• Hairer, E. & Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems（第 2 版）. Springer. 上述两卷为本领域圣经。
• Ascher, U. & Petzold, L. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM.
• Hairer, E., Lubich, C. & Wanner, G. (2006). Geometric Numerical Integration（第 2 版）. Springer. 辛方法权威参考。
• SciPy 文档：scipy.integrate.solve_ivp 与其方法类源码。

系列导航