High-order implicit Galerkin-Legendre spectral method for the two-dimensional Schrödinger equation1footnote 11footnote 1“The work of authors was partially supported by the National Natural Science Foundation of China (11271100, 11301113, 51476047, 11601103) and the Fundamental Research Funds for the Central Universities (Grant No. HIT. IBRSEM. A. 201412), Harbin Science and Technology Innovative Talents Project of Special Fund (2013RFXYJ044).”

# High-order implicit Galerkin-Legendre spectral method for the two-dimensional Schrödinger equation111“The work of authors was partially supported by the National Natural Science Foundation of China (11271100, 11301113, 51476047, 11601103) and the Fundamental Research Funds for the Central Universities (Grant No. HIT. IBRSEM. A. 201412), Harbin Science and Technology Innovative Talents Project of Special Fund (2013RFXYJ044).”

Wenjie Liu, Boying Wu Department of Mathematics, Harbin Institute of Technology, 150001, China
###### Abstract

In this paper, we propose Galerkin-Legendre spectral method with implicit Runge-Kutta method for solving the unsteady two-dimensional Schrödinger equation with nonhomogeneous Dirichlet boundary conditions and initial condition. We apply a Galerkin-Legendre spectral method for discretizing spatial derivatives, and then employ the implicit Runge-Kutta method for the time integration of the resulting linear first-order system of ordinary differential equations in complex domain. We derive the spectral rate of convergence for the proposed method in the -norm for the semidiscrete formulation. Numerical experiments show our formulation have high-order accurate.

###### keywords:
Two-dimensional Schrödinger equation, Galerkin-Legendre spectral method, Implicit Runge-Kutta method, Error estimate
journal:

## 1 Introduction

In this paper, we introduce Galerkin-Legendre spectral method for the two-dimensional Schrödinger equation

 −i∂u∂t=Δu+ψ(x,y)u, (x,y,t)∈(c,d)×(c,d)×(t0,T), (1.1a) with the initial condition u(x,y,t0)=φ(x,y), (x,y)∈Ω=(c,d)×(c,d), (1.1b) and the Dirichlet boundary conditions u(x,c,t)=g1(x,t), u(x,d,t)=g2(x,t),(x,t)∈(c,d)×(t0,T),u(c,y,t)=g3(y,t), u(d,y,t)=g4(y,t),(y,t)∈(c,d)×(t0,T), (1.1c)

where is a complex-valued function, is a smooth known real function, and are smooth known complex functions, and .

The Schrödinger equation is a famous equation used widely in many fields of physics Mohebbi2009paper (); Abdur2011paper (); Dehghan2013paper (), such as quantum mechanics, quantum dynamics calculations, optics, underwater acoustics, plasma physics and electromagnetic wave propagation.

Over the past few years, several numerical schemes have been developed for solving problem (1.1). For instance, In Subasi2002paper () Subasi used a standard finite difference method in space for solving the problem (1.1). A Crank-Nicolson method discretization for the time of the problem (1.1) was considered in in Antoine2004paper (). An implicit semi-discrete higher order compact (HOC) scheme was considered for computing the solution of the problem (1.1) in Kalita2006paper (). Dehghana and Shokriin Dehghan2007paper () proposed a numerical scheme to solve the problem (1.1) by collocation and radial basis functions in space. A meshless local boundary integral equation (LBIE) method to solve the problem (1.1) was given in Dehghan2008paper (). In Mohebbi2009paper () Mohebbi and Dehghan presented a high order method for the problem (1.1) by the compact finite difference in space and boundary value method in time.In Tian2010paper () Tian and Yu proposed a HOC-ADI method to solve the problem (1.1), which has fourth-order accuracy in space and second-order accuracy in time. Gao and Xie Gao2011paper () proposed a numerical scheme to solve the problem (1.1) by the ADI compact finite difference scheme. For the spatial discretisation of the problem (1.1) by the Chebyshev spectral collocation method considered in Abdur2011paper (). In Dehghan2013paper () Dehghan and Emami-Naeini presented Sinc-collocation and Sinc-Galerkin methods to solve the problem (1.1).

In this paper, we propose a numerical scheme for solving (1.1). We apply Galerkin-Legendre spectral method Shen2006Book (); Shen2011Book (); shj1 () for discretizing spatial derivatives, then using the implicit Runge-Kutta method for time derivatives, which is high-order accurate both in space and time.

The contents of the article is as follows. In Section 2, we introduce the implicit Runge-Kutta method for the system of ordinary differential equations in the complex domain. In Section 3, we describe the Galerkin-Legendre spectral method for the two-dimensional Schrödinger equation in space. In Section 4, we derive a priori error estimates in the -norm for the semidiscrete formulation. The analysis relies on an idea suggested by Lions et al. Lions1972Book () and Thomée Thomee2006Book (). In Section 5, we report the numerical experiments of solving two-dimensional Schrödinger equation with the new method developed in this paper, and compare the numerical results with analytical solutions and with other method in the literature Mohebbi2009paper (). We end this article with some concluding remarks in Section 6.

## 2 Implicit Runge-Kutta method

In this section we modifed the implicit Runge-Kutta (IRK) method to solve first-order linear complex ordinary differential equations. We first give the definition of Kronecker product of matrices.

###### Definition 2.1 (see, e.g., Steeb2006Book ()).

Let and be arbitrary matrices, then the matrix

 A⊗B=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝a11Ba12B…a1nBa21Ba22B…a2nB⋮⋮⋱⋮am1Bam2B…amnB⎞⎟ ⎟ ⎟ ⎟ ⎟⎠,

is called the Kronecker product of and .

###### Definition 2.2 (see, e.g., Steeb2006Book ()).

Let be any given matrix, then is defined to be a column vector of size made of the row of

 vec(A)=(a11,a12,…,a1n,a21,a22,…,am1,…,amn)T.
###### Lemma 2.1 (see, e.g., Steeb2006Book ()).

Let , , and be three given matrices. Then

 vec(ABC)=(A⊗CT)vec(B).
###### Lemma 2.2 (see, e.g., Steeb2006Book ()).

Let , . Then

 vec(ABT)=A⊗B.

We consider the following first-order initial value problem is given by

 {y′=f(t,y(t)),t0

A -stage formula of implicit Runge-Kutta method for approximating (2.1) can be written as

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩yn+1=yn+hs∑l=1blf(tn+clh,kl),kl=yn+hs∑j=1aljf(tn+cjh,kj),l=1,2,…,s, (2.2)

where is the integration step, are the internal stages and . If the method is called diagonally implicit Runge-Kutta (DIRK) method. Then the (2.2) can be written the following form

 {K=1syn+hAF(K),yn+1=yn+hBTF(K), (2.3)

where , , , and

 F(K)=[f(tn+c1h,k1),f(tn+c2h,k2),⋯,f(tn+csh,ks)]T.

If we consider the system of linear ordinary differential equations

 {−iM1y′+M2y−f(t)=0,  t0

where , and

 y(t)=[y1(t),y2(t),⋯,ym(t)]T,f(t)=[f1(t),f2(t),⋯,fm(t)]T.

Using s-stage formula of Runge-Kutta method (2.3) for approximating (2.4) can be written as

 {−iKMT1=−i1syTnMT1+hA(−KMT2+F),−iyTn+1MT1=−iyTnMT1+hBT(−KMT2+F), (2.5)

where

 K=[kT1;kT2;⋯;kTs],F=[f(tn+c1h)T;f(tn+c2h)T;⋯;f(tn+csh)T].

By using Lemma 2.1 and 2.2, we have

 (−iIs⊗M1+hA⊗M2)vec(K)=−i1s⊗(M1yn)+hA⊗Imvec(F), (2.6)

where is the unitary matrix. Thus we can obtain from (2.6) with GMRES iteration (see, e.g., (Shen2006Book, , PP. 55–58)).

We write the Runge-Kutta scheme in a tabular format known as the Butchers table, in Table 2.1, we choose a 3-stage IRK method (cf. Butcher ()).

## 3 Discretize two-dimensional Schrödinger equation in space by Galerkin-Legendre spectral method

In this section, we will present the Galerkin-Legendre spectral method to solve the unsteady two-dimensional Schrödinger equation for the space with the nonhomogeneous Dirichlet boundary conditions and initial condition. Firstly, we make the variable transformations

 x=c+~x+12(d−c), y=c+~y+12(d−c), γ=(2d−c)2, ~u(~x,~y,t)=u(c+~x+12(d−c),c+~y+12(d−c),t), ~ψ(~x,~y)=ψ(c+~x+12(d−c),c+~y+12(d−c)), ~φ(~x,~y)=φ(c+~x+12(d−c),c+~y+12(d−c)), ~g1(~y,t)=g1(c+~y+12(d−c),t), ~g2(~y,t)=g2(c+~y+12(d−c),t), ~g3(~x,t)=g3(c+~x+12(d−c),t), ~g4(~x,t)=g4(c+~x+12(d−c),t).

Then is changed to the square , and the (1.1) can be rewritten as

 −i∂~u∂t=γΔ~u+~ψ(~x,~y)~u,(~x,~y,t)∈˜Ω×(t0,T), (3.1a) with the initial condition ~u(~x,~y,t0)=~φ(~x,~y), on˜Ω, (3.1b) and the Dirichlet boundary conditions ~u(−1,y,t)=~g1(y,t), ~u(1,~y,t)=~g2(~y,t), on J=(t0,T),~u(~x,−1,t)=~g3(~x,t), ~u(~x,1,t)=~g4(~x,t), on J=(t0,T). (3.1c)

Now we recall the boundary conditions homogeneous process (see, e.g., shj1 ()). Setting

 ~u1(~x,~y,t)=~g4(~x,t)−~g3(~x,t)2~y+~g4(~x,t)+~g3(~x,t)2, ^g1(~y,t)=~g1(~y,t)−~u1(−1,~y,t), ^g2(~y,t)=~g2(~y,t)−~u1(1,~y,t), ~u2(~x,~y,t)=^g2(~y,t)−^g1(~y,t)2~x+^g2(~y,t)+^g1(~y,t)2, ^u=~u−~u1−~u2, ~f(~x,~y,t)=i(∂~u1∂t+∂~u2∂t)+γ(Δ~u1+Δ~u2)+~ψ(~u1+~u2), ^φ(~x,~y)=~φ(~x,~y)−~u1(~x,~y,t0)−~u2(~x,~y,t0),

then the (3.1) can be rewritten as

 −i∂^u∂t=γΔ^u+~ψ(~x,~y)^u+~f(~x,~y,t), (~x,~y,t)∈(−1,1)×(−1,1)×(t0,T), (3.2a) with the initial condition and homogeneous boundary value conditions ^u(~x,~y,t0)=^φ(~x,~y),on˜Ω,^u(~x,−1,t)=^u(~x,1,t)=^u(−1,~y,t)=^u(1,~y,t)=0,onJ=(t0,T). (3.2b)

We shall now discretize the equation (3.2) by using the Galerkin-Legendre spectral method in space. Let us denote the th degree Legendre polynomial (see, e.g., (Shen2006Book, , PP. 18 and 19)) and

 PN=\textupspan{L0(~x),L1(~x),⋯,LN(~x)}, VN={v∈PN:v(±1)=0}.

Then the semi-discrete Legendre-Galerkin method for (3.2) is: find such that

 {−i(∂t^uN,v)+γ(∇^uN,∇v)−(~ψ(~x,~y)^uN,v)−(^f(~x,~y,t),v)=0,∀v∈V2N,(^uN(~x,~y,t0)−^u0,v)=0,∀v∈V2N, (3.3)

where is the scalar product in , the is complex conjugate of .

The following lemma is the key to implement our algorithms.

###### Lemma 3.1.

((shj1, , Lemma 2.1)) Let us denote

 ϕk(~x)=ck(Lk(~x)−Lk+2(~x)),ck=1√4k+6,
 ~ajk=∫1−1ϕ′k(~x)ϕ′j(~x)d~x,~bjk=∫1−1ϕk(~x)ϕj(~x)d~x,

then

Obviously

 VN=\textupspan{ϕ0(~x),ϕ1(~x),⋯,ϕN−2(~x)}.

Let us setting

 ^uN(~x,~y,t)=N−2∑k,j=0αkj(t)ϕk(~x)ϕj(~y), (3.4)

then inserting (3.4) into (3.3) and taking

 v=ϕl(~x)ϕm(~y),

yields

 −i(∂t^uN,ϕl(~x)ϕm(~y)) +γ(∇^uN,∇(ϕl(~x)ϕm(~y)))−(~ψ(~x,~y)^uN,ϕl(~x)ϕm(~y)) (3.5) −(^f(~x,~y,t),ϕl(~x)ϕm(~y))=0,l,m=0,1,…,N−2.

Denote and

 ˜B=(~bkj)k,j=0,1,…,N−2,^fkj=(^f(~x,~y,t),ϕk(~x)ϕj(~y)),ˆF(t)=(^fkj)k,j=0,1,…,N−2.

The (3.5) is equivalent to the following matrix equation

 −i˜Bα′˜B+γ(˜BαIN−1+IN−1α˜B)−Wα−ˆF(t)=0,

where and

 (Wα)lm=N−2∑k,j=0αkj(t)(~ψ(~x,~y)ϕk(~x)ϕj(~y),ϕl(~x)ϕm(~y)), l,m=0,1,…,N−2.

Using Lemma 2.1, we have

 −i˜B⊗˜Bvec(α′)+(γ˜B⊗IN−1+γIN−1⊗˜B−˜W)vec(α)−vec(ˆF(t))=0,

where

Setting

 ^uN(~x,~y,t0)=N−2∑k,j=0αkj(t0)ϕk(~x)ϕj(~y),

and with

 (^uN(~x,~y,t0)−^φ(~x,~y),ϕl(~x)ϕm(~y))=0,l,m=0,1,…,N−2,

we obtian

 ˜Bα(t0)˜B=^u0,(^u0)kj=(φ(~x,~y),ϕk(~x)ϕj(~y)).

Using Lemma 2.1, we have

 ˜B⊗˜Bvec(α(t0))=vec(^u0).

Therefore, (3.3) leads to the following system of linear ordinary differential equations

 (3.6)

Hence we can use the 3-stage IRK method for (3.6).

## 4 A priori error estimate

In this section, we derive optimal a priori error bound for the semidiscrete scheme of the problem (1.1) by using the Galerkin-Legendre spectral method discretization for the space. Recall the spaces

 L2(˜Ω)={u:(u,u)<+∞},Hm(˜Ω)={u:Dku∈L2(˜Ω), 0≤|k|≤m},Hm0(˜Ω)={u∈Hm(˜Ω):Dku(∂˜Ω)=0, 0≤|k|≤m−1}.

The norms in and denoted by and , respectively, which are given as

Furthermore, we shall use to denote the semi-norm in . We now introduce the bochner space endowed with the norm

 ∥v∥Lp;X={(∫J∥v∥pXdt)1/p,1≤p<∞,\textupesssupt∈J∥v∥X,p=∞,

where or .

Let The orthogonal projection is defined by

 (Π0Nu−u,v)=0,∀v∈V2N.

The operators have the following approximation properties.

###### Lemma 4.1.

((Canuto1987Book, , PP. 309 and 310)) For any positive integer , the following estimate holds for any

 N∥u−Π0Nu∥+∥∇(u−Π0Nu)∥≤CN1−m∥u∥m,

where is independent of .

We now split the error as a sum of two terms

 ^uN−^u=(^uN−P0N^u)+(P0N^u−^u)=θ+ρ, (4.1)

where is an elliptic projection in of the exact solution , which defined by

 A(P0N^u−^u,v)=0,∀v∈V2N. (4.2)

We begin with the following auxiliary result.

###### Lemma 4.2.

((Canuto1987Book, , P. 309) Assume that and defined by (4.2). Then

 N∥P0N^u−^u∥+∥∇(P0N^u−^u)∥≤CN1−m∥^u∥m, \textupfor  t∈¯J,

where is independent of .

###### Lemma 4.3.

With defined by (4.2) and . Assume that . Then

 ∥ρ(t)∥≤CN−m∥^u∥m, \textupfor t∈¯J,
 ∥∂tρ(t)∥≤CN−m∥∂t^u∥m, \textupfor  t∈¯J,

where is independent of .

We are now ready for the -error estimate for the semidiscrete problem.

###### Theorem 4.1.

Let and be the solutions of (3.2) and (3.3), respectively. Assume that is a real function, and . Then, we have

 ∥^uN(t)−^u(t)∥≤CN−m(∥^u∥L∞(J;Hm(˜Ω))+∥∂t^u∥L∞(J;Hm(˜Ω))), \textupfor t∈J,

where is independent of .

###### Proof.

With and satisfies the following equation

 −i(∂t^uN,v)+γA(^uN,v)−(~ψ(~x,~y)^uN,v)−(^f(~x,~y,t),v)=0,∀v∈V2N,

and

 −i(∂t^u,v)+γA(^u,v)−(~ψ(~x,~y)^u,v)−(^f(~x,~y,t),v)=0,∀v∈V2N.

Hence with , we obtain

 −i(∂t(θ+ρ),v)+γA(θ+ρ,v)−(~ψ(~x,~y)(θ+ρ),v)=0,∀v∈V2N.

Taking both sides the imaginary, we obtain

 12ddt∥θ∥2=−Re(∂tρ,θ)−Im(~ψ(~x,~y)ρ,θ),

such that

 12ddt∥θ∥2≤C(∥θ∥2+∥ρ∥2+∥∂tρ∥2).

After integration, this shows

 ∥θ(t)∥2≤C∥θ(t0)∥2+C∫tt0(∥θ∥2+∥ρ∥2+∥∂sρ∥2)ds.

By using the Gronwall’s Lemma, we have

 ∥θ(t)∥2≤C∥θ(t0)∥2+C∫tt0(∥ρ∥2+∥∂sρ∥2)ds.

By Lemma 4.3, we have

 ∥θ(t)∥2≤C∥θ(t0)∥2+CN−2m(T−t0)(∥^u∥2L∞(J;Hm(˜Ω))+∥∂t^u∥2L∞(J;Hm(˜Ω)))≤C5N−2m(∥^u∥2L∞(J;Hm(˜Ω))+∥∂t^u∥2L∞(J;Hm(˜Ω))).

Using Lemmas 4.1 and 4.2, we obtain

 ∥θ(t0)∥=∥^uN(t0)−P0N^u(t0)∥≤∥Π0N^u(t0)−^u(t0)∥+∥P0N^u(t0)−^u(t0)∥≤CN−m(∥^u∥