An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions Xin Li a,b ,1 , Luming Zhang a , a b

∗,2

Department of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, Jiangsu, PR China Department of Mathematics, Anhui Science and Technology University, Fengyang, 233100, Anhui, PR China

article

info

Article history: Received 8 March 2019 Received in revised form 21 October 2019 Accepted 2 November 2019 Available online xxxx Keywords: Schrödinger equation Neumann boundary conditions Spectral-collocation Conservation laws Discrete cosine transform Fast algorithm

a b s t r a c t Based on the second-order cosine spectral differentiation matrix deduced in this paper, we present an efficient and high accurate numerical scheme for solving two-dimensional Schrödinger equation with Neumann boundary conditions. The Crank–Nicolson finite difference scheme is utilized in temporal discretization and the cosine spectral-collocation method is employed for the approximation of Laplacian operator. The new scheme conserves the mass and energy in discrete level and is implemented by a fast algorithm in terms of the relations between the cosine spectral differentiation matrix and fast cosine transform. More importantly, this strategy can be extended to address the problems with high (even)-order derivatives in space. Numerical examples and applications are listed to confirm the validity and high accuracy of the method. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction A range of physical phenomena can be described by the nonlinear Schrödinger (NLS) equation. It has been widely used in many fields, such as plasma physics, quantum mechanics, nonlinear optics, dynamics of Bose–Einstein condensates at the extremely low temperature and so on [1–6]. In this paper, we consider the following cubic NLS equation in two dimensions: iut + α ∆u + V (x)u + β|u|2 u = 0

for x ∈ Ω and 0 < t ≤ T ,

(1.1)

with Neumann boundary conditions

∂ u(x, t) = 0 for x ∈ ∂ Ω and 0 < t ≤ T , ∂n

(1.2)

and initial condition u(x, 0) = u0 (x)

¯, for x ∈ Ω

(1.3)

∗ Corresponding author. E-mail addresses: [email protected] (X. Li), [email protected] (L. Zhang). 1 He is supported by a grant KJ2018A0523 from the University Natural Science Research Key Project of Anhui Province, PR China. 2 He is supported by a grant 11571181 from the National Nature Science Foundation of China. https://doi.org/10.1016/j.camwa.2019.11.006 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

2

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

where Ω = [a, b] × [c , d] is a bounded domain, ∂ Ω is the boundary, ∆ = ∂xx + ∂yy is the two-dimensional Laplacian operator, V is a real-valued potential function, α, β are two nonzero real constants. The NLS equation (1.1) is focusing for β > 0 and defocusing for β < 0, which is a generic model for the slowly varying envelop of a wave-train in conservative, dispersive, mildly nonlinear wave phenomena. It is also obtained as the subsonic limit of the Zakharov model for Langmuir waves in plasma physics [7]. Tracing back to the literatures, extensive analytical and numerical studies have been carried out for the NLS equation. Along the mathematical front, readers can refer to [8–11] and the references therein for the derivation, well-posedness and dynamical properties of this equation. Actually, it is not difficult to check that the analytical solution of problem (1.1)–(1.3) preserves the following total mass

∫ Q (t) := Ω

|u(x, t)|2 dx = Q (0),

(1.4)

and the energy

∫ E(t) := Ω

] [ β α|∇ u(x, t)|2 − |u(x, t)|4 − v (x)|u(x, t)|2 dx = E(0). 2

(1.5)

Along the numerical front, various efficient and accurate numerical methods have been developed for the NLS equation, including the finite element method [12,13], discontinuous Galerkin method [14–16], meshless method [17,18], finite difference method [19–24], spectral/pseudo-spectral method [25–30] and so on. For comparisons between different numerical methods, we refer to [31–36] and references therein for more details. Among these approximations, the finite difference schemes are popular and widely utilized for its simplicity and flexibility in computation. Previous studies mainly focus on the Dirichlet boundary conditions and periodic boundary conditions. For the problems with Neumann boundary conditions, to the best of the authors’ knowledge, literatures are limited. Liao et al. [37] use the finite difference method and present a fourth-order compact algorithm for nonlinear reaction-diffusion equations with flux boundary conditions. The authors of [38] show a new idea and construct a high-order difference scheme at both interior and boundary points for the heat equation with Neumann boundary conditions. Then, in [39], it is generalized to the multi-dimensional heat problems with respect to the same boundary conditions. Sun [40] combines these schemes together and establishes rigorous error estimates. Theoretical results show that the convergence order of O(τ 2 + h4 ) is tough to be obtained. In general, Neumann boundary conditions are much more difficult to handle compared with the other two boundary conditions. Spectral/pseudo-spectral method is a classical technique and widely adopted to solve differential equations in recently years for its superior properties, especially the high-order accuracy. In [41,42], pseudo-spectral Chebyshev method and spectral Galerkin method are employed to investigate parabolic equations and elliptic equations with Neumann boundary conditions, respectively. Moreover, Jacobi spectral method and Chebyshev-Fourier collocation method are also applied to the problems with flux boundary conditions [43,44]. Overall, spectral/pseudo-spectral method is superior to the finite difference method in the performance of computational accuracy. Very recently, Gong and his collaborators [45] study the relations between Fourier spectral differentiation matrix and discrete Fourier transform, then successfully solve the NLS equation with periodic boundary conditions by a conservative Fourier pseudo-spectral scheme in [30]. Enlightened by this idea, we strive to extend the results to Neumann boundary conditions, and present a conservative spectral-collocation difference scheme with high accuracy for two-dimensional Schrödinger equation. Our main contribution in this paper is to deduce the second-order cosine spectral differentiation matrix (CSDM) and investigate its relationships with the high (even)-order CSDM, see Lemmas 2.2 and 2.3. It provides us a significant strategy to treat the discretizations such as the fourth-order derivative in space, and unfolds a new thinking to address the reaction–diffusion problems with the same boundary conditions. Furthermore, we also reveal the relations between the CSDM and discrete cosine transform in Lemma 2.4. It would be useful for developing the accelerated algorithm in computation. Subsequently, a numerical scheme which is relied on the Crank–Nicolson finite difference scheme in time and the cosine spectral-collocation method in space is developed. Always, we illustrate this pseudo-spectral scheme still maintains the discrete conservation laws of mass and energy as same as the finite difference schemes, see Theorem 3.1. By means of the relations described earlier, an accelerated algorithm is constructed to speed up the numerical implementation thanks to the fast cosine transform. The paper is arranged as follows. In Section 2, the second-order CSDM is derived and some essential properties are studied. We present the numerical scheme in Section 3 and then demonstrate the discrete conservation laws. In the next section, fast algorithm is discussed in detail and some examples are reported to witness the effectiveness and accuracy of the new method. Applications about simulating the soliton interactions and the dynamics of one-component Bose–Einstein condensates are displayed in Section 5. Lastly, the conclusion will be made. 2. Spectral-collocation method 2.1. Cosine spectral differentiation matrix Let Lx = b − a and Ly = d − c. For two positive integers Nx and Ny , we denote two spatial mesh sizes hx := Lx /Nx and hy := Ly /Ny in the x- and y-directions, respectively, and the maximum length h := max{hx , hy }. The time-step size is Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

3

denoted as τ := T /Nτ with Nτ is also a positive integer. We denote Ωτ = {tn |tn = nτ , n = 0, 1, . . . , Nτ } and the discrete grid points Ωh = {(xj , yk )| j = 0, 1, . . . , Nx , k = 0, 1, . . . , Ny }, where xj = a + jhx and yk = c + khy . For simplicity of notations, an index set Th = {(j, k) | j = 0, 1, . . . , Nx − 1, k = 0, 1, . . . , Ny − 1} is also introduced. π π , µy := d− and define the following To deal with the Neumann boundary conditions properly, we set µx := b− a c solution space SN := span ˜ Xl (x)˜ Ym (y), (l, m) ∈ Th ,

{

}

where ˜ Xl (x) and ˜ Ym (y) are basis functions and can be denoted by

( ) ( ) ˜ Xl (x) := cos µx l(x − a) , ˜ Ym (y) := cos µy m(y − c) .

(2.1)

Always, the collocation points are chosen in (xp+ 1 , yq+ 1 ) for (p, q) ∈ Th . Let IN be the interpolation operator onto SN , i.e. 2

2

Nx −1 Ny −1

∑∑

(IN u)(x, y) :=

˜ ulm˜ Xl (x)˜ Ym (y),

(2.2)

l=0 m=0

where the corresponding expansion coefficient reads

˜ ulm :=

Nx −1 Ny −1

4

∑∑

Nx Ny Zl Zm

u(xj+ 1 , yk+ 1 )Xl (xj+ 1 )Ym (yk+ 1 ), 2

j=0

2

2

(2.3)

2

k=0

with Zξ := 1 + δ0ξ (δ0ξ is the Kronecker delta symbol) for any integer ξ ∈ N. Recalling the basis functions in (2.1), we plug (2.3) into (2.2) and obtain the following interpolation formula Nx −1 Ny −1

∑∑

(IN u)(x, y) =

u(xj+ 1 , yk+ 1 )Xj+ 1 (x)Yk+ 1 (y), 2

j=0

2

2

(2.4)

2

k=0

with the interpolation basis functions Xj+ 1 (x) := 2

Yk+ 1 (y) := 2

Nx −1 2 ∑ 1

Nx

Zl

l=0

cos µx l(xj+ 1 − a) cos µx l(x − a) ,

m=0

(

)

(2.5)

2

Ny −1 2 ∑ 1

Ny

)

(

cos µy m(yk+ 1 − c) cos µy m(y − c) .

)

(

Zm

(

)

(2.6)

2

Based on the results mentioned above, the following lemma can be stated. Lemma 2.1. Denote θx := µx (x − xj+ 1 ), θy := µy (y − yk+ 1 ), θx′ := µx (x + xj+ 1 − 2a), θy′ := µy (y + yk+ 1 − 2c), the 2 2 2 2 interpolation basis functions can be reformulated as Xj+ 1 (x) = C(Nx , µx , θx ) + C(Nx , µx , θx′ ),

(2.7)

Yk+ 1 (y) = C(Ny , µy , θy ) + C(Ny , µy , θy ),

(2.8)

2

′

2

with C(N , µ, θ ) :=

1 2N

[sin(N θ ) cot θ2 − cos(N θ )], and they satisfy

Xj+ 1 (xp+ 1 ) = δjp , 2

Yk+ 1 (yq+ 1 ) = δkq .

2

2

(2.9)

2

Proof. For simplicity, we just demonstrate (2.7) and the first assertion of (2.9). (2.8) and the second part of (2.9) can be obtained identically. We find in (2.5) that Xj+ 1 (x) = 2

Nx −1 2 ∑

Nx

cos µx l(xj+ 1 − a) cos µx l(x − a) +

(

)

(

)

2

l=1

1 Nx

.

Applying the prosthaphaeresis formula and the definitions of θx and θx′ , it is easy to check Xj+ 1 (x) = 2

1 2Nx

Nx −1

∑[

1 + 2 cos(θx l) + 1 + 2 cos(θx′ l)

(

)]

l=1

( 1 [ ( 1 ) θx 1 ) θ′ ] = sin (Nx − )θx csc + sin (Nx − )θx′ csc x . 2Nx 2 2 2 2

(2.10)

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

4

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

Resorting to the following fact 1

(

sin (N −

2

)θ csc

)

θ 2

sin(N θ ) cos

=

θ

− cos(N θ ) sin θ2

2

sin

= sin(N θ ) cot

θ

2

θ 2

− cos(N θ ),

we obtain (2.7) directly. In view of (2.10), it is clear that Xj+ 1 (x) are the nodal basis functions which satisfy the first 2 assertion of (2.9), i.e. Xj+ 1 (xp+ 1 ) = 1 for j = p and Xj+ 1 (xp+ 1 ) = 0 for j ̸ = p. This completes the proof of Lemma 2.1. □ 2

2

2

2

In order to obtain the spatial discretization of Laplacian operator, we differentiate (2.5)–(2.6) and evaluate the expressions at the collocation points (xp+ 1 , yq+ 1 ) for (p, q) ∈ Th . Here and in what follows, with a slight abuse of notation, 2 2 we denote Ap,q as the element of the matrix A in the q + 1 column of row p + 1 for (p, q) ∈ Th . Then, the second-order CSDM is given below. Lemma 2.2. Under the definitions of Lemma 2.1, the second-order derivatives of interpolation basis functions are given as X′′

j+ 21

¨ (Nx , µx , θx ) + C¨ (Nx , µx , θx′ ), Y′′ (x) = C

k+ 21

¨ (Ny , µy , θy ) + C¨ (Ny , µy , θy′ ), (y) = C

where

C¨ (N , µ, θ ) = −

N µ2 2

θ

sin(N θ ) cot

µ2

−

2

2

θ

cos(N θ ) csc2

2

µ2

+

4N

sin(N θ )

cos

θ

2 sin3 θ2

N µ2

+

2

cos(N θ ).

Then, the elements of the second-order CSDM satisfy

⎧ (µ h ) (µ h )] µ2 [ ⎪ ⎪ ⎨ (−1)j+p x csc2 x x (j + p + 1) − csc2 x x (j − p) , 2

(C2x )j,p =

2

j ̸ = p,

2

) µ ( ⎪ 1 µ µ ⎪ ⎩ csc2 (j + )µx hx − − , 2 2 6 3 ⎧ (µ h ) (µ h )] µ2 [ ⎪ ⎪ ⎨ (−1)k+q y csc2 x x (k + q + 1) − csc2 x x (k − q) , 2 x

2 x

2

y

(C2 )k,q =

2 x

2

µ2y

⎪ ⎪ ⎩

Nx2

2

j = p, k ̸ = q,

2

csc2 (k +

(

1 2

)µy hy −

)

µ2y 6

−

Ny2 µ2y 3

,

k = q.

Proof. From (2.7)–(2.8), the second-order derivatives are straightforward. To deduce the second-order CSDM, we just consider the case of C2x in the x-direction and divide the demonstration into two stages. The other case follows the same procedure. (i) Let x = xp+ 1 ̸ = xj+ 1 . From the definitions of Lemma 2.1, one has 2

2

sin(Nx θx ) = sin(Nx θx′ ) = 0, cos(Nx θx ) = (−1)j+p , cos(Nx θx′ ) = (−1)j+p+1 . Substituting the above results into X′′

j+ 21

X′′

j+ 21

(xp+ 1 ) = (−1)j+p

µ2x ( 2

2

j+p

= (−1)

csc2

µ2x [ 2

csc

2

θx′ 2

(x), it is not difficult to check

− csc2

( µx hx 2

)

µ2x 4Nx

j+ 12

2

sin(Nx θx ) cos

θx 2

−

Nx µ2x 2

2

(j + p + 1) − csc2

(ii) Let x = xp+ 1 = xj+ 1 . The first part of X′′ 2

θx )

sin(Nx θx ) cos

( µx hx 2

(j − p)

)]

.

¨ (Nx , µx , θx ) equals to (x), i.e. C θx 2

sin2

θx 2

−

µ2x 2

cos(Nx θx ) sin

sin3 θ2x

θx 2

+

Nx µ2x 2

cos(Nx θx ).

We apply the Taylor’s expansion for each item of numerator and obtain sin(Nx θx ) cos

( )( ) ( θx )2 (Nx θx )3 = Nx θx − + ··· 1 − 2 + ··· 2 3! 2! ( θx )3 4Nx3 ( θx )3 ( ) = Nx θx − Nx − + O θx5 ,

θx

2

sin(Nx θx ) cos

θx 2

sin2

3 2 (Nx θx )3

)( )( θ )2 ( θ x )2 ( θx )3 x + ··· 1 − 2 + ··· − 2 + ··· 2 3! 2! 2 3! ( θx )3 ( 5) = 2Nx + O θx ,

θx

(

= Nx θx − 2

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

cos(Nx θx ) sin

5

( )( θ ) ( θx )3 (Nx θx )2 x = 1− + ··· − 2 + ··· 2 2! 2 3! ( ) ( ) θx 1 ( θx )3 2 θx 3 = − − 2Nx + O θx5 .

θx

2

6 2

2

Combining the above four formulas together, since θx → 0 as xp+ 1 → xj+ 1 , we find 2

C¨ (Nx , µx , θx ) → −

µ2x

j+ 12

C¨ (Nx , µx , θx′ ) = −

+

6

For the second part of X′′

µ2x 2

Nx µ2x 2

−

Nx2 µ2x 3

,

as xp+ 1 → xj+ 1 . 2

2

¨ (Nx , µx , θx′ ), due to sin (x), i.e. C θx′

(−1)j+p+1 csc2

Nx µ2x

+

2

2

2

θx′ 2

̸= 0, it follows that

(−1)j+p+1 =

µ2x

1

csc2 (j +

(

2

2

)µx hx −

)

Nx µ2x 2

.

¨ (Nx , µx , θx ) with C¨ (Nx , µx , θx′ ), for j = p, one has Combining C X′′

j+ 21

µ2x

(xp+ 1 ) =

csc2 (j +

(

2

2

1 2

)µx hx −

)

µ2x 6

−

Nx2 µ2x 3

.

(i) and (ii) complete the proof and lead to the desired assertions in Lemma 2.2.

□

2.2. Properties of the CSDM In fact, the higher-order CSDM can be deduced in the same spirit. Unfortunately, in view of the process of derivation and expansion described above, it will be tedious and expensive. Hence, the investigation of the relations between the high-order CSDM and the second-order one will be an important topic. We illustrate this result as follows. y

x and C2k as the 2k-order CSDM in the x- and y-directions, respectively. Under the definitions of Lemma 2.3. Denote C2k Lemma 2.1, there exists y

x (C2x )k = C2k

y

for k ≥ 1, k ∈ N.

and (C2 )k = C2k

(2.11)

Proof. We use the mathematical induction to demonstrate that the first part of (2.11) holds for 1 ≤ k ≤ n. Always, the second one can be obtained correspondingly. Obviously, (2.11) is valid for k = 1. Assume that (2.11) holds for 1 ≤ k ≤ n − 1, that is, x (C2x )n−1 = C2(n −1) ,

consequently, we aim to prove that (2.11) still holds for k = n. Recalling (2.5) and differentiating 2(n − 1) times, we choose x = xp+ 1 and obtain 2

x C2(n −1)

(

)

j,p

=

Nx −1 2 ∑[

Nx

(−1)n−1

(µx l)2(n−1)

(

Zl

l=0

)] .

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

)

(

2

2

(2.12)

According to the hypothesis and (2.12) with n = 2, it holds that Nx −1

(C2x )n

(

)

j,p

=

=

Nx −1

∑(

(C2x )n−1

ν=0 Nx −1{ ∑

2 ∑[ Nx

C2

ν,p

(−1)n−1

=

∑(

−

(µx ℓ)2 Zℓ

ℓ=0

x C2(n −1)

) ( x)

ν=0

(µx l)2(n−1) Zl

l=0

Nx −1

×

j,ν

Nx −1 2 ∑[

Nx

ν=0

) ( x)

j,ν

C2

ν,p

cos µx l(xj+ 1 − a) cos µx l(xν+ 1 − a)

(

)

(

2

2

} )]

cos µx ℓ(xν+ 1 − a) cos µx ℓ(xp+ 1 − a)

(

)

)]

(

2

2

.

After a simple calculation, one has (C2x )n

(

)

j,p

=

Nx −1 Nx −1{ 2 ∑ ∑ [

Nx

×

ℓ=0 Nx −1 ∑

(−1)n

(µx l)2(n−1) (µx ℓ)2 Zl

l=0

2 Nx

ν=0

Zℓ

cos µx l(xj+ 1 − a) cos µx ℓ(xp+ 1 − a)

(

2

cos µx l(xν+ 1 − a) cos µx ℓ(xν+ 1 − a)

(

)

2

)

(

2

} )

(

)]

2

.

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

6

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

Notice that Nx −1

hx

∑ ν=0

⎧ Lx , l = ℓ = 0 , ) ( ) ⎨ ( cos µx l(xν+ 1 − a) cos µx ℓ(xν+ 1 − a) = Lx /2, l = ℓ ̸ = 0, 2 2 ⎩ 0, l ̸ = ℓ,

it follows that (C2x )n

(

)

j,p

=

=

Nx −1 2 ∑[

Nx

(µx l)2(n−1) (µx l)2

(−1)

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

(

Zl

l=0

Nx −1 2 ∑[

Nx

(−1)n

(µx l)2n

2n 2

Zl

l=0

)

(

2

2

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

(

)]

)

(

2

)]

2

( x) . = C2n j,p

This implies that (2.11) is still valid for k = n. Thus, the principle of induction confirms the result and completes the proof of Lemma 2.3. □ y

Remark 1. As described earlier, with the help of Lemmas 2.2 and 2.3, we can derive C4x and C4 to discretize the fourth-order derivative in space. This feature is attractive in further developing the spectral-collocation schemes for the reaction–diffusion equations, such as Cahn–Hilliard and extended Fisher–Kolmogorov models, with the same boundary conditions. Next, we shall reveal the relations between the CSDM and discrete cosine transform in the following lemma. It will be helpful for developing the fast algorithm in computation. y

Lemma 2.4. The eigenvalues of the second-order CSDM C2x and C2 are

λC2x ,ν+1 = −ν 2 µ2x and λC y ,ω+1 = −ω2 µ2y for (ν, ω) ∈ Th .

(2.13)

2

For any k ≥ 1, k ∈ N, the CSDM have the following decomposition x C2k = CNTx (Λx )k CNx

y

C2k = CNTy (Λy )k CNy ,

and

(2.14)

where CN is the discrete cosine transform with the elements

√ (CNx )ν,p =

2

cos

Nx Zν

πν (p + 12 ) Nx

√ ,

(CNy )ω,q =

2 Ny Zω

πω(q + 12 )

cos

Ny

,

CNT (the transpose matrix of CN ) is the discrete cosine inverse transform, and

Λx = diag(λC2x ,1 , λC2x ,2 , . . . , λC2x ,Nx ),

Λy = diag(λC y ,1 , λC y ,2 , . . . , λC y ,Ny ). 2

2

2

Proof. Taking (2.12) with n = 2, we find (C2x )j,p =

Nx −1 2 ∑[

Nx

−

(µx l)2 Zl

l=0

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

(

)

(

2

2

)] .

Let

√ (CNx )ν,p = (CNTx )p,ν =

√

2 Nx Zν

( ) cos µx ν (xp+ 1 − a) = 2

2 Nx Zν

πν (p + 12 )

cos

Nx

,

then one has

√

Nx −1

(Λx CNx )ν,p =

∑

2

(Λx )ν,l (CNx )l,p = (Λx )ν,ν

Nx Zν

l=0

cos

πν (p + 21 ) Nx

.

Also, it is not difficult to find Nx −1

(CNTx Λx CNx )j,p

=

∑

Nx −1

(CNTx )j,l (Λx CNx )l,p

l=0

=

√

∑

2 Nx Zl

l=0

cos

π l(j + 12 ) Nx

√ (Λx )l,l

Nx −1

=

2 ∑ [ (Λx )l,l Nx

l=0

Zl

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

(

)

2

(

2

)]

2 Nx Zl

cos

π l(p + 12 ) Nx

= (C2x )j,p .

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Overall, we obtain the first assertion of (2.14) with k = 1. Invoking to the orthogonality (CNTx CNx )j,p

Nx −1 2 ∑[ 1

=

Nx

cos µx l(xj+ 1 − a) cos µx l(xp+ 1 − a)

(

Zl

l=0

)

(

2

)]

2

= δjp ,

Lemma 2.4 leads to x = (C2x )k = CNTx Λx CNx CNTx Λx CNx · · · CNTx Λx CNx = CNTx (Λx )k CNx . C2k

The other part of (2.14) can be proved similarly. This completes the proof of Lemma 2.4. □ 3. Numerical scheme and discrete conservation laws 3.1. Cosine spectral-collocation difference scheme With all above preparations, applying the Crank–Nicolson finite difference scheme in time and the cosine spectralcollocation approximation in space, we present the cosine spectral-collocation difference (CSCD) scheme for the original problem (1.1) 1

1

n+ 12

y

iδt+ Ujn,k + α (C2x U n+ 2 )j,k + α (U n+ 2 C2 )j,k + Vj,k Uj,k

β

+

2

2

2

n+ 12

(|Ujn,k+1 | + |Ujn,k | )Uj,k

= 0,

n+ 21

= 12 (Ujn,k+1 + Ujn,k ) for (j, k) ∈ Th and 0 ≤ n ≤ Nτ − 1. Denote } { y Vh := φ|φ = φ (xj+ 1 , yk+ 1 ) for (j, k) ∈ Th and ∆c φ := C2x φ + φ C2 , the fully-discrete scheme can be rewritten in the

where Ujn,k = U(xj+ 1 , yk+ 1 , tn ), Vj,k = V (xj+ 1 , yk+ 1 ), Uj,k 2

2

2

2

2

2

following matrix form 1

1

iδt+ U n + α ∆c U n+ 2 + VU n+ 2 +

β

2

1

2

(|U n+1 | + |U n | )U n+ 2 = 0. 2 For any φ, ψ ∈ Vh , we define the discrete inner product and (semi-)norms as Nx −1 Ny −1

⟨φ, ψ⟩ := hx hy

∑∑ j=0

√ ∥φ∥ := ⟨φ, φ⟩,

φj,k ψ j,k ,

Nx −1 Ny −1 ( ) 1p ∑∑ ∥φ∥p := hx hy |φj,k |p , j=0

k=0

∥∇c φ∥ :=

√

(3.1)

⟨−∆c φ, φ⟩,

k=0

∥φ∥∞ := max |φj,k |. (j,k)∈Th

3.2. Discrete conservation laws Corresponding to the conservation laws (1.4)–(1.5) derived by the continuous problem, the CSCD scheme preserves the mass and energy in discrete level. Before we illustrate the conservation, an auxiliary lemma is drawn. Lemma 3.1. For any numerical solution U n of the scheme, there exists 1

1

Im⟨∆c U n+ 2 , U n+ 2 ⟩ = 0,

1

Re⟨∆c U n+ 2 , δt+ U n ⟩ = −

1 (

2τ

) ∥∇c U n+1 ∥2 − ∥∇c U n ∥2 .

Proof. Resorting to the definitions of discrete inner product and (semi-)norms, we obtain 1

1

1

Im⟨∆c U n+ 2 , U n+ 2 ⟩ = −Im ∥∇c U n+ 2 ∥2 = 0.

(

)

1 Re⟨∆c U n+1 + ∆c U n , U n+1 − U n ⟩ 2τ ) 1 ( = Re⟨∆c U n+1 , U n+1 ⟩ − Re⟨∆c U n , U n ⟩ 2τ ) 1 ( =− ∥∇c U n+1 ∥2 − ∥∇c U n ∥2 . 2τ This completes the proof of Lemma 3.1. □ 1

Re⟨∆c U n+ 2 , δt+ U n ⟩ =

Based on this lemma, we state the discrete conservation law theorem below. Theorem 3.1.

The CSCD scheme satisfies the following discrete conservation laws

Q := ∥U ∥ ≡ Q 0 , n

n 2

E n := α∥∇c U n ∥2 −

(3.2)

β 2

Nx −1 Ny −1

∥U n ∥44 − hx hy

∑∑ j=0

2

Vj,k |Ujn,k | ≡ E 0 .

(3.3)

k=0

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

8

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx 1

Proof. Computing the discrete inner product of (3.1) with U n+ 2 , then taking the imaginary part, one has 1 (

) ∥U n+1 ∥2 − ∥U n ∥2 = 0,

2τ

where Lemma 3.1 was used. It leads to (3.2) immediately. Computing the discrete inner product of (3.1) with δt+ U n , then taking the real part, one has

−

) ) α( β ( n+1 4 ∥∇c U n+1 ∥2 − ∥∇c U n ∥2 + ∥U ∥4 − ∥U n ∥44 4τ Nx −1 Ny −1 ∑∑ ( ) 1 2 2 + hx hy Vj,k |Ujn,k+1 | − |Ujn,k | = 0. 2τ

2τ

j=0

k=0

It yields E n+1 = E n , where E n was defined in (3.3). This completes the proof of Theorem 3.1.

□

4. Implementation and numerical tests 4.1. Fast algorithm Based on Lemma 2.4, the fully-discrete scheme (3.1) can be decomposed (k = 1) as iU n+1(s+1) +

ατ 2

(CNTx Λx CNx U n+1(s+1) + U n+1(s+1) CNTy Λy CNy ) = Q n ,

(4.1)

where

ατ

τ

βτ

2

2

∆c U n − V (U n+1(s) + U n ) − (|U n+1(s) | + |U n | )(U n+1(s) + U n ), 2 2 4 with the initial iteration and the iterative conditions satisfied Q n = iU n −

U

n+1(0)

{ =

U n,

n=0

2U n − U n−1 ,

n≥1

,

∥U n+1(s+1) − U n+1(s) ∥∞ ≤ 10−12 .

In consideration of the fast cosine transform, we denote ˜ U = CNx UCNTy . Multiplying CNx and CNTy on both sides of (4.1), an accelerated iterative algorithm is developed by (i +

ατ 2

Λx +

ατ 2

Λy )˜ U n+1(s+1) = ˜ Q n.

(4.2)

UCNy . Calculating Eq. (4.2), we can obtain the numerical solution U n ultimately with the help of the relation U = CNTx ˜ 4.2. Numerical examples Example 1. We choose the parameters α = −1, β = 0, V (x, y) = 0 and consider the linear Schrödinger equation iut = uxx + uyy on Ω = [0, π] × [0, π]. The analytical solution is u(x, y, t) = e2it cos x cos y. To quantify the error of numerical solution, we define the error Error(τ , h) = ∥u(T ) − U Nτ ∥∞ , and the convergence order of time Order ≈ log2 Error(τ , h)/Error(τ /2, h) .

(

)

For the convergence tests, we take T = 1 and different τ with h = 0.05π to investigate the temporal convergence rate. Numerical results are reported on the top row of Table 1. Similarly, we take a fine time step τ = 0.00001 to eliminate the temporal discretization error for testing the spatial accuracy. The results are displayed on the top row of Table 2. Moreover, two tables also record the CPU times of three examples and the maximum iterations of Examples 2 and 3. Discrete conservation laws are verified with τ = 0.05, h = 0.25π until T = 20. Table 3 indicates that the CSCD scheme maintains the mass and energy in discrete level. Example 2. The 2D NLS equation iut − uxx − uyy + V (x, y)u + 4|u|2 u = 0 is concerned on Ω = [−π, π] × [−π, π] with the potential function V (x, y) = 1 − 4 cos2 2x cos2 2y. It admits a breather solution u(x, y, t) = e9it cos 2x cos 2y. We choose T = 1 and h = 0.1π for different τ to study the temporal convergence order. Numerical results are drawn in the middle row of Table 1. For spatial accuracy, we use the same time-step τ = 0.00001 for different h to avoid contamination of the temporal error and display the results in the middle row of Table 2. Table 4 shows the discrete conservation results with the parameters τ = 0.05, h = 0.25π until T = 20. Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Table 1 Temporal convergence orders of CSCD scheme for different τ .

τ

0.04

0.02

0.01

0.005

Example 1

Error Order CPU(s)

1.0591e−03 * 0.02

2.6496e−04 1.9990 0.03

6.6252e−05 1.9997 0.05

1.6564e−05 1.9999 0.10

Example 2

Error Order Iterations CPU(s)

8.6215e−02 * 10 0.17

2.1873e−02 1.9788 8 0.28

5.4882e−03 1.9947 7 0.50

1.3733e−03 1.9987 6 0.88

Example 3

Error Order Iterations CPU(s)

2.6395e−05 * 11 3.94

6.6000e−06 1.9997 9 6.23

1.6501e−06 1.9999 8 9.31

4.1253e−07 2.0000 6 15.84

Table 2 Spatial errors of CSCD scheme for different h with τ = 0.00001. h

0.2π

0.1π

0.05π

0.025π

Example 1

Error CPU(s)

6.4236e−11 18.40

7.1385e−11 39.90

7.3803e−11 55.30

8.6280e−11 168.38

Example 2

Error Iterations CPU(s)

6.0253e−09 3 110.28

5.4695e−09 3 169.07

5.8967e−09 3 422.98

5.9904e−09 3 999.44

h

0.5

0.25

0.125

0.0625

Error Iterations CPU(s)

6.7981e−05 3 391.19

9.4322e−10 3 1611.97

6.3908e−12 3 12851.46

3.2880e−12 3 49078.63

Example 3

Table 3 Discrete conservation of Example 1 at different times with τ = 0.05, h = 0.25π . t

Qn

En

t

Qn

En

2 4 6 8 10

2.467401100272361 2.467401100272376 2.467401100272388 2.467401100272405 2.467401100272427

4.934802200544725 4.934802200544754 4.934802200544778 4.934802200544813 4.934802200544856

12 14 16 18 20

2.467401100272446 2.467401100272459 2.467401100272478 2.467401100272496 2.467401100272514

4.934802200544894 4.934802200544919 4.934802200544957 4.934802200544993 4.934802200545030

Table 4 Discrete conservation of Example 2 at different times with τ = 0.05, h = 0.25π . t

Qn

En

t

Qn

En

2 4 6 8 10

9.869604401087315 9.869604401085294 9.869604401083276 9.869604401081272 9.869604401079229

83.891637409241184 83.891637409222994 83.891637409204833 83.891637409186799 83.891637409168410

12 14 16 18 20

9.869604401077211 9.869604401075218 9.869604401073204 9.869604401071182 9.869604401069164

83.891637409150221 83.891637409132329 83.891637409114182 83.891637409095992 83.891637409077845

Table 5 Discrete conservation of Example 3 at different times with τ = 0.05, h = 0.25. t

Qn

En

t

Qn

En

2 4 6 8 10

3.999999999998198 3.999999999996414 3.999999999994632 3.999999999992849 3.999999999991070

2.222222222219641 2.222222222217862 2.222222222216081 2.222222222214294 2.222222222212517

12 14 16 18 20

3.999999999989285 3.999999999987504 3.999999999985719 3.999999999983939 3.999999999982156

2.222222222210733 2.222222222208952 2.222222222207164 2.222222222205390 2.222222222203599

Example 3. Another 2D NLS equation iut + uxx + uyy + V (x, y)u + 2|u|2 u = 0 is considered with V (x, y) = 1 − 2 tanh2 x tanh2 y on Ω = [−30, 30] × [−30, 30]. We focus on the soliton wave solution u(x, y, t) = eit sech x sech y. Similar to the above two examples, we take T = 0.2 and h = 0.2 for different τ to investigate the convergence rate in time. The spatial errors for different h are computed with a sufficiently small time-step τ = 0.00001. Numerical results are also given in Tables 1 and 2. Identically, we choose τ = 0.05, h = 0.25, T = 20 to verify the discrete conservation laws and list the results in Table 5. Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

10

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. The interactions between two solitons at different times with τ = 0.1, h = 0.0625.

In summary, from Tables 1–5 and additional numerical results not shown here for brevity, the following observations are drawn: (1) The CSCD scheme is second-order of accuracy in time and spectral order of accuracy in space (the spatial errors in Table 2 are rapidly smaller and almost negligible). (2) The CSCD scheme still maintains the conservation of mass and energy in discrete level for a long time, no matter it is a linear problem or a nonlinear one. 5. Numerical applications In this section, we apply the CSCD scheme to address some practical problems so as to verify the performances. It is worthwhile to note that we always choose sufficiently large bounded rectangles as the computational domains. Since the solutions vanish at the boundary of computational domain, the Neumann boundary conditions are naturally satisfied. Example 4. The soliton interactions in optics can be modelled by the NLS equation. In order to simulate the interaction of two solitons, we numerically solve the following equation iut + uxx + uyy + V (x, y)u + 2|u|2 u = 0. The computational domain Ω = [−30, 60] × [−30, 60] and the initial condition reads u(x, y, 0) = e2i(x+y) sech x sech y + e0.05i(x+y−15) sech(x − 15) sech(y − 15). Fig. 1 describes the time evolutions of two solitons interactions with the parameters τ = 0.1, h = 0.0625 at different times t = 0, 2, 4, 6, 8, 10 Example 5. At extremely low temperatures, we use the Gross–Pitaevskii equation (also a NLS equation) to describe the dynamics of a Bose–Einstein condensates (BEC) [46], namely, 1 iut = − ∆u + V (x, y)u + β|u|2 u, 2 where the trapping potential V (x, y) =

∫ Q (t) = Ω

|u(x, t)|2 dx = 1.

1 2 (x 2

+ γy2 y2 ). The exact solution satisfies the following normalization condition (5.1)

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 2. Contour plots of |u(x, y, t)|2 of Example 5 at different times with τ = 0.01, h = 0.02.

Table 6 Relative errors of Q n at different times. t

2

4

6

8

10

Example 5 Example 6

3.3484e−13 7.6938e−13

4.6358e−12 1.0449e−12

4.9402e−12 1.1434e−12

6.2789e−12 2.0796e−12

5.0420e−12 1.3495e−12

In numerical computation, the condensate widths along two directions are approximated by

  Ny −1  N∑ x −1 ∑  2 n 2 √ (xj+ 1 − xˆ ) |Ujn,k |2 , σx = h 2

j=0

  N −1 y Nx −1  ∑ ∑  σyn = √h2 (yk+ 1 − yˆ )2 |Ujn,k |2 , 2

k=0

k=0

j=0

where Ny −1

Nx −1

xˆ = h2

∑

xj+ 1

∑

2

j=0

Ny −1

|Ujn,k |2 and yˆ = h2

k=0

∑

Nx −1

yk+ 1

∑

2

k=0

|Ujn,k |2 .

j=0

In addition, we always verify the normalization quantity numerically by Nx −1 Ny −1

Q n = h2

∑∑ j=0

|Ujn,k |2 .

k=0

We take the parameters β = 2, γy = 2 and the initial solution u(x, y, 0) =

1/4

γ √y

2π

e−0.25(x

2 +γ y2 ) y

. The computational

domain Ω = [−8, 8] × [−8, 8] and the finial time T = 10. The time evolutions of the condensate density |u(x, y, t)|2 are plotted in Fig. 2 with the time-step τ = 0.01 and spacing size h = 0.02. Fig. 3 (left) shows the condensate widths along the x- and y-directions, which are coincides with the results in [46]. We also list the relative errors of the normalization quantity Q n at different times, see the second row of Table 6. It indicates that they are corresponding to the continuous case (5.1) well. Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

12

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 3. Condensate widths of Example 5 (left) and Example 6 (right) at different times.

Example 6. We consider the 2D NLS equation with a rotational term 1 iut = − ∆u + V (x, y)u + β|u|2 u − iΥ (y∂x − x∂y )u, 2 which describes the dynamic properties of one-component BEC under external rotation (Υ > 0). The trapping potential V (x, y) = 12 (γx2 x2 + γy2 y2 ) and the initial data reads u(x, y, 0) =

( γx γy ) 14 −0.5(γ 2 x2 +γ 2 y2 ) x y e . π2

In consideration of the rotational term, the CSCD method can not be applied directly. Invoking to [47], we eliminate the rotational term by introducing the rotating Lagrangian coordinates

( A(t) =

cos(Υ t) − sin(Υ t)

sin(Υ t) . cos(Υ t)

)

The original problem becomes the usual 2D NLS equation in rotating Lagrangian coordinates 1 2 iu˜ t = − ∆u˜ + W (x˜ , y˜ , t)u˜ + β|˜u| u˜ , 2 where W (x˜ , y˜ , t) =

γx2 + γy2 4

(x˜ 2 + y˜ 2 ) +

γx2 − γy2 [ 4

(x˜ 2 − y˜ 2 ) cos(2Υ t) + 2x˜ y˜ sin(2Υ t) ,

]

and the initial condition u˜ (x˜ , y˜ , 0) = u(x, y, 0). 2 2 We here set γx = γy = 1 and the initial solution becomes u˜ (x˜ , y˜ , 0) = √1π e−0.5(x˜ +˜y ) . It is easy to check that the equation described above falls into the category of Example 5. Hence, we similarly take β = 2 and Ω = [−8, 8] × [−8, 8] for simulation. The time-step and spacing size are chosen as τ = 0.02 and h = 0.02, respectively. Fig. 4 displays the 2 time evolutions of the condensate density |˜u(x˜ , y˜ , t)| until T = 10. The condensate widths along two directions are also drawn in Fig. 3 (right), which are agreement with the results in [46]. The bottom row of Table 6 reports the relative errors of the normalization quantity Q n in this situation. Always, they are approximate to the value of 0 which confirms the normalization in (5.1). The results of all these three applications witness the reliability and validity of the CSCD method. 6. Conclusion In view of the second-order CSDM, we present a conservative CSCD scheme for 2D Schrödinger equation with Neumann boundary conditions. The derivation of high-order CSDM and its relationship with discrete cosine transform are also displayed in this work. We can use these results to discretize the high-order derivatives in space. With the help of the above relations and fast cosine transform, an accelerated iterative algorithm is developed for the computation. Numerical tests and applications are carried out to verify the performances of the new scheme. All these results confirm the effectiveness and high accuracy of the CSCD method. Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

13

2

Fig. 4. Contour plots of |˜u(x˜ , y˜ , t)| of Example 6 at different times with τ = 0.02, h = 0.02.

Acknowledgements The authors would like to thank the anonymous referees for their careful reading and many valuable comments and suggestions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

D.J. Griffiths, Introduction to Quantum Mechanics, Prentice Hall, New Jersey, 1995. A. Hasegawa, Optical Solitons in Fibers, Springer, Berlin, 1989. C.R. Menyuk, Stability of solitons in birefringent optical fibers, I: Equal propagation amplitudes, Opt. Lett. 12 (1987) 614–616. C.R. Menyuk, Stability of solitons in birefringent optical fibers, II: Arbitrary amplitudes, J. Opt. Soc. Amer. B 5 (1988) 392–402. L.P. Pitaevskii, S. Stringary, Bose–Einstein Condensation, Clarendon Press, New York, 2003. J.R. Abo-Shaeer, C. Raman, J.M. Vogels, W. Ketterle, Observation of vortex lattices in Bose–Einstein condensates, Science 292 (2001) 476. V.E. Zakharov, V.S. Synakh, The nature of self-focusing singularity, Sov. Phys. JETP 41 (1975) 465–468. C. Sulem, P.L. Sulem, The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse, Springer, New York, 1999. W.Z. Bao, Q. Du, Y.Z. Zhang, Dynamics of rotating Bose–Einstein condensates and its efficient and accurate numerical computation, SIAM J. Appl. Math. 66 (2006) 758–786. T. Cazenave, Semilinear Schrödinger Equations, in: Courant Lecture Notes in Mathematics, vol. 10, New York University, Courant Institute of Mathematical Sciences, AMS, 2003. V.G. Markowich, Dynamics of classical solitons (in non-integrable systems), Phys. Lett. C 35 (1978) 1–128. L.R.T. Gardner, G.A. Gardner, S.I. Zaki, Z. EI Sahrawi, B-spline finite element studies of the non-linear Schrödinger equation, Comput. Methods Appl. Mech. Engrg. 108 (1993) 303–318. O.A. Karakashian, G.D. Akrivis, V.A. Dougalis, On optimal order error estimates for the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 30 (1993) 377–400. G.D. Akrivis, V.A. Dougalis, O.A. Karakashian, On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schrödinger equation, Numer. Math. 59 (1991) 31–53. O.A. Karakashian, C. Makridakis, A space–time finite element method for the nonlinear Schrödinger equation: the discontinuous Galerkin method, Math. Comp. 67 (1998) 479–499. Y. Xu, C.W. Shu, Local discontinuous Galerkin methods for nonlinear Schrödinger equations, J. Comput. Phys. 205 (2005) 72–77. M. Dehghan, D. Mirzaei, Numerical solution to the unsteady two-dimensional Schrödinger equation using meshless local boundary integral equation method, Internat. J. Numer. Methods Engrg. 76 (2008) 501–520. M. Dehghan, D. Mirzaei, The meshless local Petrov–Galerkin (MLPG) method for the generalized two-dimensional non-linear Schrödinger equation, Eng. Anal. Bound. Elem. 32 (2008) 747–756. F. Zhang, V.M. Pérez-Gratciz, L. Vázquez, Numerical simulation of nonlinear Schrödinger systems: a new conservative scheme, Appl. Math. Comput. 71 (1995) 165–177. H.Q. Wang, Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations, Appl. Math. Comput. 170 (2005) 17–35.

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.

14

X. Li and L. Zhang / Computers and Mathematics with Applications xxx (xxxx) xxx

[21] M. Subasi, On the finite difference schemes for the numerical solution of two dimensional Schrödinger equations, Numer. Methods Partial Differential Equations 18 (2002) 752–758. [22] Z. Gao, S.S. Xie, Fourth-order alternating direction implicit compact finite difference schemes for two-dimensional Schrödinger equations, Appl. Numer. Math. 61 (2011) 593–614. [23] Y.Q. Xu, L.M. Zhang, Alternating direction implicit method for solving two-dimensional cubic nonlinear Schrödinger equation, Comput. Phys. Comm. 183 (2012) 1082–1093. [24] T.C. Wang, B.L. Guo, Q.B. 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. [25] M. Thalhammer, High-order exponential operator splitting methods for timedependent Schrödinger equations, SIAM J. Numer. Anal. 46 (2008) 2022–2038. [26] W.Z. Bao, J. Shen, A fourth-order time-splitting Laguerre–Hermite pseudospectral method for Bose–Einstein condensates, SIAM J. Sci. Comput. 26 (2005) 2010–2028. [27] D. Pathria, J.L. Morris, Pseudo-spectral solution of nonlinear Schrödinger equations, J. Comput. Phys. 87 (1990) 108–125. [28] M. Dehghan, A. Taleei, Numerical solution of nonlinear Schrödinger equation by using time-space pseudo-spectral method, Numer. Methods Partial Differential Equations 26 (2010) 979–992. [29] W.Z. Bao, H.L. Li, J. Shen, A generalized-Laguerre–Fourier–Hermite pseudospectral method for computing the dynamics of rotating Bose–Einstein condensates, SIAM J. Sci. Comput. 31 (2009) 3685–3711. [30] Y.Z. Gong, Q. Wang, Y.S. Wang, J.X. Cai, A conservative fourier pseudo-spectral method for the nonlinear Schrödinger equation, J. Comput. Phys. 328 (2017) 354–370. [31] 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. [32] Q.S. Chang, E.H. Jia, W. Sun, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys. 148 (1999) 397–415. [33] W.Z. 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. [34] M. Caliari, C. Neuhauser, M. Thalhammer, High-order time-splitting Hermite and Fourier spectral methods for the Gross–Pitaevskii equation, J. Comput. Phys. 228 (2009) 822–832. [35] W.Z. Bao, Y.Y. Cai, Mathematical theory and numerical methods for Bose–Einstein condensation, Kinet. Relat. Models 6 (2013) 1–135. [36] W.Z. Bao, Q.L. Tang, Z.G. Xu, Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation, J. Comput. Phys. 235 (2014) 423–445. [37] W.Y. Liao, J.P. Zhu, A.Q.M. Khaliq, A fourth-order compact algorithm for nonlinear reactiondiffusion equations with Neumann boundary conditions, Numer. Methods Partial Differential Equations 22 (2006) 600–616. [38] J. Zhao, W.Z. Dai, T.C. Niu, Fourth-order compact schemes of a heat conduction problem with Neumann boundary conditions, Numer. Methods Partial Differential Equations 23 (2007) 949–959. [39] J. Zhao, W.Z. Dai, S.Y. Zhao, Fourth-order compact schemes for solving multidimensional heat problems with Neumann boundary conditions, Numer. Methods Partial Differential Equations 24 (2008) 165–178. [40] Z.Z. Sun, Compact difference schemes for heat equation with Neumann boundary conditions, Numer. Methods Partial Differential Equations 29 (2013) 1459–1486. [41] M. Javidi, A. Golbabai, Spectral collocation method for parabolic partial differential equations with Neumann boundary conditions, Appl. Math. Sci. 1 (2007) 211–218. [42] K. Atkinson, O. Hansen, D. Chien, A spectral method for elliptic equations: the Neumann problem, Adv. Comput. Math. 34 (2011) 295–317. [43] X.H. Yu, Z.Q. Wang, Jacobi spectral method with essential imposition of Neumann boundary condition, Appl. Numer. Math. 62 (2012) 956–974. [44] B. Smith, R. Laoulache, A. Heryudono, Implementation of Neumann boundary condition with influence matrix method for viscous annular flow using pseudospectral collocation, J. Comput. Appl. Math. 285 (2015) 100–115. [45] Y.Z. Gong, J.X. Cai, Y.S. Wang, Multi-symplectic Fourier pseudospectral method for the Kawahara equation, Commun. Comput. Phys. 16 (2014) 35–55. [46] W.Z. Bao, D. Jaksch, P.A. Markowich, Numerical solution of the Gross–Pitaevskii equation for Bose–Einstein condensation, J. Comput. Phys. 187 (2003) 318–342. [47] J. Ming, Q.L. Tang, Y.Z. Zhang, An efficient spectral method for computing dynamics of rotating two-component Bose–Einstein condensates via coordinate transformation, J. Comput. Phys. 258 (2014) 538–554.

Please cite this article as: X. Li and L. Zhang, An efficient spectral-collocation difference method for two-dimensional Schrödinger equation with Neumann boundary conditions, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.006.