4. 刚性方程与隐式方法

在非线性振子数值模拟中，有一类问题会让显式方法“寸步难行”：方程本身并不复杂，但解中同时包含变化极快的瞬态分量和变化缓慢的稳态分量，两者的时间尺度相差几个数量级。这类方程称为刚性方程（stiff equation）。识别刚性、理解显式方法为何失效、掌握隐式欧拉与梯形法、隐式龙格-库塔格式，并会用 SciPy 的 solve_ivp 处理刚性问题，是数值模拟从“能跑”走向“跑得对、跑得快”的关键一步。

一、刚性现象如何识别

考虑一个带强阻尼的非线性振子示例（参数为教学假设）：

m·x'' + c·x' + k·x + β·x³ = F·sin(ωt)
取 m = 1, c = 1000, k = 1, β = 0.5, F = 0, ω = 1
化为一阶系统：y1' = y2
             y2' = -1000·y2 - y1 - 0.5·y1³

该系统的雅可比矩阵特征值约为 λ₁ ≈ -0.001、λ₂ ≈ -1000。刚性比 |λ_max|/|λ_min| ≈ 10⁶，属于典型刚性系统。物理上，y2 方向对应极快的速度衰减（时间常数约 1 ms），y1 方向对应缓慢的位移演化（时间常数约 1000 s）。

识别刚性的实用判据有三条：其一，雅可比矩阵特征值实部全为负且模长跨度极大（刚性比远大于 1，工程上常以 10³ 以上为界）；其二，用显式方法（如 RK45）求解时，步长被稳定性而非精度“卡死”，步长被迫缩到 10⁻³ 量级才能不发散；其三，解在极短瞬态后迅速进入慢变阶段，但显式方法仍被迫用小步长爬行。对上面的例子，显式 RK45 若取 h = 0.01，稳定域要求 h·|λ| ≲ 2.8，即 h ≲ 2.8×10⁻³，否则数值解会剧烈振荡甚至溢出。

二、隐式欧拉与梯形法

隐式方法的核心是把未知的 y_{n+1} 放进右端函数中，用隐式方程求解。隐式欧拉格式为：

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

对线性测试方程 y' = λy，放大因子 R(z) = 1/(1 - z)，z = hλ。当 Re(λ) < 0 时 |R| < 1 对任意 h 成立，即A-稳定，步长不再受稳定性限制。梯形法（Crank–Nicolson）为：

y_{n+1} = y_n + (h/2)·[f(t_n, y_n) + f(t_{n+1}, y_{n+1})]

其 R(z) = (1 + z/2)/(1 - z/2)，同样 A-稳定，且是二阶精度，比隐式欧拉（一阶）更精确。但梯形法对极刚性分量（z → -∞）有 R → -1，衰减不足，可能出现高频振荡；隐式欧拉则 R → 0，具有强阻尼（L-稳定）。因此对强刚性瞬态，隐式欧拉更稳健；对精度要求高的慢变阶段，梯形法更合适。实际求解器常采用 BDF 类多步法，兼顾稳定性与精度。

隐式格式每步需解非线性方程，通常用牛顿迭代：以隐式欧拉为例，令 G(y) = y - y_n - h·f(t_{n+1}, y) = 0，迭代 y^{(k+1)} = y^{(k)} - J_G⁻¹·G(y^{(k)})，其中 J_G = I - h·∂f/∂y。每步需计算雅可比并解线性方程组，单步成本高于显式方法，但步长可放大几个数量级，总成本反而更低。

三、隐式龙格-库塔格式

隐式 RK 的一般形式为 s 级：k_i = f(t_n + c_i·h, y_n + h·Σ_j a_ij·k_j)，y_{n+1} = y_n + h·Σ_i b_i·k_i。当 a_ij 对 j ≥ i 非零时为隐式。常用的是Radau IIA（3 级 5 阶，L-稳定）和SDIRK（对角隐式，每级只需解一个非线性方程，适合大规模问题）。Radau IIA 的稳定性函数在 z → -∞ 时 R → 0，兼具高阶与强阻尼，是 SciPy 中 Radau 方法的理论基础。

四、用 SciPy 的 solve_ivp 处理刚性问题

SciPy 提供 solve_ivp，通过 method 参数选择求解器。刚性问题的推荐方法是 'Radau'（5 阶，L-稳定）和 'BDF'（变阶，适合大规模）。示例代码如下（参数为教学假设）：

import numpy as np
from scipy.integrate import solve_ivp

def stiff_osc(t, y):
    y1, y2 = y
    return [y2, -1000*y2 - y1 - 0.5*y1**3]

sol = solve_ivp(stiff_osc, [0, 20], [1.0, 0.0],
                method='Radau', rtol=1e-6, atol=1e-9,
                dense_output=True)
print(sol.success, sol.nfev, sol.t.size)
# 对比显式方法
sol_rk = solve_ivp(stiff_osc, [0, 20], [1.0, 0.0],
                   method='RK45', rtol=1e-6, atol=1e-9)
print(sol_rk.nfev, sol_rk.t.size)

运行后通常可观察到：Radau 的步数在数百量级，而 RK45 因稳定性限制步数激增（可能上万步），甚至因步长过小触发 status=-1 失败。判断依据是 sol.success 为 True、sol.nfev（函数调用次数）合理、解曲线在瞬态后平滑进入慢变。若 Radau 仍报“步长过小”，应检查雅可比是否病态、初值是否落在极陡区域，可提供解析雅可比 jac=... 加速牛顿迭代。

五、实践要点与结果判读

使用 solve_ivp 处理刚性问题时，建议遵循以下做法。第一，先用 RK45 试跑并观察 nfev 与步数，若异常大或失败，即判定为刚性，切换 Radau 或 BDF。第二，设置合理的 rtol（如 1e-6）与 atol（如 1e-9），刚性问题的分量量级差异大，绝对容差过松会掩盖慢变细节。第三，尽量提供解析雅可比 jac，可显著减少牛顿迭代次数。第四，用 dense_output=True 获得连续解，便于在任意时刻插值观察瞬态到稳态的过渡。第五，结果判读看三点：success 是否为 True、nfev 是否在合理量级、解曲线是否在短暂瞬态后平滑演化且无高频伪振荡。若出现伪振荡，说明所用方法阻尼不足，应改用 L-稳定的 Radau 或隐式欧拉类格式。