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

龙格-库塔法模拟器 — RK4 与欧拉法精度对比

给定衰减率 k、步长 h、计算终点 x_end 和初值 y₀，本工具用经典 4 阶龙格-库塔法（RK4）和欧拉法同时数值积分常微分方程 dy/dx = −k·y，并实时与解析解 y = y₀·exp(−k·x) 对比。终点值和欧拉相对误差即时显示，累积绝对误差以 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 时欧拉误差爆炸式增长。

计算结果
—
RK4 终点值 y(x_end)
—
欧拉终点值 y(x_end)
—
解析值 y(x_end)
—
欧拉相对误差
解的对比（解析解 / RK4 / 欧拉）

横轴＝x∈[0, x_end]，纵轴＝y。白色粗线＝解析解 y = y₀·exp(−k·x)，蓝色＝RK4，红色＝欧拉。h 较小时 RK4 与解析解几乎重合，欧拉解整体偏下；h 增大后欧拉的偏离越发明显。

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

横轴＝x，纵轴＝log10|y_num − y_exact|。蓝色＝RK4 绝对误差，红色＝欧拉绝对误差。RK4 的误差比欧拉小数个数量级（h=0.1 时约 10⁵ 倍），且随 x 增大两者增长速率差异显著。

理论与主要公式

目标常微分方程及其解析解：

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

显式欧拉法（一阶，全局误差 O(h)）：

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

经典 4 阶龙格-库塔法（RK4，全局误差 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$ 为初值。欧拉法仅用起点斜率近似，RK4 以 Simpson 风格的权重平均 4 个斜率；步长 $h$ 减半时欧拉误差变 1/2，RK4 误差变 1/16。

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

🙋
课本上说 RK4 比欧拉法精度高，可到底高多少？我想看到具体的数字。
🎓
用本工具的默认值跑一下：k=0.5, h=0.1, x_end=4, y₀=10。在 x=4 处求 dy/dx = −0.5y 的值，解析值为 10·e⁻² = 1.3534。结果卡显示欧拉是 1.2851（相差约 −5.05 %），RK4 是 1.3534，与解析解吻合到小数点后第 4 位。下方的对数刻度误差图里，RK4 的误差在 10⁻⁵ 量级，欧拉则在 10⁻¹ 量级，相差近 10000 倍。
🙋
差一万倍？！两者都在解同一个 ODE，为什么差距这么大？
🎓
关键在于"区间内评估斜率的次数"。欧拉法只用区间起点 x_n 的斜率，Taylor 展开二阶及以上都成了误差，局部误差 O(h²)、全局误差 O(h)。RK4 在起点、两次中点和终点共 4 个位置评估斜率，并以 Simpson 风格的权重 (1, 2, 2, 1)/6 加权平均，刚好和 Taylor 展开到 4 阶吻合，因此局部误差 O(h⁵)、全局误差 O(h⁴)。h=0.1 时，欧拉是 0.1 量级，RK4 则是 0.1⁴=10⁻⁴ 量级。
🙋
我按 play 扫描 h，发现 h 一过 0.5 欧拉曲线就开始震荡，最后发散了！RK4 还安然无恙。
🎓
这是"数值稳定性"。欧拉法的稳定条件是 |1 − kh| ≤ 1，对 k=0.5 来说 h=4 是边界，但在这之前已经开始震荡。RK4 的稳定区域要大得多，大约 |λh| ≤ 2.78。对于刚性问题（k 大、尺度差异大），欧拉法需要极小的步长才能稳定，因此实际工程中大家会换成隐式方法 — 比如 Abaqus/Standard 的 Hilber-Hughes-Taylor（α 法）、LS-DYNA 的中心差分等。
🙋
那 RK4 是不是最厉害？默认就用 RK4 不就行了？
🎓
没那么简单。RK4 每步要算 4 次 f，等于欧拉法的 4 倍计算量。所以"等计算量"对比时，应当让欧拉用 h/4 与 RK4 的 h 比较 — 即便如此，RK4 仍胜出几个数量级。工程实务的标准其实是带步长控制的 Dormand-Prince RK45：MATLAB ode45、SciPy RK45、Mathematica NDSolve 都用它。RK4 是教科书经典，也是理解这些自适应方法的基石。

常见问题

欧拉法仅用区间起点 x_n 的斜率来近似区间 [x_n, x_n+h] 上的变化率，所以局部误差为 O(h²)，全局误差仅为 O(h)。而 RK4 在区间的起点、两次中点和终点共 4 个位置评估斜率，再以类似 Simpson 公式的权重 (1, 2, 2, 1)/6 加权平均，达到局部误差 O(h⁵)、全局误差 O(h⁴)。本工具的默认参数（k=0.5, h=0.1, x_end=4, y₀=10）下，欧拉法得到 y(4)=1.2851（相对解析值 1.3534 偏差 −5.05 %），而 RK4 输出 1.3534（误差小于 1e-9），相差近万倍。若要让欧拉法达到 RK4 的精度，需把 h 缩小到原来的万分之一，计算成本约为 RK4 的数千倍。
步长 h 越大，每步计算量越小，速度越快；但局部误差按 h 的幂次增长。欧拉法误差与 h 一次方成正比（一阶），RK4 与 h⁴ 成正比（四阶），所以将 h 翻倍时欧拉误差变 2 倍，RK4 误差变 16 倍。另外，对于 dy/dx = −k·y 这样的刚性问题，当 |1 − kh| > 1 时欧拉解会发散（k=0.5 时 h > 4 出现震荡和发散）。点击 play 按钮把 h 从 0.005 扫至 1.0，可以观察到 h≈0.5 附近欧拉解开始震荡，h>2 时彻底崩溃；而 RK4 在 h ≈ 2.78 之内仍稳定，因此即便在中等刚性问题上也实用。
常用的 ODE 数值解法包括：(1) 显式欧拉法（一阶，最简单）；(2) 中点法（二阶）；(3) Heun 法 / 改进欧拉法（二阶）；(4) RK4（四阶，本工具实现）；(5) Adams-Bashforth 多步法（重用过去值以节省评估）；(6) 后向欧拉、梯形法和 BDF（隐式，适合刚性问题）；(7) Dormand-Prince RK45（自适应步长，MATLAB ode45 和 SciPy RK45 的默认）。在 CAE 实务中，LS-DYNA 显式动力学使用中心差分法（二阶），Abaqus/Standard 隐式动力学使用 Hilber-Hughes-Taylor（α 法）。RK4 是教科书经典，但工程实际中带误差控制的 RK45（Cash-Karp / Dormand-Prince）才是标准。
本工具仅处理一维线性 ODE dy/dx = −k·y，不实现：(1) 耦合 ODE 系统（如 Lorenz 系统的 dx/dt, dy/dt, dz/dt）；(2) 非线性 ODE（Van der Pol、Lotka-Volterra）；(3) 偏微分方程（PDE）的时间积分（热传导、波动）；(4) 自适应步长控制（RK45）；(5) 刚性自动检测与隐式求解器切换；(6) 辛积分（Hamilton 系统的能量守恒）。实际 CAE 求解器将这些思路组合使用，例如结构动力学的 Newmark-β 法、CFD 的 Crank-Nicolson 格式、轨道计算的 Verlet 积分。本工具最适合用于直观地理解收敛阶概念。

实际应用

结构动力学（FEM 瞬态分析）：地震响应或冲击分析需对运动方程 M·ü + C·u̇ + K·u = F(t) 进行时间积分。Abaqus/Standard 使用 Hilber-Hughes-Taylor（HHT-α，隐式 2 阶）法；LS-DYNA 使用中心差分法（显式 2 阶）。显式法无需求逆刚度矩阵，特别适合碰撞分析；隐式法可使用更大步长，适合长时程响应。这正是本工具中欧拉法 vs RK4 在"精度-稳定性-计算量"之间权衡的物理写照。Newmark-β 法（β=1/4, γ=1/2）无条件稳定，是商用 CAE 中最常见的积分格式。

CFD 时间推进：OpenFOAM 的 icoFoam（非定常层流）默认采用 Euler 隐式或 Crank-Nicolson 格式，SU2 则默认提供 RK4 选项。Navier-Stokes 方程本质上是非线性、刚性的，因此 CFL（Courant-Friedrichs-Lewy）条件限制了时间步长。本工具中 h 增大时欧拉发散的现象，与 CFD 求解器在 CFL>1 时发散是完全相同的物理事件。LES（大涡模拟）通常采用 TVD-RK3（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 的仿真、卡尔曼滤波的预测步 — 都是 ODE 积分的应用。实时控制器为了满足时间约束，常将步长固定在 1 ms 左右使用欧拉法（FPGA 易实现）；而 CARLA 自动驾驶模拟、MATLAB Simulink 等离线分析默认使用 ode45（RK45）。本工具中拖动 k、h 滑块即可体会一个经验法则："控制带宽 10 倍以上的步长"（k·h ≪ 0.1）是安全的。

常见误区与注意事项

最常见的误区是"RK4 总是最精确"。RK4 的全局误差为 O(h⁴)，写成 C·h⁴ 时常数 C 依赖于解的高阶导数。对于强振荡解（高频分量大）或右端不光滑（包含 Heaviside 阶跃）的问题，C 会变得很大，RK4 的表观精度反而下降。例如含冲击的结构分析中通常不用 RK4，而选择保波动的 LeapFrog 法或 TVD-RK 方法。本工具中的平滑指数衰减最适合 RK4，但实际应用中应根据解的性质来匹配方法。

其次常见的误解是"步长 h 越小越好"。h 减小可降低离散化误差，但计算步数增加会累积浮点舍入误差。双精度浮点（IEEE 754 double）的机器精度约为 2.2e-16，h 小于 1e-8 时舍入误差将主导离散化误差，精度反而下降（最优 h 约 √(2ε) ≈ 1e-8）。本工具滑块的最小 h=0.005 在最优范围之上；若要试更小 h，请直接修改源码。

最后是"RK4 可用于刚性问题"的误解。刚性问题（时间常数相差几个数量级的 ODE 系统）下，RK4 的稳定区域 |λh| ≤ 2.78 成为限制。例如化学反应中 τ₁=1 s, τ₂=1e-6 s 的双组分系统，需要 h ≪ 1 μs 才能跟上 τ₂，相对于感兴趣的秒级尺度计算量过大。刚性问题必须用隐式方法（后向欧拉、BDF、Rosenbrock）求解 — MATLAB 的 ode15s、SciPy 的 BDF 和 Radau 都属此类。这就是为何 CAE 求解器要求选择 Implicit 或 Explicit 模式的物理原因。

相关工具

常微分方程数值方法比较 常微分方程数值方法在线比较器：免费使用Euler法、Heun法（RK2）与四阶龙格-库塔法（RK4），实时计算并可视化指数衰减、逻辑斯谛增长、谐振子等模型的数值解与全局误差。 常微分方程式相平面分析工具 探索相平面分析工具：实时可视化方向场、零线与RK4相轨迹，集成摆、捕食者-猎物等经典模型。通过Jacobian矩阵进行稳定性分析，掌握平衡点分类（稳定点、鞍点等）与线性化方法，结合R... 欧拉法与RK4对比模拟器 交互式常微分方程求解器，对比前向欧拉法、改进欧拉法（霍恩法）与四阶龙格-库塔法（RK4）的精度。通过对数坐标图直观展示O(h)、O(h²)、O(h⁴)误差收敛特性，适合CAE数值分析... 拉丁超立方采样 — LHS 与朴素蒙特卡罗对比 拉丁超立方采样（LHS）模拟器将二维分层采样与朴素蒙特卡罗并排比较，实时可视化相同样本数下积分精度提升的程度，并直观展示每行每列各放置一个采样点的机理。 数值积分法比较工具（梯形·辛普森·高斯） 比较梯形法、辛普森法、高斯求积与龙贝格积分的数值积分工具。通过实时log-log图可视化误差收敛阶次，直观展示各方法误差阶：梯形法O(h²)、辛普森法O(h⁴)、高斯n点法O(h²ⁿ... 高斯求积法模拟器 — 高斯-勒让德积分的收敛 高斯求积法模拟器可视化高斯-勒让德积分的采样点与权重，实时计算改变点数时数值积分的误差与收敛。选择被积函数和积分区间，直观学习为什么用少量点数就能获得高精度。