Accepted Manuscript A fully spectral collocation approximation for multi-dimensional fractional Schrödinger equations A.H. Bhrawy, M.A. Abdelkawy
PII: DOI: Reference:
S0021-9991(15)00224-7 http://dx.doi.org/10.1016/j.jcp.2015.03.063 YJCPH 5819
To appear in:
Journal of Computational Physics
Received date: 9 January 2015 Revised date: 25 March 2015 Accepted date: 26 March 2015
Please cite this article in press as: A.H. Bhrawy, M.A. Abdelkawy, A fully spectral collocation approximation for multi-dimensional fractional Schrödinger equations, J. Comput. Phys. (2015), http://dx.doi.org/10.1016/j.jcp.2015.03.063
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.
A fully spectral collocation approximation for multi-dimensional fractional Schr¨ odinger equations A.H. Bhrawy 1
1,2
, M.A. Abdelkawy
2
Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia E-mail:
[email protected]
2
Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt E-mail:
[email protected]
Abstract: A shifted Legendre collocation method in two consecutive steps is developed and analyzed to numerically solve time fractional one- and two-dimensional Schr¨ odinger equations (TFSEs) subject to initial-boundary and non-local conditions. The first step depends mainly on shifted Legendre Gauss-Lobatto collocation (SL-GL-C) method for spatial discretization; an expansion in a series of shifted Legendre polynomials for the approximate solution and its spatial derivatives occurring in the TFSE is investigated. In addition, the Legendre-Gauss-Lobatto quadrature rule is established to treat the nonlocal conservation conditions. Thereby, the expansion coefficients are then determined by reducing the TFSE with its nonlocal conditions to a system of fractional differential equations (SFDEs) for these coefficients. The second step is to propose a shifted Legendre Gauss-Radau collocation (SL-GR-C) scheme, for temporal discretization, to reduce such system into a system of algebraic equations which is far easier to be solved. The proposed collocation scheme, both in temporal and spatial discretizations, is successfully extended to solve the two-dimensional TFSE. Numerical results carried out to confirm the spectral accuracy and efficiency of the proposed algorithms. By selecting relatively limited Legendre Gauss-Lobatto and Gauss-Radau collocation nodes, we are able to get very accurate approximations, demonstrating the utility and high accuracy of the new approach over other numerical methods. Keyword: Fractional Schr¨odinger equations; Two-dimensional Schr¨odinger equations; Collocation method; Spectral method; Gauss-type quadrature.
1
Introduction
In recent years, spectral methods (see [1]–[7]) are often efficient and highly accurate schemes when compared with local methods. The speed of convergence is one of the great advantages of spectral methods. Besides, spectral methods have exponential rates of convergence; they also have high level
1
of accuracy. The main idea of all versions of spectral methods is to express the approximate solution of the problem as a finite sum of certain basis functions (orthogonal polynomials or combination of them) and then choose the coefficients in order to minimize the difference between the exact and approximate solutions as well as possible. The spectral collocation method is a specific type of spectral methods, that is more applicable and widely used to solve almost types of differential equations [8]-[10]. Fractional differential equations (see for instance, [11]-[17]) are used to model many phenomena in several fields such as fluid mechanics, chemistry [14, 15], biology [16], viscoelasticity [17], engineering, finance and physics [18] fields. Existence and uniqueness of solutions related to the fractional differential equations have been investigated by many authors (see for example [19]-[21]). Most of fractional differential equations do not have exact analytic solutions. Therefore, there has been a growing interest in the last decades in developing numerical methods for solving fractional differential equations. Several numerical treatment based on finite difference methods for fractional differential equations are given in [22]-[24]. Also, efficient spectral algorithms ([25]-[34]) are designed and developed for solving different kinds of fractional differential equation. The classical Schr¨ odinger equations occur in various area of physics, including nonlinear optics, plasma physics, superconductivity and quantum mechanics, and can be written in the following form
i
∂ψ(x , t) + ∇2 ψ(x , t) + γ|ψ(x , t)|2 ψ(x , t) − δR(x , t) = 0, ∂t
(1.1)
where, ψ(x , t) and R(x , t) are complex functions, at the space x , and time t. In particular equation (1.1) appears in nonlinear optics in the context of soliton propagation through optical fibers. In 3-dimensions this equation shows up in the study of optical bullets. Here, the first term represents the temporal evolution of the pulses that will be considered to be fractional in this paper. The second term is the group velocity dispersion that accounts for the dispersive effect of the solitons. The coefficient of γ represents Kerr law nonlinearity that accounts for intensity dependent refractive index of light. Finally, δ represents the coefficient of the driven term. Recently, the analytical and numerical solutions of different types of the previous classical Schr¨ odinger equations are discussed in [35, 36]. Here, we focus on the application of SL-GL-C and SL-GR-C methods in two consecutive steps for providing a high accurate numerical solution of the TFSEs
i
∂ α ψ(x , t) + ∇2 ψ(x , t) + γ|ψ(x , t)|2 ψ(x , t) − δR(x , t) = 0, ∂tα 2
(1.2)
where
∂ α ψ(x ,t) ∂tα
is the fractional derivative of order α ( 0 < α ≤ 1), in the Caputo sense. The
TFSE is equivalent to the standard one but with Hamiltonian time dependent [37]. For numerical solution of TFSEs, see for example the articles [38]–[41]. From the overview of approximation of the underline problem. The Adomian decomposition technique [38] is proposed and used to find the solution of TFSE and is extended in [39] to solve space-time fractional Schr¨odinger equations. Wei et al. [40] proposed and developed an efficient algorithm based on local discontinuous Galerkin scheme to approximate the solution of TFSE. Recently, Mohebbi et al. [41] proposed a combination of a collocation framework and radial basis functions for providing an accurate and efficient numerical solution of the TFSE, and they also developed this method to improve the accuracy of the algorithm greatly for solving nonlinear sineGordon and Klein-Gordon equations with fractional orders [42]. For other numerical schemes for fractional Schr¨odinger equations, the interested reader is referred to ([43]-[48] and the references therein). The main objective of this paper is to extend and develop the collocation method to solve one- and two-dimensional fractional Schr¨ odinger equations subject to initial-boundary and nonlocal conditions. The proposed collocation scheme is investigated for both temporal and spatial discretizations. Firstly, the SL-GL-C is proposed, with a suitable modification for treating the boundary or nonlocal conditions, for spatial discretization. This treatment improves the accuracy of the scheme greatly. Therefore, the TFSE with its nonlocal conditions is reduced to SFDEs subject to a vector of initial values. Secondly, the SL-GR-C is then investigated for temporal discretization, which is more reasonable for solving initial value problems. Thereby, the problem is reduced to a system of algebraic equations which is far easier to be solved. In addition, this algorithm is developed to numerically solve the two-dimensional TFSEs. Finally, several numerical examples with comparisons lighting the high accuracy and effectiveness of the proposed algorithm are presented. This paper is organized as follows. we present few preliminaries and some facts about shifted Legendre polynomials in Section 2. Section 3 is regulated into three subsections, the first one proposes a new collocation method for the one-dimensional TFSE subject to initial-boundary conditions, while subsections 3.2 and 3.3 are devoted to the numerical solution for the TFSE subject to non-local and mixed conditions. In Section 4, the proposed scheme is successfully extended to solve the two-dimensional version of TFSE. In Section 5, we propose the SL-GR-C scheme to solve SFDEs. Section 6 is devoted to solve five test problems. Finally, some concluding remarks are given in the last section.
3
2
Properties of shifted Legendre polynomials
There are several definitions of the fractional integration of order ν > 0, and not necessarily equivalent to each other, see [49]. Riemann-Liouville and Caputo fractional definitions are the two most used from all the other definitions of fractional calculus which have been introduced recently. Definition 2.1. The fractional integral of order ν ≥ 0 according to Riemann-Liouville sense is given by J ν f (x) =
1 Γ(ν)
x 0
(x − ζ)
ν−1
f (ζ)dζ,
ν > 0,
x > 0, (2.1)
0
J f (x) = f (x), where
∞
Γ(ν) =
xν−1 e−x dx
0
is the gamma function. The operator J ν satisfies the following properties J ν J μ f (x) = J ν+μ f (x), J ν J μ f (x) = J μ J ν f (x), J ν xβ =
(2.2)
Γ(β + 1) xβ+ν . Γ(β + 1 + ν)
Definition 2.2. The Caputo fractional derivatives of order ν is defined as Dν f (x) =
1 Γ(m − ν)
x 0
(x − ζ)
m−ν−1
dm f (ζ)dζ, dtm
m − 1 < ν ≤ m, x > 0,
(2.3)
xi , i!
(2.4)
where m is the ceiling function of ν. The Caputo fractional derivative operator satisfies J ν Dν f (x) = f (x) −
m−1
f (i) (0+ )
i=0
D ν xβ =
⎧ ⎪ ⎪ ⎪ ⎨0,
for β ∈ N0 and β < ν,
⎪ ⎪ ⎪ ⎩
for β ∈ N0 and β ≥ ν or β ∈ N and β > ν ,
Γ(β + 1) xβ−ν , Γ(β + 1 − ν)
(2.5)
where ν and ν are the ceiling and floor functions respectively, while N = {1, 2, . . .} and N0 = {0, 1, 2, . . .}. The well-known Legendre polynomials Pi (x) are defined on the interval [−1, 1]. Firstly, some properties about the standard Legendre polynomials have been recalled in this section. The Legendre polynomials Pk (x) (k = 0, 1 . . . , ) satisfy the following Rodrigue’s formula Pk (x) =
(−1)k k D ((1 − x2 )k ), 2k k! 4
(2.6)
(q)
we recall also that Pk (x) is a polynomial of degree k and therefore Pk (x) (the qth derivative of Pk (x)) is given by k−q
(q)
Pk (x) =
Cq (k, i)Pi (x),
(2.7)
i=0(k+i=even)
where 2q−1 (2i + 1)Γ( q+k−i )Γ( q+k+i+1 ) 2 2
Cq (k, i) =
Γ(q)Γ( 2−q+k−i )Γ( 3−q+k+i ) 2 2
.
The Legendre polynomials [1] satisfy the following relations P0 (x) = 1,
P1 (x) = x,
Pk+2 (x) =
2k + 3 k+1 xPk+1 w (x) − Pk (x), k+2 k+2
and the orthogonality relation 1 (Pk (x), Pl (x))w =
Pk (x)Pl (x) w(x) = hk δlk ,
(2.8)
−1
where w(x) = 1, hk =
2 2k+1 .
The Legendre-Gauss-Lobatto (SL-GL) quadrature is commonly used
to evaluate the previous integrals accurately. For any φ ∈ S2N −1 [−1, 1], we have 1 φ(x)dx =
N
N,j φ(xN,j ),
(2.9)
j=0
−1
where SN [−1, 1] is the set of polynomials of degree less than or equal to N . While the discrete inner product is given as (u, v)w =
N
u(xN,j ) v(xN,j ) N,j ,
(2.10)
j=0
where xN,j (0 ≤ j ≤ N ) and N,j (0 ≤ j ≤ N ) are used as usual the nodes and the corresponding Christoffel numbers in the interval [−1, 1], respectively. For Legendre Gauss-Lobatto, it is known that [1] xN,0 = −1,
xN,N = 1,
xN,j (j = 1, · · · , N − 1) N,j =
are the zeros of ∂x PN (x),
2 , N (N + 1)(PN (xN,j ))2
(2.11)
while the nodes and the corresponding Christoffel numbers in Legendre Gauss-Radau quadrature are given by [1] xN,j (j = 0, · · · , N ) N,0 =
2 , (N + 1)2
are the zeros of (PN (x)) + (PN +1 (x)), N,j =
1 − xj . (N + 1)2 (PN (xj ))2
(2.12)
In order to use these polynomials on the interval [0, L], we define the so-called shifted Legendre 2x −1. Let the shifted Legendre polynomials polynomials by introducing the change of variable x = L
5
2x − 1) be denoted by PL,i (x). Then PL,i (x) can be obtained with the aid of the following L recurrence formula: Pi (
(i + 1)PL,i+1 (x) = (2i + 1)(
2x − 1)PL,i (x) − iPL,i−1 (x), L
i = 1, 2, · · · .
(2.13)
The analytic form of the shifted Legendre polynomials PL,i (x) of degree i is given by PL,i (x) =
i
(−1)
i+k
k=0
(i + k)! xk , (i − k)! (k!)2 Lk
(2.14)
and the orthogonality condition is
L 0
PL,j (x)PL,k (x)wL (x)dx = L k δjk ,
(2.15)
L . 2k + 1 A function u(x), square integrable in [0, L], may be expressed in terms of shifted Legendre
where wL (x) = 1 and L k =
polynomials as u(x) =
∞
cj PL,j (x),
j=0
where the coefficients cj are given by cj =
1 L j
L 0
u(x)PL,j (x)wL (x)dx,
j = 0, 1, 2, · · · .
(2.16)
In practice, only the first (N + 1)-terms shifted Legendre polynomials are considered. Hence u(x) can be expressed in the form uN (x)
N
cj PL,j (x).
(2.17)
j=0
3
One-dimensional TFSE
In this section, we present three numerical algorithms that based on L-GL-C method to numerically solve TFSEs with different types of boundary conditions. The collocation points are selected at the (SL-GL) interpolation nodes. The core of the proposed method is to discretize the TFSE in the spatial direction along with a new treatment for the subjected conditions, to create a system of SFDEs of the unknown coefficients of spectral expansion in time direction.
3.1
TFSE with initial-boundary conditions
In particular, we consider the general initial-boundary value problem i
∂ α ψ(x, t) ∂ 2 ψ(x, t) + + γ|ψ(x, t)|2 ψ(x, t) − δR(x, t) = 0, ∂tα ∂x2
0 < α < 1,
(x, t) ∈ [0, L] × [0, T ], (3.1)
6
with the initial-boundary conditions ψ(0, t) =ζ1 (t),
ψ(L, t) = ζ2 (t),
ψ(x, 0) =ξ1 (x),
x ∈ [0, L],
t ∈ [0, T ],
(3.2)
where α represents order of fractional derivative, and 0 < α < 1. In order to begin the numerical solution, we put the complex functions (ψ(x, t), ζ1 (t), ζ2 (t), and ξ1 (x)), as proposed in [35], in their real and imaginary parts: ψ(x, t) = u(x, t) + i v(x, t), ζ2 (t) = g2 (t) + i g4 (t),
R(x, t) = f (x, t) + i g(x, t),
ζ1 (t) = g1 (t) + i g3 (t),
(3.3)
ξ1 (t) = χ1 (x) + i χ2 (x),
where u(x, t), v(x, t), f (x, t), g(x, t), g1 (t), g3 (t), g2 (t), g4 (t), χ1 (x) and χ2 (x), are real functions. In virtue of Eq. (3.3), one may write Eq. (3.1) as ∂ α u(x, t) ∂ 2 v(x, t) 2 2 + + γ(u (x, t) + v (x, t))v(x, t) − δg(x, t) i ∂tα ∂x2 (3.4) ∂ α v(x, t) ∂ 2 u(x, t) 2 2 − − γ(u (x, t) + v (x, t))u(x, t) + δf (x, t) = 0. − ∂tα ∂x2 The above equation can be separated into coupled time-dependent nonlinear fractional partial differential equations ∂ α u(x, t) ∂ 2 v(x, t) + + γ(u2 (x, t) + v 2 (x, t))v(x, t) − δg(x, t) = 0, ∂tα ∂x2 ∂ α v(x, t) ∂ 2 u(x, t) − − γ(u2 (x, t) + v 2 (x, t))u(x, t) + δf (x, t) = 0, ∂tα ∂x2 with the initial-boundary conditions, namely u(0, t) = g1 (t),
u(L, t) = g2 (t),
t ∈ [0, T ],
v(0, t) = g3 (t),
v(L, t) = g4 (t),
t ∈ [0, T ],
u(x, 0) = f1 (x),
v(x, 0) = f2 (x),
(3.5)
(3.6)
x ∈ [0, L].
The main advantage of using the node points of SL-GL quadrature is that the distribution of these points in [0, L]. In the proposed collocation scheme, unlike the most existing collocation schemes, the boundary and nonlocal conservation conditions are satisfied automatically in the collocation scheme. More specifically, no need to additional equations to enforce these conditions. Thereby, the above treatment is more reasonable. In addition, it improves the accuracy of numerical solution, see the numerical results in Section 6. Now, we outline the main step of implementing SL-GL-C scheme for reducing nonlinear coupled system (3.5)-(3.6) to SFDEs. The approximate solutions of u(x, t) and v(x, t) can be expanded, using the shifted Legendre polynomial, PL,j (x), in the forms u(x, t) =
N
aj (t)PL,j (x),
j=0
v(x, t) =
N
(3.7) bj (t)PL,j (x),
j=0
7
and in virtue of (2.15)-(2.16), we deduce that aj (t) =
L
1 L j
1 bj (t) = L j
u(x, t) wL (x) PL,j (x)dx, 0
(3.8)
L v(x, t) wL (x) PL,j (x)dx. 0
We can circumvent the need for evaluating the previous integrals, by using SL-GL quadrature formula to approximate the integrals. For any φ ∈ S2N −1 [0, L], L w(x)φ(x)dx =
N
L N,j φ(xL N,j ),
(3.9)
j=0
0
L where xL N,j (0 ≤ j ≤ N ) and N,j (0 ≤ j ≤ N ) are the nodes and the corresponding Christoffel
numbers of the SL-GL quadrature formula on the interval [0, L], respectively. Due to the property of SL-GL quadrature (3.9), the expansion coefficients aj (t) and bj (t) can be approximated by means of the solution at the SL-GL grid points, as aj (t) =
N 1 L L Pj (xL N,i )N,i u(xN,i , t), L j i=0
(3.10)
N 1 L L bj (t) = L Pj (xL N,i )N,i v(xN,i , t). j i=0
In virtue of Eqs. (3.8)-(3.10), we can write the approximate solutions as u(x, t) =
N N 1 L L PL,j (xL N,i )PL,j (x)N,i u(xN,i , t), L i=0 j=0 j
N N 1 L L v(x, t) = PL,j (xL N,i )PL,j (x)N,i v(xN,i , t). L j i=0 j=0
(3.11)
Furthermore, if we differentiate (3.11) once, the approximation of the first spatial partial derivatives of the approximate solutions, by means of the approximate solutions at the SL-GL interpolation nodes, can be expanded as ∂x u(xL N,n , t) =
N
κn,i u(xL N,i , t),
i=0
∂x v(xL N,n , t)
=
N
(3.12) κn,i v(xL N,i , t),
n = 0, 1, · · · , N,
i=0
where κn,i =
N L N,i j=0
L j
PL,j (xL N,i )
∂P
L,j (x)
∂x
x=xL N,n
.
(3.13)
Similarly, the second spatial partial derivatives for the approximate solution may also be written in the forms ∂x2 u(xL N,n , t) =
N
λn,i u(xL N,i , t),
i=0
∂x2 v(xL N,n , t)
=
N
(3.14) λn,i v(xL N,i , t),
i=0
8
n = 0, 1, · · · , N,
where N L N,i
λn,i =
j=0
L j
PL,j (xL N,i )
∂2P
L,j (x) ∂x2
x=xL N,n
.
(3.15)
Let us denote uk (t) = u(xL N,k , t),
vk (t) = v(xL N,k , t),
gk (t) = g(xL N,k , t),
fk (t) = f (xL N,k , t).
In the proposed spectral SL-GL-C scheme, the residual of (3.1) is set to be zero at the (N − 1) SL-GL interpolation nodes. Therefore, adopting (3.5)-(3.15), we may write (3.1)-(3.2) in the form Dtα un (t) = −
N
λn,i vi (t) − γ(u2n (t) + vn2 (t))vn (t) + δgn (t),
i=0
Dtα vn (t)
=
N
(3.16)
λn,i ui (t) +
γ(u2n (t)
+
vn2 (t))un (t)
− δfn (t),
n = 1, 2, · · · , N − 1,
i=0
In addition, no need for additional equations for treating the boundary conditions, because L the boundary conditions (3.2) are imposed at the two collocation points xL N,0 and xN,N in the
expansion of functions. This provides SFDEs with (2N − 2) unknown functions in time variable Dtα un (t) = −
N −1
λn,i vi (t) + λn,0 g3 (t) + λn,N g4 (t) − γ(u2n (t) + vn2 (t))vn (t) + δgn (t),
i=1
Dtα vn (t) =
N
λn,i ui (t) + λn,0 g1 (t) + λn,N g2 (t) + γ(u2n (t) + vn2 (t))un (t) − δfn (t),
(3.17)
i=0
n = 1, 2, · · · , N − 1, subject to the initial values un (0) = f1 (xL N,n ),
vn (0) = f2 (xL N,n ),
n = 1, · · · , N − 1.
(3.18)
In Section 5, we discuss the numerical approximation of such system of fractional differential equations.
3.2
TFSE with non-local conditions
Another kind of boundary condition will be discussed in this subsection. Let us consider the following TFSE: i
∂ α ψ(x, t) ∂ 2 ψ(x, t) + + γ|ψ(x, t)|2 ψ(x, t) − δR(x, t) = 0, ∂tα ∂x2
0 < α < 1,
(x, t) ∈ [0, L] × [0, T ], (3.19)
with the initial, boundary and non-local conditions l2 ψ(x, t)dx = ζ2 (t), ψ(0, t) =ζ1 (t), l1
ψ(x, 0) =ξ1 (x),
0 ≤ l1 < l2 ≤ L,
t ∈ [0, T ], (3.20)
x ∈ [0, L]. 9
As in subsection 3.1, Eqs. (3.19) and (3.20) can be rewritten as ∂ α u(x, t) ∂ 2 v(x, t) + + γ(u2 (x, t) + v 2 (x, t))v(x, t) − δg(x, t) = 0, ∂tα ∂x2 ∂ α v(x, t) ∂ 2 u(x, t) − − γ(u2 (x, t) + v 2 (x, t))u(x, t) + δf (x, t) = 0, ∂tα ∂x2
(3.21)
with the initial, boundary and non-local conditions
l2
l1 l2
u(0, t) = g1 (t), v(0, t) = g3 (t),
l1
u(x, 0) = f1 (x),
u(x, t)dx = g2 (t),
0 ≤ l1 < l2 ≤ L,
t ∈ [0, T ],
v(x, t)dx = g4 (t),
0 ≤ l1 < l2 ≤ L,
t ∈ [0, T ],
v(x, 0) = f2 (x),
(3.22)
x ∈ [0, L].
The problem now is how to deal with the non-local boundary conditions in (3.22). For this purpose, let us introduce a collocation treatment for these integral conservation conditions as
N N l2
l1
i=0
j=0
1 L L P (x )P (x) L,j L,j N,i N,i ui (t)dx = g2 (t), L j
0 ≤ l1 < l2 ≤ L,
(3.23)
the above equation may be restated as N N l2 1 L L ui (t) = g2 (t), PL,j (xN,i ) PL,j (x)dx N,i L l1 i=0 j=0 j
0 ≤ a1 < a2 ≤ L,
(3.24)
or briefly N
Ii ui (t) = g2 (t),
0 ≤ l1 < l2 ≤ L,
(3.25)
i=0
where Ii =
N l2 1 L L P (x ) P (x)dx N,i . L,j L,j N,i L l j 1 j=0
In the same way, we obtain N
Ii vi (t) = g4 (t),
0 ≤ l1 < l2 ≤ L,
(3.26)
i=0
and this in turn yields uN (t) =
N −1 N −1 1 1 g2 (t) − I0 u0 (t) − g2 (t) − I0 g1 (t) − Ii ui (t) = Ii ui (t) , IN IN i=1 i=1
N −1 N −1 1 1 vN (t) = Ii vi (t) = Ii vi (t) . g4 (t) − I0 v0 (t) − g4 (t) − I0 g3 (t) − IN IN i=1 i=1
10
(3.27)
Based on the information included in this subsection and the recent one, the problem with its conditions are reduced to the following compact form
Dtα un (t) = −
N −1
N −1 Ii vi (t) λn,N g4 (t) − i=1
λn,i vi (t) + (λn,0 − I0 )g3 (t) +
IN
i=1
− γ(u2n (t) + vn2 (t))vn (t)
+ δgn (t), Dtα vn (t) =
N
N −1 Ii ui (t) λn,N g2 (t) − i=1
λn,i ui (t) + (λn,0 − I0 )g1 (t) +
IN
i=0
− δfn (t),
+ γ(u2n (t) + vn2 (t))un (t)
n = 1, 2, · · · , N − 1, (3.28)
subject to the initial values un (0) = f1 (xL N,n ),
vn (0) = f2 (xL N,n ),
n = 1, · · · , N − 1,
(3.29)
This provides SFDEs with (2N − 2) unknown functions in time variable.
3.3
TFSE with non-local and mixed boundary conditions
This subsection is devoted to develop an algorithm based on the SL-GL-C method to numerically solve the TFSE subject to non-local and mixed boundary conditions; namely i
∂ α ψ(x, t) ∂ 2 ψ(x, t) + + γ|ψ(x, t)|2 ψ(x, t) − δR(x, t) = 0, ∂tα ∂x2
0 < α < 1,
(x, t) ∈ [0, L] × [0, T ], (3.30)
with the initial non-local and mixed boundary ψ(0, t) + ∂x ψ(0, t) = ζ1 (t), l2 ψ(x, t)dx = ζ2 (t), 0 ≤ l1 < l2 ≤ L, l1
ψ(x, 0) = ξ1 (x),
t ∈ [0, T ],
(3.31)
x ∈ [0, L].
Similar steps to that given in subsection (3.1), enable one to write the problem in the form ∂ α u(x, t) ∂ 2 v(x, t) + + γ(u2 (x, t) + v 2 (x, t))v(x, t) − δg(x, t) = 0, ∂tα ∂x2 ∂ α v(x, t) ∂ 2 u(x, t) − − γ(u2 (x, t) + v 2 (x, t))u(x, t) + δf (x, t) = 0, ∂tα ∂x2 with the initial, non-local and mixed boundary conditions l2 u(0, t) + ∂x u(0, t) = g1 (t), u(x, t)dx = g2 (t), 0 ≤ l1 < l2 ≤ L, v(0, t) + ∂x v(0, t) = g3 (t), u(x, 0) = f1 (x),
l1 l2
l1
v(x, t)dx = g4 (t),
v(x, 0) = f2 (x),
x ∈ [0, L]. 11
0 ≤ l1 < l2 ≤ L,
(3.32)
t ∈ [0, T ], t ∈ [0, T ],
(3.33)
Again, the problem now is how to deal with the non-local and mixed boundary conditions (3.33). For this purpose, let us propose collocation treatment for the above mixed and non-local conservation conditions u0 (t) +
N
κ0,i ui (t) = g1 (t),
i=0
v0 (t) +
N
(3.34) κ0,i vi (t) = g3 (t),
i=0
which may be written as N −1
(1 + κ0,0 )u0 (t) + κ0,N uN (t) +
κ0,i ui (t) = g1 (t),
i=1
(1 + κ0,0 )v0 (t) + κ0,N vN (t) +
N −1
(3.35) κ0,i vi (t) = g3 (t),
i=1
and as a consequence of the collocation treatment for the non-local conditions, we obtain N
N
Ii ui (t) = g2 (t),
i=0
Ii vi (t) = g4 (t).
(3.36)
i=0
According to (3.35) and (3.36), the values of u0 , uN , v0 and vN can be expressed in terms of ui (t) and vi (t), i = 1, · · · , N − 1 as
n−1 n−1 Ir ur (t) − IN κ0,i ur (t) + IN g1 (t) − κ0,N g2 (t) κ0,N r=1 r=1 , u0 (t) = IN (1 + κ0,0 ) − I0 κ0,N
n−1 n−1 κ0,i ur (t) − (1 + κ0,0 ) Ir ur (t) + (1 + κ0,0 )g2 (t) − I0 g1 (t) I0 r=1 r=1 , uN (t) = IN (1 + κ0,0 ) − I0 κ0,N
n−1 n−1 Ir vr (t) − IN κ0,i vr (t) + IN g3 (t) − κ0,N g4 (t) κ0,N r=1 r=1 , v0 (t) = IN (1 + κ0,0 ) − I0 κ0,N
n−1 n−1 κ0,i vr (t) − (1 + κ0,0 ) Ir vr (t) + (1 + κ0,0 )g4 (t) − I0 g3 (t) I0 r=1 r=1 . vN (t) = IN (1 + κ0,0 ) − I0 κ0,N
(3.37)
(3.38)
Based on the information included in this subsection and subsection 3.1, we obtain the following compact form of SFDEs Dtα un (t) = −
N −1
λn,i vi (t) + λn,0 v0 (t) + λn,N vN (t) − γ(u2n (t) + vn2 (t))vn (t) + δgn (t),
i=1
Dtα vn (t) =
N
λn,i ui (t) + λn,0 u0 (t) + λn,N uN (t) + γ(u2n (t) + vn2 (t))un (t) − δfn (t),
(3.39)
i=0
n = 1, 2, · · · , N − 1, subject to the initial values un (0) = f1 (xL N,n ),
vn (0) = f2 (xL N,n ), 12
n = 1, · · · , N − 1.
(3.40)
Remark 3.1. In the proposed method, unlike the most existing methods, the boundary, non-local and boundary, and mixed boundary and non-local conservation conditions are satisfied automatically in the collocation scheme. More specifically, no need to additional equations to enforce these conditions. Thereby, the above treatment is more reasonable. In addition, it provides high accurate numerical solution for the undertake problem.
4
Two-dimensional TFSE
In the present section, we provide an efficient algorithm for numerical treatment of the following two dimensional TFSEs in the following form i
∂ α ψ(x, y, t) ∂ 2 ψ(x, y, t) ∂ 2 ψ(x, y, t) = a1 +a2 + |ψ(x, y, t)|2 ψ(x, y, t) + R(x, y, t), α 2 ∂t ∂x ∂y 2 0 < α < 1,
(4.1)
(x, y, t) ∈ Ω1 × Ω2 × Ω3 ,
subject to the initial condition ψ(x, y, 0) = ζ1 (x, y),
(x, y) ∈ Ω1 × Ω2 ,
(4.2)
and the four boundary conditions ψ(0, y, t) = ζ2 (y, t),
ψ(L1 , y, t) = ζ3 (y, t),
(y, t) ∈ Ω2 × Ω3 ,
ψ(x, 0, t) = ζ4 (x, t),
ψ(x, L2 , t) = ζ5 (x, t)
(x, t) ∈ Ω1 × Ω3 ,
where Ω1 = [0, L1 ],
Ω2 = [0, L2 ]
and
(4.3)
Ω3 = [0, T ], while ζ1 (x, y), ζ2 (y, t), ζ3 (y, t), ζ4 (x, t),
ζ5 (x, t) and R(x, y, t) are given functions. Firstly, the mapping of complex to real and imaginary parts is required to the possibility of using the numerical method. To this end, suppose the transformations ψ(x, y, t) = u(x, y, t) + i v(x, y, t), ζ1 (x, y) = g1 (x, y) + i g2 (x, y), ζ3 (y, t) = g5 (y, t) + i g6 (y, t),
R(x, y, t) = f (x, y, t) + i g(x, y, t), ζ2 (y, t) = g3 (y, t) + i g4 (y, t),
(4.4)
ζ4 (x, t) = g7 (x, t) + i g8 (x, t),
ζ5 (x, t) = g9 (x, t) + i g10 (x, t), where u(x, y, t), v(x, y, t), f (x, y, t), g(x, y, t), g1 (x, y), g2 (x, y), g3 (y, t), g4 (y, t), g5 (y, t), g6 (y, t), g7 (x, t), g8 (x, t), g9 (x, t) and g10 (x, t), are real functions. Using (4.4), with the help of (4.1), we obtain a coupled system of two-dimensional fractional PFDs: ∂ 2 v(x, y, t) ∂ 2 v(x, y, t) ∂ α u(x, y, t) = a1 +a2 + ((u(x, y, t))2 + (v(x, y, t))2 ) v(x, y, t) + g(x, y, t), α 2 ∂t ∂x ∂y 2 ∂ α v(x, y, t) ∂ 2 u(x, y, t) ∂ 2 u(x, y, t) = a +a + ((u(x, y, t))2 + (v(x, y, t))2 ) u(x, y, t) + f (x, y, t), − 1 2 ∂tα ∂x2 ∂y 2 (x, y, t) ∈ Ω1 × Ω2 × Ω3 , (4.5) 13
subject to the initial conditions
u(x, y, 0) = g1 (x, y),
v(x, y, 0) = g2 (x, y),
(x, y) ∈ Ω1 × Ω2 ,
(4.6)
and the boundary conditions u(0, y, t) = g3 (y, t),
u(L1 , y, t) = g4 (y, t),
v(0, y, t) = g5 (y, t),
v(L1 , y, t) = g6 (y, t),
u(x, 0, t) = g7 (x, t),
u(x, L2 , t) = g8 (x, t),
v(x, 0, t) = g9 (x, t),
v(x, L2 , t) = g10 (x, t),
(y, t) ∈ Ω2 × Ω3 ,
(4.7)
(x, t) ∈ Ω1 × Ω3 .
Secondly, we are interested in using the SL-GL-C method to reduce the above coupled system of two-dimensional fractional PFDs with their boundary conditions into SFDEs which greatly simplify the problem. In order to do this, we outline the main steps of our algorithm based on approximating the problem using SL-GL-C method for spatial discretization. Let us expand the approximate solution in a doubly shifted Legendre series u(x, y, t) =
M N
ai,j (t)PL1 ,i (x)PL2 ,j (y),
i=0 j=0
v(x, y, t) =
M N
(4.8) bi,j (t)PL1 ,i (x)PL2 ,j (y),
i=0 j=0 L2 1 Let {xL N,i ; 0 ≤ i ≤ N } and {yM,j ; 0 ≤ j ≤ M } be the SL-GL interpolation nodes for the shifted
Legendre polynomials PL1 ,N (x) and PL2 ,M (y) respectively. Making use of SL-GL quadrature and due to the orthogonality relation, one obtains ai,j (t) =
N M 1 L2 L2 L1 L1 L1 L2 P (y ) P (x ) L ,j L ,i 2 1 M,k M,k N,l N,l u(xN,l , yM,k , t), 1 L2 L i j l=0 k=0
N M 1 L2 L2 L1 L1 L2 1 )M,k PL1 ,i (xL bi,j (t) = L1 L2 PL2 ,j (yM,k N,l )N,l v(xN,l , yM,k , t). i j l=0 k=0
(4.9)
Let us denote L2 1 u(xL N,n , yM,m , t) = un,m (t),
L2 1 v(xL N,n , yM,m , t) = vn,m (t),
then the approximate solutions (4.8) may be expressed in the form L2 L2 L1 1 M N M N PL2 ,j (yM,k )M,k PL1 ,i (xL ) N,l N,l PL1 ,i (x)PL2 ,j (y) ul,k (t), u(x, y, t) = L1 L2 i j i=0 j=0 l=0 k=0 L2 L2 L1 1 M N M N PL2 ,j (yM,k )M,k PL1 ,i (xL N,l )N,l v(x, y, t) = PL1 ,i (x)PL2 ,j (y) vl,k (t), 1 L2 L i j i=0 j=0 l=0 k=0
(4.10)
In what follows, the first-order spatial partial derivative with respect to x for the solutions L2 1 (4.10), at a specific collocation nodes xL N,n and yM,m can be written as
14
∂un,m (t) = ∂x
L2 L2 L1 1 L1 M N N M )M,k PL1 ,i (xL ) PL2 ,j (yM,k N,l N,l ∂(PL,i (xN,n )) 1 L2 L i j
l=0 k=0 i=0 j=0
n = 0, 1, · · · , N,
∂x
L2 PL2 ,j (yM,m ) ul,k (t),
m = 0, 1, · · · , M, (4.11)
for simplifying the notation, let M
N
∂un,m (t) n,m = γi,j ul,k (t), ∂x
(4.12)
l=0 k=0
n,m are the expansion coefficients of the derivative, and are obtained from where γi,j L2 L2 L1 1 L1 M N PL2 ,j (yM,k )M,k PL1 ,i (xL ) N,l N,l ∂(PL,i (xN,n )) n,m L2 PL2 ,j (yM,m = ). γi,j L1 L2 ∂x i j i=0 j=0
(4.13)
L2 1 Accordingly, the first-order derivative with respect to y for u(x, y, t), at the nodes xL N,n and yM,m ,
is
L2 L2 L1 L1 L2 N M N M ∂(PL2 ,j (yM,m )) ∂un,m (t) PL2 ,j (yM,k )M,k PL1 ,i (xN,l )N,l L1 = ul,k (t), P (x ) L,i N,n L L 1 2 ∂y ∂y i j i=0 j=0 l=0 k=0
n = 0, 1, · · · , N,
m = 0, 1, · · · , M. (4.14)
Similar argument leads to N
M
∂un,m (t) n,m = δi,j ul,k (t), ∂y
(4.15)
l=0 k=0
where n,m = δi,j
L2 L2 L1 1 M N PL2 ,j (yM,k )M,k PL1 ,i (xL ) N,l N,l i=0 j=0
1 L2 L i j
1 PL,i (xL N,n )
L2 ∂(PL2 ,j (yM,m ))
∂y
.
(4.16)
Also, the second-order spatial partial derivatives with respect to x and y are respectively given by
N
M
∂ 2 un,m (t) n,m = ξi,j ul,k (t), ∂x2 l=0 k=0 N
(4.17)
M
∂ 2 un,m (t) n,m = ηi,j ul,k (t), ∂y 2 l=0 k=0
where n,m = ξi,j
and n,m ηi,j =
L2 L2 L1 1 L1 N M 2 PL2 ,j (yM,k )M,k PL1 ,i (xL ) N,l N,l ∂ PL,i (xN,n ) i=0 j=0
1 L2 L i j
∂x2
L2 L2 L1 1 M N )M,k PL1 ,i (xL PL2 ,j (yM,k N,l )N,l i=0 j=0
1 L2 L i j
15
1 PL,i (xL N,n )
L2 PL2 ,j (yM,m ),
L2 ∂ 2 PL2 ,j (yM,m )
∂y 2
.
(4.18)
(4.19)
In a similar way, the second approximate solution for the coupled system can be obtained from N
M
∂vn,m (t) n,m = γi,j vl,k (t), ∂x l=0 k=0
2
∂ vn,m (t) = ∂x2
M N
n,m ξi,j
M
N
∂vn,m (t) n,m = δi,j vl,k (t), ∂y l=0 k=0 N
M
∂ vn,m (t) n,m = ηi,j vl,k (t). ∂y 2 2
vl,k (t),
l=0 k=0
(4.20)
l=0 k=0
In the shifted Legendre pseudo-spectral approximation for the two-dimensional version of TFSEs, the residual of (4.5) is set to be zero at (N − 1) × (M − 1) of the nodes of SL-GL interpolation. Therefore, adopting (4.10)-(4.20), enable one to write (4.5)-(4.7) as Dα un,m (t) =a1
M N
n,m ξi,j vl,k (t) + a2
l=0 k=0
M N
n,m ηi,j vl,k (t) + ((un,m (t))2 + (vn,m (t))2 )vn,m (t)
l=0 k=0
+ gn,m (t), Dα vn,m (t) = − a1
M N
n,m ξi,j ul,k (t) − a2
l=0 k=0
− fn,m (t),
M N
n,m ηi,j ul,k (t) − ((un,m (t))2 + (vn,m (t))2 )un,m (t)
l=0 k=0
n = 1, · · · , N − 1,
m = 1, · · · , M − 1, (4.21)
subject to the initial conditions L2 1 un,m (0) = g1 (xL N,n , yM,m ),
n = 1, · · · , N − 1,
L2 1 vn,m (0) = g2 (xL N,n , yM,m ),
(4.22)
m = 1, · · · , M − 1.
Moreover, the values of u0,k (t), uN,k (t), ul,0 (t) ul,N (t) v0,k (t), vN,k (t), vl,0 (t) and vl,N (t) are obtained from the relations L2 u0,k (t) = g3 (yM,k , t), 1 ul,0 (t) = g7 (xL N,l , t),
L2 v0,k (t) = g5 (yM,k , t), 1 vl,0 (t) = g9 (xL N,l , t),
L2 uN,k (t) = g4 (yM,k , t), 1 ul,N (t) = g8 (xL N,l , t),
L2 vN,k (t) = g6 (yM,k , t), 1 vl,N (t) = g10 (xL N,l , t),
k = 0, · · · , M, l = 0, · · · , N
(4.23)
k = 0, · · · , M, l = 0, · · · , N.
This provides a SFDEs in time variable with 2(N − 1) × (M − 1) unknowns. The spectral solution of this system will be discussed in details in the next section.
5
System of fractional differential equations
In this section, we propose an efficient numerical integration process for SFDEs with vector of initial values, based on SL-GR interpolation, which is easy to be implemented and possesses the spectral accuracy. This scheme has three fascinating advantages. It is easier to be implemented for nonlinear problems, it is also specially appropriate for long-time calculations and is more stable for
16
large N . This is also confirmed by the numerical results. Now, we consider the following nonlinear system of fractional initial value problems which is a generalization of the systems given in (3.17), (3.28), (3.39) and (4.21), namely Dα ur (t) = Gr (t, u1 (t), · · · , uR (t)) ,
r = 1, · · · , R,
0 < α < 1,
t ∈ [0, T ],
(5.1)
subject to r = 1, · · · , R,
ur (0) = τr , where Gr (t, u1 (t), · · · , uR (t)),
(5.2)
r = 1, · · · , R, are given nonlinear functions.
We are interested in using the SL-GR-C method to transform the previous SFDEs into a system of algebraic equations. In order to do this, we approximate the time variable using SL-GR-C method at tTN,i SL-GR interpolation nodes. Let the approximate solutions be of the form ur (t) =
K
r = 1, · · · , R.
ar,j PT,j (t),
(5.3)
j=0
Furthermore, the approximation of the time fractional derivative Dα can be computed as Dα ur (t) =
K
ar,j Dα PT,j (t),
r = 1, · · · , R.
(5.4)
j=0
The fractional derivative of shifted Legendre polynomials is given by (α)
Dα PT,i (t) = PT,i (t) =
∞
Πα (i, j)PT,j (t),
i = α, α + 1, · · · ,
(5.5)
j=0
where i (−1)i+ (2j + 1) (i + )! ( − j − α + 1)j . T α (i − )! ! Γ( − α + 1) ( − α + 1)j+1
Πα (i, j) =
=α
Consequently Dα ur (t) =
K j=0
(α)
ar,j PT,j (t),
r = 1, · · · , R.
Therefore, adopting (5.3)-(5.6), enable one to write (5.1)-(5.2) in the form: ⎛ ⎞ K K K (α) ar,j PT,j (t) = Gr ⎝t, a1,j PT,j (t), · · · , aR,j PT,j (t)⎠ , r = 1, · · · , R, j=0
j=0
(5.6)
t ∈ [0, T ], (5.7)
j=0
K
ar,j PT,j (0) = τr ,
r = 1, · · · , R.
(5.8)
j=0
In the proposed method the residual of (5.1) is set to be zero at RK collocation points. Moreover, the initial conditions in (5.2) will be collocated at the SL-GR interpolation nodes. Firstly, we obtain RK algebraic equations K j=0
(α) ar,j PT,j (tTK,s )
⎛
=Gr ⎝tTK,s ,
K
a1,j PT,j (tTK,s ), · · ·
j=0
r = 1, · · · , R,
,
K j=0
s = 1, · · · , K. 17
⎞ aK,j PT,j (tTK,s )⎠ , (5.9)
Secondly, due to initial conditions, we have other R algebraic equations K
ar,j PT,j (tTK,0 ) = τr ,
r = 1, · · · , R.
(5.10)
j=0
Finally, Eq. (5.9) guaranteed that the system (5.1) is satisfied exactly at the SL-GR interpolation nodes tTK,s ; s = 1, · · · , R. This provides RK nonlinear algebraic equations for ar,j ; r = 1, · · · , R, j = 0, · · · , K + 1. In addition, the collocation treatment of the initial values in (5.10) provides R linear algebraic equations in the unknown shifted Legendre expansion coefficients ar,j ; r = 1, · · · , R, j = 0, · · · , K + 1. The combination of these two algebraic systems constitutes a system of R(K + 1) algebraic equations which we solved it by Newton’s iteration method. Consequently, the approximate solutions (5.3) can be evaluated.
6
Numerical simulation and comparisons
This section reports several numerical examples to demonstrate the high accuracy and applicability of the proposed methods for solving one- and two-dimensional time-fractional Schr¨odinger equations. We also compare the results given from our scheme and those reported in the literature such as local discontinuous Galerkin (LDG) method [40] and the meshless technique based on collocation framework and radial basis functions [41]. The comparisons reveal that our methods are very effective and convenient.
6.1
Linear TFSE
First, we consider the time-fractional linear Schr¨odinger equation studied in [40]:
−α −α ∂αψ ∂2ψ 2t cos(2πx) 2t sin(2πx) 2 2 i α + + sin(x) + t − cos(x) , = − it ∂t ∂x2 Γ(3 − α) Γ(3 − α)
(6.1)
(x, t) ∈ [0, 2π] × [0, T ], subject to the initial-boundary conditions which can be obtained from the exact solution: ψ(x, t) = t2 (cos(x) + i sin(x)).
(6.2)
Recently, Wei et al. [40] solved this problem using the implicit fully discrete local discontinuous Galerkin (LDG) method. More recently, Bhrawy et al. [48] proposed the operational matrix formulation of Jacobi spectral collocation (JSC) method for solving such problem. In order to show that our method is accurate than the LDG [40] and JSC [48] methods. In Table 6.1, the numerical results based on maximum absolute errors obtained using the proposed algorithm, for real (M1 ) and imaginary (M2 ) parts of the solution at L = T = 1, are compared with 18
Table 6.1: Maximum absolute errors of problem (6.1) for T = 1. Our method
LDG [40] method at Δt = 1/1000
α
N =M
M1
M2
M1
M2
0.1
5
1.25 × 10−1
8.12 × 10−2
3.19 × 10−2
3.12 × 10−2
10
1.02 × 10−5
1.59 × 10−5
4.11 × 10−3
3.89 × 10−3
15
2.06 × 10−9
8.17 × 10−10
1.23 × 10−3
1.22 × 10−3
20
7.66 × 10−15
1.50 × 10−14
5.19 × 10−4
5.16 × 10−4
5
1.26 × 10−1
8.03 × 10−2
3.21 × 10−2
3.13 × 10−2
10
2.72 × 10−6
6.09 × 10−5
4.12 × 10−3
3.90 × 10−3
15
2.06 × 10−9
7.95 × 10−10
1.2 × 10−3
1.23 × 10−3
20
4.77 × 10−15
1.31 × 10−14
5.24 × 10−4
5.21 × 10−4
0.5
those presented by LDG [40] method at various N (the number of nodes of the local discontinuous Galerkin method for spatial discretization) and Δt = 1/1000 (the time step of finite difference method for temporal discretization). It is observed that the proposed method is more accurate than LDG [40] method. Moreover, In Table 6.2, we list the M1 and M2 errors. In addition, the L2 errors for real (L2,Re ) and imaginary (L2,Im ) parts, are also presented in Table 6.2, for different two choices of α. Comparing the numerical results of the present method and those achieved by JSC [48] method is displayed in Table 6.2. From this table, we observe that the present method is accurate than the proposed spectral method in [48]. The real and imaginary parts of the approximate solution obtained by the present method for α = 0.3 at N = 20 are shown in Fig. 1, for the space variable x ∈ [0, 6.5] and time value t ∈ [0, 1]. Meanwhile, absolute errors for real and imaginary solutions for the fractional derivative α = 0.3, in case of x ∈ [0, 6.5] and t ∈ [0, 1] at N = 20, are plotted in Fig. 2.
19
Table 6.2: Comparing the L2 , M1 and M2 errors of the present method and JSC [48] method for problem (6.1). Our method
0.1
0.5
JSC [48] method
N =M
L2,Re
L2,Im
M1
M2
L2,Im
M1
M2
8
1.92 × 10−5
1.84 × 10−5
1.94 × 10−4
3.34 × 10−4
6.14 × 10−5
2.07 × 10−4
8.82 × 10−4
12
3.41 × 10−7
1.29 × 10−7
1.27 × 10−6
4.96 × 10−7
1.60 × 10−6
6.45 × 10−6
2.38 × 10−6
16
5.81 × 10−11
2.19 × 10−11
2.10 × 10−10
7.87 × 10−11
...
...
...
−15
−15
−14
−15
...
...
...
20
3.50 × 10
8
4.72 × 10−4
12
1.32 × 10
1.27 × 10
5.11 × 10
1.67 × 10−4
1.94 × 10−4
8.74 × 10−4
2.16 × 10−4
2.05 × 10−4
1.03 × 10−3
3.25 × 10−7
1.179 × 10−7
1.27 × 10−6
5.15 × 10−7
3.59 × 10−5
3.65 × 10−6
5.14 × 10−5
16
5.54 × 10−11
2.02 × 10−11
2.10 × 10−10
7.66 × 10−11
...
...
...
20
3.37 × 10−15
1.25 × 10−15
1.25 × 10−14
5.11 × 10−15
...
...
...
(a) Real part
(b) Imaginary part
Figure 1: The real and imaginary parts of the approximate solution for x ∈ [0, 6.5] and t ∈ [0, 1] at N = 20 of problem (6.1).
20
(a) Real part
(b) Imaginary part
Figure 2: Absolute errors of real and imaginary parts for x ∈ [0, 6.5] and t ∈ [0, 1] at N = 20 of problem (6.1).
6.2
Attractive nonlinear TFSE with boundary conditions
We consider the time-fractional attractive nonlinear Schr¨odinger equation studied in [40, 41]
4 ∂αψ ∂2ψ 2t−α cos(2πx) 2 2 2 i α + + |ψ| ψ =t − 4π t sin(2πx) − ∂t ∂x2 Γ(3 − α)
−α sin(2πx) 2t 2 4 2 + t − 4π cos(2πx) , (x, t) ∈ [0, L] × [0, T ]. + it Γ(3 − α) (6.3) The exact solution of this problem is ψ(x, t) = t2 (sin(2πx) + i cos(2πx)).
(6.4)
The initial-boundary conditions of this problem can be extracted from the exact solution. In Table 6.3, we compare our numerical results obtained by the present method with those obtained by meshless technique [41] for two choices of the fractional derivative α = 0.1 and α = 0.3. From this table, we see that we can achieve an excellent approximation for the exact solution by using proposed method for a limited number of the collocation nodes. It is clear that our method is more accurate than the meshless technique [41]. The complete agreement of the approximate and exact solutions for the real and imaginary parts are clearly appeared in Fig. 3, for different choices of the variable t = 0.3, 0.5, 0.7 at T = 6 . In addition, the curves of real (E1 (x, t)) and imaginary (E2 (x, t)) parts of the absolute error of problem (6.3) at T = 1 are displayed in Fig. 4. In addition, to confirm the high accuracy and convergence of the present scheme, in Fig. 5, we plot the logarithmic graphs of the maximum absolute
21
Table 6.3: Maximum absolute errors of problem (6.3) for T = 1. α = 0.1
α = 0.3
N, M
M1
M2
M3
M1
M2
M3
16,16
5.73 × 10−11
2.70 × 10−10
2.74 × 10−10
6.69 × 10−11
2.66 × 10−10
2.73 × 10−10
20,20
4.39 × 10−15
2.21 × 10−14
2.21 × 10−14
6.66 × 10−15
2.38 × 10−14
6.66 × 10−15
24,24
1.01 × 10−14
5.00 × 10−15
1.12 × 10−14
4.11 × 10−15
5.77 × 10−15
6.17 × 10−15
2.18 × 10−3
3.5684 × 10−3
meshless technique [41] 2.85 × 10−3
2.18 × 10−3
3.57 × 10−3
2.86 × 10−3
0.4
0.4
0.2
0.2
0.0
v and vN
u and uN
29, 250
ux,0.3
0.0
vx,0.3 vN x,0.3
uN x,0.3 0.2
0.2
ux,0.5
vx,0.5 vN x,0.5
uN x,0.5 ux,0.7
0.4 0
1
vx,0.7
0.4
uN x,0.7 2
3
4
5
6
vN x,0.7 0
x
1
2
3
4
5
6
x
(a) Real part
(b) Imaginary part
Figure 3: Curves of real and imaginary parts of the approximate and exact solutions of problem (6.3). errors (log10 M1 and log10 M2 ) for α = 0.8 at various values of N (N = M = 4, 6, · · · , 26). This demonstrates that the proposed method provides an accurate approximation and yields exponential convergence rates.
22
1. 1013
1. 1013 8. 1014
6. 1014
E2 x,1.0
E1 x,1.0
8. 1014
4. 1014 2. 1014
6. 1014 4. 1014 2. 1014
0
0 0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
x
0.8
1.0
x
(a) Real part
(b) Imaginary part
Figure 4: Absolute errors of real and imaginary solutions of problem (6.3).
0
2
8
Log10 M 2
Log10 M 1
6
10
12
6
8
10
14 15
12
10
4
5
2
4
20
14
25
5
N
10
15
20
N
(a) Real part
(b) Imaginary part
Figure 5: Logarithmic graphs of the maximum absolute errors of problem (6.3).
23
25
Table 6.4: Maximum absolute errors of problem (6.5) at α = 0.5.
6.3
N =M
M1
M2
M3
4
7.78 × 10−2
1.75 × 10−2
7.78 × 10−2
8
2.92 × 10−5
1.69 × 10−5
3.34 × 10−5
12
1.22 × 10−9
6.84 × 10−10
1.40 × 10−9
16
1.16 × 10−13
5.73 × 10−14
1.17 × 10−13
20
3.40 × 10−13
2.68 × 10−13
3.90 × 10−13
24
1.41 × 10−13
5.93 × 10−14
1.52 × 10−13
Attractive nonlinear TFSE with boundary and non-local conditions
In order to show that our technique can be applied to TFSE with non-local conditions, we consider the time-fractional attractive nonlinear Schr¨odinger equation:
i where f (x, t) =
∂αψ ∂2ψ + + |ψ|2 ψ =f (x, t), ∂tα ∂x2
(x, t) ∈ [0, 1] × [0, T ],
et cos(πx)Γ(1 − α, t) − Γ(1 − α) sin(πx) −e2t + cos(x) + π 2 + cos(πx) Γ(1 − α)
sin(πx)Γ(1 − α, t) + cos(πx) − −e2t + cos(x) + π 2 + sin(πx) , + iet − Γ(1 − α)
(6.5)
(6.6)
with the initial, boundary and non-local conditions ψ(0, t) = i et ,
1
ψ(x, t)dx = 0
2et , π
ψ(x, 0) = (sin(πx) + i cos(πx)).
(6.7)
The exact solution is given by ψ(x, t) = et (sin(πx) + i cos(πx)).
(6.8)
Table 6.4 lists maximum absolute errors of u(x, t), v(x, t) and ψ(x, t) of problem (6.5) at α = 0.5 with various choices of N and M . In case of α = 0.5, we display the curves of real and imaginary parts of the approximate and exact solutions in Fig. 6 at N = M = 10. In Fig. 7, we display the curves of absolute errors for the real and imaginary parts of the approximate solution of problem (6.5) for α = 0.9 and N = M = 24.
6.4
1+1 TFSE with non-local and mixed boundary conditions
Now, we consider the one-dimensional TFSE:
i
∂αψ ∂2ψ + + |ψ|2 ψ =f (x, t), ∂tα ∂x2 24
(x, t) ∈ [0, 1] × [0, T ],
(6.9)
2.0
vx,0.3 vN x,0.3 vx,0.5 vN x,0.5 vx,0.7 vN x,0.7
1 v and vN
1.5 u and uN
2
ux,0.3 uN x,0.3 ux,0.5 uN x,0.5 ux,0.7 uN x,0.7
1.0
0
1
0.5
2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
x
0.6
0.8
1.0
x
(a) Real part
(b) Imaginary part
Figure 6: Curves of real and imaginary parts of the approximate and exact solutions of problem (6.5).
1. 1013 4. 1014 8. 10
14
E2 x,0.5
E1 x,0.5
3. 1014 6. 1014 4. 1014
2. 1014
1. 1014
2. 1014
0
0 0.0
0.2
0.4
0.6
0.8
1.0
0.0
x
0.2
0.4
0.6
0.8
1.0
x
(a) Real part
(b) Imaginary part
Figure 7: Absolute errors for real and imaginary parts of the approximate solution of problem (6.5).
25
Table 6.5: Maximum absolute errors of problem (6.9) for T = 1 and α = 0.5. N =M
M1
M2
M3
8
5.40 × 10−3
3.73 × 10−3
6.56 × 10−3
12
2.53 × 10−7
2.02 × 10−7
3.24 × 10−7
16
2.52 × 10−11
1.84 × 10−11
3.09 × 10−11
(a) Real part
(b) Imaginary part
Figure 8: The real and imaginary parts of the approximate solution of problem (6.9). where f (x, t) =
et cos(πx)Γ(1 − α, t) − Γ(1 − α) sin(πx) −e2t + cos(x) + π 2 + cos(πx) Γ(1 − α)
sin(πx)Γ(1 − α, t) + cos(πx) − −e2t + cos(x) + π 2 + sin(πx) , + iet − Γ(1 − α)
(6.10)
with non-local and mixed boundary conditions ψ(0, t) +
∂ψ(x, t)
= (π + i)et ,
∂x x=0
1
ψ(x, t)dx = 0
2et , π
ψ(x, 0) = (sin(πx) + i cos(πx)). (6.11)
The exact solution is ψ(x, t) = et (sin(πx) + i cos(πx)).
(6.12)
Table 6.5 lists maximum absolute errors for the real, imaginary and modulus of the approximate solutions of problem (6.9) using the proposed collocation scheme for different choices of the collocation nodes. We also plot the space graphs of the approximate solution and absolute errors for real and imaginary parts of the approximate solution in Figs. 8 and 9, respectively.
26
(a) Real part
(b) Imaginary part
Figure 9: Absolute errors of real and imaginary parts of the approximate solution of problem (6.9).
6.5
2+1 attractive TFSE
In order to confirm the high accuracy of our technique for the two-dimensional problem, we consider the 2+1 time-fractional attractive nonlinear Schr¨ odinger equation i
∂αψ ∂2ψ ∂2ψ + + + |ψ|2 ψ =f (x, y, t), ∂tα ∂x2 ∂y 2
where
(x, y, t) ∈ [0, 1] × [0, 1] × [0, 1],
4 2t−α cos(x) cos(y) 1 4 f (x, y, t) = t sin(x) sin(y) t cos(2x) cos(2y) + t − 4 − 2 Γ(3 − α) −α
2t sin(x) sin(y) 1 + i t2 + cos(x) cos(y) t4 cos(2x) cos(2y) + t4 − 4 Γ(3 − α) 2 2
(6.13)
(6.14)
with the initial-boundary conditions ψ(0, y, t) = it2 cos(y),
ψ(1, y, t) = it2 (sin(1) sin(y) + i cos(y)),
(y, t) ∈ [0, 1] × [0, 1],
ψ(x, 0, t) = it2 cos(x),
ψ(x, 1, t) = it2 (sin(1) sin(x) + i cos(x)),
(x, t) ∈ [0, 1] × [0, 1],
ψ(x, y, 0) = 0,
(6.15)
(x, y) ∈ [0, 1] × [0, 1].
The exact solution of this problem is ψ(x, y, t) = t2 (sin(x) sin(y) + i cos(x) cos(y)).
(6.16)
Table 6.6 lists the maximum absolute errors of u(x, y, t), v(x, y, t) and ψ(x, y, t) of problem (6.13) at α = 0.2 with various choices of N, M, and K. The numerical results presented in this table show that the results are very accurate for small value of N, M, and K. Fig. 10 demonstrates that the absolute errors E1 (x, y, t) and E2 (x, y, t) are very small for even the small number of grid points taken. Moreover, we see the agreement of curves of real and imaginary parts of the approximate and exact solutions in Fig. 11. 27
Table 6.6: Maximum absolute errors of problem (6.13). N =M =K
M1
M2
M3
4
6.15 × 10−4
1.36 × 10−3
1.36 × 10−3
6
1.43 × 10−6
3.14 × 10−6
3.14 × 10−6
8
1.76 × 10−9
3.84 × 10−9
3.84 × 10−9
10
7.91 × 10−11
8.84 × 10−11
1.05 × 10−10
3. 1011
8. 1011
6. 1011
E2 x,0.5,0.5
E1 x,0.5,1.0
2.5 1011
4. 1011
2. 1011 1.5 1011 1. 1011
2. 1011
5. 1012 0
0 0.0
0.2
0.4
0.6
0.8
0.0
1.0
0.2
0.4
0.6
0.8
1.0
x
x
(a) Real part
(b) Imaginary part
Figure 10: x-directional curves of real and imaginary parts of the absolute error of problem (6.13) with α = 0.2 and N = M = K = 10.
v0.3,0.3,t
u0.3,0.3,t 0.20
v N 0.3,0.3,t
0.8
uN 0.3,0.3,t
v0.4,0.4,t
u0.4,0.4,t uN 0.4,0.4,t u0.5,0.5,t uN 0.5,0.5,t
0.10
v N 0.4,0.4,t
0.6 v and v N
u and uN
0.15
v0.5,0.5,t v N 0.5,0.5,t 0.4
0.2
0.05
0.00
0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0
t
0.2
0.4
0.6
0.8
1.0
t
(a) Real part
(b) Imaginary part
Figure 11: t-directional curves of real and imaginary parts of the exact and numerical solutions of problem (6.13) with α = 0.2 and N = M = K = 10.
28
Table 6.7: Maximum absolute errors of problem (6.17).
6.6
N =M =K
M1
M2
M3
4
1.34 × 10−3
1.55 × 10−3
1.55 × 10−3
6
3.11 × 10−6
3.56 × 10−6
3.58 × 10−6
8
3.82 × 10−9
4.36 × 10−9
4.38 × 10−9
10
8.91 × 10−11
9.27 × 10−11
1.08 × 10−10
2+1 repulsive TFSE
Finally, we present the following 2+1 time-fractional nonlinear repulsive Schr¨odinger equation t3−α 6 − i t6 + 1 Γ(4 − α)tα (cos(x + y) − i sin(x + y)) ∂αψ 1 ∂2ψ 1 ∂2ψ 2 , i α + + − |ψ| ψ = ∂t 2 ∂x2 2 ∂y 2 Γ(4 − α) (x, y, t) ∈ [0, 1] × [0, 1] × [0, 1], (6.17) with the initial-boundary conditions ψ(0, y, t) = t3 (cos(y) + i sin(y)), ψ(1, y, t) = t3 (cos(1 + y) + i sin(1 + y)),
(y, t) ∈ [0, 1] × [0, 1],
ψ(x, 0, t) = t3 (cos(x) + i sin(x)),
(x, t) ∈ [0, 1] × [0, 1],
ψ(x, 1, t) = t3 (cos(x + 1) + i sin(x + 1)), ψ(x, y, 0) = 0,
(x, y) ∈ [0, 1] × [0, 1]. (6.18)
The exact solution of this problem is ψ(x, y, t) = t3 (cos(x + y) + i sin(x + y)).
(6.19)
Table 6.7 lists maximum absolute errors of u(x, y, t), v(x, y, t) and ψ(x, y, t) of problem (6.17) at α = 0.5 with various choices of N, M, and K. we observed that, a good approximation of the exact solution of the two-dimensional time-fractional nonlinear repulsive Schrodinger equation is achieved.
29
7
Conclusions
In this paper, we have proposed a new collocation algorithm to introduce an accurate numerical solution for the one-dimensional nonlinear TFSEs with different types of boundary conditions. The core of the proposed method was to discretize the TFSE in the spatial direction by SL-GL-C method, along with a new treatment for the subjected conditions, to create a system of SFDEs of the unknown coefficients of spectral expansion in time direction. An efficient numerical integration process for SFDEs was investigated based on the SL-GL-C method. The proposed method was extended to solve the two-dimensional TFSEs. The main advantage of the proposed algorithm is, adding few terms of the shifted SL-GL and SL-GR collocation nodes, a good approximation of the exact solution of the problem was achieved. Comparisons between our approximate solutions of the problems with their exact solutions and with the approximate solutions achieved by other methods were introduced to confirm the validity and accuracy of our scheme.
Acknowledgments This paper was funded by the Deanship of Scientific Research DSR, King Abdulaziz University, Jeddah, Grant no. 130-802-D1435. The authors, therefore, acknowledges with thanks DSR for the technical and financial support. The authors are very grateful to the reviewers for carefully reading the paper and for their comments and suggestions which have improved the paper.
References [1] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods: Fundamentals in Single Domains. Springer-Verlag, New York (2006) [2] W. Heinrichs, Spectral methods with sparse matrices. Numer. Math. 56, (1989) 25-41 [3] W. Heinrichs, Algebraic spectral multigrid methods. Comput. Methods Appl. Mech. Eng. 80, (1990) 281-286 [4] E.H. Doha, A.H. Bhrawy, D. Baleanu, R.M. Hafez, A new Jacobi rational-Gauss collocation method for numerical solution of generalized Pantograph equations, Appl. Numer. Math., 77 (2014) 43-54 [5] M. Tatari, M. Haghighi, A generalized Laguerre-Legendre spectral collocation method for solving initial-boundary value problems, Appl. Math. Modell., 38 (2014) 1351-1364
30
[6] M. Zayernouri, W. Cao, Z. Zhang, G.E. Karniadakis, Spectral and Discontinuous Spectral Element Methods for Fractional Delay Equations, SIAM J. Sci. Comput., 36 (2014), B904B929 [7] M. Zayernouri, G.E. Karniadakis, Fractional spectral collocation method, SIAM J. Sci. Comput., 36 (2014), A40-A62 [8] A.H. Bhrawy, An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system, Appl. Math. Comput., 247 (2014) 30-46 [9] S. Nemati, Numerical solution of Volterra-Fredholm integral equations using Legendre collocation method, J. Comput. Appl. Math., 278 (2015) 29-36 [10] A.H. Bhrawy, M.A. Zaky, Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation, Nonlinear Dyn., 80 (2015) 101-116 [11] M. Zayernouri, G.E. Karniadakis, Fractional Sturm-Liouville eigen-problems: Theory and numerical approximation, J. Comput. Phys., 252 (2013) 495-517 [12] M. Zayernouri, G.E. Karniadakis, Exponentially accurate spectral and spectral element methods for fractional ODEs, J. Comput. Phys., 257, (2014) 460-480 [13] D. Cafagna, Fractional calculus: a mathematical tool from the past for the present engineer, IEEE Transactions on Industrial Electronic (2007) 35-40 [14] J.W. Kirchner, X. Feng, C. Neal, Fractal stream chemistry and its implications for containant transport in catchments, Nature, 403 (2000) 524-526 [15] M. Giona, H.E. Roman, Fractional diffusion equation for transport phenomena in random media, Phys. A, 185 (1992) 87-97 [16] R.L. Magin, Fractional Calculus in Bioengineering, Begell House Publishers, (2006). [17] I. Podlubny, Fractional Differential Equations, in: Mathematics in Science and Engineering, Academic Press Inc., San Diego, CA, (1999) [18] R. Hilfer, Applications of Fractional Calculus in Physics, Word Scientific, Singapore, (2000) [19] M. Amairi, M. Aoun, S. Najar, M.N. Abdelkrim, A constant enclosure method for validating existence and uniqueness of the solution of an initial value problem for a fractional differential equation, Appl. Math. Comput., 217 (2010) 2162-2168
31
[20] J. Deng, L. Ma, Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations, Appl. Math. Lett., 23 (2010) 676-680 [21] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, San Diego, (2006) [22] M. M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided spacefractional partial differential equations, Appl. Numer. Math., 56(1) (2006) 80-90 [23] Z. Ding, A. Xiao, M. Li, Weighted finite difference methods for a class of space fractional partial differential equations with variable coefficients, J. Comput. Appl. Math., 233(8) (2010) 1905-1914 [24] H. Wang, N. Du, Fast alternating-direction finite difference methods for three-dimensional space-fractional diffusion equations, J. Comput. Phys., 258 (2014) 305-318 [25] E.H. Doha, A.H. Bhrawy and S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equations, Appl. Math. Model., 36 (2011) 4931-4943 [26] M. Parvizi, M. R. Eslahchi, M. Dehghan, Numerical solution of fractional advection-diffusion equation with a nonlinear source term, Numer. Algor., (2014), DOI: 10.1007/s11075-0149863-7. [27] M.R. Eslahchi, M. Dehghan, M. Parvizi, Application of the collocation method for solving nonlinear fractional integro-differential equations, J. Comput. Appl. Math., 257 (2014) 105128. [28] M. Dehghan, S. Abdi-mazraeh, M. Lakestani, Numerical solution for a class of fractional convection-diffusion equation using the flatlet oblique multiwavelets, J. Vib. Contr., 20 (2014) 913-924. [29] A.H. Bhrawy, M.A. Zaky, A method based on the Jacobi tau approximation for solving multi-term time-space fractional partial differential equations, J. Comput. Phys., 281 (2015) 876-895 [30] M. Zayernouri, M. Ainsworth, G.E. Karniadakis, A unified Petrov-Galerkin spectral method for fractional PDEs, Comput. Meth. Appl. Mech. Engineer., 283 (2015), 1545-1569 [31] V.R. Hosseini, W. Chen, Z.Avazzadeh, Numerical solution of fractional telegraph equation by using radial basis functions, Eng. Anal. Boundary Elem. 38 (2014) 31-39.
32
[32] A. Pedas, E. Tamme, Spline collocation for nonlinear fractional boundary value problems, Appl. Math. Comput., 244, (2014) 502-513 [33] M.H. Heydari, M.R. Hooshmandasl, F.M. Maalek Ghaini, C. Cattani, Wavelets method for the time fractional diffusion-wave equation, Phys. Letters A, 379 (2015) 71-76 [34] A.H. Bhrawy, E.H. Doha, D. Baleanu, and S.S. Ezz-Eldien, A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations, J. Comput. Phys., (2014) DOI: 10.1016/j.jcp.2014.03.039 [35] E.H. Doha, A.H. Bhrawy, M.A. Abdelkawy, R.A. Van Gorder, Jacobi-Gauss-Lobatto collocation method for the numerical solution of 1 + 1 nonlinear Schr¨ odinger equations, J. Comput. Phys., 261 (2014) 244-255 [36] L.W. Zhang, K.M. Liew, An element-free based solution for nonlinear Schr¨odinger equations using the ICVMLS-Ritz method, Appl. Math. Comput., 249 (2014) 333-345 [37] M. Naber, Time fractional Schr¨ oinger equation, J. Math. Phys., 45 (2004) 3339-3352. [38] S.Z. Rida, H.M. El-Sherbiny, A.A.M. Arafa, On the solution of the fractional nonlinear Schr¨ oinger equation, Phys. Lett. A, 372 (2008) 553-558 [39] M.A.E. Herzallah, Khaled A. Gepreel, Approximate solution to the time-space fractional cubic nonlinear Schr¨ odinger equation, Appl. Math. Modell., 36 (2012) 5678-5685 [40] L. Wei, Y. He, X. Zhang, S. Wang, Analysis of an implicit fully discrete local discontinuous Galerkin method for the time-fractional Schr¨odinger equation, Finite Elem Anal Des, 59 (2012) 28-34 [41] A. Mohebbi, M. Abbaszadeh, M. Dehghan, The use of a meshless technique based on collocation and radial basis functions for solving the time fractional nonlinear Schr¨odinger equation arising in quantum mechanics, Engineering Analysis with Boundary Elements, 37 (2013) 475-485 [42] M. Dehghan, M. Abbaszadeh, A. Mohebbi, An implicit RBF meshless approach for solving the time fractional nonlinear sine-Gordon and Klein-Gordon equations, Engineering Analysis with Boundary Elements, 50 (2015) 412-434 [43] F.B. Addaa, J. Cresson, Fractional differential equations and the Schr¨odinger equation, Appl. Math. Comput., 161, (2005) 323-345
33
[44] P. Wang, C. Huang, A conservative linearized difference scheme for the nonlinear fractional Schr¨ odinger equation, Numer.l Algor., 2014, doi:10.1007/s11075-014-9917-x [45] J. Dong, M. Xu, Space-time fractional Schr¨ odinger equation with time-independent potentials, J. Math. Anal. Appl., 344, (2008) 1005-1017 [46] R. Eid, S.I. Muslih, D. Baleanu, E. Rabei, On fractional Schr¨odinger
equation in α-
dimensional fractional space, Nonl. Analysis: Real World Applications, 10 (2009) 1299-1304 [47] A. Ashyralyev, B. Hicdurmaz, On the numerical solution of fractional Schr¨ odinger differential equations with the Dirichlet condition, Inter. J. Comput. Math., 89 (2012) 1927-1936. [48] A.H. Bhrawy, E.H. Doha, S.S. Ezz-Eldien, and R.A. Van Gorder, A new Jacobi spectral collocation method for solving 1+1 fractional Schr¨odinger equations and fractional coupled Schr¨ odinger systems, Eur. Phys. J. Plus, 129 (2014) 260 [49] K. Miller, B. Ross, An Introduction to the Fractional Calaulus and Fractional Differential Equations, John Wiley & Sons Inc., New York, (1993)
34