A scaling and recovering algorithm for the matrix $\varphi$-functions

TL;DR

This paper introduces a new scaling and recovering algorithm for efficiently computing matrix i-functions, crucial in exponential integrator methods, with rigorous error bounds and structure-exploiting features.

Contribution

The paper presents a novel algorithm that combines scaling, Pad approximation, and recurrence relations for accurate, efficient computation of matrix i-functions, including error analysis and structure exploitation.

Findings

Algorithm outperforms existing methods in accuracy.

Reduces computational cost via optimal scaling and Pad degrees.

Provides sharp backward error bounds for nonnormal matrices.

Abstract

A new scaling and recovering algorithm is proposed for simultaneously computing the matrix $\varphi$-functions that arise in exponential integrator methods for the numerical solution of certain first-order systems of ordinary differential equations. The algorithm initially scales the input matrix down by a nonnegative integer power of two, and then evaluates the $[m/m]$ diagonal Pad\'e approximant to $\varphi_p$, where $p$ is the largest index of interest. The remaining $[m+p{-}j/m]$ Pad\'e approximants to $\varphi_j$, $0 \le j < p$, are obtained implicitly via a recurrence relation. The effect of scaling is subsequently recovered using the double-argument formula. A rigorous backward error analysis, based on the $[m+p/m]$ Pad\'e approximant to the exponential, enables sharp bounds on the relative backward errors. These bounds are expressed in terms of the sequence $\|A^k\|^{1/k}$, which can be much smaller than $\|A\|$ for nonnormal matrices. The scaling parameter and the degrees of the Pad\'e approximants are selected to minimize the overall computational cost, which benefits from the sharp bounds and the optimal evaluation schemes for diagonal Pad\'e approximants. Furthermore, if the input matrix is (quasi-)triangular, the algorithm exploits its structure in the recovering phase. Numerical experiments demonstrate the superiority of the proposed algorithm over existing alternatives in both accuracy and efficiency.