A class of symplectic integrators with adaptive timestep for separable Hamiltonian systems

TL;DR

This paper develops adaptive-timestep symplectic integrators for separable Hamiltonians by extending the phase space with a fictitious time, enabling explicit, reversible leapfrog schemes in which the timestep depends solely on the potential energy. A key result is that, for Keplerian potentials, the γ=1 scheme can integrate Kepler orbits exactly up to a timing error, while γ=3/2 offers strong accuracy improvements over fixed-timestep methods; the Stark problem and perturbed Kepler cases further illustrate robustness and error behavior. The approach avoids common pitfalls of adaptive schemes breaking symplectic structure and shows competitive performance relative to other regularized or Wisdom-Holman-like methods, especially for test-particle or near-Keplerian dynamics. Overall, the work provides practical, structure-preserving adaptive integrators with solid long-term error control and applicability to a range of celestial-mechanical problems.

Abstract

Symplectic integration algorithms are well-suited for long-term integrations of Hamiltonian systems because they preserve the geometric structure of the Hamiltonian flow. However, this desirable property is generally lost when adaptive timestep control is added to a symplectic integrator. We describe an adaptive-timestep symplectic integrator that can be used if the Hamiltonian is the sum of kinetic and potential energy components and the required timestep depends only on the potential energy (e.g. test-particle integrations in fixed potentials). In particular, we describe an explicit, reversible, symplectic, leapfrog integrator for a test particle in a near-Keplerian potential; this integrator has timestep proportional to distance from the attracting mass and has the remarkable property of integrating orbits in an inverse-square force field with only "along-track" errors; i.e. the phase-space shape of a Keplerian orbit is reproduced exactly, but the orbital period is in error by O(1/N^2), where N is the number of steps per period.

Figures (3)

• Figure 1: Maximum fractional energy error over $2\times10^4$ orbital periods, as a function of number of steps per orbit for the Keplerian two-body problem. The curves correspond to eccentricities $e=0.9, 0.99, 0.999, 0.9999, 0.99999, 0.999999$. The integrator is standard DKD (drift-kick-drift) leapfrog with timestep $\propto r^{3/2}$, following equations (\ref{['eq:gpow']}) and (\ref{['eq:ffff']}) with $\gamma={3\over2}$. The orbits are started at pericenter. The dashed lines show analytic estimates of the energy error from equations (\ref{['eq:perierr']}) and (\ref{['eq:nstep']}). The analogous errors for DKD leapfrog with fixed timestep and $e=0.9$ and $0.99$ are shown as open circles; for larger eccentricities the fixed-timestep errors are off-scale.
• Figure 2: Fractional energy error as a function of distance from the attracting mass, for a numerical integration of the Stark problem using DKD leapfrog with $\gamma=1$, which integrates Keplerian orbits with zero energy error. The integration lasts for 1000 Keplerian periods of the initial orbit and the error is plotted every 100 timesteps. The integration parameters are $\mu=1$, $\epsilon=0.1$, the initial eccentricity is $e=0.9$, the Stark vector is $45^\circ$ from the initial line of apsides, and its magnitude is $S=\eta E^2/\mu$ where $\eta=4\times10^{-3}$. For $r\ll 1$ the points lie on a straight line, consistent with the prediction of equation (\ref{['eq:errstark']}) marked by solid line segments. The open circles show the much smaller energy errors when the initial conditions are corrected using equation (\ref{['eq:correct']}).
• Figure 3: Fractional energy error for the Stark problem, as a function of number of force evaluations per orbit. The initial eccentricity is $e=0.9$, the Stark vector is $45^\circ$ from the initial line of apsides, and the orbit starts at apocenter and is followed for $10^4$ periods. The magnitude of the Stark vector is $S=\eta E^2/\mu$ where $\eta=0.001$ (solid lines), 0.005 (dashed lines), and 0.02 (dash-dot lines). The lower curves (solid circles) represent the average of the absolute value of the energy error, and the upper curves the maximum error. The integrator is DKD leapfrog with $\gamma=1$ and initial conditions corrected using equation (\ref{['eq:correct']}); the single solid curve with open circles represents the average error for $\eta=0.001$ if no corrector is applied.