Matrix exponential

A rational Chebyshev method to compute the exponential of a symmetric matrix.

This method is based on expanding the exponential function on the negative real axis into a rational Chebyshev series. Calling rcexpmv(A,v) returns an approximation to expm(A)*v, i.e., the action of the matrix exponential of A onto a vector v (or a block of vectors). Note that this is the solution u(1) of the linear initial value problem u'(t) = A*u(t), u(0) = v. Of course, you can also compute the full matrix exponential by setting v = eye(size(A)).

The square matrix A has do be symmetric with no positive eigenvalues! If there are positive eigenvalues, just shift A sufficiently far to the left using

expm(A) = exp(s) * expm( A - s*speye(size(A)) ), s > 0.

This method is much faster than Matlab's expm if A is large. It requires the repeated solution of a single shifted linear system with A for different right-hand sides, and the shift is real. This means that

1) for real data b and A the linear system to be solved is real, and
2) a direct solver would only factor the matrix exactly once and reuse the factorization (we simply use backslash, but you can easily use your favorite linear system solver).

An error tolerance can be provided by calling rcexpmv(A,v,tol). A larger tolerance reduces the number of required linear system solves. The tolerance also influences the optimal shift, which is chosen automatically.

Included is also a routine rcexpmv2 to compute exp(M\K)*v, a problem that arises when solving a linear initial value problem with a mass matrix: M*u'(t) = K*u(t), u(0) = v.

To get started, have a look at test_rcexmpv.m or test_rcexpmv2.m.

Rational Chebyshev expansions have been first advocated in the context of spectral methods by [Grosch & Orszag 77] and [Boyd 82/87]. Old and new rational Krylov methods for approximating matrix functions are discussed in my thesis [Güttel 10] (available online).