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 highly accurate 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.
翻译:Radau IIA方法,特别是由Hairer开发的Fortran自适应阶Radau方法,被公认为是高精度求解高度刚性常微分方程(ODEs)的最先进技术。然而,传统实现专门针对特定容差范围,特别是仅支持5阶、9阶和13阶的Butcher表,且仅以双精度浮点数推导,这限制了其在需要高精度场景下的通用性。为缓解这些限制,我们实现了一种自适应时间步长与自适应阶数的Radau方法,能够动态生成任意阶数、任意精度的Radau IIA嵌入Butcher表系数。此外,我们基于Julia语言的实现包含多项现代化改进以提升性能,包括阶数自适应方案的优化和改进的线性代数集成。在与经典Fortran实现的直接基准测试中,我们证明我们的实现在一系列刚性ODE上具有约2倍的加速效果。我们将该算法与多种公认的刚性ODE数值积分器进行对比测试,在多个测试问题上展现了最先进的性能,在需要低误差容差时比常见的刚性ODE数值积分器快1.5倍。新实现的方法已在开源软件中发布,供用户免费用于刚性ODE求解。