刚性微分方程

刚性微分方程（stiff differential equation）是数值求解中这样一类初值问题：解在所关注区间可能变化得很慢，但某些快速衰减的扰动迫使显式方法采用远小于精度所需的步长。刚性描述的是方程、解的区间、容限和数值方法之间的关系，并没有一个适用于所有问题的单一阈值。

典型情形是一个化学反应或电路中同时存在很快的松弛过程与很慢的外部变化。快速过程已经衰减后，若计算器仍被它限制着步长，就可能把大量工作用在并不需要细看的一段平稳轨迹上。

慢曲线旁边的快扰动

考虑无量纲模型 $$y^{\prime }=-\mu (y-\cos t)-\sin t,y(0)=1,\mu =1000.$$ 代入即可核验 $y(t)=\cos t$。这条曲线在一单位时间内平稳变化，单看它的图像并不显得困难。现在把初值改为 $1+\epsilon$。令 $w=y-\cos t$，方程恰好化为 $$w^{\prime }=-\mu w,w(t)=\epsilon e^{-\mu t}.$$ 因此真实轨迹由一个慢部分与一个快部分组成： $$y(t)=\cos t+\epsilon e^{-1000t}.$$ 扰动的衰减时间尺度约为 $10^{-3}$，而余弦的变化时间尺度约为一。即使选 $\epsilon =0$，离散误差和舍入误差仍可能激发快分量，所以精确初值并不会消除算法的稳定性问题。Driscoll 与 Braun 用稳定性和快慢尺度解释这种计算困难。

取初始偏差 0.1，真实解为 cos(t)+0.1exp(−1000t)。放大的早期窗口显示快速层；初值恰为一的解始终为 cos(t)。

同一个问题的较长时间窗口。快速层与后续缓慢变化所需的时间分辨率并不相同。

显式 Euler 的限制来自稳定

用显式 Euler 推进时， $$y_{n+1}=y_n+h[-\mu (y_n-\cos t_n)-\sin t_n].$$ 记相对慢曲线的误差为 $e_n=y_n-\cos t_n$，两边相减得到 $$e_{n+1}=(1-\mu h)e_n+\cos t_n-h\sin t_n-\cos (t_n+h).$$ 最后三项是局部离散缺陷；即使旧误差为零，它通常也不为零。误差的齐次放大因子为 $1-\mu h$。要使扰动严格衰减，必须 $$|1-\mu h|<1,0<h<2/\mu =\mathrm{0.002.}$$ 若取看起来已经不大的 $h=0.01$，放大因子却是 $-9$：每次换号并放大九倍。曲线本身平稳，也救不了这个递推。提高浮点精度只能延迟某些误差的出现，不能把不稳定因子改为衰减因子。

步长限制也要结合观察区间判断。如果确实要分辨初始快层的形状，小步长可能是精度要求，而非浪费；在快层之后追踪慢过程时，稳定性限制才尤其显得苛刻。这是刚性与单纯“解变化很快”的区别。

后向 Euler 的衰减与误差界

把右端放在新时刻，可解出 $$y_{n+1}=\frac{y_n+\mu h\cos t_{n+1}-h\sin t_{n+1}}{1+\mu h}.$$ 其齐次误差因子为 $R=1/(1+\mu h)$，对任意 $h>0$ 都在零和一之间。更强的是，当 $\mu h$ 很大时它趋于零，能迅速去掉强衰减模式。

对本例还可以证明一个实际误差界。把精确的余弦代进更新式，剩余缺陷为 $$d_n=\cos t_n-\cos t_{n+1}-h\sin t_{n+1}.$$ 在新时刻向后作 Taylor 展开，由 $|(\cos t{)}^{″}|\le 1$，有 $|d_n|\le h^2/2$。所以 $$|e_{n+1}|\le R|e_n|+\frac{h^2}{2(1+\mu h)}.$$ 初始误差为零，等比求和得到 $$|e_n|\le \frac{h}{2\mu }(1-R^n)\le \frac{h}{2\mu }.$$ 例如 $\mu =1000,h=0.01$ 时，这个界是 $5\times 10^{-6}$。它利用了本例特殊的已知慢曲线，不能当作所有刚性问题的通用误差常数；隐式方法在更一般方程中仍须估计误差并调整步长。

取 μ=1000、h=0.01、精确初值。纵轴画绝对误差的对数；零误差的初始点不放到对数轴上。比较体现的是稳定性，而非两种方法的形式阶数不同。

A 稳定与强阻尼

在线性测试方程 $y^{\prime }=\lambda y$ 中，A 稳定要求整个左半平面都属于稳定区域。后向 Euler 的因子 $1/(1-z)$ 满足这一性质，并在负实轴远端趋零。与之相比，隐式梯形法同样 A 稳定，但因子 $$\frac{1+z/2}{1-z/2}⟶-1(z\to -\mathrm{\infty })$$ 会保留换号的快分量。取 $z=-100$，一步因子为 $-49/51$，这说明“稳定”与“迅速消掉本应衰减的分量”还不是一回事。

后向微分公式和隐式 Runge–Kutta 方法是常见刚性求解途径。真实软件还要控制误差、求解隐式非线性方程并处理 Jacobian。SUNDIALS 将非刚性 Adams 与刚性 BDF 的用途和求解环节分别说明；SciPy 的官方说明也按刚性与否区分显式方法和 Radau、BDF 等方法。

多维线性系统中，Jacobian 特征值的衰减尺度差异可以提供线索，但不能机械地把最大、最小特征值模之比当成普遍定义。复特征值、非正规矩阵的暂态放大、非线性及所选容限都可能影响实际困难。改变时间单位还会改变系数的数值，所以“系数很大”也不是独立的判断标准。

同一问题的可运行对照

以下代码专门比较本例的显式与后向 Euler，不实现通用刚性求解器。输入为正的 $\mu ,h$、非负整数步数及有限初值；返回时间与状态列表，不修改输入。两个方法每步均为常数成本，时间 $O(N)$，保存轨迹占 $O(N)$。对一般右端，隐式每步还需要额外的求根成本。

Python 3

import math


def relaxation(mu, h, steps, implicit, y0=1.0):
    if not isinstance(steps, int) or steps < 0:
        raise ValueError("steps must be a nonnegative integer")
    if not all(math.isfinite(v) for v in (mu, h, y0)) or min(mu, h) <= 0:
        raise ValueError("finite data and positive mu, h required")
    if not isinstance(implicit, bool):
        raise ValueError("implicit must be a boolean")
    z = mu * h
    if not math.isfinite(z):
        raise OverflowError("step coefficient is not representable")
    ts, ys = [0.0], [float(y0)]
    for n in range(steps):
        t, new_t, y = n * h, (n + 1) * h, ys[-1]
        if not math.isfinite(new_t):
            raise OverflowError("time outside floating-point range")
        if implicit:
            new = (y + z * math.cos(new_t) - h * math.sin(new_t)) / (1 + z)
        else:
            new = y - z * (y - math.cos(t)) - h * math.sin(t)
        if not math.isfinite(new):
            raise OverflowError("unstable growth or arithmetic overflow")
        ts.append(new_t)
        ys.append(new)
    return ts, ys


for implicit in (False, True):
    ts, ys = relaxation(1000, 0.01, 10, implicit)
    error = max(abs(y - math.cos(t)) for t, y in zip(ts, ys))
    print(implicit, error)

把初值改成 $1.1$ 时，核验用的精确解也必须改为 $\cos t+0.1e^{-1000t}$。如果仍与 $\cos t$ 比较，就把真实的初始快层误判成算法误差。选求解器之前，应先区分物理瞬态、离散误差和稳定性要求。