Heyoka Taylor Adaptive Integrator

• Heyoka Taylor Adaptive Integrator is a high-order, adaptive Taylor series method that uses automatic differentiation to solve ODEs with long-term precision and robustness.
• It employs adaptive step size and order selection to control local error tolerances while delivering dense output, efficient event detection, and SIMD-enabled performance.
• Benchmark comparisons demonstrate that Heyoka outperforms traditional integrators in speed and accuracy for complex astrodynamics and celestial mechanics applications.

The Heyoka Taylor Adaptive Integrator is a general-purpose, high-order adaptive Taylor series method for the numerical solution of ordinary differential equations (ODEs), designed for high-precision and long-term integrations in astrodynamics and celestial mechanics. Heyoka achieves both speed and accuracy superior to many domain-specialized integrators, offering robust features such as automatic differentiation–based recurrence, dense output, rigorous event detection, SIMD parallelism, and high-accuracy operation with support for extended and arbitrary-precision arithmetic. The implementation is open source, available in both C++ and Python interfaces (Biscani et al., 2021, Biscani et al., 2022).

1. Mathematical Formulation

Heyoka integrates initial-value ODE problems of the form $x'(t) = F(t,x(t)),\,x(t_0)=x_0$ by computing the truncated Taylor expansion of $x(t)$ about $t_0$ to order $p$:

$x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$

where $x^{[k]}(t) = \frac{1}{k!} \frac{d^k x}{dt^k}$ are the normalized derivatives.

Derivatives up to order $p$ are generated using forward automatic differentiation (AD) applied to the abstract syntax tree (AST) of $F$, exploiting recurrence relations (sum, product, quotient rules). The total computational complexity for all normalized derivatives is $O(p^2 N_\text{expr})$, where $N_\text{expr}$ is the number of elementary sub-expressions in $x(t)$0 (Biscani et al., 2021).

2. Adaptive Step Size and Error Control

Heyoka employs an adaptive step-size and order selection strategy to enforce a user-supplied local error tolerance $x(t)$1. The step is accepted if the local Taylor remainder satisfies:

$x(t)$2

If the step fails, $x(t)$3 is reduced according to

$x(t)$4

where $x(t)$5. The Taylor order $x(t)$6 is selected at integrator build time via $x(t)$7 and can be tuned for optimal FLOPs per unit time (Biscani et al., 2021, Biscani et al., 2022). The method is a true one-step integrator, requiring no history or warm-up.

3. Algorithmic Architecture and Features

Automatic Differentiation and Just-in-Time Compilation

The dynamical system is encoded as a C++ AST via operator overloading. During construction, common subexpressions are eliminated and LLVM IR is generated for the full AD recurrences and Taylor expansion, yielding a just-in-time (JIT) compiled stepping kernel tailored for the user’s system.

SIMD Batch Mode and Parallelism

Heyoka supports SIMD "batch mode", in which each scalar variable is replaced by a small SIMD vector (e.g., size 2 for SSE2, 4 for AVX), enabling simultaneous integration of multiple trajectories and up to 3.3x throughput improvement. Vectorized transcendental functions are handled via SLEEF (Biscani et al., 2021).

Dense Output and High-Accuracy Mode

Taylor expansion naturally provides dense output for any state or event function, accurate to $x(t)$8 at no extra runtime cost. High-accuracy mode combines compensated summation, tight tolerances, and extended-precision types to maintain invariants such as Brouwer’s law for energy drift (RMS error $x(t)$9) over $t_0$0 dynamical timescales (Biscani et al., 2021).

4. Event Detection and Hybrid Dynamics

Heyoka implements robust event detection for ODEs with event functions $t_0$1 by applying its Taylor machinery to $t_0$2 in parallel with $t_0$3. The event polynomial $t_0$4 for $t_0$5 is generated at each step and roots are located using Descartes’ rule of signs and the Collins–Akritas real-root isolation algorithm, guaranteeing complete and ordered detection of all zero crossings in the interval (Biscani et al., 2022).

Events may be flagged as terminal or non-terminal. When multiple events occur in a step, heyoka isolates all triggers, selects the earliest, and propagates to the corresponding state, applying "cool-down" intervals to avoid discontinuity sticking, e.g.,

$t_0$6

All event routines (including vectorized and in-place root refinement) are optimized for minimal additional overhead.

5. Implementation and API

C++ and Python Front Ends

The C++ interface uses symbolic variable declaration and function composition (through operator overloading) to define ODE systems; integrators are constructed and invoked via a single make_integrator<>() call. The Python front end (via pybind11) offers the same API semantics, enabling seamless workflow across both environments (Biscani et al., 2021).

Memory and Compute Optimizations

All Taylor coefficients are preallocated in flat memory arrays for state and events, and AD is performed in a single pass emulating a sparse forward-mode. No heap allocation is incurred at each step. Multiple kernel variants are supported: "default" (full loop unrolling for speed), "compact" (looped, reduced compile time at some speed cost).

6. Performance Benchmarks and Comparative Results

Two-body and N-body Problems

• For low-eccentricity Kepler orbits (double precision, $t_0$7, $t_0$8), heyoka achieves $t_0$916 steps/orbit and $p$0 energy error, whereas IAS15 requires $p$1100 steps/orbit for similar error (Biscani et al., 2021).
• In Sun+N-planet models, heyoka outperforms IAS15 by $p$235$p$3 for $p$4; the crossover where force evaluations dominate occurs at $p$5.

Long-term Solar System and Chaotic Dynamics

• 6-body (Sun+planets) integrations over $p$6 years show RMS energy error consistent with Brouwer’s law, with heyoka requiring $p$7s per $p$8 years (vs. IAS15’s $p$9s).
• In TRAPPIST-1 system stability (7-planet, near-resonant chain, $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$0 draws), 10 Myr runs confirm error growth $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$1 and $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$2.

Event-driven and Large-scale Applications

• Poincaré-section computation for the Hénon–Heiles system shows heyoka baseline performance is $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$3–$x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$4× faster than DifferentialEquations.jl and $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$5–$x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$6× faster than SciPy DOP853.
• For collision-rich N-body problems and eclipse detection on irregular small bodies, event detection adds modest overhead ($x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$720% per step in collision-prone Solar System tests), with no missed or falsely flagged events (Biscani et al., 2022).

Extended and Arbitrary Precision

When operated in 80-bit or 128-bit precision, heyoka achieves energy control at or below $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$8 (80-bit) and $x(t_0 + h) \approx \sum_{k=0}^p x^{[k]}(t_0) h^k,$9 (128-bit), with timings of $x^{[k]}(t) = \frac{1}{k!} \frac{d^k x}{dt^k}$0s and $x^{[k]}(t) = \frac{1}{k!} \frac{d^k x}{dt^k}$1s, respectively, for standard Solar System benchmarks (Biscani et al., 2021).

7. Significance and Scope

Heyoka demonstrates that a general-purpose, high-order, adaptive Taylor integrator, equipped with modern automatic differentiation, event detection, JIT compilation, and parallelism, can both rival and surpass traditional symplectic and non-symplectic methods (such as WHFast, IAS15, RKF7/8) across challenging gravitational, collisional, and chaotic regimes. It particularly excels in high-precision, multi-million year integrations, stiff/forced systems, close encounters, and event-driven hybrid dynamics. The algorithmic transparency, dense output, and open-source availability permit direct integration into computational workflows for astrodynamics, planetary science, and dynamical astronomy (Biscani et al., 2021, Biscani et al., 2022).