Computing matrix $\varphi$-functions arising in exponential integrators

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 (ODEs). The algorithm initially scales the input matrix down by a nonnegative integer power of two, then computes 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 &lt; 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<sup>k|<sup>{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 a priori sharpness of the 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.