Advances in Numerical AnalysisVolume 2012 (2012), Article ID 913429, 24 pagesdoi:10.1155/2012/913429

Research Article

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

Tingchun Wang

School of Mathematics and Statistics, Nanjing University of Information Science and Technology, Nanjing 210044, China

Received 24 May 2012; Accepted 7 August 2012

Academic Editor: Ting-Zhu Huang

Copyright © 2012 Tingchun Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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 in the discrete -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.

1. 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 norm convergence of any compact difference scheme for nonlinear equations.

In this paper, we consider the following cubic nonlinear Schrödinger (NLS) equation: subject to periodic boundary condition and initial condition where is dimensionless constant characterizing the interaction (positive for repulsive interaction and negative for the attractive interaction) between particles in BEC. 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. is the unknown periodical complex-valued function whose initial value is a given periodical complex-valued function. is the period of the .

Denote ; it is easy to show that the problem (1.1)–(1.3) conserves the total mass and the global energy 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.

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 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.

2. Notations and Compact Finite Difference Schemes

Before giving the compact difference scheme, some notations are firstly introduced. For a positive integer , choose time-step and denote time steps , where with the maximal existing time of the solution; choose mesh size , with a positive integer and denote grid points as Denote the index sets Denote to be the numerical approximation and, respectively, be the exact value of at the point for and , and denote to be the numerical solution and, respectively, be the exact solution at time . For a grid function , introduce the following finite difference operators:

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

For simplicity, we write as . Throughout the paper, we adopt the standard Sobolev spaces and their corresponding norms let denote a generic constant which is independent of mesh size and time step , and use the notation to present that there is a generic constant which is independent of time step and mesh size such that . For the exact solution of the initial-boundary value problem (1.1)–(1.3), we assume that

2.1. Some Properties of Cyclic Matrix

A matrix in the form of is called a cyclic matrix. Because the matrix is determined by the entries in the first row, the matrix can be denoted as Denote , which is called as a basic cyclic matrix. By a simple calculation, we can get that Then, any cyclic matrix can be written as a polynomial of ; that is, This implies that, for any two cyclic matrices and , . 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 we obtain where

Lemma 2.2. If a cyclic matrix is invertible, then its inverse matrix is still a cyclic matrix.

Proof. If is invertible, we assume that the inverse matrix is , then Multiplying with (2.14) yields If , we obtain the linear system . Noting that , we obtain that the system has a unique nonzero solution . This together with (2.15) gives that where is the unique nonzero solution of . It follows that the assumption is right; that is, is a cyclic matrix.

Lemma 2.3. For any a real cyclic matrix , all the eigenvalues are in the form of where and .

Proof. On the one hand, it follows from that which indicates that the eigenvalues of are On the other hand, just as (2.10), the real cyclic matrix can be written as where are real numbers. This gives that the eigenvalues of are in the form of

2.2. 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 is a -periodic function, then for the discretization of the second-order derivatives , we have the formula [31]: where , , , , , and are undetermined parameters which depend on the accuracy order constraints. After tedious calculation, we can get these constraints: According to (2.20), an eighth-order scheme is obtained with which is the unique solution of the constraint equation .

Thus, the approximation (2.20) to can be written in the form of finite-dimensional linear operator where the finite difference operator is defined as 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) is invertible and its inverse matrix also has real entries and is symmetric and positive definite;(b) can be decomposed as where also has real entries and is symmetric and positive definite.

2.3. 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 and is the local truncation error. The corresponding matrix form of the (2.28)–(2.30) is where and .

On the local truncation error , using Taylor's expansion, we obtain the follwing lemma.

Lemma 2.5. Suppose that and , there is

Omitting the small term in (2.28) yields the following difference scheme that is,

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

3. Some Useful Lemmas

Lemma 2.3 gives that the eigenvalues of the cyclic matrix are in the form of This gives It is easy to verify that , which indicates that the matrix is positive definite.

For the real, positive definite, symmetric and cyclic matrix , it follows from Lemma 2.2 that there are real numbers such that . Lemmas 2.3 and 2.4 give that is real positive definite and symmetric. And then we can introduce the following discrete norm: On the relation between and , we have the following lemma.

Lemma 3.1. The discrete norms and are equivalent; in fact, there is where , .

Proof. It follows from (3.2) that the eigenvalues of satisfy This gives the spectral radius ), and consequently Because 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 , there are

Proof. For the entry of the vector , there is where the periodicity of is used. This gives For the element of the vector , 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 ; there are

Proof. Direct calculation and noticing the periodicity of the two grid functions yield .

Lemma 3.4. For the approximation , 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 give Similarly, (3.15) can be proved. Also using Lemmas 3.2, and 3.3, definition (3.3), and the symmetry of , we obtain Similarly, (3.17) can be proved.

Lemma 3.5 (see [33]). For time sequences and , there is

4. Numerical Analysis

The corresponding matrix form of is where

4.1. Solvability of the Nonlinear Scheme

For proving the solvability, we need the following Brouwer fixed point theorem.

Lemma 4.1 (see [37]). Let be a finite-dimensional inner-product space, the associated norm, and continuous. Assume, moreover, that Then, there exists a such that and .

On the solvability of the scheme , we have the following lemma.

Lemma 4.2 (Solvability of difference scheme). Under the assumption , for any initial data , there exists a solution of for .

Proof. The HOC scheme is equivalent to the difference scheme 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 , we first prove the existence. For , rewrite (4.1) as We define a mapping : as which is obviously continue. Computing the inner product of (4.7) with and then taking the real part yield where Lemma 3.4 was used. Let , then it follows from (4.8) that . This together with Lemma 4.1 indecaites that there exists a solution satisfying . Then satisfies the scheme .

4.2. 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 conserves the discrete mass and discrete energy in the form of where

Proof. Computing the inner product of (4.1) with then taking the imaginary part yields This together with Lemma 3.4 gives (4.9).Computing the inner product of (4.1) with then taking the real part yields This together with and Lemma 3.4 gives (4.10).

4.3. A Priori Estimate

In proof of the convergence, we need the following a priori estimate.

Lemma 4.4. The difference solution of the scheme satisfies This implies the stability of the scheme .

Proof. It follows from (4.9) that Lemma 3.1 gives Using the discrete Sobolev inequality and (4.16), we obtain where is a positive constant independent of the grid parameters. Taking , then it follows from (4.16)–(4.18) and (4.10) that Utilizing the discrete Sobolev estimate and Lemma 4.3, we obtain from (4.16) and (4.19) that This completes the proof.

4.4. Unconditional Convergence in Maximum Norm

Theorem 4.5. Suppose that and is sufficiently small, then the solution of the proposed scheme converges to the solution of the problem (1.1)–(1.3) with order in the discrete -norm.

Proof. Denote the error function as Subtracting (4.1)–(4.3) from , respectively, we obtain the equation of error function as follows: where Noting that we obtain from (4.25), (2.6), and Lemma 4.4 that Computing the inner product of (4.22) with and then taking the imaginary part, we obtain where Lemma 3.4 was used.For the last two terms of (4.28), using Lemma 4.4 and (4.26) we have Substituting into (4.28) yields By virtue of Gronwall's inequality and Lemma 2.5, noting that , we obtain from (4.30) that if is sufficiently small.Computing the inner product of (4.22) with and then taking the real part, we obtain This together with Lemma 3.4 gives It follows from (4.22) that Thus, for the third term on the left side of (4.33), we have Noting (4.26), (4.31), and we obtain Similarly we obtain For the last two terms of (4.35), using the Cauchy inequality, noting that , we obtain Substituting (4.38)–(4.40) into (4.35) yields Substituting (4.41) into (4.33) yields Summing up for the superscript from to and then replacing by and noting , we obtain By virtue of Lemmas 3.5, 2.5, and (4.31), we obtain Substituting (4.44) into (4.43) and then utilizing Gronwall's inequality, we obtain if is sufficiently small. This together with (4.31) gives This completes the proof.

5. An Iterative Algorithm

Obviously, the scheme is nonlinear and implicit. In order to obtain the solution at the time level , an efficient iterative algorithm needs to be proposed. In this section, we construct an iterative algorithm to compute the scheme .

For a fixed , the difference scheme can be solved by the following iterative algorithm: with Denote then the iterative algorithm can be written as the following matrix form where is a cyclic pentadiagonal constant matrix.

Though the algebraical system generated by the proposed scheme is nonlinear and then need iteration in implementation, the coefficient matrix is cyclic pentadiagonal and invariable at different levels and then the matrix 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.

6. Numerical Experiments

In this section, numerical tests are given to support our theoretical analysis on convergence and stability. In implementation, we chose a tolerance to control the iteration; that is, .

Example 6.1. Periodic initial-boundary value problem. Consider This problem has the following plane wave solution:

Now, we compute Example 6.1 on the finite domain . In order to test the maximum norm convergence with order , we let in computation to test the eighth-order convergence in space and take a small 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:” Then we list the numerical results of Example 6.1 in Tables 1-2.

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

Error
Sorder

0.2

—

0.2/16

8.1724

8.0007

7.9989

Table 2: Errors of the numerical solution of Example 6.1 computed by the proposed scheme with at .

Error
Torder

0.2

—

0.1

2.1600

0.05

2.0140

0.025

1.9995

Example 6.2. Homogenous initial-boundary value problem.Consider The initial value problem has the one-soliton solution Because decays to zero rapidly as , so for sufficiently large number and and a finite time , the initial problem is consistent with the following initial-boundary value problem: where .

In computation of Example 6.2, we assume that the value of the exact solution is zero when and list the numerical results of Example 6.2 in Tables 3-4.

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

Error
Sorder

0.5
0.2

—

0.5/2
0.2/16

7.6658

8.0043

7.9990

Table 4: Errors of the numerical solution of Example 6.2 computed by the proposed scheme with at .

Error
Torder

0.1

—

0.05

2.1461

0.025

2.0338

0.0125

2.0064

Tables 1–4 show that the new scheme is eighth-order convergent in space and second order in time for in the discrete -norm, which verifies Theorem 4.5. Tables 2 and 4 also show that the proposed scheme is very robust.

7. 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 in the discrete -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.

Acknowledgment

This work is supported by the National Natural Science Foundation of China, nos. 11126292 and 11201239.

References

D. J. Griffiths, Introduction to Quantum Mechanics, Prentice-Hall, Englewood Cliffs, NJ, USA, 1995.

C. R. Menyuk, “Stability of solitons in birefringent optical fibers,” Journal of the Optical Society of America B, vol. 5, pp. 392–402, 1998.

M. Wadati, T. Iizuka, and M. Hisakado, “A coupled nonlinear Schrödinger equation and optical solitons,” Journal of the Physical Society of Japan, vol. 61, no. 7, pp. 2241–2245, 1992.

G. D. Akrivis, “Finite difference discretization of the cubic Schrödinger equation,” IMA Journal of Numerical Analysis, vol. 13, no. 1, pp. 115–124, 1993.

T. F. Chan and L. J. Shen, “Stability analysis of difference schemes for variable coefficient Schrödinger type equations,” SIAM Journal on Numerical Analysis, vol. 24, no. 2, pp. 336–349, 1987.

Q. Chang, E. Jia, and W. Sun, “Difference schemes for solving the generalized nonlinear Schrödinger equation,” Journal of Computational Physics, vol. 148, no. 2, pp. 397–415, 1999.

W. Z. Dai, “An unconditionally stable three-level explicit difference scheme for the Schrödinger equation with a variable coefficient,” SIAM Journal on Numerical Analysis, vol. 29, no. 1, pp. 174–181, 1992.

M. Dehghan and A. Taleei, “A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients,” Computer Physics Communications, vol. 181, no. 1, pp. 43–51, 2010.

F. Ivanauskas and M. Radžiūnas, “On convergence and stability of the explicit difference method for solution of nonlinear Schrödinger equations,” SIAM Journal on Numerical Analysis, vol. 36, no. 5, pp. 1466–1481, 1999.

P. L. Nash and L. Y. Chen, “Efficient finite difference solutions to the time-dependent Schrödinger equation,” Journal of Computational Physics, vol. 130, no. 2, pp. 266–268, 1997.

J. M. Sanz-Serna, “Methods for the numerical solution of the nonlinear Schroedinger equation,” Mathematics of Computation, vol. 43, no. 167, pp. 21–27, 1984.

Z. Z. Sun and X. Wu, “The stability and convergence of a difference scheme for the Schrödinger equation on an infinite domain by using artificial boundary conditions,” Journal of Computational Physics, vol. 214, no. 1, pp. 209–223, 2006.

L. Wu, “Dufort-Frankel-type methods for linear and nonlinear Schrödinger equations,” SIAM Journal on Numerical Analysis, vol. 33, no. 4, pp. 1526–1533, 1996.

J. M. Sanz-Serna and J. G. Verwer, “Conservative and nonconservative schemes for the solution of the nonlinear Schrödinger equation,” IMA Journal of Numerical Analysis, vol. 6, no. 1, pp. 25–42, 1986.

Z. Fei, V. M. Pérez-García, and L. Vázquez, “Numerical simulation of nonlinear Schrödinger systems: a new conservative scheme,” Applied Mathematics and Computation, vol. 71, no. 2-3, pp. 165–177, 1995.

S. Li and L. Vu-Quoc, “Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein-Gordon equation,” SIAM Journal on Numerical Analysis, vol. 32, no. 6, pp. 1839–1875, 1995.

G. D. Akrivis, V. A. Dougalis, and O. A. Karakashian, “On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schrödinger equation,” Numerische Mathematik, vol. 59, no. 1, pp. 31–53, 1991.

O. Karakashian, G. D. Akrivis, and V. A. Dougalis, “On optimal order error estimates for the nonlinear Schrödinger equation,” SIAM Journal on Numerical Analysis, vol. 30, no. 2, pp. 377–400, 1993.

Y. Tourigny, “Some pointwise estimates for the finite element solution of a radial nonlinear Schrödinger equation on a class of nonuniform grids,” Numerical Methods for Partial Differential Equations, vol. 10, no. 6, pp. 757–769, 1994.

W. Bao and D. Jaksch, “An explicit unconditionally stable numerical method for solving damped nonlinear Schrödinger equations with a focusing nonlinearity,” SIAM Journal on Numerical Analysis, vol. 41, no. 4, pp. 1406–1426, 2003.

T. Iitaka, “Solving the time-dependent Schröinger equation numerically,” Physical Review E, vol. 49, no. 5, pp. 4684–4690, 1994.

D. Kosloff and R. Kosloff, “A Fourier method solution for the time dependent Schröinger equation as a tool in molecular dynamics,” Journal of Computational Physics, vol. 52, pp. 35–53, 1983.

B. Li, G. Fairweather, and B. Bialecki, “Discrete-time orthogonal spline collocation methods for Schrödinger equations in two space variables,” SIAM Journal on Numerical Analysis, vol. 35, no. 2, pp. 453–477, 1998.

M. P. Robinson and G. Fairweather, “Orthogonal spline collocation methods for Schrödinger-type equations in one space variable,” Numerische Mathematik, vol. 68, no. 3, pp. 355–376, 1994.

G. Berikelashvili, M. M. Gupta, and M. Mirianashvili, “Convergence of fourth order compact difference schemes for three-dimensional convection-diffusion equations,” SIAM Journal on Numerical Analysis, vol. 45, no. 1, pp. 443–455, 2007.

G. Cohen, High-Order Numerical Methods for Transient Wave Equations, Springer, New York, NY, USA, 2002.

A. Gopaul and M. Bhuruth, “Analysis of a fourth-order scheme for a three-dimensional convection-diffusion model problem,” SIAM Journal on Scientific Computing, vol. 28, no. 6, pp. 2075–2094, 2006.

B. Gustafsson and E. Mossberg, “Time compact high order difference methods for wave propagation,” SIAM Journal on Scientific Computing, vol. 26, no. 1, pp. 259–271, 2004.

B. Gustafsson and P. Wahlund, “Time compact difference methods for wave propagation in discontinuous media,” SIAM Journal on Scientific Computing, vol. 26, no. 1, pp. 272–293, 2004.

J. C. Kalita, D. C. Dalal, and A. K. Dass, “A class of higher order compact schemes for the unsteady two-dimensional convection-diffusion equation with variable convection coefficients,” International Journal for Numerical Methods in Fluids, vol. 38, no. 12, pp. 1111–1131, 2002.

S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992.

H. L. Liao and Z. Z. Sun, “Maximum norm error bounds of ADI and compact ADI methods for solving parabolic equations,” Numerical Methods for Partial Differential Equations, vol. 26, no. 1, pp. 37–60, 2010.

H. L. Liao, Z. Z. Sun, and H. S. Shi, “Error estimate of fourth-order compact scheme for linear Schrödinger equations,” SIAM Journal on Numerical Analysis, vol. 47, no. 6, pp. 4381–4401, 2010.

S. Xie, G. Li, and S. Yi, “Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrödinger equation,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 9-12, pp. 1052–1060, 2009.

D. W. Zingg, “Comparison of high-accuracy finite-difference methods for linear wave propagation,” SIAM Journal on Scientific Computing, vol. 22, no. 2, pp. 476–502, 2000.

R. A. Horn and C. R. Johnson, Matrix Analysis, chapter 7, Cambridge University Press, 1985.

F. E. Browder, “Existence and uniqueness theorems for solutions of nonlinear boundary value problems,” in Application of Nonlinear Partial Differential Equations, R. Finn, Ed., Proceedings of Symposia in Applied Mathematics, pp. 24–49, American Mathematical Society, Providence, RI, USA, 1965.