1. 课程导入与刚性概念建立

在数值求解常微分方程（ODE）的实践中，很多人会遇到一种令人困惑的现象：明明选用了精度很高、稳定性理论也很成熟的显式方法，步长却被迫压到极小，计算量急剧膨胀，而解本身变化并不剧烈。这类"看起来平缓、算起来昂贵"的问题，往往就是刚性方程。本节从物理、化学、生物等真实系统出发，建立刚性的直观概念，明确其定义与特征，并给出识别与算法选型的整体流程框架。

一、从真实系统看"快慢共存"

刚性并非数学家的抽象游戏，它广泛存在于工程与自然系统中。考虑一个带小热容传感器的温度测量装置：传感器本身响应极快（时间常数约 10⁻³ s），而被测炉体升温很慢（时间常数约 10² s）。若用显式欧拉法求解该传热方程组，为保持数值稳定，步长必须小于最快模态时间常数的两倍左右，即 h ≲ 2×10⁻³ s；而为了刻画炉体升温到稳态的完整过程，总时长需要数百秒。于是步数达到十万量级，绝大部分计算都"浪费"在早已衰减的快速瞬态上。

化学动力学中同样典型。一个包含快速可逆反应与慢速分解反应的机理，其速率常数可能相差 4～6 个数量级。生物领域里，神经元的膜电位快速放电与离子浓度的慢速恢复、酶促反应中底物快速结合与产物慢速生成，都是快慢时间尺度并存的例子。这些系统的共同点是：解中包含若干衰减极快的瞬态分量，它们在很短时间后对解的贡献可忽略，但它们的存在却严格约束了显式方法的步长。

二、刚性的数学刻画

对线性常系数系统 y′ = Ay，若矩阵 A 的特征值实部均为负（系统渐近稳定），定义刚性比

S = max|Re λᵢ| / min|Re λᵢ|

当 S 远大于 1（工程上常以 S ≳ 10³ 作为经验门槛）时，称系统为刚性的。对非线性问题 y′ = f(t, y)，则用雅可比矩阵 J = ∂f/∂y 在解轨迹附近的特征值谱来局部判断。刚性比越大，快慢分离越显著，显式方法的步长限制越苛刻。

需要强调：刚性是问题本身的属性，与所用算法无关；但"是否表现为计算困难"取决于算法。刚性问题的解通常光滑、无剧烈振荡，这正是它容易被误判为"简单问题"的原因。

三、稳定性视角：为什么显式方法会"卡住"

显式方法的绝对稳定域是有界的。以经典四阶 Runge–Kutta 为例，其稳定域沿负实轴约延伸到 −2.78。若某特征值 λ 满足 |λh| 落在稳定域外，数值解就会出现虚假振荡甚至发散。因此步长被最快模态"绑架"：h ≤ 2.78/|λ_max|。而隐式方法（如后向欧拉、BDF、Radau IIA）的稳定域覆盖整个左半平面或无界区域，称为 A-稳定，允许步长由精度而非稳定性决定，从而可以用远大于快模态时间常数的步长推进。

四、识别与选型的整体流程

面对一个新问题，建议按"建模—诊断—选型—验证"的闭环推进，下图概括了各环节的输入、判断与去向。

流程中的关键判断是刚性比估计。对小型问题可直接数值求雅可比矩阵的特征值；对大型问题，可用幂法估计 |λ_max|、用反幂法估计 |λ_min|，或借助求解器自带的刚性检测（如 MATLAB 的 ode15s 在启动阶段会尝试显式步进并监测稳定性受限情况）。若无法可靠估计，可采用"试算对比"策略：分别用显式与隐式方法在相同容差下求解，比较函数评估次数与耗时，若隐式方法显著占优，则问题大概率刚性。

五、算法选型的初步原则

确认刚性后，选型应围绕三个维度权衡：稳定性需求（是否 A-稳定、是否 L-稳定以抑制快模态的虚假残留）、精度阶数（BDF 常用 1～5 阶，高阶需注意稳定性退化）、实现成本（隐式方法每步需解非线性方程组，依赖雅可比与线性求解器）。工程经验是：中小规模、要求高精度时优先 Radau IIA（3 或 5 阶，L-稳定）；大规模稀疏问题优先变阶变步长 BDF（如 CVODE、ode15s）；对雅可比难以解析提供的问题，可考虑无雅可比的 Rosenbrock 或对角隐式 RK 方法。