Symplectic integrators for orbital motion
Preserve Hamiltonian geometry while assessing energy, phase and close-encounter errors separately.
On this page
The question for today
What you’ll understand
- Construct a kick–drift–kick Verlet step.
- Distinguish symplecticity from exact conservation of the original Hamiltonian.
- Identify why fixed steps and resolved orbital time scales matter.
An orbit can look plausible for a few revolutions while its numerical phase or energy drifts steadily. Symplectic integration addresses the geometry of Hamiltonian flow, which often improves long-time behavior in conservative mechanics. It does not force the original energy to be exactly constant at every step. A useful evaluation therefore checks geometric structure, energy error, orbital phase and spatial resolution together. A tiny harmonic benchmark makes these distinctions visible before a singular gravitational potential introduces close encounters and additional modeling choices.
Build the map from exact subflows
For a separable Hamiltonian, kinetic motion changes position at fixed momentum and potential motion changes momentum at fixed position. Each subflow is canonical. Composing a half potential kick, a full kinetic drift and another half kick gives velocity Verlet in momentum form. The symmetric composition has second-order global accuracy over a fixed time interval for sufficiently smooth forces. The corresponding one-step state error is third order; those two error statements refer to different comparisons.
For a gravitational test particle, the potential is minus GMm divided by radius. The kicks then use the force at the old and new positions. Rotationally invariant pair kicks can preserve angular momentum in exact arithmetic, but energy generally oscillates around a nearby value. Arbitrarily replacing the force or applying inconsistent updates to different particles can destroy properties of the intended Hamiltonian map. One must specify the actual force calculation as carefully as the time integrator.
Preserve the two-form rather than impose an energy identity
A symplectic map preserves the canonical two-form. Its derivative obeys a matrix identity involving the canonical antisymmetric matrix J. In a single position–momentum plane this implies area preservation; in higher dimensions volume preservation alone is weaker than symplecticity. Explicit Euler applied simultaneously to position and momentum fails even the planar area test for a harmonic oscillator. Small time steps reduce its error but do not make that particular map exactly symplectic.
Backward error analysis explains why a suitable fixed-step symplectic scheme approximately follows a modified Hamiltonian. Under smoothness, bounded-motion and sufficiently small-step assumptions this can yield small energy error over long intervals. The conclusion is conditional, especially near singularities, resonances or poorly resolved fast motion. It also does not imply uniformly small phase error: a slightly incorrect orbital frequency can accumulate a substantial angular displacement while the energy remains well behaved.
Resolve the orbit before trusting long-time diagnostics
An eccentric orbit has a short dynamical time near periapsis, so a step adequate near apoapsis may miss its strongest acceleration. Ordinary state-dependent step-size substitution generally sacrifices the symplectic property; special extended-phase-space constructions require their own derivation. Collisions and close encounters may need regularization or a deliberately different integrator. Softening a force avoids a singular acceleration by changing the potential, and therefore changes the physical model being integrated.
Useful diagnostics include convergence of periapsis location, phase, conserved quantities and trajectories over a stated interval. In chaotic many-body motion, pointwise trajectories can separate even when statistical observables converge, so the target observable matters. The oscillator lab supplies the exact analytic motion of a single damped or undamped oscillator. It can benchmark expected periods, but it contains no finite time step, kick–drift update or computed orbit and cannot demonstrate a numerical energy-conservation theorem.
Work it through
One Verlet step changes the true energy slightly
Use dimensionless harmonic units m = ω = 1, so H = (p² + q²)/2. Starting at q = 1, p = 0, take one kick–drift–kick step h = 0.20. Compute the new state and its original Hamiltonian.
The first kick gives . The initial energy was one half.
The drift gives . The second kick therefore uses the new position and gives .
The original energy becomes . Its change is −0.000198, despite the symplectic construction. This single step is evidence of the distinction, not a long-time error bound for every orbit.
The new state is (q, p) = (0.98, −0.198) and H = 0.499802, slightly below the initial 0.5.
The living notebook
Explore the model
Interactive lab · oscillator
An analytic period benchmark
Reset and set BOTH Damping ratio ζ = 0 and Natural frequency ω₀ = 1; predict period 2π. Reset and set ζ = 0 and ω₀ = 2; predict period π. Compare with Play or Step over several cycles. The lab evaluates analytic cosine motion from x0 = 1, v0 = 0; it does not run Verlet, vary a numerical step or calculate gravitational orbits.
- Displacement x(t)
- 1
- Velocity v(t)
- 0
- Regime
- underdamped
Unit mass, x(0) = 1 and v(0) = 0. The analytic solution includes underdamped, critically damped and overdamped cases; time is measured in seconds. Parameter changes pause playback and recompute the model at the current elapsed time. Reset restores both parameters and returns time to zero. Playback starts paused. Adjust either parameter, then choose Play or Step.
Look forAn exact reference solution supplies a benchmark, while numerical geometric properties require examining the integrator.
Pause & reason
Check your understanding
Your turn
Try it yourself
For the same dimensionless harmonic oscillator and initial state, apply simultaneous explicit Euler with h = 0.20. Find the new energy and the determinant of its position–momentum update matrix. Does it preserve planar area?
Take a moment to make a prediction or work it out on paper. Open a hint when you need one.
Give me a hint
Euler uses qnew = q + hp and pnew = p − hq, both evaluated at the old state.
Compare with the solution
Euler gives qnew = 1 and pnew = −0.20, hence Hnew = 0.52. The update matrix is , with determinant . A small phase-space area therefore grows by four percent in one step in this harmonic test. It is not area preserving or symplectic. The larger energy here is a concrete defect for this method and data; a general comparison of integrators would also examine accumulated error, force evaluations and the requested orbital observable.
Keep the idea
Take it with you
- Verlet composes canonical kicks and drifts.
- Symplectic geometry differs from exact conservation of H.
- Phase accuracy, close encounters and step selection need separate checks.
Read further
Original sources and deeper explanations to follow at your own pace.
Optional progress, saved in this browser when available. Revisit any lesson at any time.
Lesson not marked complete.