Source-linked AI summary
IAS15: A fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits
Hanno Rein, David S. Spiegel
TL;DR
Long-term gravitational integrations must balance accuracy, timestep selection, and applicability to non-Hamiltonian forces. IAS15 addresses these needs with a 15th-order Gauß-Radau-based integrator, compensated arithmetic, and physically motivated adaptive timestepping. It maintains machine-precision behavior over billion-timescale integrations and outperforms tested integrators across broad accuracy and runtime comparisons.
Problem
Symplectic integrators face difficulties with adaptive timesteps and non-conservative velocity-dependent forces, limiting their suitability for some astrophysical systems.
Method
IAS15 combines a 15th-order Gauß-Radau scheme with compensated summation and dimensionless, physically motivated adaptive timestepping.
Results
Across diverse astrophysical tests, IAS15 is better by the evaluated metrics or within a factor of ∼2-3 of fine-tuned alternatives, while its energy error follows Brouwer’s law over at least 10^9 dynamical timescales.
Takeaways & Limitations
IAS15 provides a broadly applicable open-source integrator for conservative and non-conservative gravitational dynamics, including systems with velocity-dependent forces.
Takeaways & Limitations
Floating-point limitations can make the timestep error estimate unreliable for poorly chosen coordinate frames, motivating a potentially accelerated frame.
Abstract
from arXiv · showhide
We present IAS15, a 15th-order integrator to simulate gravitational dynamics. The integrator is based on a Gauß-Radau quadrature and can handle conservative as well as non-conservative forces. We develop a step-size control that can automatically choose an optimal timestep. The algorithm can handle close encounters and high-eccentricity orbits. The systematic errors are kept well below machine precision and long-term orbit integrations over $10^9$ orbits show that IAS15 is optimal in the sense that it follows Brouwer's law, i.e. the energy error behaves like a random walk. Our tests show that IAS15 is superior to a mixed-variable symplectic integrator (MVS) and other popular integrators, including high-order ones, in both speed and accuracy. In fact, IAS15 preserves the symplecticity of Hamiltonian systems better than the commonly-used nominally symplectic integrators to which we compared it. We provide an open-source implementation of IAS15. The package comes with several easy-to-extend examples involving resonant planetary systems, Kozai-Lidov cycles, close encounters, radiation pressure, quadrupole moment, and generic damping functions that can, among other things, be used to simulate planet-disc interactions. Other non-conservative forces can be added easily.
1 INTRODUCTION
IAS15 addresses limitations of symplectic orbit integrators by combining high-order accuracy with adaptive timestepping and support for velocity-dependent forces. It achieves machine-precision long-term integrations while following Brouwer’s law and outperforming tested alternatives across accuracy and speed.
- Motivation: Symplectic integrators complicate adaptive timestepping, require a designated central object, and are unsuitable for non-conservative velocity-dependent forces.These limitations affect systems including stellar binaries, planetary rings, and dust subject to radiation forces.
- Contribution: IAS15 is a high-order approach that does not require symplecticity and supports high-precision calculations with velocity-dependent forces.Its energy conservation is reported below double-precision machine precision.
- Results: IAS15’s energy error follows Brouwer’s law, whereas several popular integrators show linear energy growth in long simulations.Brouwer’s law describes energy error growing like a random walk and is presented as the attainable limit without extended precision.
- Results: 10^9 orbital timescales are integrated at machine precision using 100 timesteps per orbit, an order of magnitude fewer than NBI requires.NBI is identified as another integrator that achieves Brouwer’s law.
- Contribution: IAS15 suppresses systematic error below machine precision, introduces physically motivated adaptive timestepping, and controls round-off errors at machine precision.The timestep control uses no arbitrary scales, while compensated summation helps maintain symmetric round-off errors.
- Implementation: IAS15 is implemented in the freely available REBOUND package, which supports diverse astrophysical problems and provides a Python wrapper.The package is modular and can be extended to applications such as planetary rings, granular flows, and planet migration.
2 IAS15 INTEGRATOR
IAS15 builds on Gauß-Radau integration by expanding acceleration in a high-order series, evaluating forces at Gauß-Radau substeps, and solving the resulting implicit system with predictor-corrector iterations. Its convergence and timestep procedures are designed to reach machine precision while remaining adaptive and scale-independent.
- 2 IAS15 INTEGRATOR: IAS15 builds on Everhart’s 15th-order modified Runge-Kutta method and improves its accuracy and step-size control.The implementation is added to REBOUND and supplied with a Python wrapper.
- 2 IAS15 INTEGRATOR: The solved equation allows acceleration to depend on position, velocity, and time, including arbitrary velocity-dependent non-conservative forces.Such systems need not be Hamiltonian.
- 2.1 Algorithm: The acceleration is represented by a truncated high-order series whose dimensionless timestep variable rescales coefficients with acceleration units.The coefficients are rewritten in terms of the normalized variable h = t/dt.
- 2.1 Algorithm: Gauß-Radau substeps improve integral accuracy by including the timestep’s starting point and using specialized quadrature spacing.Force evaluations at the substeps determine intermediate coefficients used to recover positions and velocities.
- 2.1 Algorithm: A predictor-corrector loop iteratively updates positions, velocities, and force-dependent coefficients until machine-precision convergence.After initialization, only about two iterations are needed because previous coefficients provide predictions for subsequent timesteps.
- 2.2 Warmup procedure and convergence of the predictor-corrector loop: The convergence test terminates when the change in the highest coefficient is insignificant relative to acceleration or begins oscillating.The threshold ϵ_δb = 10^-16 ensures convergence of the implicit part rather than setting the scheme’s final order or accuracy.
- 2.3 Stepsize control: The adaptive timestep is controlled by a dimensionless smoothness parameter and remains independent of the problem’s physical scale under rescaling.This avoids the dimensional-scale dependence and code-unit pitfall described for Everhart’s implementation.
- 2.3 Stepsize control: Because IAS15 is 15th order, reducing the timestep by α reduces the error by α^16, but rapid timescale changes within one timestep can still defeat the method.The authors state they have not found a physical scenario exhibiting that failure mode.
3 ERROR ESTIMATES
This section decomposes integration error into numerical-floor, random, biased, and scheme contributions, then explains how IAS15 controls scheme error while accounting for floating-point effects. It also identifies coordinate-system and precision limits that constrain error estimation.
- Floating-point errors: IEEE 754 double precision limits relative-error accuracy to approximately 2·10^-16, defining the floor contribution Efloor.Floating-point operations also accumulate random and potentially biased errors over time.
- Floating-point errors: Compensated summation tracks lost precision when small increments are added to large values, reducing IAS15’s random error by one to two orders of magnitude.The implementation applies it to position and velocity updates during and at the ends of timesteps.
- Floating-point errors: Adaptive timestepping makes rounding errors more random by preventing identical step sizes from repeatedly producing similarly biased rounding.Numerical tests report unbiased behavior over at least 10^11 timesteps, equivalent to one billion orbits.
- Error components: The total error comprises floor, random, biased, and scheme contributions, whose largest term determines a scheme’s effective accuracy.At least one contribution may be constant, while others can grow as t^1/2, linearly, or with higher powers of time.
- IAS15 error control: IAS15 reaches machine precision when its scheme error becomes negligible, leaving floor and random errors, and possibly bias, to dominate.Its adaptive scheme keeps the scheme error negligible within machine precision for at least 10^9 dynamical timescales.
- IAS15 error control: Choosing ϵ_b ≈ 0.028 yields a one-timestep relative error at machine precision and a timestep ratio ϵ_dt ≈ 0.3.The timestep remains a fraction of the smallest characteristic dynamical time.
- Error-estimation limits: Error estimation can fail when floating-point precision degrades relative positions, especially for highly eccentric, offset particles requiring small timesteps.IAS15 excludes particles moving very little during a timestep from the relevant maxima to prevent runaway timestep reduction.
- Error-estimation limits: In precision-limited simulations, changing to a suitable accelerated frame may improve position and acceleration estimates, but the details depend on the problem.The paper identifies this as a limitation shared by all integrators rather than IAS15 alone.
4 TESTS
The tests compare IAS15 with established integrators across short-term Solar System dynamics, close encounters, Kozai-Lidov cycles, and billion-orbit integrations. IAS15 achieves high accuracy, adaptive encounter handling, and Brouwer-law energy behavior, while competing methods show narrower validity or larger errors.
- Short term simulations of the outer Solar System: IAS15 preserves phase better than WH for every timestep and better than MVS for any reasonable timestep.The phase-error test integrates the outer Solar System forward and backward for 50 Jupiter orbits.
- Jupiter-grazing comets: Adaptive timestepping detects every tested Jupiter-comet encounter, while preserving the Jacobi constant to 10^-14 or better at closest approach of approximately 10 Jupiter radii.Encounters within 1 to 10 Jupiter radii remain captured, but Jacobi-constant conservation worsens because of floating-point precision and coordinate choice.
- Kozai-Lidov cycles: In Kozai-Lidov cycles, IAS15 preserves energy to 10^-12 and angular momentum to 10^-15 across simulations, including extreme eccentricities up to e approximately 1 - 10^-10.Other tested integrators fail or are substantially less accurate in the equal-mass-binary configuration, while the canonical precision parameter is conservative for the tested single cycle.
- Long-term simulations: Over billion-orbit integrations, IAS15 follows Brouwer’s law: energy error remains near machine precision initially and later grows like √t rather than linearly.IAS15’s errors are smaller than those of MVS, RADAU, and BS throughout the comparison; competing behavior includes linear growth and strong precision-parameter dependence.
5 CONCLUSIONS
IAS15 combines high-order integration, precision safeguards, and adaptive timestepping to deliver robust long-term N-body accuracy across diverse astrophysical problems. Its energy error follows Brouwer’s law over 10^9 dynamical timescales, while broad tests show strong accuracy and runtime performance.
- Usability: The implementation provides robust default settings, an adaptive timestep algorithm, and modifiable examples for several astrophysical force models and dynamical systems.IAS15 is freely available within REBOUND, with examples including Kozai-Lidov cycles, radiation forces, quadrupole moments, and migration forces.
- Accuracy: IAS15’s compensated summation and optimized numerical constants keep bias and scheme errors subdominant over a billion dynamical timescales.These measures reduce loss of precision and support long integrations limited by floating-point precision.
- Accuracy: 10^9 dynamical timescales: IAS15’s outer-Solar-System energy error follows Brouwer’s law, reaching the limit set by floating-point precision.The paper describes this behavior as optimal without extended precision.
- Comparisons: IAS15 outperforms tested integrators across energy, phase, angular-momentum, Jacobi-constant, and runtime metrics in astrophysically interesting problems.The comparison includes symplectic, mixed-variable symplectic, Bulirsch-Stoer, and other Gauß-Radau integrators.
- Scope: IAS15 handles Hamiltonian and non-Hamiltonian systems without requiring a specific coordinate frame, including problems with velocity-dependent or non-conservative forces.This broadens its applicability beyond heliocentric planetary integrations and standard conservative dynamics.
- Scope: IAS15 is not symplectic, and fine-tuning can provide small speedups while making the integrator less robust in other situations.The paper presents this as a practical trade-off when optimizing for a specific astrophysical problem.
APPENDIX A: SIMPLE DERIVATION OF POYNTING-ROBERTSON DRAG
The appendix derives radiation pressure and Poynting-Robertson drag for a dust particle orbiting a star. It decomposes the total force into radial and azimuthal components and accounts for motion-dependent photon effects.
- Setup: The derivation assumes a dust particle absorbs stellar photons and reemits them randomly while orbiting the star.The particle’s velocity and radial position determine the radiation-force description.
- Radiation pressure: β measures stationary-particle radiation pressure relative to stellar gravity and depends on luminosity, grain density, stellar mass, and grain diameter.Because radiation pressure scales as r^-2, it effectively reduces the star’s mass by a factor involving β.
- Radial motion: Radial motion changes the radiation force through Doppler-boosted photon energy and a corresponding change in photon encounter rate.For |ṙ| ≪ c, the combined factor is approximately 1 − 2ṙ/c.
- Force components: The azimuthal force component depends on the particle’s azimuthal velocity and is identified with Poynting-Robertson drag.Radiation pressure denotes the radial component, while Poynting-Robertson drag denotes the azimuthal component.
- Force components: The total radiation force is obtained by summing its radial and azimuthal components.This decomposition provides the force used for the non-Hamiltonian example in REBOUND.