PaperlayerView exact revision on arXiv

Source-linked AI summary

WHFast: A fast and unbiased implementation of a symplectic Wisdom-Holman integrator for long term gravitational simulations

Hanno Rein, Daniel Tamayo

arXiv:1506.01084v1astro-ph.EPastro-ph.IMmath.NAnlin.CDphysics.comp-ph

TL;DR

Long-term planetary integrations need Wisdom-Holman methods that are both efficient and resistant to systematic numerical error. WHFast reimplements the method with improved Kepler solving, stable Jacobi-coordinate transformations, symplectic correctors, and variational-equation support. It is reported to be faster, to remove linear secular energy-error growth, and to follow Brouwer’s law at sufficiently small timesteps.

  • Problem

    Long-term planetary integrations require efficient symplectic methods, but floating-point and implementation errors can impair numerical stability and energy conservation.

  • Method

    WHFast reimplements the Wisdom-Holman integrator with improved Kepler solving, stable coordinate transformations, symplectic correctors, and variational-equation support.

  • Results

    WHFast is 1.5 to 5 times faster and removes linear secular energy-error growth, following Brouwer’s law over long timescales for sufficiently small timesteps.

  • Takeaways & Limitations

    WHFast provides a publicly available implementation for accurate long-term integrations and efficient calculation of chaos indicators.

  • Takeaways & Limitations

    Energy conservation depends on timestep choice: shorter steps can accumulate round-off errors, while longer steps can improve long-run energy conservation in the tested setting.

Abstract

from arXiv · show

We present WHFast, a fast and accurate implementation of a Wisdom-Holman symplectic integrator for long-term orbit integrations of planetary systems. WHFast is significantly faster and conserves energy better than all other Wisdom-Holman integrators tested. We achieve this by significantly improving the Kepler-solver and ensuring numerical stability of coordinate transformations to and from Jacobi coordinates. These refinements allow us to remove the linear secular trend in the energy error that is present in other implementations. For small enough timesteps we achieve Brouwer's law, i.e. the energy error is dominated by an unbiased random walk due to floating-point round-off errors. We implement symplectic correctors up to order eleven that significantly reduce the energy error. We also implement a symplectic tangent map for the variational equations. This allows us to efficiently calculate two widely used chaos indicators the Lyapunov characteristic number (LCN) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). WHFast is freely available as a flexible C package, as a shared library, and as an easy-to-use python module.

1 INTRODUCTION

WHFast reimplements the Wisdom-Holman approach to make long-term planetary integrations faster, more accurate, and less susceptible to systematic floating-point error. It achieves unbiased energy behavior, supports chaos analysis, and is publicly available.

  • Background: Wisdom-Holman integrators exploit a Keplerian-orbit decomposition with small planetary perturbations and are widely used for efficient celestial-mechanics integrations.Symplectic methods are favored because they typically conserve energy and angular momentum better than non-symplectic schemes.
  • Contributions: WHFast is a complete reimplementation of the Wisdom-Holman integrator focused on speed, accuracy, and floating-point stability.The authors emphasize that finite double-precision representation affects numerical stability and round-off growth.
  • Energy behavior: For sufficiently small timesteps, WHFast achieves Brouwer’s law, with energy error growing as time^1/2.The authors describe this as the behavior of an unbiased integrator whose errors are random and uncorrelated.
  • Features: WHFast improves the Kepler solver, supports variational equations and high-order symplectic correctors, and permits arbitrary units and reference frames.Correctors are implemented through order eleven, while variational equations enable chaos diagnostics.
  • Availability: The implementation is publicly available as native C99 code and an easy-to-use Python module.

2 BACKGROUND

The paper constructs a symplectic Wisdom-Holman map by expressing planetary motion in Jacobi coordinates and splitting the Hamiltonian into analytically or iteratively solvable pieces. It also introduces correctors and variational tools for improved long-term accuracy and chaos detection.

  • Jacobi Coordinates: Jacobi coordinates measure each body relative to the centre-of-mass of all lower-indexed bodies, while the central coordinate tracks the system centre of mass.For i>0, the coordinate is r′_i = r_i − R_{i−1}.
  • Hamiltonian Splitting: The Hamiltonian is split into centre-of-mass motion, Keplerian two-body terms, and interaction terms that can be evolved in suitable coordinate systems.The Keplerian term describes each body orbiting the interior mass, while interaction evolution is a kick in Cartesian coordinates.
  • Wisdom-Holman Mapping: The Wisdom-Holman map uses operator splitting with drift, kick, and drift steps, exploiting analytic or easily computed solutions for the component Hamiltonians.H0 gives straight-line motion, HInteraction changes velocities at fixed positions, and HKepler advances two-body orbits.
  • Symplectic Correctors: Symplectic correctors transform to mapping coordinates before integration and back afterward to remove high-frequency error terms.For corrector order K, the stated error is O(ϵdt^K)+O(ϵ^2dt^2), compared with O(ϵ^2dt^2) for an uncorrected second-order map.
  • Chaos Indicators: MEGNO distinguishes chaotic from quasi-periodic motion: Y(t) tends to infinity for chaos and to 2 for quasi-periodic orbits, while LCN is obtained from its time evolution.
  • Numerical Errors: Numerical error consists of a bounded integrator contribution and a finite-precision contribution caused by floating-point rounding.The bounded term is independent of time, whereas round-off accumulates through repeated operations.

3 IMPROVEMENTS

The implementation improvements target coordinate transformations and floating-point stability, where algebraically equivalent procedures can produce substantially different long-term energy behavior.

  • Floating-Point Assumptions: WHFast assumes IEEE 754 round-to-nearest, ties-to-even arithmetic and is implemented exclusively in double precision.The authors state that the approach is applicable in principle to other precisions.
  • Coordinate Choice: The interaction Hamiltonian is evolved most efficiently in Cartesian coordinates, while the Kepler Hamiltonian is easier to evolve in Jacobi coordinates.Efficient conversion between the two systems is therefore required.
  • Coordinate Transformations: Cartesian-to-Jacobi and Jacobi-to-Cartesian transformations are implemented in O(N), including systems containing test particles.
  • Coordinate Transformations: The transformations were thoroughly tested for unbiasedness because formally equivalent algorithms can differ substantially in numerical stability.
  • Coordinate Transformations: A numerically unstable coordinate-conversion variant repeatedly accesses r′_0 and subtracts in the final step, promoting error propagation.The resulting floating-point errors can be orders of magnitude larger and cause linear secular energy-error growth.

3.2 Implementation of Newton’s Method

WHFast replaces a straightforward Newton iteration with a faster convergence test and an algebraically equivalent update tailored to Kepler’s equation. The method improves accuracy and avoids indefinite floating-point cycling.

  • Baseline Method: The straightforward Newton solver stops when the relative update satisfies |dX/X| < eps, with eps typically around 10^-15.
  • Improved Method: WHFast uses an algebraically equivalent Newton update that tests convergence by comparing successive floating-point values rather than performing division and absolute-value operations.This implementation is reported to be both faster and more accurate.
  • Convergence: Tracking two previous iterates detects cycles between nearby floating-point values that otherwise might never converge.
  • Kepler Specialization: The Kepler-specific reformulation simplifies the final update using functions whose G terms depend on X and β.
  • Performance: For typical eccentricities e ≲ 0.99 and timesteps much shorter than the shortest orbital period, Newton’s method was fastest among the tested approaches.Higher-order Householder methods did not offset their greater per-iteration cost.

3.3 The Initial Guess for Kepler’s Equation: Short Timesteps

The short-timestep initial guess approximates X using the timestep and initial radius, with perturbative corrections when orbital-radius variation is sufficiently small. This approximation becomes unreliable near pericentre for highly eccentric orbits, where higher-order corrections are used.

  • Approximation regime: For low eccentricity or sufficiently short timesteps, X ≈ dt/r0 because the orbital radius varies little.The approximation uses the initial radius r0 over the timestep.
  • Approximation regime: Near pericentre, highly eccentric orbits can change radius rapidly, making X differ from dt/r0 when pericentre passages are unresolved.The radius changes by roughly a factor of 1 + e ≈ 2 toward true anomaly 90°.
  • Perturbative correction: The perturbative formulation requires three nondimensional parameters to be much smaller than unity.The variables are introduced after nondimensionalizing with X̃ = r0X/dt.
  • Perturbative correction: When ΔM ≪ (1 − e)3/2, χ, η, and ζ are small enough to extend the solution to higher order.The algorithm uses a second-order initial guess in this regime.
  • Perturbative correction: Higher-order initial guesses were slower than the chosen implementation, even for small eccentricities and timesteps.The authors attribute this to the computational efficiency of each Newton iteration.

3.4 Large Eccentricities and Timesteps

For high eccentricities or long timesteps, WHFast modifies both the Kepler solver and its initial guess. It uses stability safeguards for difficult cases and introduces a faster universal-variable guess than orbital-element conversion over a stated parameter range.

  • Solver changes: High-eccentricity and long-timestep cases require modified solvers and initial guesses.This regime is treated separately from the optimized low-eccentricity, short-timestep algorithm.
  • Solver changes: Laguerre-Conway is the most stable root-finding method tested, but bisection is sometimes required when the timestep approaches the orbital period.Long timesteps poorly sample inter-planet interactions and are more susceptible to timestep resonances.
  • Solver changes: The Newton breakout condition exits when an iteration result repeats, preventing longer floating-point cycles in this regime.The implementation stores values from each iteration to detect repetition.
  • Solver selection: A ∼1% threshold on the first Newton-iteration deviation provides a practical solver-selection balance across timestep and eccentricity values.The comparison uses a fraction of 2πβ−1/2 because it corresponds to X for a one-orbit timestep.
  • Initial guess: The usual dt/r0 guess blows up near pericentre as eccentricity increases, motivating an alternative to Danby’s orbital-element-based guess.The improved guess uses X = dt/a for a timestep of one orbital period.
  • Initial guess: Across timesteps from 0.03 to 1 orbital periods and eccentricities from 0.999 to 0.9999, the improved universal-variable guess is faster than conversion to orbital elements.The comparison covers a relevant high-eccentricity parameter range.

3.5 Implementation of c-functions

The c-function implementation prioritizes convergence, speed, and floating-point efficiency. It reduces the argument before series evaluation, stops when additional terms no longer change the result, and reconstructs the original argument using quarter-angle relations.

  • Motivation: Kepler iterations require multiple G-functions, whose c-functions therefore need optimization for both speed and accuracy.Chaos-indicator calculations require c0 through c5, while ordinary integrations require only c0 through c3.
  • Series evaluation: The algorithm repeatedly divides z by 4 until z ≤ 0.1, ensuring convergence of the c-function series while using exact power-of-two divisions.The number of reductions is tracked for later reconstruction.
  • Implementation: The implementation uses stored inverse factorials and hard-coded initial terms to avoid slower floating-point multiplication and division.The pseudocode initializes c4 and c5 with their first two terms before entering the loop.
  • Series evaluation: After the first two series terms, additional terms are added until the floating-point result stops changing.Because inverse factorials decrease quickly, this provides a fast convergence test.
  • Reconstruction: Quarter-angle relations recover c-functions for the original z after the series is evaluated at the reduced argument.The reconstructed functions include c0 through c5 as required by the calculation.
  • Implementation: The pseudocode retains previous values to detect convergence and applies the quarter-angle formula repeatedly as z is restored.The loop updates c4 and c5 before deriving c3, c2, c1, and c0.

3.6 Implementation of Gauß f and g-functions

WHFast rewrites the Gauss f and g-function calculation into an algebraically equivalent form that improves floating-point precision. The same rearrangement is applied to the variational equations.

  • Reformulation: The straightforward Gauss f and g implementation contains terms that are small relative to the leading value for timesteps below half an orbital period.The authors define new f̂ and ġ̂ functions to improve the calculation.
  • Precision: The reformulated final step sums small bracketed quantities before adding them to the larger initial value.This ordering achieves higher precision while remaining algebraically equivalent to the original equation.
  • Precision: The same numerical-stability rearrangement is used for the variational equations.This extends the precision improvement beyond the particle-state update.

3.7 A full integration in Jacobi coordinates

WHFast reduces coordinate transformations during full integrations in Jacobi coordinates, improving speed and robustness while preserving the integration workflow.

  • The Jacobi-coordinate algorithms are unbiased and fast, but WHFast aims to minimize conversions.
  • WHFast reduces each timestep to two coordinate conversions under stated assumptions about unchanged particle positions and velocities.The conversions are for positions from Jacobi to inertial coordinates and accelerations from inertial to Jacobi coordinates.
  • The integration sequence applies half-step and full-step Kepler drifts around interaction calculations, kicks, and inertial-frame updates.
  • Deferring inertial-frame velocity updates and avoiding repeated coordinate transformations reduces runtime and round-off-error buildup, making the integrator more robust.

3.8 LCN calculation

WHFast computes the Lyapunov characteristic number efficiently by updating the statistics needed for its fit at every timestep.

  • An efficient one-pass method computes the mean and covariance of Y(t) without storing all previous values.The method calculates the LCN at every timestep in O(1) and is numerically more robust.

4 NUMERICAL RESULTS

The numerical tests compare WHFast with established Wisdom-Holman and related integrators across Kepler-solver accuracy, speed, timestep scaling, and long-term energy conservation. WHFast is generally faster and less biased, with Brouwer-law behavior at sufficiently small timesteps and strong performance from symplectic correctors.

  • The experiments evaluate WHFast against publicly available integrators for speed, accuracy, and numerical stability.
  • 4.1 Two-body Kepler Solver: The Kepler-solver tests scan eccentricities from 0 to 1−10^-8 and timesteps from 0.1% to one orbital period over 100 orbits.
  • 4.1 Two-body Kepler Solver: WHFast is more accurate than standard WH across the most important parameter space, while its energy error remains random rather than biased.Biased errors produce long-term linear growth of energy error.
  • 4.1 Two-body Kepler Solver: 20%–100% speed-ups over WH are typical across the tested parameter space, even when coordinate conversions are included for fair comparison.Production runs avoid conversions at every timestep, so the reported speed and energy-conservation advantages are conservative.
  • 4.2 Short Term Energy Conservation: Second-order integrators show quadratically decreasing error as timesteps shrink from 20% to 0.1% of Jupiter’s orbital period.
  • 4.2 Short Term Energy Conservation: Below 0.1% of the shortest orbital period, floating-point round-off dominates, causing energy errors to rise and producing an integration-dependent optimum timestep.
  • 4.2 Short Term Energy Conservation: WHFast’s small-timestep error follows dt^-1/2 with increasing step count, indicating Brouwer’s law behavior.
  • 4.2 Short Term Energy Conservation: SWIFTER-TU4 also appears to follow Brouwer’s law but performs poorly at large timesteps because it lacks Keplerian splitting.

5 CONCLUSIONS

WHFast is a reimplementation of the Wisdom-Holman integrator that improves speed, accuracy, energy conservation, chaos analysis, frame and unit flexibility, and accessibility.

  • WHFast is a new implementation of a symplectic Wisdom-Holman integrator with multiple improvements over publicly available alternatives.
  • WHFast is faster by a factor of 1.5 to 5, with gains from the improved Kepler solver, combined drift steps, and fewer Jacobi-coordinate transformations.A 50% speedup comes from the improved Kepler solver.
  • The Kepler solver is more accurate and unbiased through improved Newton convergence, a Laguerre-Conway solver, accurate c and G-functions, and careful floating-point operation ordering.
  • WHFast removes the secular energy error growing linearly with integration time and follows Brouwer’s law for sufficiently small timesteps.The paper identifies unbiased Kepler solving and unbiased Jacobi-coordinate transformations as the causes.
  • Variational equations enable Lyapunov-timescale and MEGNO calculations at very low cost because they do not require solving Kepler’s equation.
  • Symplectic correctors of orders 3, 5, 7, and 11 are implemented, with three orders of magnitude improvement reported for relatively massive Solar-System-like planets.The default setting uses an 11th-order corrector, whose long-integration performance cost is negligible.
Loading 1506.01084v1…