4.3 刚性问题与隐式方法

4.3 刚性问题与隐式方法

本节摘要：刚性系统的特征是快慢时间尺度共存，显式方法的步长被绝对稳定性（而非精度）压制到快尺度量级，造成天文数字的浪费。后向欧拉与 BDF 类隐式方法以每步解一个非线性方程的代价换取无条件稳定。本节用 Van der Pol 方程完整复现刚性困境与隐式解法。

一、案发现场：一个"必须小步慢走"的方程

经典刚性试验品是刚性参数拉满的 Van der Pol 方程：

它的解在长达数百个时间单位里沿慢流形平缓滑行，只在极短瞬间快速切换。显式方法在这道题上的遭遇：

from scipy.integrate import solve_ivp
import numpy as np

def vdp_mu(mu):
    def rhs(t, y):
        return [y[1], mu*(1 - y[0]**2)*y[1] - y[0]]
    return rhs

t_end = 3000
sol_rk = solve_ivp(vdp_mu(1000), [0, t_end], [2.0, 0.0],
                   method='RK45', rtol=1e-6)
print(f"RK45 步数: {sol_rk.t.size-1}, 成功: {sol_rk.success}")
​

典型结果：几十万步，且很可能超时失败。怪异之处在于：解曲线平缓得像直线，从精度看大步长绰绰有余，RK45 却把步长压到 1e-3 以下——精度说可以走，稳定性说必须爬。

二、病理：绝对稳定域

对模型问题 y' = \lambda y（λ 为大负数，代表衰减极快的分量），显式欧拉的传播因子是 1 + h\lambda。它必须满足 |1 + h\lambda| < 1，即 h < 2/|\lambda|。λ = -1e4 时步长必须小于 2e-4——尽管这个快分量几个时间单位后就衰减到不可见，显式方法仍要在全部时间轴上为它站岗。

隐式（后向）欧拉把公式改成 y_{k+1} = y_k + h f(t_{k+1}, y_{k+1})，未知数出现在右边，每步要解方程（牛顿法，第三章工具再次上场）。它的传播因子是 1/(1 - h\lambda)，对任意大步长都小于 1——无条件稳定。代价与收益一目了然：每步从"一次函数评估"变成"若干次函数与雅可比评估加方程求解"，换来的是步长完全由精度决定。

图 4.3-1 显式与隐式欧拉的稳定域对比（示意）

三、结案演示：BDF 上场

向后微分公式（BDF）是隐式多步家族，SciPy 的 Radau、BDF 都是刚性问题的答案：

from scipy.integrate import solve_ivp
import numpy as np

def vdp_mu(mu):
    def rhs(t, y):
        return [y[1], mu*(1 - y[0]**2)*y[1] - y[0]]
    return rhs

t_end = 3000
sol_bdf = solve_ivp(vdp_mu(1000), [0, t_end], [2.0, 0.0],
                    method='BDF', rtol=1e-6)
print(f"BDF 步数: {sol_bdf.t.size-1}, 成功: {sol_rk_success if False else sol_bdf.success}")
# 步数通常是显式 RK45 的几十分之一
​

提供一个刚性的自测代码（化学动力学之外的另一个常见来源是带大系数扩散项的热方程）：

from scipy.integrate import solve_ivp
import numpy as np

# 快慢双尺度线性系统：快分量 e^(-1000t)，慢分量 e^(-t)
def rhs(t, y):
    return [-1000*y[0] - y[1], -y[1]]

for method in ['RK45', 'BDF']:
    sol = solve_ivp(rhs, [0, 10], [1.0, 1.0], method=method, rtol=1e-8)
    print(f"{method}: {sol.t.size-1} 步, success={sol.success}")
# RK45 步数是 BDF 的数十倍——这是刚性的标准指纹
​

四、诊断与处方速查

刚性没有严格单一定义，工作判据很实用：显式方法为了不发散被迫使用的步长，远小于精度所需的步长（相差百倍以上），即可判定刚性。常见来源与处方：

CFL 条件是同一病理在偏微分方程里的化身：显式格式的时间步长被限制在 h 乘以波速量级（对流问题）或 h 平方量级（扩散问题）。热传导方程显式格式的步长上限正比于网格间距的平方——网格加密一倍，允许步长缩到四分之一，计算量翻八倍。这条"隐形法规"解释了为什么工业热传导求解器清一色用隐式格式。

⚠️ 常见坑一：见 RK45 步数爆炸就收紧 rtol。南辕北辙——稳定性问题调精度参数无效，正确动作是换 method 为 BDF 或 Radau。
⚠️ 常见坑二：以为隐式方法"更高级所以更准"。隐式换的是稳定性不是精度阶：后向欧拉只有一阶精度，比 RK4 粗得多。它的赢面 solely 在步长自由。

本节要点回顾

• 刚性的工作定义：步长被稳定性而非精度压制，显式方法在为早已衰减的快分量无限站岗
• 稳定域决定论：显式欧拉的圆盘 vs 后向欧拉的整个左半平面；判据是 hλ 的位置
• 隐式的汇率：每步解一次非线性方程（通常是简化牛顿），换取步长完全交给精度决定
• BDF 家族：刚性 ODE 的工业标准，与 Radau 一起是 solve_ivp 的刚性答案
• CFL 是同案犯：PDE 显式格式的步长法规，扩散问题尤其严苛（h 平方量级）

第四章结卷。第五章离开确定性时间线，进入高维与随机的疆域——在那里，误差的主要形态从"步进累积"变成"方差"。