3. 刚性数学特征剖析

要判断一个微分方程“刚不刚”，不能只看它是否难解，也不能凭经验说“这个模型跑起来很慢”。刚性是一个可以量化的数学性质，它由系统自身的动力学结构决定，与所用求解器无关。本节围绕三个核心指标展开：刚性比、特征值分布和雅可比矩阵条件数，并说明如何把它们组合成对刚性程度的可操作描述。

一、刚性比：最直接的量化尺度

考虑线性常系数系统 y' = Ay，其中 A 的特征值为 λ₁, λ₂, …, λₙ。若所有特征值实部为负（系统稳定），定义刚性比为

R = max|Re(λᵢ)| / min|Re(λᵢ)|

直观含义是：系统中最快的衰减模态与最慢的衰减模态在时间尺度上相差多少倍。R 接近 1，说明各模态时间尺度相当，属于非刚性问题；R 达到 10³ 以上，通常已具备明显刚性；R 超过 10⁶ 时，显式方法基本不可用。

举一个教学示例（数值为合理假设）：某化学反应动力学线性化后得到特征值 λ₁ = −1.0、λ₂ = −2.0×10⁴。则 R = 2×10⁴，说明最快模态的时间常数约为 5×10⁻⁵ 秒，而最慢模态约为 1 秒。若用显式 Euler 法，为保持稳定，步长需满足 h < 2/|λ₂| ≈ 10⁻⁴ 秒；而要观察到 1 秒量级的慢过程，需要约 10⁴ 步以上。这就是刚性带来的“步长被最快模态绑架”的典型困境。

需要注意，刚性比只对稳定系统（特征值实部为负）有意义。若存在正实部特征值，系统本身不稳定，讨论刚性比会失去物理意义，此时应先处理稳定性问题。

二、特征值分布：比单一比值更丰富的信息

刚性比是一个标量，它压缩了全部特征值信息。实际工程中更值得看的是特征值在复平面上的分布形态。把全部特征值画在复平面上，观察它们相对虚轴的位置、实部跨度以及是否远离实轴，可以判断：

• 实部跨度：决定刚性比，跨度越大越刚。
• 虚部大小：虚部显著的特征值对应振荡模态，会限制显式方法的稳定域，即使刚性比不高也可能需要小步长。
• 聚集程度：若大量特征值聚集在左半平面远端，而少数靠近虚轴，说明系统存在“快变量簇”和“慢变量簇”，适合用隐式方法或指数积分方法处理。

对于非线性系统 y' = f(t, y)，特征值来自雅可比矩阵 J = ∂f/∂y 在解轨迹上的取值。由于 J 随状态变化，特征值分布也会随时间和状态漂移。因此工程上常在若干代表性时刻（如初始时刻、剧烈变化段、接近稳态段）分别计算 J 的特征值，观察刚性比是否在某个阶段突然升高。

三、雅可比矩阵条件数：衡量敏感性与数值难度

雅可比矩阵的条件数定义为

κ(J) = ‖J‖ · ‖J⁻¹‖

它刻画了 J 对输入扰动的放大能力。条件数越大，说明矩阵接近奇异，数值求解线性子问题（隐式方法每步都要解线性方程组）时误差越容易被放大。对于刚性系统，J 往往同时具有很大的范数和很小的最小奇异值，导致 κ(J) 很大。因此条件数可以作为隐式求解器数值稳定性与线性求解成本的预警指标：条件数过大时，即使方法本身是 A-稳定的，线性方程组的求解精度也可能下降，需要改用缩放、预处理或更高精度的线性代数。

四、把指标组合成刚性程度的描述

实践中可以建立一个简单的分级判断（示例阈值，需结合具体问题调整）：

下面这张图把三个指标与“快/慢模态分离”的物理图像对应起来，帮助建立从数学量到求解策略的直觉。

五、实际计算流程与注意事项

在工程中获取这些指标，通常按以下顺序操作。以 Python 为例，假设已用有限差分或自动微分得到某时刻的雅可比矩阵 J：

import numpy as np

eigvals = np.linalg.eigvals(J)
re = np.real(eigvals)
# 仅对稳定系统计算刚性比
if np.all(re < 0):
    R = np.max(np.abs(re)) / np.min(np.abs(re))
else:
    R = float('inf')  # 存在非稳定模态，先处理稳定性

kappa = np.linalg.cond(J)
print("刚性比 R =", R, " 条件数 κ =", kappa)

几个容易踩坑的地方：第一，不要用单个时刻的 J 下结论，应在解轨迹的多个阶段采样，取刚性比的最大值作为保守估计；第二，特征值实部接近零时，min|Re λ| 极小，R 会被放大到失真，此时应结合特征值分布图判断是否真的存在慢模态；第三，条件数依赖矩阵范数选择，2-范数与 ∞-范数给出的 κ 可能相差数量级，报告中应注明所用范数；第四，对于大规模稀疏系统，完整特征值分解代价高，可用幂迭代估计最大、最小特征值模，或用 scipy.sparse.linalg.eigs 只求部分特征值。

把上述指标落实到选型上：R 在 10² 以下且条件数温和，可优先考虑显式方法；R 进入 10³~10⁴ 且特征值分布呈现明显快慢分离，应转向隐式方法（如 BDF、Radau IIA）或隐式-显式组合方法；若 κ(J) 同时很大，则还需在隐式框架内加入线性求解预处理，否则每步的线性方程组求解会成为新的精度瓶颈。这样，三个数学指标就共同构成了从“识别刚性”到“选择算法”的量化桥梁。