[论文解读] Practical Implementation of High-Order Multiple Precision Fully Implicit Runge-Kutta Methods with Step Size Control Using Embedded Formula

Tomonori Kouya|arXiv (Cornell University)|Jun 11, 2013

Numerical methods for differential equations参考文献 6被引用 3

一句话总结

本论文提出了一种基于全隐式龙格-库塔（IRK）方法的高阶多精度常微分方程（ODE）求解器，结合嵌入公式与步长控制机制，利用混合精度迭代修正和基于W-变换的矩阵降阶技术，加速内层迭代过程。该实现可在刚性及大规模问题（如1D Brusselator）中达到高达29位十进制精度，其精度与效率均优于现有仅使用双精度（DP）的求解器。

ABSTRACT

We propose a practical implementation of high-order fully implicit Runge-Kutta(IRK) methods in a multiple precision floating-point environment. Although implementations based on IRK methods in an IEEE754 double precision environment have been reported as RADAU5 developed by Hairer and SPARK3 developed by Jay, they support only 3-stage IRK families. More stages and higher-order IRK formulas must be adopted in order to decrease truncation errors, which become relatively larger than round-off errors in a multiple precision environment. We show that SPARK3 type reduction based on the so-called W-transformation is more effective than the RADAU5 type one for reduction in computational time of inner iteration of a high-order IRK process, and that the mixed precision iterative refinement method is very efficient in a multiple precision floating-point environment. Finally, we show that our implementation based on high-order IRK methods with embedded formulas can derive precise numerical solutions of some ordinary differential equations.

研究动机与目标

• 开发一种实用的高阶多精度ODE求解器，能够处理超出IEEE754双精度范围的刚性与病态问题。
• 通过采用SPARK3型W-变换，替代高阶IRK方法中RADAU5型矩阵降阶方式，以提升内层迭代的计算效率，克服其低效问题。
• 通过推导任意IRK方法的嵌入公式，实现在多精度环境下的自适应步长控制。
• 利用混合精度迭代修正技术，加速内层迭代中线性系统的求解，尤其适用于大规模问题中的稀疏雅可比矩阵。
• 在基准刚性ODE问题（包括含1000个微分代数方程的1D Brusselator）上，验证求解器的高精度与高性能表现。

提出的方法

• 求解器采用高阶全隐式龙格-库塔方法，结合嵌入公式以实现多精度算术下的局部误差估计与步长控制。
• 通过W-变换将内层迭代中线性系统的系数矩阵转化为实非对称三对角形式，相比RADAU5的复对角化方式，显著提升了计算效率。
• 采用混合精度迭代修正策略，即在双精度下求解线性系统，并在多精度下进行校正，大幅缩短求解时间。
• 通过带状矩阵存储方式支持稀疏雅可比矩阵，并采用左预处理的DP-MP BiCGSTAB求解器，适用于大规模问题。
• 实现基于BNCpack（一个基于MPFR与GMP的多精度数值库）的多精度算术计算，利用CPU优化内核实现高精度运算。
• 通过与主IRK公式共存的嵌入公式实现步长自适应，确保A-稳定性和高阶精度的积分效果。

实验结果

研究问题

• RQ1 高阶IRK方法结合嵌入公式是否可在多精度环境下有效实现，从而实现比现有双精度求解器更高的精度？
• RQ2 在多精度算术下，SPARK3型W-变换是否比RADAU5型降阶方法更有效地加速高阶IRK方法的内层迭代？
• RQ3 混合精度迭代修正是否能显著降低高阶IRK方法内层迭代中大规模线性系统求解的计算成本？
• RQ4 所提出的多精度ODE求解器在刚性与大规模问题（如1D Brusselator）上的性能与精度，相较于现有求解器表现如何？
• RQ5 多精度ODE求解器的主要瓶颈是什么，特别是内层迭代中线性求解器的收敛问题，应如何缓解？

主要发现

• 在含20个阶段与20阶的1D Brusselator问题上，求解器达到最大相对误差为$1.3 \times 10^{-29}$，表明精度超过25位有效数字。
• 对于$N=500$的1D Brusselator问题，20阶段方法将计算时间缩短至3642秒（相较未预处理的BiCGSTAB求解器的64573秒），验证了左预处理的有效性。
• 混合精度迭代修正方法显著提升了线性系统求解的收敛速度，尤其在结合带状矩阵存储处理稀疏雅可比矩阵时效果更明显。
• SPARK3型W-变换在多精度环境下，相比RADAU5型降阶方式，更有效地减少了高阶IRK方法内层迭代的时间开销。
• 求解器成功计算出精度超过20位真实十进制数字的解，其结果经与100位多精度结果对比得到验证。
• 由于嵌入公式的存在，步长自适应过程稳定且高效；同时，采用左预处理的DP-MP BiCGSTAB求解器显著提升了收敛性并减少了计算时间。