Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods
Introduction
The presence of a wide range of temporal scales in a system of differential equations poses a major difficulty for their integration over long time intervals. Such stiff systems are routinely encountered in scientific applications from plasma modeling to chemical kinetics. The development of numerical techniques which provide computational savings over commonly used algorithms can allow one to solve problems faster and access previously unattainable parameter regimes.
A major difficulty in solving large stiff systems of nonlinear differential equations is choosing an efficient time integration scheme. Typically one has to make a decision whether to use an explicit or an implicit method. Explicit schemes require the least amount of computation per time step but the allowable time step is severely restricted by the stability requirements. Implicit schemes have much better stability properties and allow significantly larger time steps compared to explicit integrators. However, this advantage comes at the expense of a significant increase in the amount of computation needed at each time iteration. In particular, a typical choice for solving such problems is a Newton–Krylov implicit integrator. For large-scale stiff systems the cost of this method is dominated by the solutions of large linear systems at each Newton iteration. Experience shows that for general large-scale nonsymmetric problems GMRES is a natural choice for solving the linear systems; however, unless one can exploit the structure of the problem and develop a good preconditioner for these matrices, a Newton–Krylov method can become prohibitively expensive [1]. If it is possible to construct a good preconditioner for a particular problem the Newton–Krylov method becomes a very efficient way to solve the stiff system and is difficult to outperform. But for many problems constructing an effective preconditioner is a non-trivial task and one would like to have an alternative method that would provide savings compared to both explicit schemes and implicit Newton–Krylov integrators and would not require developing a preconditioner. These are the classes of problems where exponential propagation method can become advantageous.
The idea to use exponential time differencing (ETD) to construct an effective integrator for stiff systems has a long history and has been introduced and reintroduced in the literature many times (see Section 3). However, only when it was suggested to combine this idea with the Krylov subspace approximation of functions of matrices has it become viable to use ETD to construct time integrators for large-scale nonlinear stiff systems – the effort to construct such methods is a relatively new development. Due to the novelty of these ideas, limited understanding of their performance, and lack of well tested schemes which compare favorably to standard integrators, these methods have not yet been widely used. In fact, to our knowledge none of the previously developed exponential integrators have been clearly demonstrated to outperform the Newton–Krylov implicit integrators which are methods of choice when a large stiff nonlinear system of ODEs has to be solved.
In this paper, we introduce a new class of exponential propagators which we call exponential propagation iterative (EPI) methods. The methods are based on a key observation: if a Krylov projection is used to approximate a product f(A)b between a function f of a large stiff matrix A and a vector b then the convergence rate of the Krylov iteration will depend on the properties of the function that has to be approximated. Specifically, if a Newton–Krylov implicit integrator is used to solve a large stiff system of ODEs then at each Newton iteration the Krylov projections are used to approximate the product f(A)b = (I − A)−1b, where A is the Jacobian matrix multiplied by a time step, I is the identity matrix and b is a vector. So in this case the Krylov projections are used to approximate products of a rational function f(x) = 1/(1 − x) of a matrix and a vector. If one uses currently available exponential integrators then the Krylov projections are used to approximate products between a vector and an exponential of the matrix eAb or an expression (eA − I)A−1b. To get an accurate approximation to the solution these operations have to be performed several times per time step. We approach the construction of an exponential integrator by expressing the solution in terms of a set of special functions gk(x) or ϕγk(x), k = 0, 1, 2, …, γ (see (28), (29)). The advantage of an EPI method comes from the fact that the number of Krylov iterations needed to approximate products of these functions of a matrix with a vector decreases as the index k increases, and for each of these function the number of needed Krylov projection steps is in general smaller than the number of such steps required to approximate f(A)b with f(x) = 1/(1 − x) and f(x) = ex. Thus we are able to achieve computational savings compared to an implicit integrator while allowing much larger time steps than explicit schemes.
Below we will (i) describe the ideas behind constructing EPI methods and give an overview of other exponential integrators, (ii) introduce new EPI methods and a methodology for their construction, (iii) discuss the efficient implementation of these techniques, and (iv) based on some test problems provide guidance as to what computational savings one can expect compared to standard explicit and implicit methods. In Section 2, we present the history and the ideas behind exponential propagation methods. A procedure for constructing EPI schemes is given in Section 3. Here we introduce several new methods and discuss how the schemes should be formulated and implemented to be efficient. Finally, in Section 4 several numerical examples are used to illustrate the performance of the schemes; a discussion of their appropriate application is included.
Section snippets
History and development of exponential propagation techniques
We begin with a general overview of exponential integrators. Consider an initial-value problem for a large nonlinear autonomous system of ordinary differential equationswhere U(t) = (u1(t), u2(t), …, uN(t))T is the solution of (1) at time t and F(U) = (f1(U), f2(U), …, fN(U))T. If such system comes from a discretization of a partial differential equation, F(U) is the discrete representation of a spatial differential operator and the elements of U contain solution values at each grid
Construction of the exponential propagation iterative methods
We begin developing an exponential propagation iterative method from the integral form of the solution to (1)Recall that in this formula
is the solution of the system (1) at time tn,
is the right-hand side of the system (1) evaluated at time tn,
and
Numerical examples
To test the EPI methods and compare them with explicit and implicit schemes we consider the following three problems commonly used to test numerical methods for stiff systems. Since we are only interested in these problems from the perspective of testing the performance of our numerical methods we will discuss neither the applications associated with these differential equations nor the full spectrum of the behavior exhibited by their solutions. In fact, where possible we will set initial and
Conclusions
In this paper, we have introduced a new class of exponential integrators which we called exponential propagation iterative (EPI) methods. We have discussed the methodology for the construction of these schemes and studied their performance on several test problems. We have demonstrated that the faster convergence of Arnoldi iterations needed by EPI schemes provide computational savings compared to standard implicit Newton–Krylov integrators with no preconditioning. It is not clear at this point
Acknowledgments
The author thanks Prof. Alexandre Chorin and Prof. G.I. Barenblatt for introducing her to an inspiring and supportive scientific environment, Prof. Ole Hald and Dr. David Bernstein for helpful discussions, and Mr. John Barber for advice on Mathematica software. This work was supported in part by the NSF VIGRE Grant DMS-9819728, NSF/DOE partnership in plasma physics Grant DMS-0317511, and University of California President’s Postdoctoral Fellowship Program.
References (38)
- et al.
Jacobian-free Newton–Krylov methods: a survey of approaches and applications
J. Comput. Phys.
(2004) - et al.
A new class of time discretization schemes for the solution of nonlinear PDEs
J. Comput. Phys.
(1998) - et al.
Exponential time differencing for stiff systems
J. Comput. Phys.
(2002) Generalized integrating factor methods for stiff PDEs
J. Comput. Phys.
(2005)- et al.
Exponential Runge–Kutta methods for parabolic problems
Appl. Numer. Math.
(2005) An iterative solution method for solving f(a)x = b using Krylov subspace information obtained for the symmetric positive definite matrix a
J. Comput. Appl. Math.
(1987)- et al.
Chemical instabilities and sustained oscillations
J. Theor. Biol.
(1971) - et al.
Solving Ordinary Differential Equations
(1987) - S.M. Cox, Exponential time differencing (ETD) references. Available from:...
- J. Certaine, The solution of ordinary differential equations with large time constants, in: A. Ralston, H.S. Wilf...
An exponential method of numerical integration of ordinary differential equations
Commun. ACM
Generalized Runge–Kutta processes for stable systems with large Lipschitz constants
SIAM J. Numer. Anal.
An A-stable modification of the Adams–Bashforth methods
Lecture Notes Math.
Nineteen dubious ways to compute the exponential of a matrix
SIAM Rev.
Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later
SIAM Rev.
Efficient solution of parabolic equations by Krylov approximation methods
SIAM J. Sci. Stat. Comp.
Weighted quadrature in Krylov methods
Utilitas Mathematica
Fourth-order time stepping for stiff PDEs
SIAM J. Sci. Comp.
Iterative Methods for Sparse Linear Systems
Cited by (135)
Building blocks needed for mechanistic modeling of bioprocesses: A critical review based on protein production by CHO cells
2024, Metabolic Engineering CommunicationsBAMPHI: Matrix-free and transpose-free action of linear combinations of φ-functions from exponential integrators
2023, Journal of Computational and Applied MathematicsJacobian-free High Order Local Linearization methods for large systems of initial value problems
2023, Applied Numerical MathematicsComputing high dimensional multiple integrals involving matrix exponentials
2023, Journal of Computational and Applied MathematicsEfficient adaptive step size control for exponential integrators
2022, Computers and Mathematics with ApplicationsOn the stability of exponential integrators for non-diffusive equations
2022, Journal of Computational and Applied Mathematics