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

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.

• The paper presents a Julia implementation of a fully adaptive Radau method for stiff ODEs, dynamically adjusting order and precision.
• Benchmarking shows the Julia method provides approximately double the performance efficiency compared to older Fortran versions on standard stiff ODE test problems.
• The method's order selection algorithm and arbitrary precision handling enable robust and efficient solutions, particularly at low tolerances crucial for high-fidelity simulations.

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

This paper introduces a sophisticated approach to solving stiff ordinary differential equations (ODEs) using an advanced Radau IIA method implemented in Julia. The authors address certain limitations of established methods, particularly the adaptive order Radau method originally in Fortran by Hairer, which was constrained by predefined tableau orders and double precision floating point calculations. By introducing an adaptive-time and adaptive-order Radau method, the paper proposes improvements in computational efficiency and precision.

Key Points and Methodology

1. Adaptive Order and Precision: The method derives coefficients for the Radau IIA tableau dynamically, allowing any order and precision. This contrasts with the fixed orders in traditional implementations (5th, 9th, and 13th).
2. Integration with Modern Computing: The implementation leverages Julia, exploiting its strengths in high-level programming and numerical efficiency through arbitrary precision arithmetic and parallelization. This approach overcomes the limitations of classic Fortran implementations by utilizing features such as multithreading and SIMD (Single Instruction, Multiple Data) registers.
3. Enhanced Performance: Benchmarks indicate that the new Julia-based Radau implementation offers approximately double the performance efficiency of Hairer's Fortran version across numerous stiff ODE test problems.
4. Order and Step-size Adaptivity: The use of a novel order selection algorithm, based on convergence iteration counts, refines the adaptivity further compared to existing methods which rely primarily on divergence values.
5. Numerical Evidence: The paper provides detailed work-precision diagrams benchmarking the adaptive Radau against several well-known numerical integrators, including BDF and Rosenbrock methods. Results consistently show the Julia implementation achieving notable speed-ups, particularly at low tolerances—a distinct advantage of the method.

Practical and Theoretical Implications

The contributions in this paper hold significant implications for both computational efficiency and the precision of solutions to stiff ODEs. The ability to handle arbitrary precision and dynamically adjust both time-step and order can greatly enhance model robustness in fields like chemical kinetics, astrochemistry, and pharmacometrics.

Furthermore, the modular and open-source nature of the implementation means that it can be readily adapted and extended by others in the computational mathematics community. This extends its utility beyond specific problems to a wide range of applications requiring high-fidelity simulations.

Future Perspectives

• Broader Adoption: As computational requirements grow and more researchers adopt Julia for scientific computing, the methodologies developed here could become a cornerstone for high-precision simulations.
• Methodological Extensions: Future work could explore extending these adaptive strategies to other types of differential equations systems, potentially including partial differential equations (PDEs).
• Parallel Computing: There's also fertile ground for further research into optimizing the Radau implementation on emerging hardware architectures.

Conclusion

This paper's contribution lies in transforming the established Radau IIA methods into a versatile tool better aligned with modern computational infrastructure and practices. By tackling long-standing limitations of order flexibility and integrating modern computing capabilities, the authors present a substantial improvement on traditional approaches to solving stiff ODEs, demonstrated through comprehensive benchmarking and a robust theoretical framework.