Computing a matrix function for exponential integrators

An efficient numerical method is developed for evaluating φ(A), where A is a symmetric matrix and φ is the function defined by φ(x)=(e^x-1)/x=1+x/2+x²/6+???. This matrix function is useful in the so-called exponential integrators for differential equations. In particular, it is related to the exact solution of the ODE system dy/dt=Ay+b, where A and b are t-independent. Our method avoids the eigenvalue decomposition of the matrix A and it requires about 10n³/3 operations for a general symmetric n×n matrix. When the matrix is tridiagonal, the required number of operations is only O(n²) and it can be further reduced to O(n) if only a column of the matrix function is needed. These efficient schemes for tridiagonal matrices are particularly useful when the Lanczos method is used to calculate the product of this matrix function (for a large symmetric matrix) with a given vector. © 2003 Elsevier B.V. All rights reserved.