Optimal error estimate of a compact scheme for nonlinear Schrödinger equation

Applied Numerical Mathematics 120 (2017) 68–81

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Optimal error estimate of a compact scheme for nonlinear Schrödinger equation ✩ Jialin Hong a , Lihai Ji b , Linghua Kong c , Tingchun Wang d,∗ a

Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China b Institute of Applied Physics and Computational Mathematics, Beijing 100094, China c School of Mathematics and Information Science, Jiangxi Normal University, Nanchang, Jiangxi, 330022, China d School of Mathematics and Statistics, Nanjing University of Information Science & Technology, Nanjing, 210044, China

a r t i c l e

i n f o

Article history: Received 20 March 2016 Received in revised form 9 November 2016 Accepted 8 May 2017 Available online 15 May 2017 Keywords: Nonlinear Schrödinger equation Compact scheme Symplectic scheme Energy conservation Optimal error estimate

a b s t r a c t It has been pointed out in literature that the symplectic scheme of a nonlinear Hamiltonian system can not preserve the total energy in the discrete sense Ge and Marsden (1988) [10]. Moreover, due to the difficulty in obtaining a priori estimate of the numerical solution, it is very hard to establish the optimal error bound of the symplectic scheme without any restrictions on the grid ratios. In this paper, we develop and analyze a compact scheme for solving nonlinear Schrödinger equation. We introduce a cut-off technique for proving optimal L ∞ error estimate for the compact scheme. We show that the convergence of the compact scheme is of second order in time and of fourth order in space. Meanwhile, we define a new type of energy functional by using a recursion relationship, and then prove that the compact scheme is mass and energy-conserved, symplectic-conserved, unconditionally stable and can be computed efficiently. Numerical experiments confirm well the theoretical analysis results. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The nonlinear Schrödinger (NLS) equation describes many phenomena and has important applications in fluid dynamics, nonlinear optics and quantum mechanics [3,20]. Various kinds of numerical methods can be found for simulating solutions of the NLS equation. In [4,5,17], time-splitting pseudo-spectral methods are investigated. In [16,22,23], several important numerical methods are tested, analyzed and compared. The discontinuous Galerkin methods are considered in [13,27]. For a detailed description of various numerical methods for NLS equation as well as its implementation and applications, we refer the readers to the orthogonal spline collocation method [15,19] and the Runge–Kutta or Crank–Nicolson pseudo-spectral method [8,18] and the references therein. Based on the basic rule that numerical algorithms should preserve the intrinsic properties of the original problems as much as possible, Feng [9] first presented the concept of symplectic schemes for Hamiltonian systems and developed the structure-preserving algorithms for the general conservative dynamical systems. In [6,26], the authors propose several conservative finite difference schemes for NLS equation and present that these schemes preserve the total mass and energy

✩

*

Authors are supported by NNSFC (No. 11571181, No. 91530118, No. 11471310, No. 91630312, No. 11021101, No. 11290142, No. 11601032). Corresponding author. E-mail addresses: [email protected] (J. Hong), [email protected] (L. Ji), [email protected] (L. Kong), [email protected] (T. Wang).

http://dx.doi.org/10.1016/j.apnum.2017.05.004 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

69

in the discrete sense. Ref. [11] firstly studied the standard Crank–Nicolson finite difference scheme which is a symplectic scheme [7,21], and obtained an error estimate of the numerical solution under an rigorous condition τ /h2 ≤ o(1) with time-step τ and mesh-size h. However, the symplectic finite difference schemes for nonlinear Hamiltonian systems can not preserve the total energy in the discrete sense [10]. Furthermore, it is very difficult to obtain a priori estimate of the numerical solution of symplectic schemes in L ∞ -norm. To the best of our knowledge, there is no reference about symplectic and energy conservative numerical methods for NLS equation till now. In this work, we proposed a novel compact scheme for NLS equation. High-order compact technique is popular due to its high resolution, compactness and economy in scientific computation [14]. The proposed scheme can not only preserve the discrete symplectic structure, but also conserves the total mass and energy conservation laws. And the convergence result is analyzed, which shows that the compact scheme is of second order in time and of fourth order in space. Special techniques are utilized to achieve the optimal order of accuracy. Numerical experiments are shown to demonstrate the accuracy and capability of the proposed compact scheme. The rest of this paper is organized as follows. In Section 2, we present some preliminary results and introduce several useful lemmas. In Section 3, a novel compact scheme for NLS equation is proposed and we show that the scheme preserves the discrete symplectic conservation law, discrete mass conservation law and discrete energy evolution law. The optimal error estimate in maximum norm is presented in Section 4. In Section 5, numerical experiments are performed to testify the effectiveness and accuracy of the scheme. Finally, some conclusions are addressed in Section 6. 2. Preliminary results In order to simplify the notations, in this paper we consider one-dimensional NLS equation. The NLS equation of our interest is

iut + α u xx + V (x)u + f (|u |2 )u = 0, x ∈ , t ∈ (0, T ],

(2.1)

with an initial condition

u (x, 0) = ϕ (x)

(2.2)

and homogeneous boundary conditions. Here,  = [a, b], V (x) and f (u ) are arbitrary (smooth) nonlinear real functions and α is a real constant. It can be proved that the NLS equation (2.1) possesses the mass conservation law,

b |u (x, t )|2 dx = Constant,

Q (t ) :=

(2.3)

a

and energy conservation law

E (t ) :=

b 



α |∇ u |2 − V (x)|u |2 − F (|u |2 ) dx = Constant,

(2.4)

a

where

q F (q) =

f ( p )dp . 0

If we set u = r + is, where r, s are real-valued functions, we can separate (2.1) into the following form

rt + α sxx + V (x)s + f (r 2 + s2 )s = 0,

−st + α r xx + V (x)r + f (r 2 + s2 )r = 0. By introducing two additional new variables, p = s x , q = r x , and defining a state variable z = (r , s, p , q) T , the equation above can be transformed to the Hamiltonian system

K zt + Lzx = ∇z S (x, z), where

⎛

⎞

(2.5)

⎛

0 1 0 0 0 ⎜ −1 0 0 0 ⎟ ⎜0 ⎟ ⎜ K =⎜ ⎝ 0 0 0 0 ⎠, L = ⎝ 0 0 0 0 0 α and

0 0

α 0

0

−α 0 0

−α

⎞

0 ⎟ ⎟, 0 ⎠ 0

70

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

S (x, z) = V (x)

r 2 + s2

+α

2

p 2 + q2 2

F (r 2 + s2 )

+

2

.

Noticing that equation (2.5) has the multi-symplectic conservation law

∂t ω + ∂x κ = 0,

(2.6) 1 dz 2

1 dz 2

with ω , κ being the two forms, i.e., ω = ∧ K dz, κ = ∧ Ldz. The multi-symplectic conservation law means that symplecticity changes locally and synchronously both in time and in space directions. Now for the NLS equation (2.1), there are three conservation laws. Therefore it is natural to ask whether the numerical schemes can keep these properties. We introduce a uniform grid (x j , t n ) ∈ [a, b] × [0, T ] with mesh length h in spatial direction and mesh length τ in temporal direction, with j = 0, 1, · · · , J , n = 0, 1, · · · , N. The approximate value of the function u (x, t ) at the mesh point (x j , t n ) is denoted by unj . Denote the space

X h := {u = (u 0 , u 1 , u 2 , · · · , u J ) T | u 0 = u J = 0} ⊆ C J +1 , X h0 := {(u 1 , u 2 , · · · , u J −1 ) T | u ∈ X h } ⊆ C J −1 , where C J +1 means the ( J + 1)-dimensional complex space. For convenience, we introduce the following operators

δx+ unj = δt+ unj =

unj+1 − unj h

δx2 unj =

,

unj+1 − 2unj + unj−1 h2

,

unj +1 − unj

, un = (un1 , un2 , · · · , unJ −1 )T ,  unj−1 + 10unj + unj+1 h2 2 A h unj = 1 + δx unj =

τ

12

12

and two tridiagonal matrices

⎛

10 ⎜ 1 ⎜

1 ⎜ H= ⎜ 12 ⎜

⎞

1 10

1

..

..

⎝

.

..

.

1

2 −1 ⎜ −1 2 ⎜

⎟ ⎟ 1 ⎜ ⎟ ⎟, D = 2 ⎜ ⎟ h ⎜ ⎠ ⎝

.

10 1

⎛

1 10

..

.

⎞ −1 .. . −1

..

.

2 −1

−1

⎟ ⎟ ⎟ ⎟. ⎟ ⎠

2

Remark 1. It can be proved that H and D are two symmetric and positive definite matrices. For u , v ∈ X h , we define the discrete inner product and associated norms

u , v := h

J −1

u j v j , ||δx+ u || := h

j =1

J −1

|δx+ u j |2 ,

j =0

⎛

||u ||∞ := max |u j |, ||u || p := ⎝h 1 < j < J −1

J −1

⎞1

p

|u j |

p⎠

, 1 < p < ∞,

j =1

where v j is the conjugate of v j . If there is no confusion, we may write || v ||2 as || v ||. Throughout the paper, we adopt the standard Sobolev spaces and their corresponding norms, and use the notation w  v to present that there is a generic constant C which is independent of time step τ and mesh size h such that w ≤ C v. First of all, we introduce a useful lemma. Lemma 1. Matrices H and D hold the following properties. (a) The eigenvalues of the matrices H and D are

λ H ,k =

1



10 + 2cos

12

kπ



J

1

, λ D ,k =



h2

(b) The two matrices have the same eigenvectors, i.e.,



v k = sin

kπ J

, sin 1

2kπ J 1

, · · · , sin

(c) H −1 D = D H −1 , H −1 D 2 = D 2 H −1 .

( J − 1)kπ J

2 + 2cos

kπ J

, k = 1, · · · , J − 1.

(2.7)

T , k = 1, · · · , J − 1.

(2.8)

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

71

Define the semi-norm of u ∈ X h as follows

|||δx+ u ||| := H −1 Du, u , then it holds Lemma 2. [25] For any grid function u ∈ X h

√

||δx+ u || ≤ |||δx+ u ||| ≤

6

2

||δx+ u ||.

(2.9)

Lemma 3. [25] For any grid function un ∈ X h , n = 0, 1, 2, · · · , N

Re H −1 D (un+1 + un ), (un+1 − un ) = |||δx+ un+1 |||2 − |||δx+ un |||2 ,

(2.10)

where Re( f ) denotes the real part of f . Lemma 4. [1] For any grid function u ∈ X h , it holds

√

b − a  +  δ u  . x 2

||u ||∞ ≤

(2.11)

3. A novel compact scheme With the definition of the above difference operators and matrices, we propose the following compact scheme for (2.1) n+ 12

i δt+ unj + α A h−1 δx2 u j n+ 12

where u j

n+ 12

+ V ju j

n+ 12 2

+ f (|u j

n+ 12

| )u j

= 0,

(3.1)

= 12 (unj + unj +1 ).

3.1. Symplectic structure In order to investigate the symplectic conservation law, we rewrite (3.1) as 1

1

1

1

i δt+ un − α H −1 Dun+ 2 + Vun+ 2 + f (|un+ 2 |2 )un+ 2 = 0,

(3.2)

with



1

 n+ 1 n+ 1 T , V 2 u 2 2 , · · · , V J −1 u J −12 ,   n+ 1 n+ 1 n+ 1 n+ 1 T = f (|u 1 2 |2 )u 1 2 , · · · , f (|u J −12 |2 )u J −12 .

n+ 12

Vun+ 2 = V 1 u 1 1

1

f (|un+ 2 |2 )un+ 2

1

Assume un = rn + isn , where r, s are real-valued functions. Based on Lemma 1 and introducing new variables pn = D 2 sn , 1 2

qn = D rn , equation (3.2) can be written as the following first-order system

δt+ rn − α H −1 D 2 pn+ 2 + Vsn+ 2 + f (|rn+ 2 |2 + |sn+ 2 |2 )sn+ 2 = 0, 1

1

1

1

1

1

1

1

1

D 2 sn+ 2 = pn+ 2 , 1

1

1

1

1

1

− δt+ sn − α H −1 D 2 qn+ 2 + Vrn+ 2 + f (|rn+ 2 |2 + |sn+ 2 |2 )rn+ 2 = 0, 1

1

(3.3)

1

D 2 rn+ 2 = qn+ 2 .



Let Z n = (rn ) T , (sn ) T , (pn ) T , (qn ) T

T

∈ R4( J −1) , we can write (3.3) as

Kδt+ Z n + LD Z n+ 2 = ∇ Z S ( Z n+ 2 ), 1

1

(3.4)

where 1

S ( Z n+ 2 ) =

 1 1 1 1 1  n+ 1 2 F |r 2 | + |sn+ 2 |2 + V(|rn+ 2 |2 + |sn+ 2 |2 ) 2 2 +

α 2

1

1

1

1

(pn+ 2 )T H −1 pn+ 2 + (qn+ 2 )T H − qn+ 2



72

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

and

 K=

with

A 0

 A =

0 −I

0 0



I 0

 4 ( J −1 )

0

, L=

−B



 2 ( J −1 )

, B=

B

⎛



0

0 α H −1

4 ( J −1 )

α H −1 0

⎜ , D=⎜ ⎝

⎞

1

D2

..

. 1

D2

⎟ ⎟ ⎠ 4 ( J −1 )

. 2 ( J −1 )

It is easy to rewrite the above system into

Kδt+ Z n = ∇ Z S˜ ( Z n+ 2 ), 1

with S˜ ( Z ) = S˜ ( Z ) −

1 2

(3.5)

Z T LD Z . The variational equation for (3.5) reads

Kδt+ d Z n = ∇ 2Z S˜ ( Z n+ 2 )d Z n+ 2 , 1

1

1

then taking the wedge product with d Z n+ 2 , we get

d Z n+1 ∧ Kd Z n+1 − d Z n ∧ Kd Z n

τ

= 0.

Consequently, one obtains the following theorem: Theorem 5. The phase flow of (3.1) preserves symplectic structure

ω n +1 − ω n = 0, where ωn =

1 dZn 2

(3.6)

∧ Kd Z n .

3.2. Conservation laws and unconditionally stable This subsection investigates the global conservative properties of the compact scheme (3.1). Notice that

s2 F (s2 ) − F (s1 ) =

f (s)ds s1

and using numerical integration formula, we obtain |unj +1 |2

F (|unj +1 |2 ) = F (|unj |2 ) +



f (s)ds |unj |2 n+ 12 2

= F (|unj |2 ) + f (|u j

| )(|unj +1 |2 − |unj |2 ) + O (τ 3 ).

Define

F 0j = F (|ϕ j |2 ),

(3.7)

n+ 12 2

F nj +1 = F nj + f (|u j

| )(|unj +1 |2 − |unj |2 ).

Theorem 6. Under the homogeneous boundary conditions, the compact scheme (3.1) satisfies the discrete mass and energy conservation laws, i.e.,

Q n = ||un ||2 ≡ Q 0 , E n = − α |||δx+ un |||2 + h

J −1

j =1

n

n

V j |unj |2 + h

J −1

F nj ≡ E 0 ,

j =1

where Q , E denote the discrete mass and energy at time-step t n , respectively.

(3.8)

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

73

Proof. Multiply equation (3.2) by un+1 + un , i.e., the conjugate of un+1 + un , sum over all spatial grid points j, and take the imaginary part. Then, under the homogeneous boundary conditions

1

τ

 ||un+1 ||2 − ||un ||2 = 0,

(3.9)

we obtain the discrete mass conservation law. In the similar way, we multiply equation (3.2) by un+1 − un , sum over all spatial grid points j, and take the real part. Then, under the homogeneous boundary conditions

−α |||δx+ un+1 |||2 + α |||δx+ un |||2 + h

J −1

V j (|unj +1 |2 − |unj |2 )

j =1 J −1

+h

n+ 12 2

f (|u j

| )(|unj +1 |2 − |unj |2 ) = 0.

(3.10)

j =1

It follows from Lemma 3 and equation (3.7), we get the discrete energy conservation law. The results of this theorem are evidently consistent with (2.3) and (2.4), which means that the mass and energy conservation laws can be exactly preserved by the proposed compact scheme. Combining the above mass conservation property, we know that the compact scheme (3.1) is unconditionally stable. Corollary 1. The compact scheme for NLS equation with homogeneous boundary conditions is unconditionally stable. 4. Error estimate In this section, we derive optimal error estimate for the compact scheme under the discrete L ∞ norm. We need the following well-known result called the discrete Gronwall’s lemma (see [12]). Lemma 7. If {ak }, {bk } and {ck } are three positive sequences with {ck } being monotone such that ak + bk ≤ ck + λ and a0 + b0 ≤ c 0 then ak + bk ≤ ck exp(kλ) for all k ≥ 0. Define the local truncation error

k −1 i =0

ai with λ > 0

ηn ∈ X h of the compact scheme (3.1) as 1

1

ηnj := i δt+ u (x j , t n ) + α A h−1 δx2 u (x j , t n+ 2 ) + V j u (x j , t n+ 2 ) 1

1

+ f (|u (x j , t n+ 2 )|2 )u (x j , t n+ 2 ).

(4.1)

By Taylor expansion, we can get Lemma 8. Assume u ∈ C 6,4 ([a, b] × [0, T ]), V ∈ C ([a, b]), then the truncation errors can be bounded by

||ηn ||∞  τ 2 + h4 ,

(4.2)

||δt+ ηn ||∞  τ 2 + h4 . Based on the truncation error analysis, we arrive at the convergence analysis of our compact scheme (3.1).

Theorem 9. Assume u ∈ C 6,4 ([a, b] × [0, T ]), V ∈ C ([a, b]), then numerical solution un of the scheme (3.1) converges to the exact solution u (·, t n ) of the initial–boundary value problem (2.1)–(2.2) in the L ∞ norm at the order of O (h4 + τ 2 ), i.e.,

||u (·, t n ) − un ||∞  τ 2 + h4 , n = 0, 1, · · · , N − 1. Proof. As the nonlinear term is not global Lipschitz, it is very difficulty to establish a priori estimate of the numerical solution. To overcome this difficulty, we will use the truncation technique (see [2,24] and references therein) to obtain a truncated equation with global Lipschitz nonlinear term. Firstly, we define a cut-off function ρ ∈ C ∞ (R) as

ρ (s) =

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩

1,

| s | ≤ 1,

∈ [0, 1],

1 ≤ | s | ≤ 2,

0,

| s | ≥ 2,

(4.3)

74

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

and denote

B = ( M 0 + 1)2 , M 0 = max ||u (·, t )||∞ , f B (s) := ρ

s

t ∈[0, T ]

B

f (s), s ∈ R,

then f B ∈ C ∞ (R), and f B is globally Lipschitz continuous, i.e.,

√ √ | f B (s1 ) − f B (s2 )| ≤ C B | s1 − s2 |,

∀s1 , s2 ∈ [0, ∞).

(4.4)

ˆn

Define u ∈ X h as n+ 12

n+ 12 2

ˆ n+ 2 ) j + V j uˆ j i δt+ uˆ nj − α ( H −1 D u 1

+ f B (|uˆ j

n+ 12

| )uˆ j

= 0,

(4.5)

then uˆ nj can be viewed as another approximation of u (x j , t n ). Define the error eˆ n ∈ X h of the numerical solution of (4.5) as

eˆ nj = u (x j , t n ) − uˆ nj . And define the local truncation error σˆ n ∈ X h of (4.5) as

σˆ jn = i δt+ u (x j , t n ) − α ( H −1 Du) j (x j , t n+ 2 ) + V j u (x j , t n+ 2 ) 1

1

1

(4.6)

1

+ f B (|u (x j , t n+ 2 )|2 )u (x j , t n+ 2 ). From the definition of the function f B , we know that





1

1

f B |u (x j , t n+ 2 )|2 = f (|u (x j , t n+ 2 )|2 ), then, from Lemma 8, it holds

|σˆ jn |  τ 2 + h4 , +

|δt ˆ jn |

2

(4.7) 4

σ τ +h .

(4.8)

Subtracting (4.5) from (4.6) yields n+ 12

1

σˆ jn = i δt+ eˆ nj − α ( H −1 D uˆ n+ 2 ) j + V j eˆ j

n+ 12

+ ξˆ j

(4.9)

,

where 1

1

1

n+ n+ n+ ξˆ j 2 = f B (|u (x j , t n+ 2 )|2 )u (x j , t n+ 2 ) − f B (|uˆ j 2 |2 )uˆ j 2 . 1

1

(4.10)

Then, we have

   n+ 12   n+1   n  ξˆ   j   eˆ j  + eˆ j  .

(4.11)

For any 0 ≤ j ≤ J − 1, 0 ≤ n ≤ N and θ ∈ [0, 1], denote

u (x j , t n , θ) = θ u (x j +1 , t n ) + (1 − θ)u (x j , t n ),

unj,θ = θ unj+1 + (1 − θ)unj ,

we have n+ 12

δx+ ξˆ j

1

n+ 12 2

1

= δx+ [ f B (|u (x j , t n+ 2 )|2 )u (x j , t n+ 2 )] − δx+ [ f B (|uˆ j . =

n+ 1 Pj 2

−

n+ 12

| )uˆ j

]

n+ 1 Q j 2,

where n+ 12

Pj

=

1  1 1 1 1 ( f B (|u (x j , t n+ 2 , θ)|2 ) + f B (|u (x j , t n+ 2 , θ)|2 )|u (x j , t n+ 2 , θ)|2 )δx+ u (x j , t n+ 2 ) 0

n+ 12

Qj

 1 1 1 + f B (|u (x j , t n+ 2 , θ)|2 )(u (x j , t n+ 2 , θ))2 δx+ u (x j , t n+ 2 ) dθ,

1 

n+ 1

n+ 1

n+ 1

n+ 12

( f B (|uˆ j ,θ 2 |2 ) + f B (|uˆ j ,θ 2 |2 )|uˆ j ,θ 2 |2 )δx+ uˆ j

= 0

From the definition of f B , we know that f B ∈ C 02 ( R ) and

n+ 1

n+ 1

n+ 12

+ f B (|uˆ j ,θ 2 |2 )(uˆ j ,θ 2 )2 δx+ uˆ j

 dθ.

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

75

 1 1 1   ( f B (|u (x j , t n+ 2 , θ)|2 ) + f B (|u (x j , t n+ 2 , θ)|2 )|u (x j , t n+ 2 , θ)|2   1  n+ 1 n+ 1 n+ 1 − ( f B (|uˆ j ,θ 2 |2 ) + f B (|uˆ j ,θ 2 |2 )|uˆ j ,θ 2 |2 δx+ u (x j , t n+ 2 )    n+ 12  n+ 12   |u (x j , t , θ)| − |uˆ j ,θ | n+ 12

n+ 1

+1 | + |ˆe j +12 |  |ˆenj | + |ˆenj+1 | + |ˆenj +1 | + |ˆenj + |, 1    1 1 1  n+ 1 n+ 1    f B (|u (x j , t n+ 2 , θ)|2 )(u (x j , t n+ 2 , θ))2 − f B (|uˆ j ,θ 2 |2 )(uˆ j ,θ 2 )2 δx+ u (x j , t n+ 2 )    n+ 12  n+ 12   |u (x j , t , θ)| − |uˆ j ,θ |

 |ˆe j

n+ 12

n+ 1

+1 | + |ˆe j +12 |  |ˆenj | + |ˆenj+1 | + |ˆenj +1 | + |ˆenj + |, 1    1 1 1 1    f B (|uˆ n+ 2 |2 ) + f  (|uˆ n+ 2 |2 )|uˆ n+ 2 |2 (δ + u (x j , t n+ 12 ) − δ + uˆ n+ 2 )  |δ + eˆ n | + |δ + eˆ n+1 |, x x j x j x j B j ,θ j ,θ j ,θ  

 |ˆe j

    n+ 12 2 n+ 12 2 + n+ 12  + n+ 12  f (|uˆ ) − δx uˆ j )  |δx+ eˆ nj | + |δx+ eˆ nj +1 |.  B j ,θ | )(uˆ j ,θ ) (δx u (x j , t Thus, we have n+ 12

|δx+ ξˆ j

+1 |  |ˆenj | + |ˆenj+1 | + |ˆenj +1 | + |ˆenj + | + |δx+ eˆ nj | + |δx+ eˆ nj +1 |. 1

(4.12)

Computing the inner product of (4.9) with eˆ n+1 + eˆ n , then taking the imaginary part, we obtain

δt+ ||ˆen ||2 + Im(ξ n+ 2 , eˆ n+1 + eˆ n ) = Im(σˆ n , eˆ n+1 + eˆ n ). 1

(4.13)

By using Cauchy–Schwarz inequality and Lemma 8,

  ||ˆen+1 ||2 − ||ˆen ||2  τ ||ˆen+1 ||2 + ||ˆen ||2 + (h4 + τ 2 )2 , thus from Lemma 7, there exists sufficiently small

||ˆen ||  τ 2 + h4 ,

(4.14)

τ0 such that for τ ≤ τ0

n = 1, 2, · · · , N .

(4.15)

Computing the product of (4.9) with eˆ n+1 − eˆ n , then taking the real part, we have

|||δx+ eˆ n+1 |||2 − |||δx+ eˆ n |||2 − J −1

=

h

α

J −1 h

α

V j (|ˆenj +1 |2 − |ˆenj |2 )

j =1

n+ 1 ξ j 2 (ˆenj +1

− eˆ nj ) −

j =1

J −1 h

α

(4.16)

ˆ jn (ˆenj +1

σ

− eˆ nj ).

j =1

Summing (4.13) for n from 0 to m, and replacing m by n, it yields

|||δx+ eˆ n+1 |||2 +

J −1 h

α

V j |ˆenj +1 |2

j =1

J −1

=−

J −1 n n 1 1 l+ 1 h Re[ξ j 2 (ˆelj+1 − eˆ lj )] + h Re[σˆ lj (ˆelj+1 − eˆ lj )].

α

l =0

α

j =1

l =0

(4.17)

j =1

From (4.9), we know n+ 12

eˆ nj +1 − eˆ nj = τ [−i α ( H −1 D eˆ n+ 2 ) j + i V j eˆ j 1

1

n+ + i ξˆ j 2 − i σˆ jn ].

Substituting (4.18) into the first term of the right hand of (4.17) yields

(4.18)

76

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

J −1 n 1 l+ 1 h Re[ξˆ j 2 (ˆelj+1 − eˆ lj )]

α

=−

=−

=τ

j =1

l =0

τ α τ α

n

h

J −1

l =0

j =1

n

J −1

h

h





1



l+ Im ξˆ j 2

l+ 12

−α ( H −1 D eˆ l+ 2 ) j + V j eˆ j 1

l+ 12

−α ( H −1 D eˆ l+ 2 ) j + V j eˆ j 1



1

1

1

1

1

Im ( H − 2 D 2 ξˆ l+ 2 ) j ( H − 2 D 2 eˆ l+ 2 ) j

j =1

l =0

1

l+ 12

+ ξˆ j

 − σˆ jl

 − σˆ jl

(4.19)

j =1

l =0

J −1 n

J −1



1





τ l+ 1 l+ 1 h Im ξˆ j 2 ( V j eˆ j 2 − σˆ jl ) . α j =1 n

−



l+ Im ξˆ j 2

l =0

By using Cauchy–Schwarz inequality and Lemma 8, we obtain

   n J −1  1    ˆ l+ 2 (ˆel+1 − eˆ l )] h Re [ ξ j j j    l =0 j =1       J −1   n 1 1 1 1 1 1 − l + − h Im ( H 2 D 2 ξˆ 2 ) j ( H 2 D 2 eˆ l+ 2 ) j   τ   l =0 j =1     1  J −1   n l+ l+ 1 h Im ξˆ j 2 ( V j eˆ j 2 − σˆ jl )  + τ   l =0 j =1 τ

n 

|||δx+ ξˆ l |||2 + |||δx+ ξˆ l+1 |||2 + |||δx+ eˆ l |||2 + |||δx+ eˆ l+1 |||2

l =0

+ ||ξˆ l ||2 + ||ξˆ l+1 ||2 + ||ˆel ||2 + ||ˆel+1 ||2 + ||σˆ l ||2 τ

(4.20)



n 

 ||δx+ eˆ l ||2 + ||δx+ eˆ l+1 ||2 + (h4 + τ 2 )2 .

l =0

For the last term of (4.17), by using the summation by parts formula in time direction and Lemma 8, we obtain

     J −1 n    J −1 J −1 n −1

 l n +1    + l l +1 n n +1  n h   σˆ j (ˆe j − eˆ j ) = −τ h δt σˆ j eˆ j + h σˆ j eˆ j    j =1 l =0    j =1 j =1 l =0

(4.21)

 (h4 + τ 2 )2 . Substituting (4.20) and (4.21) into (4.17) yields

|||δx+ eˆ n+1 |||2  τ

n 

 |||δx+ eˆ l |||2 + |||δx+ eˆ l+1 |||2 + (h4 + τ 2 )2 .

(4.22)

l =0

By using Lemma 7, there exists a sufficiently small + ˆn

2

4

|||δx e |||  τ + h , for

τ0 such that

n = 1, 2, · · · , N ,

(4.23)

τ ≤ τ0 . This together with Lemma 3 gives ||δx+ eˆ n ||  τ 2 + h4 ,

n = 1, 2, · · · , N .

(4.24)

By using the discrete version of Sobolev inequality [1], we obtain

||ˆen ||∞  ||δx+ eˆ n ||  h4 + τ 2 , n = 1, 2, · · · , N , therefore,

(4.25)

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

||uˆ n ||∞ ≤ ||un ||∞ + ||ˆen ||∞ ≤ M 0 + C (h4 + τ 2 ), n = 1, 2, · · · , N . This tells us that, if h and

77

(4.26)

τ are small enough such that C (h4 + τ 2 ) ≤ 1, we have

||uˆ n ||∞ ≤ B , n = 1, 2, · · · , N .

(4.27)



n+ 12

|uˆ j

From the definition of function f B , we know that f B



n+ 12

|2 = f (|uˆ j

|2 ). At this time, the scheme (4.5) becomes the

initial scheme (3.1), and we further obtain

un = uˆ n ,

en = eˆ n .

(4.28)

This completes the proof. Remark 2. The theoretic results in this paper are also valid for the periodic boundary–initial value problem of NLS equation. Furthermore, in the case of periodic boundary condition, matrices H and D become

⎛

10 ⎜ 1 ⎜

1 ⎜ H= ⎜ 12 ⎜

1 10

1

..

..

⎝

1

.

.

1 1

..

.

10 1

⎞

⎛

2 −1 ⎜ −1 2 ⎜

⎟ ⎟ 1 ⎜ ⎟ ⎟, D = 2 ⎜ ⎟ h ⎜ ⎝ 1 ⎠

..

−1

10

.

−1 .. . −1

−1 ..

.

2 −1

⎞

⎟ ⎟ ⎟ ⎟. ⎟ −1 ⎠ 2

5. Numerical results Just as we have mentioned in the beginning of the paper and proven in the preceding sections, our new compact scheme has some nice advantages. In this section, we provide numerical examples to illustrate the accuracy and capability of the method developed in the previous sections. In each sub-interval [t n , t n+1 ], under the initial condition ϕ (x) and homogeneous boundary conditions, we write (3.1) as the form

A (n)U n+1 = B (n)U n + F (tn , tn+1 , U n , U n+1 ), n = 1, 2, · · · , where A (n), B (n) are invertible tridiagonal matrices depending on coefficients of the equation, the vector U n = (un1 , un2 , · · · , unJ −1 )T and F denotes the discretization of nonlinear term. For numerical procedure, we utilize the follow-

ing fixed-point iteration method to solve the nonlinear term

U [n0+] 1 = U n ,





U [nk++11] = A −1 B (n)U n + F (tn , tn+1 , U n , U [nk+] 1 ) , k = 1, 2, · · · . Example 5.1. Consider the following NLS equation

i ∂t u − ∂xx u + x2 u + |u |4 u = 0,

(5.1)

with the initial data



u (x, 0) = exp −

x2  2

.

In this numerical example, we take the temporal step-size τ = 0.05, the spatial mesh grid-size h = 0.05, and the time interval [0, 10], the numerical spatial domain [a, b] = [−30, 30]. As is stated in Theorem 6, the compact scheme (3.1) preserves both the discrete mass conservation law and discrete energy conservation law. We consider this phenomenon numerically in Fig. 5.1 and Fig. 5.2. From the figures, we can see that the discrete mass and energy conservation laws remain to be horizontal lines approximately. Moreover, the global residuals of the discrete mass conservation law all reach the magnitude of 10−15 and the global residuals of the discrete energy conservation law all reach the magnitude of 10−14 , respectively. We observe a good agreement with the theoretical result. Example 5.2. We investigate the convergence order of the proposed compact scheme (3.1) in this experiment. Define

e (h, τ ) = u (·, T ) − u T (·)∞ ,

(5.2)

2 where T = 1 and u |t =0 = e −0.5x . Although we do not know the explicit form of the solution to (5.1), we take the compact

scheme with the very small time step-size

τ = 0.0001 as a reference solution. We then compare it to the compact scheme

78

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

Fig. 5.1. Evolution of the discrete mass conservation law (up), and the global errors of the discrete mass conservation law (bottom).

Fig. 5.2. Evolution of the discrete energy conservation law (up), and the global errors of the discrete energy conservation law (bottom).

in order to estimate the rate of convergence. To calculate the convergence rates both in time and space of the compact scheme, we use the formula below

Orderτ =

log(e (h, 2τ )) − log(e (h, τ )) log(2τ ) − log(τ )

, Orderh =

log(e (2h, τ )) − log(e (h, τ )) log(2h) − log(h)

.

The L ∞ errors and the numerical order of accuracy in temporal and spatial directions are contained in Table 5.1 and Table 5.2, respectively. We can see that the compact scheme gives a uniform fourth-order and second-order convergence in space and in time, respectively. Example 5.3. In this example, we show the soliton propagation of the NLS equation

iut + u xx + 2|u |2 u = 0,

(5.3)

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

79

Table 5.1 Accuracy test for NLS equation (2.1) in temporal direction with h = 0.01.

τ

L ∞ error

Orderτ

0.04 0.02 0.01 0.005

1.0285e−003 2.6897e−004 6.8111e−005 1.7066e−005

– 1.93 1.98 1.99

Table 5.2 Accuracy test for NLS equation (2.1) in spatial direction with τ = 0.0001. h

L ∞ error

Orderh

0.8 0 .4 0 .2 0.1 0.05

1.4807e−002 9.7508e−004 5.1477e−005 3.1399e−006 1.9492e−007

– 3.92 4.24 4.04 4.01

Fig. 5.3. The soliton propagation of equation (5.3) in [−30, 30] until T = 5.

with the soliton solution

 

u (x, t ) = sech(x − 4t ) exp 2i x −

3  t . 2

(5.4)

Waves of this form play an important role in complex physical situations. From Fig. 5.3, it can be shown that the envelope of the modulus |u (x, t )| may be considered as a soliton. Example 5.4. In this example, we show the double soliton collision of equation (5.3) with the initial condition

u (x, 0) =

2

j =1

exp

1 2



ic j (x − x j ) sech(x − x j ),

(5.5)

80

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

Fig. 5.4. The soliton propagation of equation (5.3) with initial condition (5.5) until T = 5.

with c 1 = 4, c 2 = −4, x1 = −10, x2 = 10. The solution is computed with homogeneous boundary condition in [−30, 30] and is shown in Fig. 5.4. 6. Conclusion In this paper, we developed a novel compact scheme for NLS equation in one dimension. From the theoretical analysis and numerical experiments, our scheme can preserve the mass and energy conservation laws, and the phase flow is symplectic. By the energy method and cut-off technique, we prove that our compact scheme is of second order in time and fourth order in space. The numerical results agree excellently with the theoretical analysis. However, these results can not be extended directly to the high-dimensional cases. Building the optimal error estimates in maximum norm for high-dimensional cases is our future work. References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York–London, 1975. [2] G. Akrivis, V. Dougalis, O. Karakashian, On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schrödinger equation, Numer. Math. 59 (1991) 31–53. [3] X. Antoine, W. Bao, C. Besse, Computational methods for the dynamics of the nonlinear Schrödinger/Gross–Pitaevskii equations, Comput. Phys. Commun. 184 (2013) 2621–2633. [4] W. Bao, S. Jin, P.A. Markowich, On time-splitting spectral approximation for the Schrödinger equation in the semiclassical regime, J. Comput. Phys. 175 (2002) 487–524. [5] W. Bao, J. Shen, A fourth-order time-splitting Laguerre–Hermite pseudo-spectral method for Bose–Einstein condensates, SIAM J. Sci. Comput. 26 (2005) 2010–2028. [6] Q. Chang, E. Jia, W. Sun, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys. 148 (1999) 397–415. [7] J. Chen, M. Qin, Y. Tang, Symplectic and multi-symplectic methods for the nonlinear Schrödinger equation, Comput. Math. Appl. 43 (2002) 1095–1106. [8] M. Dehghan, A. Taleei, Numerical solution of nonlinear Schrödinger equation by using time-space pseudo-spectral method, Numer. Methods Partial Differ. Equ. 26 (2010) 979–992. [9] K. Feng, On difference schemes and symplectic geometry, in: K. Feng (Ed.), Proc. of the 1984 Beijing Symp. on Diff. Geometry and Diff. Equations, Science Press, Beijing, 1985. [10] Z. Ge, J.E. Marsden, Lie–Poisson integrators and Lie–Poisson Hamilton–Jacobi theory, Phys. Lett. A 133 (1988) 134–139. [11] B. Guo, The convergence of numerical method for nonlinear Schrödinger equation, J. Comput. Math. 4 (1986) 121–130. [12] J. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes problem, I. Regularity of solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311.

J. Hong et al. / Applied Numerical Mathematics 120 (2017) 68–81

81

[13] O. Karakashian, C. Makridakis, A space–time finite element method for the nonlinear Schrödinger equation: the discontinuous Galerkin method, Math. Comput. 67 (1998) 479–499. [14] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [15] B. Li, G. Fairweather, B. Bialecki, Discrete-time orthogonal spline collocation methods for Schrödinger equations in two space variables, SIAM J. Numer. Anal. 35 (1998) 453–477. [16] P.A. Markowich, P. Pietra, C. Pohl, Numerical approximation of quadratic observables of Schrödinger-type equations in the semi-classical limit, Numer. Math. 81 (1999) 595–630. [17] C. Neuhauser, M. Thalhammer, On the convergence of splitting methods for linear evolutionary Schrödinger equations involving an unbounded potential, BIT Numer. Math. 49 (2009) 199–215. [18] D. Pathria, J.L. Morris, Pseudo-spectral solution of nonlinear Schrödinger equations, J. Comput. Phys. 87 (1990) 108–125. [19] M.P. Robinson, G. Fairweather, Orthogonal spline collocation methods for Schrödinger-type equations in one space variable, Numer. Math. 68 (1994) 355–376. [20] C. Sulem, P.L. Sulem, The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse, Springer-Verlag, New York, 1999. [21] Y. Tang, L. Vázquez, F. Zhang, V.M. Pérez-García, Symplectic methods for the nonlinear Schrödinger equation, Comput. Math. Appl. 32 (1996) 73–83. [22] M. Thalhammer, High-order exponential operator splitting methods for time dependent Schrödinger equations, SIAM J. Numer. Anal. 46 (2008) 2022–2038. [23] M. Thalhammer, M. Caliari, C. Neuhauser, High-order time-splitting Hermite and Fourier spectral methods, J. Comput. Phys. 228 (2009) 822–832. [24] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer-Verlag, Berlin, 1997. [25] T. Wang, Optimal point-wise error estimate of a compact difference scheme for the coupled Gross–Pitaevskii equations in one dimension, J. Sci. Comput. 59 (2014) 158–186. [26] T. Wang, B. Guo, Q. Xu, Fourth-order compact and energy conservative difference schemes for the nonlinear Schrödinger equation in two dimensions, J. Comput. Phys. 243 (2013) 382–399. [27] Y. Xu, C.-W. Shu, Local discontinuous Galerkin methods for nonlinear Schrödinger equations, J. Comput. Phys. 205 (2005) 72–77.