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

A novel class of high-order linearly implicit energy-preserving exponential integrators are proposed for the nonlinear Schrödinger equation. We firstly done that the original equation is reformulated into a new form with a modified quadratic energy by the scalar auxiliary variable approach. The spatial derivatives of the system are then approximated with the standard Fourier pseudo-spectral method. Subsequently, we apply the extrapolation technique to the nonlinear term of the semi-discretized system and a linearized system is obtained. Based on the Lawson transformation, the linearized system is rewritten as an equivalent one and we further apply the symplectic Runge-Kutta method to the resulting system to gain a fully discrete scheme. We show that the proposed scheme can produce numerical solutions along which the modified energy is precisely conserved, as is the case with the analytical solution and is extremely efficient in the sense that only linear equations with constant coefficients need to be solved at every time step. Numerical results are addressed to demonstrate the remarkable superiority of the proposed schemes in comparison with other high-order structure-preserving method.

## Authors

• 10 publications
• 4 publications
• 6 publications
• 4 publications
06/29/2019

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

In this paper, we present a linearly implicit energy-preserving scheme f...
11/17/2020

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

In this paper, we present a novel strategy to systematically construct l...
05/26/2020

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

We consider the periodic initial-value problem for the Korteweg-de Vries...
08/15/2020

### An asymptotic-preserving IMEX method for nonlinear radiative transfer equation

We present an asymptotic preserving method for the radiative transfer eq...
12/15/2019

### On energy preserving high-order discretizations for nonlinear acoustics

This paper addresses the numerical solution of the Westervelt equation, ...
05/13/2020

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

This paper introduces a nonlinear acceleration technique that accelerate...
02/26/2022

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

In this paper, we present a new methodology to develop arbitrary high-or...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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

 {i∂tϕ(x,t)+Lϕ(x,t)+β|ϕ(x,t)|2ϕ(x,t)=0, (x,t)∈Ω, t>0,ϕ(x,0)=ϕ0(x), x∈Ω⊂Rd, d=1,2,3, (1.1)

where is the time variable, is the spatial variable, is the complex unit, is the complex-valued wave function, is the usual Laplace operator, is a given real constant, and is a given initial data. The nonlinear Schrödinger equation (1.1) satisfies the following two invariants, that is mass

 M(t)=∫Ω|ϕ|2dx≡M(0), t≥0, (1.2)

and Hamiltonian energy

 E(t)=∫Ω(−|∇ϕ|2+β2|ϕ|4)dx≡E(0), t≥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)=√(ϕ2,ϕ2)+C0,

where is a constant large enough to make well-defined for all . Here, is the inner product defined by where represents the conjugate of . Then, the Hamiltonian energy functional is rewritten as the following quadratic form

 E(t)=(Lϕ,ϕ)+β2p2−β2C0.

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂tϕ=i(Lϕ+β|ϕ|2ϕp√(ϕ2,ϕ2)+C0),ddtp=2Re(iLϕ,|ϕ|2ϕ√(ϕ2,ϕ2)+C0), (2.1)

with the consistent initial condition

 ϕ(x,0)=ϕ0(x), p(0)=√(ϕ20(x),ϕ20(x))+C0, (2.2)

where represents the real part of .

###### Theorem 2.1.

The SAV reformulation (2.1) preserves the following quadratic energy

 ddtE(t)=0, E(t)=(Lϕ,ϕ)+β2p2−β2C0, t≥0. (2.3)
###### Proof.

It is clear to see

 ddtE(t) =2Re(Lϕ,∂tϕ)+βp∂tp =2βpRe(Lϕ,\rm i|ϕ|2ϕ√(ϕ2,ϕ2)+C0)+βp∂tp =−βp∂tp+βp∂tp =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., in (2.1). Generalizations to or are straightforward.

### 3.1 Spatial discretization

We set computational domain and let and . Choose the mesh sizes and with and two even positive integers, and denote the grid points by for and for ; let be the numerical approximation of for , and

 ϕ:=(ϕ0,0,ϕ1,0,⋯,ϕNx−1,0,ϕ0,1,ϕ1,1,⋯,ϕNx−1,1,⋯,ϕ0,Ny−1,ϕ1,Ny−1,⋯,ϕNx−1,Ny−1)T

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

 ⟨ϕ,ψ⟩h=hxhyNx−1∑j=0Ny−1∑k=0ϕj,k¯ψj,k, ∥ϕ∥2h=⟨ϕ,ϕ⟩h.

In addition, we denote ’ as the element product of vectors and , that is

 ϕ⋅ψ= (ϕ0,0ψ0,0,⋯,ϕNx−1,0ψNx−1,0,⋯,ϕ0,Nx−1ψ0,Ny−1,⋯,ϕNx−1,Ny−1ψNx−1,Ny−1)T.

For brevity, we denote and as and , respectively.

Denote the interpolation space as

 Th=span{gj(x)gk(y), 0≤j≤Nx−1,0≤k≤Ny−1}

where and are trigonometric polynomials of degree and , given, respectively, by

 gj(x)=1NxNx/2∑l=−Nx/21aleilμx(x−xj), gk(y)=1NyNy/2∑l=−Ny/21bleilμy(y−yk),

with , and . We then define the interpolation operator as [18]:

 INϕ(x,y,t)=Nx−1∑j=0Ny−1∑k=0ϕj,k(t)gj(x)gk(y),

where .

Taking the derivative with respect to and , respectively, and then evaluating the resulting expressions at the collocation points (), we have

 ∂2INϕ(xj,yk,t)∂x2 =Nx−1∑l=0ϕl,k(t)d2gl(xj)dx2=Nx−1∑l=0(Dx2)j,lϕl,k(t), (3.1) ∂2INϕ(xj,yk,t)∂y2 =Ny−1∑l=0ϕj,ld2gl(yk)dy2=Ny−1∑l=0ϕj,l(t)(Dy2)k,l, (3.2)

where

###### Remark 3.1.

We should note that [50]

 Dw2=FHwΛwFw, w=x,y,

where

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

and

is the discrete Fourier transform (DFT) and

represents the conjugate transpose of .

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ddtϕ=i(Lhϕ+β|ϕ|2⋅ϕp√⟨ϕ2,ϕ2⟩h+C0),ddtp=2Re⟨iLhϕ,|ϕ|2⋅ϕ√⟨ϕ2,ϕ2⟩h+C0⟩h, (3.3)

where is the spectral differentiation matrix and represents the Kronecker product.

###### Theorem 3.1.

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

 ddtEh(t)=0, Eh(t)=⟨Lhϕ,ϕ⟩h+β2p2−β2C0, t≥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

 Eh(t)=⟨ϕ,Lhϕ⟩h+β2Nx−1∑j=0Ny−1∑k=0|ϕj,k|4, t≥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 and , where is the time step and are distinct real numerbers (usually ). The approximations of the function at points and are denoted by and and the approximations of the function at points and are denoted by and .

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ddt~ϕ=i(Lh~ϕ+β|ϕ∗ni|2⋅ϕ∗ni~p√⟨(ϕ∗ni)2,(ϕ∗ni)2⟩h+C0),ddt~p=2Re⟨iLh~ϕ,|ϕ∗ni|2⋅ϕ∗ni√⟨(ϕ∗ni)2,(ϕ∗ni)2⟩h+C0⟩h, (3.6)

where and are approximations of and , respectively over time interval . Here, is an (explicit) extrapolation approximation to of the order . For the more details on the construction of , 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 , and then introduce to transform (3.6) into an equivalent form

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ddtΨ=iexp(−iLht)β|ϕ∗ni|2⋅ϕ∗nip√⟨(ϕ∗ni)2,(ϕ∗ni)2⟩h+C0,ddtp=2Re⟨iexp(iLht)LhΨ,|ϕ∗ni|2⋅ϕ∗ni√⟨(ϕ∗ni)2,(ϕ∗ni)2⟩h+C0⟩h, (3.7)

with the consistent initial condition

 Ψ(x,0)=ϕ0(x), p(0)=√⟨ϕ20(x),ϕ20(x)⟩h+C0. (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 . Then are updated by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ϕn+1=exp(iLhτ)ϕn+τs∑i=1biexp(iLh(1−ci)τ)ki,pn+1=pn+τs∑i=1bili. (3.10)
###### Remark 3.3.

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

• and ;

• , 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

 ∂tz=D(Az+F′(z)), z∈Ω⊂Rd, t>0,d=1,2,3, (3.11)

with a Hamiltonian energy

 E(t)=12(z,Az)+(F(z),1), t≥0, (3.12)

where

is a self-adjoint operation and is bounded from below. By introducing a scalar auxiliary variable

 p:=p(t)=√2(F(z),1)+2C0,

where is a constant large enough to make well-defined for all . On the basis of the energy-variational principle, the system (3.11) can be reformulated into

 (3.13)

which preserves the following quadratic energy

 ddtE(t)=0, E(t)=12(z,Az)+12p2−C0, t≥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

 biaij+bjaji=bibj, ∀ i,j=1,⋯,s, (3.15)

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

 Enh=E0h, Enh=⟨Lhϕn,ϕn⟩h+β2(pn)2−β2C0, n=0,1,2,⋯. (3.16)
###### Proof.

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

 ⟨Lhϕn+1,ϕn+1⟩h−⟨Lhϕn,ϕn⟩h = 2τs∑i=1biRe⟨exp(iciτ)ϕn,Lhki⟩h+τ2s∑i,j=1bibj⟨Lhki,exp(iLh(ci−cj)τ)kj⟩h = 2τs∑i=1biRe⟨exp(iLhciτ)ϕn,Lhki⟩h+τ2s∑i,j=1(biaij+bjaji)⟨Lhki,exp(iLh(ci−cj)τ)kj⟩h = 2τs∑i=1biRe⟨exp(iLhciτ)ϕn,Lhki⟩h+2τ2s∑i,j=1biaijRe⟨exp(iLh(ci−cj)τ)kj,Lhki⟩h.

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

 2τs∑i=1biRe ⟨exp(iLhciτ)ϕn,Lhki⟩h

Thus, we can obtain

 ⟨Lhϕn+1,ϕn+1⟩h−⟨Lhϕn,ϕn⟩h=2τs∑i=1biRe⟨ϕni,Lhki⟩h. (3.17)

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

 β2[(pn+1)2−(pn)2] =βτs∑i=1bilipn+β2τ2s∑i,j=1bibjlilj =βτs∑i=1bili(pni−τs∑j=1aijlj)+β2τ2s∑i,j=1(biaij+bjaji)lilj =−2τs∑i=1biRe⟨ϕni,Lhki⟩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:

 Mnh=hxhyNx−1∑j=0Ny−1∑k=0|ϕnj,k|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])

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 as an example.

We denote

 γ∗i=i|ϕ∗ni|2⋅ϕ∗ni√⟨(ϕ∗ni)2,(ϕ∗ni)2⟩h+C0, i=1,2, (3.20)

and rewrite

 l1=−2Re⟨Lhϕn1,γ∗1⟩h, l2=−2Re⟨Lhϕn2,γ∗2⟩h. (3.21)

With (3.20), we have

 k1=βγ∗1pn−2βτa11γ∗1Re⟨Lhϕn1,γ∗1⟩h−2βτa12γ∗1Re⟨Lhϕn2,γ∗2⟩h, (3.22) k2=βγ∗2pn−2βτa21γ∗2Re⟨Lhϕn1,γ∗1⟩h−2βτa22γ∗2Re⟨Lhϕn2,γ∗2⟩h. (3.23)

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

 ϕn1 =γn1−(2βτ2a211γ∗1+2βτ2a12a21exp(iLh(c1−c2)τ)γ∗2)Re⟨Lhϕn1,γ∗1⟩h (3.24) ϕn2 =γn2−(2βτ2a21a11exp(iLh(c2−c1)τ)γ∗1+2βτ2a22a21γ∗2)Re⟨Lhϕn1,γ∗1⟩h (3.25)

where

 γn1=exp(iLhc1τ)ϕn+τβa11γ∗1pn+τβa12exp(iLh(c1−c2)τ)γ∗2pn, γn2=exp(iLhc2τ)ϕn+τβa21exp(iLh(c2−c1)τ)γ∗1pn+τβa22γ∗2pn.

Multiply both sides of (3.2) and (3.2) with and we then take discrete inner product with and , respectively, to obtain

 A11Re⟨Lhϕn1,γ∗1⟩h+A12Re⟨Lhϕn2,γ∗2⟩h=Re⟨Lhγn1,γ∗1⟩h, (3.26) A21Re⟨Lhϕn1,γ∗1⟩h+A22Re⟨Lhϕn2,γ∗2⟩h=Re⟨Lhγn2,γ∗2⟩h, (3.27)

where

 A11=1+Re⟨2βτ2a211Lhγ∗1+2βτ2a12a21exp(iLh(c1−c2)τ)Lhγ∗2,γ∗1⟩h, (3.28) A12=Re⟨2βτ2a11a12Lhγ∗1+2βτ2a12a22exp(iLh(c1−c2)τ)Lhγ∗2,γ∗1⟩h, (3.29) A21=Re⟨2βτ2a21a11exp(iLh(c2−c1)τ)Lhγ∗1+2βτ2a22a21Lhγ∗2,γ∗2⟩h, (3.30) A22=1+Re⟨2βτ2a21a12exp(iLh(c2−c1)τ)Lhγ∗1+2βτ2a222Lhγ∗2,γ∗2⟩h. (3.31)

Eqs. (3.26) and (3.27) form a linear system for the unknowns

Solving from the linear system (3.26) and (3.27), and and are then updated from (3.21)-(3.2), respectively. Subsequently, and 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: