Efficient Explicit Taylor ODE Integrators with Symbolic-Numeric Computing

Abstract: Taylor series methods show a newfound promise for the solution of non-stiff ordinary differential equations (ODEs) given the rise of new compiler-enhanced techniques for calculating high order derivatives. In this paper we detail a new Julia-based implementation that has two important techniques: (1) a general purpose higher-order automatic differentiation engine for derivative evaluation with low overhead; (2) a combined symbolic-numeric approach to generate code for recursively computing the Taylor polynomial of the ODE solution. We demonstrate that the resulting software's compiler-based tooling is transparent to the user, requiring no changes from interfaces required to use standard explicit Runge-Kutta methods, while achieving better run time performance. In addition, we also developed a comprehensive adaptive time and order algorithm that uses different step size and polynomial degree across the integration period, which makes this implementation more efficient and versatile in a broad range of dynamics. We show that for codes compatible with compiler transformations, these integrators are more efficient and robust than the traditionally used explicit Runge-Kutta methods.

• The paper introduces efficient explicit Taylor ODE integrators using a symbolic-numeric pipeline, reducing computational costs and allowing for arbitrary order and joint step size-degree adaptation.
• The method achieves Taylor-mode automatic differentiation (AD) costs of at most $O(p^2)$, making it more efficient compared to traditional AD methods such as Runge-Kutta, while requiring no user code modification.
• Numerical experiments confirm that on non-stiff benchmark problems, high-degree Taylor methods outperform traditional explicit Runge–Kutta methods in the high-accuracy regime, demonstrating a significant speedup over existing methods.

Motivation and problem statement

Taylor series methods for non-stiff ODEs offer two structural advantages over explicit Runge–Kutta schemes: the convergence order can be raised systematically by computing more terms of the solution polynomial, and dense output is available directly from the polynomial interpolant. Despite these advantages, Taylor methods have remained niche because computing high-order total derivatives of $f(u,t)$ is expensive, and because existing implementations (ATOMFT, TIDES, Jorba's taylor) require a separate build step—hand-written parsers or external computer algebra systems—that restricts how users may express their ODEs. The paper identifies a concrete complexity trap in naive approaches: if order-$k$ Taylor-mode AD costs $O(k^2)$ operations relative to $f$, then obtaining coefficients one at a time via repeated AD calls yields $O(p^3/3)$ total cost, which is prohibitive. Prior systems avoided this with special-purpose AD engines; Griewank et al.'s coefficient doubling achieves $O(p^2)$ but requires storing Jacobians $A_i = \partial u_i/\partial u_0$, giving $O(d^2)$ space that is unattractive for large systems.

Method: composed symbolic-numeric pipeline

The implementation rests on two components composed non-intrusively within Julia:

Taylor-mode automatic differentiation. The authors build on their prior engine (TaylorDiff.jl), which propagates truncated Taylor bundles $(h_0,\ldots,h_p)$ of higher-order directional derivatives through function compositions via Faà di Bruno-style pushforward rules. Pushforward rules are auto-generated from first-order rules defined in ChainRules.jl through Julia's multiple dispatch. The key guarantee is that for functions built from elementary operations with arbitrary control flow, the $p$-th order directional derivative costs at most $k$0 times evaluating $k$1 itself—in contrast to nested first-order forward-mode AD, which can be exponentially expensive ($k$2).

Compile-time symbolic recursion. Rather than calling the AD engine numerically at every step (which recomputes lower-order coefficients redundantly), the implementation executes the coefficient recursion symbolically at compile time: given target degree $k$3, it derives a directed computational graph expressing $k$4 purely as a function of $k$5, simplifies it, and JIT-compiles it into a numerical stepping function. This one-time compilation (microseconds to milliseconds) is amortized over long integrations. On the planar circular restricted three-body problem, symbolic simplification yields roughly an order-of-magnitude speedup over naive evaluation, with gains growing at higher polynomial degrees. A comparison against TaylorIntegration.jl—which relies on ad hoc macros and requires user-side code restructuring—shows the present approach is faster both before and after simplification while requiring no modification to user code. The authors concede they could not benchmark against a JAX-based Taylor solver, and infer from the JAX team's own report that coefficient doubling would perform comparably to the unsimplified naive approach here—an inference rather than a measured result.

Adaptive step size and polynomial degree

Step size control follows the standard reject-and-controller pattern: the scaled error

$k$6

is estimated from one extra derivative term (an approximation the paper states but does not bound rigorously); steps with $k$7 are rejected, otherwise an I, PI, or PID controller proposes $k$8.

The more distinctive contribution is joint adaptation of degree and step size. Since per-step work scales as $k$9 under Taylor-mode AD, the controller defines a work-per-time ratio $O(k^2)$0 for each candidate degree and selects the minimizer, restricting degree changes to $O(k^2)$1 between steps to suppress oscillation. This exploits a property unique to Taylor methods: the degree is a continuous, cheaply adjustable parameter, unlike the fixed orders of embedded Runge–Kutta pairs.

Numerical experiments

On four non-stiff benchmark problems (Lotka–Volterra, FitzHugh–Nagumo, a rigid body system, and a random 16-dimensional linear system) integrated over $O(k^2)$2, fixed-degree Taylor methods at degrees 6, 8, 10, and 12 are compared against DP5, Tsit5, Vern6, and Vern8. The reported finding is that work-precision curves are similar for matched orders (e.g., Taylor6 vs. Vern6, Taylor8 vs. Vern8), but because the Taylor degree can be raised arbitrarily, higher-degree Taylor methods dominate at tighter tolerances. This supports the claim that for compiler-transformable codes, these integrators are more efficient than traditional explicit Runge–Kutta methods specifically in the high-accuracy regime—at loose tolerances the advantage disappears.

The adaptive-degree experiment uses a logistic-type problem $O(k^2)$3 whose analytic solution has complex singularities at distance $O(k^2)$4, creating a "hard" region near the peak of $O(k^2)$5 where large degrees cannot enlarge the radius-limited step, and an "easy" region where they can. Here the adaptive-degree method (degree varying between 6 and 12) outperforms all four fixed-degree methods at lower tolerances, confirming that degree adaptivity pays off precisely when the local convergence radius varies substantially along the trajectory.

Limitations and open questions

The paper is candid that explicit Taylor methods inherit the standard limitation of all explicit integrators: on significantly stiff problems, stability forces very small steps regardless of polynomial degree, negating the efficiency advantages demonstrated here. The adaptive-degree analysis also depends on the assumption that per-step cost scales as $O(k^2)$6—valid for the current Taylor-mode AD engine but an assumption inherited from its implementation rather than a theorem about the method. Additionally, the comparison to JAX's coefficient-doubling approach is inferred rather than measured, and the error estimate used for step control approximates the Lagrange remainder by one additional derivative term without a formal error bound. Whether implicit Taylor variants can be constructed with comparable symbolic-numeric overhead remains unresolved; the authors identify this as the natural extension, alongside the observation that TaylorIntegration.jl retains an advantage for interval-arithmetic applications such as global error certification.

Conclusion

This paper demonstrates that composing a general-purpose Taylor-mode AD engine with a general-purpose computer algebra system inside a JIT-compiling language yields explicit Taylor ODE integrators that match Runge–Kutta performance at moderate accuracy, dominate at high accuracy due to arbitrarily adjustable order, require no user code transformation, and benefit further from joint step-size/degree adaptation when the local convergence radius varies. The results are confined to non-stiff problems and rest on empirical work-precision comparisons; extending the approach to implicit formulations is the principal open direction the paper itself poses.