Numerical Solutions of the Coupled Klein-Gordon ... - Semantic Scholar

Communications in Numerical Analysis 2016 No.2 (2016) 193-210

Available online at www.ispacs.com/cna Volume 2016, Issue 2, Year 2016 Article ID cna-00273, 18 Pages

doi:10.5899/2016/cna-00273 Research Article

Numerical Solutions of the Coupled Klein-Gordon-Schrödinger Equations by Differential Quadrature Methods Thoudam Roshan* Department of Mathematics, Sikkim University, 6th Mile, Samdur, Tadong, Gangtok-737102, Sikkim (India)

Copyright 2016 © Thoudam Roshan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract Numerical solutions of the coupled Klein-Gordon-Schrödinger equations is obtained by using differential quadrature methods based on polynomials and quintic B-spline functions for space discretization and RungeKutta fourth order for time discretization. Stability of the schemes are studied using matrix stability analysis. The accuracy and efficiency of the methods are shown by conducting some numerical experiments on test problems. The motion of single soliton and interaction of two solitons are simulated by the proposed methods. Keywords: Differential quadrature; Klein-Gordon-Schrödinger equations; Multidomain decomposition method; Quintic B-spline, Single soliton.

1 Introduction In this paper, we consider the standard coupled Klein-Gordon-Schrödinger (CKGS) equations 1

𝑖𝜓𝑡 + 2 𝜓𝑥𝑥 + 𝜓𝜙 = 0,

𝑥 ∈ [𝑎, 𝑏], 𝑡 ≥ 0

𝜙𝑡𝑡 − 𝜙𝑥𝑥 + 𝜙 − |𝜓|2 = 0, 𝑥 ∈ [𝑎, 𝑏], 𝑡 ≥ 0

(1.1)

which describes the interactions between conservative complex neutron field and a neutral meson Yukawa in quantum field theory, where 𝜓(𝑥, 𝑡) and 𝜙(𝑥, 𝑡) are complex and real functions respectively and 𝑖 = √−1. The CKGS equations (1.1) are supplemented by prescribing the initial and boundary conditions for 𝜓(𝑥, 𝑡) and 𝜙(𝑥, 𝑡) with 𝜓|𝑡=0 = 𝜓0 (𝑥), 𝜙|𝑡=0 = 𝜙0 (𝑥), 𝜙𝑡 |𝑡=0 = 𝜙1 (𝑥) (1.2)

lim   0,

x  

lim   0

x 

(1.3)

where 𝜓0 (𝑥), 𝜙0 (𝑥) and 𝜙1 (𝑥) are known functions. Then, we can implemented zero (or periodic) boundary conditions. The CKGS equations (1.1) have wide applications in many physical fields as cited in * Corresponding Author. Email address: [email protected]; Tel: 03592-231522 193

Communications in Numerical Analysis 2016 No. 2 (2016) 193-210 http://www.ispacs.com/journals/cna/2016/cna-00273/

194

[1-5]. Theoretical results concerning the existence of solutions, stability and asymptotic behaviour for CKGS system are found in literatures [6-10]. The numerical solutions of CKGS equations (1.1) are however important for understanding the physical behaviour of the system. Hong et al. [11] proposed the explicit multi-symplectic schemes for CKGS equations by concentrating suitable symplectic Runge-Kutta-type methods and symplectic Runge-KuttaNyström-type methods for discretizing every partial derivative in each case. A conservative finite difference method has been developed by Zhang [12] for the numerical solution of CKGS equations. Efficient and accurate numerical methods have been developed by Bao and Yang [13] based on Fourier Pseudospectral discretization for spatial derivatives and adoption of solving the ordinary differential equations in phase space analytically under appropriately chosen transmission conditions between different time intervals or Crank-Nicolson/leap-frog time derivatives. Kong et al. [14] pointed out that the CKGS equations have a natural multi-symplectic structure and proposed the multi-symplectic mid-point scheme for the discretization of the CKGS equations. A comparison of five difference schemes has been made by Hong et al. [15] for the numerical solution of the CKGS equations. Recently Dehghan and Talleil [16] have studied the CKGS equations by Chebyshev Pseudospectral multidomain method. Mittal and Bhatia [17] have also solved the CKGS system numerically by cubic B-spline collocation method. Differential quadrature method (DQM) was originally developed by simple analogy with integral quadrature by Bellman et al. [18]. This method approximates the derivatives of a function using weighted sum of the functional values at discrete points in the whole domain. The important aspect of DQM is to calculate the weighting coefficients. Various test functions have been used in the literature to calculate these weighting coefficients Legendre polynomials and spline functions by Bellman et al. [18, 19]. To improve Bellman’s approach, Quan and Chang [20], proposed and explicit formulation using Lagrange interpolation function as test functions. Shu and Xue [21] used Lagrange interpolated trigonometric functions to compute the weighting coefficients. Shu and Wu [22] presented an implicit approach using radial basis functions as the test functions. Radial basis function based differential quadrature method have used for the numerical solution of coupled Klein-Gordon-Zakharov equations by Dehghan and Nikpour [23]. Exponential cubic Bspline differential quadrature method has been developed for the numerical simulations of one dimensional advection-diffusion equation by Korkmaz and Akmaz [24]. In the recent years, the DQM has been widely popular due to its easy applicability, stability, high accuracy and adaptation with other numerical schemes and many engineering/physics problems have been solved successfully using various differential quadrature methods [25-35]. In this work, we consider an extension of the polynomial based differential quadrature (PDQ) method based on the overlapping domain decomposition algorithm and a differential quadrature method based on quintic B-spline functions for the CKGS equations. Pike and Roe’s [36] fourth-stage RK4 scheme is used to integrate the resulting first order differential equations. A linear matrix stability analysis has been conducted and it is found that the schemes are stable. Efficiency of the methods is tested on the motion of single soliton and the interactions of two solitons. Numerical results obtained by the proposed methods are compared with other results available in the literature. The article is organized as follows. In section 2, the multidomain polynomial based differential quadrature method and the quintic B-spline based differential quadrature methods are presented. The implementations of the methods are presented in section 3. In section 4, we discuss the matrix stability analysis of the schemes. We report the numerical experiments for different test problems in section 5. Finally, a brief conclusion is drawn in section 6. 2 Differential Quadrature Methods (DQMs) DQM can be defined as an approximation to a derivative of a given function by the linear sum of the functional values at specific discrete points over the whole solution domain. Let us consider the grid

International Scientific Publications and Consulting Services


195

distribution 𝑎 = 𝑥0 < 𝑥1 < ⋯ < 𝑥𝑁 = 𝑏 of the finite interval [𝑎, 𝑏]. Then, based on the DQ theory, the 𝑛𝑡ℎ order derivative value 𝑢(𝑛) (𝑥) with respect to 𝑥 at a point 𝑥𝑖 is given by N

w

𝑢(𝑛) (𝑥𝑖 ) =

(n) ij

j 0

u ( x j ) , 𝑖 = 0, 1, ⋯ , 𝑁 and 𝑛 = 1, 2, ⋯ , 𝑁 − 1

(2.4)

(𝑛)

Where 𝑤𝑖𝑗 represents the weight coefficients, and 𝑁 is the number of subintervals in which the solution domain [𝑎, 𝑏] is divided. 2.1 Polynomial based differential quadrature (PDQ) method For polynomial basis functions, a set of Lagrange interpolation functions are employed as test functions. Using the base functions 𝑔𝑘 (𝑥) = (𝑥−𝑥

𝐿(𝑥) 𝑘 )𝐿

, 𝑘 = 0, 1, ⋯ , 𝑁,

(1) (𝑥 ) 𝑘

N

where 𝐿(𝑥) =

 ( x  xl ) and 𝐿(1)(𝑥𝑘 ) = l 0

N

(x

k

 xl ) ,

l  0, l  k

Shu [37] obtained the weight coefficients of the first order derivative as: (1)

𝑤𝑖𝑗 =

𝐿 (1) (𝑥𝑖 ) , for (𝑥𝑖 −𝑥𝑗 )𝐿(1) (𝑥𝑗 )

𝑗 ≠ 𝑖, 𝑖, 𝑗 = 0, 1, ⋯ , 𝑁.

(2.5)

When the problem domain is[−1,1], we can chose Chebyshev-Gauss-Lobatto (CGL) points, 𝑖𝜋

𝑥𝑖 = − cos ( 𝑁 ), 𝑖 = 0, 1, ⋯ , 𝑁. With this choice of grid points, the formulation given in Eq. (2.5) can be modified as [30]: (1)

(−1)𝑖+𝑗 𝑐̅𝑖 , for 𝑗 (𝑥𝑖 −𝑥𝑗 )

(1)

w

𝑤𝑖𝑗 = 𝑐̅

𝑗≠𝑖

(2.6)

N

𝑤𝑖𝑖 = −

j  0, j  i

(1) ij

, for 𝑖 = 𝑗 and for all 𝑖, 𝑗 = 0, 1, ⋯ , 𝑁.

(2.7)

Where 2, 𝑖 = 0, 𝑁 𝑐̅𝑖 = { 1, 𝑖 = 1, 2, ⋯ , 𝑁 − 1.

(2.8)

Calculation of the weight coefficients of second and higher order derivatives can be computed using the following formulae [37]: (𝑛) 𝑤𝑖𝑗

=

(1) (𝑛−1) 𝑛 (𝑤𝑖𝑗 𝑤𝑖𝑖

(𝑛−1)

−

𝑤𝑖𝑗

𝑥𝑖 −𝑥𝑗

) for 𝑗 ≠ 𝑖

(2.9)

and (𝑛)

𝑤𝑖𝑖

N

=−

w

j  0, j  i

(n) ij

for 𝑖 = 𝑗,

(2.10)

where 𝑛 = 2, 3, ⋯ , 𝑁 − 1; 𝑖, 𝑗 = 0, 1, ⋯ , 𝑁. 2.1.1 Multi-domain polynomial base differential quadrature (MD-PDQ) method Following the overlapping multidomain pseudo-spectral technique [16, 38], let 𝑃𝑁 (𝐼) be the space of all algebraic polynomials of degree at most 𝑁 on the interval 𝐼 = [−1,1]. The MD-PDQ method for the differential equation 𝜃𝑢 = 𝐹, where 𝜃 is the differential operator, is to find 𝑢 ∈ 𝑃𝑁 (𝐼) such that



196

𝑗𝜋 𝑁

(𝜃𝑢 − 𝐹)(𝜂𝑗 ) = 0, 𝜂𝑗 = − cos ( ) , 𝑗 = 0, 1, ⋯ , 𝑁,

(2.11)

that satisfies the boundary conditions. In MD-PDQ method, the whole solution domain 𝐼𝑥 = [𝑎, 𝑏] is first 𝜇 𝜇 subdivided into 𝑀uniform subintervals 𝐼𝜇 = [𝑥𝐿 , 𝑥𝑅 ], 𝐼𝑥 = [𝑎, 𝑏] = ⋃𝑀 𝜇 𝐼𝜇 . Without loss of generality, we consider that each subinterval has the same length 𝑙 =

𝑏−𝑎 𝑀+(1−𝑀)(1−cos(𝜋/𝑁))/2

and the same number of CGL

collocation points (𝑁 + 1) is used in each subinterval. In this domain decomposition method, we use overlapping subintervals 𝐼𝜇 , where the first two quadrature points of the interval 𝐼𝜇+1 coincide with the last 𝜇

𝜇

𝜇

𝜇+1

𝑀 1 𝑀 𝜇 two quadrature points of the interval 𝐼𝜇 , i.e., 𝐼 = ⋃𝑀 𝜇=1 𝐼 = ⋃𝜇=1[𝑥0 , 𝑥𝑀 ], 𝑥0 = 𝑎, 𝑥𝑁 = 𝑏, 𝑥𝑁−1 = 𝑥0 𝜇 𝑥𝑁

,

𝜇+1 𝑥1 .

= In each subdomain 𝐼𝜇 , we have 𝜃𝑢𝜇 = 𝐹 where 𝑢𝜇 is the solution over the 𝜇th interval. In order to use the differential quadrature weight coefficients given in equations (2.6)-(2.10) on the arbitrary 𝜇 𝜇 interval 𝐼𝜇 = [𝑥𝐿 , 𝑥𝑅 ] we can use the mapping 𝜇

𝜇

𝑥 +𝑥

2

𝜇

𝜇

𝜂: [𝑥𝐿 , 𝑥𝑅 ] → [−1,1], 𝜂(𝑥) = 𝑥 𝜇 −𝑥𝜇 𝑥 − 𝑥𝑅𝜇 −𝑥𝐿𝜇 , 𝑅

𝐿

𝑅

𝐿

that maps the physical coordinate 𝑥 onto the collocation coordinate 𝜂. Therefore for the derivatives, we will have 𝜕𝑛 𝑢𝜇 (𝑥) 𝜕𝑥 𝑛

2

𝑛 𝑛 𝜇 𝜕 𝑢 (𝜂) 𝜕𝜂𝑛 𝐿

= (𝑥 𝜇 −𝑥𝜇 ) 𝑅

(2.12)

Thus, the weight coefficient matrix of the 𝑛𝑡ℎ order derivative is given by 2 𝑛

(𝑛)

(𝑛)

𝑾[0;𝑁,0;𝑁] = ( 𝑙 ) (𝑤𝑖𝑗 )

(𝑁+1)×(𝑁+1)

.

(2.13)

Now, to satisfy the equations globally, it requires assembling the global solution by virtue of an elementwise construction [39]. This element-wise construction is based on the summation of the local element matrices to form their global representation. Now, we interpolate the function 𝑢 to an arbitrary position in the interval 𝐼𝑥 = [𝑎, 𝑏] that the new grid points can be represented by 𝜇

𝜇+1

1 1 = 𝑥02 , 𝑥𝑁 = 𝑥12 , ⋯ , 𝑥𝑁−1 = 𝑥0 {𝑥𝑗 } = {𝑥01 , ⋯ , 𝑥𝑁−1 𝜇

𝑙

𝜇−1

where 𝑥𝑗 = 2 𝜂(𝑥𝑗

𝑙

, ⋯ , 𝑥𝑁𝑀 }

(2.14)

𝜇−1

) + 2 + 𝑥𝑁−1 . Assuming that 𝑢𝜇 = 𝑢𝜇+1 on the boundary of the overlapping domain,

i.e. 𝑥 ∈ 𝜕(𝐼𝜇 ∩ 𝐼𝜇+1 ), the functional values of 𝑢 at the above grid points can be expressed in matrix form as 𝜇 𝑀 𝑀 𝑇 ̃ = [𝑢̃10 , ⋯ , 𝑢̃1𝑁−1 , 𝑢̃12 , ⋯ , 𝑢̃𝑁−1 𝑼 , ⋯ , 𝑢̃𝑁−1 , 𝑢̃𝑁 ] 𝜇

𝜇

where 𝑢̃𝑗 = 𝑢(𝑥𝑗 ). With these definitions, the functional values of the derivative of 𝑢 at the grid points (2.14) can be expressed as ̃ (𝑛) = 𝑾 ̃(𝑛) ̃, 𝑼 𝑼 [0;𝑘,0;𝑘]

(2.15)

where,



̃(𝑛) 𝑾 [0;𝑘,0;𝑘] 1 ̃0,0 𝑊 ⋮ 1 ̃ 𝑊𝑁−1,0

1 ̃0,𝑁−1 ⋯ 𝑊 ⋱ ⋮ 1 ̃ ⋯ 𝑊𝑁−1,𝑁−1 2 ̃1,0 𝑊 ⋮ 2 ̃ 𝑊𝑁−1,0

1 ̃0,𝑁 𝑊 ⋮ 1 ̃ 𝑊𝑁−1,𝑁 2 ̃1,1 𝑊 ⋮ 2 ̃ 𝑊𝑁−1,1

⋯ ⋱ ⋯

197

2 ̃1,𝑁 𝑊 ⋮ 2 ̃ 𝑊𝑁−1,𝑁

⋱

=

𝜇 ̃1,0 𝑊 ⋮ 𝜇 ̃𝑁−1,0 𝑁

𝜇 ̃1,𝑁 ⋯ 𝑊 ⋱ ⋮ 𝜇 ̃𝑁−1,𝑁 ⋯ 𝑊

⋱

(

𝑀 ̃1,0 𝑊 ⋮ 𝑀 ̃𝑁,0 𝑊

𝑀 ̃1,𝑁 ⋯ 𝑊 ⋱ ⋮ 𝑀 ̃𝑁,𝑁 ⋯ 𝑊 )

(2.16) 𝜇

(𝑛)

̃[0;𝑘,0;𝑘] with 𝑘 = 𝑀(𝑁 − 1) + 1 are where[𝑤 ̃]𝑖𝑗 = [𝑾(𝑛) ]𝑖𝑗 and the remaining components of the matrix 𝑾 zero. 2.2 Quintic B-spline differential quadrature (QB-DQ) method Let 𝑄𝑖 (𝑥) be the quintic B-splines with knots at the point 𝑥𝑖 where the uniformly distributed grid points are chosen as 𝑎 = 𝑥0 < 𝑥1 < ⋯ < 𝑥𝑁 = 𝑏 on the ordinary real axis with ℎ = 𝑥𝑖 − 𝑥𝑖−1 , 𝑖 = 1, ⋯ , 𝑁. The quintic B-spline 𝑄𝑖 (𝑥) at the knots is given by [40]: (𝑥 − 𝑥𝑖−3 )5 , 𝑥 ∈ [𝑥𝑖−3 , 𝑥𝑖−2 ) 5 5 (𝑥 − 𝑥𝑖−3 ) − 6(𝑥 − 𝑥𝑖−2 ) , 𝑥 ∈ [𝑥𝑖−2 , 𝑥𝑖−1 ) 5 5 5 (𝑥 − 𝑥𝑖−3 ) − 6(𝑥 − 𝑥𝑖−2 ) + 15(𝑥 − 𝑥𝑖−1 ) , 𝑥 ∈ [𝑥𝑖−1 , 𝑥𝑖 ) 1 𝑄𝑖 (𝑥) = 5 (𝑥𝑖+3 − 𝑥)5 − 6(𝑥𝑖+2 − 𝑥)5 + 15(𝑥𝑖+1 − 𝑥)5 , 𝑥 ∈ [𝑥𝑖 , 𝑥𝑖+1 ) ℎ (𝑥𝑖+3 − 𝑥)5 − 6(𝑥𝑖+2 − 𝑥)5 , 𝑥 ∈ [𝑥𝑖+1 , 𝑥𝑖+2 ) (𝑥𝑖+3 − 𝑥)5 , 𝑥 ∈ [𝑥𝑖+2 , 𝑥𝑖+3 ) {0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. The quintic B-splines {𝑄−2 , 𝑄−1 , 𝑄0 , 𝑄1 , ⋯ , 𝑄𝑁−1 , 𝑄𝑁 , 𝑄𝑁+1 , 𝑄𝑁+2 } form a basis for functions defined over the domain [𝑎, 𝑏]. Each quintic B-spline covers six elements so that each element is covered by six quintic B-splines. The values of 𝑄𝑖 (𝑥) and its derivative may be tabulated as in Table 1. Using the

𝑥 𝑄𝑖 (𝑥) 𝑄𝑖′ (𝑥) 𝑄𝑖′′ (𝑥) 𝑄𝑖′′′ (𝑥) 𝑄𝑖′𝑣 (𝑥)

𝑥𝑖−3 0 0 0 0 0

Table 1: Coefficients of quintic B-spline and derivatives at points 𝑥𝑖 𝑥𝑖−2 𝑥𝑖−1 𝑥𝑖 𝑥𝑖+1 𝑥𝑖+2 𝑥𝑖+3 1 26 66 26 1 5/ℎ 50/ℎ 0 −50/ℎ −5/ℎ 20/ℎ2 40/ℎ2 −120/ℎ2 40/ℎ2 20/ℎ2 60/ℎ3 −120/ℎ3 0 120/ℎ3 60/ℎ3 4 4 4 4 120/ℎ −480/ℎ 720/ℎ −480/ℎ 120/ℎ4

0 0 0 0 0



198

quintic B-splines as test functions in the fundamental DQM, Eq. (2.4) leads to the equation 𝜕𝑛 𝑄𝑚 (𝑥𝑖 ) 𝜕𝑥 𝑛

m2

w

=

j m2

(n) ij

Qm ( x j ) , 𝑖 = 0,1, ⋯ , 𝑁, 𝑚 = −2, −1, ⋯ , 𝑁 + 1, 𝑁 + 2.

(2.17)

An arbitrary choice of 𝑖 leads to an algebraic equation system (𝑛)

𝑄−2,−4

𝑄−2,−3 𝑄−1,−3 ⋱

𝑄−2,−2 𝑄−1,−2 ⋱ 𝑄𝑁+1,𝑁−1

(

𝑄−2,−1 𝑄−1,−1 ⋱ 𝑄𝑁+1,𝑁 𝑄𝑁+2,𝑁

𝑤𝑖,−4

𝑄−2,0 𝑄−1,0 ⋱

𝑄−1,1 ⋱

𝑄𝑁+1,𝑁+1 𝑄𝑁+2,𝑁+1

𝑄𝑁+1,𝑁+2 𝑄𝑁+2,𝑁+2

(𝑛)

𝑄𝑁+1,𝑁+3 𝑄𝑁+2,𝑁+3

𝑄𝑁+2,𝑁+4 )

𝑤𝑖,−3 ⋮ (𝑛) 𝑤𝑖,𝑁+3

=𝜓

(𝑛)

(𝑤𝑖,𝑁+4 ) (2.18)

where 𝑄𝑖,𝑗 denotes 𝑄𝑖 (𝑥𝑗 ) and 𝜓 = [

𝜕𝑛 𝑄−2 (𝑥𝑖 ) 𝜕𝑛 𝑄−1 (𝑥𝑖 ) 𝜕𝑛 𝑄𝑁+1 (𝑥𝑖 ) 𝜕𝑛 𝑄𝑁+2 (𝑥𝑖 ) 𝑇 , , ⋯ , , 𝜕𝑥 𝑛 ] . 𝑛 𝑛 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝑛

The system (2.18)

consists of (𝑁 + 5) equations and (𝑁 + 9) unknowns. In order to be able to solve Eq. (2.18), it requires four more additional equations. By the addition of the equations 𝜕𝑛+1 𝑄−2 (𝑥𝑖 ) 𝜕𝑥 𝑛+1

0

=

𝜕𝑛+1 𝑄𝑁+1 (𝑥𝑖 ) 𝜕𝑥 𝑛+1

 wij( n)Q 2 ( x j ) , j 4

𝜕𝑛+1 𝑄−1 (𝑥𝑖 ) 𝜕𝑥 𝑛+1

N 3

=

 wij( n)QN 1 ( x j ) and

j  N 1

1

=

w

j  3

𝜕𝑛+1 𝑄𝑁+2 (𝑥𝑖 ) 𝜕𝑥 𝑛+1

( n) ij

Q1 ( x j ) , N 4

=

w jN

(n) ij

QN 2 ( x j ) (𝑛)

to the system (2.18) and using the values of the quintic B-splines at the grid points and eliminating 𝑤𝑖,−4, (𝑛)

(𝑛)

(𝑛)

𝑤𝑖,−3, 𝑤𝑖,𝑁+3 and 𝑤𝑖,𝑁+4 from the system, we obtain an algebraic equation system having penta-diagonal matrix of the form Φ𝑊 = Ψ

(2.19)

where (𝑛)

𝑤𝑖,−2 (𝑛)

37 82 21 8 33 18 1 1 26 66 26 1 26 66 ⋱ ⋱ 1

Φ=

(

𝑤𝑖,−1 ⋮ (𝑛) 𝑤𝑖,𝑗−2 1 26 1 and 𝑊 = ⋱ ⋱ ⋱ 26 66 26 1 1 26 66 26 1 1 18 33 8 21 82 37)

(𝑛)

𝑤𝑖,𝑗−1 (𝑛)

𝑤𝑖,𝑗

.

(𝑛)

𝑤𝑖,𝑗+1 (𝑛)

𝑤𝑖,𝑗+2 ⋮ (𝑛) 𝑤𝑖,𝑁+1 (𝑛)

(𝑤𝑖,𝑁+2 ) The entries of the load vector Ψ are given as: 1

(𝑛)

(𝑛+1)

Ψ−2 = 30 (−5𝑄−2 (𝑥𝑖 ) + ℎ𝑄−2

(𝑛) (𝑛+1) (𝑥𝑖 ) + 40𝑄−1 (𝑥𝑖 ) + 8ℎ𝑄−1 (𝑥𝑖 )),


Communications in Numerical Analysis 2016 No. 2 (2016) 193-210 http://www.ispacs.com/journals/cna/2016/cna-00273/ 1 (𝑛) (𝑛+1) (5𝑄−1 (𝑥𝑖 ) − ℎ𝑄−1 (𝑥𝑖 )), 10 (𝑛) (𝑛) (𝑛) Ψ𝑖−2 = 𝑄𝑖−2 , Ψ𝑖−1 = 𝑄𝑖−1 , Ψ𝑖 = 𝑄𝑖 , 1 (𝑛) (𝑛+1) Ψ𝑁+1 = 10 (5𝑄𝑁+1 (𝑥𝑖 ) − ℎ𝑄𝑁+1 (𝑥𝑖 )),

199

Ψ−1 =

Ψ𝑁+2 =

−1 (𝑛) ( 5𝑄𝑁+2 (𝑥𝑖 ) + 30

(𝑛)

Ψ𝑖+1 = 𝑄𝑖+1 ,

(𝑛)

Ψ𝑖+2 = 𝑄𝑖+2 ,

(𝑛+1) (𝑛) (𝑛+1) ℎ𝑄𝑁+2 (𝑥𝑖 ) − 40𝑄𝑁+1 (𝑥𝑖 ) + 8ℎ𝑄𝑁+1 (𝑥𝑖 )).

(2.20)

The system of equations (2.19) can be solved by using a variant of Thomas algorithm at each grid point (𝑛)

𝑥𝑖 , 𝑖 = 0, 1, ⋯ , 𝑁 and we can get all the weight coefficients 𝑤𝑖𝑗 , for all 𝑖, 𝑗 = 0,1, ⋯ , 𝑁; 𝑛 = 1, 2, ⋯ , 𝑁 − 1. 2.2.1 Calculation of weighting coefficients for second order derivative (2)

We have to calculate only the weighting coefficients 𝑤𝑖,𝑗 , at each grid point 𝑥𝑖 , using equations (2.19) and (2.20), as only the second order partial derivative for space is present in the CKGS equations. For example, applying the test function 𝑄𝑚 , 𝑚 = −2, −1, ⋯ , 𝑁 + 2 at the first grid point 𝑥0 by selecting 𝑖 = 0 and 𝑛 = 2 in equation (2.20), we get all the entries of the load matrix Ψ as: 1 (2) (3) (2) (3) Ψ−2 = (−5𝑄−2 (𝑥0 ) + ℎ𝑄−2 (𝑥0 ) + 40𝑄−1 (𝑥0 ) + 8ℎ𝑄−1 (𝑥0 )) 30 1

20

60

40

120

80

= 30 (−5 (ℎ2 ) + ℎ (− ℎ3 ) + 40 (ℎ2 ) + 8ℎ ( ℎ3 )) = ℎ2, Ψ−1 =

1 1 40 120 8 (2) (3) (5𝑄−1 (𝑥0 ) − ℎ𝑄−1 (𝑥0 )) = (5 ( 2 ) − ℎ ( 3 )) = 2 , 10 10 ℎ ℎ ℎ (2)

Ψ0 = 𝑄0 (𝑥0 ) = −

120 , Ψ1 ℎ2

(2)

40

20

(2)

= 𝑄1 (𝑥0 ) = ℎ2 , Ψ2 = 𝑄2 (𝑥0 ) = ℎ2,

and all Ψ𝑖 = 0, 𝑖 = 3, 4, ⋯ , 𝑁 + 2. (2)

With this load matrix, we can get the weighting coefficients 𝑤0,𝑗 , 𝑗 = −2, −1, ⋯ , 𝑁 + 2, from the following system: Φ𝑊0 = Ψx0 , where (2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

𝑊0 = [𝑤0,−2 , 𝑤0,−1 , ⋯ , 𝑤0,𝑗−2 , 𝑤0,𝑗−1 , 𝑤0,𝑗 , 𝑤0,𝑗+1 , 𝑤0,𝑗+2 , ⋯ , 𝑤0,𝑁+1 , 𝑤0,𝑁+2 ]

𝑇

And 80

8

Ψx0 = [ℎ2 , ℎ2 , −

T 120 40 20 , , , 0, ⋯ , 0] . ℎ2 ℎ2 ℎ2 (2)

With the same idea, the weighting coefficients 𝑤𝑘,𝑗 , 𝑗 = −2, , −1, ⋯ , 𝑁 + 2 at the grid points 𝑥𝑘 , 3 ≤ 𝑘 ≤ 𝑁 − 2 can be obtained from the following system: Φ𝑊𝑘 = Ψ𝑥𝑘 , where (2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

𝑊𝑘 = [𝑤𝑘,−2 , 𝑤𝑘,−1 , ⋯ , 𝑤𝑘,𝑗−2 , 𝑤𝑘,𝑗−1 , 𝑤𝑘,𝑗 , 𝑤𝑘,𝑗+1 , 𝑤𝑘,𝑗+2 , ⋯ , 𝑤𝑘,𝑁+1 , 𝑤𝑘,𝑁+2 ]

𝑇

And 20 40

Ψ𝑥𝑘 = [0, ⋯ ,0, ℎ2 , ℎ2 , −

𝑇 120 40 20 , , , 0, ⋯ ,0] . 2 2 2 ℎ ℎ ℎ



200

(2)

For the last grid point 𝑥𝑁 , the weight coefficients 𝑤𝑁,𝑗 , 𝑗 = −2, −1, ⋯ , 𝑁 + 2 can be determined from the following system: Φ𝑊𝑁 = Ψ𝑥𝑁 , Where (2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

𝑇

𝑊𝑁 = [𝑤𝑁,−2 , 𝑤𝑁,−1 , ⋯ , 𝑤𝑁,𝑗−2 , 𝑤𝑁,𝑗−1 , 𝑤𝑁,𝑗 , 𝑤𝑁,𝑗+1 , 𝑤𝑁,𝑗+2 , ⋯ , 𝑤𝑁,𝑁+1 , 𝑤𝑁,𝑁+2 ] and 20 40

Ψ𝑥𝑘 = [0, ⋯ ,0, ℎ2 , ℎ2 , −

120 8 80 𝑇 , , ] . ℎ2 ℎ2 ℎ2

3 Implementation of the DQMs to CKGS equations In this section we consider the CKGS equations (1.1) with initial conditions (1.2) and homogeneous boundary conditions. We set 𝑦(𝑥, 𝑡) = 𝜙(𝑥, 𝑡), 𝑧(𝑥, 𝑡) = 𝜙𝑡 (𝑥, 𝑡) and 𝜓(𝑥, 𝑡) = 𝑢(𝑥, 𝑡) + 𝑖𝑣(𝑥, 𝑡), where 𝑢(𝑥, 𝑡) and 𝑣(𝑥, 𝑡) are real functions. Using the above substitutions, the CKGS equations (1.1) reduce to the following form: 1

𝑢𝑡 = − 2 𝑣𝑥𝑥 − 𝑦𝑣, 1

𝑣𝑡 = 𝑢𝑥𝑥 + 𝑦𝑢, 2 𝑦𝑡 = 𝑧, 𝑧𝑡 = 𝑦𝑥𝑥 − 𝑦 + 𝑢2 + 𝑣 2 , with initial and boundary conditions given by Eqs. (1.2) & (1.3).

(3.21)

3.1 MD-PDQ method 𝜇 𝜇+1 1 Considering the collocating points 𝒙 = {𝑥01 , ⋯ , 𝑥𝑁−1 , 𝑥12 , ⋯ , 𝑥𝑁−1 , 𝑥1 , ⋯ , 𝑥𝑁𝑀 } of the solution domain [𝑎, 𝑏], let us defined the following: 𝑇

𝜇

𝑀 𝑀 𝑼(𝑡) = [𝑢10 (𝑡), ⋯ , 𝑢1𝑁−1 (𝑡), 𝑢12 (𝑡), ⋯ , 𝑢𝑁−1 (𝑡), ⋯ , 𝑢𝑁−1 (𝑡), 𝑢𝑁 (𝑡)] , 𝜇

𝑇

𝜇

𝑇

1 𝑀 𝑽(𝑡) = [𝑣01 (𝑡), ⋯ , 𝑣𝑁−1 (𝑡), 𝑣12 (𝑡), ⋯ , 𝑣𝑁−1 (𝑡), ⋯ , 𝑣𝑁−1 (𝑡), 𝑣𝑁𝑀 (𝑡)] , 1 𝑀 𝒀(𝑡) = [𝑦01 (𝑡), ⋯ , 𝑦𝑁−1 (𝑡), 𝑦12 (𝑡), ⋯ , 𝑦𝑁−1 (𝑡), ⋯ , 𝑦𝑁−1 (𝑡), 𝑦𝑁𝑀 (𝑡)] , 𝜇

1 𝑀 𝒁(𝑡) = [𝑧01 (𝑡), ⋯ , 𝑧𝑁−1 (𝑡), 𝑧12 (𝑡), ⋯ , 𝑧𝑁−1 (𝑡), ⋯ , 𝑧𝑁−1 (𝑡), 𝑧𝑁𝑀 (𝑡)] and ̃(2) 𝛀=𝑾 . (as given in Eq. (2.16))

𝑇

[0;𝑘,0;𝑘]

Discretizing the system of PDEs (3.21) at the above collocations points and using the above definitions, we obtain a system of nonlinear ODEs which can be written in the form of matrix equations as: 1

𝑼′ (𝑡) = − 2 𝛀 ∙ 𝑽(𝑡) − 𝒀(𝑡) ∘ 𝑽(𝑡), 1

𝑽′ (𝑡) = 2 𝛀 ∙ 𝑼(𝑡) + 𝒀(𝑡) ∘ 𝑼(𝑡),

(3.22)

𝒀′ (𝑡) = 𝒁(𝑡), 𝒁′ (𝑡) = 𝛀 ∙ 𝒀(𝑡) − 𝒀(𝑡) + 𝑼2 (𝑡) + 𝑽2 (𝑡), Where ′ ∙′ and ′ ∘ ′ denotes multiplication of two matrices and component by component multiplication of column matrices respectively and 𝑼2 (𝑡) = 𝑼(𝑡) ∘ 𝑼(𝑡) etc. Using the initial-boundary conditions given in (1.2) & (1.3), the above system of ODEs (3.22) is solved by Pike and Roe’s fourth-stage RK4 scheme [36].



201

3.2 QB-DQ method In this case, we consider the uniform collocation points {𝑥0 , 𝑥1 , ⋯ , 𝑥𝑁 } and define the following: 𝑼(𝑡) = [𝑢0 (𝑡), 𝑢1 (𝑡), ⋯ , 𝑢𝑁 (𝑡)]𝑇 , 𝑽(𝑡) = [𝑣0 (𝑡), 𝑣1 (𝑡), ⋯ , 𝑣𝑁 (𝑡)]𝑇 , 𝒀(𝑡) = [𝑦0 (𝑡), 𝑦1 (𝑡), ⋯ , 𝑦𝑁 (𝑡)]𝑇 and 𝒁(𝑡) = [𝑧0 (𝑡), 𝑧1 (𝑡), ⋯ , 𝑧𝑁 (𝑡)]𝑇 . (2)

By direct calculation of the weight coefficients 𝑤𝑖𝑗 , form the algebraic system (2.19), we defined weight coefficient matrix 𝑾 as: (2)

𝑾 = (𝑤𝑖𝑗 )

(𝑁+1)×(𝑁+1)

, 𝑖, 𝑗 = 0, 1, ⋯ , 𝑁.

Using the above definitions and discretizing the system of PDEs (3.21) at the collocating points {𝑥0 , 𝑥1 , ⋯ , 𝑥𝑁 }, reduced to a system of nonlinear ODEs, which can be written in matrix equation form: 1 2

𝑼′ (𝑡) = − 𝐖 ∙ 𝑽(𝑡) − 𝒀(𝑡) ∘ 𝑽(𝑡), 1 2

𝑽′ (𝑡) = 𝐖 ∙ 𝑼(𝑡) + 𝒀(𝑡) ∘ 𝑼(𝑡),

(3.23)

𝒀′ (𝑡) = 𝒁(𝑡), 𝒁′ (𝑡) = 𝐖 ∙ 𝒀(𝑡) − 𝒀(𝑡) + 𝑼2 (𝑡) + 𝑽2 (𝑡). Using initial and boundary conditions, this system of ODEs is solved by Pike and Roe’s fourth-stage RK4 scheme. 4 Stability analysis Unlike the classical finite difference and finite element methods, Von Neumann stability analysis cannot be applied in DQM discretized system, instead matrix stability or energy stability have been studied for DQMs in the literatures [41-43]. After linearization of systems (3.22) and (3.23), we can rewrite it in the following compact form: 𝑿′ (𝑡) = 𝑨 ∙ 𝑿(𝑡) (4.24) (we use homogeneous boundary conditions) Or 𝑼′ (𝑡) 𝑽′ (𝑡) ( ′ )= 𝒀 (𝑡) 𝒁′ (𝑡)

𝟎 1 𝑩 2

𝟎 (𝛼𝑢0 𝑰

1

− 2 𝑩 −𝛼𝑣0 𝑰 𝟎 𝟎 𝟎 𝛼𝑣0 𝑰

𝛼𝑢0 𝑰 𝟎 𝑩−𝑰

𝑼(𝑡) 𝟎 ∙ (𝑽(𝑡) ) 𝒀(𝑡) 𝑰 𝒁(𝑡) 𝟎)

where, 𝟎 = nullmatrix of order (𝑘 + 1) × (𝑘 + 1) for MD-PDQ and (𝑁 + 1) × (𝑁 + 1) for QB-DQ methods respectively, 𝑰 =identity matrix of order (𝑘 + 1) × (𝑘 + 1)for MD-PDQ and(𝑁 + 1) × (𝑁 + 1)for QB-DQ methods respectively, (2) ̃(2) 𝑩=𝑾 , 𝑖, 𝑗 = 0,1, 2, ⋯ , 𝑁 , for QB-DQ methods [0;𝑘,0;𝑘] , for MD-PDQ and 𝑩 = (𝑤𝑖𝑗 ) (𝑁+1)×(𝑁+1)

respectively, and 𝛼𝑢0 and 𝛼𝑣0 are absolute maximum norms of the initial vectors 𝑼(0) and 𝑽(0) and respectively.



202

Figure 1: Stability region for complex eigenvalues

The stability of a numerical scheme for numerical integration of equation (4.24) depends on the stability of numerical scheme for solving it. If the system of ODEs (4.24) is not stable, then the stable numerical scheme for temporal discretization may not produce the converged solution. The stability of equation (4.24) is dependent on the eigenvalues of coefficient matrix 𝑨, since its exact solution is directly determined by the eigenvalues of 𝑨. Let 𝜆𝑖 be the eigenvalues of the coefficient matrix 𝑨. The stable solution of 𝑿(𝑡) as 𝑡 → ∞ requires: i) if all the eigenvalues are real −2.78 < ∆𝑡𝜆𝑖 < 0 ii) if eigenvalues have only complex components, −2√2 < ∆𝑡𝜆𝑖 < 2√2 iii) if eigenvalues are complex ∆𝑡𝜆𝑖 should be in the region shown in Figure 1. We assume 𝛼𝑢0 = 1.68515 and 𝛼𝑣0 = 1.13452, which are the absolute maximum norms of the initial vectors 𝑼(0) and 𝑽(0) and respectively. For MD-PDQ method we take the following parameter values, 𝑁 = 10, 𝑀 = 40 so that 𝐾 = 362. Then we calculate the eigenvalues of the coefficient matrix 𝑨. The distribution of the eigenvalues for this case is shown in Figure 2 (a). The maximum of the absolute values of the eigenvalues is found to be 130.1762. Therefore our maximum choice of ∆𝑡 is given by ∆𝑡 < 2√2/130.1762 < 0.022. With the same choice of 𝛼𝑢0 and 𝛼𝑣0 , we calculate the eigenvalues of the coefficient matrix 𝑨, for the case of QB-DQ method, where we take 𝑁 = 400. The distribution of the eigenvalues in this case is shown in figure 2 (b). From the figure, we see that most of the eigenvalues are distributed very near to the imaginary axis and only few eigenvalues are distributed near the real axis. Here the maximum of the absolute values of the eigenvalues of the coefficient matrix is found to be 52.9733. Therefore the maximum permissible value of ∆𝑡 that satisfies the stability criteria is given by ∆𝑡 < 2√2/52.9733 < 0.0555. However in our numerical experiments, we choose smaller values of ∆𝑡 in order to get accurate results.

(a) (b) Figure 2: Distribution of eigenvalues (a) MD-PDQ method and (b) QBS-DQ method.



203

5 Numerical Results In this section, we demonstrate the effectiveness and accuracy of the proposed methods. The accuracy of the methods is measured by using 𝐿∞ −error norm which is defined by 𝐿∞ = Max uexact(i)  uapp(i) . 0i  N

The exact solution of the CKGS equations (1.1) is given by [36]: 𝜓(𝑥 − 𝑥0 , 𝑡, 𝑞) =

3√2 4√1−𝑞2 3

𝑠𝑒𝑐ℎ2

𝜙(𝑥 − 𝑥0 , 𝑡, 𝑞) = 4(1−𝑞2 ) 𝑠𝑒𝑐ℎ2

1 2√1−𝑞2 1 2√1−𝑞2

(𝑥 − 𝑞𝑡 − 𝑥0 )𝑒𝑥𝑝 (𝑖 (𝑞𝑥 +

1−𝑞2 +𝑞4 𝑡)) 2(1−𝑞2 )

(𝑥 − 𝑞𝑡 − 𝑥0 )

where |𝑞| < 1 indicates the propagating velocity of wave and 𝑥0 is the initial phase. We consider the following numerical tests: 5.1 Single soliton solution In this test the initial conditions of CKGS system are given as 3√2 1 (𝑥 − 𝑥0 )𝑒𝑥𝑝(𝑖𝑞𝑥), 𝜓0 (𝑥) = 𝑠𝑒𝑐ℎ2 4√1 − 𝑞 2 2√1 − 𝑞 2 3 1 2 (𝑥 − 𝑥0 ), 𝜙0 (𝑥) = 𝑠𝑒𝑐ℎ 4(1 − 𝑞 2 ) 2√1 − 𝑞 2 𝜙1 (𝑥) =

3𝑞 1 1 (𝑥 − 𝑥0 )𝑡𝑎𝑛ℎ (𝑥 − 𝑥0 ). 𝑠𝑒𝑐ℎ2 2 4(1 − 𝑞 2 )3/2 2√1 − 𝑞 2√1 − 𝑞 2

We take the initial phase 𝑥0 = 0 and initial velocity 𝑞 = 0.5 in the domain [−20,60] and solve the problem by the proposed two methods till 𝑡 = 40. The evolution of |𝜓𝑁 (𝑥, 𝑡)| and 𝜙𝑁 (𝑥, 𝑡) at different time levels are depicted in Figures 3 (a) and 3 (b) respectively and the corresponding 𝐿∞ -error norms are reported in Table 2. Table 3 reports the comparison of 𝐿∞ -error norms between our propose methods and various other methods found in [15, 16].

(a) (b) Figure 3: Simulation of the single soliton (a) for |𝜓| and (b) for 𝜙, at different time levels with 𝑞 = 0.5



204

Table 2: The values of 𝐿∞ -errors of different methods in the interval [−20,60] at different time levels 𝑀𝐷 − 𝑃𝐷𝑄 (𝑁 = 10, 𝑀 = 40, 𝐾 = 362, ∆𝑡 𝑄𝐵 − 𝐷𝑄(ℎ = 0.2, ∆𝑡 = 0.001) 𝑡 = 0.001) ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ 𝐶𝑃𝑈 − 𝑇𝑖𝑚𝑒 𝐶𝑃𝑈 − 𝑇𝑖𝑚𝑒 10 1.50958E-6 9.41672E-7 22.703 2.68323E-6 2.13447E-6 9.594 20 8.34757E-6 5.95820E-6 46.735 4.81529E-6 3.68966E-6 19.282 30 5.15955E-5 5.37934E-5 70.281 6.95239E-6 5.34736E-6 29.438 40 1.84882E-4 2.58408E-4 94.344 8.96068E-6 7.67140E-6 38.766

Method MD-PDQ(present)

Table 3: Comparison of the present methods and other methods [15, 16] ‖𝐸𝑟𝑟(|𝜓|)‖ ‖𝐸𝑟𝑟(𝜙)‖∞ ∆𝑡 (N-points, M-subs, K) (10,40,362) 0.001 1.55E-7 9.16E-8 (10,80,722) 0.001 3.33E-9 1.09E-8

QB-DQ (present)

(N=K=400 (h=0.2)

0.025

2.05E-6

6.81E-6

MDPS[16]

(12,12,134) (20,12,220)

0.025 0.001

1.5E-5 7.3E-9

2.1E-5 5.8E-10

SDPS[16]

N=K=220

0.00001

7.3E-8

3.5E-9

FD[15] M-M M-C T-C T-T C-C

N=K=1600(h=0.05)

0.025 3.4E-4 1.3E-4 1.4E-4 3.0E-4 1.4E-4

2.5E-5 6.1E-5 8.1E-5 1.0E-4 8.6E-5

5.2 Collision of two solitons In this example, we consider the following initial conditions for the CKGS system (1.1): 𝜓0 (𝑥) = 𝜓(𝑥 − 𝑥1 , 0, 𝑞1 ) + 𝜓(𝑥 − 𝑥2 , 0, 𝑞2 ) 𝜙0 (𝑥) = 𝜙(𝑥 − 𝑥1 , 0, 𝑞1 ) + 𝜙(𝑥 − 𝑥2 , 0, 𝑞2 ) 𝜙1 (𝑥) = 𝜙𝑡 (𝑥 − 𝑥1 , 0, 𝑞1 )|𝑡=0 + 𝜙𝑡 (𝑥 − 𝑥2 , 0, 𝑞2 )|𝑡=0 where𝑥1 , 𝑥2 are the initial phases and 𝑞1 , 𝑞2 are propagating wave velocities of two solitons. 5.2.1 Collision of symmetric solitons In this test, we study the symmetric collisions i.e. the collision of solitons with equal amplitudes and opposite velocities. We take 𝑥1 = −20, 𝑥2 = 20, 𝑞1 = 0.5, 𝑞2 = −0.5 and solve the problem over the domain [−40,40] and the results are reported in Table 4 in the time interval [0,40] for both methods. Figure 4, shows the symmetric collision of the two solitary waves by applying QB-DQ method at various time levels with space step size ℎ = 0.2 and time step size Δ𝑡 = 0.001. From the figure, we conclude that the two solitons collides at 𝑡 = 40 and propagates in their original directions after collision.



205

(a)

(b)

(c)

(d)

Figure 4: Simulation of the collision of |𝜓| symmetric solitons at different times with 𝑞1 = 0.5, 𝑞2 = −0.5.

𝑇𝑖𝑚𝑒 𝑡 10 20 30 40

Table 4: The values of 𝐿∞ -errors for the symmetric collision of solitons at different time levels 𝑀𝐷 − 𝑃𝐷𝑄 (𝑁 = 10, 𝑀 = 40, 𝐾 = 362, ∆𝑡 𝑄𝐵 − 𝐷𝑄(ℎ = 0.2, ∆𝑡 = 0.001) = 0.001) ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ 7.50073E-7 4.07215E-7 3.33500E-7 2.09654E-7 1.77095E-5 2.75110E-6 5.15559E-7 3.11319E-7 1.07122E-4 3.34164E-5 6.31913E-7 4.56692E-7 3.83663E-4 1.98931E-4 8.66583E-7 6.18802E-7

5.2.2 Collision of two asymmetric solitons Next, we consider two solitons propagating with both different velocities and different directions. We chose the parameters 𝑥1 = 15, 𝑥2 = −15, 𝑞1 = −0.8, 𝑞2 = 0.4 over the space domain [−40,40]. In Figure 5, the graphs of the collision of the solitons at different time levels are depicted with ℎ = 0.2 and Δ𝑡 = 0.001. In Table 5, the results are reported in the time interval [0,40]. Form the graph we conclude that the two solitons fully collides and emerge each other at 𝑡 = 25. But from the graph we see that the interaction is not perfectly elastic as a series of waves are accompanied after the interaction.



206

(a)

(b)

(c)

(d)

Figure 5: Simulation of the collision of |𝜓| for non-symmetric solitons at different times with 𝑞1 = −0.8, 𝑞2 = 0.4.

𝑇𝑖𝑚𝑒 𝑡 10 20 30 40

Table 5: The values of 𝐿∞ -errors for the asymmetric collision of solitons at different time levels 𝑀𝐷 − 𝑃𝐷𝑄 (𝑁 = 10, 𝑀 = 40, 𝐾 = 362, ∆𝑡 𝑄𝐵 − 𝐷𝑄(ℎ = 0.2, ∆𝑡 = 0.001) = 0.001) ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ ‖𝐸𝑟𝑟(|𝜓|)‖∞ ‖𝐸𝑟𝑟(𝜙)‖∞ 5.62934E-5 2.87534E-4 6.98856E-3 7.90144E-3 1.01608E-4 7.88700E-3 1.41011E-2 1.69570E-2 9.78951E-2 4.83740E-2 1.00919E-1 6.66452E-2 9.18392E-1 9.17535E-1 8.18325E-1 8.19442E-1

6 Conclusion In this paper, we have developed the MD-PDQ and QB-DQ methods for the finding the numerical solutions of the CKGS system. Using the multidomain scheme in spatial variable reduce memory requirements as it leads to sparse matrix and also allows for best use of parallel computers. In both the methods we used Pike and Roe’s fourth-stage RK4 scheme as the function contains no explicit dependent on 𝑡. Both the methods are tested by simulating the motion of single soliton and interactions of two solitons. In order to show the efficiency and accuracy, we compute the 𝐿∞ −norm and CPU time and compared with some earlier works. We conclude that our schemes produce better results while taking less CPU time. Hence, the combination of MD-PDQ or QB-DQ method in space and Pike and Roe’s RK4 scheme in time, gives another solution procedure in terms of accuracy and less computational cost for solving CKGS system.



207

Acknowledgement The author is very thankful to the reviewer(s) for their valuable comments and suggestions to improve the quality of the paper. This research is supported by the University Grants Commission, Government of India, under the scheme UGC-BSR Research Start-Up-Grant (No. F. 30-95/2015(BSR)). References [1] C. Bardeen, Feedback quantum control of molecular electronic population transfer, Chem. Phys. Lett, 280 (1997) 151-158. http://dx.doi.org/10.1016/S0009-2614(97)01081-6 [2] A. W. Schreiber, R. Rosenfelder, First order variational calculation of form factor in a scalar nucleonmeson theory, Necl. Phys, A 601 (1996) 397-424. http://dx.doi.org/10.1016/0375-9474(96)00016-4 [3] X. Y. Tang, W. Ding, The general Klein-Gordon-Schrödinger system: modulational instability and exact solutions, Phys. Scr. 77 (2008) 1-8. http://dx.doi.org/10.1088/0031-8949/77/01/015004 [4] Q. F. Wang, Theoretical issues of controlling nucleus in Klein-Gordon Schrödinger dynamics with perturbation in control field, Appl. Math. Comput, 206 (2008) 276-289. http://dx.doi.org/10.1016/j.amc.2008.08.057 [5] Y. P. Wang, D. F. Xia, Generalized solitary wave solutions for the Klein-Gordon-Schrödinger equations, Comput. Math. Appl, 58 (2009) 2300-2306. http://dx.doi.org/10.1016/j.camwa.2009.03.012 [6] I. Fukuda, M. Tsutsumi, On coupled Klein-Gordon-Schrödinger equations I, Bull. Sci. Engg. Res. Lab. Waseda Univ, 69 (1975) 51-62. [7] I. Fukuda, M. Tsutsumi, On coupled Klein-Gordon-Schrödinger equations II, J. Math. Anal. Appl. 66 (1978) 358-378. http://dx.doi.org/10.1016/0022-247X(78)90239-1 [8] K. N. Lee, B. X. Wang, Global attractors for the Klein-Gordon-Schrödinger equations in unbounded domain, J. Diff. Equat, 170 (2001) 281-316. http://dx.doi.org/10.1006/jdeq.2000.3827 [9] M. Ohta, Stability of stationary states for the coupled Klein-Gordon-Schrödinger equations, Nonlinear Anal, 27 (1996) 455-461. http://dx.doi.org/10.1016/0362-546X(95)00017-P [10] T. Ozawa, Y. Tsutsumi, Asymptotic behaviour of solutions for the coupled Klein-Gordon-Schrödinger equations, Adv. Stud. Pure Math, 23 (1994) 295-305.



208

[11] J. Hong, S. Jiang, C. Li, Explicit multi-symplectic method for Klein-Gordon-Schrödinger equations, J. Comput. Phys, 228 (2009) 3517-3532. http://dx.doi.org/10.1016/j.jcp.2009.02.006 [12] L. Zang, Convergence of a conservative difference scheme for a class of Klein-Gordon-Schrödinger equations in one space dimension, Appl. Math. Comput, 163 (2005) 343-355. http://dx.doi.org/10.1016/j.amc.2004.02.010 [13] W. Bao, L. Yang, Effecient and accurate numerical methods for the Klein-Gordon-Schrödinger equations, J. Comput. Phys, 225 (2007) 1863-1893. http://dx.doi.org/10.1016/j.jcp.2007.02.018 [14] L. Kong, J. Hong, R. Liu, Long term numerical simulation of the interaction between a neutron field and neutral meson field by a symplectic-preserving scheme, J. Phys. A 41 (2008) 255-207. http://dx.doi.org/10.1088/1751-8113/41/25/255207 [15] J. Hong, S. Jiang, L. Kong, C. Li, Numerical comparison of five difference schemes for coupled KleinGordon-Schrödinger equations in quantum physics, J. Phys. A Math. Theor, 40 (2007) 9125-9135. http://dx.doi.org/10.1088/1751-8113/40/30/030 [16] M. Dehghan, A. Taleei, Numerical solution of the Yukawa-coupled-Klein-Gordon-Schrödinger equations via a Chebyshev pseudospectral multidomain method, Appl. Math. Model, 36 (2012) 23402349. http://dx.doi.org/10.1016/j.apm.2011.08.030 [17] R. C. Mittal, R. Bhatia, Numerical solution of nonlinear Klein-Gordon equations by cubic B-spline collocation method, International J. Computer Mathematics, 92 (2015) 2139-2159. http://dx.doi.org/10.1080/00207160.2014.970182 [18] R. Bellman, B. G. Kashef, J. Casti, Differential quadrature: A technique for the rapid solution of nonlinear partial differential equations, J. Comput. Phys, 10 (1972) 40-52. http://dx.doi.org/10.1016/0021-9991(72)90089-7 [19] R. Bellman, B. G. Kashef, E. S. Lee, R. Vasudevan, Differential quadrature and splines, Computers & Mathematics with Applications, 1 (1975) 371-376. http://dx.doi.org/10.1016/0898-1221(75)90038-3 [20] J. R. Quan, C. Chang, New insights in solving distributed system equations by the quadrature methodI analysis, Computers & Chemical Engineering, 13 (1989) 779-788. http://dx.doi.org/10.1016/0098-1354(89)85051-3 [21] C. Shu, H. Xue, Explicit computation of weighting coefficients in the harmonic differential quadrature, Journal of Sound and Vibration, 204 (1997) 549-555. http://dx.doi.org/10.1006/jsvi.1996.0894 [22] C. Shu, Y. L. Wu, Integrated radial basis function-based differential quadrature method and its performance, International Journal for Numerical Methods in Fluids, 53 (2007) 969-987. http://dx.doi.org/10.1002/fld.1315



209

[23] M. Dehghan, A. Nikpour, The solitary wave solution of coupled Klein-Gordon-Zakharov equations via. Two different methods, Computer Physics Communications, 184 (9) (2013) 2145-2158. http://dx.doi.org/10.1016/j.cpc.2013.04.010 [24] A. Korkmaz, H. K. Akmaz, Numerical simulations for Transport of Conservative Pollutants, Selcuk Journal of Applied Mahematics, 16 (1) (2015). http://sjam.selcuk.edu.tr/sjam/article/view/435 [25] A. Başhan, S. B. G. Karakoc, T. Geyikli, Approximation of the KdVB equation by the quintic B-spline differential quadrature method, Kuwait Journal of Science, 42 (2) (2015). http://journals.ku.edu.kw/kjs/index.php/KJS/article/view/198 [26] A. Korkmaz, Numerical algorithms for solutions of Korteweg-de Vries equation, Numerical methods for partial differential equations, 26 (6) (2010) 1502-1521. http://dx.doi.org/10.1002/num.20505 [27] A. Korkmaz, I. Dag, Quartic and quintic B-spline methods for advection-diffusion equation, Applied Mathematics and Computation, 274 (2016) 208-219. http://dx.doi.org/10.1016/j.amc.2015.11.004 [28] A. Korkmaz, I. Dag, A differential quadrature algorithm for simulation of nonlinear Schrödinger equation, Computers & Mathematics with Application, 56 (9) (2008) 2222-2234. http://dx.doi.org/10.1016/j.camwa.2008.03.047 [29] B. Saka, I. Dag, Y. Dereli, A Korkmaz, Three different methods for numerical solutions of EW equation, Engineering Analysis with Boundary Elements, 32 (2008) 556-566. http://dx.doi.org/10.1016/j.enganabound.2007.11.002 [30] A. Korkmaz, I. Dag, A differential quadrature algorithm for nonlinear Schrödinger equation, Nonlinear Dynamics, 56 (1-2) (2009). http://dx.doi.org/10.1007/s11071-008-9380-0 [31] A. Korkmaz, I. Dag, Numerical simulations of Complex Modified KdV equation using polynomial differential quadrature method, J. of Mathematics and Statistics, 10 (2011) 1-13. [32] A. Korkmaz, A. M. Aksoy, I. Dag, Quartic B-spline differential quadrature method, International Journal of Nonlinear Science, 11 (4) (2011) 403-411. http://internonlinearscience.org/upload/papers/20110915034020930.pdf [33] R. C. Mittal, R. Jiwari, Differential quadrature method for two dimensional Burgers’ equations, Int. J. for Comput. Methods in Eng. Science and Mech, 10 (2009) 450-459. http://dx.doi.org/10.1080/15502280903111424 [34] R. Jiwari, S. Pandit, R. C. Mittal, Numerical simulation of two-dimensional sine-Gordon solitons by differential quadrature method, Computer Physics Communications, 183 (2012) 600-616. http://dx.doi.org/10.1016/j.cpc.2011.12.004



210

[35] R. C. Mittal, R. Bhatia, A numerical study of two dimensional hyperbolic telegraph equation by modified B-spline differential quadrature method, Appl. Math. Comput, 244 (2014) 976-997. http://dx.doi.org/10.1016/j.amc.2014.07.060 [36] J. Pike, PL. Roe, Accelerated convergence of Jameson’s finite volume Euler scheme using Van Der Houwen Integrators, Comput. Fluids, 13 (1985) 223-236. http://dx.doi.org/10.1016/0045-7930(85)90027-1 [37] C. Shu, Differential Quadrature and its Application in Engineering, Springer-Verlag London Ltd., Great Britain, (2000). http://dx.doi.org/10.1007/978-1-4471-0407-0 [38] A. Taleei, M. Dehghan, A pseudo-spectral method that uses an overlapping multidomain technique for the numerical solution of sine-Gordon equation in one and two spatial dimensions, Mathematical Methods in the Applied Sciences, 37 (13) (2014) 1909-1923. http://dx.doi.org/10.1002/mma.2943 [39] F.X. Giraldo, J. B. Perot, P. F. Fischer, A spectral element semi-Lagrange (SESL) method for the spherical shallow water equations, J. Comput. Phys, 190 (2003) 623-650. http://dx.doi.org/10.1016/S0021-9991(03)00300-0 [40] R. C. Mittal, R. K Jain, Numerical solution of General Rosenau-RLW equation using Quintic B-splines Collocation Method, Communications in Numerical Analysis, Article ID cna-00129, 2012 (2012) 19 pages. http://dx.doi.org/10.5899/2012/cna-00129 [41] S. Tomasiello, Stability and accuracy of the iterative differential quadrature method, International Journal for Numerical Methods in Engineering, 58 (6) (2003) 1277-1296. http://dx.doi.org/10.1002/nme.815 [42] S. Tomasiello, Numerical stability of DQ solutions of wave problems, Numerical Algorithms, 57 (3) (2011) 289-312. http://dx.doi.org/10.1007/s11075-010-9429-2 [43] S. Tomasiello, Stability and accuracy of DQ-based step-by-step integration methods for structural dynamics, Applied Mathematical Modelling, 37 (5) (2013) 3426-3435. http://dx.doi.org/10.1016/j.apm.2012.07.005 [44] J. Xia, M. Wang, Exact solitary solution of coupled Klein-Gordon-Schrödinger equations, Appl. Math. Mech, 23(2003) 52-57.
