Exponential integrators

Recently, exponential integrators (see a recent review by Hochbruck and Osterman published in Acta Numerica in 2010), a time integration technique that makes explicit use of the matrix exponential and related matrix functions of the Jacobian (or an approximation to it), have emerged as an efficient alternative to standard time integrators for solving stiff systems of the form (e.g. upon spatial discretizations of time-dependent PDEs)

$U'(t) =F(t, U(t)), \quad U(t_0) =U_0$.

They are fully explicit and do not suffer from the stability restriction imposed by the CFL condition for the linear part. It has also been shown that exponential integrators can take much larger time steps than implicit/IMEX methods while maintaining the same level of accuracy. Thus they can offer significant computational savings, particularly for large-scale stiff systems where no efficient preconditioner is available.

The state-of-the-art:
Non-stiff (or classical) exponential integrators have been studied intensively and derived by many authors. The term non-stiff here means that their whole construction was based on classical Taylor series expansions of matrix functions with remainder terms, however, that are not bounded independently of the stiffness. In particular, the matrix function like $e^{hA}$  is expanded as $e^{hA}=I + hA+ \ldots + \frac{(hA)^p}{p!} +\mathcal{O}(h^{p+1})$, which is only sensible for non-stiff problems where A is of small norm. This results in the so-called non-stiff order conditions (can be extended straightforwardly up to arbitrary order like explicit Runge-Kutta methods). For stiff problems, however, the Jacobian of the vector field often has a large norm or is even an unbounded operator (as stated above). Therefore, a very careful error analysis must be performed to make sure that the error terms do not contain powers of this Jacobian. Numerical experiments also showed that, depending on the problem, non-stiff integrators suffered from an order reduction when applied to stiff problems. These observations motivated the search for stiff order conditions which give rise to the construction of stiffly accurate exponential integrators. Although particular examples of exponential integrators were already discovered more than 50 years ago, the first error analysis and derivation of stiffly accurate exponential integrators had to wait another 40 years and were finally given in Hochbruck and Osterman (2005 –paper in SINUM). In particular, in that paper, the stiff order conditions for exponential Runge-Kutta methods (applied to semilinear parabolic problems) of orders up to 4 have been derived. The corresponding conditions for exponential Rosenbrock methods (making a continuous linearization of the vector field $F(t, U)$ along with the numerical solution in each step) were published in Hochbruck, Osterman, and Schweitzer (2009–paper in SINUM). However, open and important questions for exponential integrators remain: how to construct and analyze higher-order stiffly accurate methods (existing methods are up to order 4 only) or even methods of arbitrary order? How to construct parallel methods to take advantage of recent developments of numerical linear algebra in computing matrix functions or employing parallel computer architectures? How to design customized methods which can utilize preconditioners in a way that makes them clearly efficient in comparison with preconditioned implicit Newton-Krylov methods? How to build an efficient software package for the implementation of exponential integrators?

My main contributions to the field of exponential integrators:

1. Construction, analysis, and implementation of high-order exponential integrators:
In my Ph.D. research,  I focused on the construction, analysis, and implementation of high-order exponential integrators. In particular, with the help of a perturbation approach, I proposed a new and simpler approach for deriving order conditions for two popular classes of exponential integrators, namely exponential Runge-Kutta and exponential Rosenbrock methods. This allows me to give new stiff order conditions for exponential Runge-Kutta method of order up to 5 and exponential Rosenbrock methods of order up to 6 and thus construct methods of order 5. The error analysis is performed in an abstract Banach space framework of strongly continuous (or analytic) semigroups. Convergence results are proven independently of the stiffness of the problem. The implementation of new integrators and the comparison with existing methods are also carried out. See Luan (2014-Uni.Innsbruck); Luan and Ostermann (2014a-JCAM, 2014b-JCAM, 2014c HPSC12).

2. Derivation of a stiff order conditions theory for exponential integrators, “Exponential B-series“:
As mentioned above, the big and open question is how to derive stiff order conditions for exponential integrators of arbitrary order, just like explicit Runge–Kutta methods? In joint work with Prof. Ostermann, I pursued this challenging question and finally worked out a novel graph-theoretical approach to derive the order conditions for stiff problems by extending the well-known concept of B-series to exponential integrators. More precisely, we derived a stiff order condition theory based on the variation-of-constants formula. With the help of this theory, the stiff order conditions for such exponential methods of arbitrary order can be obtained in a simple and systematic way from a set of recursively defined trees. See Luan and Ostermann (2013-SINUM).

Stiff order conditions for exponential Runge-Kutta methods

Stiff order conditions for exponential Rosenbrock methods

3. Construction of parallel exponential time integration methods:
In another joint work with Prof. Ostermann, I have constructed parallel schemes based on exponential Rosenbrock methods. This was motivated by the fact that the existing exponential Rosenbrock schemes can be implemented only on serial computers (or on single processors). Our idea for constructing parallel exponential methods is to investigate the sparsity structure of the coefficient matrix of the method. In particular, we were able to construct new exponential Rosenbrock integrators of orders 4, 5, and 6 (e.g. see the left Figure below), which can be implemented on a multi-processor system or parallel computers. The new schemes of orders 4 and 5 require the same number of stages as the old schemes of the same orders of accuracy. However, while the parallel integrator of order 4 can be implemented with the same cost as a 2-stage method, the ones of orders 5 and 6 can be implemented at the cost of a 3-stage method only. This offers a significant improvement over the old schemes in terms of computational time when implemented in parallel. See Luan and Ostermann (2016-CAMWA)

4. Construction of preconditioned implicit-exponential methods:
I have constructed two new classes of time integrators for stiff PDEs: the implicit-explicit exponential (IMEXP) and the hybrid exponential methods. In contrast to the existing exponential schemes, our new methods can be used with any preconditioners and thus offer significant computational advantages. The results of our numerical simulations, in collaboration with colleagues from UCM, have verified the theoretical predictions of the performance of the new IMEXP schemes and presented them as a promising alternative to the widely used IMEX integrators for problems with very stiff nonlinearities (e.g. see the right Figure below). This leads to broader applicability of exponential methods, particularly in the context of large-scale applications such as fluid dynamics, reaction-diffusion systems, plasma physics, combustion, magnetohydrodynamic (MHD), etc. See Luan, Tokman, and Rainwater (graduate student) (paper 2017-JCP)

5. Construction of superconvergent exponential methods for time-dependent PDEs:
Among the family of fourth-order time integration schemes, the two-stage Gauss-Legendre method, which is an implicit Runge-Kutta method based on collocation, is the only superconvergent one. The computational cost of this implicit scheme for large systems, however, is very high since it requires solving a nonlinear system at every step.  In this work, I showed that one can construct and prove convergence results for exponential methods of order four which use two stages only. Specifically, I derived two new superconvergence (two-stage, fourth-order), fully explicit exponential Rosenbrock schemes for the time discretization of PDEs. It is shown that these new schemes offer great advantages over the two-stage Gauss-Legendre method as well as other fourth-order time integration schemes and stiff solvers such as ode15s (e.g. see the middle Figure below). See Luan (2017-APNUM)

An illustration of the performance of new exponential integrators vs. the standard methods such as ode15s and sBDF2 when applied to the 2D advection-diffusion-reaction equation and the Schnakenberg model from mathematical biology.

6. Construction and implementation of new efficient exponential Runge–Kuta methods:
Exponential Runge–Kutta methods have shown to be competitive for the time integration of stiff semilinear parabolic PDEs.The current construction of stiffly accurate exponential Runge–Kutta methods, however, relies on a convergence result that requires weakening many of the order conditions, resulting in schemes whose stages must be implemented in a sequential way. In this work, after showing a stronger convergence result, I derived two new families of fourth- and fifth-order exponential Runge–Kutta methods, which, in contrast to the existing methods, have multiple stages that are independent of one another and share the same format, thereby allowing them to be implemented in parallel or simultaneously, and making the methods to behave like using with much fewer stages. Moreover, all of their stages involve only one linear combination of the product of $\varphi$-functions (using the same argument) with vectors. Overall, these features make these new methods to be much more efficient to implement when compared to the existing methods of the same orders. Numerical experiments on a one-dimensional semilinear parabolic problem, a nonlinear Schrödinger equation, and a two-dimensional Gray–Scott model are given to confirm the accuracy and efficiency of the two newly constructed methods. See Luan (2021-BIT, arXiv version).

Efficiency comparison of new exponential RK methods vs. the existing methods of the same orders on a nonlinear Schr¨odinger equation.