How to apply phi-functions to matrices?

Question

I am implementing an explicit numerical integration technique, based on readings from Hochbruck and Ostermann (2010). They use the following formulation:

Excerpt from book

Here, the matrix A is a linear estimation of the Jacobian of the differential equation system, and g is the leftover nonlinear terms. It calls for the usage of a phi function, defined:

$$\phi_{1}(z)=\frac{e^{z}-1}{z}.$$

But it says to apply this to a square matrix? I understand the matrix exponential, but how do you divide by a square matrix? Should I multiply by the matrix inverse or divide term-by-term? Or is there another meaning?

Answer 1

These helper functions are defined to avoid the matrix inverse. Thus the function $\phi_1$ can also be applied to singular matrices. Insert the power series of the matrix exponential to find $$ \phi_1(A)=I+\tfrac12A+\tfrac16A^2+...+\tfrac1{n!}A^{n-1}+... $$ For efficient evaluation you would use the eigen-decomposition $A=U^{-1}JU$ and then compute $$ \phi_1(A)=U^{-1}\phi_1(J)U, $$ where the middle splits according to the Jordan blocks of $J$ and the usual rules for the evaluation of matrix functions on Jordan blocks apply.

More elaborate tricks to reduce floating-point errors and/or computation costs exist.

Answer 2

If you are interested in the action of matrix phi-functions to a vector, then you can also consider to compute the exponential of an extended matrix as shown in Theorem 1 in the paper [Expokit: A Software Package for Computing, Sidje R., ACM Trans. Math. Softw., 24(1):130-156, 1998] and others.

The following Matlab code computes $\varphi_p(tA)c\in\mathbb{C}^n$, i.e., the $p$-th phi-function of $tA\in\mathbb{C}^{n\times n}$ applied to a $c\in\mathbb{C}^n$, where $t\in\mathbb{R}$ is a possible time step, by evaluating the matrix exponential of an extended matrix $X\in\mathbb{C}^{n+p\times n+p}$:

X = A;
X(1:n,n+1) = c;
X(n+1:n+p-1,n+2:n+p)=eye(p-1);  
X(n+p,n+p)=0;
expX=expm(dt*X);
phipAc=expX(1:n,end)*dt^(-p); % = phi_p(t*A)*c

For $p=1$, a small choice of $n$, and a diagonalizable matrix $A$ with nonzero eigenvalues you can compute a reference solution using the eigenvalue decomposition of $A$ as following

[Q,lam] = eig(A,'vector');
philam = expm1(dt*lam)./(dt*lam);
phipAc_ref = Q*diag(philam)*inv(Q)*c;

(here the auxiliary $\texttt{expm1}$ computes the exponential minus 1 of the given vector and does improve the accuracy when eigenvalues of $A$ are close to zero.)

For further works on the computation of matrix phi-functions for applications in exponential integrators you can also have a look at [Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators, Higham and Al-Mohy, 2011.] or [Algorithm 919: A Krylov Subspace Algorithm for evaluating the $\phi$-Functions Appearing in Exponential Integrators, Niesen and Wright, 2012.]