JacobiâGaussâLobatto collocation method for the ... - Semantic Scholar
Journal of Computational Physics 261 (2014) 244–255
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Jacobi–Gauss–Lobatto collocation method for the numerical solution of 1 + 1 nonlinear Schrödinger equations E.H. Doha a , A.H. Bhrawy b,c , M.A. Abdelkawy c , Robert A. Van Gorder d,∗ a
Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt d Department of Mathematics, University of Central Florida, Orlando, FL 32816, USA b c
a r t i c l e
i n f o
Article history: Received 18 July 2013 Received in revised form 30 December 2013 Accepted 3 January 2014 Available online 8 January 2014 Keywords: Nonlinear complex Schrödinger equations Gross–Pitaevskii equation Collocation method Jacobi–Gauss–Lobatto quadrature Implicit Runge–Kutta method
a b s t r a c t A Jacobi–Gauss–Lobatto collocation (J-GL-C) method, used in combination with the implicit Runge–Kutta method of fourth order, is proposed as a numerical algorithm for the approximation of solutions to nonlinear Schrödinger equations (NLSE) with initialboundary data in 1 + 1 dimensions. Our procedure is implemented in two successive steps. In the first one, the J-GL-C is employed for approximating the functional dependence on the spatial variable, using ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation which depends upon two general Jacobi parameters. The resulting equations together with the two-point boundary conditions induce a system of 2( N − 1) first-order ordinary differential equations (ODEs) in time. In the second step, the implicit Runge–Kutta method of fourth order is applied to solve this temporal system. The proposed J-GL-C method, used in combination with the implicit Runge–Kutta method of fourth order, is employed to obtain highly accurate numerical approximations to four types of NLSE, including the attractive and repulsive NLSE and a Gross–Pitaevskii equation with space-periodic potential. The numerical results obtained by this algorithm have been compared with various exact solutions in order to demonstrate the accuracy and efficiency of the proposed method. Indeed, for relatively few nodes used, the absolute error in our numerical solutions is sufficiently small. © 2014 Elsevier Inc. All rights reserved.
1. Introduction Several numerical methods have been intensively investigated for the numerical solution of time dependent partial differential equations (PDEs) (see for instance, [1–5]), particularly the NLSE. The discontinuous Galerkin method was proposed in [6] for the numerical solution of single and coupled NLSE. Dehghan and Taleei [7] proposed the Chebyshev pseudospectral differentiation matrices to reduce the solution of NLSE to a system of algebraic equations and very good results were obtained concerning the behavior of the pseudo-spectral solutions. In [8], the combination of multi-domain Chebyshev collocation methods and an explicit Runge–Kutta scheme of order four for solving the coupled version of NLSE was discussed. In [9], a fourth order explicit scheme for solving the coupled version of NLSE was considered, while the Galerkin method for the same problem was considered in [10]. Recently, Xu and Zhang [11] introduced and developed alternating direction implicit methods for numerical treatment of two-dimensional cubic NLSE with studying the stability analysis. Meanwhile, the Pade approximation was considered in [12] for the solution of a higher-order NLSE with fourth order dispersion. More
*
Corresponding author. E-mail addresses: [email protected] (E.H. Doha), [email protected] (A.H. Bhrawy), [email protected] (M.A. Abdelkawy), [email protected] (R.A. Van Gorder). 0021-9991/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2014.01.003
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
245
recently, Duan and Rong [13] discussed and developed a multi-quadric quasi-interpolation scheme for solving NLSE. Many other methods have been studied in the literature; see, for example, [14–16]. Spectral methods (see [17–21]) are often an efficient and highly accurate scheme when compared with local methods. There are three primary types of spectral methods based on the choice of test functions, namely, the Galerkin, tau and collocation methods. In these methods, one reduces the original PDEs into a system of ODEs (often in the time variable) by assuming a functional dependence on the spatial part of the equation. In such algorithms, one then solves for the time evolution of the solutions. It is within the realm of nonlinear and variable coefficient problems that the advantages of these collocation methods become apparent. Indeed, the projection operators involved in the Galerkin and tau methods lead to complications in deriving the required system of ODEs, from which the final solution is obtained. In both these methods, we require that the projection of the residual onto some polynomial space be zero. In contrast, the collocation method obtains the system of ODEs by requiring that the residual vanishes at a certain set of grid points, typically some set of Gauss-type points. As a result, while the Galerkin and tau methods may become untenable for strongly nonlinear problems, the collocation method has been shown in many cases to deal with linear, variable coefficient, and nonlinear problems with equal ease [22–27]. Along these lines, a number of related methods have been proposed for solving related equations [28–34]. The Jacobi polynomials are a class of orthogonal polynomials, which satisfy an orthogonality condition on the interval [−1, 1] with respect to the weight function (1 − x)θ (l + x)ϑ . Several orthogonal polynomials can obtained as special case of the Jacobi polynomials, such as the Legendre, Chebyshev, and Gegenbauer [35]. In recent years, the use of Jacobi polynomials for solving differential equations through collocation techniques has gained increasing popularity due to the ability to obtain the solution in terms of the general Jacobi parameters θ and ϑ (see [36,37]). The main concern of this paper is to extend the J-GL-C method to the solution of various types of NLSE. In particular, we shall consider the general initial-boundary value problem
2
i ∂t ψ( y , t ) + ∂ y y ψ( y , t ) + 2γ ψ( y , t ) ψ( y , t ) − 2δ R ( y , t )ψ( y , t ) = 0,
(1.1)
where
( y , t ) ∈ D × [0, T ],
D = { y: A y B },
with the boundary-initial conditions
ψ( A , t ) = ζ1 (t ), ψ( y , 0) = ξ1 ( y ),
ψ( B , t ) = ζ2 (t ), y ∈ D.
(1.2) (1.3)
Here we have maintained a general potential R ( y , t ) and the cubic nonlinearity. This type of equation is also referred to as a Gross–Pitaevskii equation when we have non-zero R ( y , t ). The nonlinear PDEs of type (1.1) are collocated only for the space variable at ( N − 1) points, for suitable collocation points. We use the ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation, which depends upon the two general parameters θ and ϑ . These equations, together with two-point boundary conditions, constitute system of (2N − 2) ODEs in the time variable. This system can then be solved by the implicit Runge–Kutta method of fourth order. Finally, the accuracy of the proposed method is demonstrated by four test problems, which correspond to physically meaningful cases. To the best of our knowledge, there are no results on the J-GL-C method for solving NLSE. This paper is arranged as follows. In Section 2, we derive the numerical method. First, we introduce some properties of Jacobi polynomials and their associated recursion relations. Then the Gauss–Lobatto collocation technique for a general form of the NLSE obtained using the Jacobi polynomials, and the implicit Runge–Kutta method of fourth order is outlined for the resulting system of ODEs. In Section 3, the proposed method is applied to four specific NLSE. It is here where we demonstrate the usefulness of the method, in terms of efficiency and accuracy. Indeed, for relatively few nodes, we obtain rather good results. In fact, for all examples considered, we obtain errors of order 10−6 or lower using only nine nodes. Finally, some concluding remarks are given in Section 4. 2. Derivation of the numerical method In the present section, we derive the numerical method. First, we remind the reader of some useful properties of Jacobi polynomials. Then, we develop the J-GL-C method to numerically solve the NLSE (1.1). We conclude this section with the implicit Runge–Kutta method used to solve the time-dependent ODEs. 2.1. Some properties of the Jacobi polynomials We collected in this section some basic knowledge of Jacobi polynomials that are most relevant to spectral approximations [35]. It includes Legendre and Chebyshev polynomials as two special cases, so it is worthwhile work with general Jacobi polynomials. A basic property of the Jacobi polynomials is that they are the eigenfunctions to a singular Sturm– Liouville problem:
246
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
1 − x2 φ (x) + ϑ − θ + (θ + ϑ + 2)x φ (x) + n(n + θ + ϑ + 1)φ(x) = 0.
(2.1)
The following recurrence relation generate the Jacobi polynomials:
(θ,ϑ)
(θ,ϑ)
J k+1 (x) = ak (θ,ϑ)
J0
(θ,ϑ) (θ,ϑ)
x − bk (θ,ϑ)
(x) = 1,
J1
Jk
(θ,ϑ) (θ,ϑ) J k−1 (x),
(x) − ck
1
1
2
2
k 1,
(x) = (θ + ϑ + 2)x + (θ − ϑ),
where
(2k + θ + ϑ + 1)(2k + θ + ϑ + 2) , 2(k + 1)(k + θ + ϑ + 1) (ϑ 2 − θ 2 )(2k + θ + ϑ + 1) (θ,ϑ) bk = , 2(k + 1)(k + θ + ϑ + 1)(2k + θ + ϑ) (k + θ)(k + ϑ)(2k + θ + ϑ + 2) (θ,ϑ) ck = . (k + 1)(k + θ + ϑ + 1)(2k + θ + ϑ) (θ,ϑ)
ak
=
The Jacobi polynomials are satisfying the following identities: (θ,ϑ)
Jk
(θ,ϑ)
(−x) = (−1)k J k
(θ,ϑ)
(x),
Jk
(−1) =
(−1)k Γ (k + ϑ + 1) . k!Γ (ϑ + 1) (θ,ϑ)
Moreover, the q derivative of Jacobi polynomials of degree k ( J k (θ,ϑ)
D (q) J k
(x) =
(2.2)
(x)), can be obtained from:
Γ ( j + θ + ϑ + q + 1) (θ +q,ϑ+q) J (x). 2q Γ ( j + θ + ϑ + 1) k−q
(2.3)
Let w (θ,ϑ) (x) = (1 − x)θ (1 + x)ϑ , then we define the weighted space L 2w (θ,ϑ) as usual. The inner product and the norm of
L 2w (θ,ϑ)
with respect to the weight function are defined as:
1 (u , v ) w (θ,ϑ) =
u (x) v (x) w (θ,ϑ) (x) dx,
1
u w (θ,ϑ) = (u , u ) w2 (θ,ϑ) .
(2.4)
−1
The set of Jacobi polynomials forms a complete L 2w (θ,ϑ) -orthogonal system, and
(θ,ϑ) J k
w (θ,ϑ)
= hk =
2θ +ϑ+1 Γ (k + θ + 1)Γ (k + ϑ + 1)
(2k + θ + ϑ + 1)Γ (k + 1)Γ (k + θ + ϑ + 1)
(2.5)
.
2.2. Jacobi spectral collocation method The main objective of this section is to develop the J-GL-C method to numerically solve the NLSE (1.1). In order to begin the numerical solution, we put the complex functions ψ in real and imaginary parts as:
ψ( y , t ) = u ( y , t ) + i v ( y , t ),
(2.6)
where u and v are real functions, from Eqs. (1.1) and (2.6), we get
i ∂t u ( y , t ) + ∂ y y v ( y , t ) + 2γ u 2 ( y , t ) + v 2 ( y , t ) v ( y , t ) − 2δ R ( y , t ) v ( y , t )
− ∂t u ( y , t ) − ∂ y y v ( y , t ) − 2γ u 2 ( y , t ) + v 2 ( y , t ) u ( y , t ) + 2δ R ( y , t ) u ( y , t ) = 0,
(2.7)
which can be separated as coupled time-dependent nonlinear PDEs,
∂t u ( y , t ) + ∂ y y v ( y , t ) + 2γ u 2 ( y , t ) + v 2 ( y , t ) v ( y , t ) − 2δ R ( y , t ) v ( y , t ) = 0, ∂t v ( y , t ) − ∂ y y u ( y , t ) − 2γ u 2 ( y , t ) + v 2 ( y , t ) u ( y , t ) + 2δ R ( y , t ) u ( y , t ) = 0,
(2.8)
with the boundary-initial conditions
u ( A , t ) = g 1 (t ),
u ( B , t ) = g 2 (t ),
v ( A , t ) = g 3 (t ),
v ( B , t ) = g 4 (t ),
u ( y , 0) = f 1 ( y ),
v ( y , 0) = f 2 ( y ),
y ∈ D.
(2.9)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
247
Now, suppose the change of variables
x=
2 B−A
y+
A+B
u ( y , t ) = w (x, t ),
,
A−B
v ( y , t ) = z(x, t ),
which will be used to transform problem (1.1)–(1.3) into another one in the classical interval, [−1, 1], for the space variable, to directly implement collocation method based on Jacobi family defined on [−1, 1],
2
2
∂t w (x, t ) +
B−A
2
∂t z(x, t ) −
2
B−A
∂xx z(x, t ) + 2γ w 2 (x, t ) + z2 (x, t ) z(x, t ) − 2δr (x, t ) z(x, t ) = 0,
∂xx w (x, t ) − 2γ w 2 (x, t ) + z2 (x, t ) w (x, t ) + 2δr (x, t ) w (x, t ) = 0,
(2.10)
where
(x, t ) ∈ D ∗ × [0, T ],
D ∗ = {x: −1 x 1},
with the boundary-initial conditions
w (−1, t ) = g 5 (t ),
w (1, t ) = g 6 (t ),
z(−1, t ) = g 7 (t ),
z(1, t ) = g 8 (t ),
w (x, 0) = f 3 (x),
z(x, 0) = f 4 (x),
x ∈ D∗.
(2.11)
The main advantage of the use of the node points of Jacobi–Gauss–Lobatto quadrature is that the distribution of these points in [−1, 1] can be controlled by the two Jacobi parameters, θ and ϑ . Now, we outline the main step of implementing J-GL-C method for reducing nonlinear coupled system (2.10)–(2.11) to system of ODEs. The expansion of functions w (x, t ) (θ,ϑ) and z(x, t ), using the Jacobi polynomial, J j (x), takes the form
w (x, t ) =
N
(θ,ϑ)
a j (t ) J j
(x),
(2.12)
j =0
z(x, t ) =
N
(θ,ϑ)
b j (t ) J j
(x),
(2.13)
j =0
and in virtue of (2.4)–(2.5), we deduce that
a j (t ) =
b j (t ) =
1 hj 1 hj
1
(θ,ϑ)
w (x, t ) w (θ,ϑ) (x) J j
(x) dx,
−1
1
(θ,ϑ)
z(x, t ) w (θ,ϑ) (x) J j
(x) dx.
(2.14)
−1
We can circumvent the need for evaluating the previous integrals, by using Jacobi–Gauss–Lobatto quadrature formula to approximate the integrals. For any φ ∈ S 2N +1 (−1, 1),
1 w (θ,ϑ) (x)φ(x) dx =
N
(θ,ϑ) N(θ,ϑ) , j φ xN , j ,
(2.15)
j =0
−1
(θ,ϑ)
(θ,ϑ)
where S N (−1, 1) is the set of polynomials of degree less than or equal to N, x N , j (0 j N) and N , j (0 j N) are the nodes and the corresponding Christoffel numbers of the Jacobi–Gauss–Lobatto quadrature formula on the interval (−1, 1), respectively. Due to the orthogonality property of Jacobi polynomials (2.4), the expansion coefficients a j (t ), b j (t ) in terms of the solution at the Jacobi–Gauss–Lobatto grid points can be approximated by
a j (t ) =
N 1 (θ,ϑ) (θ,ϑ) (θ,ϑ) (α ,β) Jj x N ,i N ,i w x N ,i , t , hj i =0
b j (t ) =
N 1 (θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i N ,i z x N ,i , t . hj i =0
(2.16)
248
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Accordingly, (2.12)–(2.13) can be given by
w (x, t ) =
N N 1 i =0
z(x, t ) =
j =0
hj
N N 1 i =0
j =0
(θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i Jj (x) N ,i
hj
(θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i Jj (x) N ,i
(θ,ϑ)
w x N ,i , t ,
(θ,ϑ)
z x N ,i , t .
(2.17)
Furthermore, if we differentiate (2.17) once using (2.3), and evaluate it at all J-GL-C points, it is easy to compute the first spatial partial derivative in terms of the values at theses collocation points as
(θ,ϑ)
w x x N ,n , t =
N
(θ,ϑ)
A ni w x N ,i , t ,
i =0
(θ,ϑ)
zx x N ,n , t =
N
(θ,ϑ)
n = 0, 1 , . . . , N ,
B ni z x N ,i , t ,
(2.18)
i =0
where
A ni =
N j+θ +ϑ +1
2h j
j =0
B ni =
N j+θ +ϑ +1
2h j
j =0
(θ,ϑ) (θ,ϑ) (θ +1,ϑ+1) (θ,ϑ) (θ,ϑ) x N ,i J j −1 x N ,n N ,i ,
Jj
(θ,ϑ) (θ,ϑ) (θ +1,ϑ+1) (θ,ϑ) (θ,ϑ) x N ,i J j −1 x N ,n N ,i .
Jj
(2.19)
Similar steps can be applied to the second spatial partial derivative to get
(θ,ϑ)
w xx x N ,n , t =
N
(θ,ϑ)
C ni w x N ,i , t ,
i =0
(θ,ϑ)
zxx x N ,n , t
=
N
(θ,ϑ)
D ni z x N ,i , t ,
n = 0, 1 , . . . , N ,
(2.20)
i =0
where
C ni =
N ( j + θ + ϑ + 2)( j + θ + ϑ + 1)
4h j
j =0
D ni =
N ( j + θ + ϑ + 2)( j + θ + ϑ + 1)
4h j
j =0
(θ,ϑ) (θ,ϑ) (θ +2,ϑ+2) (θ,ϑ) (θ,ϑ) x N ,i J j −2 x N ,n N ,i ,
Jj
(θ,ϑ) (θ,ϑ) (θ +2,ϑ+2) (θ,ϑ) (θ,ϑ) x N ,i J j −2 x N ,n N ,i .
Jj
(2.21)
For the spectral J-GL-C method, the residual of (1.1) is set to zero at the N − 1 Jacobi–Gauss–Lobatto points. Also, the boundary conditions (1.2) can be imposed at the two collocation points −1 and 1 in the expansion of functions. Therefore, adopting (2.8)–(2.11), we may write (1.1)–(1.2) in the form
˙ n (t ) = − w
where
2 N
B−A
z˙ n (t ) =
2
2
(θ,ϑ)
i =0
2 N
B−A
D ni zi (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
C ni w i (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
n = 1, 2, . . . , N − 1,
i =0
w k (t ) = w x N ,k , t ,
(θ,ϑ)
zk (t ) = z x N ,k , t ,
k = 1, . . . , N − 1, n = 1, . . . , N − 1.
This provides two ( N − 1) system of first order ODEs in the expansion coefficients a j (t ), namely
(2.22)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
˙ n (t ) = − w z˙ n (t ) =
2 N −1
2 B−A
B−A
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
i =1
2 N −1
2
249
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
(2.23)
i =1
where
dn (t ) = D n0 g 7 (t ) + D nN g 8 (t ), dn (t ) = C n0 g 5 (t ) + C nN g 6 (t ). We get the following two system of nonlinear ODEs:
˙ n (t ) = − w z˙ n (t ) =
2 N −1
2 B−A
B−A
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
i =1
2 N −1
2
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
i =1
n = 1, 2, . . . , N − 1,
(2.24)
subject to the initial values
(θ,ϑ)
w n (0) = f 3 x N ,n
(θ,ϑ) f 4 x N ,n ,
zn (0) =
, n = 1, . . . , N − 1.
(2.25)
Finally, (2.24)–(2.25) can be rewritten into two matrix form of N − 1 first order ODEs with their vectors of initial values:
˙ (t ) = F t , w (t ), z(t ) , w z˙ (t ) = G t , w (t ), z(t ) , w(0) = f3 , z(0) = f4 ,
(2.26)
where
T
˙ (t ) = w ˙ 1 (t ), w ˙ 2 (t ), . . . , w ˙ N −1 (t ) w
, T z˙ (t ) = z˙ 1 (t ), z˙ 2 (t ), . . . , z˙ N −1 (t ) , (θ,ϑ) T (θ,ϑ) (θ,ϑ) f3 = f 3 x N ,1 , f 3 x N ,2 , . . . , f 3 x N , N −1 , (θ,ϑ) T (θ,ϑ) (θ,ϑ) f4 = f 4 x N ,1 , f 4 x N ,2 , . . . , f 4 x N , N −1 , T F t , w (t ), z(t ) = F 1 t , w (t ), z(t ) , F 1 t , w (t ), z(t ) , . . . , F N −1 t , w (t ), z(t ) ,
and
T
G t , w (t ), z(t ) = G 1 t , w (t ), z(t ) , G 1 t , w (t ), z(t ) , . . . , G N −1 t , w (t ), z(t )
.
Here,
F n t , w (t ), z(t ) = −
G n t , w (t ), z(t ) =
2 N −1
2 B−A 2
B−A
i =1
2 N −1
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ).
i =1
The two N − 1 dimensional system of ODEs (2.26) can be combined to be written in a 2N − 2 system of ODEs and then can be solved by using an implicit Runge–Kutta formula of fourth order.
250
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
2.3. Implicit Runge–Kutta (IRK) method The implicit Runge–Kutta (IRK) method is one of the suitable methods for solving a system of ODEs of first order like we have obtained above. Since the IRK methods have excellent stability properties which reduce computational costs. Consider
˙ (t ) = H t , u (t ) u u(0) = f.
(2.27)
Then, the iterative solution is
Ui , j = ui + t
m
a j ,l H (t i + cl t , Ui ,l ),
j = 1, . . . , m ,
l =1
ui +1 = ui + t
m
b j H (t i + c j t , Ui , j ),
∀j
(2.28)
j =1
here t is the step size, t i = t 0 + i t and ui is approximations to the exact solution u(t i ). Moreover, ai j , b j and c j are real coefficients, i , j = 1, 2, . . . , m, which can be represented by means of the Butcher array,
c 1 a11 · · · a1s
.. .
.. .
c s as1 b1
..
. ··· ···
.. . ass bs
3. Numerical examples This section considers four numerical examples to show the accuracy and applicability of the proposed method in the present paper. Comparison of the results obtained by various choices of Jacobi parameters θ and ϑ reveals that the present method is very effective and convenient for various choices of θ and ϑ . We focus on some well-known equations from mathematical physics in order to demonstrate the method. First we consider an attractive NLSE [38,39], and for our next example we consider a repulsive NLSE [38,39]. Third, we consider the Gross–Pitaevski equation [40,41], which describes the ground state of a quantum system of identical bosons. We assume that the potential is periodic in space, consistent with a Bose–Einstein condensate. Finally, we consider a NLSE with a linear instability added. In each of these examples, we shall highlight the accuracy and efficiency of the method. 3.1. Numerical solution for an attractive NLSE We first consider an attractive NLSE [39] in the form:
i ∂t ψ + ∂xx ψ + 2|ψ|2 ψ = 0,
( y, t ) ∈ D × R
(3.1)
with the boundary-initial conditions
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e i (t + A ) ,
ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e i (t + B ) ,
ψ(x, 0) = u (x, 0) + i v (x, 0) = e i (x) .
(3.2)
In this example, we consider D = [ A , B ] = [−1, 1]. Furthermore, we shall approximate the time domain by [0, T ] for numerical purposes. The difference between the measured or inferred value of a numerical solution and its actual value (absolute error) is given by (for real part, imaginary part, and complex modulus)
E 1 (x, t ) = u (x, t ) − u˜ (x, t ), E 2 (x, t ) = v (x, t ) − v˜ (x, t ),
˜ x, t ), E 3 (x, t ) = ψ(x, t ) − ψ(
(3.3)
˜ x, t ) = u˜ (x, t ) + i v˜ (x, t ) are the exact solution and the approximate solution at the where ψ(x, t ) = u (x, t ) + i v (x, t ) and ψ( point (x, t ), respectively. It may be shown that the error norms are related. Indeed,
E 3 (x, t ) = u (x, t ) − u˜ (x, t ) + i v (x, t ) − v˜ (x, t ) =
E 1 (x, t )2 + E 2 (x, t )2 .
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
251
Table 1 Maximum absolute errors with various choices of N, θ and ϑ for Example 1. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
3.81 × 10−1 3.62 × 10−5 4.62 × 10−8
5.73 × 10−1 15.70 × 10−5 22.59 × 10−8
6.21 × 10−1 16.04 × 10−5 22.74 × 10−8
− 12
− 12
2.30 × 10−1 5.75 × 10−5 4.93 × 10−8
4.83 × 10−1 3.65 × 10−5 8.30 × 10−8
4.83 × 10−2 5.75 × 10−5 9.21 × 10−8
2 6 10
1 2
1 2
4.40 × 10−1 5.85 × 10−5 6.41 × 10−8
6.38 × 10−1 22.07 × 10−5 33.44 × 10−8
6.67 × 10−1 22.12 × 10−5 33.64 × 10−8
1 2
0
4.13 × 10−1 7.66 × 10−5 4.86 × 10−8
5.37 × 10−1 21.36 × 10−5 25.47 × 10−8
6.78 × 10−2 22.17 × 10−5 25.58 × 10−8
Fig. 1. The absolute error (a) E 1 (x, t ) and (b) E 2 (x, t ) of problem (3.1) where θ = 0, ϑ = − 12 and N = 10. Note that while errors gradually increase with time, they remain less than 10−6 over the plotted domain.
Moreover, the maximum absolute error is given by
M 1 = Max E 1 (x, t ): ∀(x, t ) ∈ D × [0, T ] , M 2 = Max E 2 (x, t ): ∀(x, t ) ∈ D × [0, T ] , M 3 = Max E 3 (x, t ): ∀(x, t ) ∈ D × [0, T ] . (3.4) Clearly, M 3 M 12 + M 22 . The maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) related to the problem (3.1)–(3.2) are presented in Table 1 using the J-GL-C method with various choices of N, θ and ϑ . The errors are calculated through a comparison with the exact solution
ψ(x, t ) = e i (x+t ) .
(3.5)
We see in this table that the results are very accurate for even small choices of the number of nodes, N. In Fig. 1, we see that the absolute errors E 1 (x, t ) and E 2 (x, t ) are very small, despite the relatively small number of grid points used. This numerical experiment demonstrates the utility of the method.
252
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Table 2 Maximum absolute errors with various choices of N, θ and ϑ for Example 2. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
2.43 × 10−1 6.00 × 10−5 2.26 × 10−7
2.69 × 10−1 5.68 × 10−5 1.70 × 10−7
3.23 × 10−1 6.41 × 10−5 2.35 × 10−7
− 12
− 12
2.57 × 10−1 4.19 × 10−5 2.52 × 10−7
6.40 × 10−2 5.81 × 10−5 1.92 × 10−7
2.64 × 10−1 6.25 × 10−5 2.89 × 10−7
2 6 10
0
− 12
2.71 × 10−1 11.11 × 10−5 2.19 × 10−7
1.98 × 10−1 6.60 × 10−5 1.79 × 10−7
3.35 × 10−1 12.91 × 10−5 2.28 × 10−7
− 12
0
2.12 × 10−1 5.40 × 10−5 2.42 × 10−7
1.33 × 10−1 6.64 × 10−5 1.93 × 10−7
2.36 × 10−1 6.98 × 10−5 2.43 × 10−7
In the case where an exact solution is not known, one follows a similar method. The difference is that one calculates the residual errors only, since comparison to an exact solution is not possible. 3.2. Numerical solution for a repulsive NLSE Consider the repulsive NLSE [39]
i ∂t ψ + ∂xx ψ − 2γ |ψ|2 ψ = 0,
(3.6)
with the boundary-initial conditions
ψ(−1, t ) = u (−1, t ) + i v (−1, t ) = e i (−1−3t ) , ψ(1, t ) = u (1, t ) + i v (1, t ) = e i (1−3t ) , ψ(x, 0) = u (x, 0) + i v (x, 0) = e i x .
(3.7)
Note that the exact solution is given by
ψ(x, t ) = e i (x−3t ) .
(3.8)
Table 2 lists the maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) of problem (3.6)–(3.7) using the J-GL-C method for with various choices of N, θ and ϑ . The numerical results presented in this table show that the results are very accurate for small value of N. Fig. 2 demonstrates that the absolute errors E 1 (x, t ) and E 2 (x, t ) are very small for even the small number of grid points taken. 3.3. Nonlinear Gross–Pitaevskii equation with space-periodic potential Consider the complex nonlinear Gross–Pitaevskii (CNGP) equation with a space-periodic potential [39–41] given by
i ∂t ψ +
1 2
∂xx ψ − |ψ|2 ψ − cos(x)ψ = 0,
(3.9)
with the boundary-initial conditions
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e
−3it
sin( A ),
2
−3it 2
sin( B ),
ψ(x, 0) = u (x, 0) + i v (x, 0) = sin(x).
(3.10)
The exact solution of CNGP equation is given by
ψ(x, t ) = e
−3it 2
sin(x).
(3.11)
To demonstrate the accuracy of the proposed method, in Table 3 we introduce maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) of CNGP equation, using the J-GL-C method with various choices of θ and ϑ and N. We see that the approximate solutions u˜ (x, t ) ( v˜ (x, t )) and the exact solution u (x, t ) (v (x, t )) are coincided for different values of t and x. 3.4. An attractive NLSE with linear instability We consider an attractive NLSE [39] with an additional linear term
i ∂t ψ + ∂xx ψ + 2|ψ|2 ψ − 2γ ψ = 0, with the boundary-initial conditions
(3.12)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
253
Fig. 2. The absolute error (a) E 1 (x, t ) and (b) E 2 (x, t ) of problem (3.6) where θ = ϑ = − 12 and N = 10. Table 3 Maximum absolute errors with various choices of N, θ and ϑ for Example 3. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
4.29 × 10−2 3.84 × 10−6 8.32 × 10−8
4.00 × 10−2 3.28 × 10−6 7.98 × 10−8
5.87 × 10−2 3.97 × 10−6 10.49 × 10−8
− 12
− 12
5.46 × 10−2 4.90 × 10−6 3.05 × 10−7
4.00 × 10−2 5.56 × 10−6 1.42 × 10−7
5.87 × 10−2 6.84 × 10−6 3.05 × 10−7
2 6 10
1 2
0
5.25 × 10−2 3.67 × 10−6 2.11 × 10−7
4.66 × 10−2 9.26 × 10−6 1.69 × 10−7
7.02 × 10−2 9.26 × 10−6 2.43 × 10−7
1 2
1 2
4.29 × 10−2 3.97 × 10−6 1.01 × 10−7
4.00 × 10−2 5.30 × 10−6 7.13 × 10−8
5.87 × 10−2 5.40 × 10−6 1.24 × 10−7
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e i ( A +(1−2γ )t ) , ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e i ( B +(1−2γ )t ) ψ(x, 0) = u (x, 0) + i v (x, 0) = e ix .
(3.13)
Maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) related to (3.12)–(3.13) are introduced in Table 4 using J-GL-C method for with various choices of N, θ and ϑ , where γ = 1 compared with the exact solution
ψ(x, t ) = e i (x+(1−2γ )t ) .
(3.14)
Again and from Table 4, we clearly observe that the global scheme (J-GL-C) is accurate for even small values of N. 4. Conclusion We have proposed an accurate numerical algorithm based on the J-GL-C method for solving complex NLSE. Our procedure is implemented in two successive steps. First, the J-GL-C method is employed for approximating the functional dependence on the spatial variable, using ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation which depends upon two general Jacobi parameters. The resulting equations together with the two-point boundary conditions give us a system of 2( N − 1)
254
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Table 4 Maximum absolute errors with various choices of N, θ and ϑ for Example 4. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
6.12 × 10−1 1.53 × 10−4 1.82 × 10−7
2.89 × 10−1 6.80 × 10−4 1.40 × 10−7
6.21 × 10−1 1.60 × 10−4 2.15 × 10−7
− 12
− 12
4.01 × 10−1 4.19 × 10−4 2.27 × 10−7
2.69 × 10−1 5.69 × 10−4 1.04 × 10−7
4.83 × 10−1 5.76 × 10−4 2.27 × 10−7
2 6 10
0
1 2
7.43 × 10−1 2.09 × 10−4 2.07 × 10−7
5.85 × 10−2 1.26 × 10−4 1.59 × 10−7
7.43 × 10−1 2.09 × 10−4 2.35 × 10−7
1 2
1 2
6.59 × 10−1 2.18 × 10−4 2.64 × 10−7
1.60 × 10−1 7.10 × 10−4 1.96 × 10−8
6.67 × 10−1 2.21 × 10−4 3.02 × 10−7
first-order ODEs in the time variable. In the second part of the solution method, the implicit Runge–Kutta method of fourth order is applied to solve this system of 2( N − 1) ODEs in time. Note that it is possible to use other implicit Runge–Kutta formulas, say an implicit Runge–Kutta method of order greater than four, or even a second-order method. The main advantage of the use of the node points of Jacobi–Gauss–Lobatto quadrature as collocation points is that the distribution of these points in a specified interval can be tuned by the two Jacobi parameters, θ and ϑ . Through appropriate selection of these parameters, we may recover some of the well-known spectral collocation approximations. Namely, the Legendre and Chebyshev pseudo-spectral approximations can be derived as special cases from our algorithm if we suitably choose the corresponding spacial cases of Jacobi parameters θ and ϑ . From the numerical results given in Section 3, it has been concluded that the obtained results are excellent in terms of accuracy for all tested problems. For each of the examples given, we have been able to construct solutions with absolute errors between 10−7 and 10−6 , and greater accuracy may be obtained if we add more nodes. Note that the initial data considered has always been of modulus unity, so the relative errors should be of the same order of magnitude as the absolute errors. Concerning the efficiency of the method, note that we have obtained approximate numerical solutions of a desirable accuracy by taking only N = 10 (i.e., 9 nodes). Considering the size of the domains considered, this is approximately one node per 0.2 units of space. Considering the rather large average spatial gap, the results are rather good. Of course, adding more nodes will increase the accuracy further, and one should be able to employ this method in order to determine solutions of a desired accuracy. With this paper we have outlined the application of a J-GL-C method for solving complex NLSE. In principle, this algorithm may be extended to related problems, such as to coupled complex nonlinear Schrödinger equations. One might also consider higher order dispersive equations. We should note that, as a numerical method, we are restricted to solving problems over a finite domain. Hence, this method is particularly well suited for boundary value problems with finite spatial domains. Acknowledgement R.A.V. supported in part by NSF grant # 1144246. References [1] S. Islam, S. Haq, M. Uddin, A mesh free interpolation method for the numerical solution of the coupled nonlinear partial differential equations, Eng. Anal. Bound. Elem. 33 (2009) 399–409. [2] A. Mohebbi, Z. Asgari, M. Dehghan, Numerical solution of nonlinear Jaulent–Miodek and Whitham–Broer–Kaup equations, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 4602–4610. [3] R.C. Mittal, R.K. Jain, Numerical solutions of nonlinear Burgers’ equation with modified cubic B-splines collocation method, Appl. Math. Comput. 218 (2012) 7839–7855. [4] 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 Schrodinger equation arising in quantum mechanics, Eng. Anal. Bound. Elem. 37 (2013) 475–485. [5] Y. Gurefe, E. Misirli, A. Sonmezoglu, M. Ekici, Extended trial equation method to generalized nonlinear partial differential equations, Appl. Math. Comput. 219 (2013) 5253–5260. [6] H. Wang, Numerical studies on the split step finite difference method for the nonlinear Schrödinger equations, Appl. Math. Comput. 170 (2005) 17–35. [7] M. Dehghan, A. Taleei, Numerical solution of nonlinear Schrödinger equation by using time–space pseudo-spectral method, Numer. Methods Partial Differ. Equ. 26 (2010) 979–992. [8] M. Dehghan, A. Taleei, A Chebyshev pseudo-spectral multi-domain method for the soliton solution of coupled nonlinear Schrodinger equations, Comput. Phys. Commun. 182 (2011) 2519–2529. [9] M.S. Ismail, A fourth order explicit scheme for the coupled nonlinear Schrödinger equation, Appl. Math. Comput. 196 (2008) 273–284. [10] M.S. Ismail, Numerical solution of coupled nonlinear Schrödinger equation by Galerkin method, Math. Comput. Simul. 78 (2008) 532–547. [11] Y. Xu, L. Zhang, Alternating direction implicit method for solving two-dimensional cubic nonlinear Schrödinger equation, Comput. Phys. Commun. 183 (2012) 1082–1093. [12] M. Smadi, D. Bahloul, A compact split step Pade scheme for higher-order nonlinear Schrödinger equation (HNLS) with power law nonlinearity and fourth order dispersion, Comput. Phys. Commun. 182 (2011) 366–371. [13] Y. Duan, F. Rong, A numerical scheme for nonlinear Schrödinger equation by MQ quasi-interpolation, Eng. Anal. Bound. Elem. 37 (2013) 89–94. [14] H. Sakaguchi, T. Higashiuchi, Two-dimensional dark soliton in the nonlinear Schrödinger equation, Phys. Lett. A 359 (2006) 647–651. [15] Z. Kalogiratou, Th. Monovasilis, T.E. Simos, Numerical solution of the two-dimensional time independent Schrödinger equation with Numerov-type methods, J. Math. Chem. 37 (3) (2005) 271–279.
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
255
[16] M. Javidi, A. Golbabai, Numerical studies on nonlinear Schrödinger equations by spectral collocation method with preconditioning, J. Math. Anal. Appl. 333 (2007) 1119–1127. [17] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods: Fundamentals in Single Domains, Springer-Verlag, New York, 2006. [18] E.H. Doha, A.H. Bhrawy, An efficient direct solver for multidimensional elliptic Robin boundary value problems using a Legendre spectral-Galerkin method, Comput. Math. Appl. 64 (2012) 558–571. [19] S. Karimi Vanani, F. Soleymani, Tau approximate solution of weakly singular Volterra integral equations, Math. Comput. Model. 57 (3–4) (2013) 494–502. [20] W. Shao, X. Wu, Chebyshev tau meshless method based on the highest derivative for fourth order equations, Appl. Math. Model. 37 (3) (1 February 2013) 1413–1430. [21] E.H. Doha, W.M. Abd-Elhameed, M.A. Bassuony, New algorithms for solving high even-order differential equations using third and fourth Chebyshev– Galerkin methods, J. Comput. Phys. 236 (2013) 563–579. [22] C.I. Gheorghiu, Spectral Methods for Differential Problems, T. Popoviciu Institute of Numerical Analysis, Cluj-Napoca, Romaina, 2007. [23] C.I. Gheorghiu, M.E. Hochstenbach, B. Plestenjak, J. Rommes, Spectral collocation solutions to multiparameter Mathieu’s system, Appl. Math. Comput. 218 (24) (2012) 11990–12000. [24] A.H. Bhrawy, M.A. Alghamdi, A shifted Jacobi–Gauss–Lobatto collocation method for solving nonlinear factional Langevin equation involving two fractional orders in different intervals, Bound. Value Probl. 2012 (2012) 62. [25] K. Zhang, J. Li, H. Song, Collocation methods for nonlinear convolution Volterra integral equations with multiple proportional delays, Appl. Math. Comput. 218 (2012) 10848–10860. [26] A. Saadatmandi, M. Dehghan, The use of Sinc-collocation method for solving multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 593–601. [27] A.H. Bhrawy, A.S. Alofi, A Jacobi–Gauss collocation method for solving nonlinear Lane–Emden type equations, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 62–70. [28] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numer. Methods Partial Differ. Equ. 28 (2012) 74–93. [29] M.R. Eslahchi, M. Dehghan, S. Ahmadi-Asl, The general Jacobi matrix method for solving some nonlinear ordinary differential equations, Appl. Math. Model. 36 (2012) 3387–3398. [30] M.R. Eslahchi, M. Dehghan, S. Amani, The third and fourth kinds Chebyshev polynomials and best uniform approximation, Math. Comput. Model. 55 (2012) 1746–1762. [31] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Model. 53 (2011) 1865–1877. [32] M. Dehghan, R. Jazlanian, A fourth-order central Runge–Kutta scheme for hyperbolic conservation laws, Numer. Methods Partial Differ. Equ. 26 (2010) 1675–1692. [33] A. Mohebbi, M. Dehghan, The use of compact boundary value method for the solution of two-dimensional Schrödinger equation, J. Comput. Appl. Math. 225 (2009) 124–134. [34] M. Dehghan, D. Mirzaei, Numerical solution to the unsteady two-dimensional Schrödinger equation using meshless local boundary integral equation method, Int. J. Numer. Methods Eng. 76 (2008) 501–520. ˝ Orthogonal Polynomials, Colloquium Publications, vol. XXIII, American Mathematical Society, ISBN 978-0-8218-1023-1, 1939, MR 0372517. [35] G. Szego, [36] E.H. Doha, A.H. Bhrawy, Efficient spectral-Galerkin algorithms for direct solution of fourth-order differential equations using Jacobi polynomials, Appl. Numer. Math. 58 (2008) 1224–1244. [37] M. El-Kady, Jacobi discrete approximation for solving optimal control problems, J. Korean Math. Soc. 49 (2012) 99–112. [38] M.J. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 1981. [39] K. Mallory, R.A. Van Gorder, Stationary solutions for the 1 + 1 nonlinear Schrödinger equation modeling repulsive Bose–Einstein condensates in small potentials, Phys. Rev. E 88 (2013) 013205. [40] E.P. Gross, Structure of a quantized vortex in boson systems, Nuovo Cimento 20 (1961) 454–457. [41] L.P. Pitaevsk, Vortex lines in an imperfect Bose gas, Sov. Phys. JETP 13 (1961) 451–454.
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Jacobi–Gauss–Lobatto collocation method for the numerical solution of 1 + 1 nonlinear Schrödinger equations E.H. Doha a , A.H. Bhrawy b,c , M.A. Abdelkawy c , Robert A. Van Gorder d,∗ a
Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt d Department of Mathematics, University of Central Florida, Orlando, FL 32816, USA b c
a r t i c l e
i n f o
Article history: Received 18 July 2013 Received in revised form 30 December 2013 Accepted 3 January 2014 Available online 8 January 2014 Keywords: Nonlinear complex Schrödinger equations Gross–Pitaevskii equation Collocation method Jacobi–Gauss–Lobatto quadrature Implicit Runge–Kutta method
a b s t r a c t A Jacobi–Gauss–Lobatto collocation (J-GL-C) method, used in combination with the implicit Runge–Kutta method of fourth order, is proposed as a numerical algorithm for the approximation of solutions to nonlinear Schrödinger equations (NLSE) with initialboundary data in 1 + 1 dimensions. Our procedure is implemented in two successive steps. In the first one, the J-GL-C is employed for approximating the functional dependence on the spatial variable, using ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation which depends upon two general Jacobi parameters. The resulting equations together with the two-point boundary conditions induce a system of 2( N − 1) first-order ordinary differential equations (ODEs) in time. In the second step, the implicit Runge–Kutta method of fourth order is applied to solve this temporal system. The proposed J-GL-C method, used in combination with the implicit Runge–Kutta method of fourth order, is employed to obtain highly accurate numerical approximations to four types of NLSE, including the attractive and repulsive NLSE and a Gross–Pitaevskii equation with space-periodic potential. The numerical results obtained by this algorithm have been compared with various exact solutions in order to demonstrate the accuracy and efficiency of the proposed method. Indeed, for relatively few nodes used, the absolute error in our numerical solutions is sufficiently small. © 2014 Elsevier Inc. All rights reserved.
1. Introduction Several numerical methods have been intensively investigated for the numerical solution of time dependent partial differential equations (PDEs) (see for instance, [1–5]), particularly the NLSE. The discontinuous Galerkin method was proposed in [6] for the numerical solution of single and coupled NLSE. Dehghan and Taleei [7] proposed the Chebyshev pseudospectral differentiation matrices to reduce the solution of NLSE to a system of algebraic equations and very good results were obtained concerning the behavior of the pseudo-spectral solutions. In [8], the combination of multi-domain Chebyshev collocation methods and an explicit Runge–Kutta scheme of order four for solving the coupled version of NLSE was discussed. In [9], a fourth order explicit scheme for solving the coupled version of NLSE was considered, while the Galerkin method for the same problem was considered in [10]. Recently, Xu and Zhang [11] introduced and developed alternating direction implicit methods for numerical treatment of two-dimensional cubic NLSE with studying the stability analysis. Meanwhile, the Pade approximation was considered in [12] for the solution of a higher-order NLSE with fourth order dispersion. More
*
Corresponding author. E-mail addresses: [email protected] (E.H. Doha), [email protected] (A.H. Bhrawy), [email protected] (M.A. Abdelkawy), [email protected] (R.A. Van Gorder). 0021-9991/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2014.01.003
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
245
recently, Duan and Rong [13] discussed and developed a multi-quadric quasi-interpolation scheme for solving NLSE. Many other methods have been studied in the literature; see, for example, [14–16]. Spectral methods (see [17–21]) are often an efficient and highly accurate scheme when compared with local methods. There are three primary types of spectral methods based on the choice of test functions, namely, the Galerkin, tau and collocation methods. In these methods, one reduces the original PDEs into a system of ODEs (often in the time variable) by assuming a functional dependence on the spatial part of the equation. In such algorithms, one then solves for the time evolution of the solutions. It is within the realm of nonlinear and variable coefficient problems that the advantages of these collocation methods become apparent. Indeed, the projection operators involved in the Galerkin and tau methods lead to complications in deriving the required system of ODEs, from which the final solution is obtained. In both these methods, we require that the projection of the residual onto some polynomial space be zero. In contrast, the collocation method obtains the system of ODEs by requiring that the residual vanishes at a certain set of grid points, typically some set of Gauss-type points. As a result, while the Galerkin and tau methods may become untenable for strongly nonlinear problems, the collocation method has been shown in many cases to deal with linear, variable coefficient, and nonlinear problems with equal ease [22–27]. Along these lines, a number of related methods have been proposed for solving related equations [28–34]. The Jacobi polynomials are a class of orthogonal polynomials, which satisfy an orthogonality condition on the interval [−1, 1] with respect to the weight function (1 − x)θ (l + x)ϑ . Several orthogonal polynomials can obtained as special case of the Jacobi polynomials, such as the Legendre, Chebyshev, and Gegenbauer [35]. In recent years, the use of Jacobi polynomials for solving differential equations through collocation techniques has gained increasing popularity due to the ability to obtain the solution in terms of the general Jacobi parameters θ and ϑ (see [36,37]). The main concern of this paper is to extend the J-GL-C method to the solution of various types of NLSE. In particular, we shall consider the general initial-boundary value problem
2
i ∂t ψ( y , t ) + ∂ y y ψ( y , t ) + 2γ ψ( y , t ) ψ( y , t ) − 2δ R ( y , t )ψ( y , t ) = 0,
(1.1)
where
( y , t ) ∈ D × [0, T ],
D = { y: A y B },
with the boundary-initial conditions
ψ( A , t ) = ζ1 (t ), ψ( y , 0) = ξ1 ( y ),
ψ( B , t ) = ζ2 (t ), y ∈ D.
(1.2) (1.3)
Here we have maintained a general potential R ( y , t ) and the cubic nonlinearity. This type of equation is also referred to as a Gross–Pitaevskii equation when we have non-zero R ( y , t ). The nonlinear PDEs of type (1.1) are collocated only for the space variable at ( N − 1) points, for suitable collocation points. We use the ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation, which depends upon the two general parameters θ and ϑ . These equations, together with two-point boundary conditions, constitute system of (2N − 2) ODEs in the time variable. This system can then be solved by the implicit Runge–Kutta method of fourth order. Finally, the accuracy of the proposed method is demonstrated by four test problems, which correspond to physically meaningful cases. To the best of our knowledge, there are no results on the J-GL-C method for solving NLSE. This paper is arranged as follows. In Section 2, we derive the numerical method. First, we introduce some properties of Jacobi polynomials and their associated recursion relations. Then the Gauss–Lobatto collocation technique for a general form of the NLSE obtained using the Jacobi polynomials, and the implicit Runge–Kutta method of fourth order is outlined for the resulting system of ODEs. In Section 3, the proposed method is applied to four specific NLSE. It is here where we demonstrate the usefulness of the method, in terms of efficiency and accuracy. Indeed, for relatively few nodes, we obtain rather good results. In fact, for all examples considered, we obtain errors of order 10−6 or lower using only nine nodes. Finally, some concluding remarks are given in Section 4. 2. Derivation of the numerical method In the present section, we derive the numerical method. First, we remind the reader of some useful properties of Jacobi polynomials. Then, we develop the J-GL-C method to numerically solve the NLSE (1.1). We conclude this section with the implicit Runge–Kutta method used to solve the time-dependent ODEs. 2.1. Some properties of the Jacobi polynomials We collected in this section some basic knowledge of Jacobi polynomials that are most relevant to spectral approximations [35]. It includes Legendre and Chebyshev polynomials as two special cases, so it is worthwhile work with general Jacobi polynomials. A basic property of the Jacobi polynomials is that they are the eigenfunctions to a singular Sturm– Liouville problem:
246
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
1 − x2 φ (x) + ϑ − θ + (θ + ϑ + 2)x φ (x) + n(n + θ + ϑ + 1)φ(x) = 0.
(2.1)
The following recurrence relation generate the Jacobi polynomials:
(θ,ϑ)
(θ,ϑ)
J k+1 (x) = ak (θ,ϑ)
J0
(θ,ϑ) (θ,ϑ)
x − bk (θ,ϑ)
(x) = 1,
J1
Jk
(θ,ϑ) (θ,ϑ) J k−1 (x),
(x) − ck
1
1
2
2
k 1,
(x) = (θ + ϑ + 2)x + (θ − ϑ),
where
(2k + θ + ϑ + 1)(2k + θ + ϑ + 2) , 2(k + 1)(k + θ + ϑ + 1) (ϑ 2 − θ 2 )(2k + θ + ϑ + 1) (θ,ϑ) bk = , 2(k + 1)(k + θ + ϑ + 1)(2k + θ + ϑ) (k + θ)(k + ϑ)(2k + θ + ϑ + 2) (θ,ϑ) ck = . (k + 1)(k + θ + ϑ + 1)(2k + θ + ϑ) (θ,ϑ)
ak
=
The Jacobi polynomials are satisfying the following identities: (θ,ϑ)
Jk
(θ,ϑ)
(−x) = (−1)k J k
(θ,ϑ)
(x),
Jk
(−1) =
(−1)k Γ (k + ϑ + 1) . k!Γ (ϑ + 1) (θ,ϑ)
Moreover, the q derivative of Jacobi polynomials of degree k ( J k (θ,ϑ)
D (q) J k
(x) =
(2.2)
(x)), can be obtained from:
Γ ( j + θ + ϑ + q + 1) (θ +q,ϑ+q) J (x). 2q Γ ( j + θ + ϑ + 1) k−q
(2.3)
Let w (θ,ϑ) (x) = (1 − x)θ (1 + x)ϑ , then we define the weighted space L 2w (θ,ϑ) as usual. The inner product and the norm of
L 2w (θ,ϑ)
with respect to the weight function are defined as:
1 (u , v ) w (θ,ϑ) =
u (x) v (x) w (θ,ϑ) (x) dx,
1
u w (θ,ϑ) = (u , u ) w2 (θ,ϑ) .
(2.4)
−1
The set of Jacobi polynomials forms a complete L 2w (θ,ϑ) -orthogonal system, and
(θ,ϑ) J k
w (θ,ϑ)
= hk =
2θ +ϑ+1 Γ (k + θ + 1)Γ (k + ϑ + 1)
(2k + θ + ϑ + 1)Γ (k + 1)Γ (k + θ + ϑ + 1)
(2.5)
.
2.2. Jacobi spectral collocation method The main objective of this section is to develop the J-GL-C method to numerically solve the NLSE (1.1). In order to begin the numerical solution, we put the complex functions ψ in real and imaginary parts as:
ψ( y , t ) = u ( y , t ) + i v ( y , t ),
(2.6)
where u and v are real functions, from Eqs. (1.1) and (2.6), we get
i ∂t u ( y , t ) + ∂ y y v ( y , t ) + 2γ u 2 ( y , t ) + v 2 ( y , t ) v ( y , t ) − 2δ R ( y , t ) v ( y , t )
− ∂t u ( y , t ) − ∂ y y v ( y , t ) − 2γ u 2 ( y , t ) + v 2 ( y , t ) u ( y , t ) + 2δ R ( y , t ) u ( y , t ) = 0,
(2.7)
which can be separated as coupled time-dependent nonlinear PDEs,
∂t u ( y , t ) + ∂ y y v ( y , t ) + 2γ u 2 ( y , t ) + v 2 ( y , t ) v ( y , t ) − 2δ R ( y , t ) v ( y , t ) = 0, ∂t v ( y , t ) − ∂ y y u ( y , t ) − 2γ u 2 ( y , t ) + v 2 ( y , t ) u ( y , t ) + 2δ R ( y , t ) u ( y , t ) = 0,
(2.8)
with the boundary-initial conditions
u ( A , t ) = g 1 (t ),
u ( B , t ) = g 2 (t ),
v ( A , t ) = g 3 (t ),
v ( B , t ) = g 4 (t ),
u ( y , 0) = f 1 ( y ),
v ( y , 0) = f 2 ( y ),
y ∈ D.
(2.9)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
247
Now, suppose the change of variables
x=
2 B−A
y+
A+B
u ( y , t ) = w (x, t ),
,
A−B
v ( y , t ) = z(x, t ),
which will be used to transform problem (1.1)–(1.3) into another one in the classical interval, [−1, 1], for the space variable, to directly implement collocation method based on Jacobi family defined on [−1, 1],
2
2
∂t w (x, t ) +
B−A
2
∂t z(x, t ) −
2
B−A
∂xx z(x, t ) + 2γ w 2 (x, t ) + z2 (x, t ) z(x, t ) − 2δr (x, t ) z(x, t ) = 0,
∂xx w (x, t ) − 2γ w 2 (x, t ) + z2 (x, t ) w (x, t ) + 2δr (x, t ) w (x, t ) = 0,
(2.10)
where
(x, t ) ∈ D ∗ × [0, T ],
D ∗ = {x: −1 x 1},
with the boundary-initial conditions
w (−1, t ) = g 5 (t ),
w (1, t ) = g 6 (t ),
z(−1, t ) = g 7 (t ),
z(1, t ) = g 8 (t ),
w (x, 0) = f 3 (x),
z(x, 0) = f 4 (x),
x ∈ D∗.
(2.11)
The main advantage of the use of the node points of Jacobi–Gauss–Lobatto quadrature is that the distribution of these points in [−1, 1] can be controlled by the two Jacobi parameters, θ and ϑ . Now, we outline the main step of implementing J-GL-C method for reducing nonlinear coupled system (2.10)–(2.11) to system of ODEs. The expansion of functions w (x, t ) (θ,ϑ) and z(x, t ), using the Jacobi polynomial, J j (x), takes the form
w (x, t ) =
N
(θ,ϑ)
a j (t ) J j
(x),
(2.12)
j =0
z(x, t ) =
N
(θ,ϑ)
b j (t ) J j
(x),
(2.13)
j =0
and in virtue of (2.4)–(2.5), we deduce that
a j (t ) =
b j (t ) =
1 hj 1 hj
1
(θ,ϑ)
w (x, t ) w (θ,ϑ) (x) J j
(x) dx,
−1
1
(θ,ϑ)
z(x, t ) w (θ,ϑ) (x) J j
(x) dx.
(2.14)
−1
We can circumvent the need for evaluating the previous integrals, by using Jacobi–Gauss–Lobatto quadrature formula to approximate the integrals. For any φ ∈ S 2N +1 (−1, 1),
1 w (θ,ϑ) (x)φ(x) dx =
N
(θ,ϑ) N(θ,ϑ) , j φ xN , j ,
(2.15)
j =0
−1
(θ,ϑ)
(θ,ϑ)
where S N (−1, 1) is the set of polynomials of degree less than or equal to N, x N , j (0 j N) and N , j (0 j N) are the nodes and the corresponding Christoffel numbers of the Jacobi–Gauss–Lobatto quadrature formula on the interval (−1, 1), respectively. Due to the orthogonality property of Jacobi polynomials (2.4), the expansion coefficients a j (t ), b j (t ) in terms of the solution at the Jacobi–Gauss–Lobatto grid points can be approximated by
a j (t ) =
N 1 (θ,ϑ) (θ,ϑ) (θ,ϑ) (α ,β) Jj x N ,i N ,i w x N ,i , t , hj i =0
b j (t ) =
N 1 (θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i N ,i z x N ,i , t . hj i =0
(2.16)
248
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Accordingly, (2.12)–(2.13) can be given by
w (x, t ) =
N N 1 i =0
z(x, t ) =
j =0
hj
N N 1 i =0
j =0
(θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i Jj (x) N ,i
hj
(θ,ϑ) (θ,ϑ) (θ,ϑ) (θ,ϑ) Jj x N ,i Jj (x) N ,i
(θ,ϑ)
w x N ,i , t ,
(θ,ϑ)
z x N ,i , t .
(2.17)
Furthermore, if we differentiate (2.17) once using (2.3), and evaluate it at all J-GL-C points, it is easy to compute the first spatial partial derivative in terms of the values at theses collocation points as
(θ,ϑ)
w x x N ,n , t =
N
(θ,ϑ)
A ni w x N ,i , t ,
i =0
(θ,ϑ)
zx x N ,n , t =
N
(θ,ϑ)
n = 0, 1 , . . . , N ,
B ni z x N ,i , t ,
(2.18)
i =0
where
A ni =
N j+θ +ϑ +1
2h j
j =0
B ni =
N j+θ +ϑ +1
2h j
j =0
(θ,ϑ) (θ,ϑ) (θ +1,ϑ+1) (θ,ϑ) (θ,ϑ) x N ,i J j −1 x N ,n N ,i ,
Jj
(θ,ϑ) (θ,ϑ) (θ +1,ϑ+1) (θ,ϑ) (θ,ϑ) x N ,i J j −1 x N ,n N ,i .
Jj
(2.19)
Similar steps can be applied to the second spatial partial derivative to get
(θ,ϑ)
w xx x N ,n , t =
N
(θ,ϑ)
C ni w x N ,i , t ,
i =0
(θ,ϑ)
zxx x N ,n , t
=
N
(θ,ϑ)
D ni z x N ,i , t ,
n = 0, 1 , . . . , N ,
(2.20)
i =0
where
C ni =
N ( j + θ + ϑ + 2)( j + θ + ϑ + 1)
4h j
j =0
D ni =
N ( j + θ + ϑ + 2)( j + θ + ϑ + 1)
4h j
j =0
(θ,ϑ) (θ,ϑ) (θ +2,ϑ+2) (θ,ϑ) (θ,ϑ) x N ,i J j −2 x N ,n N ,i ,
Jj
(θ,ϑ) (θ,ϑ) (θ +2,ϑ+2) (θ,ϑ) (θ,ϑ) x N ,i J j −2 x N ,n N ,i .
Jj
(2.21)
For the spectral J-GL-C method, the residual of (1.1) is set to zero at the N − 1 Jacobi–Gauss–Lobatto points. Also, the boundary conditions (1.2) can be imposed at the two collocation points −1 and 1 in the expansion of functions. Therefore, adopting (2.8)–(2.11), we may write (1.1)–(1.2) in the form
˙ n (t ) = − w
where
2 N
B−A
z˙ n (t ) =
2
2
(θ,ϑ)
i =0
2 N
B−A
D ni zi (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
C ni w i (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
n = 1, 2, . . . , N − 1,
i =0
w k (t ) = w x N ,k , t ,
(θ,ϑ)
zk (t ) = z x N ,k , t ,
k = 1, . . . , N − 1, n = 1, . . . , N − 1.
This provides two ( N − 1) system of first order ODEs in the expansion coefficients a j (t ), namely
(2.22)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
˙ n (t ) = − w z˙ n (t ) =
2 N −1
2 B−A
B−A
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
i =1
2 N −1
2
249
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
(2.23)
i =1
where
dn (t ) = D n0 g 7 (t ) + D nN g 8 (t ), dn (t ) = C n0 g 5 (t ) + C nN g 6 (t ). We get the following two system of nonlinear ODEs:
˙ n (t ) = − w z˙ n (t ) =
2 N −1
2 B−A
B−A
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
i =1
2 N −1
2
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ),
i =1
n = 1, 2, . . . , N − 1,
(2.24)
subject to the initial values
(θ,ϑ)
w n (0) = f 3 x N ,n
(θ,ϑ) f 4 x N ,n ,
zn (0) =
, n = 1, . . . , N − 1.
(2.25)
Finally, (2.24)–(2.25) can be rewritten into two matrix form of N − 1 first order ODEs with their vectors of initial values:
˙ (t ) = F t , w (t ), z(t ) , w z˙ (t ) = G t , w (t ), z(t ) , w(0) = f3 , z(0) = f4 ,
(2.26)
where
T
˙ (t ) = w ˙ 1 (t ), w ˙ 2 (t ), . . . , w ˙ N −1 (t ) w
, T z˙ (t ) = z˙ 1 (t ), z˙ 2 (t ), . . . , z˙ N −1 (t ) , (θ,ϑ) T (θ,ϑ) (θ,ϑ) f3 = f 3 x N ,1 , f 3 x N ,2 , . . . , f 3 x N , N −1 , (θ,ϑ) T (θ,ϑ) (θ,ϑ) f4 = f 4 x N ,1 , f 4 x N ,2 , . . . , f 4 x N , N −1 , T F t , w (t ), z(t ) = F 1 t , w (t ), z(t ) , F 1 t , w (t ), z(t ) , . . . , F N −1 t , w (t ), z(t ) ,
and
T
G t , w (t ), z(t ) = G 1 t , w (t ), z(t ) , G 1 t , w (t ), z(t ) , . . . , G N −1 t , w (t ), z(t )
.
Here,
F n t , w (t ), z(t ) = −
G n t , w (t ), z(t ) =
2 N −1
2 B−A 2
B−A
i =1
2 N −1
D ni zi (t ) + dn (t ) − 2γ w n2 (t ) + zn2 (t ) zn (t ) + 2δrn (t ) zn (t ),
C ni w i (t ) + dn (t ) + 2γ w n2 (t ) + zn2 (t ) w n (t ) − 2δrn (t ) w n (t ).
i =1
The two N − 1 dimensional system of ODEs (2.26) can be combined to be written in a 2N − 2 system of ODEs and then can be solved by using an implicit Runge–Kutta formula of fourth order.
250
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
2.3. Implicit Runge–Kutta (IRK) method The implicit Runge–Kutta (IRK) method is one of the suitable methods for solving a system of ODEs of first order like we have obtained above. Since the IRK methods have excellent stability properties which reduce computational costs. Consider
˙ (t ) = H t , u (t ) u u(0) = f.
(2.27)
Then, the iterative solution is
Ui , j = ui + t
m
a j ,l H (t i + cl t , Ui ,l ),
j = 1, . . . , m ,
l =1
ui +1 = ui + t
m
b j H (t i + c j t , Ui , j ),
∀j
(2.28)
j =1
here t is the step size, t i = t 0 + i t and ui is approximations to the exact solution u(t i ). Moreover, ai j , b j and c j are real coefficients, i , j = 1, 2, . . . , m, which can be represented by means of the Butcher array,
c 1 a11 · · · a1s
.. .
.. .
c s as1 b1
..
. ··· ···
.. . ass bs
3. Numerical examples This section considers four numerical examples to show the accuracy and applicability of the proposed method in the present paper. Comparison of the results obtained by various choices of Jacobi parameters θ and ϑ reveals that the present method is very effective and convenient for various choices of θ and ϑ . We focus on some well-known equations from mathematical physics in order to demonstrate the method. First we consider an attractive NLSE [38,39], and for our next example we consider a repulsive NLSE [38,39]. Third, we consider the Gross–Pitaevski equation [40,41], which describes the ground state of a quantum system of identical bosons. We assume that the potential is periodic in space, consistent with a Bose–Einstein condensate. Finally, we consider a NLSE with a linear instability added. In each of these examples, we shall highlight the accuracy and efficiency of the method. 3.1. Numerical solution for an attractive NLSE We first consider an attractive NLSE [39] in the form:
i ∂t ψ + ∂xx ψ + 2|ψ|2 ψ = 0,
( y, t ) ∈ D × R
(3.1)
with the boundary-initial conditions
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e i (t + A ) ,
ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e i (t + B ) ,
ψ(x, 0) = u (x, 0) + i v (x, 0) = e i (x) .
(3.2)
In this example, we consider D = [ A , B ] = [−1, 1]. Furthermore, we shall approximate the time domain by [0, T ] for numerical purposes. The difference between the measured or inferred value of a numerical solution and its actual value (absolute error) is given by (for real part, imaginary part, and complex modulus)
E 1 (x, t ) = u (x, t ) − u˜ (x, t ), E 2 (x, t ) = v (x, t ) − v˜ (x, t ),
˜ x, t ), E 3 (x, t ) = ψ(x, t ) − ψ(
(3.3)
˜ x, t ) = u˜ (x, t ) + i v˜ (x, t ) are the exact solution and the approximate solution at the where ψ(x, t ) = u (x, t ) + i v (x, t ) and ψ( point (x, t ), respectively. It may be shown that the error norms are related. Indeed,
E 3 (x, t ) = u (x, t ) − u˜ (x, t ) + i v (x, t ) − v˜ (x, t ) =
E 1 (x, t )2 + E 2 (x, t )2 .
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
251
Table 1 Maximum absolute errors with various choices of N, θ and ϑ for Example 1. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
3.81 × 10−1 3.62 × 10−5 4.62 × 10−8
5.73 × 10−1 15.70 × 10−5 22.59 × 10−8
6.21 × 10−1 16.04 × 10−5 22.74 × 10−8
− 12
− 12
2.30 × 10−1 5.75 × 10−5 4.93 × 10−8
4.83 × 10−1 3.65 × 10−5 8.30 × 10−8
4.83 × 10−2 5.75 × 10−5 9.21 × 10−8
2 6 10
1 2
1 2
4.40 × 10−1 5.85 × 10−5 6.41 × 10−8
6.38 × 10−1 22.07 × 10−5 33.44 × 10−8
6.67 × 10−1 22.12 × 10−5 33.64 × 10−8
1 2
0
4.13 × 10−1 7.66 × 10−5 4.86 × 10−8
5.37 × 10−1 21.36 × 10−5 25.47 × 10−8
6.78 × 10−2 22.17 × 10−5 25.58 × 10−8
Fig. 1. The absolute error (a) E 1 (x, t ) and (b) E 2 (x, t ) of problem (3.1) where θ = 0, ϑ = − 12 and N = 10. Note that while errors gradually increase with time, they remain less than 10−6 over the plotted domain.
Moreover, the maximum absolute error is given by
M 1 = Max E 1 (x, t ): ∀(x, t ) ∈ D × [0, T ] , M 2 = Max E 2 (x, t ): ∀(x, t ) ∈ D × [0, T ] , M 3 = Max E 3 (x, t ): ∀(x, t ) ∈ D × [0, T ] . (3.4) Clearly, M 3 M 12 + M 22 . The maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) related to the problem (3.1)–(3.2) are presented in Table 1 using the J-GL-C method with various choices of N, θ and ϑ . The errors are calculated through a comparison with the exact solution
ψ(x, t ) = e i (x+t ) .
(3.5)
We see in this table that the results are very accurate for even small choices of the number of nodes, N. In Fig. 1, we see that the absolute errors E 1 (x, t ) and E 2 (x, t ) are very small, despite the relatively small number of grid points used. This numerical experiment demonstrates the utility of the method.
252
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Table 2 Maximum absolute errors with various choices of N, θ and ϑ for Example 2. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
2.43 × 10−1 6.00 × 10−5 2.26 × 10−7
2.69 × 10−1 5.68 × 10−5 1.70 × 10−7
3.23 × 10−1 6.41 × 10−5 2.35 × 10−7
− 12
− 12
2.57 × 10−1 4.19 × 10−5 2.52 × 10−7
6.40 × 10−2 5.81 × 10−5 1.92 × 10−7
2.64 × 10−1 6.25 × 10−5 2.89 × 10−7
2 6 10
0
− 12
2.71 × 10−1 11.11 × 10−5 2.19 × 10−7
1.98 × 10−1 6.60 × 10−5 1.79 × 10−7
3.35 × 10−1 12.91 × 10−5 2.28 × 10−7
− 12
0
2.12 × 10−1 5.40 × 10−5 2.42 × 10−7
1.33 × 10−1 6.64 × 10−5 1.93 × 10−7
2.36 × 10−1 6.98 × 10−5 2.43 × 10−7
In the case where an exact solution is not known, one follows a similar method. The difference is that one calculates the residual errors only, since comparison to an exact solution is not possible. 3.2. Numerical solution for a repulsive NLSE Consider the repulsive NLSE [39]
i ∂t ψ + ∂xx ψ − 2γ |ψ|2 ψ = 0,
(3.6)
with the boundary-initial conditions
ψ(−1, t ) = u (−1, t ) + i v (−1, t ) = e i (−1−3t ) , ψ(1, t ) = u (1, t ) + i v (1, t ) = e i (1−3t ) , ψ(x, 0) = u (x, 0) + i v (x, 0) = e i x .
(3.7)
Note that the exact solution is given by
ψ(x, t ) = e i (x−3t ) .
(3.8)
Table 2 lists the maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) of problem (3.6)–(3.7) using the J-GL-C method for with various choices of N, θ and ϑ . The numerical results presented in this table show that the results are very accurate for small value of N. Fig. 2 demonstrates that the absolute errors E 1 (x, t ) and E 2 (x, t ) are very small for even the small number of grid points taken. 3.3. Nonlinear Gross–Pitaevskii equation with space-periodic potential Consider the complex nonlinear Gross–Pitaevskii (CNGP) equation with a space-periodic potential [39–41] given by
i ∂t ψ +
1 2
∂xx ψ − |ψ|2 ψ − cos(x)ψ = 0,
(3.9)
with the boundary-initial conditions
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e
−3it
sin( A ),
2
−3it 2
sin( B ),
ψ(x, 0) = u (x, 0) + i v (x, 0) = sin(x).
(3.10)
The exact solution of CNGP equation is given by
ψ(x, t ) = e
−3it 2
sin(x).
(3.11)
To demonstrate the accuracy of the proposed method, in Table 3 we introduce maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) of CNGP equation, using the J-GL-C method with various choices of θ and ϑ and N. We see that the approximate solutions u˜ (x, t ) ( v˜ (x, t )) and the exact solution u (x, t ) (v (x, t )) are coincided for different values of t and x. 3.4. An attractive NLSE with linear instability We consider an attractive NLSE [39] with an additional linear term
i ∂t ψ + ∂xx ψ + 2|ψ|2 ψ − 2γ ψ = 0, with the boundary-initial conditions
(3.12)
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
253
Fig. 2. The absolute error (a) E 1 (x, t ) and (b) E 2 (x, t ) of problem (3.6) where θ = ϑ = − 12 and N = 10. Table 3 Maximum absolute errors with various choices of N, θ and ϑ for Example 3. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
4.29 × 10−2 3.84 × 10−6 8.32 × 10−8
4.00 × 10−2 3.28 × 10−6 7.98 × 10−8
5.87 × 10−2 3.97 × 10−6 10.49 × 10−8
− 12
− 12
5.46 × 10−2 4.90 × 10−6 3.05 × 10−7
4.00 × 10−2 5.56 × 10−6 1.42 × 10−7
5.87 × 10−2 6.84 × 10−6 3.05 × 10−7
2 6 10
1 2
0
5.25 × 10−2 3.67 × 10−6 2.11 × 10−7
4.66 × 10−2 9.26 × 10−6 1.69 × 10−7
7.02 × 10−2 9.26 × 10−6 2.43 × 10−7
1 2
1 2
4.29 × 10−2 3.97 × 10−6 1.01 × 10−7
4.00 × 10−2 5.30 × 10−6 7.13 × 10−8
5.87 × 10−2 5.40 × 10−6 1.24 × 10−7
ψ( A , t ) = u ( A , t ) + i v ( A , t ) = e i ( A +(1−2γ )t ) , ψ( B , t ) = u ( B , t ) + i v ( B , t ) = e i ( B +(1−2γ )t ) ψ(x, 0) = u (x, 0) + i v (x, 0) = e ix .
(3.13)
Maximum absolute errors of u (x, t ), v (x, t ) and ψ(x, t ) related to (3.12)–(3.13) are introduced in Table 4 using J-GL-C method for with various choices of N, θ and ϑ , where γ = 1 compared with the exact solution
ψ(x, t ) = e i (x+(1−2γ )t ) .
(3.14)
Again and from Table 4, we clearly observe that the global scheme (J-GL-C) is accurate for even small values of N. 4. Conclusion We have proposed an accurate numerical algorithm based on the J-GL-C method for solving complex NLSE. Our procedure is implemented in two successive steps. First, the J-GL-C method is employed for approximating the functional dependence on the spatial variable, using ( N − 1) nodes of the Jacobi–Gauss–Lobatto interpolation which depends upon two general Jacobi parameters. The resulting equations together with the two-point boundary conditions give us a system of 2( N − 1)
254
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
Table 4 Maximum absolute errors with various choices of N, θ and ϑ for Example 4. N
θ
ϑ
M1
M2
M3
θ
ϑ
M1
M2
M3
2 6 10
0
0
6.12 × 10−1 1.53 × 10−4 1.82 × 10−7
2.89 × 10−1 6.80 × 10−4 1.40 × 10−7
6.21 × 10−1 1.60 × 10−4 2.15 × 10−7
− 12
− 12
4.01 × 10−1 4.19 × 10−4 2.27 × 10−7
2.69 × 10−1 5.69 × 10−4 1.04 × 10−7
4.83 × 10−1 5.76 × 10−4 2.27 × 10−7
2 6 10
0
1 2
7.43 × 10−1 2.09 × 10−4 2.07 × 10−7
5.85 × 10−2 1.26 × 10−4 1.59 × 10−7
7.43 × 10−1 2.09 × 10−4 2.35 × 10−7
1 2
1 2
6.59 × 10−1 2.18 × 10−4 2.64 × 10−7
1.60 × 10−1 7.10 × 10−4 1.96 × 10−8
6.67 × 10−1 2.21 × 10−4 3.02 × 10−7
first-order ODEs in the time variable. In the second part of the solution method, the implicit Runge–Kutta method of fourth order is applied to solve this system of 2( N − 1) ODEs in time. Note that it is possible to use other implicit Runge–Kutta formulas, say an implicit Runge–Kutta method of order greater than four, or even a second-order method. The main advantage of the use of the node points of Jacobi–Gauss–Lobatto quadrature as collocation points is that the distribution of these points in a specified interval can be tuned by the two Jacobi parameters, θ and ϑ . Through appropriate selection of these parameters, we may recover some of the well-known spectral collocation approximations. Namely, the Legendre and Chebyshev pseudo-spectral approximations can be derived as special cases from our algorithm if we suitably choose the corresponding spacial cases of Jacobi parameters θ and ϑ . From the numerical results given in Section 3, it has been concluded that the obtained results are excellent in terms of accuracy for all tested problems. For each of the examples given, we have been able to construct solutions with absolute errors between 10−7 and 10−6 , and greater accuracy may be obtained if we add more nodes. Note that the initial data considered has always been of modulus unity, so the relative errors should be of the same order of magnitude as the absolute errors. Concerning the efficiency of the method, note that we have obtained approximate numerical solutions of a desirable accuracy by taking only N = 10 (i.e., 9 nodes). Considering the size of the domains considered, this is approximately one node per 0.2 units of space. Considering the rather large average spatial gap, the results are rather good. Of course, adding more nodes will increase the accuracy further, and one should be able to employ this method in order to determine solutions of a desired accuracy. With this paper we have outlined the application of a J-GL-C method for solving complex NLSE. In principle, this algorithm may be extended to related problems, such as to coupled complex nonlinear Schrödinger equations. One might also consider higher order dispersive equations. We should note that, as a numerical method, we are restricted to solving problems over a finite domain. Hence, this method is particularly well suited for boundary value problems with finite spatial domains. Acknowledgement R.A.V. supported in part by NSF grant # 1144246. References [1] S. Islam, S. Haq, M. Uddin, A mesh free interpolation method for the numerical solution of the coupled nonlinear partial differential equations, Eng. Anal. Bound. Elem. 33 (2009) 399–409. [2] A. Mohebbi, Z. Asgari, M. Dehghan, Numerical solution of nonlinear Jaulent–Miodek and Whitham–Broer–Kaup equations, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 4602–4610. [3] R.C. Mittal, R.K. Jain, Numerical solutions of nonlinear Burgers’ equation with modified cubic B-splines collocation method, Appl. Math. Comput. 218 (2012) 7839–7855. [4] 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 Schrodinger equation arising in quantum mechanics, Eng. Anal. Bound. Elem. 37 (2013) 475–485. [5] Y. Gurefe, E. Misirli, A. Sonmezoglu, M. Ekici, Extended trial equation method to generalized nonlinear partial differential equations, Appl. Math. Comput. 219 (2013) 5253–5260. [6] H. Wang, Numerical studies on the split step finite difference method for the nonlinear Schrödinger equations, Appl. Math. Comput. 170 (2005) 17–35. [7] M. Dehghan, A. Taleei, Numerical solution of nonlinear Schrödinger equation by using time–space pseudo-spectral method, Numer. Methods Partial Differ. Equ. 26 (2010) 979–992. [8] M. Dehghan, A. Taleei, A Chebyshev pseudo-spectral multi-domain method for the soliton solution of coupled nonlinear Schrodinger equations, Comput. Phys. Commun. 182 (2011) 2519–2529. [9] M.S. Ismail, A fourth order explicit scheme for the coupled nonlinear Schrödinger equation, Appl. Math. Comput. 196 (2008) 273–284. [10] M.S. Ismail, Numerical solution of coupled nonlinear Schrödinger equation by Galerkin method, Math. Comput. Simul. 78 (2008) 532–547. [11] Y. Xu, L. Zhang, Alternating direction implicit method for solving two-dimensional cubic nonlinear Schrödinger equation, Comput. Phys. Commun. 183 (2012) 1082–1093. [12] M. Smadi, D. Bahloul, A compact split step Pade scheme for higher-order nonlinear Schrödinger equation (HNLS) with power law nonlinearity and fourth order dispersion, Comput. Phys. Commun. 182 (2011) 366–371. [13] Y. Duan, F. Rong, A numerical scheme for nonlinear Schrödinger equation by MQ quasi-interpolation, Eng. Anal. Bound. Elem. 37 (2013) 89–94. [14] H. Sakaguchi, T. Higashiuchi, Two-dimensional dark soliton in the nonlinear Schrödinger equation, Phys. Lett. A 359 (2006) 647–651. [15] Z. Kalogiratou, Th. Monovasilis, T.E. Simos, Numerical solution of the two-dimensional time independent Schrödinger equation with Numerov-type methods, J. Math. Chem. 37 (3) (2005) 271–279.
E.H. Doha et al. / Journal of Computational Physics 261 (2014) 244–255
255
[16] M. Javidi, A. Golbabai, Numerical studies on nonlinear Schrödinger equations by spectral collocation method with preconditioning, J. Math. Anal. Appl. 333 (2007) 1119–1127. [17] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods: Fundamentals in Single Domains, Springer-Verlag, New York, 2006. [18] E.H. Doha, A.H. Bhrawy, An efficient direct solver for multidimensional elliptic Robin boundary value problems using a Legendre spectral-Galerkin method, Comput. Math. Appl. 64 (2012) 558–571. [19] S. Karimi Vanani, F. Soleymani, Tau approximate solution of weakly singular Volterra integral equations, Math. Comput. Model. 57 (3–4) (2013) 494–502. [20] W. Shao, X. Wu, Chebyshev tau meshless method based on the highest derivative for fourth order equations, Appl. Math. Model. 37 (3) (1 February 2013) 1413–1430. [21] E.H. Doha, W.M. Abd-Elhameed, M.A. Bassuony, New algorithms for solving high even-order differential equations using third and fourth Chebyshev– Galerkin methods, J. Comput. Phys. 236 (2013) 563–579. [22] C.I. Gheorghiu, Spectral Methods for Differential Problems, T. Popoviciu Institute of Numerical Analysis, Cluj-Napoca, Romaina, 2007. [23] C.I. Gheorghiu, M.E. Hochstenbach, B. Plestenjak, J. Rommes, Spectral collocation solutions to multiparameter Mathieu’s system, Appl. Math. Comput. 218 (24) (2012) 11990–12000. [24] A.H. Bhrawy, M.A. Alghamdi, A shifted Jacobi–Gauss–Lobatto collocation method for solving nonlinear factional Langevin equation involving two fractional orders in different intervals, Bound. Value Probl. 2012 (2012) 62. [25] K. Zhang, J. Li, H. Song, Collocation methods for nonlinear convolution Volterra integral equations with multiple proportional delays, Appl. Math. Comput. 218 (2012) 10848–10860. [26] A. Saadatmandi, M. Dehghan, The use of Sinc-collocation method for solving multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 593–601. [27] A.H. Bhrawy, A.S. Alofi, A Jacobi–Gauss collocation method for solving nonlinear Lane–Emden type equations, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 62–70. [28] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numer. Methods Partial Differ. Equ. 28 (2012) 74–93. [29] M.R. Eslahchi, M. Dehghan, S. Ahmadi-Asl, The general Jacobi matrix method for solving some nonlinear ordinary differential equations, Appl. Math. Model. 36 (2012) 3387–3398. [30] M.R. Eslahchi, M. Dehghan, S. Amani, The third and fourth kinds Chebyshev polynomials and best uniform approximation, Math. Comput. Model. 55 (2012) 1746–1762. [31] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Model. 53 (2011) 1865–1877. [32] M. Dehghan, R. Jazlanian, A fourth-order central Runge–Kutta scheme for hyperbolic conservation laws, Numer. Methods Partial Differ. Equ. 26 (2010) 1675–1692. [33] A. Mohebbi, M. Dehghan, The use of compact boundary value method for the solution of two-dimensional Schrödinger equation, J. Comput. Appl. Math. 225 (2009) 124–134. [34] M. Dehghan, D. Mirzaei, Numerical solution to the unsteady two-dimensional Schrödinger equation using meshless local boundary integral equation method, Int. J. Numer. Methods Eng. 76 (2008) 501–520. ˝ Orthogonal Polynomials, Colloquium Publications, vol. XXIII, American Mathematical Society, ISBN 978-0-8218-1023-1, 1939, MR 0372517. [35] G. Szego, [36] E.H. Doha, A.H. Bhrawy, Efficient spectral-Galerkin algorithms for direct solution of fourth-order differential equations using Jacobi polynomials, Appl. Numer. Math. 58 (2008) 1224–1244. [37] M. El-Kady, Jacobi discrete approximation for solving optimal control problems, J. Korean Math. Soc. 49 (2012) 99–112. [38] M.J. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 1981. [39] K. Mallory, R.A. Van Gorder, Stationary solutions for the 1 + 1 nonlinear Schrödinger equation modeling repulsive Bose–Einstein condensates in small potentials, Phys. Rev. E 88 (2013) 013205. [40] E.P. Gross, Structure of a quantized vortex in boson systems, Nuovo Cimento 20 (1961) 454–457. [41] L.P. Pitaevsk, Vortex lines in an imperfect Bose gas, Sov. Phys. JETP 13 (1961) 451–454.
