Hermite ...

Found Comput Math DOI 10.1007/s10208-012-9124-x

Error Analysis of the Strang Time-Splitting Laguerre–Hermite/Hermite Collocation Methods for the Gross–Pitaevskii Equation Jie Shen · Zhong-Qing Wang

Received: 20 June 2011 / Revised: 12 April 2012 / Accepted: 17 May 2012 © SFoCM 2012

Abstract The aim of this paper is to carry out a rigorous error analysis for the Strang splitting Laguerre–Hermite/Hermite collocation methods for the time-dependent Gross–Pitaevskii equation (GPE). We derive error estimates for full discretizations of the three-dimensional GPE with cylindrical symmetry by the Strang splitting Laguerre–Hermite collocation method, and for the d-dimensional GPE by the Strang splitting Hermite collocation method. Keywords Error analysis · Strang splitting · Laguerre and Hermite collocation method · Gross–Pitaevskii equation · Nonlinear Schrödinger equation Mathematics Subject Classification 65D32 · 65N35 · 35J25

Communicated by Arieh Iserles. J. Shen School of Mathematical Science, Xiamen University, Xiamen, China J. Shen () Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA e-mail: [email protected] Z.-Q. Wang Department of Mathematics, Shanghai Normal University, Shanghai, 200234, PR China e-mail: [email protected] Z.-Q. Wang Division of Computational Science of E-institute of Shanghai Universities, Shanghai, PR China

Found Comput Math

1 Introduction In this paper we consider the error analysis of fully discretized schemes for the Gross–Pitaevskii equation (GPE). The GPE, which is a nonlinear Schrödinger equation, describes Bose–Einstein condensates (BECs) in the low temperature regime (cf. [9, 19]): i∂t ψ(x, t) = −

2  2 ψ(x, t) + V (x)ψ(x, t) + N U0 ψ(x, t) ψ(x, t), 2m

(1.1)

where ψ is the condensate wave function, m is the atomic mass,  is the Planck constant, N is the number of atoms in the condensate, and V (x) is an external trapping potential. When a harmonic trap potential is considered, V (x) = m2 (ωx2 x 2 + ωy2 y 2 + ωz2 z2 ), where ωx , ωy , and ωz are the trap frequencies in the x-, y-, and z-directions, respectively. In most current experiments, the traps are cylindrically symmetric, i.e., 2 ωx = ωy . U0 = 4πm as describes the interaction between atoms in the condensate with the s-wave scattering length as (positive for repulsive interaction and negative for attractive interaction). Using the normalization    ψ(x, t)2 dx = 1, (1.2) R3

and denoting V (x) =

 1 2 2 γx x + γy2 y 2 + γz2 z2 , 2

ωm = min{ωx , ωy , ωz },

γα =

ωα , α = x, y, z, ωm

4πas N β=√ , /mωm

we arrive at the following dimensionless GPE:  i∂t ψ(x, t) = − 12 ψ(x, t) + V (x)ψ(x, t) + β|ψ(x, t)|2 ψ(x, t), ψ(x, 0) = ψ0 (x),

lim|x|→∞ ψ(x, t) = 0,

t ≥ 0,

(1.3)

which is in fact a nonlinear Schrödinger equation. Much attention has been devoted to numerical approximation of the timedependent GPE (1.3). For instance, Bao, Jaksch and Markowich [2] and Bao and Shen [1] proposed several versions of time-splitting spectral methods, Ruprecht et al. [20] used the Crank–Nicolson finite difference method, and Cerimele et al. [5] proposed a particle-inspired scheme. However, the convergence analysis of semidiscretized Strang-type splitting schemes for linear and nonlinear Schrödinger equations only became available recently. Jahnke and Lubich [16] first presented an error bound for linear Schrödinger equations; then Lubich [18] gave an error bound for nonlinear Schrödinger equations. For related analyses in this direction, we refer to [6, 13, 17, 22]. However, to the authors’ best knowledge, not much is available for the fully discretized time-splitting schemes for nonlinear Schrödinger equations. The main reason is that, unlike the error analysis for fully discrete non-splitting schemes (e.g.,

Found Comput Math

backward Euler or Crank–Nicolson schemes), the error analysis for fully discrete time-splitting schemes is much more difficult than that for their semidiscrete counterparts. Recently, Gauckler [8] performed an error analysis for a Hermite collocation Strang splitting method for the d-dimensional GPE. The aim of this paper is to carry out an error analysis for the Hermite–Laguerre collocation Strang splitting method for the three-dimensional GPE with cylindrical symmetry, and for the Hermite spectral Strang splitting method for the d-dimensional GPE. Our main contributions are twofold: (i) our error analysis for the Hermite– Laguerre collocation method, which was actually implemented in [1], is new; (ii) our results for the Hermite collocation Strang splitting method for the d-dimensional GPE significantly improve the error estimates presented in [8]. Moreover, while our analysis for the semidiscretization in time is similar to those in [8, 18], our analysis for the full discretization has some essentially different components from those in [8], and leads to improved error estimates. Nevertheless, the techniques developed in [8, 18] have been very useful for our analysis. More precisely, we first focus on the special case of (1.3) with cylindrical symmetry, i.e., γx = γy = γr and ψ0 (x, y, z) = ψ0 (r, z). Then the equation becomes: ⎧ i∂t ψ(r, z, t) = [− 2r1 ∂r (r∂r ψ(r, z, t)) + 12 γr2 r 2 ψ(r, z, t)] ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + [− 12 ∂z2 ψ(r, z, t) + 12 γz2 z2 ψ(r, z, t)] ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

r ≥ 0, z ∈ R,

+ β|ψ(r, z, t)|2 ψ(r, z, t) =: Ar ψ + Bz ψ + N (ψ),

ψ(r, z, 0) = ψ0 (r, z),

limr,|z|→∞ ψ(r, z, t) = 0,

(1.4)

t ≥ 0.

We present a full discretization scheme for (1.4) by using a Strang splitting scheme in time and a Laguerre–Hermite collocation method in space. The main results for (1.4) are summarized in Theorem 5.1 and Corollary 5.1. Note that the analysis for the axisymmetric case is much more difficult than the usual d-dimensional case due in part to the involvement of the Laguerre functions. We then consider the full discretization of the d-dimensional GPE directly by using the Strang splitting in time and a Hermite collocation method in space. The main results for the d-dimensional GPE are summarized in Theorem 6.3 and Corollary 6.1. The paper is organized as follows. In the next section, we derive some basic results for the Laguerre and Hermite approximations which will be used for the error analysis. In Sect. 3, we describe the semidiscrete Strang splitting scheme and the fully discrete Strang splitting Laguerre–Hermite collocation scheme for (1.4). The error analysis for the semidiscrete Strang splitting scheme for (1.4) is performed in Sect. 4, while that for the fully discrete Strang splitting Laguerre–Hermite collocation scheme is presented in Sect. 5. We consider the error analysis of the Strang splitting Hermite collocation scheme for the d-dimensional GPE in the last section. 2 Scaled Laguerre–Hermite–Gauss Interpolation In this section, we describe scaled Laguerre–Hermite–Gauss interpolation and derive some basic results which will be used later.

Found Comput Math

2.1 Scaled Laguerre Functions Let I = (0, ∞) and let Lm (ˆr ) be the Laguerre polynomial of degree m satisfying rˆ Lm (ˆr ) + (1 − rˆ )Lm (ˆr ) + mLm (ˆr ) = 0, 

Lm (ˆr )Ln (ˆr )e−ˆr dˆr = δmn ,

rˆ ∈ I, m ≥ 0,

m, n ≥ 0,

I

where δmn is the Kronecker delta function. For any positive integer N , we denote by PN the set of all algebraic polynomials of degree at most N . Let {ˆrj , ωˆ jr }N j =0 be the Laguerre–Gauss points and weights, and r ˆ IN : C(I ) → PN be the corresponding interpolation operator in the r-direction such that Iˆ Nr v(ˆrj ) = v(ˆrj ),

0 ≤ j ≤ N.

For any integer r ≥ 0, we define the weighted Sobolev space Hχr (I ) with the weight function χ in the usual way. In particular, L2χ (I ) = Hχ0 (I ). For any r > 0, we define the space Hχr (I ) by space interpolation as in [4]. According to Theorem 3.4 of [11], for any v ∈ Hω10 (I ) and ∂rˆ v ∈ L2ω1 (I ) with ω0 (ˆr ) = e−ˆr and ω1 (ˆr ) = rˆ e−ˆr , we have  r  Iˆ v  2 N L

1

ω0

(I )

≤ cN − 2 ∂rˆ v L2ω

0

(I )

1 + c(ln N ) 2 v L2ω

0

(I )

+ ∂rˆ v L2ω

1

(I )

 .

(2.1)

In order to determine the eigenfunctions of the linear operators Ar , Bao and Shen [1] introduced the change of variable r = γrˆr and the scaled Laguerre function

lm (r) =

γr −ˆr /2 e Lm (ˆr ) = π

γr −γr r 2 /2  2  e Lm γ r r , π

r ∈ I,

(2.2)

which satisfies Ar lm (r) := −

 1 1  ∂r r∂r lm (r) + γr2 r 2 lm (r) = μrm lm (r), 2r 2

(2.3)

μrm = γr (2m + 1), m ≥ 0,  lm (r)ln (r)r dr = δmn ,

2π

m, n ≥ 0.

I

We also denote  rˆj rj = , γr

ωjr =

π r rˆj ωˆ e , γr j

  r XN = span lm (r) : 0 ≤ m ≤ N ,

(2.4)

Found Comput Math

where rj and ωjr are the scaled Laguerre–Gauss points and weights, respectively. We recall that [1]: N 

ωjr lm (rj )ln (rj ) = δnm ,

∀0 ≤ n + m ≤ 2N + 1.

(2.5)

j =0 r , the corresponding interpolation operator, by Next, we define INr : C(I ) → XN

INr u(rj ) = u(rj ), For any u(r) = e−γr r

2 /2

0 ≤ j ≤ N.

v(γr r 2 ), we have

2

2

eγr rj /2 INr u(rj ) = eγr rj /2 u(rj ) = v(ˆrj ) = Iˆ Nr v(ˆrj ),

eγr r

INr u(r)|

r= γrˆr

∈ PN and Iˆ Nr v(ˆr ) ∈ PN . Hence,

 INr u(r)

r= γrˆr

= Iˆ Nr v(ˆr ).

2 /2

Furthermore, due to (2.2), eγr r

2 /2

0 ≤ j ≤ N.

This with (2.1) leads to        r 2   I u r dr ≤ cN −1 u(r)2 r dr + ∂r u(r)2 r −1 dr N I

I

 + c ln N ≤ cN −1  +



I

    u(r)2 r + r 3 dr +

I



  ∂r u(r)2 r dr



I

  ∂r u(r)2 r −1 dr + c ln N



I

    u(r)2 r + r 3 dr I

   ∂r u(r)2 r dr .

(2.6)

I

2.2 Scaled Hermite Functions Let Hl (z) be the standard Hermite polynomials satisfying Hl (z) − 2zHl (z) + 2lHl (z) = 0,  R

Hl (z)Hn (z)e−z dz = 2

z ∈ R, l ≥ 0,

√ l π2 l!δln ,

l, n ≥ 0.

We consider the scaled Hermite function (cf. [1]) 1

hl (z) = (γz /π) 4 e−γz z

2 /2

 √ Hl ( γz z)/ 2l l!,

z ∈ R,

(2.7)

which satisfies 1 1 Bz hl (z) := − hl (z) + γz2 z2 hl (z) = μzl hl (z), 2 2

μzl =

2l + 1 γz , l ≥ 0, 2

(2.8)

Found Comput Math

 R

hl (z)hn (z) dz = δln ,

l, n ≥ 0.

(2.9)

Next let {ˆzk , ωˆ kz }N k=0 be the Hermite–Gauss points and weights. Denote by zˆ k zk = √ , γz

ωˆ z 2 ωkz = √ k ezˆ k , γz

the scaled Hermite–Gauss points and weights, respectively. According to [1], N 

ωkz hm (zk )hn (zk ) = δnm ,

∀0 ≤ n + m ≤ 2N + 1.

(2.10)

k=0 z Let XN = span{hl (z) : 0 ≤ l ≤ N}. We define the scaled Hermite–Gauss interpolation z operator INz : C(R) → XN by

INz v(zk ) = v(zk ),

0 ≤ k ≤ N.

The following result is established in [10]:  z  I v  N

L2 (R)

  1 ≤ c v L2 (R) + N − 6 |v|H 1 (R) ,

v ∈ H 1 (R),

(2.11)

where | · |H 1 (R) denotes the seminorm of H 1 (R). We now introduce some properties about the scaled Hermite functions in R3 . Let hk (x) and hm (y) be the Hermite functions defined in (2.7) with γx and γy instead of γz , respectively. The corresponding operators in (2.8) are denoted by Bx and By , respectively, i.e., 1 1 Bx hk (x) := − hk (x) + γx2 x 2 hk (x) = μxk hk (x), 2 2 1 1 y By hm (y) := − hm (y) + γy2 y 2 hm (y) = μm hm (y), 2 2

μxk = y

2k + 1 γx , 2

μm =

2m + 1 γy . 2

(2.12)

(2.13)

We denote by L2 (R3 ) and H s (R3 ) (s > 0) the usual Sobolev spaces with the usual notation for their seminorms and norms. It can be easily shown that the linear differential operator Bx + By + Bz is positive definite and self-adjoint. Indeed, for any u and v in the domain of Bx + By + Bz , applying integration by parts leads to     (Bx + By + Bz )u, v R3 = u, (Bx + By + Bz )v R3 = a(u, v), (2.14)   (Bx + By + Bz )u, u R3 = a(u, u) > 0, if u = 0, where   1  2 2 1 γx x + γy2 y 2 + γz2 z2 u, v R3 . a(u, v) = (∇u, ∇v)R3 + 2 2 From (2.8), (2.9), (2.12), and (2.13), we have

Found Comput Math

  a hk (x)hm (y)hl (z), hk  (x)hm (y)hl  (z)   = (Bx + By + Bz )hk (x)hm (y)hl (z), hk  (x)hm (y)hl  (z) R3   y = μxk + μm + μzl δkk  δmm δll  .

(2.15)

Since Bx + By + Bz is a positive definite and self-adjoint operator in L2 (R3 ), the fractional power (Bx + By + Bz )1/2 is well defined, and the associated norms can be characterized by (see, e.g., [21, 23]):   (Bx + By + Bz )1/2 u2 2

L (R3 )

  (Bx + By + Bz )m+1/2 u2 2

L (R3 )

= a(u, u),

(2.16)

  = a (Bx + By + Bz )m u, (Bx + By + Bz )m u ,

∀m ∈ N. (2.17) For any u ∈ L2 (R3 ), we write (cf. [8]) u(x, y, z) =

∞ 

ukml hk (x)hm (y)hl (z).

k,m,l=0

We introduce the following three Sobolev spaces equipped with the norms:  u H s (R3 ) = A

∞   x s y μk + μm + μzl |ukml |2

1 2

,

k,m,l=0

 s  u H s (R3 ) = (Bx + By + Bz ) 2 uL2 (R3 ) , B

 u H s (R3 ) = C

s   2  s−k−m−l k m l 2 2  x + y 2 + z2 + 1 ∂x ∂y ∂z uL2 (R3 )

1 2

.

k+m+l=0

With a slight modification of Lemma 2.1 in [25], we can prove the following. Lemma 2.1 The previous three norms are equivalent, i.e., u H s (R3 ) = u H s (R3 ) ∼ u H s (R3 ) . A

B

C

Proof Let integer r ≥ 0. According to (2.8), (2.9), (2.12), and (2.13), we have that for s = 2r,   u 2H s (R3 ) = (Bx + By + Bz )r u, (Bx + By + Bz )r u R3 B  ∞   r y μxk + μm + μzl ukml hk hm hl , = k,m,l=0

Found Comput Math ∞   x r y μk  + μm + μzl uk  m l  hk  hm hl 



k  ,m ,l  =0

R3

∞   x 2r y μk + μm + μzl |ukml |2 = u 2H s (R3 ) .

=

A

k,m,l=0

Next, by (2.15) and (2.17), we deduce that for s = 2r + 1,   u 2H s (R3 ) = a (Bx + By + Bz )r u, (Bx + By + Bz )r u B  ∞   r y =a μxk + μm + μzl ukml hk hm hl , k,m,l=0 ∞   x r y μk  + μm + μzl uk  m l  hk  hm hl 



k  ,m ,l  =0 ∞ 

=

 x r  r y y μk + μm + μzl μxk + μm + μzl

k,m,l,k  ,m ,l  =0

× ukml uk  m l  a(hk hm hl , hk  hm hl  ) =

∞   x 2r+1 y μk + μm + μzl |ukml |2 = u 2H s (R3 ) . A

k,m,l=0

The above two estimates, together with function space interpolation as in [4], lead to the desired result u H s (R3 ) = u H s (R3 ) . Furthermore, following [25] we can verify A B readily that u H s (R3 ) ∼ u H s (R3 ) . A

C



This completes the proof. Remark 2.1 We note that Helffer [14] proved the following equivalence result:   s  s  u L2 (R3 ) + (Bx + By ) 2 uL2 (R3 ) + Bz2 uL2 (R3 )    s  ∼ u H s (R3 ) +  x 2 + y 2 2 uL2 (R3 ) + zs uL2 (R3 ) .

The above result, although very similar to Lemma 2.1, is nevertheless different. Furthermore, Abdallah, Castella, and Méhats [3] extended the above result to the more general case with nonharmonic oscillator. Lemma 2.2 We have the following inequalities: uvw L2 (R3 ) ≤ c u H 1 (R3 ) v H 1 (R3 ) w H 1 (R3 ) ,

(2.18)

uvw L2 (R3 ) ≤ c u L2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ,

(2.19)

Found Comput Math

uvw H 1 (R3 ) ≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ,

(2.20)

uvw H k (R3 ) ≤ c u H k (R3 ) v H k (R3 ) w H k (R3 ) ,

∀ integer k ≥ 2. (2.21)

Proof The results (2.18)–(2.20) and (2.21) with k = 2 are established in [18]. By induction, we can obtain the result (2.21) with integer k ≥ 3.  The previous results can be extended to the corresponding ones in HAs (R3 ), as stated below. Lemma 2.3 The following inequalities hold: uvw L2 (R3 ) ≤ c u H 1 (R3 ) v H 1 (R3 ) w H 1 (R3 ) .

(2.22)

uvw L2 (R3 ) ≤ c u L2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) .

(2.23)

uvw H 1 (R3 ) ≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) .

(2.24)

uvw H k (R3 ) ≤ c u H k (R3 ) v H k (R3 ) w H k (R3 ) ,

∀ integer k ≥ 2. (2.25)

A

A

A

A

A

A

A

A

A

A

A

A

A

Proof Cases 1 and 2. Since u H s (R3 ) ∼ u H s (R3 ) , we have u H s (R3 ) ≤ A C c u H s (R3 ) , s ≥ 0. Hence by (2.18) and (2.19) we get the results (2.22) and (2.23). A Case 3. According to the equivalence of the norms, we derive readily that   1 uvw H 1 (R3 ) ≤ c uvw H 1 (R3 ) ≤ c|uvw|H 1 (R3 ) +c x 2 +y 2 +z2 +1 2 uvw L2 (R3 ) . A

C

From (2.20) we have |uvw|H 1 (R3 ) ≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) . A

A

A

By (2.19),   2 1  x + y 2 + z2 + 1 2 uvw 

L2 (R3 )

 1  ≤ c x 2 + y 2 + z2 + 1 2 uL2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) C

C

C

≤ c u H 1 (R3 ) v H 2 (R3 ) w H 2 (R3 ) . A

A

(2.26)

A

A combination of the previous three inequalities leads to the desired result (2.24). Case 4. Obviously, uvw H 2 (R3 ) ≤ c|uvw|H 2 (R3 ) A

+c

   1  x 2 + y 2 + z2 + 1 2 ∂ k ∂ m ∂ l (uvw) x y

k+m+l=1

   + c x 2 + y 2 + z2 + 1 uvw L2 (R3 ) .

z

L2 (R3 )

Found Comput Math

From (2.21), |uvw|H 2 (R3 ) ≤ c u H 2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ≤ c u H 2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) . A

A

A

1 2

Since (x 2 + y 2 + z2 + 1) ∂x u L2 (R3 ) ≤ u H 2 (R3 ) , we can use an argument similar C to (2.26) to derive that    1  x 2 + y 2 + z2 + 1 2 ∂ k ∂ m ∂ l (uvw) 2 3 x y z L (R ) k+m+l=1

≤ c u H 2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) . A

A

A

Furthermore, by (2.19),   2   x + y 2 + z2 + 1 uvw  2 3 L (R )  2   2 2  ≤ c x + y + z + 1 uL2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) ≤ c u H 2 (R3 ) v H 2 (R3 ) w H 2 (R3 ) . A

A

A

Therefore, a combination of the previous four inequalities leads to (2.25) with k = 2. We can obtain the desired results for integer k ≥ 3 by using an argument as in the proof of (6.14) of this paper.  2.3 Approximation by the Mixed Laguerre–Hermite Functions Set Ω = I × R. In order to present the convergence of the three-dimensional GPE with cylindrical symmetry, we need some approximation results on the mixed Laguerre–Hermite functions. To this end, we define the inner product and norm of L2 (Ω) with complex-valued functions by  1 2 u(r, z)v(r, z)r dr dz, v L2 (Ω) = (v, v)Ω . (u, v)Ω = 2π Ω

We notice that the inner product introduced here is not the usual inner product on L2 (I × R), but on L2 (R3 ) using cylindrical coordinates. For any u(r, z) ∈ L2 (Ω), we write ∞  u(r, z) = uml lm (r)hl (z). m,l=0

Obviously, the linear differential operator Ar + Bz is positive definite and selfadjoint. Thus for any u and v in the domain of Ar + Bz ,     (Ar + Bz )u, v Ω = u, (Ar + Bz )v Ω = b(u, v), (2.27)   (Ar + Bz )u, u Ω = b(u, u) > 0, if u = 0, where the bilinear form   1 1  2 2 b(u, v) = (∇u, ∇v)Ω + γr r + γz2 z2 u, v Ω . 2 2

Found Comput Math

Moreover, for any functions lm (r)hl (z) and lm (r)hl  (z),     b(lm hl , lm hl  ) = (Ar + Bz )lm hl , lm hl  Ω = μrm + μzl δmm δll  .

(2.28)

The fractional power (Ar + Bz )1/2 is also well defined, and the associated norms can be characterized by   (Ar + Bz )1/2 u2 = b(u, u), (2.29) Ω     (Ar + Bz )m+1/2 u2 = b (Ar + Bz )m u, (Ar + Bz )m u , Ω

∀m ∈ N.

(2.30)

We next introduce two Sobolev spaces equipped with the norms: 

∞   r s μm + μzl |uml |2

u HAs (Ω) =

1 2

,

m,l=0

 s  u HBs (Ω) = (Ar + Bz ) 2 uL2 (Ω) . It is also easy to verify that u HAs (Ω) = u HBs (Ω) . Let XN := XN (γr , γz ) = span{lm (r)hl (z) : 0 ≤ m, l ≤ N }. According to (2.3) and (2.8), {lm (r)hl (z)} are the eigenfunctions of the operator Ar + Bz with the eigenvalues μrm + μzl . Lemma 2.4 For any φ ∈ XN and s ≥ 0, s

φ HAs (Ω) ≤ cN 2 φ L2 (Ω) . Proof Given φ ∈ XN , we write φ(r, z) =

N 

φml lm (r)hl (z).

m,l=0

For any integer s ≥ 0, N   r s   φ 2H s (Ω) = (Ar + Bz )s φ, φ Ω = μm + μzl |φml |2 A

m,l=0

≤ cN s

N 

|φml |2 = cN s φ 2L2 (Ω) .

m,l=0

This with a standard space interpolation technique [4] yields the desired result.



We now consider the orthogonal projection. For any u ∈ L2 (Ω), the orthogonal projection operator PN : L2 (Ω) → XN is defined by (u − PN u, φ)Ω = 0,

∀φ ∈ XN .

(2.31)

Found Comput Math

In particular, if u ∈ HBs (Ω) with integer s ≥ 0, we have  s s  (Ar + Bz ) 2 (u − PN u), (Ar + Bz ) 2 φ Ω   = u − PN u, (Ar + Bz )s φ Ω = 0, ∀φ ∈ XN ,

(2.32)

which means that the L2 (Ω)-orthogonal projection operator PN is also the HBs (Ω)orthogonal projection operator. Theorem 2.1 If u ∈ HAs (Ω), then for any 0 ≤ μ ≤ s, μ−s 2

u − PN u H μ (Ω) ≤ cN A

u HAs (Ω) .

Proof For any integers μ and s with 0 ≤ μ ≤ s, u − PN u 2H μ (Ω) = A

∞ N ∞     r μ  r μ μm + μzl |uml |2 + μm + μzl |uml |2 m,l=N +1

m=0 l=N +1

∞ N    r μ μm + μzl |uml |2

+

m=N +1 l=0 ∞   r s μm + μzl |uml |2 = cN μ−s u 2H s (Ω) .

≤ cN μ−s

A

(2.33)

m,l=0

This with a standard space interpolation technique leads to the desired result.



We are now in position to study the interpolation operator. The scaled Laguerre– Hermite–Gauss interpolant IN : C(Ω) → XN is determined by IN u(rj , zk ) = u(rj , zk ),

0 ≤ j, k ≤ N.

Clearly, IN u = INr INz u. Hence by (2.6), (2.11), and the equivalence of the norms, a direct calculation shows that    1 IN u 2L2 (Ω) ≤ cN −1 |∂r u|2 r −1 dr dz + N − 3 |∂z ∂r u|2 r −1 dr dz Ω

Ω



  1 |u|2 r + r 3 dr dz + N − 3

+ c ln N 

Ω

|∂r u| r dr dz + N

+

2

Ω

− 13

|∂z ∂r u| r dr dz 2

  1 ≤ cN −1 u 2H 1 (R3 ) + N − 3 u 2H 2 (R3 ) C

  1 + c ln N u 2H 1 (R3 ) + N − 3 u 2H 2 (R3 ) C

  |∂z u|2 r + r 3 dr dz

Ω

 Ω

C



C



Found Comput Math

  1 ≤ c ln N u 2H 1 (Ω) + N − 3 u 2H 2 (Ω) . A

(2.34)

A

Theorem 2.2 If u ∈ HAs (Ω), then for any 0 ≤ μ ≤ s and s ≥ 2, 1

5

u − IN u H μ (Ω) ≤ c(ln N ) 2 N 6 +

μ−s 2

A

u HAs (Ω) .

Proof From Theorem 2.1 and Lemma 2.4, for any 0 ≤ μ ≤ s,   u − IN u H μ (Ω) ≤ u − PN u H μ (Ω) + IN (u − PN u)H μ (Ω) A

A

≤ cN

μ−s 2

A

  u HAs (Ω) + cN IN (u − PN u)L2 (Ω) . μ 2

Moreover, by (2.34) and Theorem 2.1, for any s ≥ 2,    1 1−s 5 s IN (u − PN u) 2 ≤ c(ln N ) 2 N 2 u HAs (Ω) + N 6 − 2 u HAs (Ω) L (Ω) 1

5

≤ c(ln N ) 2 N 6 − 2 u HAs (Ω) . s

Therefore, 1

5

u − IN u H μ (Ω) ≤ c(ln N ) 2 N 6 + A

μ−s 2

u HAs (Ω) .



3 A Time-Splitting Laguerre–Hermite Collocation Method We now describe the time-splitting spectral method in [1] for the three-dimensional (3D) Gross–Pitaevskii equation (GPE) with cylindrical symmetry (1.4). For simplicity, we shall only consider the second-order Strang splitting scheme. It is expected that the technique presented in this paper will eventually enable us to prove error estimates for the fourth-order splitting scheme used in [1]. 3.1 Strang Splitting in Time For the semidiscretization in time, we split the 3D GPE with cylindrical symmetry (1.4) into its linear and nonlinear parts:    1 1 1 ∂r (r∂r ψ) + ∂z2 ψ + γr2 r 2 + γz2 z2 ψ, (3.1) i∂t ψ(r, z, t) = (Ar + Bz )ψ = − 2 r 2 2  i∂t ψ(r, z, t) = β ψ(r, z, t) ψ(r, z, t).

(3.2)

Equations (3.1) and (3.2) are exactly solvable since |ψ| is invariant in time along the solution of (3.2). For a given time step τ > 0, let tn = nτ , n = 0, 1, . . . , and let ψ n be the approximation of ψ(tn ). Then, the second-order Strang splitting in time for (1.4) is as follows:   −i τ2 (Ar +Bz ) n 2 τ ψ | −i τ2 (Ar +Bz ) n ψ n+1 = Φ τ ψ n := e−i 2 (Ar +Bz ) e−iτβ|e e ψ , (3.3) where ψ 0 = ψ0 .

Found Comput Math

3.2 Full Discretization Before we present the full discretization of (1.4), let us describe the semidiscretization in space using the Laguerre–Hermite collocation method. Find ψN (r, z, t) ∈ XN , i.e., N 

ψN (r, z, t) =

ψml (t)lm (r)hl (z),

(3.4)

m,l=0

such that  2 i∂t ψN (rj , zk , t) = (Ar + Bz )ψN (rj , zk , t) + β ψN (rj , zk , t) ψN (rj , zk , t), (3.5) ∀0 ≤ j, k ≤ N, where ψN (rj , zk , 0) = ψ0 (rj , zk ). The system (3.5) can also be rewritten as 

i∂t ψN (r, z, t) = (Ar + Bz )ψN (r, z, t) + βIN (|ψN (r, z, t)|2 ψN (r, z, t)), ψN (r, z, 0) = IN ψ0 (r, z).

(3.6)

We now combine the semidiscretization in time in Sect. 3.1 with the semidiscretization in space to obtain a full discretization of (1.4). To do this, we split the space-discretized equation (3.6) into its linear and nonlinear parts: i∂t ψN (r, z, t) = (Ar + Bz )ψN (r, z, t),

(3.7)

2   i∂t ψN (r, z, t) = βIN ψN (r, z, t) ψN (r, z, t) .

(3.8)

Clearly, |ψN (rj , zk , t)| is conserved in time. Thus, the fully discrete Laguerre– Hermite Strang time-splitting scheme is as follows: ψN0 (r, z) = IN ψ0 (r, z);  n   −i τ2 (Ar +Bz ) n 2 τ τ ψN | −i τ2 (Ar +Bz ) n ψNn+1 = ΦN ψN := e−i 2 (Ar +Bz ) IN e−iτβ|e e ψN , (3.9) n ≥ 0. Define the discrete inner product u, vN =

N 

ωjr ωkz u(rj , zk )v(rj , zk ).

j,k=0

According to (2.4), (2.5), (2.9), and (2.10), we have φ, ψN = (φ, ψ)Ω ,

∀φψ ∈ X2N +1 .

(3.10)

Found Comput Math

Then, we can rewrite (3.9) in a more computationally friendly algorithm: Given {ψNn (rj , xk )}, compute N 

(1)

ψN (rj , zk ) =

ml lm (rj )hl (zk ), e−i 2 (μm +μl ) U τ

r

z

m,l=0 (1)

ψN (rj , zk ) = e−iτβ|ψN (2)

N 

ψNn+1 (rj , zk ) =

(rj ,zk )|2

(1)

ψN (rj , zk ),

(3.11)

ml lm (rj )hl (zk ), e−i 2 (μm +μl ) V τ

r

z

m,l=0

where N    ml = ψ n , lm hl = ωjr ωkz ψNn (rj , zk )lm (rj )hl (zk ), U N N j,k=0 N    (2) ml = ψ (2) , lm hl = V ωjr ωkz ψN (rj , zk )lm (rj )hl (zk ). N N j,k=0

4 Error Analysis for the Semidiscrete Strang Splitting Scheme In this section, we shall derive error bounds for the semidiscretization scheme (3.3). We follow the basic procedure in [18], and generalize the error estimates to the HAk (Ω)-norms. We start by establishing a lemma which is needed for dealing with nonlinear terms. Lemma 4.1 The following inequalities hold: uvw L2 (Ω) ≤ c u H 1 (Ω) v H 1 (Ω) w H 1 (Ω) ,

(4.1)

uvw L2 (Ω) ≤ c u L2 (Ω) v H 2 (Ω) w H 2 (Ω) ,

(4.2)

uvw H 1 (Ω) ≤ c u H 1 (Ω) v H 2 (Ω) w H 2 (Ω) ,

(4.3)

A

A

A

A

A

A

A

A

A

uvw H k (Ω) ≤ c u H k (Ω) v H k (Ω) w H k (Ω) , A

A

A

A

∀ integer k ≥ 2.

(4.4)

Proof For any function  u(x, y, z) in R3 with cylindrical symmetry, we denote u(r, z) :=  u(x, y, z), (r, z) ∈ Ω. Clearly, in this case, γx = γy = γr . Thereby, u(x, y, z) = (Bx + By )

 1 2 −∂x − ∂y2 + γx2 x 2 + γy2 y 2  u(x, y, z) = Ar u(r, z), 2

and a( u, u) = b(u, u).

Found Comput Math

Accordingly, we have  u H s (R3 ) = u HBs (Ω) , B

which implies  u H s (R3 ) = u HAs (Ω) . A

We can then obtain the desired results from the above and Lemma 2.3.



Lemma 4.2 (Stability) If ψ, ϕ ∈ HA2 (Ω) ∩ HAk (Ω) with integer k ≥ 0, then cτ ( ψ 2 j + ϕ 2 j )   τ HA (Ω) HA (Ω) Φ (ψ) − Φ τ (ϕ) k ≤ e ψ − ϕ H k (Ω) , H (Ω) A

A

where j = max(k, 2). Proof We first consider the cases with 0 ≤ k ≤ 2. It is clear that the operator τ e−i 2 (Ar +Bz ) preserves the norm · HAs (Ω) for any integer s ≥ 0, since e−i 2 (Ar +Bz ) lm (r)hl (z) = e−i 2 (μm +μl ) lm (r)hl (z). τ

τ

z

r

Thereby, we only need to compare e−iτβ|ψ| ψ and e−iτβ|ϕ| ϕ, which are the solutions at time τ of the linear initial value problems: 2

2

i∂t θ = β|ψ|2 θ,

θ (0) = ψ,

(4.5)

i∂t η = β|ϕ|2 η,

η(0) = ϕ.

(4.6)

We first establish a bound for θ (t) H 2 (Ω) . By (4.5) we have A

      i (Ar + Bz )∂t θ, (Ar + Bz )θ Ω = β (Ar + Bz ) |ψ|2 θ , (Ar + Bz )θ Ω . Taking the imaginary part in the above, we obtain 2      ∂t θ (t)H 2 (Ω) ≤ 2β |ψ|2 θ (t)H 2 (Ω) θ (t)H 2 (Ω) , A

which implies that

A

A

    ∂t θ (t)H 2 (Ω) ≤ β |ψ|2 θ (t)H 2 (Ω) . A

A

On the other hand, by (4.4) we obtain     2 |ψ| θ (t) 2 ≤ c ψ 2H 2 (Ω) θ (t)H 2 (Ω) . H (Ω) A

A

A

A combination of the previous two inequalities leads to   θ (t)

HA2 (Ω)

 ≤ ψ H 2 (Ω) + c A

0

t

  ψ 2H 2 (Ω) θ (ξ )H 2 (Ω) dξ. A

A

Found Comput Math

Applying the usual Gronwall inequality, we obtain   θ (t)

≤e H 2 (Ω)

ct ψ 2

2 (Ω) HA

ψ H 2 (Ω) .

(4.7)

A

A

Next, by (4.5) and (4.6), we have i∂t (θ − η) = β(ψ − ϕ)ψθ + βϕ(ψ − ϕ)θ + β|ϕ|2 (θ − η).

(4.8)

We now proceed to treat different cases separately. (i) k = 0. By taking the inner product in (4.8) with θ − η, and then taking the imaginary part, we get   ∂t θ − η L2 (Ω) ≤ β (ψ − ϕ)ψθ + ϕ(ψ − ϕ)θ L2 (Ω) . The above with (4.2) yields   ∂t θ − η L2 (Ω) ≤ c ψ − ϕ L2 (Ω) θ H 2 (Ω) ψ H 2 (Ω) + ϕ H 2 (Ω) . A

A

A

Integrating the above from 0 to τ and using (4.7), we obtain   θ (τ ) − η(τ ) 2 ≤ ψ − ϕ L2 (Ω) L (Ω)   cτ ψ 2H 2 (Ω) A + cτ ψ − ϕ L2 (Ω) ψ 2H 2 (Ω) + ϕ 2H 2 (Ω) e A

A

cτ ψ 2 2    HA (Ω) ≤ 1 + cτ ψ 2H 2 (Ω) + ϕ 2H 2 (Ω) ψ − ϕ L2 (Ω) e A

≤e

cτ ( ψ 2

2 (Ω) HA

A

+ ϕ 2

2 (Ω) HA

)

ψ − ϕ L2 (Ω) .

(ii) k = 1, 2. By (4.8) we have   k k k i(Ar + Bz ) 2 ∂t (θ − η) = β(Ar + Bz ) 2 (ψ − ϕ)ψθ + β(Ar + Bz ) 2 ϕ(ψ − ϕ)θ  k + β(Ar + Bz ) 2 |ϕ|2 (θ − η) . (4.9) k

Take the inner product in (4.9) with (Ar + Bz ) 2 (θ − η). From (4.3) and (4.4), and using an argument similar to the case k = 0, we get   ∂t θ − η H k (Ω) ≤ c ψ − ϕ H k (Ω) θ H 2 (Ω) ψ H 2 (Ω) + ϕ H 2 (Ω) A

A

A

+ c ϕ 2H 2 (Ω) θ A

A

A

− η H k (Ω) . A

Therefore, by (4.7) and the Gronwall inequality, we obtain the desired result. (iii) k > 2. By a similar argument as before, we can establish the bound on θ (t) H k (Ω) , namely, A

  θ (t)

≤e H k (Ω) A

ct ψ 2

k (Ω) HA

ψ H k (Ω) , A

k ≥ 2.

(4.10)

Found Comput Math

Thus by (4.10), (4.4), and the Gronwall inequality, we obtain the result with integer k > 2.  We are now in position to estimate the local error. To this end, we denote (ψ) = −iβ|ψ|2 ψ. V

T(ψ) = −i(Ar + Bz )ψ,

(4.11)

Their Lie commutator (cf. [12, 15]) is as follows: ](ψ) = T (ψ)V (ψ) − V  (ψ)T(ψ) [T, V   = −β(Ar + Bz ) |ψ|2 ψ − βψ 2 (Ar + Bz )ψ + 2β|ψ|2 (Ar + Bz )ψ.

(4.12)

Lemma 4.3 For any ψ ∈ HAk+2 (Ω) with integer k ≥ 0, we have   [T, V ](ψ)

HAk (Ω)

≤ c ψ 3

HAk+2 (Ω)

(4.13)

.

If, in addition, ψ ∈ HAk+4 (Ω), then     T, [T, V ] (ψ)

HAk (Ω)

≤ c ψ 3

HAk+4 (Ω)

(4.14)

.

Proof Using (4.12) and (4.2)–(4.4), we obtain that   2   2   |ψ| ψ  k+2 ψ (Ar + Bz )ψ  k [T, V ](ψ) k ≤ c + c HA (Ω) HA (Ω) HA (Ω)  2  + c|ψ| (Ar + Bz )ψ  k HA (Ω)

≤ c ψ 3

HAk+2 (Ω)

(4.15)

.

Next, we derive a bound for the following commutator (cf. [12, 15]):     ] (ψ) = T [T, V ](ψ) − [T, V ] (ψ)T(ψ). T, [T, V By (4.11) and (4.15) we have     ](ψ)  k T [T, V

HA (Ω)

  ](ψ) k+2 = [T, V ≤ c ψ 3 H (Ω) A

HAk+4 (Ω)

.

A direct calculation shows that   ] (ψ)T(ψ) = −β(Ar + Bz ) ψ 2 T(ψ) + 2|ψ|2 T(ψ) [T, V − βψ 2 (Ar + Bz )T(ψ) + 2βψ T(ψ)(Ar + Bz )ψ  + 2β ψ T(ψ)(Ar + Bz )ψ + T(ψ)ψ(Ar + Bz )ψ  + |ψ|2 (Ar + Bz )T(ψ) .

Found Comput Math

Furthermore, by (4.4) and (4.11),      (Ar + Bz ) ψ 2 T(ψ)  k = ψ 2 T(ψ)H k+2 (Ω) ≤ c ψ 2 H (Ω) A

HAk+2 (Ω)

A

=

   T (ψ)H k+2 (Ω) A

c ψ 2 k+2 ψ H k+4 (Ω) . HA (Ω) A

By using (4.2)–(4.4) and (4.11), the same results can be derived for the estimates of ] (ψ)T(ψ). Therefore, we have the other terms in [T, V   [T, V ] (ψ)T(ψ) k ≤ c ψ 2 k+2 ψ H k+4 (Ω) . H (Ω) HA

A

(Ω)

A



A combination of the previous statements leads to the desired result.

Lemma 4.4 (Local errors) Let integer k ≥ 0. If the exact solution ψ(t) of (1.4) is in HA2 (Ω) ∩ HAk (Ω) for all 0 ≤ t ≤ τ , and ψ0 ∈ HAk+2 (Ω), then the local error of the method (3.3) is bounded by  1  ψ − ψ(τ ) k ≤ cτ 2 , (4.16) H (Ω) A

where c depends only on ψ0 H k+2 (Ω) and max0≤t≤τ ψ H j (Ω) with j = max(k, 2). A

A

If, in addition, ψ0 ∈ HAk+4 (Ω), then   1 ψ − ψ(τ )

HAk (Ω)

≤ cτ 3 ,

(4.17)

where c depends only on ψ0 H k+4 (Ω) and max0≤t≤τ ψ H j (Ω) . A

A

The proof of Lemma 4.4 is given in Appendix. The following lemma shall be used for the error analysis. Lemma 4.5 (Regularity of the numerical solution) If the exact solution of (1.4) ψ(t) ∈ HAk+2 (Ω) with integer k ≥ 2 and t ∈ [0, T ], then for small enough τ and any 1 ≤ n ≤ N0 = Tτ , we have max

 τ n    Φ ψ(j τ ) 

0≤j ≤N0 −n

HAk (Ω)

  ≤ T + max ψ(t)H k (Ω) . 0≤t≤T

A

Proof Let En = Fn =

max

0≤j ≤N0 −n

max

0≤j ≤N0 −n

 τ n    n−1     Φ ψ(j τ ) − Φ τ ψ (j + 1)τ H k (Ω) , A

 τ n    Φ ψ(j τ ) 

HAk (Ω)

,

n ≥ 1,

Then, by (4.16), we deduce that E 1 ≤ c0 τ 2 ,

n ≥ 1,

  F0 = max ψ(t)H k (Ω) . 0≤t≤T

A

Found Comput Math

where c0 depends only on max0≤t≤T ψ(t) H k+2 (Ω) . Moreover, A

F1 = ≤

max

 τ  Φ ψ(j τ ) 

max

 τ    Φ ψ(j τ ) − ψ (j + 1)τ  k + H (Ω)

0≤j ≤N0 −1

HAk (Ω)

0≤j ≤N0 −1

A

   ψ (j + 1)τ  k H (Ω)

max

0≤j ≤N0 −1

A

≤ E 1 + F0 . Next, due to Lemma 4.2, we have that for k ≥ 2,  2      2 2 E2 = max  Φ τ ψ(j τ ) − Φ τ ψ (j + 1)τ H k (Ω) ≤ E1 ecτ (F0 +F1 ) , 0≤j ≤N0 −2

A

and F2 =

 τ 2    Φ ψ(j τ ) 

max

0≤j ≤N0 −2

≤ E2 +

max

0≤j ≤N0 −2

HAk (Ω)

 τ   Φ ψ (j + 1)τ  k H (Ω) A

≤ E 2 + F1 . Finally by induction, we deduce that for k ≥ 2, ⎧ 2 +F 2 ) ⎨ En ≤ En−1 ecτ (Fn−1 n−2 , n ≥ 2, ⎩F ≤ E + F , n ≥ 1. n n n−1 Thereby, Fn ≤ Fn−1 + En−1 ecτ (Fn−1 +Fn−2 ) ≤ · · · 2

2

≤ Fn−1 + E1 ecτ (Fn−1 +2Fn−2 +···+2F1 +F0 ) 2

2

2

2

= Fn−1 + c0 τ 2 ecτ (Fn−1 +2Fn−2 +···+2F1 +F0 ) , 2

2

2

2

n ≥ 1.

Now let τ be small enough such that c0 τ e4cT F0 +4cT ≤ 1. 2

3

(4.18)

Then by induction, we derive that Fn ≤ Fn−1 + τ. Therefore, Fn ≤ F0 + nτ ≤ F0 + T ,

1 ≤ n ≤ N0 .

Thus we obtain the desired result. Remark 4.1 The condition (4.18) for τ is sufficient but not necessary.



Found Comput Math

We now present the main result for the Strang splitting scheme (3.3). Theorem 4.1 Suppose that integer k ≥ 0, τ is small enough, and the exact solution ψ(t) of (1.4) is in HA4 (Ω) ∩ HAk+2 (Ω) for all 0 ≤ t ≤ T . Then the numerical solution ψ n given by the splitting scheme (3.3) with step size τ > 0 has the following error bound:   n ψ − ψ(tn ) k ≤ cτ, tn = nτ ≤ T , H (Ω) A

where c depends only on T and max0≤t≤T ψ(t) H k+2 (Ω) . If, in addition, ψ(t) ∈ HAk+4 (Ω) for all 0 ≤ t ≤ T , then  n  ψ − ψ(tn ) k ≤ cτ 2 , H (Ω) A

A

tn = nτ ≤ T ,

where c depends only on T and max0≤t≤T ψ(t) H k+4 (Ω) . A

Proof According to Lemma 4.5, if ψ ∈ HA4 (Ω) ∩ HAk+2 (Ω) with k ≥ 0, then for m = max(k, 2),  n   max  Φ τ ψ(j τ ) H m (Ω) ≤ T + max ψ HAm (Ω) . 0≤j ≤N0 −n

0≤t≤T

A

Therefore, by Lemma 4.2, (4.16), and the previous inequality, and using the standard argument of Lady Windermere’s fan (cf. [12]), we obtain that n−1    τ n−j −1  τ   n   Φ ψ − ψ(nτ ) k Φ ψ(j τ ) ≤ H (Ω) A

j =0

  n−j −1   ψ (j + 1)τ H k (Ω) − Φτ A  τ    c1 T ≤ ne max Φ ψ(j τ ) − ψ (j + 1)τ H k (Ω) 0≤j ≤n−1

A

≤ c2 T ec1 T τ, where c1 depends only on T and max0≤t≤T ψ(t) HAm (Ω) , and c2 depends only on

max0≤t≤T ψ(t) H k+2 (Ω) . If, in addition, ψ(t) ∈ HAk+4 (Ω) for all 0 ≤ t ≤ T , then A we obtain from (4.17) that  n  ψ − ψ(nτ ) k ≤ c3 T ec1 T τ 2 , H (Ω) A

where c3 depends only on max0≤t≤T ψ(t) H k+4 (Ω) . This completes the proof. A

5 Error Analysis of the Fully Discrete Strang Splitting Laguerre–Hermite Collocation Scheme In this section, we present error bounds for the full-discretization scheme (3.9).



Found Comput Math

Lemma 5.1 We are given ψ ∈ XN . Then for ϕ ∈ HAs (Ω) with integer s ≥ 2,   τ Φ (ψ) − PN Φ τ (ϕ)2 2 N

L (Ω)

   ≤ exp cτ ψ H 2 (Ω) ϕ H 2 (Ω) exp cτ ϕ 2H 2 (Ω) ψ − PN ϕ 2L2 (Ω) A

+ cτ N 2

3 2 −s

A

A

ϕ 2H 2 (Ω) ϕ 4H s (Ω) A A

   × exp cτ ϕ 2H s (Ω) + cτ ψ H 2 (Ω) ϕ H 2 (Ω) exp cτ ϕ 2H 2 (Ω) . A

A

A

A

Proof We consider the errors in its linear and nonlinear parts, respectively. We first consider the errors of the linear parts: i∂t θ = (Ar + Bz )θ,

θ (0) = ψ,

(5.1)

i∂t η = (Ar + Bz )η,

η(0) = ϕ,

(5.2)

where θ and η correspond to the solutions of the linear parts of the full discretization and the semidiscretization in time, respectively. From (5.2) and (2.32) we further obtain     i(∂t PN η, φ)Ω = i(∂t η, φ)Ω = (Ar + Bz )η, φ Ω = (Ar + Bz )PN η, φ Ω , ∀φ ∈ XN . Therefore,     i ∂t (θ − PN η), φ Ω = (Ar + Bz )(θ − PN η), φ Ω ,

∀φ ∈ XN .

(5.3)

Taking φ = (Ar + Bz )k (θ − PN η) (∈ XN ) in the above, we obtain from its imaginary part that ∂t θ − PN η 2H k (Ω) = 0, A

whence   θ (·, t) − (PN η)(·, t) k = ψ − PN ϕ H k (Ω) , H (Ω) A

A

0 ≤ t ≤ τ.

(5.4)

Next let ψ˜ ∈ XN and consider the errors of the nonlinear parts:   ˜ i∂t θ = βIN |θ |2 θ , θ (0) = ψ,

(5.5)

i∂t η = β|η|2 η,

(5.6)

η(0) = ϕ. ˜

˜ j , zk )|, |η(r, z, t)| = From (5.5) and (5.6), one verifies readily that |θ (rj , zk , t)| = |ψ(r |ϕ(r, ˜ z)| and   i∂t (θ − PN η) = βIN |θ |2 + |PN η|2 (θ − PN η) + θ PN η(θ − PN η)

Found Comput Math

   + |PN η|2 PN η − βPN |η|2 η . Thus   i ∂t (θ − PN η), φ Ω    = β |θ |2 + |PN η|2 (θ − PN η), φ N   + β θ PN η(θ − PN η), φ N     + β |PN η|2 PN η, φ N − β |η|2 η, φ Ω ,

∀φ ∈ XN .

(5.7)

Next let M = [ N3 ]. To deal with the cubic term, we need to consider a special orthogonal projection, similar to (2.31), but with γ3r and γ3z in place of γr and γz . For M : L2 (Ω) → XM (γr /3, γz /3). More specifically, clarity, we denote it by P M u, φ)Ω = 0, (u − P

∀φ ∈ XM (γr /3, γz /3).

(5.8)

In particular, we have the following attractive property: M ψ| M η ∈ XN (γr , γz ) and ˜ 2P |P     M ψ| M ψ| M η, θ − PN η = |P M η, θ − PN η . ˜ 2P ˜ 2P |P N Ω Moreover, the same estimate in Theorem 2.1 holds for this variation. Hence, taking φ = θ − PN η (∈ XN ) in (5.7), we obtain from its imaginary part that     ∂t θ − PN η 2L2 (Ω) ≤ 2β  θ PN η(θ − PN η), θ − PN η N + |PN η|2 PN η, θ − PN η N    − |η|2 η, θ − PN η Ω  ≤ 2β

3 

|Gj |,

(5.9)

j =1

where   G1 = θ PN η(θ − PN η), θ − PN η N ,   M η|2 P M η, θ − PN η , G2 = |PN η|2 PN η − |P N   2 2  G3 = |PM η| PM η − |η| η, θ − PN η Ω . We derive from (5.9) that θ − PN η ∂t θ − PN η ≤ β

3 

|Gj |.

(5.10)

j =1

We now estimate |Gj |, j = 1, 2, 3. For any function u(r, z), (r, z) ∈ Ω, we denote  u(x, y, z) := u(r, z), (x, y, z) ∈ R3 . Then, by using a Sobolev inequality (cf., for instance, [7]) and the equivalence of norms (cf. Lemma 2.1), we obtain that for any

Found Comput Math

u ∈ HAm (Ω) with integer m > 32 , 3  3 1− 2m u L∞ (Ω) =  u L∞ (R3 ) ≤ c∂xm uL2m2 (R3 )  u L2 (R 3) 3

1−

3

3

1−

3

2m 2m 2m ≤ c  u H2mm (R3 )  u L2 (R 3 ) = c u H m (Ω) u L2 (Ω) ,

(5.11)

A

A

which implies in particular that u L∞ (Ω) ≤ C u H 2 (Ω) .

(5.12)

A

˜ j , zk )|, we have Due to |θ (rj , zk , t)| = |ψ(r ˜ L∞ (Ω) PN η L∞ (Ω) θ − PN η 2 2 |G1 | ≤ ψ L (Ω) ˜ 2 η 2 θ − PN η 2 2 . ≤ c ψ H (Ω) H (Ω) L (Ω) A

A

(5.13)

Moreover,  M η) + PN ηP M η(PN η − P M η) |G2 | =  |PN η|2 (PN η − P

  M η(PN η − P M η|PN η − P M η|2 , θ − PN η  M η) − P + PN ηP N     ≤ PN η L∞ (Ω) + PM η L∞ (Ω) PN η − PM η L∞ (Ω) PN η L2 (Ω) × θ − PN η L2 (Ω) M η 2 ∞ IN P M η L2 (Ω) θ − PN η L2 (Ω) . + PN η − P L (Ω)

By (5.11) and Theorem 2.1 we obtain that for s > 32 , M η L∞ (Ω) PN η − P 3

1−

3

M η 2s s PN η − P M η 2 2s ≤ c PN η − P H (Ω) L (Ω) A

3  M η H s (Ω) 2s ≤ c PN η − η HAs (Ω) + η − P A  3  M η L2 (Ω) 1− 2s × PN η − η L2 (Ω) + η − P 3

≤ cN 4 − 2 η HAs (Ω) . s

(5.14)

Therefore, by the above and (2.34), we can derive 3

|G2 | ≤ cN 4 − 2 η H 2 (Ω) η 2H s (Ω) θ − PN η L2 (Ω) . s

A

A

(5.15)

It is also verified readily that |G3 | ≤ cN − 2 η 2H 2 (Ω) η HAs (Ω) θ − PN η L2 (Ω) . s

A

(5.16)

Found Comput Math

Since the operator e−i 2 (Ar +Bz ) preserves the norms · HAs (Ω) , we have τ

˜ H s (Ω) = ψ H s (Ω) , ψ A A

ϕ ˜ HAs (Ω) = ϕ HAs (Ω) ,

˜ L2 (Ω) = ψ − PN ϕ L2 (Ω) , ψ˜ − PN ϕ

∀s ≥ 0.

Furthermore, by (4.10), we find   η(t)

HAs (Ω)

0 ≤ t ≤ τ,

≤e

ct ϕ 2H s (Ω)

ϕ HAs (Ω) ≤ e

A

cτ ϕ 2H s (Ω) A

ϕ HAs (Ω) ,

∀ integer s ≥ 2.

Hence, a combination of (5.10), (5.13)–(5.16), and (5.4) yields ∂t θ − PN η L2 (Ω) ≤ c ψ H 2 (Ω) ϕ H 2 (Ω) e A

3

cτ ϕ 2

2 (Ω) HA

A

+ cN 4 − 2 ϕ H 2 (Ω) ϕ 2H s (Ω) e s

A

θ − PN η L2 (Ω)

cτ ϕ 2H s (Ω) A

A

.

Finally, by using the Gronwall inequality and (5.4), we obtain the desired result.  Theorem 5.1 We assume ψ(t) ∈ HAs+2 (Ω) (0 ≤ t ≤ T ) with integer s ≥ 4. Then for N sufficiently large and τ sufficiently small (cf. (5.21) below), we have for integer s − 2 < k ≤ s,  n  3 k−s ψ − ψ(tn ) k ≤ cN 4 + 2 + cτ, N H (Ω) A

1≤n≤

T ; τ

(5.17)

1≤n≤

T . τ

(5.18)

and for integer 0 ≤ k ≤ s − 2,  n  3 k−s ψ − ψ(tn ) k ≤ cN 4 + 2 + cτ 2 , N H (Ω) A

Proof Due to ψ(t) ∈ HAs+2 (Ω) with integer s ≥ 4, we have from Lemma 4.5 that there exists a constant M0 > 0, such that for all 1 ≤ n ≤ Tτ ,  n ψ 

HAs (Ω)

≤

M0 , 2

(5.19)

provided that τ is small enough. Without loss of generality, we also assume that ψ0 H s+2 (Ω) ≤ M20 . Next, according to Theorem 2.2, we have A

1

5

ψ0 − IN ψ0 H 2 (Ω) ≤ c(ln N ) 2 N 6 − 2 ψ0 H s+2 . s

A

A

Hence, we obtain that for large N , IN ψ0 H 2 (Ω) ≤ ψ0 H 2 (Ω) + ψ0 − IN ψ0 H 2 (Ω) ≤ M0 . A

A

A

(5.20)

Found Comput Math

Now let N be large enough and τ be small enough such that 11

1

1 2 2 4 cτ M0

c2 ln NN − 6 e 2 cT M0 e

+

1 cτ M 2 1 1 2 2 1 0) 4 ≤ 1. c τ T M04 e 2 cT M0 ( 2 +e 16

(5.21)

We proceed by induction on n. (i) n = 1. By Lemma 5.1, (5.19), (5.20), and Theorems 2.1 and 2.2, we get that  1  ψ − PN ψ 1 2 2 N L (Ω) ≤e

cτ IN ψ0 H 2 (Ω) ψ0 H 2 (Ω) e A

cτ ψ0 2 2 HA (Ω)

A

IN ψ0 − PN ψ0 2L2 (Ω)

3

+ cτ 2 N 2 −s ψ0 2H 2 (Ω) ψ0 4H s (Ω) A

A

×e

cτ ψ0 2H s (Ω) +cτ IN ψ0 H 2 (Ω) ψ0 H 2 (Ω) e A

A

A

1 2 3 −s 6 1 cτ M 2 ( 1 +e 14 cτ M02 ) 0 2 cτ N 2 M0 e 2 64 1 1 cτ M 2 2 1 1 3 1 1 2 4 cτ M0 2 1 0) 4 ≤ c ln NN −s− 3 ψ0 2 s+2 e 2 cτ M0 e + cτ 2 N 2 −s M06 e 2 cτ M0 ( 2 +e HA (Ω) 64 1

1 2 2 4 cτ M0

cτ ψ0 2 2 HA (Ω)

≤ e 2 cτ M0 e

IN ψ0 − PN ψ0 2L2 (Ω) +

1 1 cτ M 2 2 1 1 3 1 1 1 2 4 cτ M0 2 1 0) 4 ≤ c ln N N −s− 3 M02 e 2 cτ M0 e + cτ 2 N 2 −s M06 e 2 cτ M0 ( 2 +e . 4 64

This with Lemma 2.4 and (5.21) gives  1  ψ − PN ψ 1 2 k N H (Ω) ≤ cN

A



k

ψN1

2 − PN ψ 1 L2 (Ω)

1 1 cτ M 2 2 1 1 3 1 1 1 2 4 cτ M0 2 1 0) 4 ≤ c2 ln N N k−s− 3 M02 e 2 cτ M0 e + c2 τ 2 N k+ 2 −s M06 e 2 cτ M0 ( 2 +e 4 64

≤

M02 k+ 3 −s N 2 . 4

(5.22)

In particular, by (5.19) and (5.22) we get that for integer s ≥ 4,  1   1   1 1 ψ  2   N HA (Ω) ≤ PN ψ HA2 (Ω) + ψN − PN ψ HA2 (Ω)     ≤ ψ 1 H 2 (Ω) + ψN1 − PN ψ 1 H 2 (Ω) ≤ M0 . A

A

(ii) Next assume that the results (5.22) and (5.23) with n = m hold, namely,  m  ψ − PN ψ m 2 N

HAk (Ω)

 2 ≤ cN k ψNm − PN ψ m L2 (Ω)

(5.23)

Found Comput Math 1 2 1 1 1 2 4 cτ M0 ≤ c2 ln N N k−s− 3 M02 e 2 mcτ M0 e 4 1 cτ M 2 3 1 1 2 1 0) 4 + mc2 τ 2 N k+ 2 −s M06 e 2 mcτ M0 ( 2 +e 64

≤

M02 k+ 3 −s N 2 , 4

(5.24)

and  m ψ 

N HA2 (Ω)

≤ M0 .

(5.25)

We now verify the results with n = m + 1. Clearly, by Lemma 5.1, (5.25), and (5.19), we derive that 2  m+1 ψ − PN ψ m+1 L2 (Ω) N ≤e

m cτ ψN ψ m H 2 (Ω) e H 2 (Ω) A

cτ ψ m 2 2 HA (Ω)

A

N

 2 3 + cτ 2 N 2 −s ψ m 

HA2 (Ω)

×e ≤e

 m  ψ − PN ψ m 2 2

L (Ω)

 m 4 ψ 

HAs (Ω)

m cτ ψ m 2H s (Ω) +cτ ψN ψ m H 2 (Ω) e H 2 (Ω) A

A

1 2 1 2 4 cτ M0 2 cτ M0 e

A

 m  ψ − PN ψ m 2 2 N

cτ ψ m 2 2 HA (Ω)

L

+ (Ω)

1 2 3 −s 6 1 cτ M 2 ( 1 +e 14 cτ M02 ) 0 2 cτ N 2 M0 e 2 . 64

The above with Lemma 2.4 and (5.24) lead to  m+1 2 ψ − PN ψ m+1 

HAk (Ω)

N

 2 ≤ cN k ψ m+1 − PN ψ m+1  2 N

1

L (Ω)

1 2 2 4 cτ M0

≤ cN k e 2 cτ M0 e +

 m  ψ − PN ψ m 2 2 N

L (Ω)

1 2 2 k+ 3 −s 6 1 cτ M 2 ( 1 +e 14 cτ M02 ) 0 2 c τ N 2 M0 e 2 64

1 2 1 1 1 2 4 cτ M0 ≤ c2 ln N N k−s− 3 M02 e 2 (m+1)cτ M0 e 4 1 cτ M 2 3 1 1 2 1 0) 4 + (m + 1)c2 τ 2 N k+ 2 −s M06 e 2 (m+1)cτ M0 ( 2 +e . 64

Hence by (5.21) and (5.26), for (m + 1)τ ≤ T ,  m+1 2 ψ − PN ψ m+1  N

HAk (Ω)

≤

M02 k+ 3 −s N 2 . 4

(5.26)

Found Comput Math

In particular, for integer s ≥ 4,  m+1      ψ  2 ≤ PN ψ m+1 H 2 (Ω) + ψNm+1 − PN ψ m+1 H 2 (Ω) N HA (Ω) A A  m+1   m+1  m+1  2  2 ≤ ψ + ψN − PN ψ ≤ M0 . (5.27) H (Ω) H (Ω) A

A

Therefore, we obtain that for integer s ≥ 4,  n  M02 k+ 3 −s ψ − PN ψ n 2 k N 2 , ≤ N HA (Ω) 4

n≤

T , τ

(5.28)

and  n ψ 

N HA2 (Ω)

≤ M0 ,

n≤

T . τ

(5.29)

Since  n      ψ − ψ(tn ) k ≤ ψNn − PN ψ n H k (Ω) + ψ n − PN ψ n H k (Ω) N HA (Ω) A A  n    + ψ − ψ(tn ) H k (Ω) , (5.30) A

we use (5.28), Theorem 2.1, (5.19), and Theorem 4.1 successively to derive the results (5.17) and (5.18).  The restriction on N in (5.21) can be removed if we set ψN0 (r, z) = PN ψ0 (r, z) instead of ψN0 (r, z) = IN ψ0 (r, z) in (3.9). More precisely, we have the following result. Corollary 5.1 Let ψN0 (r, z) = PN ψ0 (r, z) in (3.9) and ψ(t) ∈ HAs+2 (Ω) (0 ≤ t ≤ T ) with integer s ≥ 4. Then for τ sufficiently small such that 1 cτ M 2 1 1 2 2 1 0) 4 c τ T M04 e 2 cT M0 ( 2 +e ≤ 1, 16

(5.31)

we have for s − 2 < k ≤ s,  n  3 k−s ψ − ψ(tn ) k ≤ cN 4 + 2 + cτ, N H (Ω) A

1≤n≤

T ; τ

(5.32)

1≤n≤

T . τ

(5.33)

and for 0 ≤ k ≤ s − 2,  n  3 k−s ψ − ψ(tn ) k ≤ cN 4 + 2 + cτ 2 , N H (Ω) A

The proof of the above result is essentially the same as that of Theorem 5.1 without using (5.20) and with (5.31) instead of (5.21).

Found Comput Math

6 The d-Dimensional Gross–Pitaevskii Equation In this section, we shall present error bounds of the Strang splitting Hermite collocation method for the d-dimensional GPE:  i∂t ψ(x, t) = −ψ(x, t) + x2 ψ(x, t) + |ψ(x, t)|2 ψ(x, t), (6.1) ψ(x, 0) = ψ0 (x), lim|x|→∞ ψ(x, t) = 0, t ≥ 0, where x = (x1 , . . . , xd ) ∈ Rd and x2 = x12 + · · · + xd2 . 6.1 Notation and Some Basic Results For simplicity, we still denote by hl (z) the one-dimensional Hermite function as defined in (2.7) with γz ≡ 1. For l = (l1 , . . . , ld ) ∈ Nd , the d-dimensional Hermite function is defined by hl (x) = hl1 (x1 ) · · · hld (xd ), which satisfies   Lhl (x) := − + x2 hl (x) = μl hl (x),

μl = 2(l1 + · · · + ld ) + d.

We denote the spaces L2 (Rd ) and H s (Rd ), s > 0 with the inner products, seminorms, and norms as usual. For any u ∈ L2 (Rd ), we write u(x) =

∞ 

ul hl (x).

l=0

As before, the linear differential operator L is positive definite and self-adjoint. Thus for any u and v in the domain of L, (Lu, v)Rd = (u, Lv)Rd = ad (u, v), (Lu, u)Rd = ad (u, u) > 0,

if u = 0,

(6.2)

where the bilinear form   ad (u, v) = (∇u, ∇v)Rd + x2 u, v Rd . In particular, ad (hl , hl  ) = (Lhl , hl  )Rd = μl δll  .

(6.3)

The fractional power L1/2 is well defined, and the associated norms can be characterized by  1/2 2 L u d = ad (u, u), (6.4) R  m+1/2 2   L uRd = ad Lm u, Lm u ,

∀m ∈ N.

(6.5)

Found Comput Math

We also introduce the following three Sobolev spaces equipped with the norms:  u H s (Rd ) =

∞ 

A

1

 s  u H s (Rd ) = L 2 uL2 (Rd ) ,

2

μsl |ul |2

,

B

l=0

 u H s (Rd ) = C

s  l1 +···+ld =0

2  2  s−l1 −···−ld l 2  x + · · · + x2 + 1 ∂x11 · · · ∂xldd uL2 (Rd ) d 1

1 2

.

According to [25], we have the following. Lemma 6.1 The previous three norms are equivalent, i.e., u H s (Rd ) = u H s (Rd ) ∼ u H s (Rd ) . A

B

C

Remark 6.1 For d = 1, we may also refer to [24]. Hereafter, let d = 6κ + κ0 with integers κ ≥ 0 and 0 ≤ κ0 ≤ 5. For convenience, we assume that (i) if κ0 = 0, then sd = 5κ,

(ii) if κ0 = 1, then sd = 5κ + 1,

(iii) if κ0 = 2, 3, then sd = 5κ + 2, (iv) if κ0 = 4, 5, then sd = 5κ + 4.

(6.6)

Lemma 6.2 We have the following inequalities: uvw L2 (Rd ) ≤ c u

d

H 3 (Rd )

v

uvw L2 (Rd ) ≤ c u L2 (Rd ) v

d

H 3 (Rd )

H

w

d + 2 (Rd )

d

H 3 (Rd )

w

H

d + 2 (Rd )

uvw H k (Rd ) ≤ c u H k (Rd ) v H sd (Rd ) w H sd (Rd ) , ∀ integers 1 ≤ k < sd , uvw H k (Rd ) ≤ c u H k (Rd ) v H k (Rd ) w H k (Rd ) , ∀ integers k ≥ sd .

(6.7)

, ,

∀ > 0,

(6.8) (6.9)

(6.10)

Proof We proceed to treat different cases separately. d (i) The first bound follows from the Sobolev embedding H 3 (Rd ) ⊂ L6 (Rd ), and d the second bound comes from the Sobolev embedding H 2 + (Rd ) ⊂ L∞ (Rd ). (ii) We now deal with (6.9). For simplicity, we denote ∂xl u =

 l1 +···+ld =l

∂xl11 · · · ∂xldd u.

Found Comput Math

(1) 1 ≤ k < d3 . Since k is an integer number, we use (6.8) and (6.6) to deduce that for small enough  > 0, uvw H k (Rd ) ≤ u H k (Rd ) v

H

k+ d2 +

(Rd )

w

H

k+ d2 +

(Rd )

≤ u H k (Rd ) v H sd (Rd ) w H sd (Rd ) . (2)

≤ k < sd . We consider the term ∂xl u∂xm v∂xn w L2 (Rd ) . It is clear that l + m + n ≤ k. Hence (a) if max(m, n) < d3 , then by (6.8) and (6.6), for small enough  > 0,  l m n  ∂ u∂ v∂ w  2 d ≤ u l d v w n+ d + d H (R ) x x x m+ d + L (R ) d d 3

H

2

(R )

H

(R )

2

≤ u H k (Rd ) v H sd (Rd ) w H sd (Rd ) . (b) if d3 ≤ max(m, n) ≤ d2 , then l ≤ k − d3 . Since m, n are integer numbers, we use (6.7) and (6.6) to obtain that  l m n  ∂ u∂ v∂ w  2 d ≤ u d v m+ d d w n+ d d x x x l+ L (R ) d 3

H

(R )

3

H

(R )

3

H

(R )

≤ u H k (Rd ) v H sd (Rd ) w H sd (Rd ) . (c) if max(m, n) > d2 , for instance, m > d2 , then l, n < k − d2 . Therefore, by (6.8) we get that, for small enough  > 0,  l m n  ∂ u∂ v∂ w  2 d ≤ u d v H m (Rd ) w n+ d + d x x x l+ + L (R ) d H

2

(R )

H

2

(R )

≤ u H k (Rd ) v H sd (Rd ) w H sd (Rd ) . A combination of the previous statements leads to the desired result. (iii) It remains to consider (6.10). Clearly, the result (6.10) with k = sd can be proved in a similar fashion as before. Hence, by induction we obtain the result (6.10) with k > sd .  The previous results can be extended to the corresponding ones in HAk (Rd ), as stated below. Lemma 6.3 The following inequalities hold: uvw L2 (Rd ) ≤ c u

d

HA3 (Rd )

v

uvw L2 (Rd ) ≤ c u L2 (Rd ) v

d

HA3 (Rd ) d +

HA2

w

(Rd )

d

HA3 (Rd )

w

d +

HA2

, (Rd )

uvw H k (Rd ) ≤ c u H k (Rd ) v H sd (Rd ) w H sd (Rd ) , A

A

A

A

∀ integers 1 ≤ k < sd , uvw H k (Rd ) ≤ c u H k (Rd ) v H k (Rd ) w H k (Rd ) , A

∀ integers k ≥ sd .

A

A

A

(6.11)

, ∀ > 0,

(6.12)

(6.13)

(6.14)

Found Comput Math

Proof For simplicity, we only check the inequality (6.14). Clearly, uvw 2H k (Rd ) A

≤ c uvw 2H k (Rd ) C

k 

=c

2  2  k−l1 −···−ld l 2  x + · · · + x2 + 1 ∂x11 · · · ∂xldd (uvw)L2 (Rd ) d 1

l1 +···+ld =0

 2 = c∂ k (uvw) 2 x

+c

L (Rd )

k−1 

 2  k−l−m−n l m n 2 2  x + · · · + x2 + 1 ∂ u∂ v∂ w  2 1

x

d

x

x

L (Rd )

.

(6.15)

l+m+n=0 k−l−m−n

We next consider the term (x12 + · · · + xd2 + 1) 2 ∂xl u∂xm v∂xn w 2L2 (Rd ) . Without loss of generality, we assume that l = max(l, m, n). Since l + m + n ≤ k − 1, we have that m, n ≤ k−1 3 . Therefore, by (6.8) and the definition of the norm · HCk (Rd ) , a direct calculation shows that, for small enough  > 0,  2  k−l−m−n l m n  2  x + · · · + x2 + 1 ∂ u∂ v∂ w  1

x

d

x

x

L2 (Rd )

    k−l−m−n l  2 ≤  x12 + · · · + xd2 + 1 ∂x uL2 (Rd ) ∂xm v 

H

     k−l ≤  x12 + · · · + xd2 + 1 2 ∂xl uL2 (Rd ) ∂xm v 

H

d + 2 (Rd )

d + 2 (Rd )

 n  ∂ w  x

 n  ∂ w  x

H

H

d + 2 (Rd )

d + 2 (Rd )

≤ u H k (Rd ) v H k (Rd ) w H k (Rd ) . C

C

C

Further, by (6.15), (6.10), the above inequality and the equivalence of the norms, we derive that uvw H k (Rd ) ≤ c u H k (Rd ) v H k (Rd ) w H k (Rd ) . A

A

A

A

(6.16) 

Let us denote   XN = span hl1 (x1 ) · · · hld (xd ) : 0 ≤ l1 , . . . , ld ≤ N . For any u ∈ L2 (Rd ), the orthogonal projection operator PN : L2 (Rd ) → XN is defined by (u − PN u, φ)Rd = 0,

∀φ ∈ XN .

(6.17)

In particular, if u ∈ HBs (Rd ) with integer s ≥ 0, then we have    s s  L 2 (u − PN u), L 2 φ Rd = u − PN u, Ls φ Rd = 0,

∀φ ∈ XN ,

(6.18)

Found Comput Math

which means that the L2 (Rd )-orthogonal projection operator PN is also the HBs (Rd )orthogonal projection operator. We recall below an approximation result and an inverse inequality (cf. for instance [25]). Lemma 6.4 If u ∈ HAs (Rd ), then for any 0 ≤ μ ≤ s, u − PN u H μ (Rd ) ≤ cN

μ−s 2

u H s (Rd ) . A

A

Lemma 6.5 For any φ ∈ XN and s ≥ 0, s

φ H s (Rd ) ≤ cN 2 φ L2 (Rd ) . A

We now consider the interpolation operator. The Hermite–Gauss interpolant IN : C(Rd ) → XN is determined by IN u(zl1 , . . . , zld ) = u(zl1 , . . . , zld ),

0 ≤ l1 , . . . , ld ≤ N,

where {zj } are the Hermite–Gauss points. According to (2.11), we deduce readily that IN u L2 (Rd ) ≤ c

d 

N − 6 |u|H m (Rd ) , m

(6.19)

m=0

where | · |H m (Rd ) denotes the seminorm of H m (Rd ). Hence, by a similar argument as in the proof of Theorem 2.2, we can prove the following result: Theorem 6.1 If u ∈ HAs (Rd ), then for any 0 ≤ μ ≤ s and s ≥ d, u − IN u H μ (Rd ) ≤ cN 3 + d

μ−s 2

u H s (Rd ) . A

A

Proof From Lemmas 6.4 and 6.5, for any 0 ≤ μ ≤ s,   u − IN u H μ (Rd ) ≤ u − PN u H μ (Rd ) + IN (u − PN u)H μ (Rd ) A

A

≤ cN

A

  u H s (Rd ) + cN IN (u − PN u)L2 (Rd ) .

μ−s 2

μ 2

A

Moreover, by (6.19) and Lemma 6.4, for any s ≥ d,   IN (u − PN u)

L2 (Rd )

≤c

d 

N − 6 u − PN u H m (Rd ) ≤ cN 3 − 2 u H s (Rd ) . m

d

A

s

A

m=0

Therefore, u − IN u H μ (Rd ) ≤ cN 3 + d

A

μ−s 2

u H s (Rd ) . A



Found Comput Math

6.2 Error Bounds of the Semidiscretization in Time In this section, we present error bounds for the semidiscretization in time t under the HAk (Rd )-norms. The second-order Strang splitting in time for (6.1) is as follows:   −i τ L n 2 τ τ ψ n+1 = Φ τ ψ n := e−i 2 L e−iτ |e 2 ψ | e−i 2 L ψ n ,

ψ 0 = ψ0 .

(6.20)

Since the error analysis for the above scheme is essentially the same as for the scheme (3.3), we shall present the results below without proof. Lemma 6.6 (Stability) If ψ, ϕ ∈ HAsd (Rd ) ∩ HAk (Rd ) with integer k ≥ 0, then we have  τ  Φ (ψ) − Φ τ (ϕ)

≤e H k (Rd )

+ ϕ 2 j ) j HA (Rd ) HA (Rd )

cτ ( ψ 2

ψ − ϕ H k (Rd ) , A

A

(6.21)

where j = max(k, sd ). Next, we denote T(ψ) = −iLψ,

(ψ) = −i|ψ|2 ψ. V

Their Lie commutator (cf. [12, 15]) is as follows:   ](ψ) = T (ψ)V (ψ) − V  (ψ)T(ψ) = −L |ψ|2 ψ − ψ 2 Lψ + 2|ψ|2 Lψ. [T, V Lemma 6.7 If ψ ∈ HAk+2 (Rd ) ∩ HAsd (Rd ) with integer k ≥ 0, then the commutator is bounded by   [T, V ](ψ) k d ≤ c ψ 3 j , j = max(k + 2, sd ). HA (R )

HA (Rd )

If, in addition, ψ ∈ HAk+4 (Rd ), then     T, [T, V ] (ψ) k d ≤ c ψ 3 H (R )

j

HA (Rd )

A

,

j = max(k + 4, sd ).

Lemma 6.8 (Local errors) Let integer k ≥ 0. If the exact solution ψ(t) ∈ HAk (Rd ) ∩ HAsd (Rd ) for all 0 ≤ t ≤ τ , and ψ0 ∈ HAk+2 (Rd ) ∩ HAsd (Rd ), then the local errors of the method (6.20) are bounded by   1 ψ − ψ(τ ) k d ≤ cτ 2 , H (R ) A

where c depends only on ψ0 H m (Rd ) and max0≤t≤τ ψ H j (Rd ) with m = max(k + A

A

2, sd ) and j = max(k, sd ). If, in addition, ψ0 ∈ HAk+4 (Rd ), then   1 ψ − ψ(τ ) k d ≤ cτ 3 , H (R ) A

where c depends only on ψ0 H n (Rd ) and max0≤t≤τ ψ H j (Rd ) with n = max(k + A A 4, sd ).

Found Comput Math

Lemma 6.9 (Regularity of the numerical solution) If the exact solution ψ(t) ∈ HAk+2 (Rd ) with integer k ≥ sd and t ∈ [0, T ], then for small enough τ and any 1 ≤ n ≤ N0 = Tτ , we have  n   max  Φ τ ψ(j τ ) H k (Rd ) ≤ T + max ψ H k (Rd ) . 0≤j ≤N0 −n

0≤t≤T

A

A

By using the above results and the same procedure as in the proof of Theorem 4.1, we can prove the following. Theorem 6.2 Suppose that integer k ≥ 0, τ is small enough, and the exact solution s +2 ψ ∈ HAk+2 (Rd ) ∩ HAd (Rd ) for all 0 ≤ t ≤ T . Then, the solution of the scheme (6.20) satisfies  n  ψ − ψ(tn ) k d ≤ cτ, tn = nτ ≤ T , H (R ) A

where c depends only on T and max0≤t≤T ψ(t) H j (Rd ) with j = max(k + 2, sd ). If, A

in addition, ψ(t) ∈ HAk+4 (Rd ) for all 0 ≤ t ≤ T , then   n ψ − ψ(tn ) k d ≤ cτ 2 , tn = nτ ≤ T , H (R ) A

where c depends only on T and max0≤t≤T ψ(t) H m (Rd ) with m = max(k + 4, sd ). A

6.3 Error Bounds of the Full Discretization The fully discrete Strang splitting Hermite collocation scheme for (6.1) is as follows: ψN0 = IN ψ0 ;  n   −i τ2 L n 2 τ τ ψN | −i τ2 L n ψNn+1 = ΦN ψN := e−i 2 L IN e−iτ |e e ψN ,

(6.22) n ≥ 0.

Then, by using a similar procedure as in the proofs of Lemma 5.1 and Theorem 5.1, we can prove the following. Lemma 6.10 If ψ ∈ XN and ϕ ∈ HAs (Rd ) ∩ HAd (Rd ) with integer s ≥ sd , then  τ  Φ (ψ) − PN Φ τ (ϕ)2 2

L (Rd )

N

  ≤ exp cτ ψ H j (Rd ) ϕ H sd (Rd ) exp cτ ϕ 2 A

+ cτ N 2

d 2 −s

s

HAd (Rd )

A

 ψ − PN ϕ 2L2 (Rd )

ϕ 2H d (Rd ) ϕ 4H s (Rd ) A

A

  × exp cτ ϕ 2H m (Rd ) + cτ ψ H j (Rd ) ϕ H sd (Rd ) exp cτ ϕ 2 A

where j =

d 2

A

A

s

HAd (Rd )

 ,

+  and m = max(s, d).

Finally, with the above preparations and by using a similar argument as in the proof of Theorem 5.1, we can prove the following.

Found Comput Math T Theorem 6.3 Let integer s > 7d 6 − 2 and 1 ≤ n ≤ τ . If the exact solution ψ(t) ∈ HAs+2 (Rd ) for all 0 ≤ t ≤ T , then for N sufficiently large and τ sufficiently small, we have for s − 2 < k ≤ s, ⎧ d k−s ⎨ cN 4 + 2 + cτ, d < 12, s > d,   n ψ − ψ(tn ) k d ≤ (6.23) N HA (R ) d k−s ⎩ cN 3 −1++ 2 + cτ, d ≥ 12, ∀ > 0,

and for 0 ≤ k ≤ s − 2, ⎧ d k−s ⎨ cN 4 + 2 + cτ 2 ,  n  ψ − ψ(tn ) k d ≤ N HA (R ) d k−s ⎩ cN 3 −1++ 2 + cτ 2 ,

d < 12, s > d, (6.24) d ≥ 12, ∀ > 0.

As in the last section, the restriction on N in the above result can be removed if we replace ψN0 = IN ψ0 in (6.22) by ψN0 = PN ψ0 . More precisely, we can prove the following. Corollary 6.1 Let ψN0 = PN ψ0 in (6.22). Then for integer s > d, ψ(t) ∈ HAs+2 (Rd ), 0 ≤ t ≤ T , and τ sufficiently small, we have for s − 2 < k ≤ s,   n ψ − ψ(tn ) k d ≤ cN d4 + k−s 2 + cτ, N H (R ) A

1≤n≤

T , τ

(6.25)

1≤n≤

T . τ

(6.26)

and for 0 ≤ k ≤ s − 2,   n ψ − ψ(tn ) k d ≤ cN d4 + k−s 2 + cτ 2 , N H (R ) A

Remark 6.2 In Theorem 3.4 of [8], the author also presents a result under the L2 (Rd )2d norm, stated as follows: Let s > [ d+1 2 ] + 2 + 3 be an even integer. If the exact solution ψ(t) ∈ HAs+2 (Rd ) for 0 ≤ t ≤ T , then for N sufficiently large and τ sufficiently small,   n ψ − ψ(tn ) 2 d ≤ cN 1+ d3 − 2s + cτ 2 . N L (R ) It is clear that the estimates in (6.24) are better both with respect to the order of convergence in space as well as the required regularity. In addition, the order of convergence and the regularity requirements of our results are further relaxed significantly if we replace ψN0 = IN ψ0 in (6.22) by ψN0 = PN ψ0 , as indicated in Corollary 6.1. Acknowledgements The work of J. Shen was partially supported by NSF grant DMS-0915066 and AFOSR grant FA9550-11-1-0328. The work of Z.-Q. Wang was partially supported by the NSF of China, No. 11171225, the Shuguang Project of the Shanghai Education Commission, No. 08SG45, the Innovation Program of the Shanghai Municipal Education Commission, No. 12ZZ131, and the Fund for E-institute of Shanghai Universities, No. E03004.

Found Comput Math

Appendix: The Proof of Lemma 4.4

Proof The proof of Lemma 4.4 is analogous to the corresponding results for the Schrödinger–Poisson equation in [18]. For simplicity, we only verify (4.17). To this  = T + V , and denote by DH , DT , and DV the corresponding Lie derivaend, let H  , respectively. According to Sect. 4.4 of [18], tives (cf. [18]) of H , T, and V    τ 1 ψ 1 − ψ(τ ) = τf τ − f (s) ds + r2 − r1 , (A.1) 2 0 where f (s) = exp((τ − s)DT )DV exp(sDT )Id(ψ0 ), Id is the identity operator, and  τ  τ −s   exp (τ − s − σ )DH DV exp(σ DT )DV exp(sDT )Id(ψ0 ) dσ ds, r1 = 0

r2 = τ



2 0

0 1



1 (1 − θ ) exp τ DT 2

 exp(θ τ DV )DV2

  1 exp τ DT Id(ψ0 ) dθ. 2

Next we write the principal error term in the second-order Peano form:    τ  1 1 τ − τf f (s) ds = τ 3 ν(θ )f  (θ τ ) dθ 2 0 0 with the Peano kernel ν of the midpoint rule. We have (cf. Sect. 5.2 of [18])    f  (s) = exp (τ − s)DT DT , [DT , DV ] exp(sDT )Id(ψ0 )   = exp (τ − s)DT D[T,[T,V]] exp(sDT )Id(ψ0 )   ] (e−i(τ −s)(Ar +Bz ) ψ0 ), = e−is(Ar +Bz ) T, [T, V where [DT , DV ] = DT DV − DV DT . Hence by (4.14), the quadrature error is bounded in HAk (Ω) by cτ 3 ψ0 3 k+4 . Let us denote HA

(Ω)

  g(s, σ ) = exp (τ − s − σ )DT DV exp(σ DT )DV exp(sDT )Id(ψ0 ). Then the remainder term can be expressed as   τ  τ −s  1 2 1 r2 − r1 = τ g τ, 0 − g(s, σ ) dσ ds + r2 − r1 , 2 2 0 0 where 

τ

 r1 = r1 −



τ −s

g(s, σ ) dσ ds, 0

0

  1 2 1  r2 = r2 − τ g τ, 0 . 2 2

It is clear that (cf. [18])     τ  τ −s  1 2 1   τ g τ, 0 − g(s, σ ) dσ ds  2  k 2 0 0 H (Ω) A

Found Comput Math

 ≤ cτ 3 max ∂s g H k (Ω) + max ∂σ g H k (Ω) . 0≤s≤τ

0≤σ ≤τ

A

A

By using a similar argument as in Sect. 5.2 of [18], we can obtain     τ  τ −s  1 2 1   τ g τ, 0 −  g(s, σ ) dσ ds ≤ c1 τ 3 , 2  k 2 0 0 H (Ω) A

where c1 depends only on ψ0 H k+2 (Ω) . A Next, we estimate the term  r2 H k (Ω) . By the Taylor expansion, A

 exp(θ τ DV ) = I +

θτ



1

exp(ξ DV )DV dξ = I + θ τ

0

exp(θ τ ςDV )DV dς, 0

whence   r2 = τ 3 0

1 1 0

    1 1 θ (1 − θ ) exp τ DT exp(θ τ ςDV )DV3 exp τ DT Id(ψ0 ) dς dθ. 2 2

Setting φ = e−i 2 (Ar +Bz ) ψ0 and η = e−iθτ ςβ|φ| φ, a direct calculation shows that       τ 1 1 3 exp τ DT exp(θ τ ςDV )DV exp τ DT Id(ψ0 ) = iβ 3 e−i 2 (Ar +Bz ) |η|6 η . 2 2 τ

2

Hence, by (4.2)–(4.4) and (4.10), we find that for j = max(k, 2),         exp 1 τ DT exp(θ τ ςDV )D 3 exp 1 τ DT Id(ψ0 ) V   k 2 2 H (Ω) A

≤ c η 7 j HA (Ω)

≤ c2 ψ0 7 j , HA (Ω)

where c2 depends only on ψ0 H j (Ω) . Therefore,  r2 H k (Ω) ≤ c2 τ 3 . A A It remains to estimate the term  r1 H k (Ω) . By using the nonlinear variation-ofA constants formula (cf. [18]), we obtain that  τ  τ −s  τ −s−σ   exp (τ − s − σ − ξ )DH DV exp(ξ DT )DV  r1 = 0

0

0

× exp(σ DT )DV exp(sDT )Id(ψ0 ) dξ dσ ds. By (4.2)–(4.4), a direct calculation gives   exp((τ − s − σ − ξ )DH )DV exp(ξ DT )DV exp(σ DT )DV exp(sDT )Id(ψ0 ) k H (Ω) 7  ≤ cψ(τ − s − σ − ξ )H j (Ω) , A

A

j = max(k, 2).

Therefore,  r1 H k (Ω) ≤ c3 τ 3 , where c3 depends only on max0≤t≤τ ψ H j (Ω) . A A A combination of the previous statements leads to (4.17). 

Found Comput Math

References 1. W. Bao, J. Shen, A fourth-order time-splitting Laguerre–Hermite pseudo-spectral method for Bose– Einstein condensates, SIAM J. Sci. Comput. 26, 2010–2028 (2005). 2. W. Bao, D. Jaksch, P.A. Markowich, Numerical solution of the Gross–Pitaevskii equation for Bose– Einstein condensation, J. Comput. Phys. 187, 318–342 (2003). 3. N. Ben Abdallah, F. Castella, F. Méhats, Time averaging for the strongly confined nonlinear Schrödinger equation, using almost-periodicity, J. Differ. Equ. 245, 154–200 (2008). 4. J. Bergh, J. Löfström, Interpolation Spaces, An Introduction (Springer, Berlin, 1976). 5. M.M. Cerimele, M.L. Chiofalo, F. Pistella, S. Succi, M.P. Tosi, Numerical solution of the Gross– Pitaevskii equation using an explicit finite-difference scheme: an application to trapped Bose–Einstein condensates, Phys. Rev. E 62, 1382–1389 (2000). 6. E. Faou, Geometric Numerical Integration and Schrödinger Equations (European Math. Soc., Zürich, 2012). 7. A. Friedman, Partial Differential Equations (Holt, Rinehart and Winston, Inc., New York, 1969). 8. L. Gauckler, Convergence of a split-step Hermite method for the Gross–Pitaevskii equation, IMA J. Numer. Anal. 31, 396–415 (2011). 9. E. Gross, Structure of a quantized vortex in boson systems, Nuovo Cimento 20, 454–477 (1961). 10. B. Guo, J. Shen, C. Xu, Spectral and pseudospectral approximations using Hermite functions: application to the Dirac equation, Adv. Comput. Math. 19, 35–55 (2003). 11. B. Guo, L. Wang, Z. Wang, Generalized Laguerre interpolation and pseudospectral method for unbounded domains, SIAM J. Numer. Anal. 43, 2567–2589 (2006). 12. E. Hairer, S.P. Nörsett, G. Wanner, Solving Ordinary Differential Equations. I. Nonstiff Problems, 2nd edn., Springer Series in Computational Mathematics, vol. 8 (Springer, Berlin, 1993). 13. E. Hansen, A. Ostermann, Exponential splitting for unbounded operators, Math. Comput. 78, 1485– 1496 (2009). 14. B. Helffer, Théorie spectrale pour des opérateurs globalement elliptiques, Astérisque, vol. 112 (SMF, Paris, 1984). 15. W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Series in Computational Mathematics, vol. 33 (Springer, Berlin, 2003). 16. T. Jahnke, C. Lubich, Error bounds for exponential operator splittings, BIT Numer. Math. 40, 735–744 (2000). 17. O. Koch, C. Lubich, Variational-splitting time integration of the multi-configuration time-dependent Hartree–Fock equations in electron dynamics, IMA J. Numer. Anal. 31, 379–395 (2011). 18. C. Lubich, On splitting methods for Schrödinger–Poisson and cubic nonlinear Schrödinger equations, Math. Comput. 77, 2141–2153 (2008). 19. L.P. Pitaevskii, Vortex lines in an imperfect Bose gas, Sov. Phys. JETP 13, 451–454 (1961). 20. P.A. Ruprecht, M.J. Holland, K. Burrett, M. Edwards, Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms, Phys. Rev. A 51, 4704–4711 (1995). 21. R. Temam, Infinite Dimensional Dynamical System in Mechanics and Physics, Applied Mathematical Sciences, vol. 68 (Springer, New York, 1988). 22. M. Thalhammer, High-order exponential operator splitting methods for time-dependent Schrödinger equations, SIAM J. Numer. Anal. 46, 2022–2038 (2008). 23. L. Wang, Analysis of spectral approximations using prolate spheroidal wave functions, Math. Comput. 79, 807–827 (2010). 24. X. Xiang, Z. Wang, Generalized Hermite spectral method and its applications to problems in unbounded domains, SIAM J. Numer. Anal. 48, 1231–1253 (2010). 25. X. Xiang, Z. Wang, Generalized Hermite approximations and spectral method for partial differential equations in multiple dimensions, submitted.