# A Numerical Solution for Hirota-Satsuma Coupled KdV Equation

Authors: M. S. Ismail; H. A. Ashi
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A Petrov-Galerkin method and product approximation technique are used to solve numerically the Hirota-Satsuma coupled Korteweg-de Vries equation, using cubicB-splines as test functions and a linear B-spline as trial functions. The implicit midpoint rule is used to advance the solution in time. Newton’s method is used to solve the block nonlinear pentadiagonal system we have obtained. The resulting schemes are of second order accuracy in both directions, space and time. The von Neumann stability analysis of the schemes shows that the two schemes are unconditionally stable. The single soliton solution and the conserved quantities are used to assess the accuracy and to show the robustness of the schemes. The interaction of two solitons, three solitons, and birth of solitons is also discussed.


---

## Body

## 1. Introduction

In 1981, Hirota and Satsuma [1] introduced the coupled Korteweg-de Vries equation (CKdV) as follows:
(1)
∂
u
∂
t
=
a
(
∂
3
u
∂
x
3
+
6
u
∂
u
∂
x
)
+
2
b
v
∂
v
∂
x
,
∂
v
∂
t
=
-
∂
3
v
∂
x
3
-
3
u
∂
v
∂
x
,
where a and b are arbitrary constants. The CKdV equation describes interactions of two long waves with different dispersion relations. Namely, it is connected with most types of long waves with weak dispersion, internal, acoustic, and planetary waves in geophysical hydrodynamics [2, 3]. By using Hirota method [1, 4], the single solitary wave solution of this system is
(2)
u
(
x
,
t
)
=
2
λ
2
sec
h
2
(
ξ
)
,
v
(
x
,
t
)
=
1
2
w
sech
(
ξ
)
,
ξ
=
λ
(
x
-
λ
2
t
)
+
1
2
log
⁡
⁡
(
w
)
,
w
=
-
b
8
(
4
a
+
1
)
λ
4
,
where λ is an arbitrary constant. The two and three solitons solutions for a
=
1
/
2 can be found in [1]. The CKdV system has three conserved quantities
(3)
I
1
=
∫
-
∞
∞
‍
u
d
x
,
I
2
=
∫
-
∞
∞
‍
(
u
2
+
2
3
b
v
2
)
d
x
,
I
3
=
∫
-
∞
∞
‍
[
(
1
+
a
)
(
u
3
-
1
2
u
x
2
)
+
b
(
u
v
2
-
v
x
2
)
]
d
x
.The proof can be found in [1, 5, 6].The CKdV equation has been discussed analytically by many authors; Kaya and Inan [7] used Adomian decomposition method to solve this system. Fan used tanh method to find some traveling wave solution [8]. Assas [9] solved this system using variational iteration method. Abbasbandy [10] used homotopy analysis method to solve the generalized CKdV system.The numerical solutions of coupled nonlinear systems are very interesting and important in applied science, for example, the coupled nonlinear Schrödinger equation which admits soliton solution and it has many applications in communication and optical fibers; this system has been solved numerically by Ismail using finite difference and finite element methods [11–13]. The CKdV has been also considered numerically by some researchers; Halim et al. [2, 3] have introduced a numerical scheme for general CKdV systems. The scheme is valid for solving an arbitrary number of equations with arbitrary constant coefficients; the method is applied to the Hirota system and compared with its known explicit solution investigating the influence of initial conditions and grid sizes on accuracy. Wazwaz [14] produced a finite difference scheme for the periodic initial boundary problem of the CKdV system. Ismail [6] solved this system using collocation method and quintic splines; Kutluay and Ucar [15] solved this system using a quadratic B-spline Galerkin approach.In this work, we are going to solve the CKdV equation using Petrov-Galerkin method [4]. We choose the cubic B-spline as test functions and the linear B-spline as trial (basis) functions. Implicit midpoint rule is used in the time direction. Newton’s method is used to solve the nonlinear block pentadigonal system obtained from the schemes we have derived. The von Neumann stability analysis shows that the scheme is unconditionally stable. Regarding the accuracy, the scheme is of second order in space and time.The paper is organized as follows. In Section2, Petrov-Galerkin method is used to derive a numerical method for the CKdV equation; a coupled nonlinear pentadigonal system is obtained. Analysis of the method is given in Section 3. In Section 4, product approximation technique is used to derive a second method for solving the CKdV equation. Numerical results of various tests are contained in Section 5. We recap and sum up our conclusions in Section 6.


## 2. Numerical Method

To derive numerical method for the CKdV system, we consider the initial boundary value problem(4)
∂
u
∂
t
-
a
(
∂
3
u
∂
x
3
+
6
u
∂
u
∂
x
)
-
2
b
∂
v
∂
x
=
0
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
(
x
,
t
)
∈
[
x
L
,
x
R
]
×
(
0
,
T
]
,
∂
v
∂
t
+
∂
3
v
∂
x
3
+
3
u
∂
v
∂
x
=
0
,
(
x
,
t
)
∈
[
x
L
,
x
R
]
×
(
0
,
T
]
,
subject to the initial conditions
(5)
u
(
x
,
0
)
=
g
(
x
)
,
v
(
x
,
0
)
=
0
and the boundary conditions
(6)
u
(
x
L
,
t
)
=
u
(
x
R
,
t
)
=
0
,
u
x
(
x
L
,
t
)
=
u
x
(
x
R
,
t
)
=
0
.
A standard weak formulation is obtained by multiplying (4) by a twice differentiable test function ψ
(
x
) and integrating by parts to obtain
(7)
(
u
t
,
ψ
)
-
6
a
(
u
u
x
,
ψ
)
-
a
(
u
x
,
ψ
x
x
)
-
2
b
(
v
v
x
,
ψ
)
=
0
,
(
v
t
,
ψ
)
+
3
(
u
v
x
,
ψ
)
+
(
v
x
,
ψ
x
x
)
=
0
,
where (
,
) denotes the usual L
2 inner product
(8)
(
f
,
g
)
=
∫
x
L
x
R
‍
f
(
x
)
g
(
x
)
d
x
.The space interval[
x
L
,
x
R
] is discretized by uniform (
N
+
1
) grid points. Consider
(9)
x
m
=
x
L
+
m
h
,
m
=
0,1
,
…
,
N
,
where the grid spacing h is given by h
=
(
x
R
-
x
L
)
/
N. We introduce finite elements in space in (7). Approximate the exact solutions of u
(
x
,
t
) and v
(
x
,
t
) by
(10)
u
(
x
,
t
)
=
∑
m
=
0
N
‍
U
m
(
t
)
ϕ
m
(
x
)
,
v
(
x
,
t
)
=
∑
m
=
0
N
‍
V
m
(
t
)
ϕ
m
(
x
)
,
where the trial functions ϕ
m
(
x
)
=
ϕ
(
(
x
-
x
m
)
/
h
), m
=
0,1
,
…
,
N, are the piecewise linear functions. Consider
(11)
ϕ
(
x
)
=
{
1
+
x
if
-
1
<
x
≤
0
1
-
x
if
0
<
x
≤
1
0
otherwise
.The unknown functionsU
m
(
t
), V
m
(
t
), m
=
0,1
,
2
,
…
,
N, are determined from the ordinary differential system. Consider
(12)
(
U
t
,
ψ
j
)
-
6
a
(
U
U
x
,
ψ
j
)
-
a
(
U
x
,
(
ψ
j
)
x
x
)
-
2
b
(
V
V
x
,
ψ
j
)
=
0
,
(
V
t
,
,
ψ
j
)
+
3
(
U
V
x
,
ψ
j
)
+
(
V
x
,
(
ψ
j
)
x
x
)
=
0
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
j
=
0,1
,
2
,
…
,
N
,
where the test functions ψ
j
(
x
)
=
ψ
(
(
x
-
x
i
)
/
h
), j
=
1,2
,
…
,
N, are the cubic B-spline with compact support. Consider
(13)
ψ
(
x
)
=
1
4
{
(
x
+
2
)
3
if
-
2
≤
x
≤
-
1
[
(
2
+
x
)
3
-
4
(
1
+
x
)
3
]
if
-
1
<
x
≤
0
[
(
2
-
x
)
3
-
4
(
1
-
x
)
3
]
if
0
<
x
≤
1
(
2
-
x
)
3
if
1
<
x
≤
2
0
otherwise
.Direct calculation of the inner products given in (12) will lead us to the following system of first order ordinary differential system (ODES):
(14)
1
80
(
U
˙
m
-
2
+
26
U
˙
m
-
1
+
66
U
˙
m
+
26
U
˙
m
+
1
+
U
˙
m
+
2
)
-
6
a
80
h
F
(
U
m
,
U
m
)
-
2
b
80
h
F
(
V
m
,
V
m
)
-
3
a
4
h
3
(
U
m
+
2
-
2
U
m
+
1
+
2
U
m
-
1
-
U
m
-
2
)
=
0
,
1
80
(
V
˙
m
-
2
+
26
V
˙
m
-
1
+
66
V
˙
m
+
26
V
˙
m
+
1
+
V
˙
m
+
2
)
+
3
4
h
3
(
V
m
+
2
-
2
V
m
+
1
+
2
V
m
-
1
-
V
m
-
2
)
+
3
80
h
F
(
U
m
,
V
m
)
=
0
,
where
(15)
F
(
U
m
,
V
m
)
=
(
-
V
m
-
2
+
V
m
-
1
)
(
U
m
-
2
+
4
U
m
-
1
)
+
(
-
V
m
-
1
+
V
m
)
(
22
U
m
-
1
+
33
U
m
)
+
(
-
V
m
+
V
m
+
1
)
(
33
U
m
+
22
U
m
+
1
)
+
(
-
V
m
+
1
+
V
m
+
2
)
(
4
U
m
+
1
+
U
m
+
2
)
.Now to the ODES, we assumeU
m
n and V
m
n to be the fully discrete approximation to the exact solutions u
(
x
m
,
t
n
) and v
(
x
m
,
t
n
), where t
n
=
n
k and k is the time step size. By using the finite difference approximation
(16)
U
˙
m
=
U
m
n
+
1
-
U
m
n
k
,
U
m
=
U
m
n
+
1
+
U
m
n
2
=
U
*
,
V
˙
m
=
V
m
n
+
1
-
V
m
n
k
,
V
m
=
V
m
n
+
1
+
V
m
n
2
=
V
*
,
the ODES (14) will be reduced to the implicit nonlinear finite difference scheme
(17)
(
U
m
-
2
n
+
1
+
26
U
m
-
1
n
+
1
+
66
U
m
n
+
1
+
26
U
m
+
1
n
+
1
+
U
m
+
2
n
+
1
)
-
(
U
m
-
2
n
+
26
U
m
-
1
n
+
66
U
m
n
+
26
U
m
+
1
n
+
U
m
+
2
n
)
-
p
1
(
U
m
+
2
*
-
2
U
m
+
1
*
+
2
U
m
-
1
*
-
U
m
-
2
*
)
+
p
2
F
(
U
*
,
U
*
)
+
p
3
F
(
V
*
,
V
*
)
=
0
,
(18)
(
V
m
-
2
n
+
1
+
26
V
m
-
1
n
+
1
+
66
V
m
n
+
1
+
26
V
m
+
1
n
+
1
+
V
m
+
2
n
+
1
)
-
(
V
m
-
2
n
+
26
V
m
-
1
n
+
66
V
m
n
+
26
V
m
+
1
n
+
V
m
+
2
n
)
+
p
4
(
V
m
+
2
*
-
2
V
m
+
1
*
+
2
V
m
-
1
*
-
V
m
-
2
*
)
+
p
5
F
(
U
*
,
V
*
)
=
0
,
m
=
1,2
,
…
,
N
,
(19)
p
1
=
-
60
a
k
h
3
,
p
2
=
-
6
a
k
h
,
p
3
=
-
2
b
k
h
,
p
4
=
60
k
h
3
,
p
5
=
3
k
h
.The system (17)-(18) is a nonlinear block penta diagonal system in the unknown vectors U
n
+
1 and V
n
+
1 and is solved using Newton’s method. We denote this method by Scheme 1.


## 3. Analysis of the Method

In this section we will discuss the stability and the accuracy of Scheme 1.


### 3.1. Stability of the Scheme

To study the stability of Scheme 1, von Neumann stability analysis will be used and this can be only applied for linear finite difference schemes; accordingly, we consider the linear version of the proposed method(20)
(
U
m
-
2
n
+
1
+
26
U
m
-
1
n
+
1
+
66
U
m
n
+
1
+
26
U
m
+
1
n
+
1
+
U
m
+
2
n
+
1
)
-
(
U
m
-
2
n
+
26
U
m
-
1
n
+
66
U
m
n
+
26
U
m
+
1
n
+
U
m
+
2
n
)
+
p
1
(
U
m
+
2
*
-
2
U
m
+
1
*
+
2
U
m
-
1
*
-
U
m
-
2
*
)
+
p
2
F
(
U
¯
,
U
*
)
+
p
3
F
(
V
¯
,
V
*
)
=
0
,
(
V
m
-
2
n
+
1
+
26
V
m
-
1
n
+
1
+
66
V
m
n
+
1
+
26
V
m
+
1
n
+
1
+
V
m
+
2
n
+
1
)
-
(
V
m
-
2
n
+
26
V
m
-
1
n
+
66
V
m
n
+
26
V
m
+
1
n
+
V
m
+
2
n
)
+
p
4
(
V
m
+
2
*
-
2
V
m
+
1
*
+
2
V
m
-
1
*
-
V
m
-
2
*
)
+
p
5
F
(
U
¯
,
V
*
)
=
0
,
where
(21)
F
(
U
¯
,
V
)
=
U
¯
[
V
m
+
2
+
10
V
m
+
1
-
10
V
m
-
1
-
V
m
-
2
]
,p
2
′
=
5
p
2, p
3
′
=
5
p
3, and p
5
′
=
5
p
5. U
¯ and V
¯ are assumed to be constant on the whole range. We assume the solution of the linearized scheme (20) to be of the form
(22)
U
m
n
=
U
n
e
i
β
m
h
,
V
m
n
=
V
n
e
i
β
m
h
,
where β is a real constant. Direct substitution of (22) into (20) will lead us to the following system:
(23)
[
γ
1
+
i
p
1
γ
2
+
i
p
2
′
U
¯
γ
3
]
U
n
+
1
+
i
p
3
′
V
¯
γ
3
V
n
+
1
=
[
γ
1
-
i
p
1
γ
2
-
i
p
2
′
U
¯
γ
3
]
U
n
-
i
p
3
′
V
¯
γ
3
V
n
,
[
γ
1
+
i
p
4
γ
2
+
i
p
5
′
U
¯
γ
3
]
V
n
+
1
=
[
γ
1
-
i
p
4
γ
2
-
i
p
5
′
U
¯
γ
3
]
V
n
,
where
(24)
γ
1
=
66
+
2
cos
⁡
(
2
β
h
)
+
52
cos
⁡
(
β
h
)
,
γ
2
=
2
sin
(
2
β
h
)
-
4
sin
(
β
h
)
,
γ
3
=
2
sin
(
2
β
h
)
+
20
sin
(
β
h
)
.The system (23) can be written in a matrix vector form as
(25)
Ψ
n
+
1
=
B
Ψ
n
,
where Ψ
=
[
U
V
]
t and B is the (
2
×
2
) matrix. Consider
(26)
B
=
[
γ
1
+
i
c
1
c
2
0
γ
1
+
i
c
3
]
-
1
[
γ
1
-
i
c
1
-
c
2
0
γ
1
-
i
c
3
]
,
where
(27)
c
1
=
p
1
γ
2
+
p
2
′
U
¯
γ
3
,
c
2
=
p
3
′
V
¯
γ
3
,
c
3
=
p
4
γ
2
+
p
5
′
U
¯
γ
3
.
von Neumann stability condition for the system (25) states that, for the scheme to be stable, the maximum modulus of the eigenvalues of the amplification matrix B should be less than or equal to one. The eigenvalues of the matrix B are
(28)
λ
1
=
γ
1
-
i
c
1
γ
1
+
i
c
1
,
λ
2
=
γ
1
-
i
c
3
γ
1
+
i
c
3
,
which has modulus equal to one, and hence the scheme under consideration is unconditionally stable in the linearized sense.


### 3.2.  Accuracy of the Scheme

In order to study the accuracy of the scheme, we replace the numerical solutionU
m
n and V
m
n in (17) and (18) by the exact solutions u
m
n and v
m
n. By making use of the following Taylor series expansions at the point (
x
m
,
t
n
)
(29)
1
k
[
(
u
m
-
2
n
+
1
+
26
u
m
-
1
n
+
1
+
66
u
m
n
+
1
+
26
u
m
+
1
n
+
1
+
u
m
+
2
n
+
1
)
-
(
u
m
-
2
n
+
26
u
m
-
1
n
+
66
u
m
n
+
26
u
m
+
1
n
+
u
m
+
2
n
)
26
u
m
+
1
n
+
1
]
=
3
2
u
t
+
3
k
4
u
t
t
+
3
h
2
8
u
x
x
t
+
k
2
4
u
t
t
t
+
⋯
(30)
3
8
h
3
[
u
m
+
2
n
+
1
-
2
u
m
+
1
n
+
1
+
2
u
m
-
1
n
+
1
-
u
m
-
2
n
+
1
u
m
+
2
n
+
1
-
2
u
m
+
1
n
+
1
+
2
u
m
-
1
n
+
1
-
u
m
-
2
n
+
1
+
u
m
+
2
n
-
2
u
m
+
1
n
+
2
u
m
-
1
n
-
u
m
-
2
n
]
=
3
2
u
x
x
x
+
3
k
4
u
x
x
x
t
+
3
h
2
8
u
x
x
x
x
x
+
3
k
2
8
u
x
x
x
t
t
+
⋯
,
(31)
1
80
h
F
(
u
*
,
u
*
)
=
3
2
u
u
x
+
3
k
4
(
u
u
x
)
t
+
h
2
[
3
8
u
u
x
x
x
+
7
8
u
x
u
x
x
]
+
3
k
2
8
[
u
x
u
t
t
+
u
t
u
x
t
+
u
u
x
t
t
]
+
⋯
,
(32)
1
80
h
F
(
u
*
,
v
*
)
=
3
2
u
v
x
+
3
k
4
(
u
v
x
)
t
+
h
2
[
3
8
u
v
x
x
x
+
3
8
v
x
u
x
x
+
1
2
u
x
v
x
x
]
+
3
k
2
8
[
v
x
u
t
t
+
u
t
v
x
t
+
u
v
x
t
t
]
+
⋯
,
and the CKdV system (1), we obtain the local truncation error (LTE) for the proposed scheme
(33)
LTE
(
17
)
=
-
3
h
2
2
u
x
u
x
x
-
b
h
2
2
v
u
x
x
-
k
2
8
u
t
t
t
+
9
k
2
a
4
u
t
u
x
t
-
3
b
k
2
4
v
t
v
x
t
+
⋯
,
LTE
(
18
)
=
-
3
h
2
2
u
x
v
x
x
-
k
2
8
u
t
t
t
-
3
k
2
8
u
t
v
x
t
+
⋯
,
and hence the scheme is of second order in both directions space and time.


## 3.1. Stability of the Scheme

To study the stability of Scheme 1, von Neumann stability analysis will be used and this can be only applied for linear finite difference schemes; accordingly, we consider the linear version of the proposed method(20)
(
U
m
-
2
n
+
1
+
26
U
m
-
1
n
+
1
+
66
U
m
n
+
1
+
26
U
m
+
1
n
+
1
+
U
m
+
2
n
+
1
)
-
(
U
m
-
2
n
+
26
U
m
-
1
n
+
66
U
m
n
+
26
U
m
+
1
n
+
U
m
+
2
n
)
+
p
1
(
U
m
+
2
*
-
2
U
m
+
1
*
+
2
U
m
-
1
*
-
U
m
-
2
*
)
+
p
2
F
(
U
¯
,
U
*
)
+
p
3
F
(
V
¯
,
V
*
)
=
0
,
(
V
m
-
2
n
+
1
+
26
V
m
-
1
n
+
1
+
66
V
m
n
+
1
+
26
V
m
+
1
n
+
1
+
V
m
+
2
n
+
1
)
-
(
V
m
-
2
n
+
26
V
m
-
1
n
+
66
V
m
n
+
26
V
m
+
1
n
+
V
m
+
2
n
)
+
p
4
(
V
m
+
2
*
-
2
V
m
+
1
*
+
2
V
m
-
1
*
-
V
m
-
2
*
)
+
p
5
F
(
U
¯
,
V
*
)
=
0
,
where
(21)
F
(
U
¯
,
V
)
=
U
¯
[
V
m
+
2
+
10
V
m
+
1
-
10
V
m
-
1
-
V
m
-
2
]
,p
2
′
=
5
p
2, p
3
′
=
5
p
3, and p
5
′
=
5
p
5. U
¯ and V
¯ are assumed to be constant on the whole range. We assume the solution of the linearized scheme (20) to be of the form
(22)
U
m
n
=
U
n
e
i
β
m
h
,
V
m
n
=
V
n
e
i
β
m
h
,
where β is a real constant. Direct substitution of (22) into (20) will lead us to the following system:
(23)
[
γ
1
+
i
p
1
γ
2
+
i
p
2
′
U
¯
γ
3
]
U
n
+
1
+
i
p
3
′
V
¯
γ
3
V
n
+
1
=
[
γ
1
-
i
p
1
γ
2
-
i
p
2
′
U
¯
γ
3
]
U
n
-
i
p
3
′
V
¯
γ
3
V
n
,
[
γ
1
+
i
p
4
γ
2
+
i
p
5
′
U
¯
γ
3
]
V
n
+
1
=
[
γ
1
-
i
p
4
γ
2
-
i
p
5
′
U
¯
γ
3
]
V
n
,
where
(24)
γ
1
=
66
+
2
cos
⁡
(
2
β
h
)
+
52
cos
⁡
(
β
h
)
,
γ
2
=
2
sin
(
2
β
h
)
-
4
sin
(
β
h
)
,
γ
3
=
2
sin
(
2
β
h
)
+
20
sin
(
β
h
)
.The system (23) can be written in a matrix vector form as
(25)
Ψ
n
+
1
=
B
Ψ
n
,
where Ψ
=
[
U
V
]
t and B is the (
2
×
2
) matrix. Consider
(26)
B
=
[
γ
1
+
i
c
1
c
2
0
γ
1
+
i
c
3
]
-
1
[
γ
1
-
i
c
1
-
c
2
0
γ
1
-
i
c
3
]
,
where
(27)
c
1
=
p
1
γ
2
+
p
2
′
U
¯
γ
3
,
c
2
=
p
3
′
V
¯
γ
3
,
c
3
=
p
4
γ
2
+
p
5
′
U
¯
γ
3
.
von Neumann stability condition for the system (25) states that, for the scheme to be stable, the maximum modulus of the eigenvalues of the amplification matrix B should be less than or equal to one. The eigenvalues of the matrix B are
(28)
λ
1
=
γ
1
-
i
c
1
γ
1
+
i
c
1
,
λ
2
=
γ
1
-
i
c
3
γ
1
+
i
c
3
,
which has modulus equal to one, and hence the scheme under consideration is unconditionally stable in the linearized sense.


## 3.2.  Accuracy of the Scheme

In order to study the accuracy of the scheme, we replace the numerical solutionU
m
n and V
m
n in (17) and (18) by the exact solutions u
m
n and v
m
n. By making use of the following Taylor series expansions at the point (
x
m
,
t
n
)
(29)
1
k
[
(
u
m
-
2
n
+
1
+
26
u
m
-
1
n
+
1
+
66
u
m
n
+
1
+
26
u
m
+
1
n
+
1
+
u
m
+
2
n
+
1
)
-
(
u
m
-
2
n
+
26
u
m
-
1
n
+
66
u
m
n
+
26
u
m
+
1
n
+
u
m
+
2
n
)
26
u
m
+
1
n
+
1
]
=
3
2
u
t
+
3
k
4
u
t
t
+
3
h
2
8
u
x
x
t
+
k
2
4
u
t
t
t
+
⋯
(30)
3
8
h
3
[
u
m
+
2
n
+
1
-
2
u
m
+
1
n
+
1
+
2
u
m
-
1
n
+
1
-
u
m
-
2
n
+
1
u
m
+
2
n
+
1
-
2
u
m
+
1
n
+
1
+
2
u
m
-
1
n
+
1
-
u
m
-
2
n
+
1
+
u
m
+
2
n
-
2
u
m
+
1
n
+
2
u
m
-
1
n
-
u
m
-
2
n
]
=
3
2
u
x
x
x
+
3
k
4
u
x
x
x
t
+
3
h
2
8
u
x
x
x
x
x
+
3
k
2
8
u
x
x
x
t
t
+
⋯
,
(31)
1
80
h
F
(
u
*
,
u
*
)
=
3
2
u
u
x
+
3
k
4
(
u
u
x
)
t
+
h
2
[
3
8
u
u
x
x
x
+
7
8
u
x
u
x
x
]
+
3
k
2
8
[
u
x
u
t
t
+
u
t
u
x
t
+
u
u
x
t
t
]
+
⋯
,
(32)
1
80
h
F
(
u
*
,
v
*
)
=
3
2
u
v
x
+
3
k
4
(
u
v
x
)
t
+
h
2
[
3
8
u
v
x
x
x
+
3
8
v
x
u
x
x
+
1
2
u
x
v
x
x
]
+
3
k
2
8
[
v
x
u
t
t
+
u
t
v
x
t
+
u
v
x
t
t
]
+
⋯
,
and the CKdV system (1), we obtain the local truncation error (LTE) for the proposed scheme
(33)
LTE
(
17
)
=
-
3
h
2
2
u
x
u
x
x
-
b
h
2
2
v
u
x
x
-
k
2
8
u
t
t
t
+
9
k
2
a
4
u
t
u
x
t
-
3
b
k
2
4
v
t
v
x
t
+
⋯
,
LTE
(
18
)
=
-
3
h
2
2
u
x
v
x
x
-
k
2
8
u
t
t
t
-
3
k
2
8
u
t
v
x
t
+
⋯
,
and hence the scheme is of second order in both directions space and time.


## 4. Product Approximation Technique

A modified Petrov-Galerkin method for solving the CKdV system (1) can be achieved by using the product approximation technique, where we used special approximation to the nonlinear terms in the differential system. In order to apply this technique, we rewrite the CKdV in the following form:
(34)
∂
u
∂
t
=
a
(
∂
3
u
∂
x
3
+
3
∂
(
u
2
)
∂
x
)
+
b
∂
(
v
2
)
∂
x
,
(35)
∂
v
∂
t
=
-
∂
3
v
∂
x
3
-
3
u
∂
v
∂
x
.The product approximation technique is used to approximate the nonlinear terms in (34) in the following manner:
(36)
u
2
(
x
,
t
)
=
∑
m
=
1
N
‍
U
m
2
(
t
)
ϕ
m
(
x
)
,
v
2
(
x
,
t
)
=
∑
m
=
1
N
‍
V
m
2
(
t
)
ϕ
m
(
x
)
.By using the same procedure in deriving Scheme 1 and the approximation (37), we can obtain after some manipulations the following scheme:
(37)
p
1
[
(
U
m
-
2
n
+
1
+
26
U
m
-
1
n
+
1
+
66
U
m
n
+
1
+
26
U
m
+
1
n
+
1
+
U
m
+
2
n
+
1
)
-
(
U
m
-
2
n
+
26
U
m
-
1
n
+
66
U
m
n
+
26
U
m
+
1
n
+
U
m
+
2
n
)
26
U
m
+
1
n
+
1
]
-
p
2
[
(
U
m
+
2
*
)
2
+
10
(
U
m
+
1
*
)
2
-
10
(
U
m
-
1
*
)
2
-
(
U
m
-
2
*
)
2
]
-
p
3
(
U
m
+
2
*
-
2
U
m
+
1
*
+
2
U
m
-
1
*
-
U
m
-
2
*
)
-
p
4
[
(
V
m
+
2
*
)
2
+
10
(
V
m
+
1
*
)
2
-
10
(
V
m
-
1
*
)
2
-
(
V
m
-
2
*
)
2
]
=
0
,
(38)
p
1
[
(
V
m
-
2
n
+
1
+
26
V
m
-
1
n
+
1
+
66
V
m
n
+
1
+
26
V
m
+
1
n
+
1
+
V
m
+
2
n
+
1
)
-
(
V
m
-
2
n
+
26
V
m
-
1
n
+
66
V
m
n
+
26
V
m
+
1
n
+
V
m
+
2
n
)
26
U
m
+
1
n
+
1
]
+
p
5
(
V
m
+
2
*
-
2
V
m
+
1
*
+
2
V
m
-
1
*
-
V
m
-
2
*
)
+
p
6
[
-
(
U
m
-
2
*
+
4
U
m
-
1
*
)
V
m
-
2
*
l
l
l
l
l
l
l
+
(
U
m
-
2
*
-
18
U
m
-
1
*
-
33
U
m
*
)
V
m
-
1
*
l
l
l
l
l
l
l
+
22
(
U
m
-
1
*
-
U
m
+
1
*
)
V
m
*
l
l
l
l
l
l
l
+
(
33
U
m
*
+
18
U
m
+
1
*
-
U
m
+
2
*
)
V
m
+
1
*
l
l
l
l
l
l
+
(
4
U
m
+
1
*
+
U
m
+
2
*
)
V
m
+
2
*
]
=
0
,
where
(39)
U
m
*
=
U
m
n
+
1
+
U
m
n
,
V
m
*
=
V
m
n
+
1
+
V
m
n
,
p
1
=
1
80
,
p
2
=
3
a
k
64
h
,
p
3
=
3
a
k
8
h
3
,
p
4
=
b
k
64
h
,
p
5
=
3
k
8
h
3
,
p
6
=
3
k
320
h
.Again, the method is unconditionally stable according to von Neumann stability analysis, and it is of second order in both directions. The scheme produced a nonlinear block pentadiagonal system and its solution obtained using Newton’s method. We have noticed that the accuracy has been improved in the first equation (38) as we will see in the next section. We will denote the scheme obtained by using the product approximation technique by Scheme 2.


## 5. Numerical Results

To gain insight into the performance of the proposed schemes, we perform different numerical tests, like single soliton, two and three solitons interaction, and birth of solitons. The conservation properties of the proposed schemes are examined by calculatingI
1
,
I
2, and I
3 using the trapezoidal rule. The L
∞
(
u
) and L
∞
(
v
) error norms are defined as
(40)
L
∞
(
u
)
=
∥
U
n
-
u
n
∥
∞
=
max
⁡
m
|
U
m
n
-
u
m
n
|
,
L
∞
(
v
)
=
∥
V
n
-
v
n
∥
∞
=
max
⁡
m
|
V
m
n
-
v
m
n
|
.
are used to examine the accuracy of the proposed schemes.


### 5.1. Single Soliton

In the first test, we choose the initial conditions(41)
u
(
x
,
0
)
=
2
λ
2
sec
h
2
(
ξ
)
,
v
(
x
,
0
)
=
1
2
w
sech
(
ξ
)
,
ξ
=
λ
x
+
1
2
log
⁡
(
w
)
,
w
=
-
b
8
(
4
a
+
1
)
λ
4
,
which represents a single soliton solution at t
=
0. To study the behavior of numerical solution using  Scheme 1 and Scheme 2, we choose the set of parameters as h
=
0.05, k
=
0.01, tol
=
1
0
-
7, x
L
=
-
25, x
R
=
25, a
=
-
0.125, b
=
-
3, and λ
=
0.5. The conserved quantities and the error norms L
∞
(
U
), L
∞
(
V
) are displayed in Tables 1 and 2 for Scheme 1 and Scheme 2, respectively. It is clear from these tables that our schemes are highly accurate. In addition, the schemes preserve the conserved quantities exactly during the evolution of the numerical solution from t
=
0 to t
=
10. The execution time required to produce Table 1 is 2.328 second and 2.171 second to produce Table 2. We have noticed that Scheme 2 has an upper hand over Scheme 1 with respect to accuracy and CPU time. In Figures 1 and 2, we display the numerical solution of U
m
n and V
m
n for t
=
0,1
,
2
,
…
,
20.Table 1
Petrov-Galerkin method.
T
L
∞
(
U
)
L
∞
(
V
)
I
1
I
2
I
3
0
0.000000
0.000000
2.000000
0.5000000
0.112529
2
0.000015
0.000004
2.000000
0.5000000
0.112529
4
0.000021
0.000005
2.000000
0.5000000
0.112529
6
0.000023
0.000006
2.000000
0.5000000
0.112529
8
0.000024
0.000007
2.000000
0.5000000
0.112529
10
0.000025
0.000008
2.000000
0.5000000
0.112529Table 2
Product approximation technique.
T
L
∞
(
U
)
L
∞
(
V
)
I
1
I
2
I
3
0
0.000000
0.000000
2.000000
0.5000000
0.112529
2
0.000004
0.000005
2.000000
0.5000000
0.112529
4
0.000008
0.000007
2.000000
0.5000000
0.112529
6
0.000010
0.000009
2.000000
0.5000000
0.112529
8
0.000013
0.000012
2.000000
0.5000000
0.112529
10
0.000014
0.000013
2.000000
0.5000000
0.112529Figure 1
Single solitonU
:
a
=
-
0
.
125, b
=
-
3
.
0, and λ
=
0
.
5.Figure 2
Single solitonV
:
a
=
-
0
.
125, b
=
-
3
.
0, and λ
=
0
.
5.By choosing the set of valuesk
=
0.01, a
=
0.5, b
=
-
3.0, λ
=
0.5, and t
=
1, we perform a comparison of Scheme 1 and Scheme 2 with Ismail [6] and we display this in Table 3, we can easily see that the three methods produce highly accurate results with some credits for collocation method.Table 3
Comparison of three methods witha
=
0.5, b
=
-
3.0, λ
=
0.5, k
=
0.01, and t
=
1.
h
Petrov-Galerkin
Product approximation
Collocation (Ismail [6])
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
0.1
0.000051
0.000027
0.218
0.000014
0.000019
0.203
0.000000
0.000003
0.188
0.05
0.000013
0.000007
0.343
0.000009
0.000008
0.328
0.000000
0.000003
0.312
0.025
0.000003
0.000004
0.625
0.000002
0.000004
0.609
0.000000
0.000003
0.578


### 5.2.  Two Solitons Interaction

To study the interaction of two solitons, we choose the initial conditions as(42)
u
(
x
,
0
)
=
∑
j
=
1
2
u
j
(
x
,
0
)
,
v
(
x
,
0
)
=
∑
j
=
1
2
v
j
(
x
,
0
)
,
where
(43)
u
j
(
x
,
0
)
=
2
λ
j
2
sec
h
2
(
ξ
j
)
,
v
j
(
x
,
0
)
=
1
2
w
j
sech
(
ξ
j
)
,
ξ
j
=
λ
j
(
x
-
y
j
)
+
1
2
log
⁡
⁡
(
w
j
)
,
w
j
=
-
b
8
(
4
a
+
1
)
λ
j
4
,
j
=
1,2
,
which represents the sum of two single solitons, we assign the value of the parameters c
=
0, d
=
100, h
=
0.05, k
=
0.01, a
=
0.5, b
=
-
3.0, λ
1
=
1.0, λ
2
=
0.6, y
1
=
10, and y
2
=
30. In Table 4, we present the conserved quantities during the interaction scenario and show that our numerical methods achieved the goal of conserving these quantities. The interaction scenario is presented in Figures 3 and 4. The contours of the interaction process are given in Figure 5. We have noticed that the taller (faster) wave collides with the shorter (slower) wave and leaves the interaction region without any disturbance in their identities. This phenomenon indicates the interaction scenario is elastic [1].Table 4
Two solitons interaction (conserved quantities).
T
I
1
I
2
I
3
0.00
6.400000
−3.243013
−2.584118
20.0
6.399985
−3.243298
−2.586545
40.0
6.400034
−3.243016
−2.584155
60.0
6.399962
−3.243017
−2.584203
80.0
6.400234
−3.242978
−2.583492Figure 3
Numerical solutionU
n with parameters λ
1
=
1.0, λ
2
=
0
.
6, y
1
=
10, and y
2
=
30.Figure 4
Numerical solutionV
n with parameters λ
1
=
1
.
0, λ
2
=
0
.
6, y
1
=
10, and y
2
=
30.Figure 5
Contours of the numerical solutionU
n (interaction of two solitons).To examine the interaction scenario fora
≠
1
/
2, we use the set of parameters a
=
0.495, b
=
-
3, λ
1
=
1.0, λ
2
=
0.6, y
1
=
10, y
2
=
30, c
=
0, d
=
100, h
=
0.05, and k
=
0.01. Our results are demonstrated in Figures 6 and 7; the amplitudes of the two solitons have changed after the interaction, which indicates that the interaction scenario is inelastic. We have tested other values of a; in all cases we have found that the interaction is inelastic and this in agreement with [1], which they claim the CKdV equation, is integrable only for a
=
1
/
2.Figure 6
Inelastic interaction: numerical solutionU
n with a
=
0
.
495.Figure 7
Inelastic interaction: numerical solutionV
n with a
=
0
.
495.


### 5.3. Three Solitons Interaction

To study the interaction of three solitons, we choose the initial condition as a sum of three well separated single solitons in the following form:(44)
u
(
x
,
0
)
=
∑
j
=
1
3
u
j
(
x
,
0
)
,
v
(
x
,
0
)
=
∑
j
=
1
3
v
j
(
x
,
0
)
,
where
(45)
u
j
(
x
,
0
)
=
2
λ
j
2
sec
h
2
(
ξ
j
)
,
v
j
(
x
,
0
)
=
1
2
w
j
sech
(
ξ
j
)
,
ξ
j
=
λ
j
(
x
-
y
j
)
+
1
2
log
⁡
⁡
(
w
j
)
,
w
j
=
-
b
8
(
4
a
+
1
)
λ
j
4
,
j
=
1,2
,
3
.In this test we choose the set of parametersx
L
=
0, x
R
=
80, h
=
0.05, k
=
0.01, a
=
0.5, b
=
-
3.0, λ
1
=
1.0, λ
2
=
0.6, λ
3
=
0.3, y
1
=
10, y
2
=
30, and y
3
=
50.The simulation of the three solitons interaction scenario is given in Figures8 and 9, respectively, for the time duration t
=
0,4
,
…
,
80. It is very clear to see how the soliton with the largest amplitude u
1 interacts with the other two solitons and leaves the interaction region unchanged in shape. The three solitons appear after the interaction in the reverse order compared to the initial state. In Table 5, we display the conserved quantities during the interaction scenario.Table 5
Three solitons interaction: the conserved quantities.
T
I
1
I
2
I
3
0
7.600000
−3.351855
−2.610283
20
7.600042
−3.352227
−2.613149
40
7.600800
−3.352034
−2.611347
60
7.599367
−3.351547
−2.608970
80
7.600545
−3.351562
−2.602090Figure 8
Interaction of three solitons: numerical solutionU
n with parameters λ
1
=
1, λ
2
=
0
.
6, and λ
3
=
0.3.Figure 9
Interaction of three solitons: numerical solutionV
n with parameters λ
1
=
1, λ
2
=
0
.
6, and λ
3
=
0
.
3.


### 5.4.  Birth of Solitons

In this test we choose the initial condition(46)
u
(
x
,
0
)
=
exp
⁡
(
-
0.01
x
2
)
,
v
(
x
,
0
)
=
exp
⁡
(
-
0.01
x
2
)
,
with the following set of parameters a
=
0.5, b
=
-
3.0, 0
≤
x
≤
100, h
=
0.05, and k
=
0.01. We have noticed as time evolves, a birth of four solitons with different amplitudes and this can be easily seen in Figure 10. The conserved quantities are given in Table 6 which is almost conserved.Table 6
Birth of solitons witha
=
0
.
5 and b
=
-
3
.
0.
T
I
1
I
2
0.0
17.724343
−12.533142
5.0
17.723816
−12.532956
10.0
17.723352
−12.530116
15.0
17.722782
−12.529239
20.0
17.722217
−12.529013
25.0
17.721734
−12.528983Figure 10
Birth of solitons with parametersa
=
0
.
5 and b
=
-
3
.
0.


## 5.1. Single Soliton

In the first test, we choose the initial conditions(41)
u
(
x
,
0
)
=
2
λ
2
sec
h
2
(
ξ
)
,
v
(
x
,
0
)
=
1
2
w
sech
(
ξ
)
,
ξ
=
λ
x
+
1
2
log
⁡
(
w
)
,
w
=
-
b
8
(
4
a
+
1
)
λ
4
,
which represents a single soliton solution at t
=
0. To study the behavior of numerical solution using  Scheme 1 and Scheme 2, we choose the set of parameters as h
=
0.05, k
=
0.01, tol
=
1
0
-
7, x
L
=
-
25, x
R
=
25, a
=
-
0.125, b
=
-
3, and λ
=
0.5. The conserved quantities and the error norms L
∞
(
U
), L
∞
(
V
) are displayed in Tables 1 and 2 for Scheme 1 and Scheme 2, respectively. It is clear from these tables that our schemes are highly accurate. In addition, the schemes preserve the conserved quantities exactly during the evolution of the numerical solution from t
=
0 to t
=
10. The execution time required to produce Table 1 is 2.328 second and 2.171 second to produce Table 2. We have noticed that Scheme 2 has an upper hand over Scheme 1 with respect to accuracy and CPU time. In Figures 1 and 2, we display the numerical solution of U
m
n and V
m
n for t
=
0,1
,
2
,
…
,
20.Table 1
Petrov-Galerkin method.
T
L
∞
(
U
)
L
∞
(
V
)
I
1
I
2
I
3
0
0.000000
0.000000
2.000000
0.5000000
0.112529
2
0.000015
0.000004
2.000000
0.5000000
0.112529
4
0.000021
0.000005
2.000000
0.5000000
0.112529
6
0.000023
0.000006
2.000000
0.5000000
0.112529
8
0.000024
0.000007
2.000000
0.5000000
0.112529
10
0.000025
0.000008
2.000000
0.5000000
0.112529Table 2
Product approximation technique.
T
L
∞
(
U
)
L
∞
(
V
)
I
1
I
2
I
3
0
0.000000
0.000000
2.000000
0.5000000
0.112529
2
0.000004
0.000005
2.000000
0.5000000
0.112529
4
0.000008
0.000007
2.000000
0.5000000
0.112529
6
0.000010
0.000009
2.000000
0.5000000
0.112529
8
0.000013
0.000012
2.000000
0.5000000
0.112529
10
0.000014
0.000013
2.000000
0.5000000
0.112529Figure 1
Single solitonU
:
a
=
-
0
.
125, b
=
-
3
.
0, and λ
=
0
.
5.Figure 2
Single solitonV
:
a
=
-
0
.
125, b
=
-
3
.
0, and λ
=
0
.
5.By choosing the set of valuesk
=
0.01, a
=
0.5, b
=
-
3.0, λ
=
0.5, and t
=
1, we perform a comparison of Scheme 1 and Scheme 2 with Ismail [6] and we display this in Table 3, we can easily see that the three methods produce highly accurate results with some credits for collocation method.Table 3
Comparison of three methods witha
=
0.5, b
=
-
3.0, λ
=
0.5, k
=
0.01, and t
=
1.
h
Petrov-Galerkin
Product approximation
Collocation (Ismail [6])
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
L
∞
(
U
m
n
)
L
∞
(
V
m
n
)
CPU
0.1
0.000051
0.000027
0.218
0.000014
0.000019
0.203
0.000000
0.000003
0.188
0.05
0.000013
0.000007
0.343
0.000009
0.000008
0.328
0.000000
0.000003
0.312
0.025
0.000003
0.000004
0.625
0.000002
0.000004
0.609
0.000000
0.000003
0.578


## 5.2.  Two Solitons Interaction

To study the interaction of two solitons, we choose the initial conditions as(42)
u
(
x
,
0
)
=
∑
j
=
1
2
u
j
(
x
,
0
)
,
v
(
x
,
0
)
=
∑
j
=
1
2
v
j
(
x
,
0
)
,
where
(43)
u
j
(
x
,
0
)
=
2
λ
j
2
sec
h
2
(
ξ
j
)
,
v
j
(
x
,
0
)
=
1
2
w
j
sech
(
ξ
j
)
,
ξ
j
=
λ
j
(
x
-
y
j
)
+
1
2
log
⁡
⁡
(
w
j
)
,
w
j
=
-
b
8
(
4
a
+
1
)
λ
j
4
,
j
=
1,2
,
which represents the sum of two single solitons, we assign the value of the parameters c
=
0, d
=
100, h
=
0.05, k
=
0.01, a
=
0.5, b
=
-
3.0, λ
1
=
1.0, λ
2
=
0.6, y
1
=
10, and y
2
=
30. In Table 4, we present the conserved quantities during the interaction scenario and show that our numerical methods achieved the goal of conserving these quantities. The interaction scenario is presented in Figures 3 and 4. The contours of the interaction process are given in Figure 5. We have noticed that the taller (faster) wave collides with the shorter (slower) wave and leaves the interaction region without any disturbance in their identities. This phenomenon indicates the interaction scenario is elastic [1].Table 4
Two solitons interaction (conserved quantities).
T
I
1
I
2
I
3
0.00
6.400000
−3.243013
−2.584118
20.0
6.399985
−3.243298
−2.586545
40.0
6.400034
−3.243016
−2.584155
60.0
6.399962
−3.243017
−2.584203
80.0
6.400234
−3.242978
−2.583492Figure 3
Numerical solutionU
n with parameters λ
1
=
1.0, λ
2
=
0
.
6, y
1
=
10, and y
2
=
30.Figure 4
Numerical solutionV
n with parameters λ
1
=
1
.
0, λ
2
=
0
.
6, y
1
=
10, and y
2
=
30.Figure 5
Contours of the numerical solutionU
n (interaction of two solitons).To examine the interaction scenario fora
≠
1
/
2, we use the set of parameters a
=
0.495, b
=
-
3, λ
1
=
1.0, λ
2
=
0.6, y
1
=
10, y
2
=
30, c
=
0, d
=
100, h
=
0.05, and k
=
0.01. Our results are demonstrated in Figures 6 and 7; the amplitudes of the two solitons have changed after the interaction, which indicates that the interaction scenario is inelastic. We have tested other values of a; in all cases we have found that the interaction is inelastic and this in agreement with [1], which they claim the CKdV equation, is integrable only for a
=
1
/
2.Figure 6
Inelastic interaction: numerical solutionU
n with a
=
0
.
495.Figure 7
Inelastic interaction: numerical solutionV
n with a
=
0
.
495.


## 5.3. Three Solitons Interaction

To study the interaction of three solitons, we choose the initial condition as a sum of three well separated single solitons in the following form:(44)
u
(
x
,
0
)
=
∑
j
=
1
3
u
j
(
x
,
0
)
,
v
(
x
,
0
)
=
∑
j
=
1
3
v
j
(
x
,
0
)
,
where
(45)
u
j
(
x
,
0
)
=
2
λ
j
2
sec
h
2
(
ξ
j
)
,
v
j
(
x
,
0
)
=
1
2
w
j
sech
(
ξ
j
)
,
ξ
j
=
λ
j
(
x
-
y
j
)
+
1
2
log
⁡
⁡
(
w
j
)
,
w
j
=
-
b
8
(
4
a
+
1
)
λ
j
4
,
j
=
1,2
,
3
.In this test we choose the set of parametersx
L
=
0, x
R
=
80, h
=
0.05, k
=
0.01, a
=
0.5, b
=
-
3.0, λ
1
=
1.0, λ
2
=
0.6, λ
3
=
0.3, y
1
=
10, y
2
=
30, and y
3
=
50.The simulation of the three solitons interaction scenario is given in Figures8 and 9, respectively, for the time duration t
=
0,4
,
…
,
80. It is very clear to see how the soliton with the largest amplitude u
1 interacts with the other two solitons and leaves the interaction region unchanged in shape. The three solitons appear after the interaction in the reverse order compared to the initial state. In Table 5, we display the conserved quantities during the interaction scenario.Table 5
Three solitons interaction: the conserved quantities.
T
I
1
I
2
I
3
0
7.600000
−3.351855
−2.610283
20
7.600042
−3.352227
−2.613149
40
7.600800
−3.352034
−2.611347
60
7.599367
−3.351547
−2.608970
80
7.600545
−3.351562
−2.602090Figure 8
Interaction of three solitons: numerical solutionU
n with parameters λ
1
=
1, λ
2
=
0
.
6, and λ
3
=
0.3.Figure 9
Interaction of three solitons: numerical solutionV
n with parameters λ
1
=
1, λ
2
=
0
.
6, and λ
3
=
0
.
3.


## 5.4.  Birth of Solitons

In this test we choose the initial condition(46)
u
(
x
,
0
)
=
exp
⁡
(
-
0.01
x
2
)
,
v
(
x
,
0
)
=
exp
⁡
(
-
0.01
x
2
)
,
with the following set of parameters a
=
0.5, b
=
-
3.0, 0
≤
x
≤
100, h
=
0.05, and k
=
0.01. We have noticed as time evolves, a birth of four solitons with different amplitudes and this can be easily seen in Figure 10. The conserved quantities are given in Table 6 which is almost conserved.Table 6
Birth of solitons witha
=
0
.
5 and b
=
-
3
.
0.
T
I
1
I
2
0.0
17.724343
−12.533142
5.0
17.723816
−12.532956
10.0
17.723352
−12.530116
15.0
17.722782
−12.529239
20.0
17.722217
−12.529013
25.0
17.721734
−12.528983Figure 10
Birth of solitons with parametersa
=
0
.
5 and b
=
-
3
.
0.


## 6. Conclusion

In this work, we have derived two numerical schemes for solving the Hirota-Satsuma CKdV system. The resulting schemes are nonlinear, implicit, and unconditionally stable. The schemes show almost similar results. Single soliton solution and conserved quantities are used to assess the accuracy and the efficiency the derived schemes. We have noticed that the schemes accomplished the aim of preserving the conserved quantities, while maintaining small errors norm during the simulation. The features of an elastic interaction have been shown in the simulation of two and three solitons interaction using the proposed schemes fora
=
1
/
2 and inelasticity occurs for a
≠
1
/
2.To sum up, the derived methods are qualified and can be adopted for solving any CKdV like systems successfully due to their effective performance.



---

*Source: 819367-2014-08-17.xml*