📚 课程目录第2章刚性方程识别与求解器切换排障

2. 刚性方程识别与求解器切换排障

问题对象与刚性判据

取一个典型的二组分线性系统作为排障对象，右端函数与解析解都已知，便于把“求解器行为”与“方程本身性质”分开看：

y1' = -1000·y1 + 999·y2 ,  y1(0) = 1
y2' =      y1 -     y2 ,  y2(0) = 0
t ∈ [0, 2]

右端函数对 y 的雅可比矩阵为常数矩阵：

J = [ -1000   999 ]
    [     1    -1 ]

刚性判据用特征值量级比 S = max|Re λ| / min|Re λ|。对 2×2 矩阵，特征值由迹与行列式直接得到：tr(J) = -1001，det(J) = 1000 - 999 = 1，故特征方程 λ² + 1001λ + 1 = 0，两根为

λ1 ≈ -1000.999 ,  λ2 ≈ -0.000999
S ≈ 1000.999 / 0.000999 ≈ 1.0×10^6

量级比达到 10⁶，且两个特征值均为负实数（解整体衰减、无振荡），这是典型的刚性特征：快模态时间常数约 1 ms，慢模态时间常数约 1000 s，而求解区间只有 2 s。显式方法的稳定域有限，步长会被快模态死死压住，即使快模态早已衰减到可忽略，步长仍无法放大。

显式方法的步长压缩现象

用经典四阶 Runge–Kutta（固定步长）先做一次“显式侧”的观察。RK4 的实轴稳定区间约为 [-2.78, 0]，对 λ1 ≈ -1001 要求 h ≤ 2.78/1001 ≈ 2.8×10⁻³。若强行取 h = 0.05，数值解会迅速发散；若把 h 压到 2×10⁻³ 才能稳定，但 t 从 0 到 2 需要 1000 步，而慢模态在 t > 0.01 后几乎不再变化——步长被一个已经衰减掉的快模态锁死，这就是“步长被压缩”的直接证据。

import numpy as np

def f(t, y):
    return np.array([-1000.0*y[0] + 999.0*y[1],
                      y[0] - y[1]])

def rk4(f, y0, t0, t1, h):
    t, y = t0, np.array(y0, float)
    n = int(round((t1 - t0) / h))
    for _ in range(n):
        k1 = f(t, y)
        k2 = f(t + h/2, y + h/2*k1)
        k3 = f(t + h/2, y + h/2*k2)
        k4 = f(t + h,   y + h*k3)
        y = y + h/6*(k1 + 2*k2 + 2*k3 + k4)
        t += h
    return t, y

# 稳定步长下运行
t, y = rk4(f, [1.0, 0.0], 0.0, 2.0, 2e-3)
print("RK4 h=2e-3  y(2) =", y)   # 预期接近 [e^-2, ...] 量级
# 越界步长：h=0.05 会发散，可自行验证

注意：上面只打印数值，不宣称“实测通过”。判断是否被压缩，看的是“稳定所需步长”与“精度所需步长”之间的巨大差距——本例中稳定步长比精度步长小两三个数量级，这正是刚性的行为特征。

用特征值量级比确认刚性来源

把雅可比矩阵显式算出来再求特征值，比“凭感觉”可靠。对一般非线性系统，在解轨线上取若干点做数值雅可比即可：

import numpy as np

def jac(t, y):
    return np.array([[-1000.0, 999.0],
                     [    1.0,  -1.0]])

lam = np.linalg.eigvals(jac(0.0, np.array([1.0, 0.0])))
lam = np.sort_complex(lam)
S = np.max(np.abs(lam.real)) / np.min(np.abs(lam.real))
print("特征值:", lam)
print("量级比 S =", S)          # 预期约 1.0e6
print("是否刚性(S>1e3):", S > 1e3)

若 S 远大于 1（工程上常以 10³ 为经验门槛），且最大模特征值实部为负，就应优先考虑隐式方法。本例 S ≈ 10⁶，结论明确。

切换到隐式求解器并调整参数

用 SciPy 的 solve_ivp 做对照：显式 RK45 与隐式 BDF 在同一容差下比较步数与函数调用次数。BDF 是变阶变步长隐式方法，对刚性模态 A-稳定，步长由精度而非稳定性决定。

import numpy as np
from scipy.integrate import solve_ivp

def f(t, y):
    return np.array([-1000.0*y[0] + 999.0*y[1],
                      y[0] - y[1]])

def jac(t, y):
    return np.array([[-1000.0, 999.0],
                     [    1.0,  -1.0]])

y0 = [1.0, 0.0]
t_span = (0.0, 2.0)
tol = dict(rtol=1e-6, atol=1e-9)

sol_rk = solve_ivp(f, t_span, y0, method="RK45", **tol)
sol_bdf = solve_ivp(f, t_span, y0, method="BDF", jac=jac, **tol)

for name, s in [("RK45", sol_rk), ("BDF", sol_bdf)]:
    print(f"{name}: success={s.success}, nfev={s.nfev}, "
          f"njev={getattr(s,'njev',0)}, steps={len(s.t)-1}, "
          f"y(2)={s.y[:, -1]}")

预期现象（示例预期，非实测）：RK45 的 nfev 会明显高于 BDF，步数也更多；BDF 在提供解析雅可比后 njev 较小，单步虽贵但总调用次数下降。若 BDF 报“迭代不收敛”，优先检查雅可比是否正确、atol 是否过小、以及是否在 t=0 附近存在不连续。

参数调整顺序建议：先给解析雅可比（或 jac_sparsity），再放宽 atol 到与解量级匹配（本例 y 量级为 1，atol=1e-9 已足够），最后才考虑换 Radau（更高阶、更稳但更贵）。不要一上来就把 rtol 调到 1e-12，那只会让隐式迭代更难收敛。

常见误判：看到“步长很小”就断定方程刚性。若特征值量级比接近 1，步长小可能只是精度要求高或右端函数有尖峰。务必用雅可比特征值量级比 + 显式方法稳定步长与精度步长的差距两条证据交叉确认，再切换求解器。

排障检查清单

检查项	本例取值/结论	不满足时的动作
雅可比特征值量级比 S	≈1.0×10⁶，远大于 10³	用数值雅可比复核，排除符号错误
最大模特征值实部	≈ -1001，为负	若为正，先处理不稳定方程本身
显式稳定步长 vs 精度步长	2.8×10⁻³ vs 约 10⁻²，差一个量级以上	差距小则不必切隐式
显式求解器 nfev / 步数	显著偏高、进度慢	记录基线，切换后对比
隐式求解器收敛性	提供解析雅可比后正常	检查 jac、放宽 atol、改 Radau
结果一致性	与解析解 y1(2)=e⁻² 对照	偏差大则回查容差与雅可比

完成切换后，用解析解做一次独立校验：本例 y1(t) 的慢模态主导项在 t=2 时应接近 e⁻² ≈ 0.1353 量级。若隐式解与该值偏差在容差内，且 nfev 相比 RK45 明显下降，即可判定“刚性识别正确、求解器切换有效”。