› 龙格-库塔法 模拟器 返回
数值分析模拟器

龙格-库塔法 模拟器 — RK4 vs Euler 精度对比

从衰减率 k、步长 h、计算终点 x_end、初值 y₀ 出发，用经典 4 阶龙格-库塔法（RK4）与 Euler 法求解常微分方程 dy/dx = −k·y，并与精确解 y = y₀·exp(−k·x) 并排对比。实时显示终点值和 Euler 相对误差，用 log10 尺度可视化累积误差，直观理解"收敛阶数"的含义。

参数设置
衰减率 k
1/x
步长 h
—
计算终点 x_end
—
初值 y₀
—

默认值：k=0.5, h=0.1, x_end=4, y₀=10。精确解 y(4) = 10·e⁻² ≈ 1.3534。增大 h 时 Euler 误差会爆发性增长。

计算结果
—
RK4 终点值 y(x_end)
—
Euler 终点值 y(x_end)
—
精确解 y(x_end)
—
Euler 相对误差
解的对比（精确解 / RK4 / Euler）

横轴=x [0, x_end] / 纵轴=y / 黑粗线=精确解 y = y₀·exp(−k·x) / 蓝实线=RK4 / 红实线=Euler。h 越小，RK4 与精确解越接近，Euler 偏离越明显。h 越大，Euler 偏差越显著。

累积误差 |y_num − y_exact|（log10 尺度）

横轴=x / 纵轴=log10|y_num − y_exact| / 蓝色=RK4 的绝对误差 / 红色=Euler 的绝对误差。RK4 比 Euler 小数个数量级（h=0.1 时约相差 10⁵ 倍精度）。随 x 增加两者误差都增长，但增长率差异显著。

理论与主要公式

目标常微分方程与精确解：

$$\frac{dy}{dx} = -k\,y,\qquad y(0)=y_0 \;\Longrightarrow\; y(x) = y_0\,e^{-kx}$$

显式 Euler 法（1 阶，全局误差 O(h)）：

$$y_{n+1} = y_n + h\,f(x_n, y_n) = y_n(1-kh)$$

经典 4 阶龙格-库塔法（4 阶，全局误差 O(h⁴)）：

$$\begin{aligned}k_1 &= h\,f(x_n, y_n)\\ k_2 &= h\,f(x_n+\tfrac{h}{2}, y_n+\tfrac{k_1}{2})\\ k_3 &= h\,f(x_n+\tfrac{h}{2}, y_n+\tfrac{k_2}{2})\\ k_4 &= h\,f(x_n+h, y_n+k_3)\\ y_{n+1} &= y_n + \tfrac{1}{6}(k_1+2k_2+2k_3+k_4)\end{aligned}$$

$k$ 是衰减率（时间常数 $\tau=1/k$），$h$ 是步长，$y_0$ 是初值。Euler 仅用起点斜率近似，RK4 则用 4 点的 Simpson 型加权平均。将 h 减半时，Euler 误差减至 1/2，RK4 误差减至 1/16。

龙格-库塔法模拟器是什么

🙋
大学课堂上听说"龙格-库塔法的精度比 Euler 法高"，但到底高多少呢？想用数字实际体验一下。
🎓
很好的问题。试试本工具的默认值（k=0.5, h=0.1, x_end=4, y₀=10）。求解 dy/dx = −0.5y，在 x=4 时的 y 值上比较，精确解为 10·e⁻² = 1.3534，Euler 得到 1.2851，RK4 得到 1.3534。Euler 偏低约 5.05%，而 RK4 与精确解四位小数完全吻合。下面的对数尺度误差图表显示，RK4 误差为 10⁻⁵ 量级，Euler 为 10⁻¹ 量级，相差 10000 倍！
🙋
竟然相差 10000 倍！为什么同样求解一个常微分方程会有这么大的差异？
🎓
区别在于"在区间内评估几次斜率"。Euler 法仅在区间起点 x_n 评估斜率进行近似，所以 1 步前进会产生 2 阶以上的误差（局部误差 O(h²)，全局误差 O(h)）。RK4 在起点、中点（两次）、终点共 4 个位置评估斜率，再用 Simpson 法则型的权重 (1, 2, 2, 1)/6 加权平均。这样就能与 4 阶 Taylor 展开相匹配，局部误差达 O(h⁵)，全局误差达 O(h⁴)。当 h=0.1 时，Euler 误差约 0.1 阶，RK4 误差约 0.1⁴=10⁻⁴ 阶。
🙋
我用 play 按钮让 h 从小到大扫描，看到 h 超过 0.5 后 Euler 的图线开始波动，最后甚至发散了！但 RK4 似乎没事。
🎓
这就是"数值稳定性"。Euler 法的稳定条件是 |1 − kh| ≤ 1，即 h ≤ 2/k。当 k=0.5 时上限是 h=4，但更小的 h 值下也会出现振荡。而 RK4 的稳定域更宽广，约 h ≈ 2.78/k 以下都稳定。对于刚性问题（k 很大、时间尺度差异极端），Euler 必须大幅缩小步长否则爆炸，这就是为什么实际应用中不用 Euler。在 CAE 中，Abaqus/Standard 隐式求解器采用 Hilber-Hughes-Taylor 法（α 法），LS-DYNA 显式求解器采用中央差分法，都是针对各自问题特性优化的。
🙋
那 RK4 就最强了吧？我们应该总是用 RK4？
🎓
没那么简单。RK4 每步要计算 4 次函数，即 Euler 的 4 倍计算代价。在"相同计算量"的条件下，应该比较 Euler 用 h/4 与 RK4 用 h 的精度，结果还是 RK4 赢（精度差 10³ 倍）。不过在实际工程应用中，标准做法是用 Dormand-Prince RK45，它不仅是 5 阶还能自动调整步长。MATLAB 的 ode45、SciPy 的 RK45、Mathematica 的 NDSolve 都是这个。教科书上的 RK4 是"初学者的经典入门"定位。

常见问题

Euler 法在区间 [x_n, x_n+h] 内仅在起点 x_n 处计算斜率进行近似，导致局部误差为 O(h²)，全局误差为 O(h)。而 RK4 在起点、中点（2 次）、终点共 4 个位置评估斜率，用类似 Simpson 公式的权重 (1, 2, 2, 1)/6 进行加权平均，从而达到局部误差 O(h⁵)、全局误差 O(h⁴)。在本工具的默认值（k=0.5, h=0.1, x_end=4, y₀=10）下，Euler 得到 y(4)=1.2851（相对于精确值 1.3534 偏低 5.05%），而 RK4 得到 1.3534（误差小于 1e-9），精度差异达四位数量级。要用 Euler 法获得相同精度，需将 h 缩小 10000 倍，计算成本将增加数千倍。
增大步长 h 可减少每步的计算量，提高运算速度，但会导致局部误差随 h 的幂次增长。Euler 法的误差与 h 成一次比例（1 阶），RK4 与 h⁴ 成比例（4 阶），因此 h 翻倍时，Euler 误差翻倍，RK4 误差增至 16 倍。对于 dy/dx = −k·y 这样的刚性问题，当 |1 − kh| > 1 时 Euler 解会发散（k=0.5 时 h > 4 时出现振荡和发散）。通过本工具的 play 按钮将 h 从 0.005 扫描至 1.0，可观察到 h ≈ 0.5 时 Euler 开始振荡，h>2 时完全破裂，而 RK4 在 h ≈ 2.78 以下保持稳定，因此在刚性问题中更具实用性。
常见的常微分方程求解器包括：(1) 显式 Euler 法（1 阶，最简单）、(2) 中点法（2 阶）、(3) Heun 法/改进 Euler 法（2 阶）、(4) RK4（4 阶，本工具实现）、(5) Adams-Bashforth 多步法（复用过去值，降低计算量）、(6) 隐式 Euler 法、梯形法则、BDF（隐式方法，对刚性问题有利）、(7) Dormand-Prince RK45（自适应步长，MATLAB ode45 与 SciPy RK45 的标准采用）。在 CAE 实践中，结构动力学的显式求解器 LS-DYNA 采用中央差分法（2 阶），Abaqus/Standard 隐式求解器采用 Hilber-Hughes-Taylor 法（α 法）。本工具的 RK4 是教科书经典，但在工程应用中通常采用具有误差控制的 RK45（Cash-Karp/Dormand-Prince）。
本工具仅涉及 1 维线性常微分方程 dy/dx = −k·y，不处理：(1) 常微分方程组（Lorenz 系统的 dx/dt、dy/dt、dz/dt 等）、(2) 非线性常微分方程（Van der Pol、Lotka-Volterra）、(3) 偏微分方程（PDE）的时间积分（热传导、波动）、(4) 自适应步长控制（RK45）、(5) 刚性性检测与隐式方法的自动切换、(6) 辛积分（Hamilton 系统的能量守恒）。在真实 CAE 求解器中，这些方法相互结合使用，例如结构动力学的 Newmark-β 法、流体动力学的 Crank-Nicolson 格式、轨道力学的 Verlet 积分等根据具体情况选择应用。本工具最适合用于直观理解收敛阶数概念的教学用途。

现实世界中的应用

结构动力学（有限元瞬态响应分析）：地震动力分析或冲击响应分析需要对运动方程 M·ü + C·u̇ + K·u = F(t) 进行时间积分。Abaqus/Standard 采用 Hilber-Hughes-Taylor 法（HHT-α，隐式 2 阶精度），LS-DYNA 采用中央差分法（显式 2 阶精度）。显式法不需逆矩阵计算，对碰撞分析有利；隐式法用更大步长处理长时间响应 —— 这正是本工具 Euler 和 RK4 比较的"精度·稳定性·计算成本"三角拉锯战。Newmark-β 法（β=1/4, γ=1/2）因无条件稳定成为商用 CAE 最常用的时间积分方案。

流体动力学（CFD 时间积分）：OpenFOAM 的 icoFoam（非定常层流）可选 Euler 隐式或 Crank-Nicolson 法，SU2 提供 RK4 作为标准选项。Navier-Stokes 方程本质非线性且具有刚性，受 CFL 条件（Courant-Friedrichs-Lewy）约束步长。本工具中 h 增大导致 Euler 发散的现象，与 CFD 中 CFL > 1 引发求解器崩溃是完全相同的物理机制。大涡模拟（LES）采用 RK3-TVD（3 阶 Runge-Kutta with Total Variation Diminishing）在激波附近抑制振荡。

轨道力学（卫星与行星运动）：NASA 的 SPICE Toolkit 或 GMAT 用 RK7(8) 或 Bulirsch-Stoer 法计算卫星轨道。地球周围卫星需要每圈 0.1 秒精度，长期积分（10 年跨度）采用 Verlet 法或辛积分保证能量守恒。本工具的精确解 y = y₀·e⁻ᵏˣ 虽简单，但 Kepler 轨道（二体问题）也有精确解，与 RK4 比较时同样表现出"短期 RK4 高精度，长期辛方法优"的权衡。

控制工程（实时控制积分）：PID 控制的积分环节、状态空间模型 ẋ = Ax + Bu 的仿真、Kalman 滤波预测步 —— 都是常微分方程积分的应用。实时控制因计算时间受限常用步长固定的 Euler 法（FPGA 实现也更轻量）。而自动驾驶仿真器 CARLA 或 Matlab Simulink 的离线分析采用 ode45（RK45）作为标准。本工具中滑动 k 和 h 的体验表明，"控制带宽的 10 倍步长"（k·h ≪ 0.1）是安全准则。

常见误解与注意事项

最常见的误解是"使用 RK4 总能保证高精度"。RK4 的全局误差理论上为 O(h⁴)，但形式为"常数 C × h⁴"，其中常数 C 依赖于解的高阶导数。对于振荡剧烈的解（含高频分量）或右端项不连续（包含 Heaviside 函数）的情况，C 值会很大，表观精度下降。例如含冲击的结构分析中，RK4 不是最优选择，应改用波动保存的 LeapFrog 法或 TVD-RK 方案。本工具针对光滑指数衰减，RK4 最优，但面向应用问题必须根据解的性质选择合适方法。

第二常见误解是"步长 h 越小越好"。减小 h 能降低离散化误差，但增加步数会累积浮点舍入误差。IEEE 754 双精度浮点的机器精度约 2.2e-16，当 h 小于 1e-8 时舍入误差超过离散化误差，反而恶化精度（最优 h 约 √(2ε) ≈ 1e-8）。本工具的 h 滑块最小值 0.005 在最优范围内，若要尝试更小值可修改源代码。

第三常见误解是"RK4 即使对刚性问题也适用"。刚性问题（时间常数相差悬殊的方程组）中，RK4 的稳定域 |λh| ≤ 2.78 成为瓶颈。例如化学反应中 τ₁=1 s, τ₂=1e-6 s 的两分量系统，必须按 τ₂ 调整 h ≪ 1 μs，导致对秒级现象的计算量成为不可接受。刚性问题必用隐式方法（后退 Euler、BDF、Rosenbrock）。MATLAB 的 ode15s、SciPy 的 BDF/Radau 属此类。CAE 软件要求选择"隐式/显式"的物理根源正在于此。

关联工具

常微分方程·相平面 分析模拟器 常微分方程的相平面分析在线无费体验。2维自治系统的方向场、零斜线、平衡点稳定性实时可视化。振摆、Lotka-Volterra、SIR 模型等多种动态系统模拟... 拉丁超立方抽样 — LHS 与简单蒙特卡洛对比 拉丁超立方抽样（LHS）模拟器展示 2 维分层抽样与简单蒙特卡洛的比较，在相同样本数下数值积分精度改善程度的实时可视化。每行、每列各一个样本点实现拉丁超立方性质... 数值积分法对比 模拟器 数值积分法对比工具实时可视化梯形法则、Simpson 法则、Gauss 求积、Romberg 积分的误差收敛。各方法理论误差阶（如梯形法 O(h²)、Simpson 法 O(h⁴)）与实际对比... Gauss 求积法模拟器 — Gauss-Legendre 积分的收敛 Gauss 求积法模拟器可视化 Gauss-Legendre 积分的采样点与权重，实时计算点数变化时数值积分误差与收敛过程。选择函数与积分区间，体验为何少量点数可达高精度...

使用指南

  1. 设置初值条件：对微分方程 dy/dx = -ky 输入 k 值（slKVal）在 0.5～2.0 范围，初值 y0（slY0Val）设为 1.0
  2. 调整步长 h（slHVal）：分别选择 h=0.1、0.05、0.01 三个梯次进行对比，观察 RK4 与 Euler 法的误差消减过程
  3. 指定积分终点 xend（slXendVal）并运行仿真：自动计算 RK4 终点值、Euler 终点值、精确解 y(x_end)=y0×e^(-kx_end)

具体计算示例

放射性物质衰变（k=0.693、半衰期 1.0 秒）在 x_end=2.0 秒、y0=100mg、h=0.1 秒的情况下：精确解为 25.0mg，Euler 法得 26.8mg（相对误差 7.2%），RK4 法得 25.02mg（相对误差 0.08%），RK4 约高 90 倍精度。若将 h 细化至 0.05，Euler 误差降至 3.6%，但计算量翻倍。

实务中的注意点