斯托默—韦尔莱法

斯托默—韦尔莱法（Störmer–Verlet method）是求解位置和动量演化的二阶时间积分法，常用于经典力学和分子动力学。对可分离 Hamiltonian $H(q,p)=T(p)+V(q)$，它交替推进动量与位置，保持辛结构，并且具有时间可逆性。它通常使能量在一个小范围内振荡，而不是在每一步都精确保留原能量。

一个弹簧振子的三次更新

考虑无量纲系统 $q^{\prime }=p,p^{\prime }=-q$，初始位置 $q_0=1$、动量 $p_0=0$。精确轨迹为 $q(t)=\cos t,p(t)=-\sin t$。把时间步长记为 $h$。动量—位置—动量次序的 Verlet 法为 $$p_{n+1/2}=p_n-\frac{h}{2}{q}_n,q_{n+1}=q_n+h{p}_{n+1/2},p_{n+1}=p_{n+1/2}-\frac{h}{2}{q}_{n+1}.$$ 半步动量不是额外的物理粒子，而是组合两个半步作用时使用的中间值。

取 $h=0.2$，首步先得 $p_{1/2}=-0.1$，所以 $q_1=0.98$，再得 $p_1=-0.198$。第二步为 $p_{3/2}=-0.296$、$q_2=0.9208$、$p_2=-0.38808$。这些数值可与 $(\cos 0.4,-\sin 0.4)\approx (0.921060994,-0.389418342)$ 比较。两次半步都要使用各自时刻的位置；将末次力误用成旧位置，会退化成另一种方法。

从力学方程得到一般格式

对于质量 $m>0$、位置 $q\in {ℝ}^d$、动量 $p\in {ℝ}^d$，设 $T(p)=p^{T}p/(2m)$，势能 $V$ 可微。方程为 $q^{\prime }=p/m,p^{\prime }=-\mathrm{\nabla }V(q)$。保持位置不动时，动量在时间 $h/2$ 内精确改变 $-(h/2)\mathrm{\nabla }V(q)$；保持动量不动时，位置在时间 $h$ 内精确改变 $hp/m$。组合这两个可精确完成的子问题，得到 $$p_{n+1/2}=p_n-\frac{h}{2}\mathrm{\nabla }V(q_n),q_{n+1}=q_n+\frac{h}{m}{p}_{n+1/2},p_{n+1}=p_{n+1/2}-\frac{h}{2}\mathrm{\nabla }V(q_{n+1}).$$ 这也是速度 Verlet 的动量写法。一般可分离动能可将中间式改成 $q_{n+1}=q_n+h\mathrm{\nabla }T(p_{n+1/2})$。有速度依赖力、摩擦或约束时，需要重新推导算法，不能直接套这组三式。Hairer、Lubich 与 Wanner 的论文将其作为对称辛分裂的基本例子。

二阶精度、辛性与可逆性

先在一维说明精度。令 $F(q)=-V^{\prime }(q)$，并假设真解和力足够光滑。位置更新给出 $$q_{n+1}=q_n+h\frac{p_n}{m}+\frac{h^2}{2m}F(q_n),$$ 与真解的 Taylor 展开到二阶一致。把 $q_{n+1}=q_n+h{p}_n/m+O(h^2)$ 代入力的展开，得到 $$F(q_{n+1})=F(q_n)+h{F}^{\prime }(q_n)p_n/m+O(h^2).$$ 因此动量更新为 $$p_{n+1}=p_n+hF(q_n)+\frac{h^2}{2m}{F}^{\prime }(q_n)p_n+O(h^3),$$ 也与真解到二阶一致。局部状态误差为 $O(h^3)$；在固定有限时间、解留在光滑有界区域并有通常的稳定性估计时，全局误差为 $O(h^2)$。

辛性可以由每个子步直接证明。一维动量剪切的 Jacobian 为 $$K=\left(\begin{array}{cc}1& 0\\ -(h/2)V^{″}(q)& 1\end{array}\right),$$ 位置剪切的 Jacobian 为 $D=\left(\begin{array}{cc}1& h/m\\ 0& 1\end{array}\right)$。二者均满足 $M^{T}JM=J$，其中 $J=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$。组合的 Jacobian 是三个因子的乘积；逐层代入等式，就仍有 $M^{T}JM=J$。多维时用 Hessian 的对称性完成同样计算。

将 $h$ 换成 $-h$，三个子步恰好逆序抵消原来的三个子步，故 ${\Psi }_{-h}={\Psi }_h^{-1}$。这一可逆性指精确算术下的离散映射，不保证浮点舍入可以逐位倒退，也不意味着带耗散的真实过程应当时间可逆。

谐振子取 h=0.2。曲线显示原能量 H 与精确保留的二次修正量；纵轴放大后能看见原能量并非逐步守恒。

谐振子的稳定性与修正能量

消去半步变量，更新矩阵为 $$\left(\begin{array}{c}{q}_{n+1}\\ p_{n+1}\end{array}\right)=\left(\begin{array}{cc}1-h^2/2& h\\ -h(1-h^2/4)& 1-h^2/2\end{array}\right)\left(\begin{array}{c}{q}_n\\ p_n\end{array}\right).$$ 行列式等于一。其特征根在 $0<h<2$ 时为单位圆上的不同共轭根；$h>2$ 则有根模大于一，$h=2$ 的重根也可能产生线性增长。所以辛结构并不替代步长限制。

记 $c=1-h^2/2$、$b=1-h^2/4$，则 $c^2+h^{2}b=1$。由 $q_{+}=cq+hp,p_{+}=-hbq+cp$ 展开，可见交叉项相消，且 $$p_{+}^2+b{q}_{+}^2=(h^2{b}^2+b{c}^2)q^2+(c^2+b{h}^2)p^2=b{q}^2+p^2.$$ 因此 $\tilde{H}=(p^2+(1-h^2/4)q^2)/2$ 精确保留。初值下它等于 $0.495$，而首步原能量为 $(0.98^2+0.198^2)/2=0.499802$。在 $h<2$ 时，修正二次型正定并控制轨迹；到稳定边界它失去正定性，正好解释步长过大时的失败。

一个可自行重建的变式是常力 $F=2$、$m=1,q_0=p_0=0,h=0.2$。两个半步动量分别得到 $0.2,0.4$，中间位置为 $0.04$，与真解 $q=t^2,p=2t$ 完全一致。这来自常加速度的二次位置函数，不是一般力场中的精确性。

可运行的标量实现

代码输入势能导数、质量、初态、步长和步数，返回位置动量序列，不修改输入。每步调用两次势能导数；缓存相邻步共有的力可省去一次。若一次求力视为常数成本，时间为 $O(N)$，轨迹存储为 $O(N)$。多粒子模拟的主要代价往往在求力本身。

Python 3

import math


def verlet(grad_v, mass, q, p, h, steps):
    mass, q, p, h = map(float, (mass, q, p, h))
    if not all(math.isfinite(x) for x in (mass, q, p, h)) or mass <= 0 or h <= 0:
        raise ValueError("finite data, positive mass and step required")
    if not isinstance(steps, int) or steps < 0:
        raise ValueError("nonnegative integer steps required")
    states = [(q, p)]
    for _ in range(steps):
        force0 = -float(grad_v(q))
        half = p + h * force0 / 2
        q_new = q + h * half / mass
        if not all(math.isfinite(x) for x in (force0, half, q_new)):
            raise OverflowError("nonfinite intermediate state")
        p_new = half - h * float(grad_v(q_new)) / 2
        if not math.isfinite(p_new):
            raise OverflowError("nonfinite momentum")
        q, p = q_new, p_new
        states.append((q, p))
    return states


states = verlet(lambda q: q, 1, 1, 0, 0.2, 2)
print(states)  # [(1.0, 0.0), (0.98, -0.198), (0.9208, -0.38808)]
print(verlet(lambda q: -2, 1, 0, 0, 0.2, 1)[-1])  # (0.04, 0.4)

对非线性势能，通常需要用修正方程分析长期行为，不能把上面的精确二次不变量照搬。固定步长的相位误差也会累积：轨道保持有界，并不表示每个时刻的位置都接近无限长时间的真解。