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

科学和工程中，几乎所有有意思的微分方程都无法求得解析解。非线性向量场、变系数、成千上万个耦合的状态变量——纸笔在问题真正变得棘手之前就早已无能为力。数值积分是关键所在。本章将构建、评估并比较一小套算法，它们几乎能解决你遇到的任何常微分方程（ODE），同时还会提供诊断工具，帮你识别积分器何时在误导你。

为何我们需要数值方法？#

前几章介绍的解析方法——变量分离、积分因子、拉普拉斯变换、特征值展开——虽然强大，但极其脆弱。它们仅适用于特定类别的方程，一旦现实问题不再“符号友好”，这些方法便会立刻失效。以 $\frac{dy}{dx} = \sin(xy)$ 为例，它根本没有封闭形式的解。纳维-斯托克斯方程、三体问题、包含多个物种的化学反应网络、洛伦兹系统——所有这些都让符号方法束手无策。我们只能退而求其次，在一系列网格点 $x_n = x_0 + nh$ 上计算离散近似值 $y_n \approx y(x_n)$ 。

于是问题就变成了：

• 如何选择更新规则 $y_{n+1} = \Phi(y_n, h, f, \ldots)$ ？
• 误差如何随步长 $h$ 缩放？
• 当方程具有“刚性”——即小步长是为了稳定性而非精度时——方法表现如何？

本章将一一解答。

前向欧拉法：几何思想#

这就是前向欧拉法，也是所有显式单步法的原型。其直观理解是：每一步读取当前位置的斜率，沿该方向走一小段距离 $h$ ，然后重新读取斜率。就像在浓雾中爬山，始终只相信脚下的局部坡度。

精度阶#

将其与欧拉步 $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)$ 。这意味着步长减半，误差也大致减半——但计算成本却翻倍，回报率很低。

线性稳定性#

放大因子为 $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)

Heun 法、中点法与 Runge-Kutta 家族#

系数 $\{a_{ij}, b_i, c_i\}$ 通常排列成一个 Butcher 表。

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)$ 。步长减半，误差降至四分之一，已有显著改善。

经典 RK4 法，四阶#

权重 $(1, 2, 2, 1)/6$ 与辛普森积分公式一致——这并非巧合。其全局误差为 $\mathcal{O}(h^4)$ ，每步需四次函数求值，堪称科学计算中性价比极高的经典算法。对于右端光滑的非刚性问题，这一算法在过去一个多世纪里一直是工程领域的默认选择。

 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=0.1$ 时达到 RK4 的精度，你需要将步长缩小至 $h \approx 10^{-5}$ ——计算成本增加十万倍。

收敛阶的可视化呈现#

检测一个方法收敛阶最清晰的方式，是对一系列步长计算终点处的全局误差，并在双对数坐标系中绘图。此时曲线的斜率即为方法的阶数。

这里有两个实用经验：

• 阶数是一道硬性上限。无论怎么调参，$p$ 阶方法的收敛速度都不可能超过 $\mathcal{O}(h^p)$ 。若想更快，必须换用更高阶方法或采用外推技术。
• 存在真实的误差下限。在光滑问题上，RK4 在 $h = 10^{-3}$ 时已接近机器精度。继续减小 $h$ 只会导致舍入误差在更多步数中累积，反而使结果更差。

自适应步长控制#

真实解的光滑性往往不均匀：可能有平台区、尖锐瞬变和缓慢衰减的尾部。固定步长要么在平台区浪费算力，要么在瞬变区因步长过大而危险。解决方案是自适应步长：动态调整 $h$ ，使局部误差始终接近用户指定的容差。

其中 0.9 是安全因子，指数项源于 $p+1$ 阶的渐近特性。最流行的嵌入对是 Dormand-Prince RK4(5)——也就是 scipy.integrate.solve_ivp 中的 “RK45”。每步只需 6 次函数求值，即可同时得到 4 阶和 5 阶的解估计。

刚性问题：显式方法的失效模式#

某些方程在缓慢演化之上叠加了少量快速动态过程。一旦快过程衰减完毕，你自然希望采用大步长推进——但显式方法却不允许。它们要求 $h \lesssim 1/|\lambda_{\max}|$ 始终成立，仅仅是为了维持数值稳定性。这类问题被称为刚性问题。

当 $\mu = 1000$ 时，慢时间尺度约为 $\mu$ ，快时间尺度约为 $1/\mu$ ，两者相差 $10^6$ 倍。若试图用 RK45 积分此系统，结果将是灾难性的。

它在新点处计算斜率。代价是每步都需要求解关于 $y_{n+1}$ 的代数方程，通常借助牛顿法完成。但回报巨大：后向欧拉法的稳定区域覆盖整个左半复平面（即A-稳定）。无论问题多么刚性，你都可以自由选择满足精度需求的步长 $h$ 。

稳定区域图示#

对于特征值沿负实轴大幅延伸的刚性问题，显式方法的稳定区域只是一条狭窄通道。你要么将 $h$ 缩小到能容纳进去（计算成本极高），要么直接使用 A-稳定方法（完全不受 $h$ 限制）。

在实际应用中，主流的隐式方法家族包括：

• BDF（后向差分公式，1–5 阶）：工业级刚性 ODE 求解器（如 LSODE、CVODE、SciPy 的 'BDF'）的默认选择。
• Radau IIA（3、5、9 阶）：具有强稳定性的隐式 Runge-Kutta 方法；对应 SciPy 中的 'Radau'。
• 隐式中点法 / 梯形法则：A-稳定、二阶精度，且对哈密顿系统具有能量守恒特性。

实用建议：如果 solve_ivp(..., method='RK45') 运行极慢或无法满足容差要求，不妨尝试 'BDF' 或 'Radau'。若问题在刚性与非刚性状态间切换，可使用 'LSODA'，它能自动检测并切换方法。

其他值得了解的方法#

多步法：Adams-Bashforth 与 Adams-Moulton#

Adams-Bashforth（显式）和 Adams-Moulton（隐式）是经典的多步法家族。预测-校正对通常结合 AB 预测与 AM 校正。其优势在于启动后每步仅需一次函数求值；缺点是对不连续点敏感，且需要单独的启动过程。

哈密顿流的辛积分器#

对于能量守恒系统（如行星轨道、分子动力学、格点规范理论），常规积分器会导致能量漂移：守恒量因舍入误差累积而缓慢增长或衰减。辛积分器则能保持相空间的辛 2-形式结构；虽然也不能严格守恒能量，但能将能量误差限制在真实值附近的有界振荡范围内，即使积分百万个周期也不会发散。

该方法具有二阶精度、时间可逆性且辛，是 N 体天体物理的标准工具。

SciPy 实战指南#

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

from scipy.integrate import solve_ivp
import numpy as np

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)

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)

sol = solve_ivp(rhs, [0, T], y0, method='LSODA')

solve_ivp 最常用的选项包括：

• method：'RK45'（默认）、'RK23'、'DOP853'（8 阶，超高精度）、'Radau'、'BDF'、'LSODA'。
• rtol, atol：各分量的相对与绝对容差。默认值较宽松（$10^{-3}, 10^{-6}$ ）；严肃计算中应收紧至 $10^{-8}, 10^{-10}$ 或更低。
• dense_output=True：返回一个插值器 sol.sol(t)，可在任意时刻查询解，而不局限于积分器选择的步点。
• events：传入一个函数，求解器会通过根查找定位其零点——非常适合处理碰撞、反弹或阈值穿越事件。
• jac：为刚性方法提供解析雅可比矩阵，相比有限差分近似可带来显著加速。

实用可靠性检查清单#

1. 用不同容差各运行一次。若将 rtol 缩小 100 倍后结果明显变化，说明较宽松的那次结果不可靠。
2. 绘制步长历史图。若求解器频繁触及 min_step，说明存在问题；若函数求值次数 nfev 极高，很可能是在用显式方法处理刚性问题。
3. 务必检查 success 字段。solve_ivp 在失败时不会主动报错，必须手动确认。
4. 对长时间积分，使用稠密输出并在粗网格上复核。单步检查无法发现漂移，但在 $10^6$ 个周期后可能已非常明显。

方法选择速查表#

一句话建议：先用 solve_ivp(method='RK45')；若速度慢或不稳定，改用 'BDF' 再试；若问题涉及守恒量且你关心其长期行为，请考虑辛积分器。

总结#

• 欧拉方法背后的几何思想，以及为何经典 Runge-Kutta 方法几乎能免费帮你找回精度
• 收敛阶 $\mathcal{O}(h^p)$ 的清晰推导，及其在双对数图中的真实含义
• 通过嵌入式 RK 方法实现自适应步长控制（scipy.integrate.solve_ivp 的核心引擎）
• 刚性问题：它是什么、如何识别、为何显式方法会灾难性失效，以及隐式方法（BDF、Radau）如何力挽狂澜
• 复平面上的线性稳定区域，以及它们如何指导你选择步长 $h$
• 简要介绍多步法、适用于哈密顿流的辛积分器，并附上一份生产环境使用清单

前置知识：掌握 第 1–4 章 中的基本 ODE 概念，并熟悉泰勒级数展开。

练习题#

1. 阶数验证。实现 euler、heun 和 rk4。求解 $\dot y = -2y$ ，初值 $y(0)=1$ ，区间 $[0, 3]$ ，步长 $h \in \{0.5, 0.25, 0.125, 0.0625\}$ 。绘制全局误差的双对数图，并验证斜率分别为 1、2、4。

2. 稳定边界实验测定。对 $\dot y = \lambda y$ （$\lambda = -10$ ）应用欧拉法，找出欧拉法失稳的临界步长 $h_*$ ，并验证其与 $2/|\lambda|$ 一致。

3. $$ \dot y_1 = -0.04 y_1 + 10^4 y_2 y_3, $$

初值 $y(0) = (1, 0, 0)$ ，区间 $[0, 10^{11}]$ 。分别使用 'RK45' 和 'BDF' 求解，比较运行时间和步数。

4. 自适应 vs 固定步长。使用 solve_ivp 求解 $\dot y = -y + 100\,e^{-(t-1)^2/0.005}$ ，区间 $[0, 3]$ ，相对容差 rtol=1e-6。再用固定步长 RK4 实现相同精度，报告两种方法的步数。

5. 辛积分器 vs RK4 在开普勒轨道上的表现。用 leapfrog 和 RK4（相同步长）对圆形开普勒轨道积分 $10^4$ 个周期，分别绘制能量随时间的变化曲线。

下一步#

下一章换个问题：之前研究的都是初值问题（IVP），现在换成边值问题（BVP）。BVP 在工程中很常见——梁的弯曲、热传导稳态、量子力学中的本征值问题——但解的结构完全不同：可能没有解、唯一解、或无穷多解。下一章给出射击法、有限差分、配点法等求解技术。

参考文献#

• Hairer, E., Norsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed.). Springer.
• Hairer, E. & Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.). Springer. The two Hairer-Wanner volumes are the bible.
• 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 (2nd ed.). Springer. The reference for symplectic methods.
• SciPy documentation: scipy.integrate.solve_ivp and the source of its method classes.

上一章: Chapter 10: Bifurcation Theory

下一章: Chapter 12: Boundary Value Problems

本文是常微分方程系列（共 18 篇）的第 11 篇。