A Fully Adaptive Radau Method for the Efficient Solution of Stiff Ordinary Differential Equations at Low Tolerances
Shreyas Ekanathan, Oscar Smith, Christopher Rackauckas
TL;DR
The paper addresses the challenge of solving stiff ODEs at very low tolerances by enabling truly arbitrary-order Radau IIA methods. It introduces a Julia-based adaptive-order, adaptive-time Radau solver that derives tableau coefficients on the fly, supports arbitrary precision, and leverages a-transformed LU factorizations and embedded error estimation to achieve superior efficiency. Key contributions include on-demand tableau generation with caching, an embedded error estimator with a simple $\sum_{i=1}^s \tilde{b}_i c_i^{k-1} = \frac{1}{k}$ relation, and an order-adaptive controller that balances convergence and cost. Results show state-of-the-art performance on stiff benchmark problems at low tolerances, with substantial speedups over Hairer's Radau and other stiff solvers, while remaining open source and compatible with advanced Julia tooling. This work raises the standard for adaptive, high-order stiff solvers and demonstrates practical impact for high-fidelity simulations requiring extreme accuracy in stiff regimes.
Abstract
Radau IIA methods, specifically the adaptive order Radau method in Fortran due to Hairer, are known to be state-of-the-art for the high-accuracy solution of highly stiff ordinary differential equations (ODEs). However, the traditional implementation was specialized to a specific range of tolerance, in particular only supporting 5th, 9th, and 13th order versions of the tableau and only derived in double precision floating point, thus limiting the ability to be truly general purpose for high fidelity scenarios. To alleviate these constraints, we implement an adaptive-time adaptive-order Radau method which can derive the coefficients for the Radau IIA embedded tableau to any order on the fly to any precision. Additionally, our Julia-based implementation includes many modernizations to improve performance, including improvements to the order adaptation scheme and improved linear algebra integrations. In a head-to-head benchmark against the classic Fortran implementation, we demonstrate our implementation is approximately 2x across a range of stiff ODEs. We benchmark our algorithm against several well-reputed numerical integrators for stiff ODEs and find state-of-the-art performance on several test problems, with a 1.5-times speed-up over common numerical integrators for stiff ODEs when low error tolerance is required. The newly implemented method is distributed in open source software for free usage on stiff ODEs.