[Paper Review] 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 equations6 references3 citations

TL;DR

This paper presents a high-order multiple precision ODE solver using fully implicit Runge-Kutta (IRK) methods with embedded formulas and step size control, leveraging mixed-precision iterative refinement and W-transformation-based matrix reduction to accelerate inner iterations. The implementation achieves high accuracy (up to 29-digit precision) in stiff and large-scale problems like the 1D Brusselator, outperforming existing DP-only solvers in precision and efficiency.

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.

Motivation & Objective

• To develop a practical high-order multiple precision ODE solver capable of handling stiff and ill-conditioned problems beyond the limits of IEEE754 double precision.
• To overcome the inefficiency of RADAU5-type matrix reduction in high-order IRK methods by adopting SPARK3-type W-transformation for faster inner iteration.
• To enable adaptive step size control in multiple precision by deriving embedded formulas for arbitrary IRK methods.
• To accelerate linear system solves in inner iterations using mixed-precision iterative refinement, especially for sparse Jacobians in large-scale problems.
• To demonstrate the solver's high accuracy and performance on benchmark stiff ODEs, including the 1D Brusselator with 1000 DAEs.

Proposed method

• The solver uses high-order fully implicit Runge-Kutta methods with embedded formulas for local error estimation and step size control in multiple precision arithmetic.
• It applies the W-transformation to reduce the coefficient matrix of the linear system in inner iterations to real nonsymmetric tridiagonal form, improving computational efficiency over RADAU5’s complex diagonalization.
• Mixed-precision iterative refinement is employed, solving the linear system in double precision and correcting in multiple precision, drastically reducing solve time.
• The method supports sparse Jacobian matrices by using band matrix storage and left-preconditioned DP-MP BiCGSTAB solvers for large-scale problems.
• The implementation is built on BNCpack, a multiple precision numerical library using MPFR and GMP, enabling high-precision arithmetic with CPU-optimized kernels.
• Step size adaptation is achieved via embedded formulas incidental to the main IRK formula, enabling A-stable and high-order accurate integration.

Experimental results

Research questions

• RQ1 Can high-order IRK methods with embedded formulas be effectively implemented in a multiple precision environment to achieve higher accuracy than existing double-precision solvers?
• RQ2 Is the SPARK3-type W-transformation more effective than the RADAU5-type reduction for accelerating inner iterations in high-order IRK methods under multiple precision arithmetic?
• RQ3 Can mixed-precision iterative refinement significantly reduce the computational cost of solving large linear systems in the inner iteration of high-order IRK methods?
• RQ4 How does the performance and accuracy of the proposed MP ODE solver compare to existing solvers on stiff and large-scale problems like the 1D Brusselator?
• RQ5 What are the bottlenecks in MP ODE solvers, particularly regarding convergence of linear solvers in inner iterations, and how can they be mitigated?

Key findings

• The solver achieved a maximum relative error of $1.3 \times 10^{-29}$ on the 1D Brusselator problem with 20 stages and 20th order, indicating over 25 correct digits of accuracy.
• For the 1D Brusselator with $N=500$, the 20-stage method reduced computational time to 3642 seconds (vs. 64573 for unpreconditioned BiCGSTAB), showing the benefit of left preconditioning.
• The mixed-precision iterative refinement method enabled faster convergence in linear solves, especially when combined with band matrix storage for sparse Jacobians.
• The SPARK3-type W-transformation was more effective than RADAU5-type reduction in reducing inner iteration time for high-order IRK methods in multiple precision.
• The solver successfully computed solutions with over 20 true decimal digits of precision, confirmed by comparison with 100-digit MP results.
• Step size adaptation was stable and efficient due to embedded formulas, and the use of left-preconditioned DP-MP BiCGSTAB improved convergence and reduced computational time.