Geometric Exponential Integrators

TL;DR

This paper introduces geometric exponential integrators for semilinear Poisson systems that preserve structure or energy, offering improved long-term stability and computational efficiency over traditional methods.

Contribution

The paper constructs two new types of exponential integrators that preserve either Poisson structure or energy, enhancing stability and efficiency for Hamiltonian PDE discretizations.

Findings

Geometric exponential integrators outperform non-geometric ones in long-term stability.

They are more computationally efficient than symplectic and discrete gradient methods.

Numerical experiments confirm their effectiveness for semilinear Poisson systems.

Abstract

In this paper, we consider exponential integrators for semilinear Poisson systems. Two types of exponential integrators are constructed, one preserves the Poisson structure, and the other preserves energy. Numerical experiments for semilinear Possion systems obtained by semi-discretizing Hamiltonian PDEs are presented. These geometric exponential integrators exhibit better long time stability properties as compared to non-geometric integrators, and are computationally more efficient than traditional symplectic integrators and energy-preserving methods based on the discrete gradient method.

Figures

Tables

Table 1. Table 1 . DISRK

$c_1$	$\frac{b_1}{2}$	0	0	0	0
$c_2$	$b_1$	$\frac{b_2}{2}$	0	0	0
$c_3$	$b_1$	$b_2$	$\frac{b_3}{2}$	0	0
$⋮$	$⋮$	$⋮$	$⋮$	$\ddots$	$⋮$
$c_s$	$b_1$	$b_2$	$b_3$	$\cdots$	$\frac{b_s}{2}$
	$b_1$	$b_2$	$b_3$	$\cdots$	$b_s$

Table 2. Table 2 . DISEX

$e^{Ah{c}_1}$	$\frac{b_1}{2}$	0	0	0	0
$e^{Ah{c}_2}$	$b_1{e}^{Ah\left(c_2-c_1\right)}$	$\frac{b_2}{2}$	0	0	0
$e^{Ah{c}_3}$	$b_1{e}^{Ah\left(c_3-c_1\right)}$	$b_2{e}^{Ah\left(c_3-c_2\right)}$	$\frac{b_3}{2}$	0	0
$⋮$	$⋮$	$⋮$	$⋮$	$\ddots$	$⋮$
$e^{Ah{c}_s}$	$b_1{e}^{Ah\left(c_s-c_1\right)}$	$b_2{e}^{Ah\left(c_s-c_2\right)}$	$b_3{e}^{Ah\left(c_s-c_3\right)}$	$\cdots$	$\frac{b_s}{2}$
$e^{Ah}$	$b_1{e}^{Ah\left(1-c_1\right)}$	$b_2{e}^{Ah\left(1-c_2\right)}$	$b_3{e}^{Ah\left(1-c_3\right)}$	$\cdots$	$b_s{e}^{Ah\left(1-c_s\right)}$

Table 3. Table 3 . Maximum timestep for nonlinear solver to converge as a function of numerical integrator, nonlinear solver, and spatial resolution.

	midpoint		midpoint exp	discrete gradient	energy exp
$n$	fixed point	Newton	fixed point	fixed point	fixed point
11	0.02	0.1	0.1	0.02	0.1
21	0.01	0.1	0.08	0.01	0.1
41	$4\times {10}^{-3}$	0.1	0.06	$4\times {10}^{-3}$	0.1
61	$2\times {10}^{-3}$	0.1	0.04	$2\times {10}^{-3}$	0.1
81	${10}^{-3}$	0.1	0.04	${10}^{-3}$	0.1
121	$5\times {10}^{-4}$	0.1	0.04	$5\times {10}^{-4}$	0.1
161	$2\times {10}^{-4}$	0.1	0.01	$2\times {10}^{-4}$	0.1
201	$1\times {10}^{-4}$	0.1	$8\times {10}^{-3}$	$1\times {10}^{-4}$	0.1
401	$4\times {10}^{-5}$	0.1	$5\times {10}^{-3}$	$4\times {10}^{-5}$	0.1

Table 4. Table 4 . Maximum timestep for fixed point iteration to converge as a function of numerical integrator, and spatial resolution.

$n$	midpoint	midpoint exp	discrete gradient	energy exp
401	$4\times {10}^{-4}$	$8\times {10}^{-4}$	$4\times {10}^{-4}$	0.005
601	$1\times {10}^{-4}$	$6\times {10}^{-4}$	$1\times {10}^{-4}$	0.005
801	$6\times {10}^{-5}$	$6\times {10}^{-4}$	$5\times {10}^{-5}$	0.005
1001	$3\times {10}^{-5}$	$6\times {10}^{-4}$	$3\times {10}^{-5}$	0.005
1201	$1\times {10}^{-5}$	$6\times {10}^{-4}$	$1\times {10}^{-5}$	0.005
1401	$1\times {10}^{-5}$	$6\times {10}^{-4}$	$1\times {10}^{-5}$	0.005

Equations

$$\dot{q}=Aq+f(q).$$

$$[e^{-At}q(t){]}^{\prime }={\tilde{q}}^{{}^{\prime }}(t)=e^{-At}f(e^{At}\tilde{q}(t)).$$

$$q(t)=e^{A(t-t_0)}q(t_0)+{\int }_{t_0}^t{e}^{A(t-\tau )}f(q(\tau ))d\tau ,$$

$$q_{k+1}=e^{Ah}{q}_k+{\int }_0^h{e}^{A\tau }d\tau \cdot f(q_k).$$

$$\dot{q}=J(Dq+\nabla V(q))=J\nabla H(q),$$

$$\frac{\tilde{q}_{k+1}-\tilde{q}_{k}}{h}=e^{-At_{k+1/2}}f\Big{(}e^{At_{k+1/2}}\frac{\tilde{q}_{k+1}+\tilde{q}_{k}}{2}\Big{)},$$

$$t_{k+1/2}=\frac{t_k+t_{k+1}}{2},h=t_{k+1}-t_k,{\tilde{q}}_k=e^{-A{t}_k}{q}_k,{\tilde{q}}_{k+1}=e^{-A{t}_{k+1}}{q}_{k+1}.$$

$$q_{k+1}=e^{Ah}q_{k}+h\cdot e^{A\frac{h}{2}}f\Big{(}\frac{e^{A\frac{h}{2}}q_{k}+e^{-A\frac{h}{2}}q_{k+1}}{2}\Big{)}.$$

$$(\nabla \varphi )J(\nabla \varphi {)}^{\mathrm{T}}=J.$$

$$\displaystyle dq_{k+1}=e^{Ah}dq_{k}+h\cdot e^{A\frac{h}{2}}\nabla f\Big{(}\frac{1}{2}e^{A\frac{h}{2}}dq_{k}+\frac{1}{2}e^{-A\frac{h}{2}}dq_{k+1}\Big{)},$$

$$\displaystyle\Big{(}I-\frac{h}{2}\cdot e^{A\frac{h}{2}}\nabla fe^{-A\frac{h}{2}}\Big{)}dq_{k+1}=\Big{(}e^{Ah}+\frac{h}{2}\cdot e^{A\frac{h}{2}}\nabla fe^{A\frac{h}{2}}\Big{)}dq_{k}.$$

$$M$$

$$N$$

$$MJ{M}^{\mathrm{T}}$$

$$=J-\frac{h^2}{4}{e}^{A\frac{h}{2}}J{\nabla }^{2}V{e}^{-A\frac{h}{2}}J{e}^{A\frac{h}{2}}{\nabla }^{2}VJ{e}^{-A\frac{h}{2}}$$

$$=J-\frac{h^2}{4}{e}^{A\frac{h}{2}}J{\nabla }^{2}VJ{\nabla }^{2}VJ{e}^{-A\frac{h}{2}}$$

$$=NJ{N}^{\mathrm{T}}.$$

$$\{\begin{array}{rl}{Q}_i& =e^{Ah{c}_i}{q}_k+h\cdot \sum _{j=1}^i{a}_{ij}{e}^{Ah(c_i-c_j)}f(Q_j),\\ q_{k+1}& =e^{Ah}{q}_k+h\cdot \sum _{i=1}^s{b}_i{e}^{Ah(1-c_i)}f(Q_i).\end{array}$$

$$q_k\stackrel{b_{1}h}{\to }{Z}_1\stackrel{b_{2}h}{\to }{Z}_2\cdots \stackrel{b_{s}h}{\to }{Z}_s=q_{k+1}$$

$$Z_1$$

$$Z_2$$

$$⋮$$

$$Z_i$$

$$⋮$$

$$Z_s$$

$$Q_1=e^{Ah\frac{b_1}{2}}{q}_k+\frac{b_1}{2}h\cdot f(Q_1).$$

$$\{\begin{array}{rl}{c}_1& =\frac{b_1}{2},\\ c_2& =b_1+\frac{b_2}{2},\\ & ⋮\\ c_i& =b_1+b_2+\cdots b_{i-1}+\frac{b_i}{2},\\ & ⋮\\ c_s& =b_1+b_2+b_3+\cdots +b_{s-1}+\frac{b_s}{2},\\ 1& =b_1+b_2+b_3+\cdots +b_{s-1}+b_s.\end{array}$$

$$Q_2$$

$$=e^{Ah(b_1+\frac{b_2}{2})}{q}_k+b_{1}h\cdot e^{Ah(\frac{b_1}{2}+\frac{b_2}{2})}f(Q_1)+\frac{b_2}{2}h\cdot f(Q_2)$$

$$=e^{Ah{c}_2}{q}_k+b_{1}h\cdot e^{Ah(c_2-c_1)}f(Q_1)+\frac{b_2}{2}h\cdot f(Q_2).$$

$$Q_i=e^{Ah\frac{b_i}{2}}{Z}_{i-1}+\frac{b_i}{2}h\cdot f(Q_i).$$

$$Z_{i-1}$$

$$+b_{2}h\cdot e^{Ah(b_{i-1}+b_{i-2}+\cdots +\frac{b_2}{2})}f(Q_2)+\cdots b_{i-1}h\cdot e^{Ah(\frac{b_{i-1}}{2})}f(Q_{i-1}),$$

$$Q_i$$

$$+b_{2}h\cdot e^{Ah(\frac{b_i}{2}+b_{i-1}+b_{i-2}+\cdots +\frac{b_2}{2})}f(Q_2)+\cdots$$

$$+b_{i-1}h\cdot e^{Ah(\frac{b_i}{2}+\frac{b_{i-1}}{2})}f(Q_{i-1})+\frac{b_i}{2}h\cdot f(Q_i)$$

$$=e^{Ah{c}_i}{q}_k+b_{1}h\cdot e^{Ah(c_i-c_1)}f(Q_1)+\cdots$$

$$+b_{i-1}h\cdot e^{Ah(c_i-c_{i-1})}f(Q_{i-1})+\frac{b_i}{2}h\cdot f(Q_i),$$

$$q_{k+1}$$

Full text

1. Introduction

Exponential integrators [4] are a class of numerical integrators for stiff systems whose vector field can be decomposed into a linear term and a nonlinear term,

[TABLE]

Usually, the coefficient matrix $A$ has a large spectral radius, and is responsible for the stiffness of the system of differential equations, while the nonlinear term $f(q)$ is relatively smooth. There are various ways to construct an exponential integrator [7]. For example, we can perform a change of variables $\tilde{q}(t)=e^{-At}q(t)$, and transform (1) to obtain

[TABLE]

Notice that the Jacobi matrix of (2) equals $e^{-At}\nabla fe^{At}$, which has a smaller spectral radius than the Jacobi matrix $A+\nabla f$ of (1). A natural idea is to apply a classical integrator for the mollified system (2) to obtain an approximation of $\tilde{q}(t)$, then invert the change of variables to obtain an approximation of the solution $q(t)$ of (1). In Section 2, we shall demonstrate how to construct symplectic exponential integrators using this approach.

Another way of constructing exponential integrators starts from the variation-of-constants formula,

[TABLE]

which is the exact solution for (1) with initial condition $q(t_{0})=q_{0}$. Then, a computable approximation can be obtained by approximating the $f(q(\tau))$ term inside the integral. If we approximate $f(q(\tau))$ by $f(q_{k})$, we arrive at the exponential Euler method,

[TABLE]

An exponential Runge–Kutta method of collocation type [3] could also be constructed by approximating $f(q(\tau))$ with polynomials. In Section 3, we shall show how to construct energy-preserving exponential integrators from (3).

In this paper, we consider a specific form of (1) which is a Poisson system. We assume $A=JD$, $f(q)=J\nabla V(q)$, where $J^{\mathrm{T}}=-J,D^{\mathrm{T}}=D$, and $JD=DJ$, thus the coefficient matrix $A$ is also skew-symmetric. Now, the semilinear system (1) can be written as,

[TABLE]

with Hamiltonian function $H(q)=\frac{1}{2}q^{\mathrm{T}}Dq+V(q)$. Equation (5) describes a constant Poisson system, and there are at least two quantities that are preserved by the flow: the Poisson structure $J_{ij}\frac{\partial}{\partial x_{i}}\otimes\frac{\partial}{\partial x_{j}}$ and Hamiltonian $H(q)$. Geometric integrators that preserve the geometric structure and first integrals of the system typically exhibit superior qualitative properties when compared to non-geometric integrators, and they are an active area of research [2, 6, 5].

Here, we construct geometric exponential integrators that either preserve the Poisson structure or Hamiltonian of (5). They exhibit long time stability, allow for relatively larger timesteps for the stiff problem, and are computationally more efficient as they can be implemented using fixed point iterations as opposed to Newton type iterations. For the rest of the paper, Section 2 is devoted to developing symplectic exponential integrators that preserve the Poisson structure; Section 3 is devoted to developing energy preserving exponential integrators; numerical methods and experiments are presented in Section 4 and Section 5, respectively.

2. Symplectic Exponential Integrator

For constant Poisson systems (5), it was shown in [12] that the midpoint rule and diagonally implicit symplectic Runge–Kutta methods preserve the Poisson structure $J_{ij}\frac{\partial}{\partial x_{i}}\otimes\frac{\partial}{\partial x_{j}}$. We first start by constructing an exponential midpoint rule: apply the classical midpoint rule to the transformed system (2) to obtain

[TABLE]

where

[TABLE]

Transform (6) back to $q_{k}$ and $q_{k+1}$, and we obtain the exponential midpoint rule,

[TABLE]

Theorem 1**.**

The exponential midpoint rule (7) preserves the Poisson structure when applied to the semilinear Poisson system (5).

Proof.

Recall that a map $\phi$ preserves Poisson structure $J_{ij}\frac{\partial}{\partial x_{i}}\otimes\frac{\partial}{\partial x_{j}}$ iff

[TABLE]

Differentiating (7), we obtain,

[TABLE]

So the map $\phi(q_{k})=q_{k+1}$ has Jacobi matrix $\nabla\phi=M^{-1}N$, where

[TABLE]

Then, we just need to verify (8), which is equivalent to $MJM^{\mathrm{T}}=NJN^{\mathrm{T}}$,

[TABLE]

In (9), we used the property that $\nabla^{2}V$ is symmetric, that $A$ is skew-symmetric which implies that $(e^{Ah})^{\mathrm{T}}=e^{-Ah}$, and the assumptions that $JD=DJ$ and that $e^{A\frac{h}{2}}$ and $J$ commute. ∎

The exponential midpoint rule is a second-order method, and we will now develop higher-order symplectic exponential integrators. Recall that a general diagonally implicit Symplectic Runge–Kutta method (DISRK) has a Butcher tableau of the following form,

It was shown in [12] that any DISRK method can be regarded as the composition of midpoint rules with timesteps $b_{1}h,b_{2}h,b_{3}h,\cdots b_{s}h$. If we apply the DISRK method for the transformed system (2), and then convert back, we obtain the following diagonally implicit symplectic exponential (DISEX) integrator,

[TABLE]

where $a_{ij}$ are the coefficients in Table 1. This integrator can be represented in terms of the following Butcher tableau,

The DISEX method preserves the Poisson structure, as demonstrated by the following theorem.

Theorem 2**.**

The DISEX method is equivalent to the composition of exponential midpoint rules with timesteps $b_{1}h,b_{2}h,b_{3}h,\cdots b_{s}h$.

Proof.

The composition of exponential midpoint rules

[TABLE]

is represented as follows,

[TABLE]

Introduce $Q_{1}=\frac{e^{Ah\frac{b_{1}}{2}}q_{k}+e^{-Ah\frac{b_{1}}{2}}Z_{1}}{2}$ in (M.1) as an intermediate variable. Then, on both sides of (M.1), multiply by $e^{-Ah\frac{b_{1}}{2}}$, add $e^{Ah\frac{b_{1}}{2}}q_{k}$, and divide by two, which yields an equivalent form of (M.1),

[TABLE]

For the Runge–Kutta method represented by the Butcher tableau in Table 1 to be consistent, the coefficients have to satisfy,

[TABLE]

So equation (S.2) coincides with the first line of the Butcher tableau of the DISEX method. Similarly, introduce $Q_{2}=\frac{e^{Ah\frac{b_{2}}{2}}Z_{1}+e^{-Ah\frac{b_{2}}{2}}Z_{2}}{2}$. Then, on both sides of (M.2), multiply by $e^{-Ah\frac{b_{2}}{2}}$, add $e^{Ah\frac{b_{2}}{2}}Z_{1}$, then divided by two, which yields an equivalent form of (M.2):

[TABLE]

So equation (S.2) coincides with the second line of the Butcher tableau of the DISEX method. Then, as before, we introduce $Q_{i}=\frac{e^{Ah\frac{b_{i}}{2}}Z_{i-1}+e^{-Ah\frac{b_{i}}{2}}Z_{i}}{2}$, and apply the same technique to (M.i), which yields

[TABLE]

By induction,

[TABLE]

so

[TABLE]

which coincides with the $i$-th row of the Butcher tableau of the DISEX method. Finally, we have

[TABLE]

which coincides with the last row of the Butcher tableau of the DISEX method. So the composition of exponential midpoint rules with timesteps $b_{1}h,b_{2}h,b_{3}h,\cdots b_{s}h$ is equivalent to the DISEX method of Table 2. ∎

3. Energy-preserving Exponential Integrator

Though classical symplectic methods exhibit superior long time stability, it was observed that symplectic schemes are less competitive for the numerical integration of stiff systems with high frequency. In sharp contrast, energy-preserving methods perform much better [9]. A general way to construct an energy-preserving method for a Poisson system $\dot{q}=J\nabla H(q)$ is the discrete gradient method [8]. We design a discrete gradient $\overline{\nabla}H(q_{k},q_{k+1})$ that satisfies the following property,

[TABLE]

Then, the resulting discrete gradient method is given by,

[TABLE]

Multiplying $\overline{\nabla}H(q_{k},q_{k+1})$ on both sides of (12), we obtain

[TABLE]

The last equation of (13) holds simply due to the skew-symmetric property of matrix $J$, which implies that discrete gradient method (12) preserves energy. We shall combine exponential integrators with the discrete gradient method to obtain an energy-preserving exponential integrator. This approach was initially proposed in [11] for separable Hamiltonian systems using the extended discrete gradient method, and we generalize this to semilinear Poisson systems. Replace the $f(q_{k})$ term in the exponential Euler method (4) by the discrete gradient $J\overline{\nabla}V(q_{k},q_{k+1})$, which yields

[TABLE]

Theorem 3**.**

Method (14) preserves the Hamiltonian $H(q)$.

Proof.

Let $S=e^{Ah}$, $T=\int_{0}^{h}e^{A\tau}d\tau$, then $q_{k+1}=Sq_{k}+TJ\overline{\nabla}V(q_{k},q_{k+1})$, and we will show that the following properties hold:

1. (1)

$S^{\mathrm{T}}=S^{-1}$; 2. (2)

$AT=S-I$; 3. (3)

$AT^{\mathrm{T}}=I-S^{\mathrm{T}}$; 4. (4)

$S^{\mathrm{T}}T=T^{\mathrm{T}}$.

Property 1 follows from the fact that $S^{\mathrm{T}}=(e^{Ah})^{\mathrm{T}}=e^{-Ah}=S^{-1}$. Property 2 follows from

[TABLE]

Taking transposes on both sides of Property 2 and using the fact that $A$ and $T$ commute gives Property 3. Property 4 follow from

[TABLE]

Recall that the Hamiltonian of the constant Poisson system (5) is $H(q)=\frac{1}{2}q^{\mathrm{T}}Dq+V(q)$, so

[TABLE]

Applying the properties above yields the following,

[TABLE]

Substituting these into the expression for $H(q_{k+1})$ yields

[TABLE]

which proves that the Hamiltonian is preserved. ∎

One significant advantage of exponential integrators is they allow the numerical method to be implemented using fixed point iterations as opposed to the more computationally expensive Newton iterations. Recall that classical implicit Runge-Kutta methods

[TABLE]

have a form that naturally lends itself to fixed point iterations. However, when $\frac{\partial f}{\partial y}$ has a large spectral radius, the timestep is forced to be very small in order to guarantee that the fixed point iteration converges. The alternative is to use a Newton type iteration, which is time consuming since we need to perform LU decomposition ($O(n^{3})$ complexity) during each iteration. This is the problem we face for the stiff semilinear system (1) when the coefficient matrix $A$ has a large spectral radius. In contrast, in both the exponential midpoint rule

[TABLE]

and the energy-preserving exponential integrator

[TABLE]

the matrix $A$ only appears in the exponential term $e^{Ah}$. Since $A$ is skew-symmetric, this term is an orthogonal matrix, which has spectral radius 1. Thus, fixed point iterations can be used to implement the exponential integrator, regardless of the stiffness of $A$.

4. Numerical Methods

We consider two Hamilton PDEs here, namely the nonlinear Schrödinger equation,

[TABLE]

in which $\psi=u+iv$ is the wave function with real part $u$ and imaginary part $v$, and has the following equivalent form,

[TABLE]

as well as the KdV equation,

[TABLE]

To discretize the two PDEs, we impose $2\pi$ periodic boundary conditions. Given a smooth $2\pi$ periodic function $f(x)$, on the interval $[0,2\pi]$, choose $2n+1$ equispaced interpolation points $x_{j}=jh,j=0,1,2,\ldots,2n$, $h=\frac{2\pi}{2n+1}$. Given nodal values $\{v_{j}\}_{j=0}^{2n}$, there exists a unique trigonometric polynomial $v(x)$ with degree less or equal $n$, such that, $v(x_{j})=v_{j}$ (see, for example, [1]).

[TABLE]

By substituting the expression for the coefficients $\hat{v}_{k}$, we obtain

[TABLE]

where

[TABLE]

From this, we see that $\{e^{ikx}\}^{n}_{k=-n}$ and $\{\phi_{j}\}^{2n}_{j=0}$ are equivalent orthogonal bases for the trigonometric polynomial function space, and each such function can be parametrized by either the nodal values $\{v_{j}\}^{2n}_{j=0}$ or Fourier coefficients $\{\hat{v}_{k}\}^{n}_{k=-n}$. They represent the same function, but with respect to two different bases. The transformation between $\{v_{j}\}^{2n}_{j=0}$ and $\{\hat{v}_{k}\}^{n}_{k=-n}$ can be performed using the Fast Fourier transformation (FFT), which has $O(n\log n)$ complexity.

The first and second-order differentiation matrices [10] with respect to the representation in terms of nodal values $\{v_{j}\}^{2n}_{j=0}$ are given by

[TABLE]

respectively. However, with respect to the representation in terms of Fourier coefficients $\{\hat{v}_{k}\}^{n}_{k=-n}$, they are diagonal,

[TABLE]

Later, we will see that this observation is critical to a fast implementation of the product of matrix functions with vectors. We can also define a third-order differentiation matrix $D_{3}$, which has the property $D_{3}=D_{1}D_{2}=D_{2}D_{1}$, and it is diagonal with respect to the Fourier coefficients $\hat{D}_{3}=\operatorname{diag}(-ik^{3})^{n}_{k=-n}$.

4.1. Nonlinear Schrödinger equation

We perform a semi-discretization of (18) by discretizing the solution $u$, $v$ in space using their corresponding nodal values $\{q_{j}\}$ and $\{p_{j}\}$. Applying the pseudospectral method, we obtain the following system of ODEs,

[TABLE]

where the nonlinear term $(q^{2}+p^{2})p$ is computed elementwise, and represents the vector consisting of $\{(q^{2}_{j}+p^{2}_{j})p_{j}\}$ entries. We adopt this notation throughout the rest of the paper for brevity. Then, (21) can be expressed as,

[TABLE]

where

[TABLE]

[TABLE]

and $V(q,p)=\frac{1}{2}(q^{2}+p^{2})^{2}$. It is easy to verify that $J$ is skew-symmetric, $D$ is symmetric, and $JD=DJ$. Thus, (22) is a semilinear Poisson system.

To apply the exponential midpoint rule, we need to compute the product of a matrix function and a vector, which has the form $e^{Ah}\begin{pmatrix}q\\ p\end{pmatrix}$,

[TABLE]

Theorem 4**.**

Suppose that $D_{2}$, $D_{3}$ are the second and third-order differentiation matrices, respectively, $q=\{q_{j}\}_{j=0}^{2n}$ is a vector with Fourier transform $F[q]=\hat{q}=\{\hat{q}_{k}\}_{k=-n}^{n}$, and $f$ is an analytic function, then

[TABLE]

where $F^{-1}$ is the inverse Fourier transform.

Proof.

Recall that matrix $D_{2}$ is diagonalizable with eigenvalues $\lambda_{k}=-k^{2}$, and corresponding eigenvectors $e_{k}=\{e^{ikx_{j}}\}_{j=0}^{2n}$,

[TABLE]

Notice $D_{3}$ is also diagonalizable with eigenvalues $\lambda_{k}=-ik^{3}$, and corresponding eigenvectors $e_{k}$, so the property that $f(D_{3})q=F^{-1}[\operatorname{diag}(f(-ik^{3}))\hat{q}]$ can be verified in the same way. ∎

By Theorem 4, we have that

[TABLE]

In summary, the exponential midpoint rule for the nonlinear Schrödinger equation is given by

[TABLE]

where $z_{k}=\begin{pmatrix}q_{k}\\ p_{k}\end{pmatrix}$, $z_{k+1}=\begin{pmatrix}q_{k+1}\\ p_{k+1}\end{pmatrix}$, $A=\begin{pmatrix}0&-D_{2}\\ D_{2}&0\end{pmatrix}$, and $e^{Ah}z_{k}$ can be efficiently calculated using (23).

For the energy-preserving exponential integrator for the nonlinear Schrödinger equation,

[TABLE]

Here,

[TABLE]

and we can construct the discrete gradient

[TABLE]

where

[TABLE]

Notice that $\overline{\nabla}V(z_{k},z_{k+1})$ is symmetric with respect to $z_{k}$ and $z_{k+1}$, and it can be verified that it satisfies (11). A classical discrete gradient method can be constructed as follows,

[TABLE]

where

[TABLE]

The method described by (26) is very similar to classical midpoint rule, the only difference is in $\overline{\nabla}V(z_{k},z_{k+1})$, ${q^{2}}_{k+\frac{1}{2}}$ is used, while $(q_{k+\frac{1}{2}})^{2}$ is used in the midpoint rule, so (26) can be viewed as a modified midpoint rule.

4.2. KdV equation

Rewrite (19) as

[TABLE]

then apply pseudospectral semi-discretization to obtain the following system,

[TABLE]

which has the form of a semilinear Poisson system (5),

[TABLE]

where $J=-D_{1}$, $D=D_{2}$, $A=JD=-D_{3}$, $\nabla V(q)=\frac{1}{2}q^{2}$, $H(q)=\frac{1}{2}q^{\mathrm{T}}D_{2}q+\frac{1}{6}q^{3}$. The exponential midpoint rule for KdV reads as follows,

[TABLE]

and the energy preserving exponential integrator is given by,

[TABLE]

with discrete gradient $\overline{\nabla}V(q_{k},q_{k+1})=\frac{1}{6}(q_{k}^{2}+q_{k}\cdot q_{k+1}+q_{k+1}^{2})$. A related classical discrete gradient method can be constructed as follows,

[TABLE]

By Theorem 4, in each iteration, the matrix function and vector product can be implemented as

[TABLE]

and

[TABLE]

5. Numerical Experiments

5.1. Nonlinear Schrödinger equation

In Table 3 below, $n$ denotes the number of nodes we discretize the spatial domain with, and the data in the table indicates the maximum timesteps for which the nonlinear solver converges. The first two columns correspond to midpoint rule, using fixed point iteration and Newton type iteration; the third column corresponds to the exponential midpoint rule, the fourth column is discrete gradient method (26), and the last column is the energy-preserving exponential integrator, all implemented using fixed point iterations.

We observe that when the midpoint rule is implemented using fixed point iterations, there is a significant decay in the allowable timestep as the spatial resolution is increased, whereas the Newton type iteration allows a relatively large timestep that is independent of the spatial resolution. In contrast, the midpoint exponential method exhibits a slower rate of decrease in allowable timestep when using fixed point iterations. When using fixed point iterations, the discrete gradient method, which is an energy preserving method, the allowable timestep behaves similarly to the midpoint rule, and in contrast, the energy preserving exponential integrator has an allowable timestep that is independent of the spatial resolution.

In Figure 1, we observe that the allowable timestep when using fixed point iteration scale like $n^{-2}$ for the classical midpoint rule, and $n^{-1}$ for the midpoint exponential rule.

As shown in Figure 2, the exponential midpoint rule exhibits an energy error that remains small and bounded, which is consistent with it being a symplectic integrator, and the trajectory error grows linearly.

The energy preserving exponential integrator is designed to preserve energy exactly, so even when the timestep is $0.1$, we see in Figure 3 that the energy is still preserved approximately to within machine error.

We now consider the diagonally implicit symplectic exponential (DISEX) integrator with six stages:

[TABLE]

The energy and trajectory error is shown in Figure 4. Observe that the energy error is small and bounded, as expected of a symplectic integrator, and the trajectory error is small as well, as expected of a higher-order numerical integrator.

5.2. KdV

We simulate the KdV equation,

[TABLE]

where $\nu=5\times 5^{-4}$. As before, we explore how the maximum timestep for which the fixed point iteration converges depends on the choice of numerical integrator, and the spatial resolution of the semi-discretization. In Table 4, $n$ is the number of nodes we use to discretize the spatial domain, the data in the table indicates the maximum timestep for which fixed point iterations converge. For the midpoint rule, the timestep decreases like $n^{-3}$ for fixed point iterations, and a comparable timestep is required for Newton iterations which is thereby too costly to implement. The classical discrete gradient method for the KdV equation is given by (27), and it exhibits the same timestep restrictions as the midpoint rule. In contrast, both the exponential midpoint and energy preserving exponential integrator allow rather large timesteps that are independent of the spatial resolution.

In Figure 5, we observe that the exponential midpoint rule has an energy error that is small and bounded, as is typical for a symplectic integrator, and the trajectory error grows linearly. In Figure 6, the energy preserving exponential integrator has an energy that is preserved to within machine precision, and the trajectory error grows linearly.

6. Summary

For semilinear Poisson system, we have developed two types of geometric exponential integrators, one that preserves the Poisson structure, and the other preserves the energy. They allow fixed point iteration methods to be used for significantly larger timesteps as compared to non-exponential integrators, such as the classical midpoint rule and the standard discrete gradient method, which require the use of Newton type iterations to converge with comparably large timesteps. This results in substantial computational savings, particularly when applied to the simulation of semi-discretized Hamiltonian PDEs. They also exhibit long time stability and conservation of the first integrals. We notice that in practice, energy preserving exponential integrators are more stable and allow for larger timesteps than symplectic exponential integrators. In future work, we will develop higher-order energy preserving exponential integrators.

Acknowledgements

This research has been supported in part by NSF under grants DMS-1010687, CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, DMS-1345013.