A novel approach via mixed Crank–Nicolson scheme and differential quadrature method for numerical solutions of solitons of mKdV equation

(1)

https://doi.org/10.1007/s12043-019-1751-1

A novel approach via mixed Crank–Nicolson scheme

and differential quadrature method for numerical solutions of solitons of mKdV equation

ALI BA ¸SHAN

Department of Mathematics, Science and Art Faculty, Zonguldak Bulent Ecevit University, 67100 Zonguldak, Turkey

E-mail: alibashan@gmail.com

MS received 28 June 2018; revised 5 November 2018; accepted 16 November 2018;

published online 30 March 2019

Abstract. The purpose of the present study is to obtain numerical solutions of the modified Korteweg–de Vries equation (mKdV) by using mixed Crank–Nicolson scheme and differential quadrature method based on quintic B-spline basis functions. In order to control the effectiveness and accuracy of the present approximation, five well- known test problems, namely, single soliton, interaction of double solitons, interaction of triple solitons, Maxwellian initial condition and tanh initial condition, are used. Furthermore, the error norms L2andL_∞are calculated for single soliton solutions to measure the efficiency and the accuracy of the present method. At the same time, the three lowest conservation quantities are calculated and also used to test the efficiency of the method. In addition to these test tools, relative changes of the invariants are calculated and presented. After all these processes, the newly obtained numerical results are compared with results of some of the published articles.

Keywords. Partial differential equations; differential quadrature method; mKdV equation; solitons; quintic B-splines.

PACS Nos 02.30.Jr; 02.30.Mv; 02.60.x; 02.60.Cb 1. Introduction

Many physical phenomena are described via partial differential equations (PDEs). For this reason, a lot of researchers have investigated the solution of PDEs [1–5].

One of the most famous nonlinear differential equation known as the Korteweg–de Vries equation (KdV) equation in its simplest form is given by

Ut +εU Ux+μUx x x =0, (1)

where subscriptsx andtdenote partial derivatives with respect to space and time, respectively, andεandμare constant parameters.

The KdV equation stems from the study of shallow water waves [6] derived by Korteweg and de Vries to describe shallow water waves of long-wavelength and small-amplitude travelling in canals. It has been proved earlier that this equation has solitary waves as solutions, and hence it can have any number of solitons [7].

The equation has been the simplest nonlinear equation describing two important effects: nonlinearity which

is represented byU U_x and linear dispersion which is represented byUx x x.The nonlinearity ofU Ux tends to localise the wave whereas dispersion spreads the wave out. The stability of solitons is a result of the delicate equilibrium between the two effects of nonlinearity and dispersion [8–11].

One of the most important KdV-type equation is known as modified KdV (mKdV) equation which was first introduced by Miura [12] and is given as follows:

Ut +εU²Ux+μUx x x =0. (2)

The mKdV equation has many physical applications in a wide range of areas such as electrodynamics, elec- tromagnetic waves, elastic media, traffic flow [13,14], fluid dynamics [15,16] and plasma physics [17]. Various methods are used to obtain solutions of the KdV equation [18–21]. Both numerical and analytical solutions of the mKdV equation have been investigated by many researchers [22–30].

Differential quadrature method (DQM) was first introduced by Bellmanet al [31] to obtain the numerical solution of PDEs. Many researchers have developed

(2)

different types of DQMs utilising various base functions such as Legendre polynomials and spline functions [31,32], Hermite polynomials [33], radial basis functions [34], harmonic functions [35], Sinc functions [36,37], B-spline functions [38–40] and modified B-spline functions [41–43].

In this work, quintic B-spline-based Crank–Nicolson DQM (QCN-DQM) is going to be applied to obtain numerical solutions of the mKdV equation.

2. Quintic B-spline DQM

Let us take the grid distributiona = x₁ < x₂ <· · · <

xN = b of a finite interval [a,b] into consideration.

Provided that any given functionU(x)is smooth enough over the solution domain, its derivatives with respect to x at a grid point xi can be approximated by a linear summation of all the functional values in the solution domain, i.e.

U_x⁽^r⁾(xi)= d⁽^r⁾U dx⁽^r⁾|x_i =

N j=1

w_{i j}⁽^r⁾U(xj),

i=1,2, . . . ,N, r =1,2, . . . ,N −1, (3) where r denotes the order of the derivative, w_{i j}⁽^r⁾ represent the weighting coefficients of the rth-order derivative approximation and N denotes the number of grid points in the solution domain. Here, the index j

represents the fact that w⁽_{i j}^r⁾ is the corresponding weighting coefficient of the functional valueU(xj). We need the first- and the third-order derivatives of the func- tionU(x). Therefore, we are going to find the value of eq. (3) forr =1 and 3.

Let Qm(x) be the quintic B-splines with knots at points x_i where the uniformly distributed N grid points are taken as a = x1 < x2 < · · · <

xN = b on the ordinary real axis. In that case, the B-splines{Q₋₁,Q₀, . . . ,Q_N₊₂}form a basis for functions defined over[a,b]. The quintic B-splines Qm(x) are defined by the following relationships:

Q_m(x)= 1 h⁵

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

(x −xm−3)⁵, x ∈ [xm−3,xm−2], (x −xm−3)⁵−6(x−xm−2)⁵, x ∈ [xm−2,xm−1], (x −xm−3)⁵−6(x−xm−2)⁵+15(x−xm−1)⁵, x ∈ [xm−1,xm],

(x−xm−3)⁵−6(x −xm−2)⁵+15(x −xm−1)⁵

−20(x−x_m)⁵, x ∈ [xm,xm+1], (x−x_m₋₃)⁵−6(x −x_m₋₂)⁵+15(x −x_m₋₁)⁵

−20(x−xm)⁵+15(x−xm+1)⁵, x ∈ [xm+1,xm+2], (x−xm−3)⁵−6(x −xm−2)⁵+15(x −xm−1)⁵

−20(x−xm)⁵+15(x−xm+1)⁵−6(x−xm+2)⁵, x ∈ [x_m+2,x_m+3],

0, otherwise,

whereh =xm −xm−1for allm[44].

Using the quintic B-splines as test functions in the fundamental DQM, eq. (3) leads to the equation

∂⁽^r⁾Qm(xi)

∂x⁽^r⁾ =

m+2 j=m−2

w⁽_i_,^r_j⁾Q_m(x_j),

m= −1,0, . . . ,N +2,i =1,2, . . . ,N. (4) An arbitrary choice ofileads to an algebraic equation system

M1W1=1, (5)

whereQ_i,_j denotesQ_i(x_j),

M1 =

⎡

⎢⎢

⎣

Q₋1,−3 Q₋1,−2 Q₋1,−1 Q₋1,0 Q₋1,1

Q₀_,−₂ Q₀_,−₁ Q₀_,₀ Q₀_,₁ Q₀_,₂

... ... ... ... ...

Q_N_+1,N₋₁ Q_N_+1,N Q_N_+1,N₊₁ Q_N_+1,N+2 Q_N_+1,N+3

QN+2,N QN+2,N+1 QN+2,N+2 QN+2,N+3 QN+2,N+4

⎤

⎥⎥

⎦,

W1 =

w⁽_i_,−^r⁾₃ w_i⁽_,−^r⁾₂ · · · w_i⁽_,^r_N⁾₊₃ w⁽_i_,^r_N⁾₊₄T

(6)

(3)

and 1 =

∂⁽^r⁾Q₋1(xi)

∂x⁽^r⁾

∂⁽^r⁾Q0(xi)

∂x⁽^r⁾ · · ·

∂⁽^r⁾QN+1(xi)

∂x⁽^r⁾

∂⁽^r⁾QN+2(xi)

∂x⁽^r⁾ T

. (7)

The weighting coefficientsw_i⁽_,^r_j⁾ related to the ith grid points are determined by solving system (5). System (5) consists of N +8 unknowns and N +4 equations. To have a unique solution for the system, it is necessary to eliminate four unknown terms from the equation system.

By adding the following equations

∂⁽^r⁺¹⁾Q₋1(xi)

∂x⁽^r⁺¹⁾ = 1 j=−3

w⁽_i^r_,_j⁾Q₋₁(xj), (8)

∂⁽^r⁺¹⁾Q0(xi)

∂x⁽^r⁺¹⁾ = 2 j=−2

w⁽_i_,^r_j⁾Q₀(xj), (9)

∂⁽^r⁺¹⁾QN+1(xi)

∂x^(r+1) =

N+3 j=N−1

w⁽_i_,^r_j⁾Q_N₊₁(xj),

∂⁽^r+1⁾QN+2(xi)

∂x⁽^r⁺¹⁾ =

N+4 j=N

w⁽_i_,^r_j⁾Q_N₊₂(x_j), (10)

to system (5), we easily obtain

M2W1 =2, (11)

where

M2 =

⎡

⎢⎢

⎣

Q₋1,−3 Q₋1,−2 Q₋1,−1 Q₋1,0 Q₋1,1

Q_−1,−3 Q_−1,−2 Q_−1,−1 Q_−1,0 Q_−1,1

Q0,−2 Q0,−1 Q0,0 Q0,1 Q0,2

Q₀_,−₂ Q₀_,−₁ Q₀_,₀ Q₀_,₁ Q₀_,₂

Q_1,−1 Q_1,0 Q_1,1 Q_1,2 Q_1,3

... ... ... ... ...

QN+1,N−1 QN+1,N QN+1,N+1 QN+1,N+2 QN+1,N+3

Q_N₊₁_,_N₋₁ Q_N₊₁_,_N Q_N₊₁_,_N₊₁ Q_N₊₁_,_N₊₂ Q_N₊₁_,_N₊₃

QN+2,N QN+2,N+1 QN+2,N+2 QN+2,N+3 QN+2,N+4

Q_N₊₂_,_N Q_N₊₂_,_N₊₁ Q_N₊₂_,_N₊₂ Q_N₊₂_,_N₊₃ Q_N₊₂_,_N₊₄

⎤

⎥⎥

⎦ and

2=

∂⁽^r⁾Q₋1(xi)

∂x^(r)

∂⁽^r⁺¹⁾Q₋1(xi)

∂x^(r+1)

∂⁽^r⁾Q0(xi)

∂x^(r)

∂⁽^r+1⁾Q0(xi)

∂x⁽^r⁺¹⁾

∂⁽^r⁾Q1(xi)

∂x⁽^r⁾ · · · ∂⁽^r⁾QN+1(xi)

∂x⁽^r⁾

∂⁽^r+1⁾QN+1(xi)

∂x⁽^r⁺¹⁾

∂⁽^r⁾QN+2(xi)

∂x⁽^r⁾

∂⁽^r⁺¹⁾QN+2(xi)

∂x⁽^r⁺¹⁾ ^T

.

After using values of quintic B-splines at the grid points and eliminatingw_i⁽_,−^r⁾₃, w_i⁽_,−^r⁾₂, w_i⁽_,^r_N⁾₊₃andw_i⁽_,^r_N⁾₊₄from the system, we obtain an algebraic equation system having a five-banded coefficient matrix of the form

M3W2 =3, (12)

where

M3=

⎡

⎢⎢

⎣

37 82 21

8 33 18 1

1 26 66 26 1

... ... ... ... ...

1 26 66 26 1

1 18 33 8

21 82 37

⎤

⎥⎥

⎦ and

(4)

W₂ =

⎡

⎢⎢

⎢⎣ w⁽_i_,−^r⁾₁

w_i⁽_,^r₀⁾ ...

w_i⁽_,^r_i⁾₋₂ w_i⁽_,^r_i⁾₋₁ w_i^(r)_,_i w_i⁽_,^r_i⁾₊₁ w_i⁽_,^r_i⁾₊₂

...

w_i⁽_,^r_N⁾₊₁ w_i⁽_,^r_N⁾₊₂

⎤

⎥⎥

⎥⎦ .

The non-zero entries of the load vector3are given as

−1 = 1 30

−5Q^(r)₋₁(xi)+h Q^(r₋₁⁺¹⁾(xi)

+40Q⁽₀^r⁾(xi)+8h Q⁽₀^r⁺¹⁾(xi) , 0 = 1

10

5Q⁽₀^r⁾(xi)−h Q⁽₀^r⁺¹⁾(xi) , i−2= Q_i^(r)₋₂(xi),

i−1= Q_i⁽^r₋⁾₁(xi), i = Q_i⁽^r⁾(x_i), i+1= Q_i⁽^r₊⁾₁(xi), i+2= Q_i⁽^r₋⁾₂(xi), N+1= 1

10

5Q⁽_N^r⁾₊₁(xi)+h Q⁽_N^r⁺₊¹₁⁾(xi) , N+2= −1

30

−40Q⁽_N^r⁾₊₁(xi)+8h Q⁽_N^r⁺₊¹₁⁾(xi)

+5Q⁽_N^r⁾₊₂(xi)+h Q⁽_N^r⁺₊¹₂⁾(xi)

. (13)

For example, if we apply the test functions Qm,m =

−1,0, . . . ,N+2, at the first grid pointx1for first-order derivative approximation by the selection ofi =1 and r =1 in eq. (13)

−1 = 1 30

−5Q⁽¹⁾₋₁(x1)+h Q⁽²⁾₋₁(x1)

+40Q⁽₀¹⁾(x1)+8h Q⁽₀²⁾(x1) ,

−1= 1 30

−5 −5

h

+h 20

h²

+40 −50

h

+8h 40

h²

= −109 2h , 0 = 1

10

5Q⁽₀¹⁾(x1)−h Q⁽₀²⁾(x1) , 0 = 1

10

5 −50

h

−h 40

h²

= −29 h , 1 =Q⁽₁¹⁾(x1)=0,

2 =Q⁽₂¹⁾(x1)= 50 h , 3 =Q⁽₃¹⁾(x1)= 5

h, N+1 = 1

10

5Q⁽_N¹⁾₊₁(x₁)+h Q⁽_N²⁾₊₁(x₁) , N+1 = 1

10[5.0+h.0]=0, N+2 = −1

30

−40Q⁽_N^r⁾₊₁(xi)+8h Q⁽_N^r⁺₊¹₁⁾(xi)

+5Q^(r)_N₊₂(xi)+h Q^(r_N⁺¹⁾₊₂ (xi) , N+2 = −1

30 [−40.0+8h.0+5.0+h.0]=0 is obtained and written in the matrix form as

M3

⎡

⎢⎢

⎢⎣ w₁⁽¹_,−⁾₁

w⁽₁¹_,₀⁾ w₁⁽¹_,₁⁾ w_1,2⁽¹⁾ w₁⁽¹_,₃⁾ w₁⁽¹_,₄⁾ ...

w⁽₁¹_,⁾_N₊₁ w⁽₁¹_,⁾_N₊₂

⎤

⎥⎥

⎥⎦

=

⎡

⎢⎢

⎣

−109/2h

−29/h 0 50/h

5/h 0 ...

0 0

⎤

⎥⎥

⎦

. (14)

By following the same idea used before to determine the weighting coefficientsw_k⁽¹_,⁾_j, j= −1,0, . . . ,N +2, at grid pointsxk, 2 ≤ k ≤ N −1, we have obtained the following algebraic equation system:

(5)

M₃

⎡

⎢⎢

⎢⎣ w⁽_k¹_,−⁾₁

...

w⁽_k,k−3¹⁾ w⁽_k¹_,_k⁾₋₂ w⁽_k¹_,_k⁾₋₁ w⁽_k¹_,_k⁾ w⁽_k¹_,_k⁾₊₁ w⁽_k¹_,_k⁾₊₂ w⁽_k¹_,_k⁾₊₃

...

w_k⁽¹_,⁾_N₊₂

⎤

⎥⎥

⎥⎦

=

⎡

⎢⎢

⎣ 0

...

0

−5/h

−50/h 0 50/h

5/h 0

...

0

⎤

⎥⎥

⎦

. (15)

For the last grid point of the domain x_N with the same idea, we determine the weighting coefficients w⁽_N¹_,⁾_j, j = −1,0, . . . ,N+2, and obtain the algebraic equation system as

M3

⎡

⎢⎢

⎢⎣ w⁽_N¹_,−⁾ ₁

w⁽_N¹_,⁾₀ ...

w⁽_N¹_,⁾_N₋₃ w⁽_N¹_,⁾_N₋₂ w⁽_N¹_,⁾_N₋₁ w⁽_N¹_,⁾_N w⁽_N¹_,⁾_N₊₁ w⁽_N¹_,N⁾ ₊₂

⎤

⎥⎥

⎥⎦

=

⎡

⎢⎢

⎣ 0 0 ...

0

−5/h

−50/h 0 29/h 109/2h

⎤

⎥⎥

⎦

. (16)

We can obtain the third-order derivative approximations with similar calculation. Hence, system (12) is solved by the pentadiagonal Thomas algorithm.

3. Numerical discretisation

We have discretised eq. (2) using the forward finite dif- ference and Crank–Nicolson-type schemes. First, eq. (2) is discretised as

Uⁿ⁺¹−Uⁿ

t +μU_3xⁿ⁺¹+U_3xⁿ 2 +ε

U²Ux

n+1

+ U²Ux

n

2 =0. (17)

Equation (17) is rewritten as follows:

2Uⁿ⁺¹+t

μU_3xⁿ⁺¹+ε U²Ux

n+1

=2Uⁿ+t

−μU_3xⁿ −ε

U²U_xn

. (18)

Then, the Rubin and Graves-type linearisation tech- nique [45] is used on the left-hand side of eq. (18) to linearise the nonlinear terms as given below:

(U Ux)ⁿ⁺¹=

Uⁿ⁺¹U_xⁿ+UⁿU_xⁿ⁺¹−UⁿU_xⁿ

, (19)

(U Ux)ⁿ =UⁿU_xⁿ. (20)

Accordingly, we have obtained 2Uⁿ⁺¹

+t

μU_3xⁿ⁺¹+ε U²n

U_xⁿ⁺¹+2UⁿU_xⁿUⁿ⁺¹

=2Uⁿ+t

−μU_3xⁿ +ε U²n

U_xⁿ

. (21)

Let us define some terms used in eq. (21) as Aⁿ_i =

N j=1

w_{i j}⁽¹⁾Uⁿ_j =U_xⁿ_i,

B_iⁿ = N

j=1

w⁽_{i j}³⁾Uⁿ_j =U_3xⁿ_i, (22) whereAⁿ_i andB_iⁿare the first- and third-order derivative approximations of theU functions at thenth time level on pointsx_i, respectively. By substituting definition (22) in eq. (21), we obtain

2U_iⁿ⁺¹+t

⎡

⎣μ N

j=1

w⁽_{i j}³⁾Uⁿ_j⁺¹

+ε

⎛

⎝U_iⁿ2

N j=1

w⁽_{i j}¹⁾Uⁿ_j⁺¹+2U_iⁿAⁿ_iU_iⁿ⁺¹

⎞

⎠

⎤

⎦

=φ_iⁿ, (23)

where

φiⁿ =2U_iⁿ+t

−μB_iⁿ+ε U_iⁿ2

Aⁿ_i , for i =1(1)N.

Then we have reorganised eq. (23) for each grid point as follows:

2+t

μw⁽_{i i}³⁾+ε U_iⁿ2

w⁽¹⁾_{i i} +2U_iⁿAⁿ_i

U_iⁿ⁺¹

+

⎡

⎣ ^N

j=1,i=j

t

μw⁽_{i j}³⁾+ε U_iⁿ2

w_{i j}⁽¹⁾ Uⁿ_j⁺¹

⎤

⎦

=φ_iⁿ. (24)

(6)

By implementing the system of eq. (24) on xi, i = 1(1)N grid points, N equations consisting of N unknowns which are denoted byUⁿ⁺¹will be obtained.

The equation system is shown in the matrix form below:

⎡

⎢⎢

⎣

K₁_,₁ K₁_,₂ · · · K₁_,_N K2,1 K2,2 · · · K2,N

... ... ... ...

KN−1,1 KN−1,2 · · · KN−1,N

K_N_,1 K_N,2 · · · K_N,N

⎤

⎥⎥

⎦

⎡

⎢⎢

⎢⎣ U₁ⁿ⁺¹ U₂ⁿ⁺¹

...

U_N−1ⁿ⁺¹ U_Nⁿ⁺¹

⎤

⎥⎥

⎥⎦

=

⎡

⎢⎢

⎣ φⁿ₁ φⁿ₂ ...

φⁿ_N₋₁ φⁿ_N

⎤

⎥⎥

⎦

. (25)

Then the boundary conditions are applied to the system of eq. (25) and the first and last equations are eliminated from the systems. Hence,

⎡

⎢⎢

⎢⎣

K2,2 K2,3 · · · K2,N−1

K3,2 K3,3 · · · K3,N−1

... ... ... ...

KN−1,2 KN−1,3 · · · KN−1,N−1

⎤

⎥⎥

⎥⎦

⎡

⎢⎢

⎣ U₂ⁿ⁺¹ U₃ⁿ⁺¹

...

U_Nⁿ⁺₋¹₁

⎤

⎥⎥

⎦

=

⎡

⎢⎢

⎢⎣

φ₂ⁿ−K_2,1U₁ⁿ⁺¹−K_2,NU_Nⁿ⁺¹ φ₃ⁿ−K3,1U₁ⁿ⁺¹−K3,NU_Nⁿ⁺¹

...

φⁿ_N₋₁−KN−1,1U₁ⁿ⁺¹−KN−1,NU_Nⁿ⁺¹

⎤

⎥⎥

⎥⎦ (26)

is obtained and solved by the Gauss elimination method easily.

4. Numerical examples

In this section, the five well-known test problems are investigated. The accuracy of the numerical method is checked by using the error norms L2 and L_∞, respectively:

L2= h

N J=1

U^exact_j −(UN)j², L_∞=max

j

U_j^exact−(UN)j. (27)

Moreover, the following lowest three invariants corresponding to the conservation of mass, momentum and energy are computed:

I1= _b

a

Udx, I2 = _b

a

U²dx, I3=

_b

a

U⁴−6μ ε

U2

dx. (28)

4.1 Single soliton

The mKdV equation has an analytic solution given in the following form:

U(x,t)=kpsech

kx−kx0−k³μt

, (29)

where p=

6μ ε

1/2

, (30)

which represents a single soliton originally located at x0 moving to the right with velocity k²μ. Solitons may have positive or negative amplitudes depending on the sign of k but all of them have positive velocities.

We take eq. (29) as initial condition at t = 0 of the form

U(x,0)=kpsech(kx−kx₀), (31) and to allow comparison with earlier works [29,30], we use ε = 3, μ = 1, kp = c = 1.3, x0 = 15 and 0 ≤ x ≤ 200. For the present case, the obtained solution is going to move towards the right, having a constant speed with unchanged amplitude. We have plotted the graphs of the numerical solution of a single soliton with t = 0.025 and N = 1001 fromt = 0 to 100 in figure 1. To make a quantitative comparison, the error norms L2 and L_∞ have been computed and compared with earlier works [29,30] in table 1 until t =10, respectively. It is clearly seen from table1that by using the same parameters (t = 0.025) and less number of grid points(N = 761),the present results are superior. Besides this, by decreasing the time step size fromt=0.025 to 0.001, the error normsL2and L_∞decrease to 1.6×10⁻⁵and 1.0×10⁻⁵, respectively att =10. After that, three lowest invariants,I1, I2 and I3, are computed with the same parameterst=0.025 andN =1001 and compared with earlier works [29,30]

in table2untilt =100.It is seen from table2that the present results are superior, again. It is obviously seen from table2that the three lowest invariants,I₁,I₂andI₃, are changed by less than 2.7×10⁻⁶,−4.4×10⁻⁵and

−1.3×10⁻⁴, respectively, with respect to their original values during the long run t = 100 and so are small enough to accept. To show the accuracy of the present method for a long time, i.e.t =100 with a small time