High-order linearly implicit structure-preserving exponential integrators for the nonlinear Schrödinger equation

1 Introduction

We consider the nonlinear Schrödinger equation (NLSE) given as follows

{\begin{matrix} i \partial_{t} ϕ (x, t) + L ϕ (x, t) + β | ϕ (x, t) |^{2} ϕ (x, t) = 0, (x, t) \in Ω, t > 0, ϕ (x, 0) = ϕ_{0} (x), x \in Ω \subset R^{d}, d = 1, 2, 3, \end{matrix}

(1.1)

where $t$ is the time variable, $x$ is the spatial variable, $i = \sqrt{- 1}$ is the complex unit, $ϕ := ϕ (x, t)$ is the complex-valued wave function, $L$ is the usual Laplace operator, $β$ is a given real constant, and $ϕ_{0} (x)$ is a given initial data. The nonlinear Schrödinger equation (1.1) satisfies the following two invariants, that is mass

M (t) = \int_{Ω} | ϕ |^{2} d x \equiv M (0), t \geq 0,

(1.2)

and Hamiltonian energy

E (t) = \int_{Ω} (- | \nabla ϕ |^{2} + \frac{β}{} 2 | ϕ |^{4}) d x \equiv E (0), t \geq 0.

(1.3)

A numerical scheme that preserves one or more invariants of original systems is called an energy-preserving method. The earlier attempts to develop energy-preserving methods for NLSE can date back to 1981, when Delfour, Fortin, and Payre [22] constructed a two-level finite difference scheme (also called Crank-Nicolson finite difference (CNFD) scheme) which can satisfy the discrete analogue of conservation laws (1.2) and (1.3). Some comparisons with nonconservation schemes are investigated by Sanz-Serna and Verwer [49], which show the non-energy-preserving schemes may easily show nonlinear blow-up. Furthermore, it can be proven rigorously in mathematics that the CNFD scheme is second-order accurate in time and space [6, 7, 54]. In particular, we note that the discrete conservative laws play a crucial role in numerical analyses for numerical schemes of NLSE. However, it is fully implicit and at every time step, a nonlinear equation shall be solved by using a nonlinear iterative method and thus it may be time consuming. Zhang et al. [58] firstly proposed a linearly implicit Crank-Nicolson finite difference (LI-CNFD) scheme in which a linear system is to be solved at every time step. Thus it is computationally much cheaper than that of the CNFD scheme. The LI-CNFD scheme can satisfy a new discrete analogue of conservation laws (1.2) and (1.3) and its stability and convergence were analyzed in [54]. Up to now, the energy-preserving Crank-Nicolson schemes for NLSE are extensively extended and analyzed [3, 2, 17, 23, 26, 41, 56]. Another popular linearly implicit energy-preserving method is the relaxation finite difference scheme [9]

. In addition, the energy-preserving schemes for NLSE provided by the averaged vector field method

[47] and the discrete variational derivative method [21, 24] can be founded in [16, 38, 45, 25]. More recently, the scalar auxiliary variable (SAV) approach [51, 52] have been a particular powerful tool for the design of linearly implicit energy-preserving numerical schemes for NLSE [5]. Nevertheless, to our best knowledge, all of the schemes can achieve at most second-order in time.

It has been proven that no Runge-Kutta (RK) method can preserve arbitrary polynomial invariants of degree 3 or higher of arbitrary vector fields [14]. Thus, over that last decades, how to develop high-order energy-preserving methods for general conservation systems attracts a lot of attention. The notable ones include, but are not limited to, high-order averaged vector field (AVF) method [44, 47], Hamiltonian Boundary Value Methods (HBVMs) [11, 12] and energy-preserving continuous stage Runge-Kutta methods (CSRKs) [19, 29, 46, 55]. Actually, the HBVMs and CSRKs have been shown to be an efficient method to develop high-order energy-preserving schemes for NLSE [8, 43], however, the proposed schemes are fully implicit. In fact, at every time step, one needs to solve a large fully nonlinear system and thus it might lead to high computational costs.

Due to the high computational cost of the high-order fully implicit schemes, in the literatures, ones are devoted to construct energy-preserving explicit schemes for NLSE. Based on the the invariant energy quadratization (IEQ) approach [57] or the SAV approach and the projection method [13, 30, 28], the authors proposed some explicit energy-preserving schemes for NLSE [34, 59]. Such scheme is extremely easy to implement, however, it requires very small time step sizes. Especially, this limitation is more serious in 2D and 3D. The linearly implicit one which needs solve a linear equations at every time step, can remove this limitation. Li et al. proposed a class of linearly implicit schemes for NLSE, which can preserve the discrete analogue of conservation law (1.2) [42]. It involves solving linear equations with complicated variable coefficients at every time step and thus it may be very time consuming. Another way to achieve this goal is to combine the SAV approach with the extrapolation technique or the prediction-correction strategy [10].

In the past few decades, many energy-preserving exponential integrators have been done for conservative systems. The exponential integrator often involves exact integration of the linear part of the target equation, and thus it can achieve high accuracy, stability for a very stiff differential equation such as highly oscillatory ODEs and semidiscrete time-dependent PDEs. For more details on the exponential integrator, please refer to the excellent review article provided by Hochbruck and Ostermann [32]. Based on the projection approach, Celledoni et al. developed some symmetry- and energy-preserving implicit exponential schemes for the cubic Schrödinger equation [15]. By combining the exponential integrator with the AVF method, Li and Wu constructed a class of second-order implicit energy-preserving exponential schemes for canonical Hamiltonian systems and successfully applied it to solve NLES [44]. Further analysis and generalization is investigated by Shen and Leok [53]. More recently, Jiang et al. showed that the SAV approach is also an efficient approach to develop second-order linearly implicit energy-preserving exponential scheme for NLSE [33]. Overall, there exist very few works devoted to development of high-order linearly implicit exponential schemes with energy-preserving property for NLSE, which motivates this paper.

In this paper, following the idea of the SAV approach, we firstly reformulate the original system (1.1) into a new system by introducing a new auxiliary variable, which satisfies a modified energy. The spatial discretization is then performed with the standard Fourier pseudo-spectral method [18, 50]. Subsequently, the extrapolation technique is employed to the nonlinear term of the semi-discrete system and a linearized system is obtained. Based on the Lawson transformation, the semi-discrete system is rewritten as an equivalent form and the fully discrete scheme is further obtained by using the symplectic RK method in time. We show that the proposed scheme can preserve the modified energy in the discretized level and at every time step, only a linear equations with constant coefficients is solved for which there is no existing references [10, 40] considering this issue.

The rest of this paper is organized as follows. In Section 2, we use the idea of the SAV approach to reformulate the system (1.1) into a new reformulation, which satisfies a modified energy. In Section 3, the linearly implicit high-order energy-preserving exponential scheme is proposed and its energy-preservation is discussed. In Section 4, several numerical examples are investigated to illustrate the efficiency of the proposed scheme. We draw some conclusions in Section 5.

2 Model reformulation

In this section, we employ the SAV approach to recast the NLS equation (1.1) into a new reformulation which satisfies a quadratic energy conservation law. The resulting reformulation provides an elegant platform for developing high-order linearly implicit energy-preserving schemes.

Based on the idea of the SAV approach [51, 52], we introduce an auxiliary variable

p := p (t) = \sqrt{(ϕ^{2}, ϕ^{2}) + C_{0}},

where $C_{0}$ is a constant large enough to make $p$ well-defined for all $ϕ$ . Here, $(v, w)$ is the inner product defined by $\int_{Ω} v ¯ w d x$ where $¯ w$ represents the conjugate of $w$ . Then, the Hamiltonian energy functional is rewritten as the following quadratic form

E (t) = (L ϕ, ϕ) + \frac{β}{2} p^{2} - \frac{β}{2} C_{0} .

According to the energy variational principle, we obtain the following SAV reformulated system

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \partial_{t} ϕ = i (L ϕ + \frac{β | ϕ |^{2} ϕ p}{\sqrt{(ϕ^{2}, ϕ^{2}) + C_{0}}}), \frac{d}{d t} p = 2 R e (i L ϕ, \frac{| ϕ |^{2} ϕ}{\sqrt{(ϕ^{2}, ϕ^{2}) + C_{0}}}), \end{matrix}

(2.1)

with the consistent initial condition

ϕ (x, 0) = ϕ_{0} (x), p (0) = \sqrt{(ϕ_{0}^{2} (x), ϕ_{0}^{2} (x)) + C_{0}},

(2.2)

where $R e (u)$ represents the real part of $u$ .

Theorem 2.1.

The SAV reformulation (2.1) preserves the following quadratic energy

\frac{d}{d t} E (t) = 0, E (t) = (L ϕ, ϕ) + \frac{β}{2} p^{2} - \frac{β}{2} C_{0}, t \geq 0.

(2.3)

Proof.

It is clear to see

	$\frac{d}{d t} E (t)$	$= 2 Re (L ϕ, \partial_{t} ϕ) + β p \partial_{t} p$
		$=2βpRe(Lϕ,\rm i\|ϕ\|2ϕ√(ϕ2,ϕ2)+C0)+βp∂tp$
		$= - β p \partial_{t} p + β p \partial_{t} p$
		$= 0.$

This completes the proof. ∎

3 The construction of the high-order linearly implicit exponential integrator

In this section, a class of high-order linearly implicit energy-preserving integrators are proposed for the SAV reformulated system (2.1). For simplicity, in this paper, we shall introduce our schemes in two space dimension, i.e., $d = 2$ in (2.1). Generalizations to $d = 1$ or $3$ are straightforward.

3.1 Spatial discretization

We set computational domain $Ω = [a, b] \times [c, d]$ and let $l_{x} = b - a$ and $l_{y} = d - c$ . Choose the mesh sizes $h_{x} = l_{x} / N_{x}$ and $h_{y} = l_{y} / N_{y}$ with $N_{x}$ and $N_{y}$ two even positive integers, and denote the grid points by $x_{j} = j h_{x}$ for $j = 0, 1, 2, \dots, N_{x}$ and $y_{k} = k h_{y}$ for $k = 0, 1, 2, \dots, N_{y}$ ; let $ϕ_{j, k}$ be the numerical approximation of $ϕ (x_{j}, y_{k}, t)$ for $j = 0, 1, \dots, N_{x}, k = 0, 1, 2, \dots, N_{y}$ , and

ϕ := (ϕ_{0, 0}, ϕ_{1, 0}, \dots, ϕ_{N_{x} - 1, 0}, ϕ_{0, 1}, ϕ_{1, 1}, \dots, ϕ_{N_{x} - 1, 1}, \dots, ϕ_{0, N_{y} - 1}, ϕ_{1, N_{y} - 1}, \dots, ϕ_{N_{x} - 1, N_{y} - 1})^{T}

be the solution vector; we also define discrete inner product and norm as

⟨ ϕ, ψ ⟩_{h} = h_{x} h_{y} N_{x} - 1 \sum j = 0 N_{y} - 1 \sum k = 0 ϕ_{j, k} {¯ ψ}_{j, k}, ∥ ϕ ∥_{h}^{2} = ⟨ ϕ, ϕ ⟩_{h} .

In addition, we denote $‘ \cdot$ ’ as the element product of vectors $ϕ$ and $ψ$ , that is

ϕ \cdot ψ =

(ϕ_{0, 0} ψ_{0, 0}, \dots, ϕ_{N_{x} - 1, 0} ψ_{N_{x} - 1, 0}, \dots, ϕ_{0, N_{x} - 1} ψ_{0, N_{y} - 1}, \dots, ϕ_{N_{x} - 1, N_{y} - 1} ψ_{N_{x} - 1, N_{y} - 1})^{T} .

For brevity, we denote $ϕ \cdot ϕ$ and $ϕ \cdot ¯ ϕ$ as $ϕ^{2}$ and $| ϕ |^{2}$ , respectively.

Denote the interpolation space as

T_{h} = span {g_{j} (x) g_{k} (y), 0 \leq j \leq N_{x} - 1, 0 \leq k \leq N_{y} - 1}

where $g_{j} (x)$ and $g_{k} (y)$ are trigonometric polynomials of degree $N_{x} / 2$ and $N_{y} / 2$ , given, respectively, by

g_{j} (x) = \frac{1}{N_{x}} N_{x} / 2 \sum l = - N_{x} / 2 \frac{1}{a_{l}} e^{i l μ_{x} (x - x_{j})}, g_{k} (y) = \frac{1}{N_{y}} N_{y} / 2 \sum l = - N_{y} / 2 \frac{1}{b_{l}} e^{i l μ_{y} (y - y_{k})},

with $a_{l} = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} 1, | l | < \frac{N_{x}}{2}, 2, | l | = \frac{N_{x}}{2}, \end{matrix}, b_{l} = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} 1, | l | < \frac{N_{y}}{2}, 2, | l | = \frac{N_{y}}{2}, \end{matrix}$ , $μ_{x} = \frac{2 π}{l_{x}}$ and $μ_{y} = \frac{2 π}{l_{y}}$ . We then define the interpolation operator $I_{N} : C (Ω) \to T_{h}$ as [18]:

I_{N} ϕ (x, y, t) = N_{x} - 1 \sum j = 0 N_{y} - 1 \sum k = 0 ϕ_{j, k} (t) g_{j} (x) g_{k} (y),

where $ϕ_{j, k} (t) = ϕ (x_{j}, y_{k}, t)$ .

Taking the derivative with respect to $x$ and $y$ , respectively, and then evaluating the resulting expressions at the collocation points ( $x_{j}, y_{k}$ ), we have

	$\frac{\partial^{2} I_{N} ϕ (x_{j}, y_{k}, t)}{\partial x^{2}}$	$= N_{x} - 1 \sum l = 0 ϕ_{l, k} (t) \frac{d^{2} g_{l} (x_{j})}{d x^{2}} = N_{x} - 1 \sum l = 0 (D_{2}^{x})_{j, l} ϕ_{l, k} (t),$		(3.1)
	$\frac{\partial^{2} I_{N} ϕ (x_{j}, y_{k}, t)}{\partial y^{2}}$	$= N_{y} - 1 \sum l = 0 ϕ_{j, l} \frac{d^{2} g_{l} (y_{k})}{d y^{2}} = N_{y} - 1 \sum l = 0 ϕ_{j, l} (t) (D_{2}^{y})_{k, l},$		(3.2)

where

Remark 3.1.

We should note that [50]

D_{2}^{w} = F_{w}^{H} Λ_{w} F_{w}, w = x, y,

where

Λw=−(2πlw)2\rm diag[02,12,⋯,( Nw2)2,(−Nw2+1)2,⋯,(−2)2,(−1)2],

and $F_{w}$

is the discrete Fourier transform (DFT) and

F_{w}^{H}

represents the conjugate transpose of

F_{w}

Then, we use the standard Fourier pseudo-spectral method to solve (2.1)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{d}{d t} ϕ = i (L_{h} ϕ + \frac{β | ϕ |^{2} \cdot ϕ p}{\sqrt{⟨ ϕ^{2}, ϕ^{2} ⟩_{h} + C_{0}}}), \frac{d}{d t} p = 2 R e ⟨ i L_{h} ϕ, \frac{| ϕ |^{2} \cdot ϕ}{\sqrt{⟨ ϕ^{2}, ϕ^{2} ⟩_{h} + C_{0}}} ⟩_{h}, \end{matrix}

(3.3)

where $L_{h} = I_{N_{y}} \otimes D_{2}^{x} + D_{2}^{y} \otimes I_{N_{x}}$ is the spectral differentiation matrix and $\otimes$ represents the Kronecker product.

Theorem 3.1.

The semi-discrete system (3.3) preserves the following semi-discrete quadratic energy

\frac{d}{d t} E_{h} (t) = 0, E_{h} (t) = ⟨ L_{h} ϕ, ϕ ⟩_{h} + \frac{β}{2} p^{2} - \frac{β}{2} C_{0}, t \geq 0.

(3.4)

Proof.

The proof is similar to Theorem 2.1, thus, for brevity, we omit it.

Remark 3.2.

When the standard Fourier pseudo-spectral method is employed to the system (1.1) for spatial discretizations, we can obtain semi-discrete Hamiltonian energy as

E_{h} (t) = ⟨ ϕ, L_{h} ϕ ⟩_{h} + \frac{β}{2} N_{x} - 1 \sum j = 0 N_{y} - 1 \sum k = 0 | ϕ_{j, k} |^{4}, t \geq 0.

(3.5)

We note that however, the quadratic energy (3.4) is only equivalent to the Hamiltonian energy (3.5) in the continuous sense, but not for the discrete sense.

3.2 Temporal exponential integration

Denote $t_{n} = n τ,$ and $t_{n i} = t_{n} + c_{i} τ, i = 1, 2, \dots, s$ , $n = 0, 1, 2, \dots,$ where $τ$ is the time step and $c_{1}, c_{2}, \dots, c_{s}$ are distinct real numerbers (usually $0 \leq c_{i} \leq 1$ ). The approximations of the function $ϕ (x, y, t)$ at points $(x_{j}, y_{k}, t_{n})$ and $(x_{j}, y_{k}, t_{n i})$ are denoted by $ϕ_{j, k}^{n}$ and $(ϕ_{n i})_{j, k}$ and the approximations of the function $p (t)$ at points $t_{n}$ and $t_{n i}$ are denoted by $p^{n}$ and $p_{n i}$ .

We first apply the extrapolation technique to the nonlinear term of (3.3) and a linearized system is obtained, as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{d}{d t} ~ ϕ = i (L_{h} ~ ϕ + \frac{β | ϕ_{n i}^{*} |^{2} \cdot ϕ_{n i}^{*} ~ p}{\sqrt{⟨ (ϕ_{n i}^{*})^{2}, (ϕ_{n i}^{*})^{2} ⟩_{h} + C_{0}}}), \frac{d}{d t} ~ p = 2 R e ⟨ i L_{h} ~ ϕ, \frac{| ϕ_{n i}^{*} |^{2} \cdot ϕ_{n i}^{*}}{\sqrt{⟨ (ϕ_{n i}^{*})^{2}, (ϕ_{n i}^{*})^{2} ⟩_{h} + C_{0}}} ⟩_{h}, \end{matrix}

(3.6)

where $~ ϕ := ~ ϕ (t)$ and $~ p := ~ p (t)$ are approximations of $ϕ (t)$ and $p (t)$ , respectively over time interval $[t_{n}, t_{n + 1}], n = 0, 1, 2, \dots,$ . Here, $ϕ_{n i}^{*}$ is an (explicit) extrapolation approximation to $ϕ (t_{n} + c_{i} τ), i = 1, 2, \dots, s$ of the order $O (τ^{s + 1})$ . For the more details on the construction of $ϕ_{n i}^{*}$ , please refer to Refs. [27, 40].

By using the Lawson transformation [39], we multiply both sides of the first equation of (3.6) by the operator $e x p (i L_{h} t)$ , and then introduce $Ψ = e x p (- i L_{h} t) ϕ$ to transform (3.6) into an equivalent form

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{d}{d t} Ψ = i exp (- i L_{h} t) \frac{β | ϕ_{n i}^{*} |^{2} \cdot ϕ_{n i}^{*} p}{\sqrt{⟨ (ϕ_{n i}^{*})^{2}, (ϕ_{n i}^{*})^{2} ⟩_{h} + C_{0}}}, \frac{d}{d t} p = 2 R e ⟨ i exp (i L_{h} t) L_{h} Ψ, \frac{| ϕ_{n i}^{*} |^{2} \cdot ϕ_{n i}^{*}}{\sqrt{⟨ (ϕ_{n i}^{*})^{2}, (ϕ_{n i}^{*})^{2} ⟩_{h} + C_{0}}} ⟩_{h}, \end{matrix}

(3.7)

with the consistent initial condition

Ψ (x, 0) = ϕ_{0} (x), p (0) = \sqrt{⟨ ϕ_{0}^{2} (x), ϕ_{0}^{2} (x) ⟩_{h} + C_{0}} .

(3.8)

We first apply an RK method to the linearized system (3.7), and then the discretization is rewritten in terms of the original variables to give a class of linearly implicit exponential integrations (LI-EIs) for solving (1.1)

(3.9)

where $i = 1, 2, \dots, s$ . Then $(ϕ^{n + 1}, p^{n + 1})$ are updated by

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} ϕ^{n + 1} = exp (i L_{h} τ) ϕ^{n} + τ s \sum i = 1 b_{i} exp (i L_{h} (1 - c_{i}) τ) k_{i}, p^{n + 1} = p^{n} + τ s \sum i = 1 b_{i} l_{i} . \end{matrix}

(3.10)

Remark 3.3.

According to the definition of matrix-valued function [31], it holds

$exp (i L_{h} τ)^{H} = exp (- i L_{h} τ)$ and $exp (i L_{h} τ) L_{h} = L_{h} exp (i L_{h} τ)$ ;
$exp (i L_{h} τ) = F_{N_{y}}^{H} \otimes F_{N_{x}}^{H} exp (i (I_{N_{y}} \otimes Λ_{x} + Λ_{y} \otimes I_{N_{x}}) τ) F_{N_{y}} \otimes F_{N_{x}}$ , which can be efficiently implemented by the matlab functions fftn.m and ifftn.m.

Remark 3.4.

We consider a Hamiltonian PDEs system given by

\partial_{t} z = D (A z + F^{^{'}} (z)), z \in Ω \subset R^{d}, t > 0, d = 1, 2, 3,

(3.11)

with a Hamiltonian energy

E (t) = \frac{1}{2} (z, A z) + (F (z), 1), t \geq 0,

(3.12)

where $z := z (x, t),$ $D$

is a skew-adjoint operation and

A

is a self-adjoint operation and

\int_{Ω} F (z (x, t)) d x

is bounded from below. By introducing a scalar auxiliary variable

p := p (t) = \sqrt{2 (F (z), 1) + 2 C_{0}},

where $C_{0}$ is a constant large enough to make $p$ well-defined for all $z (x, t), x \in Ω$ . On the basis of the energy-variational principle, the system (3.11) can be reformulated into

(3.13)

which preserves the following quadratic energy

\frac{d}{d t} E (t) = 0, E (t) = \frac{1}{2} (z, A z) + \frac{1}{2} p^{2} - C_{0}, t \geq 0.

(3.14)

The proposed linearly implicit schemes can be easily generalized to solve the above system.

Theorem 3.2.

If the coefficients of (3.7) satisfy

b_{i} a_{i j} + b_{j} a_{j i} = b_{i} b_{j}, \forall i, j = 1, \dots, s,

(3.15)

the proposed scheme (3.9)-(3.10) can preserve the discrete quadratic energy, as follows:

E_{h}^{n} = E_{h}^{0}, E_{h}^{n} = ⟨ L_{h} ϕ^{n}, ϕ^{n} ⟩_{h} + \frac{β}{2} (p^{n})^{2} - \frac{β}{2} C_{0}, n = 0, 1, 2, \dots .

(3.16)

Proof.

We obtain from the first equality of (3.10) and (3.15) that

		$⟨ L_{h} ϕ^{n + 1}, ϕ^{n + 1} ⟩_{h} - ⟨ L_{h} ϕ^{n}, ϕ^{n} ⟩_{h}$
	$=$	$2 τ s \sum i = 1 b_{i} R e ⟨ exp (i c_{i} τ) ϕ^{n}, L_{h} k_{i} ⟩_{h} + τ^{2} s \sum i, j = 1 b_{i} b_{j} ⟨ L_{h} k_{i}, exp (i L_{h} (c_{i} - c_{j}) τ) k_{j} ⟩_{h}$
	$=$	$2 τ s \sum i = 1 b_{i} R e ⟨ exp (i L_{h} c_{i} τ) ϕ^{n}, L_{h} k_{i} ⟩_{h} + τ^{2} s \sum i, j = 1 (b_{i} a_{i j} + b_{j} a_{j i}) ⟨ L_{h} k_{i}, exp (i L_{h} (c_{i} - c_{j}) τ) k_{j} ⟩_{h}$
	$=$	$2 τ s \sum i = 1 b_{i} R e ⟨ exp (i L_{h} c_{i} τ) ϕ^{n}, L_{h} k_{i} ⟩_{h} + 2 τ^{2} s \sum i, j = 1 b_{i} a_{i j} R e ⟨ exp (i L_{h} (c_{i} - c_{j}) τ) k_{j}, L_{h} k_{i} ⟩_{h} .$

We then can deduce from the first equality of (3.9) that

	$2 τ s \sum i = 1 b_{i} R e$	$⟨ exp (i L_{h} c_{i} τ) ϕ^{n}, L_{h} k_{i} ⟩_{h}$

Thus, we can obtain

⟨ L_{h} ϕ^{n + 1}, ϕ^{n + 1} ⟩_{h} - ⟨ L_{h} ϕ^{n}, ϕ^{n} ⟩_{h} = 2 τ s \sum i = 1 b_{i} R e ⟨ ϕ_{n i}, L_{h} k_{i} ⟩_{h} .

(3.17)

Moreover, we derive from (3.9) and (3.10) that

$\frac{β}{2} [(p^{n + 1})^{2} - (p^{n})^{2}]$	$= β τ s \sum i = 1 b_{i} l_{i} p^{n} + \frac{β}{2} τ^{2} s \sum i, j = 1 b_{i} b_{j} l_{i} l_{j}$
	$= β τ s \sum i = 1 b_{i} l_{i} (p_{n i} - τ s \sum j = 1 a_{i j} l_{j}) + \frac{β}{2} τ^{2} s \sum i, j = 1 (b_{i} a_{i j} + b_{j} a_{j i}) l_{i} l_{j}$

	$= - 2 τ s \sum i = 1 b_{i} R e ⟨ ϕ_{n i}, L_{h} k_{i} ⟩_{h} .$	(3.18)

It follows from (3.17) and (3.2), we completes the proof.∎

Remark 3.5.

It follows from Remark 3.2 that the proposed scheme cannot preserve the discrete Hamiltonian energy (3.5). In addition, the nonlinear terms of the semi-discrete system (3.3) are explicitly discretized, thus our scheme cannot preserve the discrete mass, as follows:

M_{h}^{n} = h_{x} h_{y} N_{x} - 1 \sum j = 0 N_{y} - 1 \sum k = 0 | ϕ_{j, k}^{n} |^{2} .

(3.19)

Remark 3.6.

The Gauss collocation method is symplectic (see Refs. [30, 48] and references therein), thus, it can preserve the discrete quadratic energy (see (3.16)). In particular, the coefficients of Gauss collocation methods of order 4 and 6 can be given explicitly by (see [30])

$\frac{1}{2} - \frac{\sqrt{3}}{6}$	$\frac{1}{4}$	$\frac{1}{4} - \frac{\sqrt{3}}{6}$
$\frac{1}{2} + \frac{\sqrt{3}}{6}$	$\frac{1}{4} + \frac{\sqrt{3}}{6}$	$\frac{1}{4}$
	$\frac{1}{2}$	$\frac{1}{2}$

, $\frac{1}{2} - \frac{\sqrt{15}}{10}$ $\frac{5}{36}$ $\frac{2}{9} - \frac{\sqrt{15}}{15}$ $\frac{5}{36} - \frac{\sqrt{15}}{30}$ $\frac{1}{2}$ $\frac{5}{36} + \frac{\sqrt{15}}{24}$ $\frac{2}{9}$ $\frac{5}{36} - \frac{\sqrt{15}}{24}$ $\frac{1}{2} + \frac{\sqrt{15}}{10}$ $\frac{5}{36} + \frac{\sqrt{15}}{30}$ $\frac{2}{9} + \frac{\sqrt{15}}{15}$ $\frac{5}{36}$ $\frac{5}{18}$ $\frac{4}{9}$ $\frac{5}{18}$ .

Table. 1: RK coefficients of Gaussian collocation methods of order 4 (left) and 6 (right).

Besides its energy-preserving property, a most remarkable thing about the above scheme (3.9)-(3.10) is that it can be solved efficiently. For simplicity, we take $s = 2$ as an example.

We denote

γ_{i}^{*} = i \frac{| ϕ_{n i}^{*} |^{2} \cdot ϕ_{n i}^{*}}{\sqrt{⟨ (ϕ_{n i}^{*})^{2}, (ϕ_{n i}^{*})^{2} ⟩_{h} + C_{0}}}, i = 1, 2,

(3.20)

and rewrite

l_{1} = - 2 Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{*} ⟩_{h}, l_{2} = - 2 Re ⟨ L_{h} ϕ_{n 2}, γ_{2}^{*} ⟩_{h} .

(3.21)

With (3.20), we have

	$k_{1} = β γ_{1}^{} p^{n} - 2 β τ a_{11} γ_{1}^{} Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{} ⟩_{h} - 2 β τ a_{12} γ_{1}^{} Re ⟨ L_{h} ϕ_{n 2}, γ_{2}^{*} ⟩_{h},$		(3.22)
	$k_{2} = β γ_{2}^{} p^{n} - 2 β τ a_{21} γ_{2}^{} Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{} ⟩_{h} - 2 β τ a_{22} γ_{2}^{} Re ⟨ L_{h} ϕ_{n 2}, γ_{2}^{*} ⟩_{h} .$		(3.23)

Then it follows from the first equality of (3.9) and (3.22)-(3.23) that

$ϕ_{n 1}$	$= γ_{1}^{n} - (2 β τ^{2} a_{11}^{2} γ_{1}^{} + 2 β τ^{2} a_{12} a_{21} exp (i L_{h} (c_{1} - c_{2}) τ) γ_{2}^{}) Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{*} ⟩_{h}$
		(3.24)
$ϕ_{n 2}$	$= γ_{2}^{n} - (2 β τ^{2} a_{21} a_{11} exp (i L_{h} (c_{2} - c_{1}) τ) γ_{1}^{} + 2 β τ^{2} a_{22} a_{21} γ_{2}^{}) Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{*} ⟩_{h}$
		(3.25)

where

	$γ_{1}^{n} = exp (i L_{h} c_{1} τ) ϕ^{n} + τ β a_{11} γ_{1}^{} p^{n} + τ β a_{12} exp (i L_{h} (c_{1} - c_{2}) τ) γ_{2}^{} p^{n},$
	$γ_{2}^{n} = exp (i L_{h} c_{2} τ) ϕ^{n} + τ β a_{21} exp (i L_{h} (c_{2} - c_{1}) τ) γ_{1}^{} p^{n} + τ β a_{22} γ_{2}^{} p^{n} .$

Multiply both sides of (3.2) and (3.2) with $L_{h}$ and we then take discrete inner product with $γ_{1}^{*}$ and $γ_{2}^{*}$ , respectively, to obtain

	$A_{11} Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{} ⟩_{h} + A_{12} Re ⟨ L_{h} ϕ_{n 2}, γ_{2}^{} ⟩_{h} = Re ⟨ L_{h} γ_{1}^{n}, γ_{1}^{*} ⟩_{h},$		(3.26)
	$A_{21} Re ⟨ L_{h} ϕ_{n 1}, γ_{1}^{} ⟩_{h} + A_{22} Re ⟨ L_{h} ϕ_{n 2}, γ_{2}^{} ⟩_{h} = Re ⟨ L_{h} γ_{2}^{n}, γ_{2}^{*} ⟩_{h},$		(3.27)

where

	$A_{11} = 1 + Re ⟨ 2 β τ^{2} a_{11}^{2} L_{h} γ_{1}^{} + 2 β τ^{2} a_{12} a_{21} exp (i L_{h} (c_{1} - c_{2}) τ) L_{h} γ_{2}^{}, γ_{1}^{*} ⟩_{h},$		(3.28)
	$A_{12} = Re ⟨ 2 β τ^{2} a_{11} a_{12} L_{h} γ_{1}^{} + 2 β τ^{2} a_{12} a_{22} exp (i L_{h} (c_{1} - c_{2}) τ) L_{h} γ_{2}^{}, γ_{1}^{*} ⟩_{h},$		(3.29)
	$A_{21} = Re ⟨ 2 β τ^{2} a_{21} a_{11} exp (i L_{h} (c_{2} - c_{1}) τ) L_{h} γ_{1}^{} + 2 β τ^{2} a_{22} a_{21} L_{h} γ_{2}^{}, γ_{2}^{*} ⟩_{h},$		(3.30)
	$A_{22} = 1 + Re ⟨ 2 β τ^{2} a_{21} a_{12} exp (i L_{h} (c_{2} - c_{1}) τ) L_{h} γ_{1}^{} + 2 β τ^{2} a_{22}^{2} L_{h} γ_{2}^{}, γ_{2}^{*} ⟩_{h} .$		(3.31)

Eqs. (3.26) and (3.27) form a $2 \times 2$ linear system for the unknowns

Solving from the $2 \times 2$ linear system (3.26) and (3.27), and $l_{i}, k_{i}$ and $ϕ_{n i}, i = 1, 2$ are then updated from (3.21)-(3.2), respectively. Subsequently, $ϕ^{n + 1}$ and $p^{n + 1}$ are obtained by (3.10).

To summarize, the scheme (3.9)-(3.10) of order 3 and 4 can be efficiently implemented, respectively, as follows:

Algorithm 1: Linearly implicit exponential integration (3.9)-(3.10) of order 3 (denote by 3rd-LI-EI)

Initialization: apply the fourth-order exponential integration [20] to compute

ϕ_{01}, ϕ_{02}, ϕ^{1}

and

p^{1}

For n=1,2,

\dots

, do

Compute the extrapolation approximation of order 3 [27]

ϕ_{n 1}^{*} = (

High-order linearly implicit structure-preserving exponential integrators for the nonlinear Schrödinger equation

Authors

A linearly implicit structure-preserving scheme for the Camassa-Holm equation based on multiple scalar auxiliary variables approach

Arbitrary high-order linearly implicit energy-preserving algorithms for Hamiltonian PDEs

A high order fully discrete scheme for the Korteweg-de vries equation with a time-stepping procedure of Runge-Kutta-composition type

An asymptotic-preserving IMEX method for nonlinear radiative transfer equation

On energy preserving high-order discretizations for nonlinear acoustics

Nonlinear Fokker-Planck Acceleration for Forward-Peaked Transport Problems in Slab Geometry

Arbitrary high-order structure-preserving methods for the quantum Zakharov system

1 Introduction

2 Model reformulation

Theorem 2.1.

Proof.

3 The construction of the high-order linearly implicit exponential integrator

3.1 Spatial discretization

Remark 3.1.

Theorem 3.1.

Proof.

Remark 3.2.

3.2 Temporal exponential integration

Remark 3.3.

Remark 3.4.

Theorem 3.2.

Proof.

Remark 3.5.

Remark 3.6.

High-order linearly implicit structure-preserving exponential integrators for the nonlinear Schrödinger equation

Authors

Related Research

A linearly implicit structure-preserving scheme for the Camassa-Holm equation based on multiple scalar auxiliary variables approach

Arbitrary high-order linearly implicit energy-preserving algorithms for Hamiltonian PDEs

A high order fully discrete scheme for the Korteweg-de vries equation with a time-stepping procedure of Runge-Kutta-composition type

An asymptotic-preserving IMEX method for nonlinear radiative transfer equation

On energy preserving high-order discretizations for nonlinear acoustics

Nonlinear Fokker-Planck Acceleration for Forward-Peaked Transport Problems in Slab Geometry

Arbitrary high-order structure-preserving methods for the quantum Zakharov system

This week in AI

1 Introduction

2 Model reformulation

Theorem 2.1.

Proof.

3 The construction of the high-order linearly implicit exponential integrator

3.1 Spatial discretization

Remark 3.1.

Theorem 3.1.

Proof.

Remark 3.2.

3.2 Temporal exponential integration

Remark 3.3.

Remark 3.4.

Theorem 3.2.

Proof.

Remark 3.5.

Remark 3.6.