Computing the Wave-Kernel Matrix Functions

The wave-kernel functions $\cosh{\sqrt{A}}$ and $\mathrm{sinhc}{\sqrt{A}}$ arise in the solution of second order differential equations such as $u''(t) - Au(t) = b(t)$ with initial conditions at $t=0$. Here, $A$ is an arbitrary square matrix and $\mathrm{sinhc}{z} = \sinh(z)/z$. The square root in these formulas is illusory, as both functions can be expressed as power series in $A$, so there are no questions about existence of the functions.

How can these functions be computed efficiently? In Computing the Wave-Kernel Matrix Functions (SIAM J. Sci. Comput., 2018) Prashanth Nadukandi and I develop an algorithm based on Padé approximation and the use of double angle formulas. The amount of scaling and the degree of the Padé approximant are chosen to minimize the computational cost subject to achieving backward stability for $\cosh{\sqrt{A}}$ in exact arithmetic.

In the derivation we show that the backward error of any approximation to $\cosh{\sqrt{A}}$ can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for $\|A^k\|^{1/k}$ in terms of the norms of lower powers of $A$; this bound is sharper than one previously obtained by Al-Mohy and Higham.

Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur–Parlett algorithm applied to these functions.

The fundamental regions of the

function cosh(sqrt(z)), needed for the backward error analysis

underlying the algorithm.