# An Efficient Approach to Numerical Study of the MRLW Equation with B-Spline Collocation Method

Authors: Seydi Battal Gazi Karakoç; Turgut Ak; Halil Zeybek
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A septic B-spline collocation method is implemented to find the numerical solution of the modified regularized long wave (MRLW) equation. Three test problems including the single soliton and interaction of two and three solitons are studied to validate the proposed method by calculating the error normsL
2 and L
∞ and the invariants  I
1, I
2, and I
3. Also, we have studied the Maxwellian initial condition pulse. The numerical results obtained by the method show that the present method is accurate and efficient. Results are compared with some earlier results given in the literature. A linear stability analysis of the method is also investigated.


---

## Body

## 1. Introduction

The generalized regularized long wave equation is given by(1)
U
t
+
U
x
+
δ
U
p
U
x
-
μ
U
x
x
t
=
0
,
where p is a positive integer and δ and μ are positive constants. This equation is one of the most important nonlinear wave equation used a model for small amplitude long waves on the surface of water in a channel [1, 2]. A few authors solved the equation numerically: among others, Zhang [3] used a finite difference method for a Cauchy problem and Kaya [4] applied the Adomian decomposition method and a quasilinearization method based on finite differences was used by Ramos [5]. Roshan [6] implemented the Petrov-Galerkin method using a linear hat function as the trial function and a quintic B-spline function as the test function. A mesh-free technique for the numerical solution of the equation has been presented by Mokhtari and Mohammadi [7]. For   p
=
1,
(2)
U
t
+
U
x
+
δ
U
U
x
-
μ
U
x
x
t
=
0
.
Equation (2) is known as regularized long wave equation, originally introduced to describe the behavior of the undular bore by Peregrine [1] and later widely studied by Benjamin et al. [8]. The RLW equation has been solved numerically by finite element methods [9–22], finite difference methods [23–26], Fourier pseudospectral [27], and mesh-free method [28]. For  p
=
2,
(3)
U
t
+
U
x
+
6
U
2
U
x
-
μ
U
x
x
t
=
0
.
Another particular case of (2) is called modified regularized long wave (MRLW) equation. Like RLW equation, the MRLW equation has been solved numerically by various methods. Among many others, a collocation solution to the equation using quintic B-spline finite element method is developed by Gardner et al. [29]. Khalifa et al. [30, 31] obtained the numerical solutions of the equation using finite difference method and cubic B-spline collocation finite element method. Solutions based on collocation method with quadratic B-spline finite elements and the central finite difference method for time are investigated by Raslan [32]. The equation was solved with a collocation finite element method using quadratic, cubic, quartic, and quintic B-splines to obtain the numerical solutions of the single solitary wave by Raslan and Hassan [33]. Haq et al. [34] have designed a numerical scheme based on quartic B-spline collocation method for the numerical solution of the equation. Ali [35] has formulated a classical radial basis functions (RBFs) collocation method for solving the equation. Karakoc et al. [36] have obtained a type of the quintic B-spline collocation procedure in which nonlinear term in the equation is linearized by using the form introduced by the Rubin and Graves [37] to solve the equation. A Petrov-Galerkin method using cubic B-spline function as trial function and a quadratic B-spline function as the test function is set up to solve the equation by Karakoc and Geyikli [38]. A homotopy analysis method has been employed to obtain approximate numerical solution of the modified regularized long wave (MRLW) equation with some specified initial conditions by Khan et al. [39].In the present paper, a numerical scheme based on the septic B-spline collocation method has been set up for solving the MRLW equation with a variant of both initial and boundary conditions. This paper is set out as follows. In Section2, septic B-spline collocation scheme is presented. Also stability analysis is considered. In Section 3, test problems including single, two, and three solitary waves and Maxwellian initial condition are discussed. Finally in Section 4, a summary is given at the end of the paper.


## 2. Septic B-Spline Finite Element Solution

Consider the MRLW Equation (3) given with the following boundary conditions,
(4)
U
(
a
,
t
)
=
0
,
U
(
b
,
t
)
=
0
,
U
x
(
a
,
t
)
=
U
x
(
b
,
t
)
=
0
,
U
x
x
(
a
,
t
)
=
U
x
x
(
b
,
t
)
=
0
,
U
x
x
x
(
a
,
t
)
=
U
x
x
x
(
b
,
t
)
=
0
,
t
>
0
,
and the initial condition
(5)
U
(
x
,
0
)
=
f
(
x
)
a
≤
x
≤
b
.In order to be able to apply the numerical method, the solution region of the problem is restricted over an intervala
≤
x
≤
b. Space interval [
a
,
b
] is partitioned into uniformly sized finite elements of length h by the nodes x
m such that a
=
x
0
<
x
1
<
⋯
<
x
N
-
1
<
x
N
=
b and h
=
x
m
+
1
-
x
m
,
m
=
-
3
,
-
2
,
…
,
N
+
2
,
N
+
3. The set of septic B-spline functions {
ϕ
-
3
(
x
)
,
ϕ
-
2
(
x
)
,
ϕ
-
1
(
x
)
,
…
,
ϕ
N
+
1
(
x
)
,
ϕ
N
+
2
(
x
)
,
ϕ
N
+
3
(
x
)
} forms a basis over the problem domain [
a
,
b
]. A global approximation U
N
(
x
,
t
) is expressed in terms of septic B-splines as
(6)
U
N
(
x
,
t
)
=
∑
i
=
-
3
N
+
3
‍
ϕ
i
(
x
)
δ
i
(
t
)
,
where δ
i
(
t
)s are time dependent parameters to be determined from the initial, boundary, and collocation conditions.Septic B-splinesϕ
m
(
x
), (m
=
-
3
(
1
)
N
+
3), at the knots x
m are defined over the interval [
a
,
b
] by [40]
(7)
ϕ
m
(
x
)
=
1
h
7
{
(
x
-
x
m
-
4
)
7
[
x
m
-
4
,
x
m
-
3
]
(
x
-
x
m
-
4
)
7
-
8
(
x
-
x
m
-
3
)
7
hhhhhhhlhhh
[
x
m
-
3
,
x
m
-
2
]
(
x
-
x
m
-
4
)
7
-
8
(
x
-
x
m
-
3
)
7
+
28
(
x
-
x
m
-
2
)
7
hhhhhhhlhhh
[
x
m
-
2
,
x
m
-
1
]
(
x
-
x
m
-
4
)
7
-
8
(
x
-
x
m
-
3
)
7
+
28
(
x
-
x
m
-
2
)
7
-
56
(
x
-
x
m
-
1
)
7
hhhhhhhhlhh
[
x
m
-
1
,
x
m
]
(
x
m
+
4
-
x
)
7
-
8
(
x
m
+
3
-
x
)
7
+
28
(
x
m
+
2
-
x
)
7
-
56
(
x
m
+
1
-
x
)
7
hhhhhhhhlhh
[
x
m
,
x
m
+
1
]
(
x
m
+
4
-
x
)
7
-
8
(
x
m
+
3
-
x
)
7
+
28
(
x
m
+
2
-
x
)
7
hhhhhhhlhhh
[
x
m
+
1
,
x
m
+
2
]
(
x
m
+
4
-
x
)
7
-
8
(
x
m
+
3
-
x
)
7
hhhhhhlhhhh
[
x
m
+
2
,
x
m
+
3
]
(
x
m
+
4
-
x
)
7
[
x
m
+
3
,
x
m
+
4
]
0
hhhhhhhhhh
otherwise
.Using expansion (6) and trial function (7), the nodal values (
U
N
)
m and their first, second, and third derivatives (
U
N
)
x
,
(
U
N
)
x
x
,
(
U
N
)
x
x
x can be calculated at the nodal points x
m in terms of nodal parameters δ
m by the following set of equations:
(8)
(
U
N
)
m
=
δ
m
-
3
+
120
δ
m
-
2
+
1191
δ
m
-
1
+
2416
δ
m
(
U
N
x
)
m
=
+
1191
δ
m
+
1
+
120
δ
m
+
2
+
δ
m
+
3
,
(
U
N
x
)
m
=
7
h
(
-
δ
m
-
3
-
56
δ
m
-
2
-
245
δ
m
-
1
(
U
N
x
)
m
hhl
=
+
245
δ
m
+
1
+
56
δ
m
+
2
+
δ
m
+
3
)
,
(
U
N
x
x
)
m
=
42
h
2
(
δ
m
-
3
+
24
δ
m
-
2
+
15
δ
m
-
1
-
80
δ
m
(
U
N
x
)
m
hhhh
=
+
15
δ
m
+
1
+
24
δ
m
+
2
+
δ
m
+
3
)
,
(
U
N
x
x
x
)
m
=
210
h
3
(
-
δ
m
-
3
-
8
δ
m
-
2
+
19
δ
m
-
1
(
U
N
x
x
x
)
m
hhhhl
=
-
19
δ
m
+
1
+
8
δ
m
+
2
+
δ
m
+
3
)
.
The splines ϕ
m
(
x
)  and their six principle derivatives vanish outside the interval [
x
m
-
4
,
x
m
+
4
].To apply the proposed method, Crank-Nicolson approximation for the space derivativesU
x and U
x
x and usual first order forward difference formula for the time derivative of the U in (3) have been used, which lead to
(9)
U
n
+
1
-
U
n
Δ
t
+
U
x
n
+
1
+
U
x
n
2
+
6
(
U
2
U
x
)
n
+
1
+
(
U
2
U
x
)
n
2
-
μ
U
x
x
n
+
1
-
U
x
x
n
Δ
t
=
0
.In order to linearize the nonlinear term (U
2
U
x
)
n
+
1, we can write the term as follows,
(10)
U
2
U
x
=
U
U
U
x
,
and apply the linearization form introduced by Rubin and Graves [37]
(11)
(
U
U
x
)
n
+
1
=
U
n
+
1
U
x
n
+
U
n
U
x
n
+
1
-
U
n
U
x
n
to (9),  and we get
(12)
U
n
+
1
+
Δ
t
2
U
x
n
+
1
+
3
Δ
t
U
n
+
1
×
(
U
n
+
1
U
x
n
+
U
n
U
x
n
+
1
-
U
n
U
x
n
)
-
μ
U
x
x
n
+
1
=
U
n
-
Δ
t
2
U
x
n
-
3
Δ
t
U
n
(
U
U
x
)
n
-
μ
U
x
x
n
.Substituting the approximate solutionU
N and putting the nodal values of U and its derivatives given by (8) into (12) one obtains the following iterative system for m
=
0,1
,
…
,
N
:
(13)
δ
m
-
3
n
+
1
(
1
+
α
1
-
α
2
-
α
3
)
+
δ
m
-
2
n
+
1
(
120
+
120
α
1
-
56
α
2
-
24
α
3
)
+
δ
m
-
1
n
+
1
(
1191
+
1191
α
1
-
245
α
2
-
15
α
3
)
+
δ
m
n
+
1
(
2416
+
2416
α
1
+
80
α
3
)
+
δ
m
+
1
n
+
1
(
1191
+
1191
α
1
+
245
α
2
-
15
α
3
)
+
δ
m
+
2
n
+
1
(
120
+
120
α
1
+
56
α
2
-
24
α
3
)
+
δ
m
+
3
n
+
1
(
1
+
α
1
+
α
2
-
α
3
)
=
δ
m
-
3
n
(
1
-
α
1
+
α
2
-
α
3
)
+
δ
m
-
2
n
(
120
-
120
α
1
+
56
α
2
-
24
α
3
)
+
δ
m
-
1
n
(
1191
-
1191
α
1
+
245
α
2
-
15
α
3
)
+
δ
m
n
(
2416
-
2416
α
1
+
80
α
3
)
+
δ
m
+
1
n
(
1191
-
1191
α
1
-
245
α
2
-
15
α
3
)
+
δ
m
+
2
n
(
120
-
120
α
1
-
56
α
2
-
24
α
3
)
+
δ
m
+
3
n
(
1
-
α
1
-
α
2
-
α
3
)
,
where
(14)
α
1
=
3
Δ
t
,
α
2
=
7
Δ
t
2
h
,
α
3
=
42
μ
h
2
.
The newly obtained iterative system (13) consists of N
+
1 linear equation in N
+
7 unknowns (
δ
-
3
,
δ
-
2
,
δ
-
1
,
δ
0
,
…
,
δ
N
,
δ
N
+
1
,
δ
N
+
2
,
δ
N
+
3
)
T. To obtain a unique solution of this system, six additional constraints are required. Applying the boundary conditions (4) and using the values of (8), these constraints are used and this enables us to eliminate the unknowns δ
-
3
,
δ
-
2
,
δ
-
1
,
δ
0
,
…
,
δ
N
,
δ
N
+
1
,
δ
N
+
2
,
δ
N
+
3 from system (13). So system (13) is reduced to a septa-diagonal system of (
N
+
1
) linear equations in (
N
+
1
) unknowns given by A
δ
n
+
1
=
B
δ
n, where δ
n
=
(
δ
0
,
δ
1
,
δ
2
,
…
,
δ
N
-
2
,
δ
N
-
1
,
δ
N
)
T. The coefficient matrixes are given by(15)
A
=
[
a
11
a
12
a
13
a
14
a
21
a
22
a
23
a
24
a
25
a
31
a
32
a
33
a
34
a
35
a
36
⋱
a
m
,
m
-
3
a
m
,
m
-
2
a
m
,
m
-
1
a
m
,
m
a
m
,
m
+
1
a
m
,
m
+
2
a
m
,
m
+
3
⋱
a
n
-
1
,
n
-
4
a
n
-
1
,
n
-
3
a
n
-
1
,
n
-
2
a
n
-
1
,
n
-
1
a
n
-
1
,
n
a
n
-
1
,
n
+
1
a
n
,
n
-
3
a
n
,
n
-
2
a
n
,
n
-
1
a
n
,
n
a
n
,
n
+
1
a
n
+
1
,
n
-
2
a
n
+
1
,
n
-
1
a
n
+
1
,
n
a
n
+
1
,
n
+
1
]
,
B
=
[
b
11
b
12
b
13
b
14
b
21
b
22
b
23
b
24
b
25
b
31
b
32
b
33
b
34
b
35
b
36
⋱
b
m
,
m
-
3
b
m
,
m
-
2
b
m
,
m
-
1
b
m
,
m
b
m
,
m
+
1
b
m
,
m
+
2
b
m
,
m
+
3
⋱
b
n
-
1
,
n
-
4
b
n
-
1
,
n
-
3
b
n
-
1
,
n
-
2
b
n
-
1
,
n
-
1
b
n
-
1
,
n
b
n
-
1
,
n
+
1
b
n
,
n
-
3
b
n
,
n
-
2
b
n
,
n
-
1
b
n
,
n
b
n
,
n
+
1
b
n
+
1
,
n
-
2
b
n
+
1
,
n
-
1
b
n
+
1
,
n
b
n
+
1
,
n
+
1
]
,where
(16)
a
11
=
-
41472
297
α
3
,
b
11
=
-
41472
297
α
3
,
a
12
=
-
73224
297
α
3
,
b
12
=
-
73224
297
α
3
,
a
13
=
-
20736
297
α
3
,
b
13
=
-
20736
297
α
3
,
a
14
=
-
648
297
α
3
,
b
14
=
-
648
297
α
3
,
a
21
=
215411
297
+
215411
297
α
1
-
11761
297
α
2
+
17893
297
α
3
,
b
21
=
215411
297
-
215411
297
α
1
+
11761
297
α
2
+
17893
297
α
3
,
a
22
=
616822.5
297
+
616822.5
297
α
1
+
43513.5
297
α
2
+
38665.5
297
α
3
,
b
22
=
616822.5
297
-
616822.5
297
α
1
-
43513.5
297
α
2
+
38665.5
297
α
3
,
a
23
=
340284
297
+
340284
297
α
1
+
78528
297
α
2
-
2532
297
α
3
,
b
23
=
340284
297
-
340284
297
α
1
-
78528
297
α
2
-
2532
297
α
3
,
a
24
=
35525.5
297
+
35525.5
297
α
1
+
16682.5
297
α
2
-
7109.5
297
α
3
,
b
24
=
35525.5
297
+
35525.5
297
α
1
+
16682.5
297
α
2
-
7109.5
297
α
3
,
a
25
=
1
+
α
1
+
α
2
-
α
3
,
b
25
=
1
-
α
1
-
α
2
-
α
3
,
a
31
=
34432
297
+
34432
297
α
1
-
15424
297
α
2
-
5920
297
α
3
,
b
31
=
34432
297
-
34432
297
α
1
+
15424
297
α
2
-
5920
297
α
3
,
a
32
=
352833
297
+
352833
297
α
1
-
71871
297
α
2
-
3561
297
α
3
,
b
32
=
352833
297
-
352833
297
α
1
+
71871
297
α
2
-
3561
297
α
3
,
a
33
=
717432
297
+
717432
297
α
1
+
120
297
α
2
+
23880
297
α
3
,
b
33
=
717432
297
-
717432
297
α
1
-
120
297
α
2
+
23880
297
α
3
,
a
34
=
353726
297
+
353726
297
α
1
+
72766
297
α
2
-
4454
297
α
3
,
b
34
=
353726
297
-
353726
297
α
1
-
72766
297
α
2
-
4454
297
α
3
,
a
35
=
120
+
120
α
1
+
56
α
2
-
24
α
3
,
b
35
=
120
-
120
α
1
-
56
α
2
-
24
α
3
,
a
36
=
1
+
α
1
+
α
2
-
α
3
,
b
36
=
1
+
α
1
+
α
2
-
α
3
,
a
m
,
m
-
3
=
1
+
α
1
-
α
2
-
α
3
,
b
m
,
m
-
3
=
1
-
α
1
+
α
2
-
α
3
,
a
m
,
m
-
2
=
120
+
120
α
1
-
56
α
2
-
24
α
3
,
b
m
,
m
-
2
=
120
-
120
α
1
+
56
α
2
-
24
α
3
,
a
m
,
m
-
1
=
1191
+
1191
α
1
-
245
α
2
-
15
α
3
,
b
m
,
m
-
1
=
1191
-
1191
α
1
+
245
α
2
-
15
α
3
,
a
m
,
m
=
2416
+
2416
α
1
+
80
α
3
,
b
m
,
m
=
2416
-
2416
α
1
+
80
α
3
,
a
m
,
m
+
1
=
1191
+
1191
α
1
+
245
α
2
-
15
α
3
,
b
m
,
m
+
1
=
1191
-
1191
α
1
-
245
α
2
-
15
α
3
,
a
m
,
m
+
2
=
120
+
120
α
1
+
56
α
2
-
24
α
3
,
b
m
,
m
+
2
=
120
-
120
α
1
-
56
α
2
-
24
α
3
,
a
m
,
m
+
3
=
1
+
α
1
+
α
2
-
α
3
,
b
m
,
m
+
3
=
1
-
α
1
-
α
2
-
α
3
,
h
h
h
h
h
h
h
h
l
m
=
4
(
1
)
n
-
2
,
a
n
-
1
,
n
-
4
=
1
+
α
1
-
α
2
-
α
3
,
b
n
-
1
,
n
-
4
=
1
-
α
1
+
α
2
-
α
3
,
a
n
-
1
,
n
-
3
=
120
+
120
α
1
-
56
α
2
-
24
α
3
,
b
n
-
1
,
n
-
3
=
120
-
120
α
1
+
56
α
2
-
24
α
3
,
a
n
-
1
,
n
-
2
=
353726
297
+
353726
297
α
1
-
72766
297
α
2
-
4454
297
α
3
,
b
n
-
1
,
n
-
2
=
353726
297
-
353726
297
α
1
+
72766
297
α
2
-
4454
297
α
3
,
a
n
-
1
,
n
-
1
=
717432
297
+
717432
297
α
1
-
120
297
α
2
+
23880
297
α
3
,
b
n
-
1
,
n
-
1
=
717432
297
-
717432
297
α
1
+
120
297
α
2
+
23880
297
α
3
,
a
n
-
1
,
n
=
352833
297
+
352833
297
α
1
+
71871
297
α
2
-
3561
297
α
3
,
b
n
-
1
,
n
=
352833
297
-
352833
297
α
1
-
71871
297
α
2
-
3561
297
α
3
,
a
n
-
1
,
n
+
1
=
34432
297
+
34432
297
α
1
+
15424
297
α
2
-
5920
297
α
3
,
b
n
-
1
,
n
+
1
=
34432
297
-
34432
297
α
1
-
15424
297
α
2
-
5920
297
α
3
,
a
n
,
n
-
3
=
1
+
α
1
-
α
2
-
α
3
,
b
n
,
n
-
3
=
1
-
α
1
+
α
2
-
α
3
,
a
n
,
n
-
2
=
35525.5
297
+
35525.5
297
α
1
-
16682.5
297
α
2
-
7109.5
297
α
3
,
b
n
,
n
-
2
=
35525.5
297
-
35525.5
297
α
1
+
16682.5
297
α
2
-
7109.5
297
α
3
,
a
n
,
n
-
1
=
340284
297
+
340284
297
α
1
-
78528
297
α
2
-
2532
297
α
3
,
b
n
,
n
-
1
=
340284
297
-
340284
297
α
1
+
78528
297
α
2
-
2532
297
α
3
,
a
n
,
n
=
616822.5
297
+
616822.5
297
α
1
-
43513.5
297
α
2
+
38665.5
297
α
3
,
b
n
,
n
=
616822.5
297
-
616822.5
297
α
1
+
43513.5
297
α
2
+
38665.5
297
α
3
,
a
n
,
n
+
1
=
215411
297
+
215411
297
α
1
+
11761
297
α
2
+
17893
297
α
3
,
b
n
,
n
+
1
=
215411
297
-
215411
297
α
1
-
11761
297
α
2
+
17893
297
α
3
,
a
n
+
1
,
n
-
2
=
-
648
297
α
3
,
b
n
+
1
,
n
-
2
=
-
648
297
α
3
,
a
n
+
1
,
n
-
1
=
-
20736
297
α
3
,
b
n
+
1
,
n
-
1
=
-
20736
297
α
3
,
a
n
+
1
,
n
=
-
73224
297
α
3
,
b
n
+
1
,
n
=
-
73224
297
α
3
,
a
n
+
1
,
n
+
1
=
-
41472
297
α
3
,
b
n
+
1
,
n
+
1
=
-
41472
297
α
3
.Before starting the solution process, initial parametersδ
0 must be determined by using the initial condition and the following derivatives at the boundaries:
(17)
U
N
(
x
,
0
)
=
U
(
x
m
,
0
)
m
=
0,1
,
2
,
…
,
N
,
(
U
N
)
x
(
a
,
0
)
=
0
,
(
U
N
)
x
(
b
,
0
)
=
0
,
(
U
N
)
x
x
(
a
,
0
)
=
0
,
(
U
N
)
x
x
(
b
,
0
)
=
0
,
(
U
N
)
x
x
x
(
a
,
0
)
=
0
,
(
U
N
)
x
x
x
(
b
,
0
)
=
0
.
So we have the following matrix form for the initial vector δ
0:
(18)
M
δ
0
=
k
,
where(19)
M
=
[
1536
2712
768
24
82731
81
210568.5
81
104796
81
10063.5
81
1
9600
81
96597
81
195768
81
96474
81
120
1
⋱
1
120
1191
2416
1191
120
1
1
120
96474
81
195768
81
96597
81
9600
81
1
10063.5
81
104796
81
210568.5
81
82731
81
24
768
2712
1536
]
,
δ
0
=
(
δ
0
,
δ
1
,
δ
2
,
…
,
δ
N
-
2
,
δ
N
-
1
,
δ
N
)
T
,and k
=
(
U
(
x
0
,
0
)
,
U
(
x
1
,
0
)
,
…
,
U
(
x
N
-
1
,
0
)
,
U
(
x
N
,
0
)
)
T.


### 2.1. A Linear Stability Analysis

We have investigated stability analysis by applying the von-Neumann approach in which the growth factor of typical Fourier mode is given by(20)
δ
j
n
=
ξ
n
e
i
j
k
h
,
(
i
=
-
1
)
,
where k is a mode number and h is the element size. To apply this method, we have linearized the nonlinear term U
2
U
x by considering U
2 as a constant such as Z
i in (9). If we substitute (20) into the iterative system (13) we obtain the following equation:
(21)
ξ
n
+
1
=
g
ξ
n
,
where g is the growth factor. We have identified the collocation points with the nodes and used (8) to evaluate U
m and its space derivatives in (3). This leads to a set of ordinary differential equations in the following form:
(22)
δ
˙
m
-
3
+
120
δ
˙
m
-
2
+
1191
δ
˙
m
-
1
+
2416
δ
˙
m
+
1191
δ
˙
m
+
1
+
120
δ
˙
m
+
2
+
δ
˙
m
+
3
+
7
h
(
1
+
6
Z
m
)
×
(
-
δ
m
-
3
-
56
δ
m
-
2
-
245
δ
m
-
1
hhhhhh
+
245
δ
m
+
1
+
56
δ
m
+
2
+
δ
m
+
3
)
-
42
μ
h
2
(
δ
˙
m
-
3
+
24
δ
˙
m
-
2
+
15
δ
˙
m
-
1
-
80
δ
˙
m
h
h
h
h
h
h
h
h
h
h
h
+
15
δ
˙
m
+
1
+
24
δ
˙
m
+
2
+
δ
˙
m
+
3
)
=
0
,
where Z
m
=
(
δ
m
-
3
+
120
δ
m
-
2
+
1191
δ
m
-
1
+
2416
δ
m
+
1191
δ
m
+
1
+
120
δ
m
+
2
+
δ
m
+
3
)
2. Here · denotes derivative with respect to time. If the parameters δ
i’s and their time derivatives in (22) are discretized by the Crank-Nicolson formula and usual forward finite difference approximation, respectively,
(23)
δ
i
=
δ
n
+
δ
n
+
1
2
,
δ
˙
i
=
δ
n
+
1
-
δ
n
Δ
t
,
we obtain a recurrence relationship between two time levels n and n
+
1 relating two unknown parameters δ
i
n
+
1, δ
i
n for i
=
m
-
3
,
m
-
2
,
…
,
m
+
2
,
m
+
3,
(24)
γ
1
δ
m
-
3
n
+
1
+
γ
2
δ
m
-
2
n
+
1
+
γ
3
δ
m
-
1
n
+
1
+
γ
4
δ
m
n
+
1
+
γ
5
δ
m
+
1
n
+
1
+
γ
6
δ
m
+
2
n
+
1
+
γ
7
δ
m
+
3
n
+
1
=
γ
7
δ
m
-
3
n
+
γ
6
δ
m
-
2
n
+
γ
5
δ
m
-
1
n
+
γ
4
δ
m
n
+
γ
3
δ
m
+
1
n
+
γ
2
δ
m
+
2
n
+
γ
1
δ
m
+
3
n
,
where
(25)
γ
1
=
1
-
E
-
M
,
γ
2
=
120
-
56
E
-
24
M
,
γ
3
=
1191
-
245
E
-
15
M
,
γ
4
=
2416
+
80
M
,
γ
5
=
1191
+
245
E
-
15
M
,
γ
6
=
120
+
56
E
-
24
M
,
γ
7
=
1
+
E
-
M
,
m
=
0,1
,
…
,
N
,
E
=
7
Δ
t
2
h
(
1
+
6
Z
m
)
,
M
=
42
μ
h
2
.Substituting the Fourier mode (20) into (24) leads to the growth factor of the form
(26)
ξ
=
a
-
i
b
a
+
i
b
,
where
(27)
a
=
1208
+
40
M
+
(
1191
-
15
M
)
cos
⁡
[
h
k
]
+
(
120
-
24
M
)
cos
⁡
[
2
h
k
]
+
(
1
-
M
)
cos
⁡
[
3
h
k
]
,
b
=
245
E
sin
[
h
k
]
+
56
E
sin
[
2
h
k
]
+
E
sin
[
3
h
k
]
,
|
ξ
|
2
=
1; therefore the linearized scheme is unconditionally stable.


## 2.1. A Linear Stability Analysis

We have investigated stability analysis by applying the von-Neumann approach in which the growth factor of typical Fourier mode is given by(20)
δ
j
n
=
ξ
n
e
i
j
k
h
,
(
i
=
-
1
)
,
where k is a mode number and h is the element size. To apply this method, we have linearized the nonlinear term U
2
U
x by considering U
2 as a constant such as Z
i in (9). If we substitute (20) into the iterative system (13) we obtain the following equation:
(21)
ξ
n
+
1
=
g
ξ
n
,
where g is the growth factor. We have identified the collocation points with the nodes and used (8) to evaluate U
m and its space derivatives in (3). This leads to a set of ordinary differential equations in the following form:
(22)
δ
˙
m
-
3
+
120
δ
˙
m
-
2
+
1191
δ
˙
m
-
1
+
2416
δ
˙
m
+
1191
δ
˙
m
+
1
+
120
δ
˙
m
+
2
+
δ
˙
m
+
3
+
7
h
(
1
+
6
Z
m
)
×
(
-
δ
m
-
3
-
56
δ
m
-
2
-
245
δ
m
-
1
hhhhhh
+
245
δ
m
+
1
+
56
δ
m
+
2
+
δ
m
+
3
)
-
42
μ
h
2
(
δ
˙
m
-
3
+
24
δ
˙
m
-
2
+
15
δ
˙
m
-
1
-
80
δ
˙
m
h
h
h
h
h
h
h
h
h
h
h
+
15
δ
˙
m
+
1
+
24
δ
˙
m
+
2
+
δ
˙
m
+
3
)
=
0
,
where Z
m
=
(
δ
m
-
3
+
120
δ
m
-
2
+
1191
δ
m
-
1
+
2416
δ
m
+
1191
δ
m
+
1
+
120
δ
m
+
2
+
δ
m
+
3
)
2. Here · denotes derivative with respect to time. If the parameters δ
i’s and their time derivatives in (22) are discretized by the Crank-Nicolson formula and usual forward finite difference approximation, respectively,
(23)
δ
i
=
δ
n
+
δ
n
+
1
2
,
δ
˙
i
=
δ
n
+
1
-
δ
n
Δ
t
,
we obtain a recurrence relationship between two time levels n and n
+
1 relating two unknown parameters δ
i
n
+
1, δ
i
n for i
=
m
-
3
,
m
-
2
,
…
,
m
+
2
,
m
+
3,
(24)
γ
1
δ
m
-
3
n
+
1
+
γ
2
δ
m
-
2
n
+
1
+
γ
3
δ
m
-
1
n
+
1
+
γ
4
δ
m
n
+
1
+
γ
5
δ
m
+
1
n
+
1
+
γ
6
δ
m
+
2
n
+
1
+
γ
7
δ
m
+
3
n
+
1
=
γ
7
δ
m
-
3
n
+
γ
6
δ
m
-
2
n
+
γ
5
δ
m
-
1
n
+
γ
4
δ
m
n
+
γ
3
δ
m
+
1
n
+
γ
2
δ
m
+
2
n
+
γ
1
δ
m
+
3
n
,
where
(25)
γ
1
=
1
-
E
-
M
,
γ
2
=
120
-
56
E
-
24
M
,
γ
3
=
1191
-
245
E
-
15
M
,
γ
4
=
2416
+
80
M
,
γ
5
=
1191
+
245
E
-
15
M
,
γ
6
=
120
+
56
E
-
24
M
,
γ
7
=
1
+
E
-
M
,
m
=
0,1
,
…
,
N
,
E
=
7
Δ
t
2
h
(
1
+
6
Z
m
)
,
M
=
42
μ
h
2
.Substituting the Fourier mode (20) into (24) leads to the growth factor of the form
(26)
ξ
=
a
-
i
b
a
+
i
b
,
where
(27)
a
=
1208
+
40
M
+
(
1191
-
15
M
)
cos
⁡
[
h
k
]
+
(
120
-
24
M
)
cos
⁡
[
2
h
k
]
+
(
1
-
M
)
cos
⁡
[
3
h
k
]
,
b
=
245
E
sin
[
h
k
]
+
56
E
sin
[
2
h
k
]
+
E
sin
[
3
h
k
]
,
|
ξ
|
2
=
1; therefore the linearized scheme is unconditionally stable.


## 3. Numerical Examples and Results

In this section, we have obtained numerical solution of the MRLW equation for motion of single solitary wave, interaction of two and three solitary waves, and development of the Maxwellian initial condition into solitary waves. Accuracy of the method is measured by using the following error norms:(28)
L
2
=
∥
U
exact
-
U
N
∥
2
≃
h
∑
J
=
1
N
‍
|
U
j
exact
-
(
U
N
)
j
|
2
,
L
∞
=
∥
U
exact
-
U
N
∥
∞
≃
max
⁡
j
|
U
j
exact
-
(
U
N
)
j
|
.
The discrete conservation properties of the MRLW equation corresponding to mass, momentum, and energy are determined by finding the following three invariants [41]:
(29)
I
1
=
∫
a
b
‍
U
d
x
≃
h
∑
J
=
1
N
‍
U
j
n
,
I
2
=
∫
a
b
‍
[
U
2
+
μ
(
U
x
)
2
]
d
x
≃
h
∑
J
=
1
N
‍
[
(
U
j
n
)
2
+
μ
(
U
x
)
j
n
]
,
I
3
=
∫
a
b
‍
(
U
4
-
μ
U
x
2
)
d
x
≃
h
∑
J
=
1
N
‍
[
(
U
j
n
)
4
-
μ
(
U
x
)
j
n
]
.


### 3.1. The Motion of Single Solitary Wave

As a first problem, MRLW equation (3) is considered with the boundary conditions U
→
0 as x
→
±
∞ and the initial condition
(30)
U
(
x
,
0
)
=
c
sec
⁡
h
[
p
(
x
-
x
0
)
]
.
Single solitary wave solution of the MRLW equation has an analytical solution of the form
(31)
U
(
x
,
t
)
=
c
sec
⁡
h
[
p
(
x
-
(
c
+
1
)
t
-
x
0
)
]
,
where p
=
c
/
μ
(
c
+
1
) and x
0 and c are arbitrary constants. This solution corresponds to motion of single solitary wave with amplitude c, initially centered at x
0 and with wave velocity 1
+
μ
c. For this problem the analytical values of the invariants are [29]
(32)
I
1
=
π
c
p
,
I
2
=
2
c
p
+
2
μ
p
c
3
,
I
3
=
4
c
2
3
p
-
2
μ
p
c
3
.For the computational work, two sets of parameters have been chosen and discussed. First of all, we have taken the parametersc
=
1
,
μ
=
1
,
h
=
0.2
,
x
0
=
40
,
k
=
0.025 over the interval [
0,100
] to compare our results with [6, 29, 30, 35, 36]. Thus, the solitary wave has an amplitude 1.0 and the computations are done up to time t
=
10 to obtain the invariants and error norms L
2 and L
∞. Values of the three invariants and error norms are reported in Table 1. It is clearly seen from the table that the error norms are satisfactorily small enough and the computed values of invariants are in good agreement with their analytical values I
1
=
4.4428829
,
I
2
=
3.2998316, and I
3
=
1.4142135. The percentage of the relative error of the conserved quantities I
1, I
2, and I
3 is calculated with respect to the conserved quantities at t
=
0. Percentage of relative changes of I
1, I
2, and I
3 is found to be 0.001
×
1
0
-
3
%
,
3.261
×
1
0
-
3
%, and 7.609
×
1
0
-
3
%, respectively. Thus, the quantities I
1, I
2, and I
3 remain constant during the computer run. For this case in Table 2, we compare the values of the invariants and error norms obtained by using the present method and some earlier methods [6, 29, 30, 35, 36]. From the table, we observed that the error norms obtained by the present method are less than those of other methods [6, 29, 30, 35, 36]. Figure 1 shows the motion of solitary wave with c
=
1
,
h
=
0.2
,
k
=
0.025 at different time levels. It is observed that the solitary wave moves to the right with constant velocity and amplitude. At t
=
0, the amplitude is 1.0 which is located at x
=
40, while it is 0.999950 at t
=
10 located at x
=
60. The absolute difference in amplitudes at times t
=
0 and t
=
10 is found to be 5
×
1
0
-
5, so there is a little change between the amplitudes.Table 1
Invariants and error norms for single solitary wave withc
=
1, h
=
0.2, k
=
0.025, 0
≤
x
≤
100.
t
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
4.4428661
3.2998227
1.4142046
0.00000000
0.00000000
1
4.4428661
3.2998085
1.4142188
0.28537793
0.16594258
2
4.4428661
3.2997808
1.4142465
0.56248008
0.31854916
3
4.4428661
3.2997573
1.4142700
0.82836630
0.45535369
4
4.4428661
3.2997415
1.4142858
1.08566992
0.58528925
5
4.44286601
3.2997313
1.4142960
1.33772774
0.71261445
6
4.4428661
3.2997248
1.4143025
1.58675627
0.83879372
7
4.4428661
3.2997207
1.4143067
1.83403948
0.96441682
8
4.4428661
3.2997180
1.4143093
2.08032250
1.08975930
9
4.4428661
3.2997162
1.4143111
2.32602024
1.21494581
10
4.4428661
3.2997151
1.4143122
2.57148152
1.34021078Table 2
Errors and invariants for single solitary wave withc
=
1, h
=
0.2, k
=
0.025, 0
≤
x
≤
100, at t
=
10.
Method
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
Analytical
4.4428829
3.2998316
1.4142135
0
0
Present
4.4428661
3.2997151
1.4143122
2.57148
1.334021
[6]
4.44288
3.29981
1.41416
3.00533
1.68749
Cubic B-splines coll-CN [29]
4.442
3.299
1.413
16.39
9.24
Cubic B-splines coll-PA-CN+ [29]
4.440
3.296
1.411
20.3
11.2
Cubic B-splines coll [30]
4.44288
3.29983
1.41420
9.30196
5.43718
MQ [35]
4.4428829
3.29978
1.414163
3.914
2.019
IMQ [35]
4.4428611
3.29978
1.414163
3.914
2.019
IQ [35]
4.4428794
3.29978
1.414163
3.914
2.019
GA [35]
4.4428829
3.29978
1.414163
3.914
2.019
TPS [35]
4.4428821
3.29972
1.414104
4.428
2.306
Quintic B-splines coll [36]
4.4428661
3.2997108
1.4143165
2.58891
1.35164Figure 1
Single solitary wave withc
=
1
,
h
=
0.2
,
Δ
t
=
0.025,0
≤
x
≤
100 at times t
=
0,2
,
4,6
,
8, and 10.Secondly, we have taken the parametersμ
=
1
,
c
=
0.3
,
h
=
0.1
,
k
=
0.01 and x
0
=
40 with range [
0,100
] to enable comparison with [6, 30, 32, 34–36]. So the solitary wave has amplitude 0.547723. Simulations are run up to time t
=
20. Error norms L
2 and L
∞ and conserved quantities are tabulated in Table 3, together with the results obtained in [6, 30, 32, 34–36]. These results show high degree of accuracy and efficiency of the method. The invariants I
1, I
2, and I
3 have changed by less than 0.001
×
1
0
-
3
%
,
0.189
×
1
0
-
3
%, and 1.656
×
1
0
-
3
% percent, respectively. Moreover, the variation of the invariants I
2 and I
3 from the initial variants is less than 2.5
×
1
0
-
6 and 2.6
×
1
0
-
6, respectively, whereas the changes in the invariant I
1 approach to zero throughout the interval. The perspective views of traveling solitary wave at different time levels are shown in Figure 2. The distributions of the errors at time t
=
10 and t
=
20 are shown graphically for solitary waves amplitudes 1 and 0.3 in Figure 3. It is seen that the maximum errors are close to the tip of the solitary waves and between -
6
×
1
0
-
3 and 6
×
1
0
-
3, -
2
×
1
0
-
4, and 2
×
1
0
-
4, respectively. The CPU times for an Intel(R) Core i5, 2.53 GHz, are also given in Tables 4 and 5.Table 3
Invariants and error norms for single solitary wave withc
=
0.3, h
=
0.1, k
=
0.01, 0
≤
x
≤
100.
t
I
1
I
2
I
3
L
2
×
1
0
4
L
∞
×
1
0
4
0
3.5819531
1.3450721
0.1537217
0.0000000
0.0000000
2
3.5819531
1.3450719
0.1537219
0.0373696
0.0211791
4
3.5819531
1.3450715
0.1537223
0.0711480
0.0387624
6
3.5819531
1.3450711
0.1537227
0.1001141
0.0515117
8
3.5819531
1.3450708
0.1537231
0.1249329
0.0614203
10
3.5819531
1.3450705
0.1537234
0.1466243
0.0700260
12
3.5819531
1.3450702
0.1537236
0.1659668
0.0775889
14
3.5819531
1.3450700
0.1537238
0.1833628
0.0844911
16
3.5819531
1.3450698
0.1537240
0.2015361
0.0909663
18
3.5819531
1.3450697
0.1537241
0.2560750
0.0993420
20
3.5819531
1.3450696
0.1537243
0.3585031
0.1702101
20 [6]
3.58197
1.34508
0.153723
0.645295
0.301923
20 [30]
3.58197
1.34508
0.153723
6.06885
2.96650
20 [34]
3.581967
1.345076
0.153723
0.508927
0.222284
20 [35] MQ
3.5819665
1.3450764
0.153723
0.51498
0.22551
20 [35] IMQ
3.5819664
1.3450764
0.153723
0.51498
0.22551
20 [35] IQ
3.5819654
1.3450764
0.153723
0.51498
0.22551
20 [35] GA
3.5819665
1.3450764
0.153723
0.51498
0.22551
20 [35] TPS
3.5819663
1.3450759
0.153723
0.51498
0.26605
20 [36]
3.5820204
1.3450974
0.1537250
0.8112594
0.3569076Table 4
Invariants and error norms for single solitary wave withc
=
0.6, h
=
0.1, k
=
0.1, -
40
≤
x
≤
60.
t
I
1
I
2
I
3
L
2
×
1
0
4
L
∞
×
1
0
4
CPU Time
4
8.070902
4.100549
14.361115
2.716636
1.323585
0.437 s
8
8.070925
4.100534
14.361111
5.240855
2.458097
0.908 s
12
8.070943
4.100519
14.361103
7.543043
3.372256
1.435 s
16
8.070949
4.100503
14.361089
9.662561
4.154427
1.774 s
20
8.070921
4.100489
14.361074
11.660256
4.862191
2.319 sTable 5
Invariants and error norms for single solitary wave withc
=
0.18, h
=
0.1, k
=
0.1, -
80
≤
x
≤
120.
t
I
1
I
2
I
3
L
2
×
1
0
5
L
∞
×
1
0
5
CPU Time
4
7.809873
2.129887
7.031111
3.077451
1.138290
0.801 s
8
7.809875
2.129887
7.031112
6.124574
2.313797
1.571 s
12
7.809877
2.129887
7.031112
9.106266
3.455915
2.470 s
16
7.809879
2.129887
7.031113
11.967534
4.530312
3.124 s
20
7.809880
2.129887
7.031115
14.731556
5.538799
3.789 sFigure 2
Single solitary wave withc
=
0.3
,
h
=
0.1
,
Δ
t
=
0.01
,
0
≤
x
≤
100 at times t
=
0,5
,
10,15, and 20.Error with (a)c
=
1
,
h
=
0.2
,
Δ
t
=
0.025
,
t
=
10
,
0
≤
x
≤
100 and (b) c
=
0.3
,
h
=
0.1
,
Δ
t
=
0.01
,
t
=
20,0
≤
x
≤
100.
(a)
(b)


### 3.2. Interaction of Two Solitary Waves

As a second problem, interaction of two well separated solitary waves having different amplitudes and traveling in the same direction is considered by using the initial condition(33)
U
(
x
,
0
)
=
∑
j
=
1
2
‍
A
j
sec
⁡
h
(
p
j
[
x
-
x
j
]
)
,
where A
j
=
c
j, p
j
=
c
j
/
μ
(
c
j
+
1
), j
=
1,2, and c
j and x
j are arbitrary constants. The analytical values of the invariants can be found as [29]
(34)
I
1
=
∑
j
=
1
2
‍
π
c
j
p
j
=
11.467698
,
I
2
=
∑
j
=
1
2
‍
(
2
c
j
p
j
+
2
μ
p
j
c
j
3
)
=
14.629243
,
I
3
=
∑
j
=
1
2
‍
(
4
c
j
2
3
p
j
-
2
μ
p
j
c
j
3
)
=
22.880466
.We have studied the interaction of two positive solitary waves having the parametersμ
=
1
,
h
=
0.2
,
k
=
0.025
,
c
1
=
4
,
c
2
=
1
,
x
1
=
25
,
x
2
=
55 through the interval 0
≤
x
≤
250 to coincide with those used by [6, 30, 34–36]. The simulations are maintained up to t
=
20. Constant values I
1
,
I
2, and I
3 at various time steps together with equivalent results for the previous methods are shown in Table 6. It is seen that the numerical values of the invariants remain almost constant during the computer run. The interaction of two solitary waves is shown in Figure 4. It can be seen from the figure that, at t
=
0, the wave with larger amplitude is to the left of the second wave with smaller amplitude. Since the taller wave moves faster than the shorter one, it catches up and collides with the shorter one at t
=
8 and then moves away from the shorter one as time increases. At t
=
20, the amplitude of larger waves is 2.001102 at the point x
=
127.4, whereas the amplitude of the smaller one is 0.996403 at the point x
=
92. It is found that the absolute difference in amplitude is 3.59
×
1
0
-
3 for the smaller wave and 1.10
×
1
0
-
3 for the larger wave for this case.Table 6
Comparison of the invariants for the interaction of two solitary waves with results from [30] with h
=
0.2, k
=
0.025 in the region 0
≤
x
≤
250.
Present method
[30]
t
I
1
I
2
I
3
CPU Time
I
1
I
2
I
3
0
11.4676
14.6292
22.8803
—
11.4677
14.6291
22.8806
2
11.4678
14.6282
22.8812
0.209 s
11.4677
14.6292
22.8807
4
11.4679
14.6282
22.8812
0.413 s
11.4677
14.6292
22.8807
6
11.4681
14.6181
22.8914
1.101 s
11.4677
14.6295
22.8806
8
11.4675
14.1393
23.3702
1.219 s
11.4677
14.6451
22.8454
10
11.4674
14.0502
23.4592
1.421 s
11.4677
14.5963
22.8913
12
11.4685
14.6817
22.8278
2.027 s
11.4677
14.6287
22.8814
14
11.4687
14.6649
22.8446
2.228 s
11.4677
14.6295
22.8807
16
11.4688
14.6459
22.8635
2.432 s
11.4677
14.6294
22.8808
18
11.4690
14.6370
22.8725
3.004 s
11.4677
14.6293
22.8809
20
11.4691
14.6331
22.8763
3.239 s
11.4677
14.6292
22.8809
20 [6]
11.4677
14.6299
22.8806
20 [34]
11.4677
14.5830
22.6965
20 [35] MQ
11.4676
14.5830
22.6965
20 [35] IMQ
11.4676
14.5830
22.6965
20 [35] IQ
11.4676
14.5830
22.6965
20 [35] GA
11.4676
14.5830
22.6965
20 [35] TPS
11.4677
14.5824
22.6942
20 [36]
11.4691
14.6331
22.8764Figure 4
Interaction of two solitary waves witht
=
0
,
4
,
8
,
10
,
14, and  20.


### 3.3. Interaction of Three Solitary Waves

As a third problem, interaction of three solitary waves having different amplitudes and travelling in the same direction is studied. We consider (3) with initial conditions given by
(35)
U
(
x
,
0
)
=
∑
j
=
1
3
‍
A
j
sec
⁡
h
(
p
j
[
x
-
x
j
]
)
,
where A
j
=
c
j, p
j
=
c
j
/
μ
(
c
j
+
1
), j
=
1,2
,
3, c
j and x
j are arbitrary constants. For this problem the analytical values of the invariants are found from (32) as
(36)
I
1
=
∑
j
=
1
3
‍
π
c
j
p
j
=
14.9801
,
I
2
=
∑
j
=
1
3
‍
(
2
c
j
p
j
+
2
μ
p
j
c
j
3
)
=
15.8218
,
I
3
=
∑
j
=
1
3
‍
(
4
c
j
2
3
p
j
-
2
μ
p
j
c
j
3
)
=
22.9923
.In order to be able to compare with the previous works, computations are performed for the parametersμ
=
1
,
h
=
0.2
,
k
=
0.025
,
c
1
=
4
,
c
2
=
1
,
c
3
=
0.25
,
x
1
=
15
,
x
2
=
45
,
x
3
=
60 over the interval 0
≤
x
≤
250. The simulation is run from t
=
0 to t
=
45. Table 7 shows a comparison of the values of the invariants obtained by the present method with those obtained in [30, 34–36]. It is evident from the table that the obtained values of the invariants remain almost constant during the computer run which are all in good agreement with their analytical values given by (36). Figure 5 depicts the interaction of three solitary waves at different times. It is observed from Figure 5 that interaction started about time t
=
10, overlapping processes occurred between time t
=
15 and t
=
40, and waves started to resume their original shapes after the time t
=
40.Table 7
Comparison of invariants for the interaction of three solitary waves with results from [30] with h
=
0.2, k
=
0.025 in the region 0
≤
x
≤
250.
Present method
[30]
t
I
1
I
2
I
3
CPU Time
I
1
I
2
I
3
0
14.9800
15.8374
23.0081
—
13.6891
15.4549
22.8816
5
14.9804
15.7991
23.0465
0.536 s
13.6891
15.3109
22.6939
10
14.9783
14.3275
24.5180
1.442 s
13.6891
15.6514
22.8388
15
14.9809
15.5048
23.3408
2.349 s
13.6891
15.6548
22.9347
20
14.9814
15.7504
23.0952
3.251 s
13.6891
15.6557
22.9330
25
14.9819
15.8496
22.9960
4.162 s
13.6892
15.6559
22.9336
30
14.9823
15.8626
22.9829
5.063 s
13.6894
15.6559
22.9348
35
14.9828
15.8622
22.9833
5.573 s
13.6913
15.6564
22.9343
40
14.9832
15.8605
22.9851
6.481 s
13.7015
15.6566
22.9335
45
14.9461
15.8613
22.9843
7.388 s
13.7043
15.6563
22.9303
45 [34]
13.7043
15.6563
22.9303
45 [35] MQ
14.9681
15.7343
22.5966
45 [35] IMQ
14.9680
15.7343
22.5966
45 [35] IQ
14.9681
15.7343
22.5966
45 [35] GA
14.9681
15.7343
22.5966
45 [35] TPS
14.9682
15.7337
22.5944
45 [36]
14.7145
15.4927
23.3529Figure 5
Interaction of three solitary waves witht
=
0
,
5
,
8
,
15
,
20
,
and
40.


### 3.4. The Maxwellian Initial Condition

As our last problem, the development of the Maxwellian initial condition,(37)
U
(
x
,
0
)
=
exp
⁡
(
-
(
x
-
40
)
2
)
,
into a train of solitary waves is considered. For the Maxwellian initial condition, behavior of the solution depends on the values of μ. So we take μ
=
0.1
,
μ
=
0.04
,
μ
=
0.015, and μ
=
0.01. The numerical values of the invariants quantities during the simulations are given in Table 8. For μ
=
0.1, only a single soliton occurred as depicted in Figure 6(a). When μ
=
0.015 and μ
=
0.01, two and three stable solitons occurred, respectively, as depicted in Figures 6(b) and 6(c). For μ
=
0.04, Maxwellian initial condition has decayed into four solitary waves as depicted in Figure 6(d). As is seen, when μ is reduced, more solitary waves occurred. All figures are drawn up at time t
=
14.5. The peaks of the well developed wave lie on a straight line so that their velocities are linearly dependent on their amplitudes and also we observe a small oscillating tail appearing behind the last wave in all Maxwellian figures.Table 8
Invariants of MRLW equation using the Maxwellian initial condition.
t
μ
I
1
I
2
I
3
μ
I
1
I
2
I
3
0
1.7724809
1.3786628
0.7609109
1.7724809
1.2721327
0.8674410
3
1.7725930
1.4731551
0.6664186
1.7753954
1.4487372
0.6908365
6
0.1
1.7725929
1.4727329
0.6668408
0.015
1.7770146
1.4518203
0.6877534
9
1.7725935
1.4725552
0.6670185
1.7786793
1.4534252
0.6861485
12
1.7725942
1.4724619
0.6671118
1.7804784
1.4552349
0.6843388
15
1.77259499
1.4724042
0.6671695
1.7824228
1.4573443
0.6822294
0
1.7724809
1.2658662
0.8737075
1.7724809
1.3034651
0.8361087
3
1.7796208
1.4490883
0.6904854
1.7730165
1.4543782
0.6851955
6
0.01
1.7863091
1.4571676
0.6824061
0.04
1.7730562
1.4540759
0.6854978
9
1.7952456
1.4658449
0.6737288
1.7730962
1.4538204
0.6857534
12
1.8089221
1.4802019
0.6593718
1.7731365
1.4537359
0.6858378
15
1.8364987
1.5120445
0.6275292
1.7731770
1.4537234
0.6858503Maxwellian initial condition att
=
14.5 with (a) μ
=
0.1, (b) μ
=
0.015, (c) μ
=
0.01, and (d) μ
=
0.04.
(a)
(b)
(c)
(d)


## 3.1. The Motion of Single Solitary Wave

As a first problem, MRLW equation (3) is considered with the boundary conditions U
→
0 as x
→
±
∞ and the initial condition
(30)
U
(
x
,
0
)
=
c
sec
⁡
h
[
p
(
x
-
x
0
)
]
.
Single solitary wave solution of the MRLW equation has an analytical solution of the form
(31)
U
(
x
,
t
)
=
c
sec
⁡
h
[
p
(
x
-
(
c
+
1
)
t
-
x
0
)
]
,
where p
=
c
/
μ
(
c
+
1
) and x
0 and c are arbitrary constants. This solution corresponds to motion of single solitary wave with amplitude c, initially centered at x
0 and with wave velocity 1
+
μ
c. For this problem the analytical values of the invariants are [29]
(32)
I
1
=
π
c
p
,
I
2
=
2
c
p
+
2
μ
p
c
3
,
I
3
=
4
c
2
3
p
-
2
μ
p
c
3
.For the computational work, two sets of parameters have been chosen and discussed. First of all, we have taken the parametersc
=
1
,
μ
=
1
,
h
=
0.2
,
x
0
=
40
,
k
=
0.025 over the interval [
0,100
] to compare our results with [6, 29, 30, 35, 36]. Thus, the solitary wave has an amplitude 1.0 and the computations are done up to time t
=
10 to obtain the invariants and error norms L
2 and L
∞. Values of the three invariants and error norms are reported in Table 1. It is clearly seen from the table that the error norms are satisfactorily small enough and the computed values of invariants are in good agreement with their analytical values I
1
=
4.4428829
,
I
2
=
3.2998316, and I
3
=
1.4142135. The percentage of the relative error of the conserved quantities I
1, I
2, and I
3 is calculated with respect to the conserved quantities at t
=
0. Percentage of relative changes of I
1, I
2, and I
3 is found to be 0.001
×
1
0
-
3
%
,
3.261
×
1
0
-
3
%, and 7.609
×
1
0
-
3
%, respectively. Thus, the quantities I
1, I
2, and I
3 remain constant during the computer run. For this case in Table 2, we compare the values of the invariants and error norms obtained by using the present method and some earlier methods [6, 29, 30, 35, 36]. From the table, we observed that the error norms obtained by the present method are less than those of other methods [6, 29, 30, 35, 36]. Figure 1 shows the motion of solitary wave with c
=
1
,
h
=
0.2
,
k
=
0.025 at different time levels. It is observed that the solitary wave moves to the right with constant velocity and amplitude. At t
=
0, the amplitude is 1.0 which is located at x
=
40, while it is 0.999950 at t
=
10 located at x
=
60. The absolute difference in amplitudes at times t
=
0 and t
=
10 is found to be 5
×
1
0
-
5, so there is a little change between the amplitudes.Table 1
Invariants and error norms for single solitary wave withc
=
1, h
=
0.2, k
=
0.025, 0
≤
x
≤
100.
t
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
0
4.4428661
3.2998227
1.4142046
0.00000000
0.00000000
1
4.4428661
3.2998085
1.4142188
0.28537793
0.16594258
2
4.4428661
3.2997808
1.4142465
0.56248008
0.31854916
3
4.4428661
3.2997573
1.4142700
0.82836630
0.45535369
4
4.4428661
3.2997415
1.4142858
1.08566992
0.58528925
5
4.44286601
3.2997313
1.4142960
1.33772774
0.71261445
6
4.4428661
3.2997248
1.4143025
1.58675627
0.83879372
7
4.4428661
3.2997207
1.4143067
1.83403948
0.96441682
8
4.4428661
3.2997180
1.4143093
2.08032250
1.08975930
9
4.4428661
3.2997162
1.4143111
2.32602024
1.21494581
10
4.4428661
3.2997151
1.4143122
2.57148152
1.34021078Table 2
Errors and invariants for single solitary wave withc
=
1, h
=
0.2, k
=
0.025, 0
≤
x
≤
100, at t
=
10.
Method
I
1
I
2
I
3
L
2
×
1
0
3
L
∞
×
1
0
3
Analytical
4.4428829
3.2998316
1.4142135
0
0
Present
4.4428661
3.2997151
1.4143122
2.57148
1.334021
[6]
4.44288
3.29981
1.41416
3.00533
1.68749
Cubic B-splines coll-CN [29]
4.442
3.299
1.413
16.39
9.24
Cubic B-splines coll-PA-CN+ [29]
4.440
3.296
1.411
20.3
11.2
Cubic B-splines coll [30]
4.44288
3.29983
1.41420
9.30196
5.43718
MQ [35]
4.4428829
3.29978
1.414163
3.914
2.019
IMQ [35]
4.4428611
3.29978
1.414163
3.914
2.019
IQ [35]
4.4428794
3.29978
1.414163
3.914
2.019
GA [35]
4.4428829
3.29978
1.414163
3.914
2.019
TPS [35]
4.4428821
3.29972
1.414104
4.428
2.306
Quintic B-splines coll [36]
4.4428661
3.2997108
1.4143165
2.58891
1.35164Figure 1
Single solitary wave withc
=
1
,
h
=
0.2
,
Δ
t
=
0.025,0
≤
x
≤
100 at times t
=
0,2
,
4,6
,
8, and 10.Secondly, we have taken the parametersμ
=
1
,
c
=
0.3
,
h
=
0.1
,
k
=
0.01 and x
0
=
40 with range [
0,100
] to enable comparison with [6, 30, 32, 34–36]. So the solitary wave has amplitude 0.547723. Simulations are run up to time t
=
20. Error norms L
2 and L
∞ and conserved quantities are tabulated in Table 3, together with the results obtained in [6, 30, 32, 34–36]. These results show high degree of accuracy and efficiency of the method. The invariants I
1, I
2, and I
3 have changed by less than 0.001
×
1
0
-
3
%
,
0.189
×
1
0
-
3
%, and 1.656
×
1
0
-
3
% percent, respectively. Moreover, the variation of the invariants I
2 and I
3 from the initial variants is less than 2.5
×
1
0
-
6 and 2.6
×
1
0
-
6, respectively, whereas the changes in the invariant I
1 approach to zero throughout the interval. The perspective views of traveling solitary wave at different time levels are shown in Figure 2. The distributions of the errors at time t
=
10 and t
=
20 are shown graphically for solitary waves amplitudes 1 and 0.3 in Figure 3. It is seen that the maximum errors are close to the tip of the solitary waves and between -
6
×
1
0
-
3 and 6
×
1
0
-
3, -
2
×
1
0
-
4, and 2
×
1
0
-
4, respectively. The CPU times for an Intel(R) Core i5, 2.53 GHz, are also given in Tables 4 and 5.Table 3
Invariants and error norms for single solitary wave withc
=
0.3, h
=
0.1, k
=
0.01, 0
≤
x
≤
100.
t
I
1
I
2
I
3
L
2
×
1
0
4
L
∞
×
1
0
4
0
3.5819531
1.3450721
0.1537217
0.0000000
0.0000000
2
3.5819531
1.3450719
0.1537219
0.0373696
0.0211791
4
3.5819531
1.3450715
0.1537223
0.0711480
0.0387624
6
3.5819531
1.3450711
0.1537227
0.1001141
0.0515117
8
3.5819531
1.3450708
0.1537231
0.1249329
0.0614203
10
3.5819531
1.3450705
0.1537234
0.1466243
0.0700260
12
3.5819531
1.3450702
0.1537236
0.1659668
0.0775889
14
3.5819531
1.3450700
0.1537238
0.1833628
0.0844911
16
3.5819531
1.3450698
0.1537240
0.2015361
0.0909663
18
3.5819531
1.3450697
0.1537241
0.2560750
0.0993420
20
3.5819531
1.3450696
0.1537243
0.3585031
0.1702101
20 [6]
3.58197
1.34508
0.153723
0.645295
0.301923
20 [30]
3.58197
1.34508
0.153723
6.06885
2.96650
20 [34]
3.581967
1.345076
0.153723
0.508927
0.222284
20 [35] MQ
3.5819665
1.3450764
0.153723
0.51498
0.22551
20 [35] IMQ
3.5819664
1.3450764
0.153723
0.51498
0.22551
20 [35] IQ
3.5819654
1.3450764
0.153723
0.51498
0.22551
20 [35] GA
3.5819665
1.3450764
0.153723
0.51498
0.22551
20 [35] TPS
3.5819663
1.3450759
0.153723
0.51498
0.26605
20 [36]
3.5820204
1.3450974
0.1537250
0.8112594
0.3569076Table 4
Invariants and error norms for single solitary wave withc
=
0.6, h
=
0.1, k
=
0.1, -
40
≤
x
≤
60.
t
I
1
I
2
I
3
L
2
×
1
0
4
L
∞
×
1
0
4
CPU Time
4
8.070902
4.100549
14.361115
2.716636
1.323585
0.437 s
8
8.070925
4.100534
14.361111
5.240855
2.458097
0.908 s
12
8.070943
4.100519
14.361103
7.543043
3.372256
1.435 s
16
8.070949
4.100503
14.361089
9.662561
4.154427
1.774 s
20
8.070921
4.100489
14.361074
11.660256
4.862191
2.319 sTable 5
Invariants and error norms for single solitary wave withc
=
0.18, h
=
0.1, k
=
0.1, -
80
≤
x
≤
120.
t
I
1
I
2
I
3
L
2
×
1
0
5
L
∞
×
1
0
5
CPU Time
4
7.809873
2.129887
7.031111
3.077451
1.138290
0.801 s
8
7.809875
2.129887
7.031112
6.124574
2.313797
1.571 s
12
7.809877
2.129887
7.031112
9.106266
3.455915
2.470 s
16
7.809879
2.129887
7.031113
11.967534
4.530312
3.124 s
20
7.809880
2.129887
7.031115
14.731556
5.538799
3.789 sFigure 2
Single solitary wave withc
=
0.3
,
h
=
0.1
,
Δ
t
=
0.01
,
0
≤
x
≤
100 at times t
=
0,5
,
10,15, and 20.Error with (a)c
=
1
,
h
=
0.2
,
Δ
t
=
0.025
,
t
=
10
,
0
≤
x
≤
100 and (b) c
=
0.3
,
h
=
0.1
,
Δ
t
=
0.01
,
t
=
20,0
≤
x
≤
100.
(a)
(b)


## 3.2. Interaction of Two Solitary Waves

As a second problem, interaction of two well separated solitary waves having different amplitudes and traveling in the same direction is considered by using the initial condition(33)
U
(
x
,
0
)
=
∑
j
=
1
2
‍
A
j
sec
⁡
h
(
p
j
[
x
-
x
j
]
)
,
where A
j
=
c
j, p
j
=
c
j
/
μ
(
c
j
+
1
), j
=
1,2, and c
j and x
j are arbitrary constants. The analytical values of the invariants can be found as [29]
(34)
I
1
=
∑
j
=
1
2
‍
π
c
j
p
j
=
11.467698
,
I
2
=
∑
j
=
1
2
‍
(
2
c
j
p
j
+
2
μ
p
j
c
j
3
)
=
14.629243
,
I
3
=
∑
j
=
1
2
‍
(
4
c
j
2
3
p
j
-
2
μ
p
j
c
j
3
)
=
22.880466
.We have studied the interaction of two positive solitary waves having the parametersμ
=
1
,
h
=
0.2
,
k
=
0.025
,
c
1
=
4
,
c
2
=
1
,
x
1
=
25
,
x
2
=
55 through the interval 0
≤
x
≤
250 to coincide with those used by [6, 30, 34–36]. The simulations are maintained up to t
=
20. Constant values I
1
,
I
2, and I
3 at various time steps together with equivalent results for the previous methods are shown in Table 6. It is seen that the numerical values of the invariants remain almost constant during the computer run. The interaction of two solitary waves is shown in Figure 4. It can be seen from the figure that, at t
=
0, the wave with larger amplitude is to the left of the second wave with smaller amplitude. Since the taller wave moves faster than the shorter one, it catches up and collides with the shorter one at t
=
8 and then moves away from the shorter one as time increases. At t
=
20, the amplitude of larger waves is 2.001102 at the point x
=
127.4, whereas the amplitude of the smaller one is 0.996403 at the point x
=
92. It is found that the absolute difference in amplitude is 3.59
×
1
0
-
3 for the smaller wave and 1.10
×
1
0
-
3 for the larger wave for this case.Table 6
Comparison of the invariants for the interaction of two solitary waves with results from [30] with h
=
0.2, k
=
0.025 in the region 0
≤
x
≤
250.
Present method
[30]
t
I
1
I
2
I
3
CPU Time
I
1
I
2
I
3
0
11.4676
14.6292
22.8803
—
11.4677
14.6291
22.8806
2
11.4678
14.6282
22.8812
0.209 s
11.4677
14.6292
22.8807
4
11.4679
14.6282
22.8812
0.413 s
11.4677
14.6292
22.8807
6
11.4681
14.6181
22.8914
1.101 s
11.4677
14.6295
22.8806
8
11.4675
14.1393
23.3702
1.219 s
11.4677
14.6451
22.8454
10
11.4674
14.0502
23.4592
1.421 s
11.4677
14.5963
22.8913
12
11.4685
14.6817
22.8278
2.027 s
11.4677
14.6287
22.8814
14
11.4687
14.6649
22.8446
2.228 s
11.4677
14.6295
22.8807
16
11.4688
14.6459
22.8635
2.432 s
11.4677
14.6294
22.8808
18
11.4690
14.6370
22.8725
3.004 s
11.4677
14.6293
22.8809
20
11.4691
14.6331
22.8763
3.239 s
11.4677
14.6292
22.8809
20 [6]
11.4677
14.6299
22.8806
20 [34]
11.4677
14.5830
22.6965
20 [35] MQ
11.4676
14.5830
22.6965
20 [35] IMQ
11.4676
14.5830
22.6965
20 [35] IQ
11.4676
14.5830
22.6965
20 [35] GA
11.4676
14.5830
22.6965
20 [35] TPS
11.4677
14.5824
22.6942
20 [36]
11.4691
14.6331
22.8764Figure 4
Interaction of two solitary waves witht
=
0
,
4
,
8
,
10
,
14, and  20.


## 3.3. Interaction of Three Solitary Waves

As a third problem, interaction of three solitary waves having different amplitudes and travelling in the same direction is studied. We consider (3) with initial conditions given by
(35)
U
(
x
,
0
)
=
∑
j
=
1
3
‍
A
j
sec
⁡
h
(
p
j
[
x
-
x
j
]
)
,
where A
j
=
c
j, p
j
=
c
j
/
μ
(
c
j
+
1
), j
=
1,2
,
3, c
j and x
j are arbitrary constants. For this problem the analytical values of the invariants are found from (32) as
(36)
I
1
=
∑
j
=
1
3
‍
π
c
j
p
j
=
14.9801
,
I
2
=
∑
j
=
1
3
‍
(
2
c
j
p
j
+
2
μ
p
j
c
j
3
)
=
15.8218
,
I
3
=
∑
j
=
1
3
‍
(
4
c
j
2
3
p
j
-
2
μ
p
j
c
j
3
)
=
22.9923
.In order to be able to compare with the previous works, computations are performed for the parametersμ
=
1
,
h
=
0.2
,
k
=
0.025
,
c
1
=
4
,
c
2
=
1
,
c
3
=
0.25
,
x
1
=
15
,
x
2
=
45
,
x
3
=
60 over the interval 0
≤
x
≤
250. The simulation is run from t
=
0 to t
=
45. Table 7 shows a comparison of the values of the invariants obtained by the present method with those obtained in [30, 34–36]. It is evident from the table that the obtained values of the invariants remain almost constant during the computer run which are all in good agreement with their analytical values given by (36). Figure 5 depicts the interaction of three solitary waves at different times. It is observed from Figure 5 that interaction started about time t
=
10, overlapping processes occurred between time t
=
15 and t
=
40, and waves started to resume their original shapes after the time t
=
40.Table 7
Comparison of invariants for the interaction of three solitary waves with results from [30] with h
=
0.2, k
=
0.025 in the region 0
≤
x
≤
250.
Present method
[30]
t
I
1
I
2
I
3
CPU Time
I
1
I
2
I
3
0
14.9800
15.8374
23.0081
—
13.6891
15.4549
22.8816
5
14.9804
15.7991
23.0465
0.536 s
13.6891
15.3109
22.6939
10
14.9783
14.3275
24.5180
1.442 s
13.6891
15.6514
22.8388
15
14.9809
15.5048
23.3408
2.349 s
13.6891
15.6548
22.9347
20
14.9814
15.7504
23.0952
3.251 s
13.6891
15.6557
22.9330
25
14.9819
15.8496
22.9960
4.162 s
13.6892
15.6559
22.9336
30
14.9823
15.8626
22.9829
5.063 s
13.6894
15.6559
22.9348
35
14.9828
15.8622
22.9833
5.573 s
13.6913
15.6564
22.9343
40
14.9832
15.8605
22.9851
6.481 s
13.7015
15.6566
22.9335
45
14.9461
15.8613
22.9843
7.388 s
13.7043
15.6563
22.9303
45 [34]
13.7043
15.6563
22.9303
45 [35] MQ
14.9681
15.7343
22.5966
45 [35] IMQ
14.9680
15.7343
22.5966
45 [35] IQ
14.9681
15.7343
22.5966
45 [35] GA
14.9681
15.7343
22.5966
45 [35] TPS
14.9682
15.7337
22.5944
45 [36]
14.7145
15.4927
23.3529Figure 5
Interaction of three solitary waves witht
=
0
,
5
,
8
,
15
,
20
,
and
40.


## 3.4. The Maxwellian Initial Condition

As our last problem, the development of the Maxwellian initial condition,(37)
U
(
x
,
0
)
=
exp
⁡
(
-
(
x
-
40
)
2
)
,
into a train of solitary waves is considered. For the Maxwellian initial condition, behavior of the solution depends on the values of μ. So we take μ
=
0.1
,
μ
=
0.04
,
μ
=
0.015, and μ
=
0.01. The numerical values of the invariants quantities during the simulations are given in Table 8. For μ
=
0.1, only a single soliton occurred as depicted in Figure 6(a). When μ
=
0.015 and μ
=
0.01, two and three stable solitons occurred, respectively, as depicted in Figures 6(b) and 6(c). For μ
=
0.04, Maxwellian initial condition has decayed into four solitary waves as depicted in Figure 6(d). As is seen, when μ is reduced, more solitary waves occurred. All figures are drawn up at time t
=
14.5. The peaks of the well developed wave lie on a straight line so that their velocities are linearly dependent on their amplitudes and also we observe a small oscillating tail appearing behind the last wave in all Maxwellian figures.Table 8
Invariants of MRLW equation using the Maxwellian initial condition.
t
μ
I
1
I
2
I
3
μ
I
1
I
2
I
3
0
1.7724809
1.3786628
0.7609109
1.7724809
1.2721327
0.8674410
3
1.7725930
1.4731551
0.6664186
1.7753954
1.4487372
0.6908365
6
0.1
1.7725929
1.4727329
0.6668408
0.015
1.7770146
1.4518203
0.6877534
9
1.7725935
1.4725552
0.6670185
1.7786793
1.4534252
0.6861485
12
1.7725942
1.4724619
0.6671118
1.7804784
1.4552349
0.6843388
15
1.77259499
1.4724042
0.6671695
1.7824228
1.4573443
0.6822294
0
1.7724809
1.2658662
0.8737075
1.7724809
1.3034651
0.8361087
3
1.7796208
1.4490883
0.6904854
1.7730165
1.4543782
0.6851955
6
0.01
1.7863091
1.4571676
0.6824061
0.04
1.7730562
1.4540759
0.6854978
9
1.7952456
1.4658449
0.6737288
1.7730962
1.4538204
0.6857534
12
1.8089221
1.4802019
0.6593718
1.7731365
1.4537359
0.6858378
15
1.8364987
1.5120445
0.6275292
1.7731770
1.4537234
0.6858503Maxwellian initial condition att
=
14.5 with (a) μ
=
0.1, (b) μ
=
0.015, (c) μ
=
0.01, and (d) μ
=
0.04.
(a)
(b)
(c)
(d)


## 4. Conclusion

In this paper, a numerical treatment of the MRLW equation has been introduced using septic B-spline collocation finite element method. The nonlinear term in the equation has been linearized by using a form given in the paper [37]. To examine the performance of the scheme, four test problems have been studied. The performance and accuracy of the method have been tested by calculating the error norms L
2 and L
∞ and the invariant quantities I
1
,
I
2, and I
3. Linear stability analysis proved that the present scheme is unconditionally stable. The experimental results of the algorithm are much satisfactory in comparison with the previous results. Our method can successfully be used to model the motion and interaction of the solitary waves. Thus, we can assert that our scheme is efficient and reliable for obtaining the numerical solutions of the other physically important nonlinear partial differential equations.



---

*Source: 596406-2014-07-22.xml*