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

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.