Symplectic Integrators and the Shape of Energy Drift

Integrate a planetary orbit with a fourth-order Runge-Kutta method at a small step and the energy decreases monotonically. The orbit spirals inward, slowly at first, and the inward drift continues as long as the integration runs, because each step introduces a small dissipative error of the same sign. Use the second-order Störmer-Verlet method at the same step and the energy oscillates around its initial value without a trend. That difference is not a matter of accuracy. Verlet has a larger error per step than the fourth-order method, and it still produces the qualitatively correct long-term behavior while the more accurate method does not.

The structure that the equation possesses

A Hamiltonian system has the form

$$ \dot q = \frac{\partial H}{\partial p}, \qquad \dot p = -\frac{\partial H}{\partial q} $$

for a scalar function $H(q,p)$, the Hamiltonian, usually the total energy. The flow map $\varphi_t$ generated by the system preserves two structures: it is symplectic, meaning $(D\varphi_t)^T J (D\varphi_t) = J$ with $J$ the block matrix $\begin{pmatrix}0 & I\ -I & 0\end{pmatrix}$, and it preserves the level sets of $H$. The two are not independent. For a time-independent Hamiltonian, preservation of $H$ follows from symplecticity and the equivariance of the flow under time translation.

The practical consequence of symplecticity is a conservation law about areas in phase space. For a one-degree-of-freedom system the condition is that the map preserves areas, which rules out the monotone inward spiral: a flow that spirals inward must shrink areas, and a symplectic map cannot.

Verlet as a composition

The Störmer-Verlet method is defined by a splitting. Write the Hamiltonian as $H = T(p) + V(q)$ with $T$ depending only on momenta and $V$ only on positions. Each part is exactly integrable: the $T$ part moves $q$ at constant $p$, and the $V$ part changes $p$ at constant $q$. Verlet composes the three flows

$$ p_{n+1/2} = p_n - \frac{h}{2}\nabla V(q_n), \qquad q_{n+1} = q_n + h,\nabla T(p_{n+1/2}), \qquad p_{n+1} = p_{n+1/2} - \frac{h}{2}\nabla V(q_{n+1}) $$

This is a symmetric composition, and each of its factors is the exact flow of a Hamiltonian. The composition of exact Hamiltonian flows is symplectic, so Verlet is symplectic by construction rather than by coefficient design. For a separable Hamiltonian with $T = \frac{1}{2}p^T M^{-1} p$ and $V$ independent of $p$, the method needs one force evaluation per step and no storage beyond the current state.

Composition is a general recipe for building methods. Composing a second-order symmetric scheme with its own adjoint in the pattern $\varphi_{h/2} \circ \varphi^*_{h/2}$ gives fourth order, at the cost of several evaluations per step; the pattern of the composition determines the order, and there is a theory of composition methods based on the Baker-Campbell-Hausdorff expansion. The catch is that negative substeps appear in some of the higher-order compositions, and negative time steps behave badly for dissipative perturbations of the system.

Implementing the leapfrog form, which is algebraically equivalent when $V$ does not depend on $p$, needs three arrays:

/* Störmer-Verlet (velocity Verlet) for q'' = -grad V(q). */
for (size_t n = 0; n < steps; n++) {
    for (int i = 0; i < dim; i++) {
        p[i] -= 0.5 * h * gradV(q, i);     /* half kick */
        q[i] += h * p[i] / m[i];           /* drift */
        /* positions updated; recompute forces for the second half kick */
    }
    for (int i = 0; i < dim; i++)
        p[i] -= 0.5 * h * gradV(q, i);
}

The inner half-kick requires a force evaluation per half step; combining the second half-kick of step $n$ with the first half-kick of step $n+1$ halves the number of evaluations for a separable potential, which is why the leapfrog form stores the momenta at half-integer times.

The modified Hamiltonian and why the energy stays bounded

A symplectic method is not the exact flow of the original Hamiltonian, but it is the exact flow of a different one. Backward error analysis constructs a modified Hamiltonian $H_h = H + h H_1 + h^2 H_2 + \cdots$ whose exact flow, sampled at the step points, coincides with the method’s iterates up to an error that is exponentially small in $1/h$ over times that grow with the step.

The consequence for energy is direct: the method exactly conserves $H_h$, and near a regular orbit $H_h$ differs from $H$ by $O(h^p)$ where $p$ is the method’s order. So the original energy satisfies $H(t) = H(0) + O(h^p)$ along the trajectory, with oscillations of amplitude $O(h^p)$ but no secular drift. The energy error is bounded on exponentially long time intervals, not merely on short ones.

The bound comes with preconditions that are easy to violate. It requires the step to be small relative to the fastest frequency and the orbit to stay away from resonance; near a resonance the denominators in the expansion become small and the bound degrades, and the energy begins to wander. It also requires a regular orbit: chaotic trajectories have bounded energy error but their phase-space position diverges from the true one at a rate set by the Lyapunov exponent, so the trajectory is not accurate even though the energy is.

Why non-symplectic methods fail differently

A standard fourth-order Runge-Kutta method conserves neither energy nor a modified energy. Its local error has a component along the gradient of $H$, which is the direction that changes energy, and this component accumulates with a definite sign on the average. The result is linear drift: the energy error grows like $t$ rather than oscillating. For an orbit integrated over many periods, the drift eventually becomes comparable to the system’s energy scale and the trajectory is qualitatively wrong.

The accurate characterization is not that Runge-Kutta is bad and symplectic methods are good. A fourth-order method’s energy error per step can be far smaller than a second-order method’s oscillation amplitude; the difference is that one is cumulative and the other is bounded. On a short integration the fourth-order method may be better by every measure. The distinction appears when the integration length is long compared with the number of steps required to accumulate the drift, which for astronomically long simulations is the normal regime.

Partitioned Runge-Kutta methods, which use separate tableaux for positions and momenta, can be symplectic, and the corresponding condition on the two tableaux was characterized by Sanz-Serna and by Suris working independently. The condition is an algebraic relation among the tableau coefficients, not a property of the order. There is also a barrier analogous to the order barrier discussed in Butcher tableaux: no explicit Runge-Kutta method of the usual form is symplectic, and symplectic general linear methods have their own order ceiling.

Where the constraint bites the step size

Symplectic methods are not A-stable, and they should not be. Verlet’s stability region is an interval on the imaginary axis, and it requires $h\omega < 2$ for the fastest frequency $\omega$ of the system; for a molecular dynamics system the fastest frequency is a stiff bond vibration, and the step is bounded by it. The constraint is the same in kind as the parabolic step limit derived in explicit diffusion and the parabolic step limit, but the failure mode is different: exceeding it does not slow the solution’s drift, it makes the resonance with the bond frequency cause an instability that grows.

The standard responses are to constrain the fast bonds, replacing them with holonomic constraints solved by an iterative procedure, and to use an integrator that handles the constrained system while preserving the symplectic structure. A multiple-timestepping scheme applies a small step to the fast forces and a large step to the slow ones, which is a composition whose accuracy rests on the separation of time scales.

Symplecticity also does not survive a step-size change. A variable-step symplectic integrator is symplectic only in the extended phase space where the time step is an additional coordinate; the naive adjustment of the step from an error estimate destroys the conservation property. Adaptive schemes for Hamiltonian problems either adopt the extended formulation or accept the loss and use a non-symplectic method with a good error estimate, which is the trade made explicit in adaptive step size.

What this says about choosing a method

The relevant question is what property of the solution must survive a long integration. If the answer is the value of a conserved quantity, a symplectic or more generally a structure-preserving method is the default, and the step should be chosen for accuracy of the trajectory rather than for a global error tolerance. If the answer is the trajectory’s phase at a fixed final time, and the integration is short, a high-order non-symplectic method with step control is the better instrument.

There is a further class of methods that conserve the energy exactly by solving an additional constraint, the projection methods, and those are appropriate when the energy is the quantity of interest and the method needs to be general. And there is the reverse case: for a dissipative system, symplecticity is irrelevant and the structure of interest is the monotone decrease of a Lyapunov function, for which a method with a dissipative stability function is preferable. The geometric view supplies the criterion; the choice still depends on which invariant the application cares about.

One consequence deserves emphasis because it is easy to get wrong when reading energy plots. A symplectic method’s bounded oscillation and a non-symplectic method’s slow drift are both small on a short run, and a plot of energy against time over a few hundred steps may show no difference. The distinction emerges only when the integration is long enough for the drift to accumulate past the oscillation amplitude, which is a property of the method and the step, not of the problem. Comparing two integrators on a short run is therefore the wrong test; extrapolating the observed trend is the right one.