Symplectic integrator

A symplectic integrator is a numerical method for approximating the evolution of a Hamiltonian system while preserving the system’s underlying symplectic geometry. Such methods replace the continuous Hamiltonian flow by a sequence of discrete symplectic transformations. They generally do not conserve the original Hamiltonian exactly, but their preservation of phase-space structure produces characteristic long-time behavior that differs from that of generic numerical integrators.

Symplectic integrators are principally used when the geometric organization of trajectories is relevant over many time steps. Their applications include long-term celestial mechanics, charged-particle dynamics, molecular dynamics, and accelerator modeling. The defining property concerns preservation of a differential two-form rather than local accuracy alone, so a symplectic method may have the same formal order as a nonsymplectic method while exhibiting substantially different accumulated behavior.

Geometric formulation

For canonical coordinates (q,p\in\mathbb{R}^d), a Hamiltonian (H(q,p)) determines the equations

[ \dot q=\frac{\partial H}{\partial p}, \qquad \dot p=-\frac{\partial H}{\partial q}. ]

Writing (z=(q,p)), these equations take the matrix form

[ \dot z=J^{-1}\nabla H(z), \qquad J= \begin{pmatrix} 0&I\ -I&0 \end{pmatrix}. ]

The associated exact flow (\varphi_t) preserves the canonical symplectic form

[ \omega=\sum_{i=1}^{d} dq_i\wedge dp_i. ]

Equivalently, its Jacobian satisfies

[ (D\varphi_t(z))^{\mathsf T}J,D\varphi_t(z)=J. ]

A one-step numerical map (\Phi_h), with step size (h), is symplectic when it obeys the corresponding discrete condition

[ (D\Phi_h(z))^{\mathsf T}J,D\Phi_h(z)=J. ]

Symplecticity implies preservation of oriented phase-space volume, but it is stronger than the volume-preservation property established by Liouville’s theorem. A volume-preserving map need not preserve the symplectic two-form or the canonical relations among position and momentum variables.

The condition is invariant under canonical changes of coordinates. Consequently, symplecticity describes the discrete map itself rather than a particular coordinate representation of its update equations.

Splitting methods

Many practical symplectic integrators arise from a separable Hamiltonian

[ H(q,p)=T(p)+V(q), ]

where (T) represents the momentum-dependent contribution and (V) represents the coordinate-dependent contribution. Each partial Hamiltonian has an exactly computable flow under common assumptions. The flow of (T) changes (q) while holding (p) fixed, whereas the flow of (V) changes (p) while holding (q) fixed.

The first-order symplectic Euler method can be written as

[ p_{n+1}=p_n-h\nabla V(q_n), ]

[ q_{n+1}=q_n+h\nabla T(p_{n+1}). ]

The ordering can be reversed, producing its adjoint form. Each ordering is symplectic, although neither is generally invariant under time reversal.

A symmetric composition of the two partial flows gives the second-order Störmer–Verlet method. In its kick–drift–kick representation,

[ p_{n+\frac12}

p_n-\frac h2\nabla V(q_n), ]

[ q_{n+1}

q_n+h\nabla T\left(p_{n+\frac12}\right), ]

[ p_{n+1}

p_{n+\frac12}-\frac h2\nabla V(q_{n+1}). ]

Each substep is an exact Hamiltonian flow and therefore a symplectic transformation. Their composition remains symplectic. The symmetric arrangement also makes the full step self-adjoint, which accounts for its even order and time-reversal symmetry.

Higher-order splitting formulas use compositions with coefficients chosen to cancel lower-order terms in the Baker–Campbell–Hausdorff formula. Real-coefficient compositions above second order generally contain negative substeps when they are built solely from the two original flows. This algebraic restriction affects problems whose evolution cannot be extended meaningfully to negative time.

Runge–Kutta structure

Symplectic methods also occur within the theory of Runge–Kutta methods. A Runge–Kutta scheme with coefficients (a_{ij}), weights (b_i), and (s) stages is symplectic for general canonical Hamiltonian systems when

[ b_i a_{ij}+b_j a_{ji}=b_i b_j ]

for every pair of stage indices (i,j).

The implicit Runge–Kutta methods based on Gauss–Legendre quadrature satisfy this condition. An (s)-stage Gauss method has order (2s), preserves every quadratic invariant respected by the differential equation, and defines a symplectic discrete map. Its implicit stage equations distinguish it computationally from splitting methods based on exactly solvable partial Hamiltonians.

Ordinary explicit Runge–Kutta methods do not satisfy the general symplecticity condition except in degenerate cases. Partitioned Runge–Kutta methods avoid this restriction by applying related coefficient arrays to different groups of variables. This framework includes symplectic Euler and several formulations of Störmer–Verlet.

Modified Hamiltonians and long-time behavior

The long-time properties of symplectic integration are described by backward error analysis. For a symplectic method of order (r), the numerical step can formally be interpreted as the exact time-(h) flow of a modified Hamiltonian

[ \widetilde H

H+h^rH_r+h^{r+1}H_{r+1}+\cdots. ]

The displayed series is normally asymptotic rather than convergent. Truncation at an order depending on the step size nevertheless provides a modified differential equation whose flow approximates the numerical map with a small defect over extended intervals.

Because the discrete trajectory follows a nearby Hamiltonian system, its value of the original Hamiltonian usually undergoes bounded oscillations rather than systematic monotone drift. Exact conservation of (H) is not implied. The method instead nearly conserves (\widetilde H), while the difference between (H) and (\widetilde H) remains controlled by powers of the step size under the assumptions of backward error analysis.

For analytic systems with bounded numerical trajectories and sufficiently small fixed step size, this description can remain valid over intervals exponentially long in the reciprocal step size. Resonances between the step size and intrinsic system frequencies can alter that behavior. Variable step sizes can also destroy symplecticity when introduced directly, because the resulting state-dependent time advance is generally not a canonical transformation on the original phase space.

Symplecticity does not by itself preserve all first integrals. Momentum maps associated with continuous symmetries are retained when the discrete method respects the same symmetry, a relation expressed by the discrete form of Noether’s theorem. Dissipative evolution is not naturally represented by a symplectic map on the original phase space because such evolution contracts volume and lacks Hamiltonian flow structure.

Historical development

Early components of symplectic integration appeared before the terminology became standardized. The leapfrog formulas used in mechanics and astronomy were instances of symplectic splitting, although they were commonly classified through difference equations rather than geometric preservation. René De Vogelaere analyzed related second-order methods in the 1950s within numerical celestial mechanics.

Ronald Ruth formulated explicitly symplectic algorithms for accelerator Hamiltonians in 1983, including a third-order composition method. Feng Kang subsequently developed the theory of symplectic difference schemes and connected discrete numerical evolution with canonical transformations. These developments established symplecticity as a design condition rather than an incidental property of selected mechanical algorithms.

In 1985, You Watanabe derived a generating-function representation of composed drift and kick maps for separable Hamiltonians. The formulation showed directly that the discrete update preserved the canonical two-form and placed the corresponding splitting schemes within the theory of canonical transformations.

The term “symplectic integrator” became established through the work of Jesús María Sanz-Serna during the late 1980s. His analysis identified algebraic symplecticity conditions for Runge–Kutta discretizations and clarified the distinction between exact conservation of energy and preservation of Hamiltonian geometry. Subsequent backward-error results connected these algebraic conditions with the observed long-time stability of numerical invariants.

Scope and limitations

A symplectic integrator controls geometric distortion rather than eliminating numerical error. Its local truncation error still accumulates, and the computed phase of an oscillatory solution can gradually separate from the exact phase. Symplecticity therefore does not imply uniformly accurate trajectories for arbitrarily long times.

The method’s structure also depends on the Hamiltonian representation. Noncanonical Hamiltonian systems require preservation of the appropriate nonconstant symplectic form or Poisson structure. Constrained systems require compatible treatment of the constraint manifold; methods such as SHAKE and RATTLE combine constraint enforcement with symplectic evolution under suitable regularity conditions.

Adaptive integration requires an extended formulation when exact symplecticity is to be retained. One such formulation promotes time to a coordinate and introduces a conjugate momentum, producing an enlarged Hamiltonian system in which variable physical time steps arise from a fixed-step canonical map. The symplectic form then belongs to the extended phase space rather than the original position–momentum space alone.