Algorithm 919

Abstract

We develop an algorithm for computing the solution of a large system of linear ordinary differential equations (ODEs) with polynomial inhomogeneity. This is equivalent to computing the action of a certain matrix function on the vector representing the initial condition. The matrix function is a linear combination of the matrix exponential and other functions related to the exponential (the so-called ϕ-functions). Such computations are the major computational burden in the implementation of exponential integrators, which can solve general ODEs. Our approach is to compute the action of the matrix function by constructing a Krylov subspace using Arnoldi or Lanczos iteration and projecting the function on this subspace. This is combined with time-stepping to prevent the Krylov subspace from growing too large. The algorithm is fully adaptive: it varies both the size of the time steps and the dimension of the Krylov subspace to reach the required accuracy. We implement this algorithm in the MATLAB function phipm and we give instructions on how to obtain and use this function. Various numerical experiments show that the phipm function is often significantly more efficient than the state-of-the-art.

On the use of exponential time integration methods in atmospheric models

“…All of the methods presented require the evaluation of expressions of the form Niesen and Wright (2012) have developed the phipm algorithm to efficiently carry out this task. This is achieved by using an adaptive time-stepping on the non-autonomous linear ODE…”

“…This is done using a Krylov projection as outlined in the previous section. The full description of the algorithm is given in the paper of Niesen and Wright (2012). For a specified accuracy, the efficiency is determined both by the size of the Krylov subspaces used and the number of steps taken.…”

“…This appears in eqs. (17) and (18) of Niesen and Wright (2012) and relates to the initial choice of step size in the algorithm and the adaptive procedure for this.…”

“…The default value used by Niesen and Wright (2012) is 10 (7 . This parameter determines the accuracy of the evaluation of the 8-functions by controlling the dimension of the Krylov subspace and the length of the internal pseudo-time-step.…”

“…In particular, we discuss the adaptive Krylov subspace algorithm of Niesen and Wright (2012), which we will use for the efficient handling of matrix exponentials. In Section 3, we introduce the shallow water model and describe the implementation of the exponential methods.…”