Crank–Nicolson difference scheme for the coupled nonlinear Schrödinger equations with the Riesz space fractional derivative

Accepted Manuscript Crank-Nicolson difference scheme for the coupled nonlinear Schrödinger equa‐ tions with the Riesz space fractional derivative Dongling Wang, Aiguo Xiao, Wei Yang PII: DOI: Reference:

S0021-9991(13)00156-3 http://dx.doi.org/10.1016/j.jcp.2013.02.037 YJCPH 4491

To appear in:

Journal of Computational Physics

Received Date: Revised Date: Accepted Date:

6 May 2012 16 February 2013 21 February 2013

Please cite this article as: D. Wang, A. Xiao, W. Yang, Crank-Nicolson difference scheme for the coupled nonlinear Schrödinger equations with the Riesz space fractional derivative, Journal of Computational Physics (2013), doi: http://dx.doi.org/10.1016/j.jcp.2013.02.037

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Crank-Nicolson difference scheme for the coupled nonlinear Schr¨ odinger equations with the Riesz space fractional derivative Dongling Wang

Aiguo Xiao∗ Wei Yang

Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Key Laboratory of Intelligent Computing & Information Processing of Ministry of Education, Xiangtan University, Xiangtan, Hunan 411105, P.R.China

Abstract In this paper, the Crank-Nicolson (CN) difference scheme for the coupled nonlinear Schr¨odinger equations with the Riesz space fractional derivative is studied. The existence of this difference solution is proved by the Brouwer fixed point theorem. The stability and convergence of the CN scheme are discussed in the L2 norm. When the fractional order is two, all those results are in accord with the difference scheme developed for the classical non-fractional coupled nonlinear Schr¨ odinger equations. Some numerical examples are also presented. Keywords: Fractional Schr¨odinger equation; Crank-Nicolson scheme; fractional centered difference

1

Introduction

The Schr¨odinger equation is one of the most important equations in quantum mechanics. Feynman and Hibbs derived the standard Schr¨ odinger equation from the path integrals over Brownian paths. A natural generalization of the Brownian motion is L´evy stochastic process. Recently, Laskin [1, 2] derived the space fractional Schr¨ odinger equation by extending the Feynman path integral to L´evy one. The fractional Schr¨ odinger equation contains the standard Schr¨odinger equation as its special case. Some authors [3, 4] discussed the physical applications of fractional Schr¨ odinger equation and obtained the exact solutions with several kinds of potentials. However, the exact solutions of fractional differential equations often contains some special functions, such as Fox functions, which are not studied very well. The numerical methods for fractional differential equations become important tools to understand the behaviors of the equations. Amore et.al. [5] developed the collocation method for fractional quantum mechanics. ∗

Corresponding author. E-mail addresses: [email protected] (Dongling Wang), [email protected] (Aiguo Xiao), [email protected] (Wei Yang)

1

The coupled nonlinear Schr¨odinger equations arise in many physical fields, especially in optics and hydrodynamics [6, 7]. In recent years, more and more researchers [8, 9, 10] found that fractional partial differential equations can describe some real physical phenomena better than integral partial differential equations. Guo, Han and Xin [11] proved the existence of the global smooth solution of fractional nonlinear Schr¨ odinger equation with periodic boundary conditions. Hu, Xin and Hu [12] further studied the existence and uniqueness of the global solution of fractional coupled nonlinear Schr¨ odinger equations. For solving the coupled nonlinear Schr¨ odinger equations, some authors developed the finite difference methods [13-17], pseudo-spectral methods [18], and multi-symplectic methods [19, 20]. The authors [13, 14, 17] discussed the stability and convergence of some difference schemes for coupled nonlinear Schr¨ odinger equations. Sun and Zhao [15] studied a conservative difference scheme and obtained the second order convergence result in L∞ norm. As far as we know, there exist few studies on numerical methods for fractional coupled nonlinear Schr¨odinger equations. In this paper, we propose the Crank-Nicolson (CN) difference scheme for fractional coupled nonlinear Schr¨ odinger equations and discuss its properties. This paper is organized as follows. In Section 2, the fractional coupled Schr¨ odinger equations are introduced. In Section 3, the fractional centered difference approximation to the Riesz fractional derivative is introduced and the CN difference scheme is proposed. The stability and convergence are studied. Some examples are analyzed in Section 4.

2

Space fractional Schr¨ odinger equations

We briefly present some definitions and properties of fractional derivatives which are used later. The left and right Riemann-Liouville fractional derivative of order α (n−1 < α < n) on infinite domain are defined as Z x 1 ∂n α (x − τ )n−1−α u(τ, t)dτ, −∞ Dx u(x, t) = Γ(n − α) ∂xn −∞ (1) Z +∞ 1 ∂n α n−1−α D u(x, t) = (τ − x) u(τ, t)dτ, x +∞ Γ(n − α) ∂xn x where Γ(x) is the gamma function. The Riesz fractional derivative is defined as £ α ∂α 1 α u(x, t) = −(−∆) 2 u(x, t) = − πα −∞ Dx u(x, t) + α ∂|x| 2cos 2

α x D+∞ u(x, t)

¤

.

(2)

The Riesz fractional derivative can also be characterized as α

−(−∆) 2 u(x, t) = −F −1 (|ξ|α u b(ξ, t)) ,

(3)

where F is the Fourier transform. If u(x, t) is defined on the finite interval [a, b] and satisfies the boundary conditions u(a, t) = u(b, t) = 0, we can extend the function by taking u(x, t) ≡ 0 for x ≤ a and x ≥ b. 2

Furthermore, if u(x, t) satisfies that u0 (a, t) = u0 (b, t) = 0, by the Fourier transform (3), it is shown in [21] that the Riesz fractional derivative on the finite interval [a, b] can be defined as α ∂α 1 u(x, t) = −(−∆) 2 u(x, t) = − [a Dxα u(x, t) + x Dbα u(x, t)] , (4) α ∂|x| 2cos πα 2 where (1 < α ≤ 2). Some other definitions, for example, Gr¨ unwald-Letnikov (G-L), Weyl and Jumaire fractional derivatives, can be found in [8, 9, 10]. Our purpose in this paper is to consider the numerical solution of the coupled nonlinear Schr¨odinger equations with the Riesz space fractional derivative (1 < α ≤ 2) α

iut − γ(−∆) 2 u + ρ(|u|2 + β|v|2 )u = 0, α 2

2

2

ivt − γ(−∆) v + ρ(|v| + β|u| )v = 0,

a < x < b, 0 < t ≤ T, a < x < b, 0 < t ≤ T

(5)

with the initial conditions u(x, 0) = u0 (x), v(x, 0) = v0 (x),

(6)

and the Dirichlet boundary conditions u(a, t) = u(b, t) = 0, v(a, t) = v(b, t) = 0,

(7)

where i2 = −1, the parameters γ, ρ and β are some real constants. The initial conditions u0 (x) and v0 (x) are some given smooth functions vanishing at the end points x = a and x = b. When α = 2 this system becomes the classical coupled nonlinear Schr¨ odinger equations. When ρ = 0, this system is decoupled and becomes the fractional Schr¨ odinger equation of free particles solved in [10]. When β = 0, it reduces to the single fractional nonlinear Schr¨odinger equation considered in [11].

3 3.1

Crank-Nicolson difference scheme Numerical approximation

There are several numerical methods to approximate the Riesz fractional derivative. Yang, Liu and Turner [22] proposed the standard and shifted Gr¨ unwald formula and the L1 and L2 approximations. For homogeneous boundary conditions, Ili´c, Liu, Turner and Anh [23] proposed the matrix transform method (MTM) and this method is further studied in [24, 25]. A finite element method is presented in [26]. Some other efficient numerical methods about the Riesz fractional differential equations and some physical applications can be founded in [29-32]. We now introduce the fractional centered difference proposed by Ortigueira [27]. Suppressing the time variable and setting u(x) = u(x, t), the fractional centered difference for α > −1 is defined as ∆αh u(x)

=

+∞ X k=−∞

(−1)k Γ(α + 1) u(x − kh). Γ(α/2 − k + 1)Γ(α/2 + k + 1) 3

(8)

It is shown that for 1 < α ≤ 2, one have +∞ X ∆αh u(x) ∂α ck u(x) = − lim = − lim u(x − kh), α α h→0 h→0 ∂|x| h hα

(9)

k=−∞

k

(−1) Γ(α+1) where the coefficients ck = Γ(α/2−k+1)Γ(α/2+k+1) . C ¸ elik and Duman [28] discussed the errors of this approximation by the Fourier transformation method and applied it to fractional diffusion equations. Lemma 3.1[28] Let u ∈ C 5 (R), all derivatives up to order five belong to L1 (R) and the fractional centered difference ∆αh u(x) be given in (8). Then for 1 < α ≤ 2, we have c0 ≥ 0, ck = c−k ≤ 0 for all k ≥ 1, and that

−h−α ∆αh u(x) =

∂α u(x) + O(h2 ). ∂|x|α

(10)

We now discretize the equations (5) on [a, b] × [0, T ]. Let τ = T /N and h = (b − a)/M , and define Ωh = {xj = jh, j = 0, 1, · · · , M }, Ωτ = {tn = nτ, n = 0, 1, · · · , N }, Ωhτ = Ωh × Ωτ . Suppose that Ω = {wjn ; j = 0, 1, ..., M, n = 0, 1, ..., N } is a discrete function on Ωhτ and Ωh,0 = {w|w = (w0 , w1 , ..., wM ), w0 = wM = 0} be the space complex grid n ). Introduce the following notations: functions on Ωh and ω n = (ω0n , ω1n , ..., ωM n+ 12

wj

=

wjn+1 + wjn 2

(wjn )t =

,

||wn ||2 = hwn , wn i,

wjn+1 − wjn

hun , wn i = h

τ M −1 X

,

unj wjn .

j=1

Let Ujk and Vjk denote the numerical approximations to u(xj , tk ) and v(xj , tk ). The CN difference scheme for the fractional Schr¨ odinger equations (5) is ï ¯ ¯ ¯ ! M −1 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 ¡ n¢ γ X n+ 12 n+ 1 i Uj t − α cj−l Ul + ρ ¯¯Uj 2 ¯¯ + β ¯¯Vj 2 ¯¯ Uj 2 = 0, h l=1

M −1 ¡ ¢ γ X n+ 1 i Vjn t − α cj−l Vl 2 h

1 ≤ j ≤ M − 1, 1 ≤ n ≤ N − 1, ï ¯ ¯ ¯ ! ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 n+ 1 + ρ ¯¯Vj 2 ¯¯ + β ¯¯Uj 2 ¯¯ Vj 2 = 0,

(11)

l=1

1 ≤ j ≤ M − 1, Uj0 = u0 (xj ), U0n

=

n UM

= 0,

Vj0 = v0 (xj ), V0n

=

n VM

4

= 0,

1 ≤ n ≤ N − 1,

1 ≤ j ≤ M − 1, 1 ≤ n ≤ N − 1.

(12)

3.2

Numerical analysis of the Crank-Nicolson difference scheme

In this section, we analyze the solvability, stability and convergence of the CN difference scheme proposed in (11)-(12). Lemma 3.2 The formula   M −1 M −1 X X 1 1 n+ n+ Im  cj−l Ul 2 Uj 2  = 0 (13) j=1 l=1

holds, where “Im” means taking the imaginary part. Proof: Calculating directly, we have M −1 M −1 X X

n+ 1 cj−l Ul 2

n+ 1 Uj 2

= c0

M −1 X

j=1 l=1

n+ 12

Uj

n+ 12

Uj

+

l=1

M −2 M −1 X X

µ ¶ n+ 12 n+ 12 n+ 12 n+ 12 cj−l Ul Uj + cl−j Uj Ul .

(14)

l=1 j=l+1 n+ 12

Since cl−j = cj−l and Ul

n+ 12

Uj

n+ 12

= Uj c0

M −1 X

n+ 12

Ul

, we get that

n+ 1 n+ 1 Uj 2 Uj 2

j=1 M −2 M −1 X X

° c0 ° ° n+ 12 °2 = °U ° , h

µ ¶ n+ 12 n+ 12 n+ 12 n+ 12 cj−l Ul Uj + cl−j Uj Ul

(15)

l=1 j=l+1

= 2Re

M −2 M −1 X X

µ ¶ n+ 12 n+ 12 cj−l Ul Uj ,

l=1 j=l+1

which imply the result.

¤

Theorem 3.1 The scheme (11)-(12) is conservative in the sense Qn1 = Qn−1 = ... = Q01 , 1 Qn2 = Qn−1 = ... = Q02 , 2

(16)

° °2 ° °2 where Ql1 = °U l ° , Ql2 = °V l ° , l = 0, 1, 2, ..., n, 0 ≤ n ≤ N . 1 Proof: Computing the inner product of (11) with U n+ 2 , we get that ¿ n+1 À ¶ M −1 M −1 µ X X U − U n U n+1 + U n n+ 1 n+ 1 i , − γh1−α cj−l Ul 2 Uj 2 τ 2 j=1 l=1

+h

M −1 X

µ ¶ n+ 1 n+ 1 aj Uj 2 Uj 2 = 0,

j=1

5

(17)

ï ¯ ¯ ¯ ! ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 n where the real coefficients aj = ρ ¯¯Uj 2 ¯¯ + β ¯¯Vj 2 ¯¯ , and U n = (U1n , U2n , ..., UM −1 ), n V n = (V1n , V2n , ..., VM −1 ). Then taking the imaginary part and using Lemma 3.2, we have

° n+1 °2 °U ° = kU n k2 .

(18)

° n+1 °2 °V ° = kV n k2 .

(19)

Similarly, we can obtain

¤ When α = 2, this theorem is in accords with the results given in [13]. and can be seen as the discrete masses for the fractional coupled Schr¨ odinger equations. Theorem 3.1 implies that the scheme (11)-(12) is unconditionally stable for initial values. Moreover, if u0 (x), v0 (x) ∈ L2 [a, b], it follows from Theorem 3.1 that the numerical solution of (11)-(12) is bounded, i.e., there exists some constant C > 0, such that Qn1

kU n k ≤ C,

kV n k ≤ C,

Qn2

n = 0, 1, 2, ....

(20)

We now prove the existence of numerical solution. Lemma 3.3 (Brouwer fixed point theorem) Let (H, h·, ·i) be a finite dimensional inner product space, k · k be the associated norm, and g : H → H be continuous. Assume that there exists the real constant δ > 0 such that Rehg(x), xi ≥ 0, for ∀x ∈ H, kxk = δ. Then, there exists an element x∗ ∈ H and kx∗ k ≤ δ such that g(x∗ ) = 0. Let Z4 = {s|s = (s1 , s2 ), s1 , s2 ∈ Ωh,0 }, and define that hs, s0 i = hs1 , s01 i + hs2 , s02 i, ksk2 = ks1 k2 + ks2 k2 . Theorem 3.2 The solution of the difference scheme (11)-(12) exists. Proof. We prove this theorem by induction. It follows from (11)-(12) that (u0 , v 0 ) ∈ Z4 exists and satisfies the difference scheme. Assume that there exist (u0 , v 0 ), (u1 , v 1 ), ..., (un , v n ) ∈ Z4 which satisfy (11)-(12) for n ≤ N − 1. We now try to prove that there exists (un+1 , v n+1 ) ∈ Z4 which also satisfies this difference scheme. n+ 12

Since Ujn+1 = 2Uj (11)-(12) in the form n+ 1 Uj 2 n+ 1 Vj 2

=

Ujn

n+ 12

− Ujn and Vjn+1 = 2Vj

− Vjn , for fixed n, we rewrite the

ï ¯ ¯ ¯ ! 1 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 2 2 ¯U ¯ ¯ ¯ U n+ 2 , j ¯ j ¯ + β ¯Vj ¯ ï ¯ ¯ ¯ ! ¯ n+ 1 ¯2 iτ ρ ¯¯ n+ 12 ¯¯2 n+ 1 + Vj + β ¯¯Uj 2 ¯¯ Vj 2 . ¯ ¯ 2

M −1 iτ γ X iτ ρ n+ 1 − α cj−l Ul 2 + 2h 2

iτ γ = Vjn − α 2h

l=1 M −1 X l=1

n+ 1 cj−l Ul 2

6

(21)

n+ 12

Let (s1 )j = Uj

n+ 12

, (s2 )j = Vj

(g1 (s))j = (s1 )j − (U n )j +

and define the mapping g = (g1 , g2 ) on Z4 by

M −1 ´ iτ γ X iτ ρ ³ 2 2 c (s ) + |(s ) | + β |(s ) | (s1 )j , 1 j 2 j j−l 1 l 2hα 2

iτ γ (g2 (s))j = (s2 )j − (V n )j + α 2h

l=1 M −1 X l=1

´ iτ ρ ³ cj−l (s2 )l + |(s2 )j |2 + β |(s1 )j |2 (s2 )j , 2

(22)

where s = (s1 , s2 ). Using Lemma 3.2, we have hg1 (s), s1 i =ks1 k2 − hU n , s1 i +

+ ih

M −1 X





iτ γ Re  hα−1

cj−l (s1 )l (s1 )j 

M −1 M −1 X X j=1 l=1

(23)

bj (s1 )j (s1 )j ,

j=1

where the real coefficients bj =

τρ 2

³ ´ |(s1 )j |2 + β |(s2 )j |2 . Therefore,

Rehg1 (s), s1 i =ks1 k2 − hU n , s1 i. Similarly, we have Rehg2 (s), s2 i =ks2 k2 − hV n , s2 i. We have that Rehg(s), si = Rehg1 (s), s1 i + Rehg2 (s), s2 i = ks1 k2 − hU n , s1 i + ks2 k2 − hV n , s2 i ¢ 1¡ ≥ ks1 k2 + ks2 k2 − ks1 k2 + ks2 k2 + kU n k2 + kV n k2 2 ¢ 1¡ = k(s1 , s2 )k2 − k(U n , V n k2 . 2

(24)

p Hence, taking δ = k(U n , V n )k2 + 1, for ∀s : ||s|| = δ, we have Rehg(s), si ≥ 12 . By Lemma 3.3, there exists an element s∗ ∈ Z4 such that g(s∗ ) = 0, which completes the proof. ¤ To prove the convergence, we need the following lemmas. Lemma 3.4 (Young Inequality) If a ≥ 0, b ≥ 0, then there is the inequality ab ≤

ap bq + , p q

or ab ≤

(εa)p (b/ε)q + , p q

where p > 1, q > 1, ε > 0, and p1 + 1q = 1. Lemma 3.5 (Gronwall Inequality) Assume that {Gn |n ≥ 0} is a nonnegative sequence, and satisfies Gn+1 ≤ (1 + cτ )Gn + τ σ, 7

n = 0, 1, 2, ...,

where c and σ are nonnegative constants. Then Gn satisfies ¡ ¢ Gn+1 ≤ ecnτ G0 + σ/c . Lemma 3.6 [15] For any complex functions U, V, u, v, one has ¯ 2 ¯ ¯|U | V − |u|2 v ¯ ≤ (max{|U |, |V |, |u|, |v|})2 · (2|U − u| + |V − v|) . Theorem 3.3 Suppose that the original problem (5)-(7) has a smooth solution. Then the the numerical solution (U n , V n ) of the difference scheme (11)-(12) is convergent to the ture solution (u(x, t), v(x, t)) with the error O(τ 2 + h2 ) in the L2 norm when 0 < τ ≤ τ0 , τ0 is a constant. Proof. Let unj = u(xj , tn ), vjn = v(xj , tn ) and enj = unj − Ujn ,

fjn = vjn − Vjn ,

j = 0, 1, ..., M, n = 0, 1, ..., N,

and ï ¯ ¯ ¯ ! M −1 X 1 1 ¯2 1 ¯2 ¯ ¯ ¡ ¢ γ n+ n+ n+ n+ 1 pnj = i unj t − α cj−l ul 2 + ρ ¯¯uj 2 ¯¯ + β ¯¯vj 2 ¯¯ uj 2 , h l=1 ï ¯ ¯ ¯ ! M −1 X 1 1 ¯2 1 ¯2 ¯ ¯ ¡ ¢ γ n+ n+ n+ n+ 1 qjn = i vjn t − α cj−l vl 2 + ρ ¯¯vj 2 ¯¯ + β ¯¯uj 2 ¯¯ vj 2 . h

(25)

l=1

It follows from Lemma 3.1 and Taylor’s expansion that there exists constant C1 > 0 such that |pnj | ≤ C1 (τ 2 + h2 ),

|qjn | ≤ C1 (τ 2 + h2 ).

(26)

From (11) and (24), we have M −1 ¡ ¢ γ X n+ 1 n+ 1 pnj = i enj t − α cj−l el 2 + Gj 2 , h

¡ ¢ γ qjn = i fjn t − α h

l=1 M −1 X

(27) n+ 1 cj−l fl 2

+

n+ 1 Fj 2 ,

l=1

where n+ 1 Gj 2 n+ 1 Fj 2

ï ï ¯ ¯ ¯ ! ¯ ¯ ¯ ! 1 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 n+ 12 2 2 2 2 ¯ U n+ 2 , = ρ ¯¯uj ¯¯ + β ¯¯vj ¯¯ uj − ρ ¯¯Uj ¯¯ + β ¯¯Vj j ¯ ï ï ¯ ¯ ¯ ! ¯ ¯ ¯ ! 1 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 ¯ n+ 1 ¯2 n+ n+ 1 = ρ ¯¯vj 2 ¯¯ + β ¯¯uj 2 ¯¯ vj 2 − ρ ¯¯Vj 2 ¯¯ + β ¯¯Uj 2 ¯¯ Vj 2 .

8

(28)

° ° ° n+ 1 °2 2° ° We now estimate °Gj ° and

° ° ° n+ 1 °2 °F 2 ° . By Lemma 3.6, we have ° j °

ï ! ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ n+ 1 ¯ ¯ n+ 1 ¯2 n+ 1 ¯ n+ 1 ¯2 n+ 1 ¯ n+ 1 ¯2 n+ 1 ¯ n+ 1 ¯2 n+ 1 ¯G 2 ¯ = |ρ| ¯u 2 ¯ u 2 − ¯U 2 ¯ U 2 + β ¯v 2 ¯ u 2 − β ¯V 2 ¯ U 2 ¯ j ¯ ¯ j ¯ j ¯ j ¯ j ¯ j ¯ j ¯ j ¯ j à µ ¯ ½ ¾¶2 ¯ ¯ n+ 1 ¯ n+ 12 n+ 12 2¯ ¯ ≤ |ρ| 3 · max |Uj |, |uj | ¯ej ¯ + ¯ ¯ ¯¶! µ ½ ¾¶2 µ ¯ ¯ n+ 1 ¯ ¯ n+ 1 ¯ n+ 12 n+ 12 n+ 12 n+ 12 2 2 |β| · max |Vj |, |Uj |, |uj |, |vj | · 2 ¯¯fj ¯¯ + ¯¯ej ¯¯ . (29) It follows from (20) and (29) that there exists constant C2 > 0 such that ¯ ¯ ¯ ¯ ¯¶ µ¯ ¯ n+ 1 ¯ ¯ n+ 1 ¯ ¯ n+ 1 ¯ 2 2 2 ¯G ¯ ¯ ¯ ¯ ¯ ¯ j ¯ ≤ C2 ¯ej ¯ + ¯fj ¯ , which implies that there exists constant C3 > 0 such that µ° ° ° ° ° ° ¶ ° n+ 12 °2 ° n+ 12 °2 ° n+ 12 °2 °G ° ≤ C3 °e ° + °f ° .

(30)

In the same way, we have that µ° ° ° ° ° ° ¶ ° n+ 12 °2 ° n+ 12 °2 ° n+ 12 °2 °F ° ≤ C3 °e ° + °f ° .

(31)

D E Computing the inner product pn , en+1 + en and taking the imaginary part, and using Lemma 3.2, we have D E 1¡ D E ¢ 1 1 Im pn , en+1 + en = ken+1 k2 − ken k2 + Im Gn+ 2 , 2en+ 2 , τ 1

1

1

1

where en = (en1 , en2 , ..., enM −1 ), and en+1 , pn , en+ 2 , Gn+ 2 F n+ 2 and f n+ 2 are defined similarly. Hence, D E D E ¢ 1 1 1 1 ¡ n+1 2 ke k − ken k2 = Im pn , 2en+ 2 + Im −Gn+ 2 , 2en+ 2 , τ µ ° ° ¶ µ° ° ° ° ¶ ° n+ 12 °2 ° n+ 12 °2 ° n+ 12 °2 n 2 ≤ kp k + °e ° + °G ° + °e ° .

(32)

In the same way, we have that D E D E ¢ 1 1 1 1 ¡ n+1 2 kf k − kf n k2 = Im q n , 2f n+ 2 + Im F n+ 2 , 2f n+ 2 , τ µ ° ° ¶ µ° ° ° ° ¶ 1 °2 1 °2 1 °2 ° ° ° ≤ kq n k2 + °f n+ 2 ° + °F n+ 2 ° + °f n+ 2 ° .

9

(33)

Adding (32) and (33) leads to ¢ 1 ¡ n+1 2 (ke k + kf n+1 k2 ) − (ken k2 + kf n k2 ) τ µ° ° ° ° ¶ µ° ° ° ° ¶ ³ ´ ° n+ 12 °2 ° n+ 12 °2 ° n+ 12 °2 ° n+ 12 °2 n 2 n 2 ≤ kp k + kq k + 2 °e ° + °f ° + °F ° + °G ° µ° ° ° ° ¶ ¡ ¢2 1 °2 1 °2 ° ° ≤ (2 + 2C3 ) °en+ 2 ° + °f n+ 2 ° + 2(b − a)C12 τ 2 + h2

(34)

¡ ¢ ¡ ¢2 ≤ (1 + C3 ) ken+1 k2 + kf n+1 k2 + ken k2 + kf n k2 + 2(b − a)C12 τ 2 + h2 . Therefore, there exists the constant τ0 : 0 < τ0 < 1/(1 + C3 ) such that when 0 < τ ≤ τ0 , we have µ ¶ ¡ n 2 ¢ ¢2 2τ (1 + C3 ) 2τ (b − a)C12 ¡ 2 n+1 2 n+1 2 ke k + kf k ≤ 1+ ke k + kf n k2 + τ + h2 . 1 − τ0 (1 + C3 ) 1 − τ0 (1 + C3 ) Let C4 =

2(1+C3 ) 1−τ0 (1+C3 )

and C5 =

2(b−a)C12 1−τ0 (1+C3 ) ,

we have

¡ ¢ ¡ ¢2 ken+1 k2 + kf n+1 k2 ≤ (1 + C4 τ ) ken k2 + kf n k2 + C5 τ τ 2 + h2 .

(35)

By the Gronwall inequality, we obtain that ken+1 k2 + kf n+1 k2 ≤

¢2 C5 C 4 T ¡ 2 e τ + h2 , C4

n = 0, 1, 2, ..., N.

(36)

This completes the proof. ¤ From the proofs of Theorems 3.1, 3.2 and 3.3, we know that the Lemma 3.1 plays an important role for fractional coupled nonlinear Schr¨ odinger equations. The symmetry of the coefficients of numerical approximation to fractional derivatives, i.e., ck = c−k , k ≥ 1, is often crucial for numerical analysis.

4

Examples and numerical results Example 4.1 Let β = 0. Then the system (5) is decoupled and becomes α

iut − γ(−∆) 2 u + ρ|u|2 u = 0

(37)

u(x, 0) = sech(x) · exp(2ix).

(38)

subject to the initial value

For this problem, we take the parameters γ = 1, ρ = 2. When α = 2, the exact solution is given by u(x, t) = sech(x − 4t) · exp(i(2x − 3t)).

10

Figure 1. Numerical solution for α = 1.4

Figure 2. Numerical solution for α = 1.6

Figure 3. Numerical solution for α = 1.8

Figure 4. Numerical solution for α = 1.999

1.4

1.2

1

0.8

0.6

0.4

t=1 0.2

t=2 t=3

0 −20

−15

−10

−5

0

5

10

15

20

Figure 5. Numerical solution for α = 2 Figure 6. Numerical solution of time t = 1, 2, 3 for α = 1.8 The initial-boundary value problem (5)-(7) can be considered to the initial value problem (37)-(38) and for a ¿ 0 and b À 0. Since the initial value u(x, 0) exponentially decays 11

Table 1. The errors of ||u − uh ||L2 and the order for α = 2 h 0.1 0.05 0.025 0.0125

dt 0.1 0.05 0.025 0.0125

||u − uh || 0.2321895347117 0.0588407677959 0.0146381849150 0.0036626158838

order 1.9804 2.0071 1.9988

¯ t ¯ ¯ n t0 ¯ Table 2. The errors of ¯ Q Q−Q ¯ for different values of α and h = 0.05, τ = 0.05 t0

α=1.4 α=1.6 α=1.8 α=1.99 α=2

t=1 9.8014e-08 1.4949e-07 7.6308e-08 5.0251e-10 5.1593e-10

t=2 3.7045e-07 1.9680e-07 7.8807e-08 1.2179e-09 1.1654e-09

t=3 1.5930e-07 9.6761e-08 8.0580e-08 1.8880e-09 1.8535e-09

to zero with the variable x away from the origin, the wave function can be negligible outside the interval [a, b] and we can set u(a, t) = u(b, t) = 0 for a ¿ 0 and b À 0. In this example, we chose a = −20 and b = 20. In Figures 1 to 6, we take h = τ = 0.05. The figures present the numerical solutions for different values of order α. We observe that the order α will affect the shape of the soliton. It will not only change the height and width of the solitary wave solution but produce two turning points. See the details in Figure 6. When α becomes smaller, the shape of the soliton will change more quickly. This property of the fractional Schr¨ odinger equations can be used in physics to modify the shape of wave without change of the nonlinearity and dispersion effects. When α tends to 2, the numerical solutions of the fractional equation are convergent to the solutions of the usual classical non-fractional equation. Table 1 shows that the method has second order accuracy for α = 2, as expected. Table 2 shows that the method preserve the discrete masses for different values of α. Example 4.2 For the coupled system β 6= 0, α

iut − γ(−∆) 2 u + ρ(|u|2 + β|v|2 )u = 0, α

ivt − γ(−∆) 2 v + ρ(|v|2 + β|u|2 )v = 0,

(39)

we take the initial conditions in the form u(x, 0) = sech(x + D0 ) · exp(iV0 x), v(x, 0) = sech(x − D0 ) · exp(−iV0 x).

12

(40)

When α = 2 and β = 1, the system is the Manakov model which is completely integrable and hence the collision is elastic, see Figure 7. In this example, we fix the parameters D0 = 10, V0 = 3, γ = 1, ρ = 2 and take h = 0.2, τ = 0.1, a = −20, b = 20.

1 |u|, |v| 0.5 0 5

t 4

3

2

1

0

20

15

10

5

0 x

−5

−10

−15

−20

Figure 7. Numerical solution for α = 2, β = 1

Figure 8. Numerical solution for α = 1.4, β = 1

Figure 9. Numerical solution for α = 1.5, β = 3

Figure 10. Numerical solution for α = 1.5, β = 2/3

¯ ¯ ¯ Qtn −Qt0 ¯ 1 1 ¯ ¯ Table 3. The errors of ¯ t ¯ for different values of α and β Q10

α=1.4, β=1 α=1.5, β = 23 α=1.5, β=3

t=2 4.9609e-07 3.4084e-07 3.4084e-07

t=4 1.4035e-06 9.0982e-07 4.7427e-07

t=6 9.3659e-07 1.9468e-06 6.9829e-07

t=8 1.8240e-06 2.2889e-06 1.5004e-06

t = 10 3.2615e-06 2.9502e-06 2.1962e-06

As in Example 4.1, the initial-boundary value problem (5)-(7) can be considered to initial value problem (39)-(40) and we can set u(a, t) = u(b, t) = 0 for a ¿ 0 and b À 0. 13

The Figures 7 to 12 present the numerical solutions for different values of order α and β. We can find that that the order α will affect the shape of the solitons. When β = 1 and α 6= 2, the collision of solitons are not elastic any more. When the α becomes smaller, the shape of the solitons will change more quickly. When α tends to 2, the numerical solutions of the fractional equation are also convergent to the solutions of the usual classical nonfractional equation. Table 3 shows that the method preserve the discrete masses for different values of α and β.

5

Acknowledgements

The authors would like to thank the anonymous referees for very useful comments and suggestions. This work is supported by projects NSF of China (11271311,10971175), Program for Changjiang Scholars and Innovative Research Team in University of China (IRT1179), NSF of Hunan Province (10JJ7001), the Aid Program for Science and Technology, Innovative Research Team in Higher Educational Institutions of Hunan Province of China, and Hunan Province Innovation Foundation for Postgraduate (CX2011B244).

References [1] N. Laskin, Fractional quantum mechanics, Phys. Rev. E 62(2000)3135-3145. [2] N. Laskin, Fractional quantum mechanics and L´evy path integrals, Phys. Lett. A 268(4-6)(2000)298-305. [3] X. Guo, M. Xu, Some physical applications of fractional Schr¨ odinger equation, J. Math. Phys. 47(2006)082-104. [4] J. Fujioka, A. Espinosa, R.F. Rodr´ıguez, Fractional optical solitons, Phy. Lett. A 374(2010)1126-1134. [5] P. Amore, F. M. Fern´andez, C. P. Hofmann, and R. A. S´a enz, Collocation method for fractional quantum mechanics, J. Math. Phys. 51(2010)122101. [6] S. Leble and B. Reichel, Coupled nonlinear Schr¨ odinger equations in optic fibers theory: From general to solitonic aspects, Eur. Phys. J. Special Topics 173(2009)555. [7] W. Bao, D. Jaksch, An explicit unconditionally stable numerical method for solving damped nonlinear Schr¨odinger equations with a focusing nonlinearity, SIAM J. Numer. Anal. 41(2003)1406-1426. [8] V.E. Tarasov, Fractional dynamics: Application of fractional calculus to dynamics of particles, fields and media, Higher Education Press and Springer, 2010.

14

[9] R. Hilfer, Applications of fractional calculus in physics, World Scientific, Singapore, 1999. [10] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and applications of fractional differential equations, Elsevier, Amsterdam, 2006. [11] B. Guo, Y. Han, J. Xin, Existence of the global smooth solution to the period boundary value problem of fractional nonlinear Schr¨ odinger equation, Appl. Math. Comput. 204(2008)468-477. [12] J. Hu, J. Xin, H. Lu, The global solution for a class of systems of fractional nonlinear Schr¨odinger equations with periodic boundary condition, Computers Math. Appl. 62(2011)1510-1521. [13] T. Wang, T. Nie, L. Zhang, Analysis of a symplectic difference scheme for a coupled nonlinear Schr¨odinger system, J. Comput. Appl. Math. 231(2)(2009)745-759. [14] T. Wang, B. Guo, L. Zhang, New conservative difference schemes for a coupled nonlinear Schr¨odinger system, Appl. Math. Comput. 217(2010)1604-1619. [15] Z. Sun, D. Zhao, On the L∞ convergence of a difference secheme for coupled nonlinear Schr¨odinger equations, Computers Math. Appli. 59(2010)3286-3300. [16] M.S. Ismail, A fourth-order explicit schemes for the coupled nonlinear Schr¨ odinger equation, Appl. Math. Comput. 196(2008)273-284. [17] F. Ivanauskas, M. Radˇzi¯ unas, On convergence and stability of the explicit difference method for solution of nonlinear Schr¨ odinger equations, SIAM J. Numer. Anal. 36(1999)1466-1481. [18] L. Gauckler, C. Lubich, Nonlinear Schr¨ odinger equations and their spectral discretizations over long times, Found. Comput. Math. 10(2010)141-169. [19] J. Sun, M. Qin, Multi-symplectic methods for the coupled 1D nonlinear Schr¨ odinger system, Comput. Phys. Commun. 155(2003)221-235. [20] A. Aydin, B. Karas¨ozen, Multi-symplectic integration of coupled nonlinear Schr¨odinger system with soliton solutions, Int. J. Computer Math. 86(2009)864-882. [21] P. Zhuang, F. Liu, V. Anh, I. Turner, Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term, SIAM J. Numer. Anal. 47(2009)1760-1781. [22] Q. Yang, F. Liu, I. Turner, Numerical methods for fractional partial differential equations with Riesz space fractional derivatives, Appl. Math. Model. 34(1)(2010)200-218. [23] M. Ili´c, F. Liu, I. Turner and V. Anh, Numerical approximation of a fractional-inspace diffusion equation, I. Fractional Calculus & Applied Analysis 8(3)(2005)323-341. 15

[24] S. Shen, F. Liu, V. Anh, I. Turner and J. Chen, A novel numerical approximation for the Riesz space fractional advection-dispersion equation, IMA J. Appl. Math. (2012), in press, doi:10.1093/imamat/hxs073. [25] Q. Yang, I. Turner, F. Liu and Milos Ilis, Novel numerical methods for solving the time-space fractional diffusion equation in 2D, SIAM J. Sci. Comput. 33(2011)11591180. [26] H. Zhang, F. Liu, V. Anh, Galerkin finite element approximations of symmetric spacefractional partial differential equations, Appl. Math. Comp. 217(6)(2010), 2534-2545. [27] M.D. Ortigueira, Riesz potential operators and inverses via fractional centred derivatives, Int. J. Math. Math. Sci. (2006)1-12. [28] C. C ¸ elik, M. Duman, Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys. 231(2012)1743-1750. [29] Q. Yu, F. Liu, I. Turner and K. Burrage, A computationally effective alternating direction method for the space and time fractional Bloch-Torrey equation in 3-D, Appl. Math. Comp. 219(2012)4082-4095. [30] S. Shen, F. Liu and V. Anh, Numerical approximations and solution techniques for the space-time Riesz-Caputo fractional advection-diffusion equation, Numerical Algorithm, 56(3)(2011)383-404. [31] H. Zhang and F. Liu, Numerical simulation of the Riesz fractional diffusion equation with a nonlinear source term, J. Appl. Math. Informatics, 26(1-2)(2008)1-14. [32] S. Shen, F. Liu, V. Anh and I. Turner, The fundamental solution and numerical solution of the Riesz fractional advection-dispersion equation, IMA J. Appl. Math. 73(2008)850-872.

16