- Email: [email protected]
Numerical solution of coupled nonlinear Schrödinger equation by Galerkin method
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 78 (2008) 532–547
Numerical solution of coupled nonlinear Schr¨odinger equation by Galerkin method M.S. Ismail ∗ Department of Mathematics, College of Science, P.O. Box 80203, King Abdulaziz University, Jeddah, Saudi Arabia Received 19 November 2006; received in revised form 15 June 2007; accepted 31 July 2007 Available online 7 August 2007
Abstract The coupled nonlinear Schr¨odinger equation models several interesting physical phenomena presents a model equation for optical fiber with linear birefringence. In this paper we derive a finite element scheme to solve this equation, we test this method for stability and accuracy, many numerical tests have been conducted. The scheme is quite accurate and describe the interaction picture clearly. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Galerkin method; Coupled nonlinear Schr¨odinger equation; Soliton interaction
1. Introduction The coupled nonlinear Schr¨odinger equations (CNLSE) arise in a great variety of physical situations. In fiber communication system, such equations have been shown to govern pulse propagation along orthogonal polarization axes in nonlinear optical fibers and in wavelength-division-multiplexed systems [10,19,24]. These equations also model beam propagation inside crystals or photorefractives as well as water wave interactions. Solitary waves in these equations are often called vector solitons in the literature as they generally contain two components. In all the above physical situations, collision of vector solitons is an important issue. This system has been studied intensively in the past 10 years. It has been shown that, in addition to passing-through collision, vector solitons can also bounce off each other or trap each other. In this paper we will consider the coupled nonlinear Schr¨odinger equations (CNLSE): i
∂Ψ1 1 ∂ 2 Ψ1 + + (|Ψ1 |2 + e|Ψ2 |2 )Ψ1 = 0, ∂t 2 ∂x2
−∞ < x < ∞,
(1)
i
∂Ψ2 1 ∂ 2 Ψ2 + + (e|Ψ1 |2 + |Ψ2 |2 )Ψ2 = 0, ∂t 2 ∂x2
−∞ < x < ∞,
(2)
where Ψ1 and Ψ2 are the wave amplitudes in two polarizations and e is the cross-phase modulation coefficient, with initial conditions: Ψ1 (x, 0) = g1 (x), ∗
Ψ2 (x, 0) = g2 (x)
Tel.: +966 506608949. E-mail address: [email protected].
0378-4754/$32.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.07.003
(3)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
533
and boundary conditions: ∂Ψ1 (x, t) ∂Ψ2 (x, t) = =0 ∂x ∂x
as |x| → ∞.
(4)
For linearly birefringent fibers e = 2/3, and for elliptically birefringent fibers e can take other positive values. When e = 0 this system becomes two decoupled nonlinear Schr¨odinger equations (NLS) and when e = 1, the system is known as Manakov equations. In both cases, the system is integrable. In particular, solitons of one polarization should pass through pulses of the opposite polarization without creating shadows and so the collision is elastic. For other values of e the system is nonintegrable, the collision may be nontrivial and various complex phenomena like reflection, transmission, trapping and creation of new solitary waves can occur. Solitons in optical fibers can be defined as nonlinear pulses that propagate nearly distortion-free for long distances and that undergo elastic collision. Some new methods are used to derive the analytical solution of the CNLS system, Sweilam and Al-Bar [16] used the variational iteration method to derive some approximations to the exact solution of the CNLSE. This method is a powerful method and many researchers are currently used this technique [21,22]. Following the discussion of Wadati et al. [19], we derived the exact solution of the CNLSE which can be given in a compact form as 2 √ 2α v Ψ1 (x, t) = sech( 2α(x − vt)) exp i vx − −α t , (5) 1+e 2 2 √ 2α v Ψ2 (x, t) = ± sech( 2α(x − vt)) exp i vx − −α t , (6) 1+e 2 with the conserved quantities (i) Mass conservations: ∞ I1 = |Ψ1 |2 dx, −∞
I2 =
∞
−∞
|Ψ2 |2 dx.
(ii) Momentum conservation:
∞ 2 ∂Ψj ∂Ψ¯ j I3 = Ψ¯ j i − Ψj dx. ∂x ∂x −∞
(7)
(8)
j=1
(iii) Energy conservation: ⎧ ⎫ 2 ∞ ⎨ 2 2 ⎬ 1 ∂Ψj 1 I4 = − |Ψj |4 − e|Ψ1 |2 |Ψ2 |2 dx. ⎭ 2 ∂x 2 −∞ ⎩ j=1
(9)
j=1
Many numerical methods have been developed for the NLS equation [9,11–14], comparing to few numerical methods for the coupled nonlinear Schr¨odinger equation, most of these methods using finite difference method [5–8,15], recently a multi-symplectic methods are used for solving the CNLS system [1,2,17,18] and a very good results obtained concerning the behavior of the solution and the conserved quantities. Wang [20] presented a numerical solution of the single and coupled nonlinear Schr¨odinger equation using split-step finite difference method. Xu and Shu [23], they solve the NLS and CNLSE using the local discontinuous Galerkin method, they analyzed and test this method, the results are in a very good agreement with the previous existing numerical works. Concerning the finite element method many papers have been published for nonlinear Schr¨odinger equation [3,4]. In this work we are going to solve the CNLSE using one of the most popular finite element method which is called Galerkin method. This method will be analyzed for accuracy and stability. We have found that the scheme is second order in both directions space and time, it is unconditionally stable in the linearized sense. Newton’s method is used to solve the nonlinear system obtained from the full discertization of the coupled nonlinear Schr¨odinger equation. The paper is organized as follows: In Section 2, the derivation of the method and the analysis of the properties of the scheme, like stability, accuracy and the solution of the nonlinear system are given. Numerical tests with various initial
534
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
conditions of the CNLSE like single soliton, interaction of three solitons, interaction of arbitrary initial conditions are reported in Section 3. Conclusions are contained in Section 4. 2. Numerical method We assume that the solution of the system (1) and (2) is negligibly small outside the interval [xL , xR ], and so we consider the coupled nonlinear Schr¨odinger: i
∂Ψ1 1 ∂ 2 Ψ1 + (|Ψ1 |2 + e|Ψ2 |2 )Ψ1 = 0, + ∂t 2 ∂x2
xL < x < xR ,
(10)
i
∂Ψ2 1 ∂ 2 Ψ2 + + (e|Ψ1 |2 + |Ψ2 |2 )Ψ2 = 0, ∂t 2 ∂x2
xL < x < xR ,
(11)
where Ψ1 and Ψ2 are the wave amplitudes in two polarizations and e is the cross-phase modulation coefficient, with initial conditions: Ψ1 (x, 0) = g1 (x),
Ψ2 (x, 0) = g2 (x)
(12)
and boundary conditions: ∂Ψ1 (x, t) ∂Ψ2 (x, t) = = 0, ∂x ∂x
x = xL
and
x = xR .
(13)
For our numerical we decompose the complex functions Ψ1 and Ψ2 in the CNLSE into its real and imaginary parts by writing: Ψ1 (x, t) = u1 + iu2 ,
Ψ2 (x, t) = u3 + iu4
(14)
where (ui , i = 1, . . . , 4) are real functions [4]. On substituting these in Eqs. (1) and (2), the following system is obtained: ∂u1 1 ∂ 2 u2 + z1 u2 = 0, + ∂t 2 ∂x2
∂u2 1 ∂ 2 u1 − z1 u1 = 0, − ∂t 2 ∂x2
∂u3 1 ∂ 2 u4 + z2 u4 = 0, + ∂t 2 ∂x2
∂u4 1 ∂ 2 u3 − − z 2 u3 = 0 ∂t 2 ∂x2
(15)
where z1 = u21 + u22 + e(u23 + u24 ),
z2 = e(u21 + u22 ) + u23 + u24 .
(16)
The above system can be written in a matrix–vector form as ∂u 1 ∂2 u + A 2 + F (u)u = 0, ∂t 2 ∂x where
⎡
u1
⎤
⎢u ⎥ ⎢ 2⎥ u = ⎢ ⎥, ⎣ u3 ⎦ u4
⎡
0
(17)
1
⎢ −1 0 ⎢ A=⎢ ⎣ 0 0 0 0
0
0
⎤
0 0⎥ ⎥ ⎥, 0 1⎦ −1 0
⎡
0
⎢ −z ⎢ 1 F (u) = ⎢ ⎣ 0 0
z1
0
0 0
0 0
0
−z2
0
⎤
0⎥ ⎥ ⎥. z2 ⎦ 0
The finite element method [4,8] for the resulting system (17) can be described as follows. The weak form of (17) is obtained by taking the L2 -inner product with elements of the space H 1 (xL , xR ). After integration by parts, and using the boundary conditions we obtain
∂u 1 ∂u ∂ψ (18) ,ψ − A , + (G(u), ψ) = 0, ∀ψ ∈ H 1 (xL , xR ) ∂t 2 ∂x ∂x
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
535
where G(u) = F (u)u and the L2 -inner product (·, ·) is interpreted pointwise and denoted by xR (f, g) = f (x)g(x) dx. xL
The interval [xL , xR ] is discretized by uniform (N + 1) grid points: xm = xL + mh,
m = 0, 1, . . . , N,
where the grid spacing h is given by h = (xR − xL )/N. Let S h be an N-dimensional subspace of H 1 . We choose the elements of S h to be the piecewise linear hat functions and denoted by {ψm }N m=1 . We assume the approximation of the exact solution u(x, t) is given by U(x, t) =
N
Um (t)ψm (x),
(19)
m=0
where Um , m = 0, 1, . . . , N, are time dependent coefficients which are to be determined. Then the discrete form of (18) is given by N m=0
N
N
m=0
m=0
˙ m (ψm , ψj ) − 1 A Um (ψm (x), ψj ) + (G(U), ψj ) = 0, U 2
j = 0, 1, . . . , N,
(20)
˙ m = ∂Um /∂t, Um = ∂Um /∂x. In order to deal with the nonlinear term in (20), we use product approximation where U [4] by which (G(U), ψm ) =
N
G(Um )(ψm , ψm )
(21)
m=0
Since the trial and the test function are chosen to be the piecewise linear functions: ⎧ ⎪ ⎨ (x − xm−1 )/ h if xm−1 < x ≤ xm , ψm (x) = (xm+1 − x)/ h if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise. We note that the solution (19) is zero outside the interval [xm−1 , xm+1 ], and according to this we can write it as ⎧ ⎪ ⎨ Um−1 (t)ψm−1 (x) + Um (t)ψm (x) if xm−1 < x ≤ xm , U(x) = Um (t)ψm (x) + Um+1 (t)ψm+1 (x) if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise and the nonlinear term as ⎧ ⎪ ⎨ G(Um−1 (t))ψm−1 (x) + G(Um (t))ψm (x) if xm−1 < x ≤ xm , G(U(x, t)) = G(Um (t))ψm (x) + G(Um+1 (t))ψm+1 (x) if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise. By using these definitions, the inner products in (20) can be calculated as follows: xm+1 xm xim+1 ˙ m−1 + 4U ˙m +U ˙ m+1 ), ˙ ψm dx = ˙ ψm dx + ˙ ψm dx = 1 (U U U U 6 xm−1 xm−1 xm xm+1 xm xm+1 1 Ux (ψx )m dx = Ux (ψx )m dx + Ux (ψx )m dx = 2 (Um−1 − 2Um + Um+1 ) h xm−1 xm−1 xm
536
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
and
xm+1
xm−1
G(U)ψm dx =
xm
xm−1
G(U)ψm dx +
xm+1
xm
G(U)ψm dx =
1 (G(Um−1 ) + 4G(Um ) + G(Um+1 )). 6
By making use of these integrals in (20), this will lead us to the first order ordinary differential system:
1 2 ˙ 1 1 2 2 1 + δx Um + 2 Aδx Um + 1 + δx G(Um ) = 0, 6 2h 6
(22)
where δ2x Um = Um−1 − 2Um + Um+1 . On expansion of the central difference operator, and imposing the homogeneous Neumann boundary conditions, we will get ˙ = [S + PD(U)]U, PU
(23)
where each of D(u), P and S are a block-diagonal skew symmetric, a block-tridiagonal symmetric and a blocktridiagonal skew symmetric matrices, respectively, and the structure of these matrices are ⎤ ⎤ ⎡ ⎡ −A A 0 0 ··· 0 2I I 0 0 ··· 0 ⎢ I 4I ⎢ A −2A A 0 ··· 0 ⎥ I 0 ··· 0 ⎥ ⎥ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ .. ⎥ .. ⎥ .. .. .. .. .. .. ⎥ ⎢ ⎢ . . . . . . 0 0 0 . ⎥ . ⎥ −1 ⎢ 1⎢ 0 ⎥, ⎥, ⎢ S = P= ⎢ ⎥ ⎥ .. .. .. 6⎢ 2h2 ⎢ ⎢ ... · · · . . . . . . . . . 0 ⎥ ⎢ ... . . . ··· 0 ⎥ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎦ ⎣ 0 ··· ⎣ 0 ··· 0 A −2A A ⎦ 0 I 4I I 0 ··· ··· 0 A −A 0 ··· ··· 0 I 2I D(u) = −diag[B1 (u1 ), B2 (u2 ), . . . , BN (uN )], where
⎡
0
1
⎢ −1 0 ⎢ A=⎢ ⎣ 0 0 0 0
0
0
⎤
0 0⎥ ⎥ ⎥, 0 1⎦ −1 0
⎡
0 ⎢ −(z ) 1 i ⎢ Bi (ui ) = ⎢ ⎣ 0 0
(z1 )i 0 0
0 0 0
0
−(z2 )i
⎤ 0 0 ⎥ ⎥ ⎥, (z2 )i ⎦
z1 = u21 + u22 + e(u23 + u24 ),
0
z2 = e(u21 + u22 ) + u23 + u24 , and 0 and I are (4 × 4) zero and identity matrices, respectively. 2.1. Full discertization In order to solve the system (23), we first discretize the time interval by the grid points t = nk, n = 0, 1, 2, . . ., n is the approximation of the exact solution u(x , t ). Making use of where k is the time step size, we assume that Um m n the substitution: n+1 n ˙ m = Um − Um U k
and Um =
n Un+1 m + Um , 2
which is equivalent to the implicit midpoint rule. This will lead us to the implicit nonlinear scheme:
n+1 1 2 Um + Unm 1 2 n+1 n 2 n+1 n 1 + δx (Um − Um ) + rAδx (Um + Um ) + k 1 + δx G = 0, 6 6 2
(24)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
where r = k/4h2 . The above scheme can be displayed in a matrix vector form as n+1 n n+1 n U + U + U U P(Un+1 − Un ) = k S + PD . 2 2
537
(25)
The resulting system (25) is a block nonlinear tridiagonal one of 4N nonlinear algebraic equations, in the unknown vector Un+1 , and this can be solved by using Newton’s method, or by a fixed point method. 2.2. Newton’s method The solution of the nonlinear system obtained in (25) can be solved by Newton’s method. To describe Newton’s method for our system we do the following. First we rewrite the nonlinear system (25) as
w + Un w + Un n Q(w) = P(w − U ) − k S + PD = 0, (26) 2 2 where w = Un+1 . Newton’s method for this system can be given as w(l+1) = w(l) − J −1 (w(l) )Q(w(l) ),
l = 0, 1, . . . ,
(27)
where J(w) is the Jacobian matrix contains the partial derivatives of Q(w), it is an (N × N) square matrix whose elements are blocks of size (4 × 4). The initial guess vector is assumed as w(0) =Un . Newton’s method is generally expected to give quadratic convergence, provided that a sufficiently accurate starting value is known and J −1 (w) exists. The iterative scheme in (27) is usually executed in two steps: (i) Solve the linear block tridiagonal system: J(w(l) )y = Q(w(l) ) by Crout’s method for the intermediate vector y. (ii) Update your solution by using: w(l+1) = w(l) − y. The two steps (i) and (ii) will be applied iteratively till the following condition is satisfied: ||w(l+1) − w(l) ||∞ ≤ ∈ where ∈ is a prescribed tolerance. 2.3. Conserved quantities In this section we will study the conservation properties of the CNLSE and the analogues property of the numerical method [6,13,14]. In this work we derive the conservation laws: ∞ ∞ |Ψ1 |2 dx = constant and |Ψ2 |2 dx = constant. (28) −∞
−∞
The other quantities will be calculated using a numerical quadrature. To show the system in Eq. (1) satisfies (28), we multiply Eq. (1) by Ψ¯ 1 (the complex conjugate of Ψ1 ) and its complex conjugate by Ψ1 , and then subtract the later
538
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
from the former to obtain
∂Ψ1 ∂Ψ¯ 1 1 ∂ 2 Ψ1 1 ∂2 Ψ¯ 1 2 2 2 2 ¯ Ψ¯ 1 i − Ψ i =0 + + (|Ψ | + e|Ψ | )Ψ + (|Ψ | + e|Ψ | ) Ψ + 1 1 2 1 1 2 1 ∂t 2 ∂x2 ∂t 2 ∂x2
(29)
and this gives ∂Ψ1 ∂Ψ¯ 1 1 ¯ ∂ 2 Ψ1 ∂2 Ψ¯ 1 i Ψ¯ 1 + Ψ1 Ψ1 2 − Ψ1 2 = 0, + ∂t 2 ∂t ∂x ∂x which is equivalent to ∂ 1 ∂ i (Ψ1 Ψ¯ 1 ) + ∂t 2 ∂x
∂Ψ1 ∂Ψ¯ 1 Ψ¯ 1 − Ψ1 ∂x ∂x
= 0.
(30)
By integrating Eq. (30) with respect to x from −∞ to ∞ one can get ∞ ∂ ∞ ∂Ψ¯ 1 1 ¯ ∂Ψ1 2 i |Ψ1 | dx + Ψ1 =0 − Ψ1 ∂t −∞ 2 ∂x ∂x −∞
(31)
and by using the homogenous Neumann boundary conditions ∂Ψ1 /∂x = 0 as |x| → ∞ we find that ∂ ∞ |Ψ1 |2 dx = 0 ∂t −∞ and hence the first conserved quantity: ∞ I1 = |Ψ1 |2 dx = constant
(32)
is obtained. By adopting the same criteria, the second conserved quantity: ∞ I2 = |Ψ2 |2 dx = constant
(33)
−∞
−∞
can be obtained. The exact values of these quantities using (5) and (6) are 2 √ I1 = I2 = 2α. 1+e
(34)
Conservation of I1 and I2 implies L2 boundedness of the solution, and hence no blow up is expecting. To prove the discrete analogue of the conservation of (28), we apply the derived scheme for the first component of the CNLSE to get P v˙ = [S + PD(u)]v, 0 ˜ = where v = [u1 , u2 ]t , A −1
(35)
1 , and D(u) is a block-diagonal matrix: 0
˜ σ2 A, ˜ . . . , σN A] ˜ D(u) = −diag[σ1 A,
with σj = u21,j + u22,j + e(u23,j + u24,j )
where in this case u consists of the first two components of u a way from the nonlinear terms. The corresponding implicit scheme for (35) is Un+1 + Un Vn+1 + Vn n+1 n P(V − V ) = k S + PD , (36) 2 2 which can be written as n+1
V
n
−V =k P
−1
S+D
Un+1 + Un 2
Vn+1 + Vn 2
.
(37)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
The quantity ||v||2 = v, v is conserved in the evolution in time, where ∀v, α ∈ R2N : ⎡ ⎤ N−1 1 1
v, α = h ⎣ vt1 α1 + vtj αj + vtN αN ⎦ . 2 2
539
(38)
j=2
From the definitions of P, S and D(u), we conclude that
P −1 Sv, v = 0,
¯ v = 0
D(u)v,
(39)
From this we conclude:
˙v, v =
1 d
v, v = 0. 2 dt
(40)
Hence, the semidiscretized system (35) problem is conservative. Also we can prove the same results for the numerical solution Vn of the proposed scheme (36) by observing that
P −1 S(Vn+1 + Vn ), (Vn+1 + Vn ) = 0,
n = 0, 1, 2, . . . ,
n+1 ¯
D(U)(V + Vn ), (Vn+1 + Vn ) = 0,
n = 0, 1, 2, . . . , n
n+1
∀U ,U
(41)
¯ ∈ R2N . Recall (37) we find immediately ,U
(Vn+1 − Vn ), (Vn+1 + Vn ) = ||Vn+1 ||22 − ||Vn ||22 = 0
(42)
and hence ||Vn+1 || = ||Vn ||,
n = 0, 1, 2, . . . ,
(43)
this proves that the derived scheme posses the conservation property, which ensures the boundedness of the approximations Vnm as n → ∞, ruling out the occurrence of the nonlinear blowup [6,12–14]. 2.4. Stability To study the stability of the numerical method, that is, the sensitivity of the numerical solution to perturbations in the initial data. The von Neumann stability analysis will be used. This method only applicable for linear problems. So we consider the linear version of the CNLS equation: ∂u 1 ∂2 u + A 2 + αAu = 0 ∂t 2 ∂x
(44)
where α = max{z1 , z2 }. The linear version of the proposed scheme is
1 2 1 2 1 n+1 n 2 n+1 n n 1 + δx (Um − Um ) + rAδx (Um + Um ) + kαA 1 + δx (Un+1 m + Um ) = 0, 6 2 6 we assume that Unm = Gn U0 eiβmh ,
where i =
√
−1
(45)
(46)
be the test function, where β ∈ R, U0 ∈ R4 , and G ∈ R4×4 being the amplification matrix. Substituting the test function (46) into (45), we get after some manipulation the amplification matrix G, which is given implicitly by the equation: [γI − ωA]G − [γI + ωA] = 0,
(47)
where μ = sin2
βh , 2
γ =1−
μ 3
k and ω = μ − γα 2
540
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
The von Neumann necessary condition for the stability of a system is maxj |λj | ≤ 1,
j = 1, 2, 3, 4,
where λj are the eigenvalues of G. Recall the skew symmetric property of A, it is easy to see the matrix [γI − ωA] is nonsingular and shares the same set of eigenvalues {γ + ωi, γ − ωi, γ + ωi, γ − ωi} with the matrix [γI + ωA]. Thus the maximal module of the eigenvalues of G is one. Hence, the linearized scheme is nondissipative and the scheme is unconditionally stable. Which means that no restriction on the time step size k can be imposed, but we must choose it in such a way that we get accurate results. 2.5. Accuracy of the scheme The accuracy of the derived scheme (25) can be studied by replacing the numerical solution Unm by the exact solution Doing this the proposed scheme will be of the form
n+1
n+1
1 um − unm um + unm 1 2 1 n 1 + δ2x g + u ) + 1 + + 2 Aδ2x (un+1 δ = 0, (48) m m 6 k 4h 6 x 2
unm .
where g(u) = B(u)u. Taylor’s series expansion of all terms in Eq. (48) about unm can be displayed as
n+1
1 2 um − unm ∂u k ∂2 u k2 ∂2 u h2 ∂3 u 1 + δx + + = + + ··· 6 k ∂t 2 ∂t 2 6 ∂t 2 12 ∂x2 ∂t 1 2 n+1 1 ∂2 u k ∂3 u h2 ∂ 4 u n δ (u + u ) = + + ··· + x m m 4h2 2 ∂x2 4 ∂x2 ∂t 24 ∂x4
n+1
1 um + unm k ∂g h2 ∂2 g kh2 ∂3 g 1 + δ2x g + + ···. =g+ + 6 2 2 ∂t 12 ∂x2 24 ∂t∂x2 Substituting these expansions into (48), this will lead us to the following:
n
n h2 ∂ 4 u ∂u 1 ∂2 u k ∂ ∂u 1 ∂2 u k 2 ∂ 2 u h2 ∂ 3 u n + A 2 + g(u) + A 2 + g(u) + + + + Tm = ∂t 2 ∂x 2 ∂t ∂t 2 ∂x 6 ∂t 2 12 ∂x2 ∂t 24 ∂x4 m m +
h2 ∂ 2 g + ···. 12 ∂x2
By using the differential system (17) all terms inside the brackets are equal to zero, and so we get Tmn =
h 2 ∂ 4 u h2 ∂ 2 g k 2 ∂ 2 u h2 ∂ 3 u + + + + ···, 2 2 6 ∂t 12 ∂x ∂t 24 ∂x4 12 ∂x2
we deduce that the proposed scheme is second order in space and time, it is consistent since the local truncation Tmn tends to zero as h and k approaches zero. According to the Lax theorem, the proposed scheme is convergent since it is consistent and unconditionally stable. 3. Numerical results In this section, we will present numerical results for the proposed scheme (25). The accuracy of this method is tested by looking at their conservation properties and the error using infinity norm. Trapezoidal rule is used to calculate the conserved quantities.
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
541
Table 1 Single soliton T
ER (Crank Nicolson)
ER (Galerkin method)
10 20 30 40
0.037770 0.073790 0.109349 0.145004
0.036107 0.070403 0.104821 0.138414
Fig. 1. Single soliton.
In order to get a reliable solution, discrete conservation laws are important feature in computing smooth soliton solution of the CNLS equations. The accuracy is measured by using the L∞ error norm n n + iU2,m ||)|} ||ER||∞ = max {|(||Ψ1 (xm , tn )|| − ||U1,m 1
We will consider the following problems to highlight the properties of the derived scheme. 3.1. Single soliton In this test we choose the initial condition √ 2α Ψ1 (x, 0) = sech( 2αx) exp i{vx}, 1+e
Ψ2 (x, 0) =
√ 2α sech( 2αx) exp i{vx} 1+e
(49)
where α and e and v are constants. To compute the numerical solution. The following parameters are used xL = −20,
xR = 60,
h = 0.1,
k = 0.01,
v = 1.0,
α = 1.0,
e=1
In Table 1, a comparison with Crank-Nicolson method [7] has been made, we found that the finite element method is more accurate than the Crank- Nicolson method. In Table 2, we display the conserved quantities obtained from the proposed scheme, it is very easy to see that the first conserved quantity is exactly conserved, while the other two are nearly conserved. Fig. 1 shows the evolution of single soliton moving to the right with velocity v = 1. Table 2 Single soliton (conserved quantities) Time
I1
I3
I4
8 16 24 32 40
1.189207 1.189207 1.189207 1.189207 1.189207
0.463807 0.463777 0.463694 0.463727 0.463791
−2.814566 −2.814536 −2.814437 −2.814409 −2.814422
542
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Table 3 Rate of convergence h
L∞
Order
0.2 0.1 0.05 0.025
0.06477 0.01722 0.00436 0.00108
– 1.9112 1.9817 2.0093
Table 4 Three solitons (conserved quantities) T
I1
I2
I3
10 20 30 40
1.985222 1.985222 1.985222 1.985222
0.784914 0.749783 0.785296 0.784572
−1.785707 −1.775935 −1.785222 −1.785310
To test whether the proposed numerical scheme exhibit the expected convergence rate in space we perform some numerical experiments for different values of h and fixed time t = 4 with k = 0.01 to ensure the temporal error is negligible. The rate of convergence is calculated using the formula: rate of convergence ≈
ln(E(h2 )/E(h1 )) ln(h2 / h1 )
where E(h) is the L∞ -error. The error and the convergence rate at t = 4 is shown in Table 3. 3.2. Three solitons interaction To study the interaction of three solitons we take as an initial conditions: Ψ1 (x, 0) =
3 j=1
2αj sech( 2αj xj ) exp i{(vj + xj )} 1+e
Fig. 2. Three solitons interaction.
(50)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
543
Fig. 3. Interaction of two solitons (Manakov equation).
and Ψ2 (x, 0) =
3 j=1
2αj sech( 2αj xj ) exp i{(vj + xj )} 1+e
(51)
which represents three waves. The total sum of the masses asymptotically can be as 2 2αj . 1+e 3
I1 =
j=1
In this test we choose the parameters xL = −20, xR = 60, x1 = 0, x2 = 25, x3 = 50, v1 = 1.0, v2 = 0.1, v3 = −1.0, α1 = 1, α2 = 0.6, α3 = 0.3, e = 2/3. We found that three iterations at most are needed for convergence with tolerance 10−6 . This method is exactly conserves I1 and nearly conserves I3 and I4 , see Table 4. In Fig. 2. we display the interaction scenario of the three waves two of them are moving in the same direction with different velocities, while the third one is moving in the opposite direction. Its is observed that, the three waves interacting and leaving the interaction region unchanged in shape and velocities.
Fig. 4. Reflection of two solitons: v = 0.4, e = 2.
544
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Table 5 Interaction of two solitons (Manakov) velocity and conserved quantities T
v
I1
I4
0 10 20 30 40 50 60
0.2 0.2 0.19 0.34 0.2 0.2 0.2
1.760223 1.760223 1.760223 1.760223 1.760223 1.760223 1.760223
−2.087838 −2.087547 −2.087529 −2.087473 −2.087473 −2.087485 −2.087471
Fig. 5. Transmission scenario: |Ψ1 | and v = 1.6, e = 2/3.
3.3. Soliton collisions of different initial conditions We consider now the initial conditions [24]: Ψ1 (x, 0) = 2α1 sech( 2α1 x + x0 ) exp(iv1 x),
Ψ2 (x, 0) =
2α2 sech( 2α2 x − x0 ) exp(iv2 x)
Fig. 6. Transmission scenario: |Ψ2 | and v = 1.6, e = 2/3.
(52)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
545
Fig. 7. Fussion of two solitons: v = 0.4, e = 0.3.
where amplitudes 2αj , velocities vj and initial phase constant x0 with v1 = −v2 = v/2 (> 0). So the approaching velocity of the two solitons is v. We also fix α1 = 1.2 and α2 = 1. The choices for the initial position parameter x0 can be totally arbitrary. They will not affect the collision outcome as long as x0 is large enough. We fix x0 = 10 in our calculations. The only free parameters left are the cross-phase modulational coefficient e, and the collision velocity v, which will be used as control parameters. In all calculations we choose h = 0.1 and k = 0.01. We will consider different tests. In the first test, we choose e = 1 and v = 0.4. This is a Manakov model which is completely integrable and hence we expect the interaction of the two solitons will be elastic and this is indeed the case, see Fig. 3. The velocity and the conserved quantities before and after the interaction are given in Table 5. In the second test we select e = 2/3 which corresponds to linearly birefringent fibers, and v = 0.4. During collision we observe that the velocity of the right moving soliton steadily decreases, and become negative when it emerges from the collision. This means that this soliton is reflected back by the collision. The same thing happens to the other soliton, it initially moves to the left but turn around after the collision. This reflection scenario has been reported in [24]. The amplitudes of the two solitons also changed after the collision, the larger soliton get even larger and the smaller one gets even smaller. We have observed a small pulses that split off from the solitary wave and propagate along beside it but in the other mode, these waves are called daughter waves. This scenario is displayed in Fig. 4.
Fig. 8. Creation scenario for |Ψ1 | with e = 2 and v = 0.8.
546
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Fig. 9. Creation scenario: |Ψ2 | with e = 2 and v = 0.8.
We repeat the test by taking v = 1.6 and e = 2/3. We have noticed that the two waves pass through each other with some reshaping and radiation shedding and daughter waves are generated. The transmission scenario is displayed in Figs. 5 and 6. In the third test we choose e = 0.3 and v = 0.4, we have noticed in this case that the two solitons fusion into one soliton, see Fig. 7. Finally, we choose e = 2 and v = 0.8, a creation of new soliton has been observed, this scenario is displayed in Figs. 8 and 9. 4. Conclusions We have developed a finite element method for solving the coupled nonlinear Schr¨odinger equations. The developed scheme is second order in both directions space and time, and it is unconditionally stable. The accuracy of the numerical results highlighted the efficiency of the method. The method simulates the interaction picture clearly in the case of three solitons interaction, interaction of two orthogonal solitons [24], we have noticed different scenario for different values of e and v. The derived method can be easily generalized to solve N coupled nonlinear Schr¨odinger equation. Acknowledgments The author thanks the referees for their valuable comments and suggestions. References [1] A. Aydin, B. Karasozen, Multi-symplectic integration of coupled nonlinear Schr¨odinger system with soliton solutions, Int. J. Comput. Math., in press. [2] A. Aydin, B. Karasozen, Symplectic and multisymplectic methods for the coupled nonlinear Schr¨odinger equations with periodic solutions, Comput. Phys. Commun. 177 (7) (2007) 566–583. [3] J. de Frutos, J.M. Sanz-Serna, An easily implementable fourth order method for the time integration of waves problems, Report 1991/2, Universidad De Valladolid, Spain, 1991. [4] D.F. Griffiths, A.R. Mitchell, J. Li Morris, A numerical study of the nonlinear Schr¨odinger equation, NA/52, University of Dundee, 1982. [5] M.S. Ismail, T.R. Taha, A linearly implicit conservative scheme for the coupled nonlinear Schr¨odinger equation, Math. Comp. Simul. 74 (2007) 302–311. [6] M.S. Ismail, S.Z. Alamri, ’Highly accurate finite difference method for coupled nonlinear Schr¨odinger equation, Int. J. Comput. Math. 81 (3) (2004) 333–351. [7] M.S. Ismail, T.R. Taha, Numerical simulation of coupled nonlinear Schr¨odinger equation, Math. Comput. Simul. 56 (2001) 547–562. [8] M.S. Ismail, T. Taha, A finite element solution for the coupled Schr¨odinger equation, in: Proceedings of the 16th IMACS World Congress on Scientific Computation, Lausanne, 2000. [9] M.S. Ismail, Finite difference method with cubic spline for solving nonlinear Schr¨odinger equation, Int. J. Comput. Math. 62 (1996) 101–112. [10] C.R. Menyuk, Stability of solitons in birefringent optical fibers, J. Opt. Soc. Am. B 5 (1998) 392–402. [11] G.M. Muslu, H.A. Erbay, Higher-order split-step Fourier schemes for the generalized nonlinear Schr¨odinger equation, Math. Comput. Simul. 67 (2005) 581–595.
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
547
[12] J.M. Sanz-Serna, J.G. Verwer, Conservative and nonconservative Schr¨odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. [13] A.B. Shamerdan, The numerical treatment of the nonlinear Schr¨odinger equation, Comput. Math. Appl. 19 (1990) 67–73. [14] Q. Sheng, A.Q.M. Khaliq, E.A. Al-Said, Solving the generalized nonlinear Schr¨odinger equation via quartic spline approximation, J. Comput. Phys. 166 (2001) 400–417. [15] W.J. Somnnier, C.I. Christov, Strong coupling of Schr¨odinger equations conservative scheme approach, Math. Comput. Simul. 69 (2005) 314–325. [16] N.H. Sweilam, R.F. Al-Bar, Variational iteration method for coupled nonlinear Schr¨odinger equations, Comput. Math. Appl., in press. [17] J.Q. Sun, X.Y. Gu, Z.Q. Ma, Numerical study of the soliton waves of the coupled nonlinear Schr¨odinger system, Physica D 196 (2004) 311–328. [18] J.Q. Sun, M.Z. Qin, Multi-symplectic methods for the coupled 1D nonlinear Schr¨odinger system, Comput. Phys. Commun. 155 (2003) 221–235. [19] M. Wadati, T. Izuka, M. Hisakado, A coupled nonlinear Schr¨odirnger equation and optical solitons, J. Phys. Soc. Jpn. 61 (7) (1992) 2241–2245. [20] H. Wang, Numerical studies on the split step finite difference method for the nonlinear Schr¨odinger equations, Appl. Math. Comput. 170 (2005) 17–35. [21] A.M. Wazwaz, The variational method for rational solutions for KdV, K(2,2), Burger, and Cubic Boussinesq equation, J. Comput. Appl. Math. 207 (2007) 18–23. [22] A.M. Wazwaz, A study on linear and nonlinear Schr¨odinger equations by the variational iteration method, Chaos, Solitons & Fractals, in press. [23] Y. Xu, C. Shu, Local discontinuous Galerkin method for nonlinear Schr¨odinger equations, J. Comput. Phys. 205 (2005) 72–97. [24] J. Yang, Multisoliton perturbation theory for the Manakov equations and its applications to nonlinear optics, Phys. Rev. E 59 (2) (1999) 2393–2405.
Mathematics and Computers in Simulation 78 (2008) 532–547
Numerical solution of coupled nonlinear Schr¨odinger equation by Galerkin method M.S. Ismail ∗ Department of Mathematics, College of Science, P.O. Box 80203, King Abdulaziz University, Jeddah, Saudi Arabia Received 19 November 2006; received in revised form 15 June 2007; accepted 31 July 2007 Available online 7 August 2007
Abstract The coupled nonlinear Schr¨odinger equation models several interesting physical phenomena presents a model equation for optical fiber with linear birefringence. In this paper we derive a finite element scheme to solve this equation, we test this method for stability and accuracy, many numerical tests have been conducted. The scheme is quite accurate and describe the interaction picture clearly. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Galerkin method; Coupled nonlinear Schr¨odinger equation; Soliton interaction
1. Introduction The coupled nonlinear Schr¨odinger equations (CNLSE) arise in a great variety of physical situations. In fiber communication system, such equations have been shown to govern pulse propagation along orthogonal polarization axes in nonlinear optical fibers and in wavelength-division-multiplexed systems [10,19,24]. These equations also model beam propagation inside crystals or photorefractives as well as water wave interactions. Solitary waves in these equations are often called vector solitons in the literature as they generally contain two components. In all the above physical situations, collision of vector solitons is an important issue. This system has been studied intensively in the past 10 years. It has been shown that, in addition to passing-through collision, vector solitons can also bounce off each other or trap each other. In this paper we will consider the coupled nonlinear Schr¨odinger equations (CNLSE): i
∂Ψ1 1 ∂ 2 Ψ1 + + (|Ψ1 |2 + e|Ψ2 |2 )Ψ1 = 0, ∂t 2 ∂x2
−∞ < x < ∞,
(1)
i
∂Ψ2 1 ∂ 2 Ψ2 + + (e|Ψ1 |2 + |Ψ2 |2 )Ψ2 = 0, ∂t 2 ∂x2
−∞ < x < ∞,
(2)
where Ψ1 and Ψ2 are the wave amplitudes in two polarizations and e is the cross-phase modulation coefficient, with initial conditions: Ψ1 (x, 0) = g1 (x), ∗
Ψ2 (x, 0) = g2 (x)
Tel.: +966 506608949. E-mail address: [email protected].
0378-4754/$32.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.07.003
(3)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
533
and boundary conditions: ∂Ψ1 (x, t) ∂Ψ2 (x, t) = =0 ∂x ∂x
as |x| → ∞.
(4)
For linearly birefringent fibers e = 2/3, and for elliptically birefringent fibers e can take other positive values. When e = 0 this system becomes two decoupled nonlinear Schr¨odinger equations (NLS) and when e = 1, the system is known as Manakov equations. In both cases, the system is integrable. In particular, solitons of one polarization should pass through pulses of the opposite polarization without creating shadows and so the collision is elastic. For other values of e the system is nonintegrable, the collision may be nontrivial and various complex phenomena like reflection, transmission, trapping and creation of new solitary waves can occur. Solitons in optical fibers can be defined as nonlinear pulses that propagate nearly distortion-free for long distances and that undergo elastic collision. Some new methods are used to derive the analytical solution of the CNLS system, Sweilam and Al-Bar [16] used the variational iteration method to derive some approximations to the exact solution of the CNLSE. This method is a powerful method and many researchers are currently used this technique [21,22]. Following the discussion of Wadati et al. [19], we derived the exact solution of the CNLSE which can be given in a compact form as 2 √ 2α v Ψ1 (x, t) = sech( 2α(x − vt)) exp i vx − −α t , (5) 1+e 2 2 √ 2α v Ψ2 (x, t) = ± sech( 2α(x − vt)) exp i vx − −α t , (6) 1+e 2 with the conserved quantities (i) Mass conservations: ∞ I1 = |Ψ1 |2 dx, −∞
I2 =
∞
−∞
|Ψ2 |2 dx.
(ii) Momentum conservation:
∞ 2 ∂Ψj ∂Ψ¯ j I3 = Ψ¯ j i − Ψj dx. ∂x ∂x −∞
(7)
(8)
j=1
(iii) Energy conservation: ⎧ ⎫ 2 ∞ ⎨ 2 2 ⎬ 1 ∂Ψj 1 I4 = − |Ψj |4 − e|Ψ1 |2 |Ψ2 |2 dx. ⎭ 2 ∂x 2 −∞ ⎩ j=1
(9)
j=1
Many numerical methods have been developed for the NLS equation [9,11–14], comparing to few numerical methods for the coupled nonlinear Schr¨odinger equation, most of these methods using finite difference method [5–8,15], recently a multi-symplectic methods are used for solving the CNLS system [1,2,17,18] and a very good results obtained concerning the behavior of the solution and the conserved quantities. Wang [20] presented a numerical solution of the single and coupled nonlinear Schr¨odinger equation using split-step finite difference method. Xu and Shu [23], they solve the NLS and CNLSE using the local discontinuous Galerkin method, they analyzed and test this method, the results are in a very good agreement with the previous existing numerical works. Concerning the finite element method many papers have been published for nonlinear Schr¨odinger equation [3,4]. In this work we are going to solve the CNLSE using one of the most popular finite element method which is called Galerkin method. This method will be analyzed for accuracy and stability. We have found that the scheme is second order in both directions space and time, it is unconditionally stable in the linearized sense. Newton’s method is used to solve the nonlinear system obtained from the full discertization of the coupled nonlinear Schr¨odinger equation. The paper is organized as follows: In Section 2, the derivation of the method and the analysis of the properties of the scheme, like stability, accuracy and the solution of the nonlinear system are given. Numerical tests with various initial
534
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
conditions of the CNLSE like single soliton, interaction of three solitons, interaction of arbitrary initial conditions are reported in Section 3. Conclusions are contained in Section 4. 2. Numerical method We assume that the solution of the system (1) and (2) is negligibly small outside the interval [xL , xR ], and so we consider the coupled nonlinear Schr¨odinger: i
∂Ψ1 1 ∂ 2 Ψ1 + (|Ψ1 |2 + e|Ψ2 |2 )Ψ1 = 0, + ∂t 2 ∂x2
xL < x < xR ,
(10)
i
∂Ψ2 1 ∂ 2 Ψ2 + + (e|Ψ1 |2 + |Ψ2 |2 )Ψ2 = 0, ∂t 2 ∂x2
xL < x < xR ,
(11)
where Ψ1 and Ψ2 are the wave amplitudes in two polarizations and e is the cross-phase modulation coefficient, with initial conditions: Ψ1 (x, 0) = g1 (x),
Ψ2 (x, 0) = g2 (x)
(12)
and boundary conditions: ∂Ψ1 (x, t) ∂Ψ2 (x, t) = = 0, ∂x ∂x
x = xL
and
x = xR .
(13)
For our numerical we decompose the complex functions Ψ1 and Ψ2 in the CNLSE into its real and imaginary parts by writing: Ψ1 (x, t) = u1 + iu2 ,
Ψ2 (x, t) = u3 + iu4
(14)
where (ui , i = 1, . . . , 4) are real functions [4]. On substituting these in Eqs. (1) and (2), the following system is obtained: ∂u1 1 ∂ 2 u2 + z1 u2 = 0, + ∂t 2 ∂x2
∂u2 1 ∂ 2 u1 − z1 u1 = 0, − ∂t 2 ∂x2
∂u3 1 ∂ 2 u4 + z2 u4 = 0, + ∂t 2 ∂x2
∂u4 1 ∂ 2 u3 − − z 2 u3 = 0 ∂t 2 ∂x2
(15)
where z1 = u21 + u22 + e(u23 + u24 ),
z2 = e(u21 + u22 ) + u23 + u24 .
(16)
The above system can be written in a matrix–vector form as ∂u 1 ∂2 u + A 2 + F (u)u = 0, ∂t 2 ∂x where
⎡
u1
⎤
⎢u ⎥ ⎢ 2⎥ u = ⎢ ⎥, ⎣ u3 ⎦ u4
⎡
0
(17)
1
⎢ −1 0 ⎢ A=⎢ ⎣ 0 0 0 0
0
0
⎤
0 0⎥ ⎥ ⎥, 0 1⎦ −1 0
⎡
0
⎢ −z ⎢ 1 F (u) = ⎢ ⎣ 0 0
z1
0
0 0
0 0
0
−z2
0
⎤
0⎥ ⎥ ⎥. z2 ⎦ 0
The finite element method [4,8] for the resulting system (17) can be described as follows. The weak form of (17) is obtained by taking the L2 -inner product with elements of the space H 1 (xL , xR ). After integration by parts, and using the boundary conditions we obtain
∂u 1 ∂u ∂ψ (18) ,ψ − A , + (G(u), ψ) = 0, ∀ψ ∈ H 1 (xL , xR ) ∂t 2 ∂x ∂x
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
535
where G(u) = F (u)u and the L2 -inner product (·, ·) is interpreted pointwise and denoted by xR (f, g) = f (x)g(x) dx. xL
The interval [xL , xR ] is discretized by uniform (N + 1) grid points: xm = xL + mh,
m = 0, 1, . . . , N,
where the grid spacing h is given by h = (xR − xL )/N. Let S h be an N-dimensional subspace of H 1 . We choose the elements of S h to be the piecewise linear hat functions and denoted by {ψm }N m=1 . We assume the approximation of the exact solution u(x, t) is given by U(x, t) =
N
Um (t)ψm (x),
(19)
m=0
where Um , m = 0, 1, . . . , N, are time dependent coefficients which are to be determined. Then the discrete form of (18) is given by N m=0
N
N
m=0
m=0
˙ m (ψm , ψj ) − 1 A Um (ψm (x), ψj ) + (G(U), ψj ) = 0, U 2
j = 0, 1, . . . , N,
(20)
˙ m = ∂Um /∂t, Um = ∂Um /∂x. In order to deal with the nonlinear term in (20), we use product approximation where U [4] by which (G(U), ψm ) =
N
G(Um )(ψm , ψm )
(21)
m=0
Since the trial and the test function are chosen to be the piecewise linear functions: ⎧ ⎪ ⎨ (x − xm−1 )/ h if xm−1 < x ≤ xm , ψm (x) = (xm+1 − x)/ h if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise. We note that the solution (19) is zero outside the interval [xm−1 , xm+1 ], and according to this we can write it as ⎧ ⎪ ⎨ Um−1 (t)ψm−1 (x) + Um (t)ψm (x) if xm−1 < x ≤ xm , U(x) = Um (t)ψm (x) + Um+1 (t)ψm+1 (x) if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise and the nonlinear term as ⎧ ⎪ ⎨ G(Um−1 (t))ψm−1 (x) + G(Um (t))ψm (x) if xm−1 < x ≤ xm , G(U(x, t)) = G(Um (t))ψm (x) + G(Um+1 (t))ψm+1 (x) if xm < x ≤ xm+1 , ⎪ ⎩ 0 otherwise. By using these definitions, the inner products in (20) can be calculated as follows: xm+1 xm xim+1 ˙ m−1 + 4U ˙m +U ˙ m+1 ), ˙ ψm dx = ˙ ψm dx + ˙ ψm dx = 1 (U U U U 6 xm−1 xm−1 xm xm+1 xm xm+1 1 Ux (ψx )m dx = Ux (ψx )m dx + Ux (ψx )m dx = 2 (Um−1 − 2Um + Um+1 ) h xm−1 xm−1 xm
536
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
and
xm+1
xm−1
G(U)ψm dx =
xm
xm−1
G(U)ψm dx +
xm+1
xm
G(U)ψm dx =
1 (G(Um−1 ) + 4G(Um ) + G(Um+1 )). 6
By making use of these integrals in (20), this will lead us to the first order ordinary differential system:
1 2 ˙ 1 1 2 2 1 + δx Um + 2 Aδx Um + 1 + δx G(Um ) = 0, 6 2h 6
(22)
where δ2x Um = Um−1 − 2Um + Um+1 . On expansion of the central difference operator, and imposing the homogeneous Neumann boundary conditions, we will get ˙ = [S + PD(U)]U, PU
(23)
where each of D(u), P and S are a block-diagonal skew symmetric, a block-tridiagonal symmetric and a blocktridiagonal skew symmetric matrices, respectively, and the structure of these matrices are ⎤ ⎤ ⎡ ⎡ −A A 0 0 ··· 0 2I I 0 0 ··· 0 ⎢ I 4I ⎢ A −2A A 0 ··· 0 ⎥ I 0 ··· 0 ⎥ ⎥ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ .. ⎥ .. ⎥ .. .. .. .. .. .. ⎥ ⎢ ⎢ . . . . . . 0 0 0 . ⎥ . ⎥ −1 ⎢ 1⎢ 0 ⎥, ⎥, ⎢ S = P= ⎢ ⎥ ⎥ .. .. .. 6⎢ 2h2 ⎢ ⎢ ... · · · . . . . . . . . . 0 ⎥ ⎢ ... . . . ··· 0 ⎥ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎦ ⎣ 0 ··· ⎣ 0 ··· 0 A −2A A ⎦ 0 I 4I I 0 ··· ··· 0 A −A 0 ··· ··· 0 I 2I D(u) = −diag[B1 (u1 ), B2 (u2 ), . . . , BN (uN )], where
⎡
0
1
⎢ −1 0 ⎢ A=⎢ ⎣ 0 0 0 0
0
0
⎤
0 0⎥ ⎥ ⎥, 0 1⎦ −1 0
⎡
0 ⎢ −(z ) 1 i ⎢ Bi (ui ) = ⎢ ⎣ 0 0
(z1 )i 0 0
0 0 0
0
−(z2 )i
⎤ 0 0 ⎥ ⎥ ⎥, (z2 )i ⎦
z1 = u21 + u22 + e(u23 + u24 ),
0
z2 = e(u21 + u22 ) + u23 + u24 , and 0 and I are (4 × 4) zero and identity matrices, respectively. 2.1. Full discertization In order to solve the system (23), we first discretize the time interval by the grid points t = nk, n = 0, 1, 2, . . ., n is the approximation of the exact solution u(x , t ). Making use of where k is the time step size, we assume that Um m n the substitution: n+1 n ˙ m = Um − Um U k
and Um =
n Un+1 m + Um , 2
which is equivalent to the implicit midpoint rule. This will lead us to the implicit nonlinear scheme:
n+1 1 2 Um + Unm 1 2 n+1 n 2 n+1 n 1 + δx (Um − Um ) + rAδx (Um + Um ) + k 1 + δx G = 0, 6 6 2
(24)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
where r = k/4h2 . The above scheme can be displayed in a matrix vector form as n+1 n n+1 n U + U + U U P(Un+1 − Un ) = k S + PD . 2 2
537
(25)
The resulting system (25) is a block nonlinear tridiagonal one of 4N nonlinear algebraic equations, in the unknown vector Un+1 , and this can be solved by using Newton’s method, or by a fixed point method. 2.2. Newton’s method The solution of the nonlinear system obtained in (25) can be solved by Newton’s method. To describe Newton’s method for our system we do the following. First we rewrite the nonlinear system (25) as
w + Un w + Un n Q(w) = P(w − U ) − k S + PD = 0, (26) 2 2 where w = Un+1 . Newton’s method for this system can be given as w(l+1) = w(l) − J −1 (w(l) )Q(w(l) ),
l = 0, 1, . . . ,
(27)
where J(w) is the Jacobian matrix contains the partial derivatives of Q(w), it is an (N × N) square matrix whose elements are blocks of size (4 × 4). The initial guess vector is assumed as w(0) =Un . Newton’s method is generally expected to give quadratic convergence, provided that a sufficiently accurate starting value is known and J −1 (w) exists. The iterative scheme in (27) is usually executed in two steps: (i) Solve the linear block tridiagonal system: J(w(l) )y = Q(w(l) ) by Crout’s method for the intermediate vector y. (ii) Update your solution by using: w(l+1) = w(l) − y. The two steps (i) and (ii) will be applied iteratively till the following condition is satisfied: ||w(l+1) − w(l) ||∞ ≤ ∈ where ∈ is a prescribed tolerance. 2.3. Conserved quantities In this section we will study the conservation properties of the CNLSE and the analogues property of the numerical method [6,13,14]. In this work we derive the conservation laws: ∞ ∞ |Ψ1 |2 dx = constant and |Ψ2 |2 dx = constant. (28) −∞
−∞
The other quantities will be calculated using a numerical quadrature. To show the system in Eq. (1) satisfies (28), we multiply Eq. (1) by Ψ¯ 1 (the complex conjugate of Ψ1 ) and its complex conjugate by Ψ1 , and then subtract the later
538
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
from the former to obtain
∂Ψ1 ∂Ψ¯ 1 1 ∂ 2 Ψ1 1 ∂2 Ψ¯ 1 2 2 2 2 ¯ Ψ¯ 1 i − Ψ i =0 + + (|Ψ | + e|Ψ | )Ψ + (|Ψ | + e|Ψ | ) Ψ + 1 1 2 1 1 2 1 ∂t 2 ∂x2 ∂t 2 ∂x2
(29)
and this gives ∂Ψ1 ∂Ψ¯ 1 1 ¯ ∂ 2 Ψ1 ∂2 Ψ¯ 1 i Ψ¯ 1 + Ψ1 Ψ1 2 − Ψ1 2 = 0, + ∂t 2 ∂t ∂x ∂x which is equivalent to ∂ 1 ∂ i (Ψ1 Ψ¯ 1 ) + ∂t 2 ∂x
∂Ψ1 ∂Ψ¯ 1 Ψ¯ 1 − Ψ1 ∂x ∂x
= 0.
(30)
By integrating Eq. (30) with respect to x from −∞ to ∞ one can get ∞ ∂ ∞ ∂Ψ¯ 1 1 ¯ ∂Ψ1 2 i |Ψ1 | dx + Ψ1 =0 − Ψ1 ∂t −∞ 2 ∂x ∂x −∞
(31)
and by using the homogenous Neumann boundary conditions ∂Ψ1 /∂x = 0 as |x| → ∞ we find that ∂ ∞ |Ψ1 |2 dx = 0 ∂t −∞ and hence the first conserved quantity: ∞ I1 = |Ψ1 |2 dx = constant
(32)
is obtained. By adopting the same criteria, the second conserved quantity: ∞ I2 = |Ψ2 |2 dx = constant
(33)
−∞
−∞
can be obtained. The exact values of these quantities using (5) and (6) are 2 √ I1 = I2 = 2α. 1+e
(34)
Conservation of I1 and I2 implies L2 boundedness of the solution, and hence no blow up is expecting. To prove the discrete analogue of the conservation of (28), we apply the derived scheme for the first component of the CNLSE to get P v˙ = [S + PD(u)]v, 0 ˜ = where v = [u1 , u2 ]t , A −1
(35)
1 , and D(u) is a block-diagonal matrix: 0
˜ σ2 A, ˜ . . . , σN A] ˜ D(u) = −diag[σ1 A,
with σj = u21,j + u22,j + e(u23,j + u24,j )
where in this case u consists of the first two components of u a way from the nonlinear terms. The corresponding implicit scheme for (35) is Un+1 + Un Vn+1 + Vn n+1 n P(V − V ) = k S + PD , (36) 2 2 which can be written as n+1
V
n
−V =k P
−1
S+D
Un+1 + Un 2
Vn+1 + Vn 2
.
(37)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
The quantity ||v||2 = v, v is conserved in the evolution in time, where ∀v, α ∈ R2N : ⎡ ⎤ N−1 1 1
v, α = h ⎣ vt1 α1 + vtj αj + vtN αN ⎦ . 2 2
539
(38)
j=2
From the definitions of P, S and D(u), we conclude that
P −1 Sv, v = 0,
¯ v = 0
D(u)v,
(39)
From this we conclude:
˙v, v =
1 d
v, v = 0. 2 dt
(40)
Hence, the semidiscretized system (35) problem is conservative. Also we can prove the same results for the numerical solution Vn of the proposed scheme (36) by observing that
P −1 S(Vn+1 + Vn ), (Vn+1 + Vn ) = 0,
n = 0, 1, 2, . . . ,
n+1 ¯
D(U)(V + Vn ), (Vn+1 + Vn ) = 0,
n = 0, 1, 2, . . . , n
n+1
∀U ,U
(41)
¯ ∈ R2N . Recall (37) we find immediately ,U
(Vn+1 − Vn ), (Vn+1 + Vn ) = ||Vn+1 ||22 − ||Vn ||22 = 0
(42)
and hence ||Vn+1 || = ||Vn ||,
n = 0, 1, 2, . . . ,
(43)
this proves that the derived scheme posses the conservation property, which ensures the boundedness of the approximations Vnm as n → ∞, ruling out the occurrence of the nonlinear blowup [6,12–14]. 2.4. Stability To study the stability of the numerical method, that is, the sensitivity of the numerical solution to perturbations in the initial data. The von Neumann stability analysis will be used. This method only applicable for linear problems. So we consider the linear version of the CNLS equation: ∂u 1 ∂2 u + A 2 + αAu = 0 ∂t 2 ∂x
(44)
where α = max{z1 , z2 }. The linear version of the proposed scheme is
1 2 1 2 1 n+1 n 2 n+1 n n 1 + δx (Um − Um ) + rAδx (Um + Um ) + kαA 1 + δx (Un+1 m + Um ) = 0, 6 2 6 we assume that Unm = Gn U0 eiβmh ,
where i =
√
−1
(45)
(46)
be the test function, where β ∈ R, U0 ∈ R4 , and G ∈ R4×4 being the amplification matrix. Substituting the test function (46) into (45), we get after some manipulation the amplification matrix G, which is given implicitly by the equation: [γI − ωA]G − [γI + ωA] = 0,
(47)
where μ = sin2
βh , 2
γ =1−
μ 3
k and ω = μ − γα 2
540
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
The von Neumann necessary condition for the stability of a system is maxj |λj | ≤ 1,
j = 1, 2, 3, 4,
where λj are the eigenvalues of G. Recall the skew symmetric property of A, it is easy to see the matrix [γI − ωA] is nonsingular and shares the same set of eigenvalues {γ + ωi, γ − ωi, γ + ωi, γ − ωi} with the matrix [γI + ωA]. Thus the maximal module of the eigenvalues of G is one. Hence, the linearized scheme is nondissipative and the scheme is unconditionally stable. Which means that no restriction on the time step size k can be imposed, but we must choose it in such a way that we get accurate results. 2.5. Accuracy of the scheme The accuracy of the derived scheme (25) can be studied by replacing the numerical solution Unm by the exact solution Doing this the proposed scheme will be of the form
n+1
n+1
1 um − unm um + unm 1 2 1 n 1 + δ2x g + u ) + 1 + + 2 Aδ2x (un+1 δ = 0, (48) m m 6 k 4h 6 x 2
unm .
where g(u) = B(u)u. Taylor’s series expansion of all terms in Eq. (48) about unm can be displayed as
n+1
1 2 um − unm ∂u k ∂2 u k2 ∂2 u h2 ∂3 u 1 + δx + + = + + ··· 6 k ∂t 2 ∂t 2 6 ∂t 2 12 ∂x2 ∂t 1 2 n+1 1 ∂2 u k ∂3 u h2 ∂ 4 u n δ (u + u ) = + + ··· + x m m 4h2 2 ∂x2 4 ∂x2 ∂t 24 ∂x4
n+1
1 um + unm k ∂g h2 ∂2 g kh2 ∂3 g 1 + δ2x g + + ···. =g+ + 6 2 2 ∂t 12 ∂x2 24 ∂t∂x2 Substituting these expansions into (48), this will lead us to the following:
n
n h2 ∂ 4 u ∂u 1 ∂2 u k ∂ ∂u 1 ∂2 u k 2 ∂ 2 u h2 ∂ 3 u n + A 2 + g(u) + A 2 + g(u) + + + + Tm = ∂t 2 ∂x 2 ∂t ∂t 2 ∂x 6 ∂t 2 12 ∂x2 ∂t 24 ∂x4 m m +
h2 ∂ 2 g + ···. 12 ∂x2
By using the differential system (17) all terms inside the brackets are equal to zero, and so we get Tmn =
h 2 ∂ 4 u h2 ∂ 2 g k 2 ∂ 2 u h2 ∂ 3 u + + + + ···, 2 2 6 ∂t 12 ∂x ∂t 24 ∂x4 12 ∂x2
we deduce that the proposed scheme is second order in space and time, it is consistent since the local truncation Tmn tends to zero as h and k approaches zero. According to the Lax theorem, the proposed scheme is convergent since it is consistent and unconditionally stable. 3. Numerical results In this section, we will present numerical results for the proposed scheme (25). The accuracy of this method is tested by looking at their conservation properties and the error using infinity norm. Trapezoidal rule is used to calculate the conserved quantities.
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
541
Table 1 Single soliton T
ER (Crank Nicolson)
ER (Galerkin method)
10 20 30 40
0.037770 0.073790 0.109349 0.145004
0.036107 0.070403 0.104821 0.138414
Fig. 1. Single soliton.
In order to get a reliable solution, discrete conservation laws are important feature in computing smooth soliton solution of the CNLS equations. The accuracy is measured by using the L∞ error norm n n + iU2,m ||)|} ||ER||∞ = max {|(||Ψ1 (xm , tn )|| − ||U1,m 1
We will consider the following problems to highlight the properties of the derived scheme. 3.1. Single soliton In this test we choose the initial condition √ 2α Ψ1 (x, 0) = sech( 2αx) exp i{vx}, 1+e
Ψ2 (x, 0) =
√ 2α sech( 2αx) exp i{vx} 1+e
(49)
where α and e and v are constants. To compute the numerical solution. The following parameters are used xL = −20,
xR = 60,
h = 0.1,
k = 0.01,
v = 1.0,
α = 1.0,
e=1
In Table 1, a comparison with Crank-Nicolson method [7] has been made, we found that the finite element method is more accurate than the Crank- Nicolson method. In Table 2, we display the conserved quantities obtained from the proposed scheme, it is very easy to see that the first conserved quantity is exactly conserved, while the other two are nearly conserved. Fig. 1 shows the evolution of single soliton moving to the right with velocity v = 1. Table 2 Single soliton (conserved quantities) Time
I1
I3
I4
8 16 24 32 40
1.189207 1.189207 1.189207 1.189207 1.189207
0.463807 0.463777 0.463694 0.463727 0.463791
−2.814566 −2.814536 −2.814437 −2.814409 −2.814422
542
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Table 3 Rate of convergence h
L∞
Order
0.2 0.1 0.05 0.025
0.06477 0.01722 0.00436 0.00108
– 1.9112 1.9817 2.0093
Table 4 Three solitons (conserved quantities) T
I1
I2
I3
10 20 30 40
1.985222 1.985222 1.985222 1.985222
0.784914 0.749783 0.785296 0.784572
−1.785707 −1.775935 −1.785222 −1.785310
To test whether the proposed numerical scheme exhibit the expected convergence rate in space we perform some numerical experiments for different values of h and fixed time t = 4 with k = 0.01 to ensure the temporal error is negligible. The rate of convergence is calculated using the formula: rate of convergence ≈
ln(E(h2 )/E(h1 )) ln(h2 / h1 )
where E(h) is the L∞ -error. The error and the convergence rate at t = 4 is shown in Table 3. 3.2. Three solitons interaction To study the interaction of three solitons we take as an initial conditions: Ψ1 (x, 0) =
3 j=1
2αj sech( 2αj xj ) exp i{(vj + xj )} 1+e
Fig. 2. Three solitons interaction.
(50)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
543
Fig. 3. Interaction of two solitons (Manakov equation).
and Ψ2 (x, 0) =
3 j=1
2αj sech( 2αj xj ) exp i{(vj + xj )} 1+e
(51)
which represents three waves. The total sum of the masses asymptotically can be as 2 2αj . 1+e 3
I1 =
j=1
In this test we choose the parameters xL = −20, xR = 60, x1 = 0, x2 = 25, x3 = 50, v1 = 1.0, v2 = 0.1, v3 = −1.0, α1 = 1, α2 = 0.6, α3 = 0.3, e = 2/3. We found that three iterations at most are needed for convergence with tolerance 10−6 . This method is exactly conserves I1 and nearly conserves I3 and I4 , see Table 4. In Fig. 2. we display the interaction scenario of the three waves two of them are moving in the same direction with different velocities, while the third one is moving in the opposite direction. Its is observed that, the three waves interacting and leaving the interaction region unchanged in shape and velocities.
Fig. 4. Reflection of two solitons: v = 0.4, e = 2.
544
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Table 5 Interaction of two solitons (Manakov) velocity and conserved quantities T
v
I1
I4
0 10 20 30 40 50 60
0.2 0.2 0.19 0.34 0.2 0.2 0.2
1.760223 1.760223 1.760223 1.760223 1.760223 1.760223 1.760223
−2.087838 −2.087547 −2.087529 −2.087473 −2.087473 −2.087485 −2.087471
Fig. 5. Transmission scenario: |Ψ1 | and v = 1.6, e = 2/3.
3.3. Soliton collisions of different initial conditions We consider now the initial conditions [24]: Ψ1 (x, 0) = 2α1 sech( 2α1 x + x0 ) exp(iv1 x),
Ψ2 (x, 0) =
2α2 sech( 2α2 x − x0 ) exp(iv2 x)
Fig. 6. Transmission scenario: |Ψ2 | and v = 1.6, e = 2/3.
(52)
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
545
Fig. 7. Fussion of two solitons: v = 0.4, e = 0.3.
where amplitudes 2αj , velocities vj and initial phase constant x0 with v1 = −v2 = v/2 (> 0). So the approaching velocity of the two solitons is v. We also fix α1 = 1.2 and α2 = 1. The choices for the initial position parameter x0 can be totally arbitrary. They will not affect the collision outcome as long as x0 is large enough. We fix x0 = 10 in our calculations. The only free parameters left are the cross-phase modulational coefficient e, and the collision velocity v, which will be used as control parameters. In all calculations we choose h = 0.1 and k = 0.01. We will consider different tests. In the first test, we choose e = 1 and v = 0.4. This is a Manakov model which is completely integrable and hence we expect the interaction of the two solitons will be elastic and this is indeed the case, see Fig. 3. The velocity and the conserved quantities before and after the interaction are given in Table 5. In the second test we select e = 2/3 which corresponds to linearly birefringent fibers, and v = 0.4. During collision we observe that the velocity of the right moving soliton steadily decreases, and become negative when it emerges from the collision. This means that this soliton is reflected back by the collision. The same thing happens to the other soliton, it initially moves to the left but turn around after the collision. This reflection scenario has been reported in [24]. The amplitudes of the two solitons also changed after the collision, the larger soliton get even larger and the smaller one gets even smaller. We have observed a small pulses that split off from the solitary wave and propagate along beside it but in the other mode, these waves are called daughter waves. This scenario is displayed in Fig. 4.
Fig. 8. Creation scenario for |Ψ1 | with e = 2 and v = 0.8.
546
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
Fig. 9. Creation scenario: |Ψ2 | with e = 2 and v = 0.8.
We repeat the test by taking v = 1.6 and e = 2/3. We have noticed that the two waves pass through each other with some reshaping and radiation shedding and daughter waves are generated. The transmission scenario is displayed in Figs. 5 and 6. In the third test we choose e = 0.3 and v = 0.4, we have noticed in this case that the two solitons fusion into one soliton, see Fig. 7. Finally, we choose e = 2 and v = 0.8, a creation of new soliton has been observed, this scenario is displayed in Figs. 8 and 9. 4. Conclusions We have developed a finite element method for solving the coupled nonlinear Schr¨odinger equations. The developed scheme is second order in both directions space and time, and it is unconditionally stable. The accuracy of the numerical results highlighted the efficiency of the method. The method simulates the interaction picture clearly in the case of three solitons interaction, interaction of two orthogonal solitons [24], we have noticed different scenario for different values of e and v. The derived method can be easily generalized to solve N coupled nonlinear Schr¨odinger equation. Acknowledgments The author thanks the referees for their valuable comments and suggestions. References [1] A. Aydin, B. Karasozen, Multi-symplectic integration of coupled nonlinear Schr¨odinger system with soliton solutions, Int. J. Comput. Math., in press. [2] A. Aydin, B. Karasozen, Symplectic and multisymplectic methods for the coupled nonlinear Schr¨odinger equations with periodic solutions, Comput. Phys. Commun. 177 (7) (2007) 566–583. [3] J. de Frutos, J.M. Sanz-Serna, An easily implementable fourth order method for the time integration of waves problems, Report 1991/2, Universidad De Valladolid, Spain, 1991. [4] D.F. Griffiths, A.R. Mitchell, J. Li Morris, A numerical study of the nonlinear Schr¨odinger equation, NA/52, University of Dundee, 1982. [5] M.S. Ismail, T.R. Taha, A linearly implicit conservative scheme for the coupled nonlinear Schr¨odinger equation, Math. Comp. Simul. 74 (2007) 302–311. [6] M.S. Ismail, S.Z. Alamri, ’Highly accurate finite difference method for coupled nonlinear Schr¨odinger equation, Int. J. Comput. Math. 81 (3) (2004) 333–351. [7] M.S. Ismail, T.R. Taha, Numerical simulation of coupled nonlinear Schr¨odinger equation, Math. Comput. Simul. 56 (2001) 547–562. [8] M.S. Ismail, T. Taha, A finite element solution for the coupled Schr¨odinger equation, in: Proceedings of the 16th IMACS World Congress on Scientific Computation, Lausanne, 2000. [9] M.S. Ismail, Finite difference method with cubic spline for solving nonlinear Schr¨odinger equation, Int. J. Comput. Math. 62 (1996) 101–112. [10] C.R. Menyuk, Stability of solitons in birefringent optical fibers, J. Opt. Soc. Am. B 5 (1998) 392–402. [11] G.M. Muslu, H.A. Erbay, Higher-order split-step Fourier schemes for the generalized nonlinear Schr¨odinger equation, Math. Comput. Simul. 67 (2005) 581–595.
M.S. Ismail / Mathematics and Computers in Simulation 78 (2008) 532–547
547
[12] J.M. Sanz-Serna, J.G. Verwer, Conservative and nonconservative Schr¨odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. [13] A.B. Shamerdan, The numerical treatment of the nonlinear Schr¨odinger equation, Comput. Math. Appl. 19 (1990) 67–73. [14] Q. Sheng, A.Q.M. Khaliq, E.A. Al-Said, Solving the generalized nonlinear Schr¨odinger equation via quartic spline approximation, J. Comput. Phys. 166 (2001) 400–417. [15] W.J. Somnnier, C.I. Christov, Strong coupling of Schr¨odinger equations conservative scheme approach, Math. Comput. Simul. 69 (2005) 314–325. [16] N.H. Sweilam, R.F. Al-Bar, Variational iteration method for coupled nonlinear Schr¨odinger equations, Comput. Math. Appl., in press. [17] J.Q. Sun, X.Y. Gu, Z.Q. Ma, Numerical study of the soliton waves of the coupled nonlinear Schr¨odinger system, Physica D 196 (2004) 311–328. [18] J.Q. Sun, M.Z. Qin, Multi-symplectic methods for the coupled 1D nonlinear Schr¨odinger system, Comput. Phys. Commun. 155 (2003) 221–235. [19] M. Wadati, T. Izuka, M. Hisakado, A coupled nonlinear Schr¨odirnger equation and optical solitons, J. Phys. Soc. Jpn. 61 (7) (1992) 2241–2245. [20] H. Wang, Numerical studies on the split step finite difference method for the nonlinear Schr¨odinger equations, Appl. Math. Comput. 170 (2005) 17–35. [21] A.M. Wazwaz, The variational method for rational solutions for KdV, K(2,2), Burger, and Cubic Boussinesq equation, J. Comput. Appl. Math. 207 (2007) 18–23. [22] A.M. Wazwaz, A study on linear and nonlinear Schr¨odinger equations by the variational iteration method, Chaos, Solitons & Fractals, in press. [23] Y. Xu, C. Shu, Local discontinuous Galerkin method for nonlinear Schr¨odinger equations, J. Comput. Phys. 205 (2005) 72–97. [24] J. Yang, Multisoliton perturbation theory for the Manakov equations and its applications to nonlinear optics, Phys. Rev. E 59 (2) (1999) 2393–2405.