SIMD-vectorized implicit symplectic integrators can outperform explicit ones

TL;DR

The paper presents a $16$-th order, $8$-stage Gauss–Legendre collocation IRK method (IRKGL16) implemented with SIMD acceleration (IRKGL16-SIMD) in Julia to solve non-stiff Hamiltonian ODEs at high precision. It leverages a symplecticly correct floating-point reformulation and a vectorized fixed-point iteration with an energy-stable stopping criterion, achieving near-transparent SIMD parallelism across $8$-wide vectors. Through comprehensive benchmarks against state-of-the-art explicit symplectic integrators on problems including Schwarzschild-like first-order systems and second-order models such as the outer Solar System and Hénon–Heiles, the authors show that IRKGL16-SIMD can surpass explicit methods in accuracy for a given CPU time at high precision. The work demonstrates that implicit schemes, when properly vectorized, offer practical advantages for long-time, high-accuracy Hamiltonian integrations, expanding the toolbox for structure-preserving numerical integration.

Abstract

The main purpose of this work is to present a SIMD-vectorized implementation of the symplectic 16th-order 8-stage implicit Runge-Kutta integrator based on collocation with Gauss-Legendre nodes (IRKGL16-SIMD), and to show that it can outperform state-of-the-art symplectic explicit integrators for high-precision numerical integrations (in double-precision floating-point arithmetic) of non-stiff Hamiltonian ODE systems. Our IRKGL16-SIMD integrator leverages Single Instruction Multiple Data (SIMD) based parallelism (in a way that is transparent to the user) to significantly enhance the performance of the sequential IRKGL16 implementation. We present numerical experiments comparing IRKGL16-SIMD with state-of-the-art high-order explicit symplectic methods for the numerical integration of several Hamiltonian systems in double-precision floating-point arithmetic.

Figures (6)

• Figure 1: Work-precision diagrams (maximum local error in the Hamiltonian) for the Schwarzschild black hole problem using the IRKGL16 method and symmetric composition schemes of orders $r = 2, 4, 6, 8$, and $10$. In the right plot, the sequential IRKGL16 implementation is denoted as IRKGL16-SEQ
• Figure 2: Left: global error in the Hamiltonian. Right: relative error in the radial coordinate $r$ for the Schwarzschild black hole problem. Results are shown for IRKGL16-SIMD and symmetric composition methods of orders $r = 2, 4, 6, 8, 10$.
• Figure 3: Work-precision diagrams (maximum local Hamiltonian error) for the 6-body outer Solar System and the Hénon-Heiles problem using IRKGL16, symmetric compositions (orders $r = 2, 4, 6, 8, 10$), and symplectic RKN methods of orders 6 and 8.
• Figure 4: Time evolution of the global Hamiltonian error for the 6-body outer Solar System (left) and the Hénon-Heiles problem (right), using IRKGL16, symmetric compositions, and RKN methods.
• Figure 5: Evolution of position error for each planet in the 6-body outer Solar System using IRKGL16 and SS05 ($r=10$).