M084 · Geometry of numerical integrators

Hamiltonian Phase Space & Integrators

When a computer simulates an orbit, can a small energy error hide a wrong trajectory? It can. A symplectic integrator, one built to respect the structure of Hamilton's equations, keeps the energy error small and bounded for a very long time. The computed body can still drift away from the true trajectory, running ahead of or behind it along the orbit. Energy behavior, accuracy, symplectic structure and time reversibility are separate properties; check each one.

Run Euler, RK4, symplectic Euler and Verlet side by side. Compare how closely each follows the true phase or orbit over long times, how its energy drifts, and whether its step map keeps the local canonical geometry and runs backward correctly. Then halve the step and take on a Kepler periapsis where the motion changes fast.

Reviewed scaled time $\text{Pending}$

Actual constant step $\text{Pending}$

Total steps per method $\text{Pending}$

Selected method $\text{Pending}$

Current relative energy error $\text{Pending}$

Maximum absolute energy error $\text{Pending}$

Current state error $\text{Pending}$

Last retained state error $\text{Pending}$

Local full-map determinant $\text{Pending}$

Local symplectic defect $\text{Pending}$

One-step reversal error $\text{Pending}$

Reference construction $\text{Pending}$

Reference refinement disagreement $\text{Pending}$

Selected method record $\text{Pending}$

Retained record end time $\text{Pending}$

Canonical dimension $\text{Pending}$

Change the setup and compare trajectories

When a computer simulates an orbit, can a small energy error hide a wrong trajectory?

System: oscillator / pendulum / Kepler

\mathcal M

0.00\,

Initial coordinate

q_0

1.00\,

Initial momentum

p_0

0.00\,

Initial Kepler eccentricity

e

0.40\,

Requested step fraction

h/T_0

0.02\,

Nominal cycles

T/T_0

20.00\,

Method: Euler / RK4 / symplectic Euler / Verlet

\mathcal I

3.00\,

Canonical patch radius

\rho

0.03\,

Reviewed time

t/T

0.00\,\%

Drag the orange initial-state ring; the matching coordinate and momentum sliders work from the keyboard. In Kepler mode the ring sets the eccentricity, and the oscillator controls have no effect. Every record is computed in advance, so playback only replays stored samples. The lower-right view plots each patch relative to its own center. The Kepler patch there is only a projection; the determinant and the symplectic error use the full four-dimensional canonical map. A record that diverges or nears a singularity stops and keeps its label.

The three systems are conservative, dimensionless and of unit mass: a harmonic oscillator, an ideal pendulum and a point mass on a planar Kepler orbit. All four methods use the same constant step. Nominal cycles count the harmonic or Kepler reference period, not the period of the pendulum at finite amplitude. The oscillator and Kepler references are analytic solutions; the pendulum reference uses RK4 at two finer steps. The difference between those two refinements is an estimate, not a rigorous error bound. There is no adaptive step size, damping or physical collision model.

Start with the essentials

Focus question: When a computer simulates an orbit, can a small energy error hide a wrong trajectory?
One-sentence intuition: Accuracy, energy behavior, symplecticity and reversibility are related, but each one measures something different.

Core mathematical model

Hamilton's equations for the three systems

\dot{\mathbf q}=\mathbf p,\qquad\dot{\mathbf p}=-\nabla V,\qquad H=\frac{|\mathbf p|^2}{2}+V

Unit mass and scaled units let the three conservative systems be compared directly. The harmonic frequency and the Kepler gravity parameter are both one.

A kick followed by a drift

\mathbf p_{j+1}=\mathbf p_j-h\nabla V(\mathbf q_j),\qquad\mathbf q_{j+1}=\mathbf q_j+h\mathbf p_{j+1}

Symplectic Euler is first order. It differs from explicit Euler because the drift uses the momentum that was just updated.

Symmetric Verlet steps

\mathbf p_{j+1/2}=\mathbf p_j-\frac h2\nabla V(\mathbf q_j),\quad\mathbf q_{j+1}=\mathbf q_j+h\mathbf p_{j+1/2},\quad\mathbf p_{j+1}=\mathbf p_{j+1/2}-\frac h2\nabla V(\mathbf q_{j+1})

Splitting the kick into two halves gives a second-order, time-reversible map. That does not conserve the true Hamiltonian exactly, and it does not make a periapsis taken in too few steps accurate.

Checking that the step map is symplectic

D\Phi_h^{\mathsf T}J D\Phi_h=J,\qquad J=\begin{pmatrix}0&I\\-I&0\end{pmatrix}

The readout gives the largest deviation of any matrix entry, estimated with central differences. In four dimensions, a determinant of one is a weaker test on its own. The Kepler patch is a two-dimensional projection and never measures the full phase-space volume.

Energy error and state error are measured separately

\epsilon_H=\frac{H_j-H_0}{\max(|H_0|,10^{-6})},\qquad\epsilon_z=\|\mathbf z_j-\mathbf z_{\rm ref}(t_j)\|_2

The maximum energy error uses every step, but the plots keep at most about 600 samples per method. A steady energy curve can hide a phase error that keeps growing. A record that stops is not extended with a prediction.

How the logarithmic axes work

y_H=\operatorname{sgn}(\epsilon_H)\log_{10}\left(1+\frac{|\epsilon_H|}{10^{-8}}\right),\qquad y_z=\log_{10}\left(1+\frac{\epsilon_z}{10^{-8}}\right)

One shared logarithmic scale shows tiny and large errors together and keeps the sign of the energy drift. Nominal cycles use a reference period of two pi in scaled time. The pendulum reference compares runs with 32 and 64 RK4 substeps in each coarse step.

Common difficulties

Good energy does not mean a good orbit

Typical misconceptionIf the energy error is small, the orbit must be right.

Better mental modelPhase error can build up while the energy stays nearly constant. Compare with the reference state at the same times.

Symplectic does not mean the energy is exact

Typical misconceptionA symplectic method conserves the exact Hamiltonian at every step.

Better mental modelWhen the fixed step resolves the motion, the energy error usually oscillates within bounds, but the true energy is still not conserved exactly.

Periapsis still needs small steps

Typical misconceptionPreserving the geometry makes up for a large time step.

Better mental modelThe fast pass close to the centre still needs a short step. Halve the step and compare the orbit error, not just the determinant.

A projected area is not phase-space volume

Typical misconceptionIf the projected Kepler patch keeps its area, the method is symplectic.

Better mental modelKepler phase space has four dimensions. The full matrix test gives the relevant local evidence.

Run the experiment

01
Compare all methods
Review a full harmonic record and switch the selected method.
What to observe: Every color uses the same step, yet the energy behaves differently: it grows, drifts or oscillates within bounds.
02
Halve the step
Use the long RK4 record, then halve the fixed step.
What to observe: Both the state error and the energy drift change with the step size. A high local order is not the only thing that matters.
03
Leave the small-angle limit
Compare the large-amplitude libration and rotating pendulum presets.
What to observe: The reference here is a refined numerical solution, not a harmonic approximation.
04
Take on a Kepler periapsis
Compare a resolved and an unresolved periapsis, check the full-map defect, then halve the step.
What to observe: A symplectic map can keep its geometry and still give a poor trajectory.