A Fully Adaptive Radau Method for the Efficient Solution of Stiff Ordinary Differential Equations at Low Tolerances

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.5x 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.