Jitse Niesen
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.
Supplemental Material
Available for Download
Software for A Krylov Subspace Algorithm for Evaluating the ϕ-Functions Appearing in Exponential Integrators
- Al-Mohy, A. H. and Higham, N. J. 2011. Computing the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput. 33, 488--511. Google Scholar
Digital Library
- Butcher, J. C. 2008. Numerical Methods for Ordinary Differential Equations 2nd Ed. John Wiley & Sons Ltd., Chichester.Google Scholar
- Caliari, M. and Ostermann, A. 2009. Implementation of exponential Rosenbrock-type integrators. Appl. Numer. Math. 59, 568--581. Google ScholarDigital Library
- Crank, J. and Nicolson, P. 1947. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Proc. Cambridge Philos. Soc. 43, 50--67.Google ScholarCross Ref
- Davis, T. A. 2007. The University of Florida sparse matrix collection. Tech. rep. REP-2007-298, CISE Department, University of Florida.Google Scholar
- Douglas Jr., J. and Rachford Jr., H. H. 1956. On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc. 82, 421--439.Google ScholarCross Ref
- Duff, I. S., Grimes, R. G., and Lewis, J. G. 1989. Sparse matrix test problems. ACM Trans. Math. Softw. 15, 1, 1--14. Google ScholarDigital Library
- Friesner, R. A., Tuckerman, L. S., Dornblaser, B. C., and Russo, T. V. 1989. A method for the exponential propagation of large stiff nonlinear differential equations. J. Sci. Comput. 4, 4, 327--354. Google ScholarDigital Library
- Gallopoulos, E. and Saad, Y. 1992. Efficient solution of parabolic equations by Krylov approximation methods. SIAM J. Sci. Statist. Comput. 13, 1236--1264. Google ScholarDigital Library
- Hairer, E., Nørsett, S. P., and Wanner, G. 1993. Solving Ordinary Differential Equations. I 2nd Ed. Springer Series in Computational Mathematics, vol. 8, Springer-Verlag, Berlin. Google ScholarDigital Library
- Heston, S. L. 1993. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Rev. Finan. Stud. 6, 2, 327--43.Google ScholarCross Ref
- Higham, N. J. 2005. The scaling and squaring method for the matrix exponential revisited. SIAM J. Matrix Anal. Appl. 26, 4, 1179--1193. Google ScholarDigital Library
- Higham, N. J. 2008. Functions of Matrices. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA. Google ScholarDigital Library
- Hochbruck, M. and Lubich, C. 1997. On Krylov subspace approximations to the matrix exponential operator. SIAM J. Numer. Anal. 34, 5, 1911--1925. Google ScholarDigital Library
- Hochbruck, M., Lubich, C., and Selhofer, H. 1998. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput. 19, 5, 1552--1574. Google ScholarDigital Library
- Horn, R. A. and Johnson, C. R. 1991. Topics in Matrix Analysis. Cambridge University Press, Cambridge, UK. Google ScholarDigital Library
- Hundsdorfer, W. and Verwer, J. 2003. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer Series in Computational Mathematics, vol. 33, Springer-Verlag, Berlin.Google Scholar
- In ’t Hout, K. J. 2007. ADI schemes in the numerical solution of the Heston PDE. Numer. Anal. Appl. Math. 936, 10--14.Google ScholarCross Ref
- In ’t Hout, K. J. and Foulon, S. 2010. ADI finite difference schemes for option pricing in the Heston model with correlation. Int. J. Numer. Anal. Model. 7, 2, 303--320.Google Scholar
- In ’t Hout, K. J. and Welfert, B. D. 2009. Unconditional stability of second-order ADI schemes applied to multi-dimensional diffusion equations with mixed derivative terms. Appl. Numer. Math. 59, 4, 677--692. Google ScholarDigital Library
- Koikari, S. 2007. An error analysis of the modified scaling and squaring method. Computers Math. Applic. 53, 8, 1293--1305. Google ScholarDigital Library
- Koikari, S. 2009. On a block Schur--Parlett algorithm for ϕ-functions based on the sep-inverse estimate. ACM Trans. on Math. Soft. 38, 2, Article 12. Google ScholarDigital Library
- López-Fernández, M. 2010. A quadrature based method for evaluating exponential-type functions for exponential methods. BIT 50, 3, 631--655.Google ScholarCross Ref
- Minchev, B. and Wright, W. M . 2005. A review of exponential integrators for semilinear problems. Tech. rep. 2/05, Department of Mathematical Sciences, NTNU, Norway. http://www.math.ntnu.no/preprint/.Google Scholar
- Moler, C. B. and Van Loan, C. F. 2003. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review 45, 1, 3--49.Google ScholarDigital Library
- Moret, I. 2007. On RD-rational Krylov approximations to the core-functions of exponential integrators. Numer. Linear Algebra Appl. 14, 5, 445--457.Google ScholarCross Ref
- Niesen, J. and Wright, W. M. 2011. A Krylov subspace method for option pricing. http://ssrn.com/abstract=1799124.Google Scholar
- Saad, Y. 1992. Analysis of some Krylov subspace approximations to the matrix exponential operator. SIAM J. Numer. Anal. 29, 1, 209--228. Google ScholarDigital Library
- Schmelzer, T. and Trefethen, L. N. 2007. Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals. Electron. Trans. Numer. Anal. 29, 1--18.Google Scholar
- Sidje, R. 1998. EXPOKIT: Software package for computing matrix exponentials. ACM Trans. Math. Soft. 24, 1, 130--156. http://www.maths.uq.edu.au/expokit/ for latest software. Google ScholarDigital Library
- Skaflestad, B. and Wright, W. M. 2009. The scaling and modified squaring method for matrix functions related to the exponential. Appl. Numer. Math. 59, 3--4, 783--799. Google ScholarDigital Library
- Sofroniou, M. and Spaletta, G. 2007. Efficient computation of Krylov approximations to the matrix exponential and related functions. Personal communication.Google Scholar
- Tavella, D. and Randall, C. 2000. Pricing Financial Instruments. John Wiley & Sons Ltd., New York.Google Scholar
- Trefethen, L. N. and Bau III, D. 1997. Numerical Linear Algebra. SIAM, Philadelphia, PA.Google Scholar
Index Terms
- Algorithm 919: A Krylov Subspace Algorithm for Evaluating the ϕ-Functions Appearing in Exponential Integrators
Recommendations
Exponential Integrators for Large Systems of Differential Equations
We study the numerical integration of large stiff systems of differential equations by methods that use matrix--vector products with the exponential or a related function of the Jacobian. For large problems, these can be approximated by Krylov subspace ...
The scaling and modified squaring method for matrix functions related to the exponential
In recent years, there has been a resurgence in the construction and implementation of exponential integrators, which are numerical methods specifically designed for the numerical solution of spatially discretized semi-linear partial differential ...
Approximation of matrix operators applied to multiple vectors
In this paper we propose a numerical method for approximating the product of a matrix function with multiple vectors by Krylov subspace methods combined with a QR decomposition of these vectors. This problem arises in the implementation of exponential ...
Comments