A Krylov subspace algorithm for evaluating the phi-functions appearing in exponential integrators

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 phi-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.Comment: 20 pages, 3 colour figures, code available from http://www.maths.leeds.ac.uk/~jitse/software.html . v2: Various changes to improve presentation as suggested by the refere

Exponential Rosenbrock methods without order reduction when integrating nonlinear initial value problems

A technique is described in this paper to avoid order reduction when integrating reaction-diffusion initial boundary value problems with explicit exponential Rosenbrock methods. The technique is valid for any Rosenbrock method, without having to impose any stiff order conditions, and for general time-dependent boundary values. An analysis on the global error is thoroughly performed and some numerical experiments are shown which corroborate the theoretical results, and in which a big gain in efficiency with respect to applying the standard method of lines can be observed.Comment: 27 pages, 2 figure

A robust spectral method for pricing of American put options on zero-coupon bonds

American put options on a zero-coupon bond problem is reformulated as a linear complementarity problem of the option value and approximated by a nonlinear partial differential equation. The equation is solved by an exponential time differencing method combined with a barycentric Legendre interpolation and the Krylov projection algorithm. Numerical examples shows the stability and good accuracy of the method. A bond is a financial instrument which allows an investor to loan money to an entity (a corporate or governmental) that borrows the funds for a period of time at a fixed interest rate (the coupon) and agrees to pay a fixed amount (the principal) to the investor at maturity. A zero-coupon bond is a bond that makes no periodic interest payments

Adaptive rational Krylov methods for exponential Runge--Kutta integrators

We consider the solution of large stiff systems of ordinary differential equations with explicit exponential Runge--Kutta integrators. These problems arise from semi-discretized semi-linear parabolic partial differential equations on continuous domains or on inherently discrete graph domains. A series of results reduces the requirement of computing linear combinations of $\varphi$-functions in exponential integrators to the approximation of the action of a smaller number of matrix exponentials on certain vectors. State-of-the-art computational methods use polynomial Krylov subspaces of adaptive size for this task. They have the drawback that the required Krylov subspace iteration numbers to obtain a desired tolerance increase drastically with the spectral radius of the discrete linear differential operator, e.g., the problem size. We present an approach that leverages rational Krylov subspace methods promising superior approximation qualities. We prove a novel a-posteriori error estimate of rational Krylov approximations to the action of the matrix exponential on vectors for single time points, which allows for an adaptive approach similar to existing polynomial Krylov techniques. We discuss pole selection and the efficient solution of the arising sequences of shifted linear systems by direct and preconditioned iterative solvers. Numerical experiments show that our method outperforms the state of the art for sufficiently large spectral radii of the discrete linear differential operators. The key to this are approximately constant rational Krylov iteration numbers, which enable a near-linear scaling of the runtime with respect to the problem size