跳至主要內容

数值积分入门：前向 / 后向欧拉、RK4 与韦尔莱

StonyBear2026/7/29大约 6 分钟图形学数值积分物理模拟游戏开发

物理引擎与动画里，微分方程几乎无处不在：弹簧、刚体、布料、简单 AI 运动……核心往往是：已知 y˙=f(t,y)\dot y = f(t,y)（或二阶 x¨=a(x,v)\ddot x = a(x,v)），用离散时间步 Δt\Delta t 往前推状态。本文介绍四种常用方法：前向欧拉、后向欧拉、四阶龙格–库塔（RK4）、韦尔莱（Verlet），并在同一谐振子上对比它们的行为。

问题设定

一阶初值问题：

y˙=f(t,y),y(t0)=y0.\dot y = f(t,y),\qquad y(t_0)=y_0.

游戏里更常见的是二阶力学：

x¨=a(x,v,t),\ddot x = a(x,v,t),

可写成一阶系统 x˙=v, v˙=a\dot x = v,\ \dot v = a。下面四种方法里，前三种是通用一阶积分器；韦尔莱专为二阶力学设计，在保守力场里特别有用。

对比实验取无阻尼谐振子（ω=1\omega=1）：

x¨=−x,x(0)=1, x˙(0)=0,\ddot x = -x,\qquad x(0)=1,\ \dot x(0)=0,

解析解 x(t)=cos⁡tx(t)=\cos t。步长刻意取偏大的 Δt=0.35\Delta t = 0.35，好让稳定性差异一眼可见：

谐振子对比（Δt = 0.35）

图中可以看出：前向欧拉振幅爆炸；后向欧拉过度阻尼、振幅塌缩；RK4 几乎贴住解析解；Verlet 相位会漂一点，但振幅长期可控。

前向欧拉（显式欧拉）

yn+1=yn+Δt f(tn,yn).y_{n+1} = y_n + \Delta t\, f(t_n, y_n).

力学形式（半隐式变体很常见，但标准显式是）：

vn+1=vn+Δt a(xn,vn),xn+1=xn+Δt vn.\begin{aligned} v_{n+1} &= v_n + \Delta t\, a(x_n, v_n),\\ x_{n+1} &= x_n + \Delta t\, v_n. \end{aligned}

（若用 xn+1=xn+Δt vn+1x_{n+1}=x_n+\Delta t\, v_{n+1}，则是 semi-implicit Euler / symplectic Euler，稳定性更好，别和纯前向欧拉混为一谈。）

特点

阶局部截断 O(Δt2)O(\Delta t^2)，全局 O(Δt)O(\Delta t)
代价每步一次 ff，最便宜
稳定条件稳定；步长太大时谐振子能量单调增加
适用原型、弱刚性、对精度不敏感的效果

直觉：每步用当前斜率沿切线走出一整步，弯曲轨道上会系统性“冲出去”。

x := x + dt * v
v := v + dt * a(x_old, v)   // 注意：若 a 依赖新 x，顺序要写清楚

后向欧拉（隐式欧拉）

yn+1=yn+Δt f(tn+1,yn+1).y_{n+1} = y_n + \Delta t\, f(t_{n+1}, y_{n+1}).

未知量出现在等式两边，一般要解非线性方程（牛顿迭代）；对线性系统可直接求逆。上面谐振子例子里，把 (x,v)(x,v) 联立后得到闭式：

xn+1=xn+Δt vn1+ω2Δt2,vn+1=vn−ω2Δt xn+1.x_{n+1} = \frac{x_n + \Delta t\, v_n}{1 + \omega^2 \Delta t^2},\quad v_{n+1} = v_n - \omega^2 \Delta t\, x_{n+1}.

特点

阶仍是一阶，但 A-稳定：大步长通常不炸
代价每步要求解；刚性系统里往往值得
行为人为耗散强，振荡会被“捏死”
适用刚性约束、隐式约束求解、需要无条件稳定的场合

游戏里“完全后向欧拉 + 大规模非线性求解”成本高，更常见的是约束投影、XPBD、或半隐式方法；但理解后向欧拉有助于理解为什么隐式方法偏阻尼。

龙格–库塔四阶（RK4）

在 [tn,tn+Δt][t_n,t_n+\Delta t] 上取四个斜率再加权：

k1=f(tn,yn),k2=f ⁣(tn+Δt2, yn+Δt2k1),k3=f ⁣(tn+Δt2, yn+Δt2k2),k4=f ⁣(tn+Δt, yn+Δt k3),yn+1=yn+Δt6(k1+2k2+2k3+k4).\begin{aligned} k_1 &= f(t_n, y_n),\\ k_2 &= f\!\left(t_n+\tfrac{\Delta t}{2},\, y_n+\tfrac{\Delta t}{2}k_1\right),\\ k_3 &= f\!\left(t_n+\tfrac{\Delta t}{2},\, y_n+\tfrac{\Delta t}{2}k_2\right),\\ k_4 &= f\!\left(t_n+\Delta t,\, y_n+\Delta t\, k_3\right),\\ y_{n+1} &= y_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{aligned}

对状态向量 (x,v)(x,v)，把 ff 换成系统右端即可（上文图表即如此）。

特点

阶全局 O(Δt4)O(\Delta t^4)，同样 Δt\Delta t 通常远准于欧拉
代价每步 4 次 ff
稳定显式方法，稳定域比欧拉大，但不是无条件稳定
适用轨迹精度重要、力场光滑、可负担 4 倍力评估

经典“教科书精度”；实时物理里若力很贵（复杂碰撞、布料），会改用更便宜的 symplectic / Verlet 族。

韦尔莱积分（Verlet）

针对 x¨=a(x)\ddot x = a(x)（加速度主要依赖位置，速度阻尼需另处理）。基本（Störmer–Verlet）形式：

xn+1=2xn−xn−1+Δt2 a(xn).x_{n+1} = 2x_n - x_{n-1} + \Delta t^2\, a(x_n).

只需当前位置与上一步位置，不必显式存速度。速度可由中心差分事后估计：

vn≈xn+1−xn−12Δt.v_n \approx \frac{x_{n+1}-x_{n-1}}{2\Delta t}.

Velocity Verlet（与位置 Verlet 同族、同阶）更方便同时推进 x,vx,v：

xn+1=xn+Δt vn+12Δt2an,an+1=a(xn+1),vn+1=vn+12Δt (an+an+1).\begin{aligned} x_{n+1} &= x_n + \Delta t\, v_n + \tfrac12 \Delta t^2 a_n,\\ a_{n+1} &= a(x_{n+1}),\\ v_{n+1} &= v_n + \tfrac12 \Delta t\,(a_n + a_{n+1}). \end{aligned}

特点

阶位置全局约 O(Δt2)O(\Delta t^2)
结构辛（symplectic）/ 时间可逆（在理想设定下）→ 长期能量误差有界振荡，不易像前向欧拉那样发散
代价每步约一次 a(x)a(x)（Velocity Verlet 巧妙复用）
适用分子动力学、布料/软体、弹道、需要长时间稳定的保守系统

启动时若只有 x0,v0x_0,v_0，可用

x−1=x0−Δt v0+12Δt2a(x0)x_{-1} = x_0 - \Delta t\, v_0 + \tfrac12 \Delta t^2 a(x_0)

生成“上一步”，再进入递推。

怎么选

场景更合适
快速原型、视觉凑合前向欧拉 / 半隐式欧拉
刚性大、绝不能炸后向欧拉或其它隐式
要准、力不太贵RK4
长期力学、能量不能漂飞Verlet / Velocity Verlet
有强阻尼、非光滑碰撞单纯积分器不够，需配合约束/冲击求解

经验口诀：

  • 显式便宜，步长受稳定性限制；
  • 隐式稳，但耗散 + 求解成本；
  • 高阶 RK 准，但不保辛结构；
  • Verlet 为力学二阶方程“量身”，长时间表现往往优于同阶的朴素欧拉。

小结

  1. 前向欧拉：y+Δtfny+\Delta t f_n，简单，大步长易不稳定。
  2. 后向欧拉：y+Δtfn+1y+\Delta t f_{n+1}，稳但偏阻尼，要解方程。
  3. RK4：四次采样加权，精度高，力评估 ×4。
  4. Verlet：二阶力学专用，长时间振幅/能量更老实。

同一谐振子、同一大步长下，四条曲线已经把性格画清楚：会炸的、会瘫的、会贴解析解的、会慢慢相位漂移但仍“弹着”的。把这张图记在脑子里，读引擎 integrators 源码时会快很多。