Stability of the first-order exponential integrator method

Question

The question is about the first-order exponential integration method described in this article.

Consider a system of ordinary differential equations

$$y'(t) = -A\,y(t) + \mathcal{N}(t, y), \qquad y(0) = y_0,$$

where $\mathcal{N}$ is a nonlinear function. The solution can be approximated using the exponential (explicit) Euler method given by

$$y_{n+1} = e^{-Ah}y_n + A^{-1}(1-e^{-Ah}) \mathcal{N}(t_n, \, y(t_n)).$$

The question is: What can be said about the stability of the above integration scheme provided that A is a real, symmetric matrix? What if A is also positive semi-definite? I am aware that the ordinary (explicit) Euler scheme is known to be unstable, but what about its exponential analogue?

Thank you!

Regards, Ivan

Answer 1

It depends on your nonlinear function $\mathcal{N}.$ The stability means the solution continuously depends on the data. In particular, for discretised ODEs, it means if you give a perturbation to $y_n,$ the perturbation propagating to $y_{n+m}$ is bounded uniformly w.r.t. $m.$

In your case, we have $$ |\tilde{y}_{n+1}-y_{n+1}|\leq |e^{-Ah}|\cdot|\tilde{y}_n-y_n| + (h + \frac{1}{2}|A|h^2 + ..)|\mathcal{N}(t_n,\tilde{y}_n)-\mathcal{N}(t_n,{y}_n)|. $$ Assuming $\mathcal{N}$ is Lipschitz w.r.t. the second argument, you will need something like $$ |e^{-Ah}| + Ch = 1 + \mathcal{O}(h) $$ for asymptotic stability (stable as $h\rightarrow 0$) or if you want absolute stability $$ |e^{-Ah}| + Ch<1, $$ where $C$ depends on $|A|.$

Answer 2

Exponential time integrators are exact for linear systems of ordinary differential equations with constant coefficients and are thus trivially A-stable.

See [1] and the numerous references therein.

For a general system of ODEs, the basic idea is to split up the terms in such a way that the largest stable time step for the exponential integrator is significantly larger than for a corresponding explicit discretization, thus reducing the computer time required to solve the equations.

To achieve this goal without sacrificing accuracy, three issues must be addressed:

1. all terms required for the fastest solution components should be treated in the linear part $A$,
2. the energy in those fast components should be insignificant, and
3. the resulting matrix exponential should be relatively easy to solve.

Points 1 and 3 should be obvious. Point 2 is not always clearly stated, but it is necessary: contrary to implicit method, where the fastest waves are slowed, solution modes treated in the exponential part are not distorted by the time discretization, and thus an accurate solution can be obtained only if there is little energy in these modes. The slowest modes are usually treated in $N$

[1] Tokman, M. 2006. Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods. J. Comput. Phys. 213, 748–776.