Arbitrary-order energy-preserving exponential integrators for the cubic Schrödinger equation

# Arbitrary-order energy-preserving exponential integrators for the cubic Schrödinger equation

Bin Wang  111School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R. China; Mathematisches Institut, University of Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany. The research is supported in part by the Alexander von Humboldt Foundation and by the Natural Science Foundation of Shandong Province (Outstanding Youth Foundation) under Grant ZR2017JL003. E-mail: wang@na.uni-tuebingen.de    Xinyuan Wu School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R. China; Department of Mathematics, Nanjing University, Nanjing 210093, P.R. China. The research is supported in part by the National Natural Science Foundation of China under Grant 11671200. E-mail: xywu@nju.edu.cn
###### Abstract

In this paper we derive and analyse new and efficient energy-preserving exponential integrators of arbitrarily high order to solve the cubic Schrödinger Cauchy problem on a -dimensional torus. Energy preservation is a key feature of the cubic Schrödinger equation. It is proved that the novel integrators can be of arbitrarily high order which exactly preserve the continuous energy of the original continuous system. The existence and uniqueness, regularity, global convergence, nonlinear stability of the new integrators are studied in detail. One of the new energy-preserving exponential integrators is constructed and two numerical experiments are included. The numerical results illustrate the efficiency of the new integrator in comparison with existing numerical methods in the literature.

Keywords: Cubic Schrödinge equationEnergy preservationExponential integrators

MSC: 65P1035Q5565M1265M70

## 1 Introduction

It is well known that one of the cornerstones of quantum physics is the Schrödinger equation. The nonlinear Schrödinger equation has long been used to approximately describe the dynamics of complicated systems, such as Maxwell¡¯s equations for the description of nonlinear optics or the equations describing surface water waves, including rogue waves which appear from nowhere and disappears without a trace (see, e.g. [1, 2, 9]). This paper is devoted to designing and analysing novel numerical integrators to preserve the continuous energy of the cubic Schrödinger equation

 {iut+△u=λ|u|2u,  (t,x)∈[0,T]×Td,u(0,x)=u0(x),      x∈Td, (1)

where is denoted the one-dimensional torus and for the given positive integer , denotes the -dimensional torus. In this paper, we consider this equation as a Cauchy problem in time (no space discretisation is made). The solutions of this equation have conservation of the following energy

 H[u]=12(2π)d∫Td(|∇u|2+λ2|u|4)dx. (2)

It is of great interest to devise numerical schemes which can conserve the continuous version of this important invariant, and the aim of this paper is to formulate a novel kind of exponential integrators with this fundamental property.

Schrödinger equations frequently arise in a wide variety of applications including several areas of physics, fiber optics, quantum transport and other applied sciences (see, e.g. [3, 24, 35, 43, 53]). Many efficient and effective methods have been proposed for the numerical integration of Schrödinger equations, such as splitting methods (see, e.g. [5, 10, 23, 26, 44, 46, 54]), exponential-type integrators (see, e.g. [14, 15, 19, 21, 49]), multi-symplectic methods (see, e.g. [6, 50]) and other efficient methods (see, e.g. [8, 27, 34, 36]).

In recent decades, structure-preserving algorithms of differential equations have been received much attention and for the related work, we refer the reader to [25, 29, 39, 55, 57, 59, 56, 58, 60, 62] and references therein. It is well known that structure-preserving algorithms are able to exactly preserve some structural properties of the underlying continuous system. Amongst the typical subjects of structure-preserving algorithms are energy-preserving schemes, which exactly preserve the energy of the underlying system. There have been a lot of studies on this topic for Hamiltonian partial differential equations (PDEs). In [52], finite element methods were introduced systematically for numerical solution of PDEs. The authors in [20] researched discrete gradient methods for PDEs. The work in [16] investigated the average vector field (AVF) method for discretising Hamiltonian PDEs. Hamiltonian Boundary Value Methods (HBVMs) were studied for the semilinear wave equation in [12] and were recently researched for nonlinear Schrödinger equations with wave operator in [13]. The adapted AVF method for Hamiltonian wave equations was analysed in [41, 42]. Other related work is referred to [11, 18, 22, 37, 40, 45, 47, 51]. On the other hand, exponential integrators have been widely introduced and developed for solving first-order ODEs, and we refer the reader to [29, 30, 31, 32, 48] for example. This kind of methods has also been studied in the numerical integration of Schrödinger equations (see, e.g. [14, 15, 19, 21, 49]). However, until now, energy preserving exponential integrators for Schrödinger equations that exactly preserve the continuous energy, do not seem to have been studied in the literature, which motivates this paper.

With this premise, this paper is mainly concerned with arbitrary-order energy-preserving exponential integrators for solving cubic Schrödinger equations. The remainder of this paper is organized as follows. We first present some notations and preliminaries in Section 2. Then the scheme of energy-preserving exponential integrators is formulated in Section 3. Section 4 discusses the implementation issues. In Section 5, we analyse the existence, uniqueness and smoothness of the integrators. Section 6 pays attention to the regularity. The convergence of the integrators is studied in Section 7 and the nonlinear stability is discussed in Section 8. Section 9 is devoted to constructing one of practical exponential integrators in the light of the approach proposed in this paper, and reporting two numerical experiments to demonstrate the excellent qualitative behavior of the new exponential integrator. Section 10 focuses on the concluding remarks.

## 2 Notations and preliminaries

In this paper, we use the following notations.

• For , we denote For and , denotes the function on defined by

• We denote by (or simply ) the set of (classes of) complex functions on such that , endowed with the norm

• For all functions and all , denote by the Fourier coefficient

• For , denote by (or simply ) the space of (classes of) complex functions such that endowed with the norm Note that with the same norm.

• For all , we call a linear operator acting on as diagonal, if is in the domain of and there exists a such that . For example, the Laplace operator acting on is diagonal since for all , Another example is the identity operator on which will be denoted by .

• For all such operator and all functions from to itself, denote by the diagonal linear operator acting on whose domain is the set of linear combinations of functions defined for all by

• For all and all linear operators from to itself, we denote

The following result given in [21] will be useful for the analysis of this paper.

###### Proposition 1

(See [21]) With the above notations, if is a function from to bounded by on , then for all and , is a bounded linear operator from to itself with . For example, for all , it is true that

In order to derive the energy-preserving exponential integrators, we will use the idea of continuous time finite element methods in a generalised function space. To this end, we first present the following three definitions.

###### Definition 1

Define the generalised function space =span on as follows:

 Y={w:w(t,x)=r−1∑i=0φi(t)Ui(x), Ui(x):=(U1i(x)U2i(x)), U1i(x),U2i(x)∈Hα(Td)},

where is an integer satisfying and are supposed to be linearly independent on and sufficiently smooth. In this paper, we choose and as follows:

 Y=span{φ0(t),…,φr−1(t)}, X=% span{1,∫t0φ0(s)ds,…,∫t0φr−2(s)ds}

with for Then choose a time stepsize and define and on as

 Yh=span{~φ0(τ),…,~φr−1(τ)}, Xh=span{1,∫τ0~φ0(s)ds,…,∫τ0~φr−2(s)ds}, (3)

where for . It is noted that throughout this paper, the notations and are referred to as and for all the functions, respectively.

###### Definition 2

The inner product for the time is defined by

 ⟨~w1,~w2⟩=⟨~w1(τ,x),~w2(τ,x)⟩τ=∫10~w1(τ,x)⋅~w2(τ,x)dτ,

where and are two integrable functions for (scalar-valued or vector-valued) on , and if they are both vector-valued functions, ‘’ denotes the entrywise multiplication operation.

###### Definition 3

A projection onto is defined as

 ⟨~v(τ,x),Ph~w(τ,x)⟩=⟨~v(τ,x),~w(τ,x)⟩,for any  ~v(τ,x)∈Yh, (4)

where is a continuous two-dimensional vector function for .

With regard to the projection operation , the following property is important.

###### Lemma 1

The projection can be explicitly expressed as

 Ph~w(τ,x)=⟨Pτ,σ,~w(σ,x)⟩σ,

where and is a standard orthonormal basis of under the inner product given by

 ~ψj(t)=(−1)j√2j+1j∑k=0(jk)(j+kk)(−t)k

for

Proof The standard orthonormal basis of is obtained immediately from the choice of . Since , it can be expressed as By taking in (4) for and , we obtain

 ⟨~ψi(τ)ej,r−1∑k=0~ψk(τ)Uk(x)⟩=r−1∑k=0⟨~ψi(τ)ej,~ψk(τ)Uk(x)⟩ = r−1∑k=0⟨~ψi(τ),~ψk(τ)⟩(ej⋅Uk(x))=⟨~ψi(τ)ej,~w(τ,x)⟩,

which gives By the standard orthonormal basis , this result can be formulated as

 ⎛⎜ ⎜⎝U0(x)⋮Ur−1(x)⎞⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝⟨~ψ0(τ),~w(τ,x)⟩⋮⟨~ψr−1(τ),~w(τ,x)⟩⎞⎟ ⎟ ⎟⎠.

Then one has

 Ph~w(τ,x)=(~ψ0(τ),…,~ψr−1(τ))⊗I2×2⎛⎜ ⎜⎝U0(x)⋮Ur−1(x)⎞⎟ ⎟⎠ = (~ψ0(τ),…,~ψr−1(τ))⊗I2×2⎛⎜ ⎜ ⎜⎝⟨~ψ0(τ),~w(τ,x)⟩⋮⟨~ψr−1(τ),~w(τ,x)⟩⎞⎟ ⎟ ⎟⎠=⟨Pτ,σ,~w(σ,x)⟩σ,

which proves the result.

###### Remark 1

It is noted that the above three definitions and Lemma 1 can be considered as the generalised version of those presented in [39]. We also remark that one can make different choices of and , which will produce different methods by taking the methodology proposed in this paper.

## 3 Formulation of energy-preserving schemes

Define the linear differential operator by

 (Au)(x,t)=Δu(x,t)

and let . Then the system (1) is identical to

 ut=iAu+if(u),    u(0,x)=u0(x). (5)

The solutions of this equation satisfy Duhamel’s formula

 u(tn+h)=eihAu(tn)+i∫h0ei(h−ξ)Af(u(tn+ξ))dξ. (6)

Let and then the equation (1) can be rewritten as a pair of real initial-value problems

 (ptqt)=(0−AA0)(pq)+(λ(p2+q2)q−λ(p2+q2)p), (7) (p(0,x)q(0,x))=(\textmdRe(u0(x))\textmdIm(u0(x))).

In this case, the energy of this system is expressed by

 (8)

Accordingly, the system (7) can be formulated as the following infinite-dimensional real Hamiltonian system

 yt=Ky+g(y)=J−1δHδy,    y0(x)=(p(0,x)q(0,x)), (9)

where and

With the preliminaries described above, we first derive the energy-preserving exponential integrators for solving the real-valued equation (9) and then present the integrators for solving the cubic Schrödinger equation (1).

Find with , satisfying that

 ⟨~v(τ,x),~zt(τ,x)⟩−⟨~v(τ,x),K~z(τ,x)⟩=⟨~v(τ,x),g(~z(τ,x))⟩ (10)

for any , where From the definition (4) of the operator , it follows that

 ~zt(τ,x)−K~z(τ,x)=Ph(g(~z(τ,x))),

which yields the following result by considering Lemma 1

 ~zt(τ,x)=K~z(τ,x)+⟨Pτ,σ,g(~z(σ,x))⟩σ. (11)

Applying Duhamel’s formula (6) to (11) leads to

 ~z(τ,x)=eτhKy0(x)+τh∫10e(1−ξ)τhK⟨Pξτ,σ,g(~z(σ,x))⟩σdξ =eτhKy0(x)+τh∫10e(1−ξ)τhK∫10r−1∑i=0~ψi(ξτ)~ψi(σ)g(~z(σ,x))dσdξ =eτhKy0(x)+τh∫10r−1∑i=0∫10e(1−ξ)τhK~ψi(ξτ)dξ~ψi(σ)g(~z(σ,x))dσ.

This yields the following definition of energy-preserving exponential integrators for (9) in this paper.

###### Definition 4

The energy-preserving exponential integrator for solving the real Hamiltonian initial-value problem (9) is defined by

 ⎧⎪⎨⎪⎩~z(τ,x)=eτhKy0(x)+τh∫10¯Aτ,σ(K)g(~z(σ,x))dσ,0≤τ≤1,y1(x)=~z(1,x), (12)

where is a time stepsize and

 ¯Aτ,σ(K)=∫10e(1−ξ)τhKPξτ,σdξ=r−1∑i=0∫10e(1−ξ)τhK~ψi(ξτ)dξ~ψi(σ). (13)
###### Theorem 1

The continuous energy determined by (8) is preserved exactly by the integrator (12), i.e.,

 H(y1(x))=H(y0(x)).

Proof For each function , it follows from (10) that

 ∫10~w(τ,x)i~zt(τ,x)idτ=∫10~w(τ,x)i(K~z(τ,x))idτ+∫10~w(τ,x)ig(~z(τ,x))idτ,

for , where is the th vector of units and denotes the th entry of a vector. Then, one arrives at

 ∫10~w(τ,x)⊺~zt(τ,x)dτ=2∑i=1∫10~w(τ,x)i~zt(τ,x)idτ (14) = 2∑i=1[∫10~w(τ,x)i(K~z(τ,x))idτ+∫10~w(τ,x)ig(~z(τ,x))idτ] = ∫10~w(τ,x)⊺(K~z(τ,x)+g(~z(τ,x)))dτ.

Since , we obtain that and . Letting in (14) gives

 ∫10~zt(τ,x)⊺J⊺~zt(τ,x)dτ=∫10~zt(τ,x)⊺J⊺(K~z(τ,x)+g(~z(τ,x)))dτ.

Because is skew symmetric, we have

 0=∫10~zt(τ,x)⊺J⊺~zt(τ,x)dτ=−∫10~zt(τ,x)⊺J(K~z(τ,x)+g(~z(τ,x)))dτ.

Therefore, it is true that

 H(y1(x))−H(y0(x))=∫10∂∂τH(~z(τ,x))dτ=h∫10~zt(τ,x)⊺δH(~z)δ~zdτ =h∫10~zt(τ,x)⊺J(K~z(τ,x)+g(~z(τ,x)))dτ=h⋅0=0.

On the basis of the analysis stated above, we now present the energy-preserving scheme for solving the cubic Schrödinger equation (1).

###### Definition 5

The energy-preserving exponential integrator (denoted as EPEIr) with a time stepsize for the cubic Schrödinger equation (1) is defined by

 ⎧⎪⎨⎪⎩~u(τ,x)=eτhiAu0(x)+τh∫10¯Aτ,σ(iA)if(~u(σ,x))dσ,0≤τ≤1,u1(x)=~u(1,x), (15)

where

 ¯Aτ,σ(iA)=∫10e(1−ξ)τhiAPξτ,σdξ=r−1∑j=0∫10e(1−ξ)τhiA~ψj(ξτ)dξ~ψj(σ). (16)
###### Theorem 2

The integrator (15) can exactly preserve the continuous energy of the cubic Schrödinger equation (1).

Proof  In what follows, we prove that the schemes (15) and (12) share the same result for solving (1). With careful calculations, we obtain

 eiAu0(x)= [Id+∞∑k=0(−A4k+2(4k+2)!+A4k+4(4k+4)!)]\textmdRe(u0(x)) −∞∑k=0(A4k+1(4k+1)!+−A4k+3(4k+3)!)\textmdIm(u0(x)) +i[Id+∞∑k=0(−A4k+2(4k+2)!+A4k+4(4k+4)!)]\textmdIm(u0(x)) +i∞∑k=0(A4k+1(4k+1)!+−A4k+3(4k+3)!)\textmdRe(u0(x)).

On the other hand, we have

 eK(\textmdRe(u0(x))\textmdIm(u0(x))) =

These formulae show that the above two different expressions share the same result for solving (1). Moreover, it is clear that and then and yield the same result. According to the above analysis and looking closer at the two schemes (15) and (12), it can be confirmed that they are the same integrator when applied to (1).

## 4 Implementation issues

Introduce

 li(τ)=r∏j=1,j≠iτ−djdi−dj,   i=1,2…,r

with respect to distinct points . Then is a basis of satisfying It is assumed that since . Choosing for and using , one gets which leads to Let and denote . The integrator (15) now becomes

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩yσ=r∑i=1ydili(σ),ydi=edihiAy0(x)+dih∫10¯Adi,σ(iA)if(yσ)dσ,i=1,…,r,y1(x)=ehiAy0(x)+h∫10¯A1,σ(iA)if(yσ)dσ. (17)

It is clear that in (17), inserting the appearing in the two integrals by the first equation yields the scheme of EPEIr integrators. Consequently, we have the following definition for energy-preserving exponential integrators in practice.

###### Definition 6

The practical EPEIr integrator for integrating (1) is defined by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Yn,i=edihiAyn+diih∫10¯Adi,σ(iA)f(r∑j=1Yn,jlj(σ))dσ,i=1,…,r,yn+1=ehiAyn+ih∫10¯A1,σ(iA)f(r∑j=1Yn,jlj(σ))dσ, (18)

where is time stepsize.

## 5 The existence, uniqueness and smoothness

In the remainder of this paper, we use the following assumptions on the exact solution of (1) and on the non-linearity .

###### Assumption 1

It is assumed that the Schrödinger equation (1) admits an exact solution which is sufficiently smooth. In particular, there exists such that for all , where

 B(¯u0,R)={u∈Hα(Td):||u−¯u0||Hα≤R}

and .

###### Assumption 2

We assume that the mapping is sufficiently smooth.

###### Assumption 3

The th-order derivative of is assumed that . And we denote for

Note that these assumptions are fulfilled in the case of the cubic non-linear Schrödinger equation (see [4] for example, Chapter II, Proposition 2.2 ), at least when and .

According to Proposition 1, one gets that the coefficients and of our integrators for and are uniformly bounded. Hence, we let

 Mk=maxτ,σ,h∈[0,1]∥∥ ∥∥∂k¯Aτ,σ∂hk∥∥ ∥∥Hα→Hα, Ck=maxτ,h∈[0,1]∥∥∥∂keτhiA∂hku0(x)∥∥∥Hα, k=0,1. (19)
###### Theorem 3

Under the above assumptions, if satisfies

 0≤h≤κ

then the EPEIr integrator (15) admits a unique solution which smoothly depends on .

Proof  By setting and defining

 ~un+1(τ,x)=eτhiAu0(x)+τh∫10¯Aτ,σ(iA)ig(~un(σ,x))dσ,n=0,1,…, (21)

we get a function series If is uniformly convergent, is a solution of the EPEIr integrator (15).

By induction, it follows from (20) and (21) that for Using (21), we have

 ||~un+1(τ,x)−~un(τ,x)||Hα≤τh∫10M0D1||~un(σ,x)−~un−1(σ,x)||Hαdσ ≤h∫10M0D1||~un(σ,x)−~un−1(σ,x)||Hαdσ≤β||~un−~un−1||c,β=κM0D1,

where is the maximum norm for continuous functions defined as for a continuous function with on . Hence, one arrives at

 ||~un+1(τ,x)−~un(τ,x)||c≤β||