Convergence of an Eighth-Order Compact Difference Scheme for the Nonlinear Schrödinger Equation

A new compact difference scheme is proposed for solving the nonlinear Schrödinger equation. The scheme is proved to conserve the total mass and the total energy and the optimal convergent rate, without any restriction on the grid ratio, at the order of O h8 τ2 in the discrete L∞-norm with time step τ and mesh size h. In numerical analysis, beside the standard techniques of the energy method, a new technique named “regression of compactness” and some lemmas are proposed to prove the high-order convergence. For computing the nonlinear algebraical systems generated by the nonlinear compact scheme, an efficient iterative algorithm is constructed. Numerical examples are given to support the theoretical analysis.


Introduction
The time-dependent Schr ödinger equation is one of the most important equations in quantum mechanics 1 such as the Bose-Einstein condensate BEC , which is also used widely in many other different fields 2, 3 such as plasma physics, nonlinear optics, water waves, and bimolecular dynamics.There are many studies on numerical approaches, including finite difference 4-16 , finite element 17-19 , and polynomial approximation methods 20-24 , of the initial or initial-boundary value problems of the Schr ödinger equations.
Recently, there has been growing interest in high-order compact methods for solving partial differential equations PDEs 25-35 .It was shown that the high-order difference methods play an important role in the simulation of high-frequency wave phenomena.However, there is few proof of unconditional H 1 norm convergence of any compact difference scheme for nonlinear equations.

Advances in Numerical Analysis
In this paper, we consider the following cubic nonlinear Schr ödinger NLS equation: subject to periodic boundary condition u x, t u x L, t , x ∈ R, t > 0, 1.2 and initial condition where q is dimensionless constant characterizing the interaction positive for repulsive interaction and negative for the attractive interaction between particles in BEC.V x is a real function corresponding to the external trap potential and it is often chosen as a harmonic potential, that is, a quadratic polynomial, in most experiments.u : u x, t is the unknown periodical complex-valued function whose initial value ϕ x is a given periodical complexvalued function.L > 0 is the period of the u x, t .Denote Ω 0, L ; it is easy to show that the problem 1.1 -1.3 conserves the total mass and the global energy x u q|u| 4 dx ≡ E ϕ , t ≥ 0. 1.5 Fei et al. pointed out in 15 that the nonconservative schemes may easily show nonlinear blowup, and they presented a new conservative linear difference scheme for nonlinear Schr ödinger equation.In 16 , Li and Vu-Quoc also said, ". . . in some areas, the ability to preserve some invariant properties of the original differential equation is a criterion to judge the success of a numerical simulation." In this paper, a high-order compact HOC difference scheme is proposed for solving the problem 1.1 -1.3 .On the proposed HOC scheme, we obtain the following novel results.
i The proposed scheme is eighth order in space and second order in time.
ii The new scheme conserves the discrete total mass and the discrete global energy corresponding to 1.4 and 1.5 , respectively.
iii A new technique named "regression of compactness" is introduced to analyze the convergence of the proposed.
iv An efficient iterative algorithm is constructed to compute the proposed scheme.

Advances in Numerical Analysis
The remainder of this paper is arranged as follows.In Section 2, some notations are introduced, and a HOC difference scheme is constructed for the NLS equation.In Section 3, some useful lemmas are introduced or proved.In Section 4, discrete conservation laws of the proposed scheme are studied, then the convergence rate is proved at, without any restriction on the mesh ratio, the order of O h 8 τ 2 in the discrete maximum norm.In Section 5, an efficient iterative algorithm is designed to solve the nonlinear scheme.Lastly, numerical experiments are presented in Section 6 and some remarks are given in the concluding section.

Notations and Compact Finite Difference Schemes
Before giving the compact difference scheme, some notations are firstly introduced.For a positive integer N, choose time-step τ T/N and denote time steps t n nτ, where 0 < T < T max with T max the maximal existing time of the solution; choose mesh size h L/J, with J a positive integer and denote grid points as x j jh, j 0, 1, . . ., J.

2.1
Denote the index sets Denote U n j to be the numerical approximation and, respectively, u n j be the exact value of u x, t at the point x j , t n for j ∈ T e J and n 0, 1, 2, . . ., N, and denote U n ∈ C J 3 to be the numerical solution and, respectively, u n ∈ C J 3 be the exact solution at time t t n .For a grid function u {u n j | j ∈ T J , n 0, 1, 2, . . ., N}, introduce the following finite difference operators:

2.3
We denote the space and define seminorms and discrete inner product over X J as

Advances in Numerical Analysis
For simplicity, we write v 2 as v .Throughout the paper, we adopt the standard Sobolev spaces and their corresponding norms let C denote a generic constant which is independent of mesh size h and time step τ, 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 Cv.For the exact solution of the initial-boundary value problem 1.1 -1.3 , we assume that

Some Properties of Cyclic Matrix
A matrix in the form of is called a cyclic matrix.Because the matrix A is determined by the entries in the first row, the matrix can be denoted as A C a 0 , a 1 , . . ., a n−1 .

2.8
Denote G C 0, 1, 0, . . ., 0 , which is called as a basic cyclic matrix.By a simple calculation, we can get that G n C 1, 0, . . ., 0 I, I is the unit matrix.

2.9
Then, any cyclic matrix A C a 0 , a 1 , . . ., a n−1 can be written as a polynomial of G; that is,

2.10
This implies that, for any two cyclic matrices A and B, AB BA.Now we list some useful lemmas as follows.
Lemma 2.1.If A and B are two cyclic matrices, AB is still a cyclic matrix.
Proof.For any two cyclic matrices by using the multiplication of polynomials and the properties 2.9 we obtain where

2.13
Lemma 2.2.If a cyclic matrix A is invertible, then its inverse matrix A −1 is still a cyclic matrix.

Lemma 2.3. For any a real cyclic matrix
where Proof.On the one hand, it follows from |εI − G| 0 that ε n 1 which indicates that the eigenvalues of G are 2.17

Advances in Numerical Analysis
On the other hand, just as 2.10 , the real cyclic matrix A can be written as where a 0 , a 1 , . . ., a n−1 are n real numbers.This gives that the eigenvalues of A are in the form of

Compact Approximation of Second Derivative
The key point of the HOC method is to discretize the derivative with the fewest nodes to get the expected accuracy.In order to achieve higher accuracy, we adopt implicit compact scheme which equals a combination of the nodal derivatives to a combination of the nodal values of the function.Next, by Taylor's formula, expanding all the nodal derivatives and nodal values of the function at the same node, we can determine the coefficients of the combinations under the constrains of accuracy order.The nodal derivatives are implicitly evaluated herein by solving linear algebraic equations in strip.The stencil and the weighted factors are restricted to be symmetric so that the resulting schemes are nondissipative.Suppose f x is a L-periodic function, then for the discretization of the second-order derivatives f xx x , we have the formula 31 : where f j f x j , f xx j f xx x j , a, b, α, and β are undetermined parameters which depend on the accuracy order constraints.After tedious calculation, we can get these constraints:

2.21
According to 2.20 , an eighth-order scheme is obtained with which is the unique solution of the constraint equation 2.21 .
Thus, the approximation 2.20 to f xx j can be written in the form of finitedimensional linear operator where the finite difference operator A h is defined as

2.24
The corresponding matrix form is where is a positive, symmetric, and cyclic pentadiagonal matrix.

Lemma 2.4 see 36 . If a matrix H has real entries and is symmetric and positive definite, then
a H is invertible and its inverse matrix A −1 also has real entries and is symmetric and positive definite; b H can be decomposed as where D also has real entries and is symmetric and positive definite.

High-Order Compact Scheme
Imposing the approximations given in Section 2.2 on the spatial direction and the Crank-Nicolson discretization on the temporal direction of the NLS equation 1.1 gives where V j V x j and ξ n j is the local truncation error.The corresponding matrix form of the 2.28 -2.30 is u n ∈ X J , 0 n N,

2.31
where On the local truncation error ξ n j , using Taylor's expansion, we obtain the follwing lemma.

2.32
Omitting the small term ξ n j in 2.28 yields the following difference scheme Advances in Numerical Analysis 9 that is,

2.34
Remark 2.6.Though 2.33 and 2.34 are two different forms of a same scheme, the latter is suitable for computing in implementation.

Some Useful Lemmas
Lemma 2.3 gives that the eigenvalues of the cyclic matrix A C 1, α, β, 0, . . ., 0, β, α are in the form of This gives It is easy to verify that 1 − 2α − 2β > 0, which indicates that the matrix A is positive definite.
For the real, positive definite, symmetric and cyclic matrix A C 1, α, β, 0, . . ., 0, β, α , it follows from Lemma 2.2 that there are J real numbers m 0 , m 1 , . . ., m J−1 such that M C m 0 , m 1 , . . ., m J−1 .Lemmas 2.3 and 2.4 give that M is real positive definite and symmetric.And then we can introduce the following discrete norm:

3.3
On the relation between ||| • ||| and • , we have the following lemma.Proof.It follows from 3.2 that the eigenvalues of M satisfy , j ∈ T J .

Advances in Numerical Analysis
This gives the spectral radius ρ M 1/ 1 − 2α − 2β , and consequently then 3.5 follows from 3.3 and 3.7 .
The following lemmas will be used frequently in this paper.
Lemma 3.2.For any a grid function u ∈ X J , there are Mu Mδ 3.9 Proof.For the kth entry of the vector Mu, there is where the periodicity of u is used.This gives For the kth element of the vector Mδ x u, there is Then 3.8 is gotten from 3.11 and 3.12 .Similarly, 3.9 can be proved.
Lemma 3.3.For any two grid functions u, v ∈ X J ; there are

3.13
Advances in Numerical Analysis 11 Proof.Direct calculation and noticing the periodicity of the two grid functions yield 3.13 .
Lemma 3.4.For the approximation u n , u n 1 ∈ X J , there are where " • " and " • " mean taking the imaginary part and the real part, respectively.
Proof.Lemmas 3.2, 3.3, and definition 3.3 together with the symmetry of M give Mδ x u n u n 1 , δ x u n u n 1 δ x u n u n 1 2 0.

3.18
Similarly, 3.15 can be proved.Also using Lemmas 3.2, and 3.3, definition 3.3 , and the symmetry of M, we obtain

Numerical Analysis
The corresponding matrix form of 2.33 is where

Solvability of the Nonlinear Scheme
For proving the solvability, we need the following Brouwer fixed point theorem.Then, there exists a z * ∈ H such that g z * 0 and z * α.
On the solvability of the scheme 2.34 , we have the following lemma.
Lemma 4.2 Solvability of difference scheme .Under the assumption (A), for any initial data u 0 ∈ X J , there exists a solution U n ∈ X J of 2.34 for n ≥ 1.
Proof.The HOC scheme 2.33 is equivalent to the difference scheme 2.34 whose matrix form is 4.1 -4.3 , so we here just prove the unique solvability of the difference scheme 4.1 -4.3 .In 4.1 , for given U n ∈ X J n ≥ 1 , we first prove the existence.For j ∈ T J , rewrite 4.1 as

4.6
Advances in Numerical Analysis 13 We define a mapping F: which is obviously continue.Computing the inner product of 4.7 with u * and then taking the real part yield where Lemma 3.4 was used.Let U * U n 1, then it follows from 4.8 that F U * , U * ≥ 0. This together with Lemma 4.1 indecaites that there exists a solution U n 1 ∈ X J satisfying F U n 1 0. Then U n 1 satisfies the scheme 2.34 .

Discrete Conservation Laws
Corresponding to 1.4 -1.5 , the proposed scheme also possesses two counterpart discrete conservation laws.
Lemma 4.3.The difference solution of scheme 2.34 conserves the discrete mass and discrete energy in the form of where

4.11
Proof.Computing the inner product of 4.1 with U n U n 1 then taking the imaginary part yields

4.12
This together with Lemma 3.4 gives 4.9 .
Computing the inner product of 4.1 with δ t U n then taking the real part yields

4.13
This together with

A Priori Estimate
In proof of the convergence, we need the following a priori estimate.
Lemma 4.4.The difference solution of the scheme 2.34 satisfies

4.15
This implies the stability of the scheme 2.34 .
Proof.It follows from 4.9 that Lemma 3.1 gives 4.17 Using the discrete Sobolev inequality and 4.16 , we obtain

4.18
Advances in Numerical Analysis 15 where c 0 is a positive constant independent of the grid parameters.Taking ε ac

4.20
This completes the proof.

Unconditional Convergence in Maximum Norm
Theorem 4.5.Suppose that u x, t ∈ C 10,4 R × 0, T and τ is sufficiently small, then the solution u n j of the proposed scheme 2.34 converges to the solution u x j , t n of the problem V j e n e n 1 q 4 g n 1 ξ n , n 0, 1, . . ., N − 1, 4.22 e n ∈ X J , n 0, 1, 2, . . ., N.

Advances in Numerical Analysis
Noting that Computing the inner product of 4.22 with e n e n 1 and then taking the imaginary part, we obtain e n e n 1 , δ t e n a 4 Mδ 2 x e n e n 1 , δ t e n q 4 g n 1 , δ t e n ξ n , δ t e n , n 1, 2, . . ., N − 1.

4.32
This together with Lemma 3.4 gives

4.33
It follows from 4.22 that

4.34
Thus, for the third term on the left side of 4.33 , we have

4.42
Summing up for the superscript n from 0 to m and then replacing m by n and noting |δ

4.46
This completes the proof.

An Iterative Algorithm
Obviously, the scheme 2.34 is nonlinear and implicit.In order to obtain the solution U n 1 j at the time level n 1, an efficient iterative algorithm needs to be proposed.In this section, we construct an iterative algorithm to compute the scheme 2.34 .
For a fixed n, the difference scheme 2.34 can be solved by the following iterative algorithm: .

Advances in Numerical Analysis
Denote then the iterative algorithm can be written as the following matrix form where H C c, b, a, 0, . . ., 0, a, b 5.5 is a cyclic pentadiagonal constant matrix.Though the algebraical system Hu n 1 s 1 d n 1 s generated by the proposed scheme is nonlinear and then need iteration in implementation, the coefficient matrix H is cyclic pentadiagonal and invariable at different levels and then the matrix H should not be computed repeatedly.This together with the proposed iterative algorithm and the Thomas algorithm can improve greatly the computational efficiency.To some extent, the proposed scheme is more efficient than most linearized scheme.

Numerical Experiments
In this section, numerical tests are given to support our theoretical analysis on convergence and stability.In implementation, we chose a tolerance ε 10 −10 to control the iteration; that is, U n 1 s 1 − U n 1 s ∞ 10 −10 .Example 6.1.Periodic initial-boundary value problem.

6.1
This problem has the following plane wave solution: Now, we compute Example 6.1 on the finite domain 0, π .In order to test the maximum norm convergence with order O h 8 τ 2 , we let τ O h 4 in computation to test the eighth-order convergence in space and take a small h in computation to test second-order convergence in time.In order to quantify the numerical results, we introduce the following "error" functions and "convergence rate:"  where Ω −40, 40 .
In computation of Example 6.2, we assume that the value of the exact solution is zero when x ∈ Ω c : R \ Ω and list the numerical results of Example 6.2 in Tables 3-4.
Tables 1-4 show that the new scheme is eighth-order convergent in space and second order in time for U in the discrete L ∞ -norm, which verifies Theorem 4.5.Tables 2 and 4 also show that the proposed scheme is very robust.

Conclusion
In this paper, we construct an eighth-order compact different scheme for numerically solving the NLS equation, use the discrete energy method together with a new technique named "regression of compactness" to prove that the proposed HOC scheme is convergent, without any restriction on the grid ratio, at the order of O h 8 τ 2 in the discrete L ∞ -norm.What's better, the HOC scheme conserves the discrete total mass and energy, which changed the viewpoint that the compact difference scheme of a nonlinear equation cannot preserve the discrete conservation laws.In order to implement the nonlinear compact difference scheme, we propose an iterative algorithm which is numerically proved to be very reliable.

Proof.
Denote the error function as e n u n − U n , n 0, 1, 2, . . ., N. 4.21 Subtracting 4.1 -4.3 from 2.31 , respectively, we obtain the equation of error function as follows: 4.25 , 2.6 , and Lemma 4.4 that Computing the inner product of 4.22 with δ t e n and then taking the real part, we obtain

Table 1 :
Errors of the numerical solution of Example 6.1 computed by the proposed scheme at t 1.

Table 2 :
Errors of the numerical solution of Example 6.1 computed by the proposed scheme with h π/40 at t 1.
, τThen we list the numerical results of Example 6.1 in Tables 1-2.

Table 3 :
Errors of the numerical solution of Example 6.2 computed by the proposed scheme at t 1.

Table 4 :
Errors of the numerical solution of Example 6.2 computed by the proposed scheme with h 0.15 at t 1.Because u x, t decays to zero rapidly as |x| → ∞, so for sufficiently large number −x L and x R and a finite time T , the initial problem is consistent with the following initial-boundary value problem: 2 x, t 2|u x, t | 2 u x, t 0, x, t ∈ R × 0, T , u x, 0 e 2ix , x ∈ R, u x, t −→ 0, as |x| −→ 0.