Efficient exponential splitting spectral methods for linear Schrödinger equation in the semiclassical regime

Journal Pre-proof Efficient exponential splitting spectral methods for linear Schrödinger equation in the semiclassical regime

Wansheng Wang, Jiao Tang

PII:

S0168-9274(20)30041-6

DOI:

https://doi.org/10.1016/j.apnum.2020.02.006

Reference:

APNUM 3771

To appear in:

Applied Numerical Mathematics

Received date:

10 March 2019

Revised date:

28 January 2020

Accepted date:

9 February 2020

Please cite this article as: W. Wang, J. Tang, Efficient exponential splitting spectral methods for linear Schrödinger equation in the semiclassical regime, Appl. Numer. Math. (2020), doi: https://doi.org/10.1016/j.apnum.2020.02.006.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.

EFFICIENT EXPONENTIAL SPLITTING SPECTRAL METHODS ¨ FOR LINEAR SCHRODINGER EQUATION IN THE SEMICLASSICAL REGIME∗ WANSHENG WANG∗† AND JIAO TANG‡ Abstract. The design of efficient numerical methods, which produce an accurate approximation of the solutions, for solving time-dependent Schr¨ odinger equation in the semiclassical regime, where the Planck constant ε is small, is a formidable mathematical challenge. In this paper a new method is shown to construct exponential splitting schemes for linear time-dependent Schr¨ odinger equation with a linear potential. The local discretization error of the two time-splitting methods constructed here is O(max{Δt3 , Δt5 /ε}), while the well-known Lie-Trotter splitting scheme and the Strang splitting scheme are O(Δt2 /ε) and O(Δt3 /ε), respectively, where Δt is the time step-size. The global error estimates of new exponential splitting schemes with spectral discretization suggests that larger time step-size is admissible for obtaining high accuracy approximation of the solutions. Numerical studies verify our theoretical results and reveal that the new methods are especially efficient for linear semiclassical Schr¨ odinger equation with a quadratic potential. Key words. Time-splitting methods; time-dependent Schr¨ odinger equation; semiclassical regime; exponential splitting; error estimates; spectral methods AMS subject classifications. 65M70, 35Q41, 65L05

1. Introduction. The goal of this paper is to construct new time-splitting methods for linear semiclassical Schr¨ odinger equation whose computation presents numerous enduring challenges [27]. This type of equation plays central role in a wide range of applications and is the fundamental model of quantum mechanics [23]. We consider the time-dependent Schr¨odinger equation in a single space variable i

∂u 2 ∂ 2 u ˜ − V (˜ x)u, =− 2m ∂ x ˜2 ∂ t˜

x ˜ ∈ [−(2m) ˜ −1/2 , (2m) ˜ −1/2 ],

0 ≤ t˜ ≤ T˜,

(1.1)

with an initial condition and periodic boundary conditions, where u = u(x, t) is the wave function, the periodic function V˜ is a smooth electrostatic potential and m ˜ is the mass of the underlying particle, a small quantity, but substantially larger than , the reduced Planck constant,  ≈ 1.05457168 × 10−34 . After the scale transformation √ ˜ x/ε and t = t˜/ε, (1.1) is replaced with x = 2m˜ iε

∂uε ∂ 2 uε = −ε2 − V (x)uε , ∂t ∂x2

x ∈ Ω = [−1, 1],

0 ≤ t ≤ T,

(1.2)

where 10−8 ≤ ε ≤ 10−2 is a small dimensionless parameter. The wave function uε (x, t) is used to compute primary physical quantities, usually referred to as observables (see, e.g., [4, 22]), such as the position density, the current density and the energy density: nε = |uε (x, t)|2 , J ε = 2εIm(uε (x, t)∇uε (x, t)),

position density; current density;

eε = 2ε|∇uε (x, t)|2 − 2V (x)|uε (x, t)|2 , energy density.

(1.3) (1.4) (1.5)

∗ This work was supported by a grant from the National Natural Science Foundation of China (Grant No. 11771060, 11371074). † Department of Mathematics, Shanghai Normal University, Shanghai, 200234, China ([email protected]). ‡ School of Mathematics and Statistics, Changsha University of Science & Technology, 410076, Hunan, China.

1

2

W. WANG AND J. TANG

The small size of ε poses a huge challenge for numerical computations of (1.2), since uε features oscillations with frequency 1/ε in space and time, which strain computational resources when a naive approach is applied to numerically solve (1.2) in the semiclassical regime ε  1. In [30] and [31], Markowich et al. utilised the Wigner measure to study the finite difference approximation to the Schr¨ odinger equation with small ε and showed that to guarantee good approximations to all observables, one needs the following constraint h = o(ε),

Δt = o(ε),

(1.6)

for the best combination of the time and space discretizations, where h and Δt are the space and time step-sizes, respectively. To obtain an accurate L2 -approximation of the wave function itself much more restrictive conditions are needed. This motivates the researchers to pursue alternative approaches such as the exponential splitting methods (see, e.g., [21, 32, 29, 24]), WKB (Wentzel-KramersBrillouin)-based splitting schemes (see, e.g., [7, 11, 12, 15, 16]), et. al. To explore the exponential splitting methods, we rewrite (1.2) as an abstract form uε (t) = Auε (t) + Buε (t),

0 ≤ t ≤ T,

(1.7)

with unbounded linear operator A comprising the second derivative operator, i.e., A = iε∂x2 , and multiplication operator B defined by the potential, i.e., B = iε−1 V (x). Two well-known exponential splitting methods are the Lie-Trotter splitting eΔtB eΔtA

(1.8)

and Strang splitting 1

1

e 2 ΔtA eΔtB e 2 ΔtA .

(1.9)

The local splitting errors of Lie-Trotter splitting and Strang splitting ware shown in [4] to be O(Δt2 /ε) and O(Δt3 /ε), respectively. The following global error of the splitting spectral methods, i.e., the Lie-Trotter splitting or the Strang splitting for time discretization and spectral method for space discretization, is also provided in [4]:  m h T CT Δtγ , (1.10) + uε (tn ) − IM U ε,n L2 ≤ Gm Δt 2ε ε under assumption (4.7) (see Section 4), where IM U ε,n is the continuous approximation given by the splitting spectral method, m is an any positive integer which is independent of the three small quantities, ε, h, and Δt, the positive constants C and Gm are also independent of the three small quantities, ε, h, and Δt. Here γ = 1 for the Lie-Trotter splitting spectral method and γ = 2 for the Strang splitting spectral method. The theoretical results (1.10) and numerical experiments carried out in [4] (see also [5] for nonlinear Schr¨odinger equation) confirm the meshing strategy h = O(ε),

Δt = o(ε),

(1.11)

giving L2 -approximation of the wave function uε for the Lie-Trotter splitting spectral method. For the Strang splitting spectral method, the authors omitted the proof because of more complicated calculations and believed that one can establish an estimate at least as good as the one for the Lie-Trotter splitting spectral method. From

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

3

(1.10) and numerical results presented in [4], we know that to obtain an accurate L2 approximation of the wave function, the grid scales for the Strang splitting spectral method should satisfy h = O(ε),

Δt = O(ε).

(1.12)

It is apparent that the Strang splitting spectral method can be comparable to the Lie-Trotter splitting spectral method, and has a high computational efficiency. We also observe that the error bound in the first term on the right hand side of (1.10) is obviously deteriorated as the time-step becomes small. Later on Descombes and Thalhammer [18, 19] gave an exact local error representation of the Lie-Trotter splitting spectral method. For nonlinear Schr¨odinger equation in the semiclassical regime, Carles in [11] proved a global error estimate for a Lie-Trotter splitting operator. The error bounds of high order exponential splitting methods for linear and nonlinear semiclassical Schr¨ odinger equation were also obtained in [35] and [36]. The a posteriori error estimates for Lie-Trotter and Strang splitting spectral methods were investigated in [28, 1] and adaptive splitting methods based on a posteriori local error estimators are proposed in [1]. The splitting schemes for semiclassical Schr¨ odinger equation on an unbounded domain were also investigated in [38]. Li and Xiao [37] studied time-splitting finite difference method for semiclassical Gross-Pitaevskii equation in supercritical case. Multiscale finite element methods were used to solve the semiclassical Schr¨ odinger equation with time-dependent potentials in [17]. The error bounds of the Lie-Trotter splitting and Strang splitting for the Dirac equation in the nonrelativistic limit regime in the absence of external magnetic potentials, with a small parameter, was established in [6]. Boscarino, Pareschi, and Russo recently studied asymptotic properties of the Implicit-Explicit Runge-Kutta methods for hyperbolic systems with a critical parameter (see, e.g., [8, 9]). The standard generalisation of the Strang splitting bears the form eα1 ΔtA eβ1 ΔtB eα2 ΔtA · · · eαs ΔtA eβs ΔtB eαs ΔtA · · · eα2 ΔtA eβ1 ΔtB eα1 ΔtA ,

(1.13)

with parameters αi , βi , i = 1, 2, · · · , s. Observe that to attain significant order an inordinately large number of exponentials is required. For example, the simplest means toward a high-order splitting, the Yo˘ sida method (see, e.g., [32, 39]), needs to compute 2 · 3p−1 + 1 exponentials to attain order 2p. Recently, based on the idea of asymptotic splitting, Bader et. al. [2] introduced and analyzed exponential splitting of the form eΔt(A+B) = eR0 eR1 · · · eRs eTs+1 eRs eR1 · · · eR0 + O(ε2s+2 ),

(1.14)

with R0 = R0 (Δt, ε, A, B) = O(ε0 ), Rk = Rk (Δt, ε, A, B) = O(ε2k−2 ), Ts+1 = Ts+1 (Δt, ε, A, B) = O(ε2s ).

k = 1, 2, · · · , s,

Note that the terms in Rk are the same order as ε2k−2 and the splitting can be obtained by symmetric Baker-Campbell-Hausdorff formula (usually known in abbreviate form as the symmetric BCH formula [24]). The solution of the semi-classical Schr¨ odinger equation subject with explicitly time-dependent potentials was further investigated based on a combination of the Zassenhaus decomposition with the Magnus expansion of the time-dependent Hamiltonian (see, e.g., [3, 25, 26]).

4

W. WANG AND J. TANG

In this paper, we consider another type of splitting based on symmetric BCH formula. More precisely, we construct high-order exponential splitting approximation eΔt[

s j=0

αj Δtj fj (B)]A Δt[

e

s j=0

βj Δtj gj (B)]B Δt[

e

s j=0

αj Δtj fj (B)]A

,

(1.15)

,

(1.16)

or eΔt[

s j=0

βj Δtj gj (B)]B Δt[

e

s j=0

αj Δtj fj (B)]A Δt[

e

s j=0

βj Δtj gj (B)]B

where fj and gj are operators acting on B. More specifically, we introduce and analyse the following splitting schemes 1

2

eΔt(A+B) ≈ e 2 Δt[1+βΔt

2 ∂x V ]A iε−1 Δt[V −αΔt2 (∂x V )2 ]

e

1

2

−1

Δt[V −αΔt2 (∂x V )2 ]

e 2 Δt[1+βΔt

2 ∂x V ]A

(1.17)

and 1

eΔt(A+B) ≈ e 2 iε

−1

2 Δt[V −αΔt2 (∂x V )2 ] Δt[1+βΔt2 ∂x V ]A

e

1

e 2 iε

.(1.18)

It is shown that the local splitting error of (1.17) is O(max{Δt3 , Δt5 /ε}) whenever 1 and β = 13 . It should be α = − 16 and β = − 16 , and so is (1.18) whenever α = 12 3 pointed out that to attain the local splitting error of O(Δt ), the splitting (1.14) needs to compute at least five exponentials but three for ours when we solve linear problems with a harmonic potential. The paper is organized as follows. In Section 2 we introduce and analyze our exponential splitting scheme in an infinite-dimensional operator setting, and show their local splitting errors based on symmetric BCH formula. In Section 3, with the exponential splitting for time discretization and spectral method for space discretization, the full discretization approximations of linear Schr¨odinger equation are considered and two new exponential splitting spectral schemes are presented. Section 4 is devoted to the stability analysis and error estimates of the two new exponential splitting spectral schemes. To avoid the term T /Δt, which deteriorates the error bound in (1.10), we use Lay Windermere’s fan and show that for both our exponential splitting spectral schemes, when the space step size h satisfies h = O(ε), the following global error estimate holds   m  4 h ε ε,n 2 Δt , (1.19) + max Δt , u (tn ) − IM U L2 ≤ C 2ε ε where C is independent of ε, h, and Δt. In Section 5 numerical results of the new exponential splitting spectral methods are presented and compared with the Strang splitting spectral methods for a linear problem with a quadratic potential. The last section, Section 6, will contain a few concluding remarks. 2. New exponential splitting for semi-discretization approximations. In this section we introduce our new exponential splitting approximations and analyse their local error. For this purpose, we need the well-known symmetric BCH formula [24]. 2.1. The symmetric BCH formula. Let X and Y be two terms in a Lie algebra g. Unless X and Y commute, it is in general not true that et(X+Y ) = etX etY . It is well known that the commutator [X, Y ] of operators X and Y action on u is defined by [X, Y ]u := X(Y u) − Y (Xu).

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

5

Here for simplicity’s sake we assume u is periodic function in C ∞ [−1, 1], although the definition can be extended to functions of lower smoothness in a straightforward manner. Then the symmetric BCH formula reads 1

1

e 2 tX etY e 2 tX = esBCH(tX,tY ) ,

(2.1)

sBCH(tX, tY ) = tS1 + t3 S3 + t5 S5 + · · ·

(2.2)

where

with S1 = X + Y ; 1 1 [Y, [Y, X]] ; S3 = − [X, [X, Y ]] + 24 12 1 7 [X, [X, [X, [X, Y ]]]] − [Y, [Y, [Y, [Y, X]]]] S5 = 5760 720 1 1 [X, [Y, [Y, [Y, X]]]] + [Y, [X, [X, [X, Y ]]]] + 360 360 1 1 [X, [X, [Y, [Y, X]]]] + [Y, [Y, [X, [X, Y ]]]] . − 480 120

(2.3)

Because (2.1) is palindromic, only odd powers of t feature in the expansion. The expansion (2.1) can be computed to an arbitrary power of t using an algorithm from Casas and Murua [13]. It is a common practice that let X = B and Y = A such that S = O(ε−1 ). Instead of such simple choice, we take 1 = A + B. In this case, S3  s s j X = [ j=0 βj t gj (B)]B and Y = [ j=0 αj tj fj (B)]A with fj and gj being operators acting on B such that S1 = A + B, S3 = O(1) and S5 = O(ε−1 ). Thus, high order splitting schemes can be obtained. This leads to the following choices. 2.2. Two new exponential splitting schemes. Since the vector field in the linear Schr¨ odinger equation (1.2) is a linear combination of the action of two operators, ∂x2 and the multiplication by the interaction potential V , to compute commutators generated by ∂x2 and V , we first need to compute [V, ∂x2 ]u = V (∂x2 u) − ∂x2 (V u) = −(∂x2 V )u − 2(∂x V )∂x u, which implies that [V, ∂x2 ] = −(∂x2 V ) − 2(∂x V )∂x . Let the time step be tn = nΔt, n = 0, 1, 2, · · · , with the time stepsize Δt > 0. We set   Y = iε(1 + βΔt2 ∂x2 V )∂x2 , (2.4) X = iε−1 V (x) − αΔt2 (∂x V )2 , where α =

1 12

and β = 13 . Then we have

 [X, Y ] = (1 + βΔt2 ∂x2 V ) ∂x2 V + 2∂x V ∂x

 −2αΔt2 (∂x2 V )2 − 2αΔt2 ∂x V ∂x3 V − 4αΔt2 ∂x V ∂x2 V ∂x .

(2.5)

Furthermore, one can obtain [X, [X, Y ]] = 2iε−1 (1 + βΔt2 ∂x2 V )   × −(∂x V )2 + 4αΔt2 (∂x V )2 ∂x2 V − 4α2 Δt4 (∂x V )2 (∂x2 V )2 (2.6)

6

W. WANG AND J. TANG

and  [Y, [Y, X]] = iε(1 + βΔt2 ∂x2 V )2 −4∂x2 V ∂x2 − 4∂x3 V ∂x − ∂x4 V +6αΔt2 (∂x3 V )2 + 8αΔt2 ∂x2 V ∂x4 V + 2αΔt2 ∂x V ∂x5 V  +24αΔt2 ∂x2 V ∂x3 V ∂x + 8αΔt2 (∂x2 V )2 ∂x2 + 8αΔt2 ∂x V ∂x3 V ∂x2  −iεβΔt2 (1 + βΔt2 ∂x2 V )∂x4 V ∂x2 V + 2∂x V ∂x  −2αΔt2 (∂x2 V )2 − 2αΔt2 ∂x V ∂x3 V − 4αΔt2 ∂x V ∂x2 V ∂x  −2iεβΔt2 (1 + βΔt2 ∂x2 V )∂x3 V ∂x3 V + 3∂x2 V ∂x + 2∂x V ∂x2 − 6αΔt2 ∂x2 V ∂x3 V − 2αΔt2 ∂x V ∂x4 V  −6αΔt2 (∂x2 V )2 ∂x − 6αΔt2 ∂x V ∂x3 V ∂x − 4αΔt2 ∂x V ∂x2 V ∂x2 . (2.7) Substituting (2.5), (2.6), and (2.7) into the symmetric BCH formula (2.1) to get sBCH(ΔtX, ΔtY )  

= Δt iε−1 V (x) − αΔt2 (∂x V )2 + iε(1 + βΔt2 ∂x2 V )∂x2 1 − iε−1 Δt3 (1 + βΔt2 ∂x2 V ) 12   × −(∂x V )2 + 4αΔt2 (∂x V )2 ∂x2 V − 4α2 Δt4 (∂x V )2 (∂x2 V )2 1 + iεΔt3 (1 + βΔt2 ∂x2 V )(−4∂x2 V ∂x2 − 4∂x3 V ∂x − ∂x4 V + · · · ) 12    1 = Δt iε−1 V (x) + iε∂x2 + − α iε−1 Δt3 (∂x V )2 12   1 1 iεΔt3 ∂x2 V ∂x2 − (4α − β)iε−1 Δt5 ∂x2 V (∂x V )2 + β − 12 3   1 1 1 − iεΔt3 ∂x3 V ∂x − iεΔt3 ∂x4 V − iεΔt5 8β(∂x2 V )2 ∂x2 + · · · . 3 12 12

(2.8)

Classical WKB analysis of the semicalssical Schr¨odinger equation [27] shows that the solution develops spatial oscillations of order O(ε−1 ), provided the initial datum in itself is ε-oscillatory, that is, uε0 =

ε iS /ε n0 e 0 .

As a consequence, each ∂xd scales as O(ε−d ). Taking this into account, we set α = and β = 13 in (2.8) and obtain the exponential splitting scheme −1

eΔt(A+B) ≈ e 2 iε Δt[V − 12 Δt (∂x V ) ] eΔt[1+ 3 Δt 1 1 = e 2 ΔtX eΔtY e 2 ΔtX 1

1

2

2

1

2

2 ∂x V ]A

1

e 2 iε

−1

1 12

1 Δt[V − 12 Δt2 (∂x V )2 ]

(2.9)

with approximation error O(max{Δt3 , Δt5 /ε}). Note that the approximation error scales as Δt3 whenever Δt = O(εσ ), 1/2 ≤ σ ≤ 1. For simplicity, the scheme (2.9) will be referred to as BAB type of schemes.

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

7

Let us now consider the splitting of the form ABA. Analogously, we have sBCH(ΔtY, ΔtX)  

= Δt iε−1 V (x) − αΔt2 (∂x V )2 + iε(1 + βΔt2 ∂x2 V )∂x2 1 + iε−1 Δt3 (1 + βΔt2 ∂x2 V ) 6  × −(∂x V )2 + 4αΔt2 (∂x V )2 ∂x2 V − 4α2 Δt4 (∂x V )2 (∂x2 V )2 1 − iεΔt3 (1 + βΔt2 ∂x2 V )(−4∂x2 V ∂x2 − 4∂x3 V ∂x − ∂x4 V + · · · ) 24     1 = Δt iε−1 V (x) + iε∂x2 − + α iε−1 Δt3 (∂x V )2 6   1 1 −1 5 2 2 iεΔt3 ∂x2 V ∂x2 + (4α − β)iε Δt (∂x V ) ∂x V + β + 6 6   1 1 1 + iεΔt3 ∂x3 V ∂x + iεΔt3 ∂x4 V + iεΔt5 8β(∂x2 V )2 ∂x2 + · · · . 6 24 24

(2.10)

Considering the fact that the solution develops spatial oscillations of order O(ε−1 ), provided the initial datum in itself is ε-oscillatory, we set α = − 16 and β = − 16 and obtain the following splitting scheme: 1

1

eΔt(A+B) ≈ e 2 Δt[1− 6 Δt

2

2 ∂x V ]A iε−1 Δt[V + 16 Δt2 (∂x V )2 ]

e

1

1

e 2 Δt[1− 6 Δt

2

2 ∂x V ]A

,

(2.11)

with approximation error O(max{Δt3 , Δt5 /ε}). Similarly, the approximation error scales as Δt3 whenever Δt = O(εσ ), 1/2 ≤ σ ≤ 1. Observe that if we set α = β = 0 in (2.4), the splittings (1.18) and (1.17) reduce respectively to the Strang splittings 1

eΔt(A+B) ≈ e 2 iε

−1

ΔtV ΔtA

e

1

e 2 iε

−1

ΔtV

(2.12)

and 1

eΔt(A+B) ≈ e 2 ΔtA eiε

−1

ΔtV

1

e 2 ΔtA .

(2.13)

It should be pointed out that for constant potential V (x) ≡ V = constant, (2.9) and (2.11) reduce to (2.12) and (2.13), respectively, just as the Strang splittings reduce to the Lie-Trotter splitting in this case. 2.3. Local splitting error of exponential splitting algorithms. Let  · L2 be the usual L2 -norm on the interval Ω = [−1, 1] which is defined as  u2L2

1

= −1

|u(x)|2 dx,

u ∈ L2 (Ω).

(2.14)

Now we are ready to give the local splitting error of (2.9) and (2.11). Let w be the solution obtained from the exponential splitting schemes after one time step with the exact initial data at tn without spatial discretization, that is 1

w = e 2 iε =e

−1

2 Δt[V −αΔt2 (∂x V )2 ] Δt[1+βΔt2 ∂x V ]A

1 2 ΔtX

e

e

ΔtY

e

1 2 ΔtX

ε

u (x, tn ).

1

e 2 iε

−1

Δt[V −αΔt2 (∂x V )2 ] ε

u (x, tn ) (2.15)

8

W. WANG AND J. TANG

On the other hand, the exact solution uε (x, tn+1 ) is given by uε (x, tn+1 ) = eΔt(A+B) uε (x, tn ). and β = 13 ), we have  sBCH(ΔtX, ΔtY ) = Δt (A + B) + O max{Δt3 , Δt5 /ε} ,

Then for the new splitting considered here (α =

(2.16)

1 12

(2.17)

and therefore,  uε (tn+1 ) − wL2 = O max{Δt3 , Δt5 /ε} .

(2.18)

Furthermore, if we take Δt = O(εσ ), 1/2 ≤ σ ≤ 1, we have  sBCH(ΔtX, ΔtY ) = Δt (A + B) + O Δt3 ,

(2.19)

and  uε (tn+1 ) − wL2 = O Δt3 .

(2.20)

But for the Strang splitting (2.12) (α = 0 and β = 0), it is easy to get from (2.8) that  3  −1 Δt 2 sBCH(ΔtX, ΔtY ) = Δt iε V (x) + iε∂x + O ε  3 Δt = Δt (A + B) + O . (2.21) ε Consequently, the following holds  uε (tn+1 ) − wL2 = O

Δt3 ε

 .

(2.22)

For the new splitting (2.11), we also have (2.18) and (2.20) whenever Δt = O(εσ ), 1/2 ≤ σ ≤ 1. For the Strang splitting (2.13) of the form ABA, we simi1 larly get (2.22). It is noteworthy that if we take Δt = O(ε 2 ) for the Strang splittings (2.12) and (2.13), the order of the splitting error is only 1. 3. Exponential splitting spectral methods. The splittings (2.9) and (2.11) are expressed in operatorial terms. To obtain the full discretization algorithms we must consider the space discretization and replace ∂x2 with an appropriate differentiation matrix, acting on an M -dimensional space. The finite difference, the spectral collocation, and the pseudo-spectral methods are considered in the literatures. Faced with this embarrassment of riches, we opt for the spectral methods. It is common to use spectral discretization in the numerical solution of the Schr¨ odinger equation (see, e.g., [21, 27]), since they exhibit spectral convergence: for sufficiently large M the error decays faster than M −r for any r > 0. 3.1. The exponential splitting spectral method of the first kind (BAB). We choose the spatial mesh size h = 2/M with M being an even positive integer. Let Ujε,n be the approximation of uε (xj , tn ) and U ε,n be the solution vector with ˆ ε of a vector U ε are defined as components Ujε,n . The Fourier coefficients U l ˆlε = U

M −1  j=0

Ulε e−iμl (xj +1) ,

l=−

M M ,··· , − 1. 2 2

(3.1)

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

9

For a given initial function uε (x, 0) = uε0 (x), we set Ujε,0 = uε (xj , 0) = uε0 (xj ),

j = 0, 1, 2, · · · , M.

(3.2)

Considering spectral discretization for the space approximation, the algorithm (2.9) to compute the numbers Ujε,n+1 from the collection of numbers Ujε,n then reads: 1

1

Ujε,∗ = e 2ε iΔt(V (xj )− 12 Δt Ujε,∗∗ =

1 M

Ujε,n+1 = e



2

(∂x V (xj ))2 )

Ujε,n ,

j = 0, 1, · · · , M − 1,

M/2−1

e−iεΔt(1+ 3 Δt 1

2

2 ∂x V )μ2l

ˆ ε,∗ eiμl (xj +1) , j = 0, 1, · · · , M − 1, U l

l=−M/2

1 2ε iΔt(V

1 (xj )− 12 Δt2 (∂x V (xj ))2 )

Ujε,∗∗ ,

j = 0, 1, · · · , M − 1,

(3.3)

M with μl = πl, l = − M 2 , · · · , 2 − 1. The above algorithm can be rewritten in the abstract formula 1

1

U ε,n+1 = e 2 ΔtXh eΔtYh e 2 ΔtXh U ε,n .

(3.4)

Observe that the only time discretization error of this method is the splitting error, which is third order in Δt for any small 0 < ε  1. 3.2. The exponential splitting spectral method of the second kind (ABA). Using the spectral method for the space disretization, we split the Schr¨ odinger equation (1.2) via the splitting (2.11) Ujε,∗ Ujε,∗∗

1 = M =e

Ujε,n+1 =



M/2−1

e− 2 iεΔt(1− 6 Δt 1

1

2

2 ∂x V )μ2l

l=−M/2

ε−1 iΔt(V (xj )+ 16 Δt2 (∂x V (xj ))2 )

1 M

ˆ ε,n eiμl (xj +1) , j = 0, 1, 2, · · · , M − 1, U l



Ujε,∗ ,

j = 0, 1, 2, · · · , M − 1,

(3.5)

M/2−1

e− 2 iεΔt(1− 6 Δt 1

1

2

2 ∂x V )μ2l

ˆ ε,∗∗ eiμl (xj +1) , j = 0, 1, 2, · · · , M − 1. U l

l=−M/2

3.3. The Strang splitting spectral methods (SP2). For comparison, we also introduce the Strang splitting spectral methods. The Strang splitting spectral algorithm of the form BAB to compute the numbers Ujε,n+1 from the collection of numbers Ujε,n reads: Ujε,∗ = eiΔtV (xj )/(2ε) Ujε,n , Ujε,∗∗ =

1 M



j = 0, 1, 2, · · · , M − 1,

M/2−1 2 ˆ ε,∗ eiμl (xj +1) , e−iεΔtμl U l

j = 0, 1, 2, · · · , M − 1,

(3.6)

l=−M/2

Ujε,n+1 = eiΔtV (xj )/(2ε) Ujε,∗∗ ,

j = 0, 1, 2, · · · , M − 1.

Similarly, we can define Strang splitting spectral algorithm of the form ABA. 4. Stability and error estimates. To study the stability and error estimates for the proposed methods, we introduce the discrete l2 -norm  · l2 on the interval Ω which is defined as U 2l2 =

M −1 2  |Uj |2 , M j=0

U = (U1 , U2 , · · · , UM −1 )T .

(4.1)

10

W. WANG AND J. TANG

It is useful to note the following relation for every periodic function f IM f 2L2

=

f 2l2

M −1 2  2 = |f (xj )| M j=0

(4.2)

where IM f stands for the trigonometric interpolant of f on {x0 , x1 , · · · , xM }, defined by

IM f (x) =

1 M



M/2−1

fˆl eiμl (x+1) ,

fˆl =

M −1 

f (xj )e−iμl (xj +1) ,

j=−

j=0

l=−M/2

M M , · · · , −1. 2 2

To measure the errors between IM f and f in Sobolev spaces, we denote Hpr (Ω) the subspace of H r (Ω), which consists of functions with derivatives up to r − 1 being 2-periodic. The norm and semi-norm of Hpr (Ω) are equivalent to  f r =

∞ 



 12 (1 + l2 )r |fˆl |2

,

|f |r =

l=−∞

 12

∞ 

l2r |fˆl |2

.

l=−∞

In view of the property of spectral methods (see, e.g. [34]), for u ∈ H r with r > 0, we have the following    1 1 1 1    IM e 2 ΔtXh eΔtYh e 2 ΔtXh − e 2 ΔtX eΔtY e 2 ΔtX f 

L2

≤ CΔtM −r |f |r .

We also have the following lemma. Lemma 4.1 (see, e.g., [33, 34]). Suppose that f ∈ H r (Ω) with r > following estimate holds IM f − f k ≤ CM k−r |f |r ,

0 ≤ k ≤ r.

1 2,

(4.3)

then the

(4.4)

4.1. Stability. We are ready to prove the following theorem, which shows that the total charge is conserved. Theorem 4.2. The exponential splitting spectral schemes (3.3) and (3.5) are unconditionally stable, that is, under any mesh size h and time step Δt, U ε,n l2 = uε0 l2 ,

n = 1, 2, · · · ,

(4.5)

and consequently IM U ε,n L2 = IM U ε,0 L2 ,

n = 1, 2, · · · .

(4.6)

Here, IM U ε,n stands for the trigonometric polynomial interpolating {(x0 , U0ε,n ), (x1 , U1ε,n ), ε,n )}. · · · , (xM , UM Proof. In view of the orthogonal property of {ϕl (x) := eiμl (x+1) }, the proof is standard (see, e.g., [4, 20]) and we omit it.

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

11

4.2. Error estimates. In this subsection we establish error estimates for our two schemes. To this end, we assume that the solution uε = uε (x, t) of (1.2) and the potential V (x) in (1.2) are C ∞ (R) and 2-periodic. In classical WKB analysis, it is assumed that the initial datum in itself is ε-oscillatory. Then the initial value problem for the linear Schr¨ odinger equation propagates one ε-scale of oscillations. Therefore, we here assume that the solution oscillates in space and time with wavelength ε, that is, there are positive constants Cm > 0, Dm > 0 independent of ε, x, t such that  m +m  m    ∂ 1 2 ε  d  Cm1 +m2    ≤ m1 +m2 , V ≤ Dm ,  ∂xm1 ∂tm2 u    ∞ m ε dx C([0,T ];L2 (Ω)) L (Ω) for all m, m1 , m2 ∈ N ∪ {0}. Let XM := span{ϕl (x) : − M 2 ≤ l ≤ L2 -orthogonal projection, defined by

M 2

(4.7)

− 1} and let PM : L2 (Ω) → XM be the

(u − PM u, v)L2 = 0,

∀v ∈ XM .

Then PM can be written as 

M/2−1

(PM u)(x, t) =

u ˜l (t)ϕl (x),

(4.8)

l=−M/2

 1 M where u ˜l (t) = 2L u(x)ϕ¯l (x)dx, − M 2 ≤ l ≤ 2 − 1. The following lemma bounds Ω the error between PM u and u. Lemma 4.3 (see, e.g., [34]). Suppose that u ∈ H r (Ω) with r > 0, then the following estimate holds PM u − uL2 ≤ CM −r |u|r .

(4.9)

The following lemma provides an upper bound for the error between the projection operator PM and the interpolation operator IM . Lemma 4.4 (see, e.g., [34]). Suppose that u ∈ H r (Ω) with r > 0, then the following estimate holds PM u − IM uL2 ≤ CM −r |u|r .

(4.10)

With these preparations, we have the following error estimate. Theorem 4.5. Let uε = uε (x, t) be the exact solution of (1.2), and U ε,n be the discrete approximation given by (3.3) or (3.5). Assuming h = O(ε), then under assumption (4.7), for any Δt > 0, we have for all positive integers m ≥ 1 which is independent of ε, h, Δt, and tn ∈ [0, T ] that   m  h Δt4 , (4.11) + max Δt2 , uε (tn ) − IM U ε,n L2 ≤ C 2ε ε where C is independent of ε, h, and Δt. Proof. We only show the conclusion for the exponential splitting spectral method (3.3) of the first kind (BAB). The proof for the algorithm (3.5) is similar.

12

W. WANG AND J. TANG

We first consider the error splitting uε (tn ) − IM U ε,n = uε (tn ) − un + un − PM un + PM un − IM U ε,n ,

(4.12)

n

where u is the splitting solution at t = tn satisfying exactly the operator splitting scheme with the exact initial data uε0 = uε (0), that is, n  1 1 (4.13) un = e 2 ΔtX eΔtY e 2 ΔtX uε0 . From Lemma 4.3, we have un − PM un L2 ≤ CM −m |un |r .

(4.14)

Now we direct our attention to estimate u (tn ) − u and PM u − IM U we represent them as, by use of standard Lay Windermere’s fan, ε

n

n

ε,n

. To do this,

uε (tn ) − un n  n−j    1 1 1 1 e 2 ΔtX eΔtY e 2 ΔtX eΔt(A+B) − e 2 ΔtX eΔtY e 2 ΔtX uε (tj−1 ) = j=1

and PM un − IM U ε,n n  n−j  1 1 IM e 2 ΔtXh eΔtYh e 2 ΔtXh = j=1

  1 1 1 1 × PM e 2 ΔtX eΔtY e 2 ΔtX − IM e 2 ΔtXh eΔtYh e 2 ΔtXh uj−1 . Thanks to (2.18), uε (tn ) − un can be estimated by uε (tn ) − un L2 n     1 1   ≤ eC(tn −tj )  eΔt(A+B) − e 2 ΔtX eΔtY e 2 ΔtX uε (tj−1 ) j=1 n 



Δt5 ≤ e C max Δt , ε j=1   Δt4 ≤ C max Δt2 , . ε C(tn −tj )

L2



3

(4.15)

For PM un − IM U ε,n , using Lemma 4.4 and (4.3), we can bound it as PM un − IM U ε,n L2    1 1 1 1   ≤  PM e 2 ΔtX eΔtY e 2 ΔtX − IM e 2 ΔtXh eΔtYh e 2 ΔtXh un−1   n−1  n−j  1 1 IM e 2 ΔtXh eΔtYh e 2 ΔtXh +   j=1    1 1 1 1  × PM e 2 ΔtX eΔtY e 2 ΔtX − IM e 2 ΔtXh eΔtYh e 2 ΔtXh uj−1  ≤ CM −m |un−1 |m +

n 

eC(tn −tj ) ΔtCM −m

j=1

≤ CM

−m

sup 0≤j≤n−1

|uj |m .

sup

|uj |m

0≤j≤n−2

(4.16)

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

13

Since uj has the same space regularity as uε (tj ), we have |un |m ≤ Cε−m with C being independent of ε. Then combine (4.12), (4.14), (4.15) and (4.16) to get   Δt4 + C(M ε)−m , (4.17) uε (tn ) − IM U ε,n L2 ≤ C max Δt2 , ε which completes the proof. For Strang splitting spectral methods, we also have the following estimate. Theorem 4.6. Let uε = uε (x, t) be the exact solution of (1.2), and U ε,n be the discrete approximation given by (3.6). Under assumption (4.7), and assuming h = O(ε), we have for all positive integers m ≥ 1 which is independent of ε, h, Δt and tn ∈ [0, T ] that   m h Δt2 uε (tn ) − IM U ε,n L2 ≤ C , (4.18) + 2ε ε where C is independent of ε, h, and Δt. We conclude this section with a few remarks about our results. Our first remark is that compared to Theorem 4.6, the theoretical results obtained in Theorem 4.5 reveal an interesting property of our new exponential splitting spectral schemes: they have convergence order 4 (high order) in time for large time step-size 1 (Δt = O(ε 2 )) and have uniform convergence order 2 (independent of ε) for small time step-size. The second remark is about the factor T /Δt, which deteriorates the bound in (1.10) as the time-step becomes small. We observe that it does not appear in the error bounds in Theorems 4.5 and 4.6. Finally, it should be emphasized that to guarantee an accurate L2 -approximation of the wave function, the spatial discretization parameter h must scale with the semiclassical parameter ε for both our exponential splitting spectral methods. 5. Numerical Examples. In this section, we illustrate the behaviour of the schemes (3.3) and (3.5) introduced in Section 3 and compare their properties to those of the Strang splitting spectral method (3.6) for (1.2). Following the numerical observations of Example 3 in [4] for the Strang splitting, we illustrate the high accuracy approximation of our new splitting methods with respect to the critical parameter 0 < ε  1 when applied to the time-dependent linear Schr¨ odinger equation iε

∂uε 1 1 = − ε2 Δuε + x2 uε , −2 ≤ x ≤ 2, ∂t 2 2 ε −25(x−1/2)2 i(1+x)/ε e , u0 = e

0 ≤ t ≤ T,

(5.1) (5.2)

subject to periodic boundary conditions. We set T = 0.64 as in [4]. Although it seems that the form of equation (5.1) is different from the form (1.2) considered here, they can transform each other. The values of the critical parameter ε are chosen in the range 6.25 × 10−5 to 4 × 10−2 . For each fixed ε, a reference solution uεr is computed by using the Strang splitting spectral method (2.12) with a very fine mesh, 1 , and a very small time step, e.g., Δt = 0.00001. Tables 5.1-5.3 show e.g., h = 32768 the errors uεr − uε,h,Δt l2 produced by our new schemes at t = 0.64 for different combinations of ε, h, and Δt. For the purpose of comparison, the errors produced by the Strang splitting spectral methods, which have been shown in [4], are also

14

W. WANG AND J. TANG

Table 5.1 The errors uεr − uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with ε = 0.04 at t = 0.64 and their convergence orders with respect to Δt.

Methods

Strang BAB

Δt

New ABA

h=

1 16

10−3

h=

1 64

1 256

7.555963 ·

Order −

10−3

7.556521 ·

0.04

0.709733

4.800879 · 10−4

4.697152 · 10−4

4.697153 · 10−4

2.003879

0.01

0.709797

1.041984 · 10−4

2.934735 · 10−5

2.934743 · 10−5

2.000240

0.709801

1.000429 ·

10−4

1.834209 ·

10−6

1.834283 ·

10−6

1.999973

6.038058 ·

10−3

6.036899 ·

10−3

6.036899 ·

10−3

−

10−4

3.751863 ·

10−4

3.751861 ·

10−4

2.004065

0.711782

7.555963 ·

10−3

h=

0.708577

0.16

New BAB

1 4

0.16

0.0025

Strang ABA

h=

0.04

0.709919

3.884670 ·

0.01

0.709808

1.027791 · 10−4

2.344107 · 10−5

2.344089 · 10−5

2.000254

0.0025

0.709802

1.000424 ·

10−4

10−6

10−6

2.000009

0.16

0.708768

1.031813 · 10−4

2.998464 · 10−5

2.998468 · 10−5

0.709745

9.979304 ·

10−5

1.173653 ·

10−7

1.173743 ·

10−7

3.998482

10−4

4.141435 ·

10−9

3.717781 ·

10−9

2.490266

0.04

1.465222 ·

1.465038 ·

−

0.01

0.709798

1.000141 ·

0.0025

0.709801

1.000282 · 10−4

3.656301 · 10−9

3.390904 · 10−9

0.066386

0.16

0.711667

1.179610 · 10−4

2.996040 · 10−5

2.996045 · 10−5

−

0.04

0.709912

1.004810 · 10−4

1.170153 · 10−7

1.169691 · 10−7

4.000394

0.709808

1.000563 ·

10−4

4.216853 ·

10−9

3.675497 ·

10−9

2.496022

1.000308 ·

10−4

3.666237 ·

10−9

3.395276 ·

10−9

0.057205

0.01 0.0025

0.709802

listed here. In these tables, we also present the convergence orders of these numerical 1 in Tables 5.1 and methods, which are computed at the finest spatial mesh (h = 256 1 5.2 and h = 1024 in Tables 5.3). From Tables 5.1-5.3, we observe that the splitting schemes of the form ABA have the same error behaviour as those of the form BAB for both the Strang and our new splitting schemes. As a consequence, we only list the numerical results produced by the schemes of the form BAB in the following numerical experiments. Table 5.4 shows the errors uεr − uε,h,Δt l2 produced by the schemes of the form BAB at t = 0.64 and their convergence orders with respect to Δt 1 and ε = 6.25 × 10−4 . when h = 4096 To show the computational efficiency, we present the CPU times of the schemes for the Schr¨odinger equation with a much smaller semiclassical parameter ε = 6.25×10−5 in Table 5.5. These numerical results on the approximation errors of the wave function uε (x, t) are summarized in Fig. 5.1. The approximation errors of the position density nεr = |uεr |2 are presented in Fig. 5.2. From the theoretical analysis given in this paper and the numerical results shown in Tables 5.1-5.5 and Fig. 5.1, we come to the following remarks: (i) Our two schemes introduced in this paper have much higher accuracy than the frequently used Strang splitting schemes for this special linear problem. In fact, we observe that with the same space step-sizes, the errors of our new methods with Δt = 0.04 are smaller than those of the Strang splitting schemes with Δt = 0.0025. It is illustrated by CPU times in Table 5.5 that

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

15

Table 5.2 The errors uεr − uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with ε = 0.01 at t = 0.64 and their convergence orders with respect to Δt. Methods

Δt

Strang BAB

10−2

h=

1 256

Order −

0.707554

1.454233 · 10−3

1.454233 · 10−3

2.004042

0.01

0.707586

9.085735 · 10−5

9.085736 · 10−5

2.000256

0.707588

5.678705 ·

10−6

5.678711 ·

10−6

1.999984

2.457835 ·

10−2

2.457836 ·

10−2

−

10−3

1.527747 ·

10−3

2.003956

0.708492

2.339846 ·

10−2

0.04

0.04

0.707655

1.527747 ·

0.01

0.707592

9.545035 · 10−5

9.545091 · 10−5

2.000251

0.0025

0.707588

5.965372 ·

10−6

10−6

1.999970

0.16

0.707407

1.116185 · 10−4

1.116359 · 10−4

0.707552

4.343070 ·

10−7

4.339900 ·

10−7

4.003347

10−9

5.404860 ·

10−9

3.163630

5.965933 ·

−

0.01

0.707586

5.368670 ·

0.0025

0.707588

5.111160 · 10−9

5.252668 · 10−9

0.020603

0.16

0.708555

1.116362 · 10−4

1.116181 · 10−4

−

0.04

0.707659

4.343114 · 10−7

4.339943 · 10−7

4.003340

0.707592

5.368672 ·

10−9

5.404860 ·

10−9

3.163638

5.111159 ·

10−9

5.252668 ·

10−9

0.020603

0.01 0.0025

0.707588

h=1/1024

0

1 64

2.339846 ·

0.04

New ABA

h=

0.707427

0.16

New BAB

1 16

0.16

0.0025

Strang ABA

h=

h=1/1024

−2

10

10

−1

10

−4

10

−2

10

−6

10 −3

Error

Error

10

−4

−8

10

10

−10

−5

ε=0.000625 ε=0.0025 ε=0.01 ε=0.04 slope 2

10

−6

10

−7

10

10

ε=0.000625 ε=0.0025 ε=0.01 ε=0.04 slope 4

−12

10

−14

−2

−1

10

10 time step

10

−2

−1

10

10 time step

Fig. 5.1. The errors uεr −uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with different ε at t = 0.64. Left: Strang BAB splitting spectral method; Right: New BAB splitting spectral method.

the new schemes are more efficient than the Strang splitting spectral schemes. (ii) Our error analysis and numerical results confirm that to give L2 -approximation of the wave function, the space step h should scale as O(ε) for all splitting spectral methods presented here. We also find that all numerical data of the exponential splitting spectral methods with a space step-size h satisfying h 2ε < 1 are almost the same as those of these methods with a smaller space step-size h. This implies that the exponential splitting spectral methods have

16

W. WANG AND J. TANG

Table 5.3 The errors uεr − uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with ε = 0.0025 at t = 0.64 and their convergence orders with respect to Δt. Methods

Strang BAB

Δt

New ABA

h=

1 256

10−2

h=

1 1024

Order −

9.152470 ·

0.04

0.707961

5.698667 · 10−3

5.698667 · 10−3

2.002732

0.01

0.708013

3.560410 · 10−4

3.560409 · 10−4

2.000255

0.708015

2.225185 ·

10−5

2.251728 ·

10−5

2.000027

9.815456 ·

10−2

9.815456 ·

10−2

−

10−3

6.118257 ·

10−3

2.001931

0.708046

9.152470 ·

10−2

0.708042

0.16

New BAB

1 64

0.16

0.0025

Strang ABA

h=

0.04

0.708099

6.118257 ·

0.01

0.708018

3.822600 · 10−4

3.822600 · 10−4

2.000247

0.0025

0.708015

2.389042 ·

10−5

10−5

2.000024

0.16

0.708052

4.443512 · 10−4

4.443513 · 10−4

0.707961

1.729432 ·

10−6

1.729473 ·

10−6

4.002612

10−8

1.003947 ·

10−8

3.714253

0.04

2.389046 ·

−

0.01

0.708013

1.015861 ·

0.0025

0.708015

7.314576 · 10−9

7.108638 · 10−9

0.249019

0.16

0.708151

4.444520 · 10−4

4.444521 · 10−4

−

0.04

0.708099

1.729456 · 10−6

1.729497 · 10−6

4.002766

0.708018

1.015862 ·

10−8

1.003947 ·

10−8

3.714263

7.314576 ·

10−9

7.108638 ·

10−9

0.249019

0.01 0.0025

0.708015

Table 5.4 The errors uεr − uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with ε = 0.000625 at t = 0.64 and their convergence orders with respect to Δt. Methods

Strang BAB

Δt

1 256

h=

1 1024

10−1

h=

1 4096

Order −

0.708048

3.549128 ·

0.04

0.708068

2.276224 · 10−2

2.276224 · 10−2

1.981375

0.708046

1.422300 ·

10−3

1.422300 ·

10−3

2.000172

8.889057 ·

10−5

8.889046 ·

10−5

2.000027

10−3

1.776873 ·

10−3

−

0.0025

0.708044

3.549128 ·

10−1

0.16

0.01

New BAB

h=

0.16

0.708057

1.776873 ·

0.04

0.708069

6.915572 · 10−6

6.915672 · 10−6

4.002628

0.01

0.708046

2.907784 · 10−8

2.905954 · 10−8

3.947357

0.708044

10−9

10−9

1.044202

0.0025

7.292578 ·

6.833080 ·

spectral accuracy for spatial discretization. The main error of the exponential h < 1. splitting spectral methods is the time discretization error whenever 2ε (iii) It is confirmed that for any fixed ε, the Strang splitting spectral methods have order 2 for time approximation. For our new methods, from the theoretical analysis, we know that for any fixed ε, if Δt2 = O(ε), the approximation order of our new methods is of 4, and if Δt2 = o(ε), the approximation order is uniformly of 2. This is verified by our numerical results and particularly evident for a small parameter ε (see Table 5.5). It should be pointed out that

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

17

Table 5.5 The errors uεr − uε,h,Δt l2 produced by different methods for (5.1)-(5.2) with ε = 6.25 × 10−5 at t = 0.64, the CPU times (s), and their convergence orders with respect to Δt. Methods

Strang BAB

New BAB

Δt

h=

1 4096

CPU time

h=

1 16384

CPU time

Order

10−1

0.213107

−

0.04

0.708047

0.060967

2.248578 ·

0.01

0.708033

0.203797

1.422110 · 10−2

0.782239

1.991454

0.0025

0.708056

0.587193

8.888276 · 10−4

3.050981

1.999993

10−5

11.318841

2.000174

0.000625

0.708044

2.118856

5.553831 ·

0.16

0.708109

0.011110

1.776743 · 10−2

0.071795

−

0.052502

6.915551 ·

10−5

0.197560

4.002568

10−7

0.894927

3.966176

2.931852

2.065649

0.04

0.708056

0.01

0.708033

0.135650

2.831070 ·

0.0025

0.708056

0.840875

1.615497 · 10−8

h=1/1024

−2

h=1/1024

−5

10

10

−6

10

−3

10

−7

10 −4

Error

Error

10

−8

10

−5

10

−9

ε=0.000625 ε=0.0025 ε=0.01 ε=0.04 slope 2

−6

10

−7

10

10

ε=0.000625 ε=0.0025 ε=0.01 ε=0.04 slope 4

−10

10

−11

−2

−1

10

10 time step

10

−2

−1

10

10 time step

Fig. 5.2. The errors nεr − nε,h,Δt l2 of the position density produced by different methods for (5.1)-(5.2) with different ε at t = 0.64. Left: Strang BAB splitting spectral method; Right: New BAB splitting spectral method.

since uεr − uε,h,Δt l2 ≤ uεr − uε (t)l2 + uε (t) − uε,h,Δt l2 , the accuracy of our new methods with Δt = 0.0025 is limited by the accuracy of the reference solution uεr − uε (t)l2 , which seems to be around 10−9 . To illustrate the effect of the Planck constant ε on the errors, we list the orders of these numerical methods with respect to ε in Table 5.6. From Table 5.6, we observe that the errors of the Strang splitting spectral methods scale as 1/ε for any fixed Δt, while the errors of our new methods scale as 1/ε only when Δt2 ≥ Cε. When Δt2 = o(ε), it become evident that our schemes are uniformly convergent, which is independent of the Planck constant ε. 6. Concluding remarks. The computation of the semiclassical Schr¨odinger equation presents major challenges because of the presence of the small parameter ε. In this paper we proposed two exponential splitting scheme based on the symmetric BCH formula for time disctrization of semiclassical Schr¨odinger equation. Theoretical analysis and numerical experiments reveal that the new exponential splitting spectral methods have much higher accuracy and order of convergence than the Strang splitting spectral methods for any small parameter ε when they are applied to a linear

18

W. WANG AND J. TANG Table 5.6 The orders of different methods with respect to ε, where t = 0.64 and h = Methods

Strang BAB

New BAB

1 . 1024

ε

Δt = 0.16

Δt = 0.04

Δt = 0.01

Δt = 0.0025

0.04

−

−

−

−

0.01

−0.815363

−0.815200

−0.815184

−0.815167

0.0025

−0.983874

−0.985183

−0.985185

−0.985169

0.000625

−0.977616

−0.998973

−0.999056

−0.999057

0.04

−

−

−

−

0.01

−0.948138

−0.943555

−0.390194

−0.434190

0.0025

−0.996562

−0.996824

−0.331681

−0.106306

0.000625

−0.999784

−0.999758

−0.767119

−0.018428

problem with a quadratic potential. We only considered a linear time-independent potential semiclassical Schr¨odinger equation to which our new time-splitting spectral methods apply. It is natural to ask whether the new time-splitting spectral methods can be applied to the (nonautonomous) semiclassical Schr¨ odinger equation with time-dependent potential V (x, t) and nonlinear semiclassical Schr¨odinger equation. These will be matters for further research. Acknowledgments. The first author thanks Arieh Iserles for suggesting this problem and for stimulating and inspiring discussions during the author’s visit to the Department of Applied Mathematics and Mathematical Physics, University of Cambridge, as a postdoctoral research fellow. He also thanks Karolina Kropielnicha and Pranav Singh for helpful discussions. The authors would like to thank the associate editor and the anonymous referees for the valuable comments that lead to great improvements in the presentation of this paper. REFERENCES [1] W. Auzinger, T. Kassebacher, O. Koch and M. Thalhammer, Adaptive splitting methods for nonlinear Schr¨ odinger equations in the semiclassical regime, Numer. Algor., 72 (2016), pp. 1-35. [2] P. Bader, A. Iserles, K. Kropielnicka and P. Singh, Efficient approximation for the semiclassical Schr¨ odinger equation, FoCM, 14 (2014), pp. 689-720. [3] P. Bader, A. Iserles, K. Kropielnicka and P. Singh, Efficient methods for linear schr¨ odinger equation in the semiclassical regime with time-dependent potential, Proc. R. Soc. A., 472: 20150733, https://doi.org/10.1098/rspa.2015.0733 [4] W. Bao, S. Jin and P. A. Markowich, On time-splitting spectral approximations for the Schr¨ odinger equation in the semiclassical regime, J. Comput. Phys., 175 (2002), pp. 487– 524. [5] W. Bao, S. Jin and P. A. Markowich, Numerical study of time-splitting spectral discretizations of nonlinear Schr¨ odinger equations in the semiclassical regimes, SIAM J. Sci. Comput., 25 (2003), pp. 27-64. [6] W. Bao, Y. Cai and J. Yin, Uniform error bounds of time-splitting methods for the nonlinear Dirac equation in the nonrelativistic limit regime, arXiv:1906.11101v1. ´hats, An asymptotic preserving scheme based on a new [7] C. Besse, R. Carles and F. Me formulation for NLS in the semiclassical limit, Multiscale Model. Simul., 11 (2013), pp. 1228–1260.

¨ ESSP FOR SEMICLASSICAL SCHRODINGER EQUATION

19

[8] S. Boscarino, L. Pareschi and G. Russo, Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, SIAM J. Sci. Comput., 35 (2013), pp. A22–A51. [9] S. Boscarino, L. Pareschi and G. Russo, A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation, SIAM J. Numer. Anal., 55 (2017), pp. 2085-2109. [10] R. Carles, Semi-classical analysis for nonlinear Schr¨ odinger equations, World Scientific, New Jersey, 2008. [11] R. Carles, On Fourier time-splitting methods for nonlinear Schr¨ odinger equations in the semiclassical limit, SIAM J. Numer. Anal., 51 (2013), pp. 3232-3258. [12] R. Carles and C. Gallo, On Fourier time-splitting methods for nonlinear Schr¨ odinger equations in the semi-classical limit II. Analytic regularity, Numer. Math., 136 (2017), pp. 315-342. [13] F. Casas and A. Murua, An efficient algorithm for computing the Baker-Campbell-Hausdorff series and some of its applications, J. Math. Phys., 50 (2009), electronic. ´hats, An efficient algorithm for computing the Baker[14] P. Chartier, L. Le Treust and F. Me Campbell-Hausdorff series and some of its applications, J. Math. Phys., 50 (2009), electronic. ´hats, Uniformly accurate time-splitting methods for the [15] P. Chartier, L. Le Treust and F. Me semiclassical Schr¨ odinger equation. Part 1: Construction of the schemes and simulations, arXiv:1605.03446v1, (2016). ´hats, Uniformly accurate time-splitting methods for [16] P. Chartier, L. Le Treust and F. Me the semiclassical Schr¨ odinger equation. Part 2: Numerical analysis of the linear case, arXiv:1601.04825v1, (2016). [17] J. Chen, S. Li and Z. Zhang, Efficient multiscale methods for the semiclassical Schr¨ odinger equation with time-dependent potentials, arXiv: 1909.07203v1. [18] S. Descombes and M. Thalhammer, An exact local error representation of exponential operator splitting methods for evolutionary problems and applications to linear Schr¨ odinger equations in the semiclassical regime, BIT, 50 (2010), pp. 729–749. [19] S. Descombes and M. Thalhammer, The Lie-Trotter splitting for nonlinear evolutionary problems with critical parameters: a compact local error representation and application to nonlinear Schr¨ odinger equations in the semiclassical regime, IMA J. Numer. Anal., 33 (2013), pp. 722–745. [20] S. Duo and Y. Zhang, Mass-conservative Fourier spectral methods for solving the fraction nonlinear Schr¨ odinger equation, Comput. Math. Appl., 71 (2016), pp. 2257-2271. odinger equations, Zurich Lectures in Ad[21] E. Faou, Geometric numerical integration and Schr¨ vanced Mathematics, Europ. Math. Soc., Z¨ urich, 2012. ´rard, P. A. Markowich, N. J. Mauser, and F. Poupaud, Homogenization limits and [22] P. Ge Wigner transforms, Comm. Pure Appl. Math., 50 (1997), pp. 323-379. [23] D. J. Griffiths, Introduction to quantum mechnics, 2nd edn, Prentice Hall, Upper Saddle River, NJ, 2004. [24] E. Hairer, C. Lubich and G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, 2nd edn, Springer Verlag, Berlin, 2006. [25] A. Iserles, K. Kropielnicka and P. Singh, Commutator-free Magnus-Lanczos methods for the linear Schr¨ odinger equation, SIAM J Num. Anal. 56 (2018), pp. 1547-1569. [26] A. Iserles, K. Kropielnicka and P. Singh, Magnus-Lanczos methods with simplified commutators for the Schr¨ odinger equation with a time-dependent potential, J. Comput. Phys. 376 (2018), pp. 564-584. [27] S. Jin, P. A. Markowich and C. Sparber, Mathematical and computational methods for semicalssical Schr¨ odinger equations, Acta Numer., 20 (2011), pp. 121-209. [28] I. Kyza, C. Makridakis and M. Plexousakis, Error control for time-splitting spectral approximations of the semiclassical Schr¨ odinger equation, IMA J. Numer. Anal., 31 (2011), pp. 416-441. [29] C. Lubich, From quantum to classical molecular dynamics: reduced models and numerical analysis, Zurich Lectures in Advanced Mathematics, Europ. Math. Soc., Z¨ urich, 2008. [30] P. A. Markowich, P. Pietra and C. Pohl, Numerical approximation of quadratic observables of Schr¨ odinger-type equations in the semi-classical limit, Numer. Math., 81 (1999), 595630. [31] P. A. Markowich, P. Pietra, C. Pohl and H. Stimming, A Wigner-measure analysis of the Dufort-Frankel scheme for the Schr¨ odinger equation, SIAM J. Numer. Anal., 40 (2002) pp. 1281-1310. [32] R. I. McLachlan and G. R. W. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–434.

20

W. WANG AND J. TANG

[33] J. E. Pasciak, Spectral and pseudo-spectral methods for advection equations, Math. Comput., 35 (1980), pp. 1081–1092. [34] J. Shen, T. Tang and L. L. Wang, Spectral methods: algorithms, analyses and applications, Springer, Berlin, 2011. [35] M. Thalhammer, High-order exponential operator splitting methods for time-dependent Schr¨ odinger equations, SIAM J. Numer. Anal., 46 (2008), pp. 2022-2038. [36] M. Thalhammer, Convergence analysise of high-order time-splitting pseudospectral methods for nonlinear Schr¨ odinger equations, SIAM J. Numer. Anal., 50 (2012), pp. 3231-3258. [37] X. Li and A. Xiao, Time-splitting finite difference method with the wavelet-adaptive grids for semiclasseicalssical Gross-Pitaevskii equation in supercritical case, J. Comput. Phys., 267 (2014), pp. 146-161. [38] X. Yang and J. W. Zhang, Computation of the Schr¨ odinger equation in the semiclassical regime on an unbounded domain, SIAM J. Numer. Anal., 52 (2014), pp. 808-831. [39] H. Yo˘ sida, Construction of higher order symplectic integrators, Phys. Lett., 150 (1990), pp. 262–268.