4. 数值积分方法入门

在粒子运动的数值模拟中，我们通常已知粒子的初始位置和速度，以及作用在粒子上的力（或加速度），希望求出之后每个时刻的位置与速度。牛顿第二定律给出的是微分方程：a = F/m，即 dv/dt = a、dx/dt = v。数值积分要做的，就是把连续的时间切成一个个小步长 Δt，用离散的递推公式近似这个连续过程。本节聚焦最基础也最常用的两种方法：显式欧拉法与半隐式欧拉法，并通过一个自由落体加空气阻力的例子，观察它们的实现方式与误差表现。

一、显式欧拉法：用当前状态直接外推

显式欧拉法（Forward Euler）的思想最直接：用当前时刻的速度更新位置，用当前时刻的加速度更新速度。设第 n 步的位置为 xₙ、速度为 vₙ，则

vₙ₊₁ = vₙ + a(xₙ, vₙ) · Δt
xₙ₊₁ = xₙ + vₙ · Δt

注意这里位置更新用的是旧速度 vₙ，而不是更新后的 vₙ₊₁。这正是显式欧拉法的关键特征，也是它误差偏大的根源。从泰勒展开看，真实解为 x(t+Δt) = x(t) + v(t)Δt + ½a(t)Δt² + …，而欧拉法只保留了前两项，截断误差为 O(Δt²)（单步），全局误差为 O(Δt)。这意味着步长减半，误差大约也减半，收敛较慢。

更严重的问题出现在振动或轨道类运动中。以简谐振子为例，显式欧拉法每一步都会让系统的“数值能量”单调增加，粒子轨迹呈向外发散的螺旋，即使真实运动是等幅振荡。这种数值不稳定性使得显式欧拉法很少用于长时间模拟。

二、半隐式欧拉法：先更新速度，再用新速度更新位置

半隐式欧拉法（Semi-implicit Euler，也称辛欧拉法）只调整了更新顺序：

vₙ₊₁ = vₙ + a(xₙ, vₙ) · Δt
xₙ₊₁ = xₙ + vₙ₊₁ · Δt

差别仅在于位置更新改用新速度 vₙ₊₁。这个看似微小的改动带来了显著改善：对简谐振子，半隐式欧拉法保持数值能量在一个有界范围内振荡，不会系统性发散，因此特别适合力学、分子动力学和游戏物理中的粒子模拟。它的单步截断误差同样是 O(Δt²)，但具有更好的长期稳定性（属于辛积分器的一阶形式）。

三、两种方法的更新顺序与误差来源对比

下面的图示对比了两种方法在一个时间步内的状态更新路径，以及误差累积的定性趋势。可以看到，显式欧拉的位置更新滞后于速度变化，而半隐式欧拉让位置“跟上”了新速度，从而在振动问题中避免了能量漂移。

四、简单粒子运动模拟实现

考虑一个竖直方向运动的粒子：质量 m = 0.5 kg，初始高度 y₀ = 10 m，初速度 v₀ = 0，重力加速度 g = 9.8 m/s²，并加入与速度成正比的空气阻力 F_d = −k·v，取 k = 0.2 kg/s。规定竖直向上为正方向，则加速度为 a = −g − (k/m)·v。下面用半隐式欧拉法实现，步长 Δt = 0.01 s，模拟 2 s。

m, g, k = 0.5, 9.8, 0.2
y, v = 10.0, 0.0
dt = 0.01
t = 0.0
while t < 2.0:
    a = -g - (k/m)*v
    v = v + a*dt          # 先更新速度
    y = y + v*dt          # 再用新速度更新位置
    t = t + dt
    # 每 0.5 s 记录一次 (t, y, v)

若把最后两行换成“先 y = y + v·dt，再 v = v + a·dt”，就是显式欧拉法。两种方法在短时间内的结果接近，但随步数增加，差异会逐渐显现。判断模拟是否可信，可以从三个角度检查：

• 步长收敛性：把 Δt 从 0.01 s 减到 0.005 s，若结果变化明显，说明步长仍偏大；
• 能量或守恒量：对无阻尼振动，观察数值能量是否漂移；有阻尼时检查速度是否单调趋近终端速度 v_t = −mg/k ≈ −24.5 m/s；
• 与解析解对照：有阻力落体的速度解析解为 v(t) = −(mg/k)(1 − e^(−kt/m))，可在若干时刻对比数值结果。

五、误差分析与方法选择

两种方法的局部截断误差均为 O(Δt²)，全局误差为 O(Δt)。显式欧拉的误差主要来自位置更新滞后，在振动和轨道问题中表现为能量漂移；半隐式欧拉通过“速度先行”修正了这一点，在同等计算量下获得更好的长期行为。因此，对于粒子运动模拟，优先选择半隐式欧拉法；若需要更高精度，可进一步采用 Velocity-Verlet 等二阶方法。理解这两种入门方法，是后续学习更高级积分器的必要基础。