以半顯隱Runge-Kutta及擬譜法求解薛丁格方程

國

立

交

通

大

學

應用數學系數學建模與科學計算碩士班

碩

士

論

文

以半顯隱

Runge-Kutta

及擬譜法求解薛丁格方程

Implicit-Explicit Runge-Kutta and pseudospectral methods for Schrodinger

equation

研 究 生：林沛沅

指導教授：賴明治 博士

共同指導教授：鄧君豪 博士

以半顯隱

Runge-Kutta

及擬譜法求解薛丁格方程

Implicit-Explicit Runge-Kutta and pseudospectral methods for Schrodinger

equation

研 究 生：林沛沅 Student：Pei-Yuan Lin

指導教授：賴明治 Advisor：Ming-Chih Lai

共同指導教授：鄧君豪 Co-Advisor：Chun-Hao Teng

國 立 交 通 大 學

應 用 數 學 系

數 學 建 模 與 科 學 計 算 碩 士 班

碩 士 論 文

Submitted to Department of Applied Mathematics College of Science Institute of Mathematical Modeling and Scientific Computing

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Applied Mathematics January 2013

Hsinchu, Taiwan, Republic of China

以半顯隱 Runge-Kutta 及擬譜法求解薛丁格方程

學生：林沛沅

共同

指導教授：賴明治

共同指導教授：鄧君豪

國立交通大學應用數學系數學建模與科學計算碩士班

摘

要

本論文使用半顯隱 Runge-Kutta 和擬譜法建立計算格式以求解薛丁格方

程。利用補償法將邊界條件加入格式中，藉由離散的能量估計，訂立適當的

懲罰參數。我們應用 Legendre-Gauss-Lobatto 及 Chebyshev-Gauss-Lobatto 這

兩組不同的網格點進行計算，並從幾個數值實驗來驗證此格式。

Implicit-Explicit Runge-Kutta and pseudospectral methods

for Schrodinger equation

Student：Pei-Yuan Lin

Advisor：Ming-Chih Lai

Co-Advisor：Chun-Hao Teng

Institute of Mathematic Modeling and Scientific Computing

Department of Applied Mathematics

National Chiao Tung University

Abstract

In this paper, we present a scheme for solving the Schrödinger equation based on

Implicit-Explicit Runge-Kutta and pseudospectral method. The boundary conditions

are imposed to the scheme through the penalty methodology. By conducting the

energy estimate, we determine the values of penalty parameters. We apply

Legendre-Gauss-Lobatto and Chebyshev-Gauss-Lobatto grid points for numerical

computations. Several numerical experiments are shown to validate the scheme.

誌

謝

目

錄

中文提要

………

i

英文提要

………

ii

誌謝

………

iii

目錄

………

iv

一、

Introduction………

1

二、

Formulation………

2

2.1

Model problem and energy estimate ………

2

2.2

Basic concepts of the pseudospectral method………

3

2.3

Semi-discrete schemes………

6

2.3.1 One-dimensional problem ………

6

2.3.2 Two-dimensional problem ………

11

2.4

Time integration ………

14

三、

Numerical results………

16

3.1

One-dimensional problem………

16

3.2

Two-dimensional problem ………

21

四、

Concluding remarks………

25

參考文獻

………

25

1

Introduction

The Sch¨ordinger equation is a partial differential equation (PDE) which describes

the behavior of a particle in quantum mechanics, formulated by the Austrian physicist

Erwin Sch¨ordinger [1]. It predicts the probability of observing the particle in a particular

position, and the equation is applied wildly in the fields of nuclear physics and quantum chemistry.

Due to the repaid development of computers, people have started to solve PDEs by numerical computations. To solve a problem numerically, it is necessary to construct a computational scheme. Therefore, how to obtain an efficient scheme and whether the numerical solution from the scheme converges to the exact solution become important.

The finite difference method is often used to discretize the Sch¨ordinger equation [6, 11].

However, it requires many grid points to obtain an accurate results in heavy computations. To reduce the computational loading, we introduce the pseudospectral methods [5] to solve problems.

The Lax-Richtmyer theorem [4] provides us a simple way to examine the convergence of a scheme for a linear problem. It states that a consistent scheme for a well-posed linear initial value problem is convergent if and only if it is stable. Therefore, we can ensure the convergence of a scheme by examining the consistency and stability of the scheme. A procedure proposed by von Neumann is commonly used to check the stability for partial differential equations. But in this paper, we establish the stability of proposed schemes by the energy method.

In this study, we present a pseudospectral scheme for the Sch¨ordinger equation defined

on the square domain subject to different types of boundary conditions. The Legendre-Gauss-Lobatto and Chebyshev-Legendre-Gauss-Lobatto grid points are introduced to discrete the space. The boundary conditions are imposed to the scheme through the penalty method-ology [10]. We pay attention to the stability of the scheme by conducting the energy

estimate. Through the discrete energy estimate, we determine the values of penalty

parameters to ensure the stable computations. For time discretization, we use the Crank-Nicolson and implicit-explicit Runge-Kutta methods. Because these methods are implicit,

it is necessary to invert matrices. Here, we adopt the eigen-decomposition approach [9] to conduct the matrix inversion.

This thesis is organized as follows. Section 2 states the initial boundary value problem and examines the well-posedness by conducting an energy estimate. In section 3, the concepts of pseudospectral methods are introduced. Then we propose the pseudospectral scheme and analyze the stability of the scheme. Section 4 presents the numerical results with several experiments. The concluding remarks are drawn in Section 5.

2

Formulation

2.1

Model problem and energy estimate

We consider the space domain Ω ⊂ R2 _{and denote the space and time coordinates by}

x = (x, y) and t, respectively. Let u = u(x, t) be a complex-valued function satisfying the initial boundary value problem (IBVP):

i∂u(x, t)

∂t = −ρ∇

2_{u(x, t), x ∈ Ω, t ≥ 0, (1a)}

u(x, 0) = f (x), x ∈ Ω, (1b)

Bu(x, t) = α(x)u(x, t) + β(x)∂u(x, t)

∂n = g(t), x ∈ ∂Ω, t ≥ 0, (1c)

where ρ is a positive constant, ∇2 is the Laplace operator and i =√−1 is the imaginary

unit, f is the initial data of u and Bu = g is the boundary condition imposed at the boundary domain ∂Ω, and B is the boundary operator parameterized by non-negative

functions α(x) and β(x) which satisfy the constrains α2_{(x) + β}2_{(x) 6= 0 on ∂Ω.}

We consider the homogeneous boundary conditions and assume that there is an unique

solution to the IBVP. Multiplying −iρ−1u∗ and iρ−1u to Eq. (1a) and its complex

conju-gate, respectively, and summing the resultants, we have 1 ρ  u∗∂u ∂t + u ∂u∗ ∂t  = iu∗∇2_{u − iu∇}2_u∗ , (2)

with the symbol ∗ denoting the complex conjugate. Integrating Eq. (2) over the domain

Ω, it becomes 1 ρ Z Ω  u∗∂u ∂t + u ∂u∗ ∂t  dx = i Z Ω u∗∇2_{u
dx − i} Z Ω u∇2u∗
dx. (3)

The left-hand side of Eq. (3) can be simplified as 1 ρ Z Ω  u∗∂u ∂t + u ∂u∗ ∂t  dx = 1 ρ d dt Z Ω |u|2dx.

Invoking the divergence theorem, the terms on the right-hand side of Eq. (3) become

i Z Ω u∗∇2_{u
dx = −} Z Ω (∇u∗· ∇u) dx + i I ∂Ω u∗(n · ∇u)dx, −i Z Ω u∇2u∗
dx = Z Ω (∇u · ∇u∗) dx − i I ∂Ω u(n · ∇u∗)dx,

where H () · dx denotes the surface integration. Applying the boundary condition, we

obtain the energy rate equation 1 ρ d dt Z Ω |u|2_{dx = −2} I ∂Ω |u|2_Im α β  dx = 0. (4)

It leads to an energy bound for u as Z Ω |u(x, t)|2_{dx =} Z Ω |f (x)|2_{dx, ∀t > 0}

Thus, this problem is well-posed.

2.2

Basic concepts of the pseudospectral method

Let N be a positive integer. The Legendre-Gauss-Lobatto (LGL) grid points xj are

the roots of the polynomial (1 − x2_)P0

N(x), where the prime denotes the differentiation

and PN(x) is the N -th degree Legendre polynomial defined by

PN(x) = 1 2N_{N !} dN dxN(x 2_{− 1)}N_{. (5)}

These points are arranged in ascending order in the interval I = [−1, 1]. In addition to the LGL grid points, we introduce another set of collocation points, the

Chebyshev-Gauss-Lobatto (CGL) grid points. The CGL points are the zeros of the polynomial (1−x2_)T0

N(x)

with TN(x) being the N -th degree Chebyshev polynomial

and the CGL points can be defined explicitly as xj = − cos  jπ N  , j = 0, 1, . . . , N. (7)

Pseudospectral methods are commonly based on interpolations at the LGL or CGL points. Based on a set of collocation points, we approximate a function f (x) defined on I as f (x) ≈ fN(x) = N X j=0 lj(x)f (xj), (8)

with lj(x) being the Lagrange basis polynomials given as

lj(x) = Y 0≤m≤N m6=j x − xm xj− xm , j = 0, 1, . . . , N.

Then the p-th derivative of f (x) can be also approximated as

f(p)(x) ≈ f_N(p)(x) = N X j=0 dplj(x) dxp f (xj).

Through a matrix-vector multiplication, the numerical derivatives at the grid points can be evaluated as

f(p) = Dpf ,

where f and f(p) _{are vectors given by}

f = [fN(x0), fN(x1), · · · , fN(xN)]T , f(p) =

h

f_N(p)(x0), f_N(p)(x1), · · · , f_N(p)(xN)

iT

, with the superscript T being the vector transpose, and D is called the differential matrix

with the entries Djk = lk0(xj). Notice that the differential matrix D varies when we use

different types of collocation points. The Legendre differentiation matrix DL has entries

DL_ji =    (−1)i+j 2  sin (i − j)π N + 1 −1 , i 6= j, 0, i = j.

And the entries of the Chebyshev differentiation matrix DC _are

DC_ij =                        ci cj (−1)i+j xi − xj , i 6= j −xj 2(1 − x2 j) , 2 ≤ i = j ≤ N, 2n2_{+ 1} 6 , i = j = 1, −2n 2_{+ 1} 6 , i = j = N + 1,

where

ci =

 2, i = 1 or N + 1, 1, 2 ≤ i = j ≤ N.

Associated with the LGL points especially, we have the quadrature formula

N X j=0 ωjf (xj) = Z 1 −1 f (x)dx, (9)

with f (x) being a polynomial of degree at most 2N − 1, and the quadrature weights ωj

are given by ωj =      − 2 N + 1[PN(xj)PN −1(xj)] −1 , j = 1, 2, . . . , N − 1, 2 N (N + 1), j = 0, N.

For further use, we have the following rules based on Eq. (9). Let u(x) and v(x) be polynomials of degree at most N . We have

N

X

j=0

ωju(xj)v0(xj) = u(xN)v(xN) − u(x0)v(x0) −

N X j=0 ωju0(xj)v(xj), (10) N X j=0 ωju(xj) (v(x)lj(x)) 0  xk = u(xN)v(xN)lN(xk) − u(x0)v(x0)l0(xk) − ωku0(xk)v(xk). (11) The above concepts can be extended for problems defined on two dimensional space.

Given {xj}Mj=0and {yk}Nk=0, the sets of grid points on [-1,1] along the x-axis and the y-axis,

respectively. Define the two-dimensional collection points (xj, yk). Based on these points,

we construct the two-dimensional Lagrange basis polynomials as

Ljk(x, y) = lxj(x)l y k(y), where lx j(x), l y

k(y) are the one-dimensional Lagrange interpolation polynomials based on

{xj}Mj=0 and {yk}Nk=0, respectively. Then we approximate a function f (x, y) defined on

I2 _{= [−1, 1]}2 _as f (x, y) ≈ fM N(x, y) = M X j=0 N X k=0 Ljk(x, y)f (xj, yk).

The partial derivatives of f at the grid point (xj, yk) are approximated as follows, ∂f (xj, yk) ∂x ≈ ∂fM N(xj, yk) ∂x = M X j0₌₀ N X k0₌₀ ∂Lj0_k0(x_j, y_k) ∂x f (xj0, yk0), ∂f (xj, yk) ∂y ≈ ∂fM N(xj, yk) ∂y = M X j0₌₀ N X k0₌₀ ∂Lj0_k0(x_j, y_k) ∂y f (xj0, yk0).

The numerical partial derivatives can also be calculated through a matrix-vector multi-plication as

Fx = DxF , Fy = F DyT,

where F is an (M + 1) × (N + 1) matrix with entries being Fjk = f (xj, yk), Dx and Dy

are the (M + 1) × (M + 1) and (N + 1) × (N + 1) differentiation matrices in the x- and the

y-directions, respectively. Fx and Fy are the matrices whose elements are the numerical

partial derivatives ∂fM N(xj, yk)/∂x and ∂fM N(xj, yk)/∂y, respectively.

Let u(x, y) and v(x, y) both be polynomials of degree at most M and N in x and y,

respectively. Denote ujk = u(xj, yk) and vjk = v(xj, yk). We have

M X j=0 N X k=0  ωu∂v ∂x     jk = N X k=0 ωy_k[(uv)|M k − (uv)|0k] − M X j=0 N X k=0  ω∂u ∂xv     jk, (12a) M X j=0 N X k=0  ωu∂v ∂y     jk = M X j=0 ωx_j[(uv)|jN − (uv)|j0] − M X j=0 N X k=0  ω∂u ∂yv     jk . (12b)

2.3

Semi-discrete schemes

Consider the problem on the interval I = [−1, 1].

i∂u ∂t = −ρ ∂2_u ∂x2, x ∈ I, t ≥ 0, (13a) u(x, 0) = f (x), x ∈ I, (13b) B−u(−1, t) = α−u(−1, t) − β− ∂u(−1, t) ∂x = g−(t), t ≥ 0, (13c) B+u(+1, t) = α+u(+1, t) + β+ ∂u(+1, t) ∂x = g+(t), t ≥ 0. (13d)

To solve the problem numerically, we collocate N + 1 LGL points xj for j = 0, 1, . . . , N

the form v(x, t) = N X j=0 lj(x)vj(t),

satisfying the collocation equations

i∂v(xj, t) ∂t = −ρ ∂F (xj, t) ∂x , j = 0, 1, . . . , N, (14a) v(xj, 0) = f (xj), j = 0, 1, . . . , N, (14b) where F (x, t) = ∂v(x, t) ∂x + τ−S−(x)(B−v0− g−(t)) − τ+S−(x)(B+vN − g+(t)), (14c) with S−(x) = (1 − x)P_N0 (x) 2P_N0 (−1) , S+(x) = (1 + x)P_N0 (x) 2P_N0 (1) , B−v0 = α−v(−1, t) − β− ∂v(−1, t) ∂x , B+vN = α+v(+1, t) + β+ ∂v(+1, t) ∂x .

The symbols τ− and τ+ are called the penalty parameters, and their values will be

de-termined later such that the scheme is stable. It is also noticed that S−(x) and S+(x)

coincide with the l0(x) and lN(x), respectively. The purpose of introducing S−(x) and

S+(x) is to avoid confusion when we use CGL grid points which will be shown later.

For stability analysis, we consider the homogeneous boundary conditions, namely,

g±(t) = 0. Consider Eq. (14a) and its complex conjugate version:

i∂v ∂t    j = −ρ ∂ ∂x  ∂v ∂x    j + τ−S−(xj)(B−v0) − τ+S+(xj)(B+vN)  , j = 0, 1, . . . , N, (15a) −i∂v ∗ ∂t    j = −ρ ∂ ∂x  ∂v∗ ∂x    j+ τ ∗ −S−(xj)(B−v∗0) − τ ∗ +S+(xj)(B+v∗N)  , j = 0, 1, . . . , N. (15b)

Multiplying Eq. (15a) and Eq. (15b) by −iρ−1ωjvj∗ and iρ

−1_ω

sum-ming the resultants over the index j = 0 to N , we obtain 1 ρ N X j=0  ωv∗∂v ∂t    j =i N X j=0  ωv∗∂ 2_q ∂x2    j + iτ− N X j=0 ωjvj∗ ∂ ∂x(S−(xj)B−v0) − iτ+ N X j=0 ωjvj∗ ∂ ∂x(S+(xj)B+vN) (16a) 1 ρ N X j=0  ωv∂v ∗ ∂t     j = − i N X j=0  ωv∂ 2_v∗ ∂x2     j − iτ ∗ − N X j=0 ωjvj ∂ ∂x(S−(xj)B−v ∗ 0) + iτ₊∗ N X j=0 ωjvj ∂ ∂x(S+(xj)B+v ∗ N) (16b)

Applying Eq. (10) to the first summation term on the right-hand side of Eq. (16a), we obtain i N X j=0  ωv∗∂ 2_v ∂x2    j = iv∗∂v ∂x    N − iv∗∂v ∂x    0 − i N X j=0  ω∂v ∗ ∂x ∂v ∂x    j . (17)

The other terms on the right-hand side of Eq. (16a) can be evaluated by invoking Eq. (11), and the results are given as follows,

iτ− N X j=0 ωjv∗j ∂ ∂x(S−(xj)B−v0) = −iτ−(v ∗ + ω0 ∂v∗ ∂x)(α−v − β− ∂v ∂x)    0 (18) −iτ+ N X j=0 ωjv∗j ∂ ∂x(S+(xj)B+vN) = −iτ+(v ∗_{− ω} N ∂v∗ ∂x)(α+v + β+ ∂v ∂x)    N (19)

With similar arguments, for the terms on the right-hand side of Eq. (16b), we have the results −i N X j=0  ωv∂ 2_v∗ ∂x2     j = −iv ∂v∗ ∂x    N + iv ∂v∗ ∂x    0+ i N X j=0  ω∂v ∂x ∂v∗ ∂x     j, (20) −iτ₋∗ N X j=0 ωjvj ∂ ∂x(S−(xj)B−v ∗ 0) = iτ ∗ −(v + ω0 ∂v ∂x)(α−v ∗_{− β} − ∂v∗ ∂x)    0 , (21) iτ₊∗ N X j=0 ωjvj ∂ ∂x(S+(xj)B+v ∗ N) = iτ ∗ +(v − ωN ∂v ∂x)(α+v ∗ + β+ ∂v∗ ∂x)    N . (22)

Define a vector function r and a matrix function A as

r(r1, r2) = [r1, r2]T, A(α, β, τ, ω) =  −τ α + τ∗_{α 1 − τ β − ωτ}∗_α −1 + τ∗_{β + ωτ α ωτ β − ωτ}∗_β  (23)

Adding Eq. (16a) to Eq. (16b) and substituting Eqs. (17)-(22) into the resultants, we have an energy rate equation

1 ρ d dt N X j=0 ωj|vj|2 = ir∗−A−r−+ ir∗₊A+r+,

where r−, r+ are vectors given as

r−(t) = r(v(x0, t), −∂v(x0, t)/∂x), r+(t) = r(v(xN, t), ∂v(xN, t)/∂x), (24)

and

A± = A(α±, β±, τ±, ¯ω), ω =¯

2

N (N + 1). (25)

Taking the values of τ− and τ+ as

τ− = 1 ¯ ωα−+ β− , τ+= 1 ¯ ωα++ β+ , (26)

we have A− and A+ being zero matrices. This leads to

1 ρ d dt N X j=0 ωj|vj|2 = 0,

implying the stability of the scheme.

The semi-discrete scheme Eq. (14) can be represented in the form of matrices and vectors. We introduce the matrices

E− = eeT0, E+ = eeTN, (27)

with e0, eN, and e being the vectors of length N + 1, given by

e0 = [1, 0, . . . , 0]T, eN = [0, . . . , 0, 1]T, e = [1, . . . , 1]T.

Define the solution vector

v(t) = [v(x0, t), v(x1, t), . . . , v(xN, t)]T.

Then the semi-discrete scheme Eq. (14) can be expressed as ∂v

∂t = Lv + g(t), (28a)

where L is the matrix operator and g(t) collects the terms caused by the time-dependent boundary conditions, given as

L = iρ(D(D + τ−S−E−(α−I0− β−I0D) − τ+S+E+(α+IN + β+IND))), (29)

g(t) = iρ(−τ−g−(t)DS−E−e0+ τ+g+(t)DS+E+eN). (30)

f = [f (x0), f (x1), · · · , f (xN)]T is a vector concerned with initial data. In the above

expression, I0, IN, S−, and S+ are (N + 1) × (N + 1) matrices defined by

S± = diag(S±(x0), S±(x1), . . . , S±(xN)), (31)

I0 = diag(1, 0, . . . , 0), IN = diag(0, . . . , 0, 1), (32)

and D is the differentiation matrix corresponding to the grid points.

We have constructed a stable scheme based on the LGL grid points. The constructed scheme can be slightly modified for computations based on CGL grid points, also known as the Chebyshev-Legendre method [2]. Here we briefly summarized the modification. We seek a numerical solution of the form

v(x, t) = N X j=0 l_jc(x)vj (33) where lc

j(x) are the Lagrange basis polynomials based on the CGL grid points and vj are

the field values collocated at the CGL grid points. We require the solution satisfy the scheme Eq. (14) at the CGL grid points,

i∂v(x c j, t) ∂t = −ρ ∂F (xc j, t) ∂x , j = 0, 1, . . . , N, (34a) v(xc_j, 0) = f (xc_j), j = 0, 1, . . . , N, (34b) where F (x, t) = ∂v(x, t) ∂x + τ−S−(x)(B−v0− g−(t)) − τ+S−(x)(B+vN − g+(t)), (34c) with S−(x) = (1 − x)P_N0 (x) 2P_N0 (−1) , S+(x) = (1 + x)P_N0 (x) 2P_N0 (1) ,

Notice that numerical solution v(x, t) in fact satisfy the partial differential equation:

i∂v(x, t)

∂t = −ρ

∂F (x, t) ∂x

Hence, we can establish the stability of the scheme based on the Legendre integration quadrature rule as shown before. We conclude the scheme is stable provided that the penalty parameters are given in Eq. (26).

2.3.2 Two-dimensional problem

Let us consider the IBVP problem defined on I2 _{= [−1, 1] × [−1, 1]:}

i∂u(x, y, t)

∂t = −ρ∇

2_{u(x, y, t), (x, y) ∈ I}2_{, t ≥ 0, (35a)}

u(x, y, 0) = f (x, y), (x, y) ∈ I2, (35b)

B(a)_{u(−1, y, t) = g}

−(y, t), B(a) = α(a)− β(a)

∂ ∂x, y ∈ I, t ≥ 0, (35c) B(b)u(+1, y, t) = g+(y, t), B(b) = α(b)+ β(b) ∂ ∂x, y ∈ I, t ≥ 0, (35d) B(c)u(x, −1, t) = h−(x, t), B(c) = α(c)− β(c) ∂ ∂y, x ∈ I, t ≥ 0, (35e) B(d)_{u(x, +1, t) = h} +(x, t), B(d) = α(d)+ β(d) ∂ ∂y, x ∈ I, t ≥ 0. (35f)

For γ ∈ {a, b, c, d}, B(γ) are the boundary operators defined on the edges of the domain.

Each boundary operator is parameterized by two non-negative constants α(γ) and β(γ)

satisfying the constrain (α(γ)₎2_{+ (β}(γ)₎2 _{6= 0.}

Introduce the two-dimensional LGL grid points (xj, yk) based on the sets of LGL

points {xj}Mj=0 and {yk}Nk=0. We seek a numerical solution of the form

v(x, y, t) = M X j=0 N X k=0 Ljk(x, y)v(xj, yk, t).

satisfying the scheme:

i∂v(xj, yk, t) ∂t = −ρ ∂Fx(xj, yk, t) ∂x − ρ ∂Fy(xj, yk, t) ∂y , 0 ≤ j ≤ M, 0 ≤ k ≤ N, (36a) v(xj, yk, 0) = f (xj, yk), 0 ≤ j ≤ M, 0 ≤ k ≤ N, (36b)

with Fx(x, y, t) = ∂v(x, y, t) ∂x + N X k0₌₀

L0k0(x, y)τ(a)(B(a)v_0k0 − g₋(y, t))

− N X k0₌₀ LM k0(x, y)τ(b)(B(b)v_{M k}0 − g₊(y, t)), (36c) Fy(x, y, t) = ∂v(x, y, t) ∂y + M X j0₌₀ Lj0₀(x, y)τ(c)(B(c)v_j0₀− h₋(x, t)) − M X j0₌₀ Lj0_N(x, y)τ(d)(B(d)v_j0_N − h₊(x, t)). (36d)

τ(a), τ(b), τ(c), and τ(d) are the penalty parameters associated with the edges and their

values will be determined through conducting an energy estimate for the scheme.

To conduct the stability analysis, we consider the problem subject to the homogeneous

boundary conditions, that is, g±(y, t) = 0 and h±(x, t) = 0. Define ωjk = ωxjω

y

k with ωjx

and ω_ky being the quadrature weights associated with LGL points xj and yk, respectively.

Multiplying −iρ−1ωjkvjk∗ and iρ

−1_ω

jkvjk to Eq. (36a) and its complex conjugate,

respec-tively, summing the resultants over the indices j = 0 to M and k = 0 to N , and adding them together, we have the energy rate equation

1 ρ d dt M X j=0 N X k=0 ω|v|2
  jk = M X j=0 N X k=0  iωv∗∂Fx ∂x     jk + M X j=0 N X k=0  iωv∗∂Fy ∂y     jk + M X j=0 N X k=0  −iωv∂F ∗ x ∂x     jk+ M X j=0 N X k=0  −iωv∂F ∗ y ∂y     jk. (37)

Invoking Eqs. (12a)-(12b), we can evaluate the summation terms on the right-hand side of Eq. (37) and have

1 ρ d dt N X j=0 M X k=0 ω|v|2
   jk = i N X k=0 ωy_kr_k(a) ∗ A(a)r_k(a)+ i N X k=0 ωy_kr_k(b) ∗ A(b)r_k(b) + i M X j=0 ω_jxr_j(c) ∗ A(c)r(c)_j + i M X j=0 ω_jxr(d)_j  ∗ A(d)r_j(d)

where through Eq. (23)

r(a)_k = r(v0k, −∂v0k/∂x), r (b) k = r(vM k, ∂vM k/∂x), r_j(c) = r(vj0, −∂vj0/∂y), r (d) j = r(vjN, ∂vjN/∂y),

and

A(γ)= A(α(γ), β(γ), τ(γ), ω(γ)), for γ = a, b, c, and d.

ω(γ)=

( ₂

M (M +1) if γ = a, b 2

N (N +1) if γ = c, d

To ensure the stability of the scheme, we request

τ(γ) = 1

ω(γ)_α(γ)_{+ β}(γ),

so that A(γ) _{are all zero matrices, and we have}

1 ρ d dt M X j=0 N X k=0 ω|v|2
  jk = 0.

The semi-discrete scheme has a matrix-vector representation. Let v(t) be a solution

matrix with the entries vjk = v(xj, yk, t). Dx and Dy are the differentiation matrices

respect to x− and y− directions. The scheme (36) can be represented as ∂v

∂t = Lv + vR + G(t), (38a)

v(0) = f , (38b)

where L ∈ C(M +1)×(M +1) and R ∈ C(N +1)×(N +1) are matrix operators

L = iρ(Dx(Dx+ τ(a)S−xE−x(α(a)I0x− β

(a)_Ix 0Dx) − τ(b)S+xE x +(α (b)_Ix M + β (b)_Ix MDx))), R = iρ((D_yT + τ(c)(α(c)I₀y − β(c)_DT yI y 0)(S y −E y −)T − τ(d)(α(d)I y N + β (d)_DT yI y N)(S y +E y +) T_)DT y),

and G(t) consists of the terms from discrete boundary conditions

G(t) = iρ(−τ(a)DxS−xE−xI0x(exg−(t)) + τ(b)DxS+xE x +I x M(exg+(t))) + iρ(−τ(c)(h−(t)eTy)(DyS y −E y −I y 0) T + τ(d)(h+(t)eTy)(DyS y +E y +I y N) T ),

with ex and ey being the vectors of length N + 1 and M + 1 that all the components are

equal to 1, and

g±(t) = [g±(y0, t), g±(y1, t), . . . , g±(yN, t)],

h±(t) = [h±(x0, t), h±(x1, t), . . . , h±(xM, t)]T.

For ν ∈ {x, y}, the matrices Sν

2.4

Time integration

To march the numerical solution in time, we adopt the Crank-Nicolson (CN) method [7] and Implicit-explicit Runge-Kutta (IMEX-RK) method [3, 8]. Denote the time step

by ∆t and the n-th time level by tn _{= n∆t. Let u}n _{be the numerical solution at time t}n

for the following differential equation du

dt = F (t, u(t)), (39a)

u(0) = u0. (39b)

The Crank-Nicolson algorithm solves Eqs. (39) as

un+1_{− u}n ∆t = 1 2(F (tn+1, u n+1 ) + F (tn, un)). (40)

Then we have the fully-discrete version for Eqs. (14) as

iv n+1 j − vnj ∆t = − ρ 2  ∂F (xj, tn+1) ∂x + ∂F (xj, tn) ∂x  , j = 0, 1, . . . , N, (41a) v_j0 = f (xj), j = 0, 1, . . . , N. (41b)

To analyze the stability of Eqs. (41), we consider the homogeneous boundary conditions.

For convenience, we denote (v_jn+1+ vn

j)/2 by v n+1/2 j . Multiplying −iρ −1_ω j  vn+1/2_j  ∗ and iρ−1ωjv n+1/2

j to Eq. (41a) and its complex conjugate, respectively, summing the resultants,

and following a similar approach in semi-discrete scheme, we obtain 1 2ρ∆t N X j=0 ωj(|vn+1j | 2_{− |v}n j| 2_{) = ir}∗ −A−r−+ ir∗₊A+r+,

where A± are given in Eq. (25) and

r− = r(v0n+1/2, −∂v n+1/2 0 /∂x), r+ = r(v n+1/2 N , ∂v n+1/2 N /∂x).

For τ± given in Eq. (26), A± are zero matrices. Thus, we have

N X j=0 ωj|vjn+1| 2 ₌ N X j=0 ωj|vjn|2 = · · · = N X j=0 ωj|f (xj)|2,

In addition, we have the fully-discrete scheme for Eq. (38a) with vnbeing the numerical

solution matrix at time tn

vn+1− vn ∆t = 1 2(Lv n+1_{+ v}n+1_{R + G(t} n+1) + Lvn+ vnR + G(tn)). (42)

The Eq. (42) can be rewritten as

Avn+1+ vn+1B = F , (43) where A = 1 2I (M )₋ ∆t 2 L, B = 1 2I (N )₋ ∆t 2 R, F = vn₊ ∆t 2 (Lv n_{+ v}n_{R + +G(t}n_{) + G(t} n+1)) .

with I(M ) _{∈ R}(M +1)×(M +1) _{and I}(N ) _{∈ R}(N +1)×(N +1) _{being identity matrices. Assuming}

that A and B are diagonalizable, we can directly solve vn+1 _{by the eigenvector}

decom-position method [9].

For Eqs. (39), the IMEX-RK method splits the flux function F into two parts as

F (t, u(t)) = F[im](t, u(t)) + F[ex](t, u(t)),

where F[im]_{and F}[ex] _{are flux functions to be treated in implicit and explicit ways,}

respec-tively. For a s-stage IMEX-RK method, we solve Eqs. (39) numerically by the following steps: u(i) = un+ ∆t s X j=1 a[im]_ij F[im](tn+ cj∆t, u(j)) + ∆t s X j=1 a[ex]_ij F[ex](tn+ cj∆t, u(j)), 1 ≤ i ≤ s un+1= un+ ∆t s X i=1

b[im]_i F[im](tn+ ci∆t, u(i)) + ∆t

s

X

i=1

b[ex]_i F[ex](tn+ ci∆t, u(i)).

The coefficients of a[im]_ij , a[ex]_ij , b[im]_i , b[ex]_i , cj can be found in [3]. We consider Lv+vR+G(t)

in Eq. (38a) to be the implicit part, and the explicit part is regarded as zero. Then we have a fully-discrete scheme with IMEX-RK

v(i) = vn+ ∆t s X j=1 a[im]_ij (Lv(j)+ v(j)R + G(tn+ cj∆t)), 1 ≤ i ≤ s (44a) vn+1= vn+ ∆t s X

We can rewrite Eq. (44a) as

Av(i)_{+ v}(i)_{B = F}(i)_{, 1 ≤ i ≤ s, (45)}

where A = 1 2I (M )₋∆t 2 L, B = 1 2I (N )₋ ∆t 2 R, F(i) ₌        vn_{+ ∆t a}[im] 11 G(tn+ c1∆t), if i = 1, vn_{+ ∆t} i−1 X j=1 a[im]_ij F(j)_{+ ∆t a}[im] ii G(t n_{+ c} i∆t), if i 6= 1.

Thus, we solve Eqs. (45) directly to obtain v(i)_{, and then v}n+1 _{can be computed by}

Eqs. (44b).

3

Numerical results

In this section, we illustrate the proposed methods by several examples. Denote ∆x =

2N−1 which is the mean distance for a set of N + 1 grid points. The time step ∆t is

computed adaptively as

∆t = CF L∆x,

where CF L is referred to the Courant-Friedrichs-Lewy number. The convergence order is calculated as

q = log(e(N1)/e(N2))

log(N2/N1)

,

where e(N ) = ku − uNk∞is the maximum error with u and uN being the analytic solution

and the numerical solution corresponding to polynomial degree N , respectively.

3.1

One-dimensional problem

Example 1. Let I = [−1, 1]. Consider q(x, t) = ei(x−t) _{satisfying the problem:}

i∂q

∂t = −

∂2q

∂x2, x ∈ I, t ≥ 0,

q(x, 0) = eix, x ∈ I,

Table 1-3 shows the results subject to different types of boundary conditions imposed at x = ±1 applying Crank-Nicolson method. For each terminal time T , the error decreases as N increases and the rate of convergence is approximately a two. Compare the results from different collocation points, there is little difference between them and the rate of convergence is similar. Thus, we know that the scheme can achieve a convergent result.

Table 1: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by Crank-Nicolson method. α± = 1, β±= 0, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 7.10e-05 - 7.10e-05 - 7.22e-05 - 7.22e-05

-12 3.12e-05 2.03 3.12e-05 2.03 3.21e-05 2.00 3.21e-05 2.00

16 1.76e-05 1.99 1.76e-05 1.99 1.81e-05 1.98 1.81e-05 1.98

20 1.13e-05 1.98 1.13e-05 1.98 1.18e-05 1.94 1.18e-05 1.94

Table 2: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by Crank-Nicolson method. α± = 0, β±= 1, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.05e-04 - 1.05e-04 - 6.53e-05 - 6.53e-05

-12 4.64e-05 2.02 4.64e-05 2.02 2.86e-05 2.03 2.86e-05 2.03

16 2.61e-05 2.00 2.61e-05 2.00 1.64e-05 1.94 1.64e-05 1.94

20 1.67e-05 2.00 1.67e-05 2.00 1.06e-05 1.97 1.06e-05 1.97

The IMEX-RK method can be applied to solve Example 1, and the results are shown in Table 4-6. We see the convergence order is approximately a three in each case. So the scheme produced by IMEX-RK is still convergent whether we use CGL or LGL points. Fig. 1 reveals that the discrete energy is conserved, where ∆E is defined as

∆E = N X j=1 wj|vj|2− N X j=1 wj|f (xj)|2.

Table 3: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by Crank-Nicolson method. α± = 1.5, β± = 0.5, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 2.24e-04 - 2.24e-04 - 2.05e-04 - 2.05e-04

-12 9.97e-05 2.00 9.97e-05 2.00 9.07e-05 2.01 9.07e-05 2.01

16 5.60e-05 2.00 5.60e-05 2.00 5.10e-05 2.00 5.10e-05 2.00

20 3.59e-05 2.00 3.59e-05 2.00 3.26e-05 2.01 3.26e-05 2.01

Table 4: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. α±= 1, β±= 0, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.43e-06 - 1.43e-06 - 1.33e-06 - 1.34e-06

-12 4.05e-07 3.10 4.07e-07 3.10 4.13e-07 2.88 4.14e-07 2.91

16 1.72e-07 2.98 1.72e-07 2.99 1.75e-07 2.98 1.75e-07 2.99

20 8.84e-08 2.98 8.84e-08 2.98 8.86e-08 3.05 8.86e-08 3.05

Table 5: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. α±= 0, β±= 1, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.81e-06 - 1.81e-06 - 1.16e-06 - 1.16e-06

-12 2.68e-07 4.71 2.68e-07 4.71 2.79e-07 3.50 2.79e-07 3.50

16 1.14e-07 2.98 1.14e-07 2.98 1.20e-07 2.95 1.20e-07 2.95

20 5.87e-08 2.97 5.87e-08 2.97 6.18e-08 2.96 6.18e-08 2.96

Table 6: The rate of convergence for Example 1 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. α±= 1.5, β± = 0.5, CF L = 0.1

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 2.42e-06 - 2.45e-06 - 2.17e-06 - 2.17e-06

-12 7.21e-07 2.99 7.21e-07 3.02 6.66e-07 2.92 6.66e-07 2.92

16 3.04e-07 3.00 3.04e-07 3.00 2.82e-07 2.99 2.82e-07 2.99

0 2 4 6 8 10 0 0.5 1 1.5 2 2.5 3 3.5x 10 −5 T ∆ E 0 2 4 6 8 10 −1.5 −1 −0.5 0 0.5 1 1.5x 10 −7 T ∆ E

Figure 1: The difference of discrete energy for Example 1 subject to Dirichlet boundary condition. Left: Time integration by Crank-Nicolson method. Left: Time integration by third-order IMEX-RK.

Example 2. Let I = [−1, 1]. Consider q(x, t) = sin(πx) sin(πt) satisfying the problem:

i∂q

∂t = −

∂2_q

∂x2 + iπ sin(πx) cos(πt) − π

2_{sin(πx) sin(πt), x ∈ I, t ≥ 0,}

q(x, 0) = 0, x ∈ I,

B±q(±1, t) = α±sin(±π) sin(πt) ± β±π cos(±π) sin(πt), t ≥ 0,

In this experiment, the partial differential equation contains a source term. To solve Example 2 with the IMEX-RK method, the source term is treated in the explicit flux. The results are shown in Table 7- 8. We observe an order reduction occurs when the boundary conditions are not Dirichlet-type. However, if we apply the Crank-Nicolson method to solve the same problem, the convergence order is approximately a two for every types of boundary conditions. Since the Dirichlet boundary condition in this problem is exactly zero, we guess that the order reduction occurs if we use IMEX-RK in time to solve a problem which contains a source term and nonzero boundary conditions.

Table 7: Convergence order for Example 2. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson at T = 1 with CGL points. CF L = 0.1

N

α±= 1, β± = 0 α±= 0, β± = 1

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.72e-04 - 8.14e-05 - 1.01e-03 - 5.98e-04

-12 5.38e-05 2.87 3.40e-05 2.16 4.02e-04 2.27 2.65e-04 2.01

16 2.24e-05 3.04 1.92e-05 1.98 1.94e-04 2.54 1.48e-04 2.02

Table 8: Convergence order for Example 2. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson at T = 1 with LGL points. CF L = 0.1

N

α±= 1, β± = 0 α±= 0, β± = 1

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.77e-04 - 7.85e-05 - 1.00e-03 - 6.08e-04

-12 5.37e-05 2.94 3.39e-05 2.07 4.02e-04 2.25 2.64e-04 2.05

16 2.26e-05 3.01 1.94e-05 1.95 1.94e-04 2.54 1.48e-04 2.02

20 1.15e-05 3.03 1.24e-05 1.99 1.11e-04 2.48 9.41e-05 2.03

Example 3. Let I = [−1, 1]. Consider q(x, t) = ei(x+t) _{satisfying the problem:}

i∂q ∂t = − ∂2_q ∂x2 − 2e i(x+t)_{, x ∈ I, t ≥ 0,} q(x, 0) = eix, x ∈ I,

B±q(±1, t) = α±ei(±1+t)± iβ±ei(±1+t), t ≥ 0,

We construct Example 3 whose boundary conditions are all nonhomogeneous. Table 9-10 show the results as we solve this problem with IMEX-RK and Crank-Nicolson method, respectively. The order reduction still occurs in the case of IMEX-RK. Meanwhile, the results produced by Crank-Nicolson method present a second-order convergence as we expect.

Table 9: Convergence order for Example 3. The numerical results are computed by third-order IMEX-RK and Crank-Nicolson with CGL points at T = 1. CF L = 0.1

N

α±= 1, β± = 0 α±= 0, β± = 1

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 3.25e-04 - 3.33e-05 - 3.27e-05 - 6.37e-05

-12 1.56e-04 1.81 1.48e-05 1.99 1.20e-05 2.46 2.83e-05 2.00

16 8.28e-05 2.20 8.48e-06 1.95 5.47e-06 2.74 1.59e-05 2.00

Table 10: Convergence order for Example 3. The numerical results are computed by third-order IMEX-RK and Crank-Nicolson with LGL points at T = 1. CF L = 0.1

N

α±= 1, β± = 0 α±= 0, β± = 1

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 3.25e-04 - 3.33e-05 - 3.27e-05 - 6.35e-05

-12 1.56e-04 1.81 1.48e-05 1.99 1.20e-05 2.46 2.83e-05 2.00

16 8.28e-05 2.20 8.49e-06 1.94 5.47e-06 2.74 1.59e-05 2.00

20 5.16e-05 2.12 5.46e-06 1.98 3.05e-06 2.62 1.02e-05 2.00

3.2

Two-dimensional problem

Example 4. Let I2 = [−1, 1]×[−1, 1]. Consider q(x, y, t) = e−2iπ2tcos(π(x+y)) satisfying

the problem: i∂q ∂t = − ∂2_q ∂x2 − ∂2_q ∂y2, (x, y) ∈ I 2_{, t ≥ 0,} q(x, y, 0) = cos(π(x + y)), (x, y) ∈ I2,

B(a)_{q(−1, y, t) = α}(a)_e−2iπ2_t

cos(π(y − 1)) + β(a)πe−2iπ2tsin(π(y − 1)), y ∈ I, t ≥ 0,

B(b)q(+1, y, t) = α(b)e−2iπ2tcos(π(y + 1)) − β(b)πe−2iπ2tsin(π(y + 1)), y ∈ I, t ≥ 0,

B(c)_{q(x, −1, t) = α}(c)_e−2iπ2_t

cos(π(x − 1)) + β(c)πe−2iπ2tsin(π(x − 1)), x ∈ I, t ≥ 0,

B(d)_{q(x, +1, t) = α}(d)_e−2iπ2_t

cos(π(x + 1)) − β(d)πe−2iπ2tsin(π(x + 1)), x ∈ I, t ≥ 0.

We consider Example 4 to test the scheme for two-dimension problems. For simpli-fication, the grid revolutions N and M are set to be equal. We apply Crank-Nicolson method in time to solve this problem, and show the convergence results in Table 11-13. We use different collocation points to compute the numerical solutions. For each terminal time T , The error decreases as N increases with a second-order convergence. The scheme is stable after a long-time computation (T = 100). Moveover, for a fixed N , the error grows linearly in time.

Table 14-16 present the results computed by IMEX-RK in time. For each terminal time T , the the rate of convergence decays rapidly in the beginning, and then goes down to third order. We see the error grows linearly for a fixed grid revolution N . It shows

Table 11: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) _{= 1, β}(γ) _{= 0.}

CF L = 0.01

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 3.58e-02 - 3.37e-02 - 3.50e-01 - 3.32e-01

-12 1.79e-02 1.71 1.79e-02 1.57 1.78e-01 1.67 1.77e-01 1.55

16 9.79e-03 2.10 9.92e-03 2.04 9.71e-02 2.10 9.84e-02 2.05

20 6.30e-03 1.97 6.24e-03 2.08 6.28e-02 1.96 6.21e-02 2.06

Table 12: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by third-order Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) ₌

0, β(γ)_{= 1. CF L = 0.01}

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 4.81e-02 - 4.84e-02 - 3.95e-01 - 3.95e-01

-12 1.80e-02 2.42 1.80e-02 2.43 1.78e-01 1.97 1.78e-01 1.97

16 1.02e-02 2.00 1.02e-02 2.00 1.00e-01 2.00 1.00e-01 2.00

20 6.50e-03 2.00 6.50e-03 2.00 6.41e-02 2.00 6.41e-02 2.00

Table 13: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by third-order Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) ₌

1.5, β(γ)_{= 0.5. CF L = 0.01}

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.79e-02 - 1.80e-02 - 1.68e-02 - 1.69e-02

-12 4.76e-04 8.95 4.92e-04 8.88 6.58e-04 7.99 6.75e-04 7.94

16 2.76e-04 1.88 2.80e-04 1.95 3.88e-04 1.84 3.91e-04 1.90

Table 14: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ) _{= 1, β}(γ) _{= 0.}

CF L = 0.01

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 4.58e-03 - 4.82e-03 - 2.26e-02 - 2.33e-02

-12 1.84e-04 7.93 1.83e-04 8.06 1.82e-03 6.21 1.82e-03 6.29

16 7.52e-05 3.11 7.61e-05 3.05 7.44e-04 3.11 7.55e-04 3.05

20 3.87e-05 2.98 3.83e-05 3.08 3.85e-04 2.96 3.81e-04 3.06

Table 15: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ) _{= 0, β}(γ) _{= 1.}

CF L = 0.01

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.05e-02 - 1.06e-02 - 9.26e-03 - 9.18e-03

-12 1.85e-04 9.95 1.84e-04 10.00 1.82e-03 4.02 1.82e-03 3.99

16 7.84e-05 2.99 7.84e-05 2.96 7.68e-04 3.00 7.68e-04 3.00

20 4.02e-05 3.00 4.02e-05 3.00 3.93e-04 3.00 3.93e-04 3.00

Table 16: The rate of convergence for Example 4 with CGL and LGL points at different terminal time. The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ)_{= 1.5, β}(γ)_{= 0.5.}

CF L = 0.01

N

T = 10 T = 100

CGL LGL CGL LGL

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.62e-02 - 1.62e-02 - 1.41e-02 - 1.41e-02

-12 1.21e-05 17.74 1.22e-05 17.74 8.97e-06 18.15 9.13e-06 18.11

16 2.33e-06 5.73 2.34e-06 5.73 3.25e-06 3.53 3.25e-06 3.59

Example 5. Let I2 = [−1, 1] × [−1, 1]. Consider q(x, y, t) = e−itsin(πx) sin(πy) satisfying the problem: i∂q ∂t = − ∂2_q ∂x2 − ∂2_q ∂y2 + (1 − 2π

2_)e−it_{sin(πx) sin(πy), (x, y) ∈ I}2_{, t ≥ 0,}

q(x, y, 0) = sin(πx) sin(πy), (x, y) ∈ I2,

B(a)_{q(−1, y, t) = α}(a)_e−it

sin(−π) sin(πy) − β(a)πe−itcos(−π) sin(πy), y ∈ I, t ≥ 0,

B(b)_{q(+1, y, t) = α}(b)_e−it_{sin(+π) sin(πy) + β}(b)_πe−it_{cos(+π) sin(πy), y ∈ I, t ≥ 0,}

B(c)_{q(x, −1, t) = α}(c)_e−it

sin(πx) sin(−π) − β(c)πe−itsin(πx) cos(−π), x ∈ I, t ≥ 0,

B(d)_{q(x, +1, t) = α}(d)_e−it

sin(πx) sin(+π) + β(d)πe−itsin(πx) cos(+π), x ∈ I, t ≥ 0.

We can use the scheme to solve Example 5 which contains a source term. Table 17-18 present the results subject to different types of boundary conditions. We see the results computed by the Crank-Nicolson method achieve a convergence of second order indeed. The order reduction still occurs while we use IMEX-RK to solve this problem with the nonhomogeneous boundary conditions.

It has been observed that if we use pseudospectral penalty method to construct a scheme, the order reduction may occurs for Runge-Kutta methods. The paper [2] provide us a procedure to deal with explicit Runge-Kutta. However, there is no way for implicit Runge-Kutta up to now.

Table 17: Convergence order for Example 5 at T = 1. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson with LGL points. CF L = 0.1

N

α(γ)_{= 1, β}(γ)_{= 0 α}(γ) _{= 0, β}(γ) _{= 1}

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.12e-03 - 1.07e-03 - 1.68e-03 - 3.20e-03

-12 1.69e-05 10.36 4.59e-07 19.13 2.51e-04 4.69 6.98e-06 15.11

16 6.96e-06 3.07 1.43e-07 4.05 1.24e-04 2.44 3.76e-06 2.15

Table 18: Convergence order for Example 5 at T = 1. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson with CGL points. CF L = 0.1

N

α(γ)= 1, β(γ)= 0 α(γ) = 0, β(γ) = 1

IMEX-RK CN IMEX-RK CN

e(N ) q e(N ) q e(N ) q e(N ) q

8 1.13e-03 - 1.27e-03 - 1.62e-03 - 3.17e-03

-12 1.69e-05 10.36 4.66e-07 19.52 2.55e-04 4.56 6.99e-06 15.08

16 6.86e-06 3.13 1.41e-07 4.15 1.25e-04 2.47 3.76e-06 2.15

20 3.54e-06 2.96 7.57e-08 2.79 6.99e-05 2.61 2.45e-06 1.92

4

Concluding remarks

In this study, we proposed a numerical scheme for solving the Sch¨ordinger equation

based on pseudospectral penalty method. For stable computations, we determine the penalty parameters subject to different types of boundary conditions by conducting the energy estimate. Although we establish the stability through LGL points, the scheme still works when we employ CGL points. We apply Crank-Nicolson method and

IMEX-RK for time-discretization. Several numerical experiments are shown to validate the

scheme, and we observe the expected convergence rate in most cases. An order reduction occurs only when we use IMEX-RK to solve the problems that contain a source term and nonhomogeneous boundary conditions.

The present method well solves the problems in a simple domain with scalar boundary parameters. In the future, we can conduct some numerical experiments for the problems whose boundary conditions are parameterized by functions. The method may be

gener-alized to solve the non-linear Sch¨ordinger equations defined on complicated domains.

References

[1] E. Schr¨odinger, An Undulatory Theory of the Mechanics of Atoms and Molecules.

Phys. Rev. Vol. 28, No. 6 (1926) pp. 1049-1070

[2] W.S. Don, D. Gottlieb, The Chebyshev-Legendre method: Implementing Legendre methods on Chebyshev points. SIAM J. Numer. Anal., 31 (1994), pp. 1519-1534.

[3] C.A. Kennedy, M.H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations. Appl. Numer. Math., 44 (2003), pp. 139-181.

[4] P. D. Lax, R. D. Richtmyer, Survey of the stability of linear finite difference equations. Communications on Pure and Applied Mathematics, 9 (1956) 267-293.

[5] J.S. Hesthaven, S. Gottlieb, D. Gottlieb, Spectral Methods for Time-Dependent Problems. Cambridge Universuty Press. Prentice Hall, (2007).

[6] Anna Nissen, Gunilla Kreiss, Margot Gerritsen, High Order Stable Finite Difference

Methods for the Schr¨odinger Equation. J. Sci. Comput. 55(1): 173-199 (2013)

[7] J. Crank, P. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Camb. Phil. Soc., 43, 1 (1947) 50-67.

[8] A. Kanevsky, M.H. Carpenter, D. Gottlieb, J.S. Hesthaven, Application of implic-itVexplicit high order RungeVKutta methods to discontinuous-Galerkin schemes. J. Comput. Phys., 225 (2007), pp. 1753-1781.

[9] Robert E. Lynch, John R. Rice, Donald H. Thomas, Direct solution of partial differ-ential equations by tensor product methods. Numer. Math., 6 (1964) 185-199. [10] D. Funaro, D. Gottlieb, A new method of imposing boundary conditions in

pseu-dospectral approximations of hyperbolic equations. Math. Comput. 51 (1988) 599-613.

[11] Delfour, M, Fortin, M, Payr, G, Finite-difference solutions of a non-linear Schr¨odinger