› Runge-Kutta法 稳定性模拟器 返回
数值分析

Runge-Kutta法 稳定性模拟器

常微分方程求解的数值积分法（欧拉法·RK2·RK4）的「绝对稳定性」可视化工具。改变步长 h 和特征值 λ，即可实时了解放大率 R(z) 和复平面的稳定域，以及数值解是否发散。

参数设置
数值积分法
确定稳定函数 R(z) 的阶数
测试方程的 λ
y'=λy 的特征值。λ<0 时真解衰减
步长 h
每步前进的时间幅度
计算时间 t_end
s
进行积分的总时间
计算结果
—
z = hλ
—
|R(z)| 放大率
—
稳定判定
—
数值解最终值
—
精确解最终值
—
最大误差
—
复平面的绝对稳定域 |R(z)| ≤ 1

横轴为 z 的实部，纵轴为虚部。填充区域为所选方法的稳定域。当前工作点 z=hλ（在实轴上）用绿色（稳定）或红色（不稳定）表示。

数值解 vs 精确解 y(t)
放大率 |R(z)| 与步长 h 的关系
理论·主要公式

$$y'=\lambda y,\qquad y_{n+1}=R(z)\,y_n,\quad z=h\lambda$$

测试方程 y'=λy（精确解 y(t)=y₀e^(λt)）。数值解在每步乘以稳定函数 R(z)，n 步后为 y_n = R(z)ⁿ·y₀。

$$R_{\text{Euler}}=1+z,\qquad R_{\text{RK4}}=1+z+\tfrac{z^{2}}{2}+\tfrac{z^{3}}{6}+\tfrac{z^{4}}{24}$$

各方法的稳定函数。RK2 为 $R=1+z+z^{2}/2$。阶数越高，越接近指数函数 eᶻ，稳定域越大。

$$|R(z)|\le 1 \;\Longleftrightarrow\; \text{绝对稳定}$$

绝对稳定条件。满足此条件的 z 集合称为「绝对稳定域」。对衰减系统（λ<0），选择步长使 z=hλ 在此域内。

Runge-Kutta法的稳定性基础

🙋
我学过把微分方程用计算机求解时，步长 h 越小越精确。那么干脆把 h 设得很小不就行了？
🎓
精度角度确实是这样。但「精度」和「稳定性」是两个不同的概念。稳定性是指，当 h 太大时，数值解会像失控一样发散，结果完全错误。比如真解明明在平稳衰减到零，但数值解却在振荡增长走向无穷——这就是不稳定。h 太小当然能修复它，但那样的话计算工作量会成倍增加。你需要知道「h 最多能取多大而不出问题」，才能有效利用计算资源。
🙋
那发散和不发散是什么决定的？我试着把左边的 λ 拖成非常负的数，突然就变成红的「不稳定」了。
🎓
很好的观察！关键量是 z = hλ 这一个无量纲数。在求解测试方程 y'=λy 时，一步后数值解总是 y_{n+1} = R(z)·y_n 的形式。这里 R(z) 叫「稳定函数」。走 n 步，数值解就变成 y_n = R(z)ⁿ。所以当 |R(z)| 大于1时，每步都在放大，最后指数爆炸。|R(z)| ≤ 1 时才能稳定衰减。把 λ 弄成大的负数，z=hλ 就跑出了稳定区域，所以你看到红色了。
🙋
哦，那上面那个复平面的图——欧拉法和RK4的稳定域形状不一样，这就是原因吗？
🎓
完全正确。前进欧拉法的稳定函数是 R(z)=1+z，稳定域只是一个小圆盘（圆心 -1、半径1），在实轴上只能到 z∈[-2,0]。RK4 是 R(z)=1+z+z²/2+z³/6+z⁴/24，稳定域是个更大的叶片形，实轴上能伸到约 z≈-2.78。这意味着面对同一个衰减率 λ 的问题，RK4 能用更大的步长 h。从欧拉法的 h < 2/|λ| 升级到 RK4 的 h < 2.78/|λ|，计算效率会显著提升。
🙋
那这样的话，我们不应该总是选 RK4、把 h 拉到稳定边界吗？这样岂不是最快？
🎓
听起来合理，但这是个陷阱。在稳定域的边缘附近，|R(z)| 非常接近1，虽然确实不会爆炸，但精度变得很糟。本来应该快速衰减的解，由于 |R(z)| ≈ 1，衰减速度被拖得很慢，留下很多数值噪声。稳定性说的是「不发散的最低条件」，精度是「答案有多准」，这两个是独立的。实务中，我们通常在稳定上限的一半甚至更低的地方取 h，确保既稳定又精确。特别是遇到刚性方程时，这个制约会变得要命。
🙋
刚性方程是什么？好像很可怕的样子。
🎓
刚性方程（stiff equation）是快速衰减和缓慢演变混在一起的系统。表面上看解很温和地变化，但真实的物理过程中混有快速响应和慢速漂移。比如化学反应，有的反应在毫秒内达到平衡，有的反应要进行几小时。虽然最后我们只关心慢的那部分，但要稳定求解快速那部分，h 就得被压到 2/|λ_max| 以下。结果呢，就算解看着变化不快，你也得用细密的时间网格，计算量爆炸。显式Runge-Kutta法（欧拉法、RK4）稳定域有限，对付刚性很费力。标准做法是改用后向欧拉法这样的隐式方法，它们的稳定域能覆盖整个左半平面，不再受 h 的束缚。这工具中如果把 λ 设到 -30，你就能亲身体验到显式方法有多难受。

常见问题

绝对稳定性是指，在数值求解测试方程 y'=λy（λ<0 时真解衰减）时，数值解不发散而持续衰减的性质。每步的放大率称为稳定函数 R(z)（z=hλ），n 步后数值解为 y_n = R(z)^n。当 |R(z)| ≤ 1 时，重复迭代也不会使解变大，称为绝对稳定。当 |R(z)| > 1 时，尽管真解在衰减，数值解却指数级发散。
在实轴上，前进欧拉法 R(z)=1+z 仅在 z∈[-2,0] 时稳定（复平面上为圆心 -1、半径1的单位圆）。RK4 R(z)=1+z+z²/2+z³/6+z⁴/24 在实轴上约可稳定到 z∈[-2.78,0]，稳定域呈更大的叶状。也就是说，对同一个衰减率 λ 的问题，RK4 可以取更大的步长 h。从 h < 2/|λ|（欧拉法）扩展到 h < 2.78/|λ|（RK4）。
从稳定性角度，应选择步长使得 z=hλ 落在所选方法的绝对稳定域内。对于衰减系统（λ<0），欧拉法要求 h < 2/|λ|，RK4 约为 h < 2.78/|λ|。但稳定仅保证不发散，精度另当别论。实务中应在稳定上限之下再取更小的 h，以满足所需精度（局部截断误差）。稳定性是「不发散的最低条件」，精度是「结果有多准确」，两者独立。
刚性方程是指衰减很快的分量（|λ| 很大的特征值）和缓慢变化的分量混在一起的系统。虽然解看起来很温和，但为了稳定求解快速衰减分量，需要将 h 压到 2/|λ_max| 以下，导致计算效率极低。显式Runge-Kutta法（如欧拉法、RK4）的稳定域有限，对刚性问题无力，通常改用后向欧拉法等 A-稳定的隐式方法（稳定域涵盖左半平面）。本工具中将 λ 设为 -30 等很负的值时，可以亲身体验显式方法的局限。

实际应用

结构振动·过渡响应分析：建筑物或机械的动态响应通过时间积分求解时，高振动模态的特征值 |λ| 更大，对欧拉法的稳定条件 h<2/|λ| 要求更严。细化网格会导致最高固有振动频率增大，所需时间步也必须相应缩小。在冲撞、爆炸等显式分析中，稳定性限制以「Courant条件」的形式出现，网格尺寸直接支配时间步长。

化学反应·燃烧模拟：多个反应速率常数相差数个量级的反应方程是典型的刚性方程。如果用步长去适应最快反应，为了跟踪最慢反应，需要的总步数会大得离谱。实务中不使用欧拉法或RK4这样的显式方法，而改用后向微分公式（BDF）等隐式求解器，消除稳定性对步长的制约。

电路·控制系统数值分析：RC电路、运放电路的过渡分析中，寄生电容的时常数小，产生大的 |λ|。SPICE类电路模拟器采用梯形法或后向欧拉法，正是为了应对这种刚性。控制系统的状态方程离散化时，采样周期 h 和特征值的乘积 z=hλ 必须在稳定域内，是系统设计的前提。

CAE求解器设置确认：显式动力学分析中结果突然发散，往往是因为时间步超过了稳定限界。理解本工具展示的稳定域知识，你就能诊断「细化网格后发散了」或「材料变硬后发散了」的现象：这是因为 λ 增大，使 z=hλ 超出了稳定区域。反过来，如果改用隐式方法，就不必为了稳定而人为缩小 h。

常见误解及注意事项

最大的误解是「高阶方法必然稳定」的想法。RK4 确实比欧拉法精度高，在实轴方向稳定域更宽（-2 对比 -2.78），但「高阶=无条件稳定」是错的。显式Runge-Kutta法不管多高阶，稳定域都是有限的，永远无法覆盖左半平面全体。面对刚性问题，升高阶数没有本质帮助，必须切换到A-稳定的隐式方法。很多人吃过「用RK4还是爆炸了」的亏，根子就在这里。

另一个误解是「稳定就是正确」。|R(z)| ≤ 1 只保证解不爆炸，和精度无关。当 z 在稳定域边缘时，|R(z)| 接近1，本应迅速衰减的分量被人为拖延，甚至产生虚假振荡。稳定通过检验不代表数值解准确。标准做法是：稳定性通过后，还要进行网格/步长独立性检验，逐步缩小 h 看答案有没有变化，直到收敛为止。

最后一个容易被忽视的是「测试方程太简单，和真实问题无关」的想法。y'=λy 看似玩具，但任何常微分方程在平衡点附近线性化，ヤコビ行列的每个特征值都扮演一个 λ 的角色。也就是说，稳定性分析实际是在检查「对所有特征值 λ_i，z=hλ_i 都落在稳定域内吗」。最大模的特征值支配着允许的时间步——这正是刚性的本质。测试方程的稳定理论，直接对应到实际问题中ヤコビ矩阵的分析，用处很大。

使用指南

  1. 在左侧设置特征值 λ（lamNum/lamRange，范围 -5 到 +5）和步长 h（hNum/hRange，范围 0.01 到 0.5）
  2. 模拟器计算复平面上 z=hλ 的位置，以及各Runge-Kutta法（欧拉法·RK2·RK4）的放大率 |R(z)|
  3. 判断稳定域（|R(z)|≤1）或非稳定，比较常微分方程 dy/dt=λy 的数值解和精确解 exp(λt)，在设定积分时间内（teNum/teRange）确认最大误差

具体计算例子

当 λ=-2.0、h=0.1、积分时间 t=1.0 时，z=hλ=-0.2。RK4法中 z=-0.2 在稳定域内（复平面左半平面），|R(z)|≈0.98，精确解 e^(-2.0)≈0.135，数值解≈0.134，误差 <0.001。相比之下，λ=-3.0、h=0.15 时 z=-0.45，欧拉法不稳定（|R(z)|=1.45>1），解发散；但RK4法仍保持 |R(z)|≈0.997，处于稳定状态。

实务中的注意事项

  1. 遇到包含负大特征值的刚性问题时，欧拉法的稳定域限制在实轴 z<-2，此时应改用隐式Runge-Kutta法
  2. 对于具有复特征值 λ=α±βi 的耦合系统，要确保 z=h(α±βi) 完全落在复平面的稳定域内，步长 h 由最大模特征值支配
  3. RK4法虽然稳定域广，但每步需4次函数求值；当精度要求不高时应考虑RK2（2次求值）以节省计算