刚性问题与非刚性问题

刚性是stiff这个词的翻译，其实并不容易建立直观的物理直觉。但如果把非刚性解释成柔性，那么理解起来会比较方便，即系统随时间的变化比较柔和，即为非刚性问题。这样一来，刚性问题就是，在时间尺度上，突然出现某些迅速而剧烈的变化。

所谓“迅速而剧烈”，是一个相对性的词汇，一个更加规范的描述是，系统中同时存在极快和极慢的变化过程，这个特性导致我们在根据时间划分步长的时候，很难找到一个平衡点，如果迁就缓变区域，步长设得很长，将捕捉不到快速变化的过程，导致求解失败；反之，如果从快速变化的过程出发，使用极短的步长，那么缓变区域将浪费大量的时间。

在sciML中，对于非刚性问题，首推Tsit5求解器，对于刚性问题，则推荐Rodas5。

Tsit5是一种显示龙格库塔法，全称是Tsitouras 5/4 Runge-Kutta，主解阶数为5阶，4阶误差估计，其下一步状态$u_{n+1}$直接由当前状态和若干中间斜率组合得出

$$\begin{array}{rl}{k}_i& =f(t_n+c_{i}h,u_n+h\sum _{j=1}^{i-1}{a}_{ij}{k}_j)\\ u_{n+1}& =u_n+h\sum _{i+1}^s{b}_i{k}_i\end{array}$$

Rodas5算法则属于隐式Rosenbrock算法，同样5阶主解、4阶误差估计，但在计算过程中，每一步都需要求解线性方程组

$$\begin{array}{rl}(I-\gamma hJ)k_i& =f(t_n+c_{i}h,u_n+h\sum _{j-1}^{i-1}{a}_{ij}{k}_j)+hJ\sum _{j=1}^{i-1}{c}_{ij}{k}_j\\ u_{n+1}& =u_n+h\sum _{i=1}^s{b}_i{k}_i\end{array}$$

非刚性问题求解

Lotka-Volterra 模型是一个典型的非刚性问题，用于描述捕食者和猎物在封闭环境中的行为。比如草原上只有两种动物，狐狸和兔子，二者的种群密度分别是$y(t),x(t)$，则二者随时间演化如下所示

$$\begin{array}{r}\frac{\mathrm{d}x}{\mathrm{d}t}=\alpha x-\beta xy\\ \frac{\mathrm{d}y}{\mathrm{d}t}=\delta xy-\gamma y\end{array}$$

式中，$\alpha$为兔子$x$的自然增长率，$\gamma$为狐狸$y$的固有死亡率，$\beta$为捕食率，$\delta$为捕食转化效率。该模型中，所有参数量级一致，其雅可比矩阵特征值的量级也接近，所以是非刚性问题。

此问题可用Tsit5来解决，得到两个物种的变化关系为

using DifferentialEquations, Plots

# 定义方程
function lotka(du, u, p, t)
    α, β, δ, γ = p
    du[1] = α*u[1] - β*u[1]*u[2]  # 兔子
    du[2] = δ*u[1]*u[2] - γ*u[2]  # 狐狸
end

u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
p = (1.5, 1.0, 3.0, 1.0)
prob = ODEProblem(lotka, u0, tspan, p)

# ✅ 推荐：Tsit5() 速度快
sol1 = solve(prob, Tsit5())
plot(sol1, title="Tsit5()解决非刚性问题", label=["兔子" "狐狸"])

刚性问题求解

刚性问题在化学反应中十分常见，例如，反应物$A$在生成$C$的过程中，有一中间产物$B$，其中$A\to B$这个过程相对比较缓慢，但$B$生成之后，会马上将$C$转换为$A$，并且以更快的速度，将自己变成$C$，写成伪化学方程式如下

$$\begin{array}{rlr}& A\stackrel{k_1}{\to }B& （慢，一级反应）\\ & B+C\stackrel{k_2}{\to }A+C& （快，二级，催化循环）\\ & B+B\stackrel{k_3}{\to }C+B& （极快，二级，自催化消耗）\end{array}$$

这个模型叫Robertson模型。当$k_i$之间存在明显的量级差别时，自然就出现了步长到底迁就谁的问题，即出现了刚性，其数学表达式为

$$
\begin{aligned}
\dot y_1&=-k_1y_1+k_2y_2y_3\
\dot y_2&=k_1y_1-k_2y_2y_3-k_3y_2^2\
\dot y_3&=k_3y_2^2

\end{aligned}
$$

刚性问题的首选求解器是Rodas，求解结果如下

using DifferentialEquations, Plots

# 定义方程 (Robertson 问题)
function robertson(du, u, p, t)
    du[1] = -0.04*u[1] + 1e4*u[2]*u[3]
    du[2] =  0.04*u[1] - 1e4*u[2]*u[3] - 3e7*u[2]^2
    du[3] =  3e7*u[2]^2
    nothing
end

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 1e5)  # 注意：时间跨度很大
prob = ODEProblem(robertson, u0, tspan)

# sol1 = solve(prob, Tsit5())  求解结果不对

# ✅ 正确示范：Rodas5() 轻松求解
sol2 = solve(prob, Rodas5(), saveat=10.0) # saveat 避免输出太多点

plot(sol2, title="刚性问题：必须用 Rodas5()", label=["A" "B" "C"], xscale=:log10)
plot(sol1, title="刚性问题：必须用 Rodas5()", label=["A" "B" "C"], xscale=:log10)