# A Galerkin Finite Element Method for Numerical Solutions of the Modified Regularized Long Wave Equation

Authors: Liquan Mei; Yali Gao; Zhangxin Chen
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A Galerkin method for a modified regularized long wave equation is studied using finite elements in space, the Crank-Nicolson scheme, and the Runge-Kutta scheme in time. In addition, an extrapolation technique is used to transform a nonlinear system into a linear system in order to improve the time accuracy of this method. A Fourier stability analysis for the method is shown to be marginally stable. Three invariants of motion are investigated. Numerical experiments are presented to check the theoretical study of this method.


---

## Body

## 1. Introduction

Much physical phenomena are described by nonlinear partial differential equations. Most of these equations do not have an analytical solution, or it is extremely difficult and expensive to compute their analytical solutions. Hence a numerical study of these nonlinear partial differential equations is important in practice. The regularized long wave (RLW) equation(1)
u
t
+
u
x
+
u
u
x
-
μ
u
x
x
t
=
0
,
where μ is a positive constant, is a nonlinear evolution equation, which was originally introduced by Peregrine [1] in describing the behavior of an undular bore and studied later by Benjamin et al. [2]. This equation plays an important role in describing physical phenomena in various disciplines, such as the nonlinear transverse waves in shallow water, ion-acoustic waves in plasma, magnetohydrodynamics waves in plasma, longitudinal dispersive waves in elastic rods, and pressure waves in liquid’s gas bubbles. Many numerical methods for the RLW equation have been proposed, such as finite difference methods [3, 4], the Galerkin finite element method [5–8], the least squares method [9–11], various collocation methods with quadratic B-splines [12], cubic B-splines [13] and septic splines [14], meshfree method [15, 16], and an explicit multistep method [17].The RLW equation is a special case of the generalized regularized long wave (GRLW) equation(2)
u
t
+
u
x
+
δ
u
p
u
x
-
μ
u
x
x
t
=
0
,
where δ and μ are positive constants and p is a positive integer. Some numerical methods [18–24] for the GRLW equation have been also presented, such as a finite difference method [18], a decomposition method [20], and a sinc-collocation method [23]. Another special case of the GRLW equation is called the modified regularized long wave (MRLW) equation in which p
=
2 [25]. Some authors have studied the MRLW equation using various numerical methods, such as a finite difference method [26], a collocation method [27], a spline method [28, 29], and the Adomian decomposition method [30].In this paper, we study a Galerkin method for the MRLW equation by using linear finite elements in space and extrapolation to remove the nonlinear term. We discuss the properties of this method and compare its accuracy with previous studies. The interaction of two and three solitons is also studied. Moreover, the propagation of the Maxwellian initial condition is simulated.The outline of this paper is as follows. In the next section, the governing equation, its analytical solution, and three invariants are given. In Section3, we propose the numerical method, including semidiscretization in space and full discretization in space and time. In Section 4, we present a stability analysis. In Section 5, some numerical experiments are presented. Finally, we give a brief conclusion in Section 6.


## 2. Governing Equation

We assumeΩ
=
(
a
,
b
) to be an open, bounded subset of R and T
>
0 to be a final time. We set Ω
T
=
Ω
×
(
0
,
T
), J
=
(
0
,
T
] and consider the following MRLW equation:
(3)
u
t
+
u
x
+
6
u
2
u
x
-
μ
u
x
x
t
=
0
,
in
Ω
T
,
u
=
0
,
on
⁡
∂
Ω
×
[
0
,
T
]
,
u
=
u
0
,
on
⁡
Ω
×
{
t
=
0
}
,
where u
0
:
Ω
→
R is the initial datum, μ is a positive constant, and the subscripts x and t denote differentiation in space and time, respectively. The physical boundary condition requires u
→
0 as x
→
±
∞.The MRLW equation (3) has the exact solution [25]:
(4)
u
(
x
,
t
)
=
c
sech
(
p
(
x
-
(
c
+
1
)
t
-
x
0
)
)
,
where p
=
c
/
μ
(
c
+
1
) and x
0 and c are arbitrary constants. Furthermore, (3) possesses three invariants of motion corresponding to conservation of mass, momentum, and energy [25]:
(5)
I
1
=
∫
a
b
‍
u
d
x
,
I
2
=
∫
a
b
‍
(
u
2
+
μ
u
x
2
)
d
x
,
I
3
=
∫
a
b
‍
(
u
4
-
μ
u
x
2
)
d
x
.
These invariants are used to check the conservative properties of a numerical method, especially for problems without an analytical solution and during collision of solitons.


## 3. Numerical Methods

SetV
=
H
0
1
(
Ω
)
=
{
v
∈
H
1
(
Ω
)
:
v
|
∂
Ω
=
0
}. Applying Green’s formula to problem (3) and using the boundary condition in the definition of V, we derive the variational form of problem (3). Find u
:
J
→
V such that
(6)
(
∂
u
∂
t
,
v
)
+
B
(
u
,
v
)
+
C
(
u
,
v
)
+
D
(
u
,
v
)
=
0
,
∀
v
∈
V
,
t
∈
J
,
(
u
(
x
,
0
)
,
v
)
=
(
u
0
(
x
)
,
v
)
,
∀
v
∈
V
,
x
∈
Ω
,
where
(7)
B
(
u
,
v
)
=
∫
a
b
‍
μ
u
x
t
·
v
x
d
x
,
C
(
u
,
v
)
=
∫
a
b
‍
u
x
·
v
d
x
,
D
(
u
,
v
)
=
∫
a
b
‍
6
u
2
u
x
·
v
d
x
.


### 3.1. Semidiscretization in Space

Let us consider a uniform1
D mesh with the mesh size h
=
x
i
+
1
-
x
i, i
=
0,1
,
…
,
N
-
1, which consists of N
+
1 points a
=
x
0
<
x
1
<
⋯
<
x
N
-
1
<
x
N
=
b. Then we obtain
(8)
∑
i
=
0
N
-
1
‍
(
∫
x
i
x
i
+
1
‍
∂
u
∂
t
·
v
d
x
+
∫
x
i
x
i
+
1
‍
μ
u
x
t
·
v
x
d
x
+
∫
x
i
x
i
+
1
‍
u
x
·
v
d
x
+
∫
x
i
x
i
+
1
‍
6
u
2
u
x
·
v
d
x
)
=
0
.
We make the transformation x
=
x
i
+
η
h, 0
≤
η
≤
1, i
=
0,1
,
…
,
N
-
1, in order to transfer an element into a standard interval. From (8), we derive the following equation:
(9)
∫
0
1
‍
∂
u
∂
t
·
v
d
η
+
μ
h
2
∫
0
1
‍
u
η
t
·
v
η
d
η
+
1
h
∫
0
1
‍
u
η
·
v
d
η
+
6
h
∫
0
1
‍
u
2
u
η
·
v
d
η
=
0
.Now, we define the finite dimensional subspaceV
h
⊂
V, V
h
=
span
⁡
{
L
1
,
L
2
}, where
(10)
L
1
=
η
,
L
2
=
1
-
η
are the linear basis functions on each element. Then the semidiscrete scheme for problem (3) is formulated as follows. Find u
h
:
J
→
V
h such that
(11)
(
∂
u
h
∂
t
,
v
)
+
B
(
u
h
,
v
)
+
C
(
u
h
,
v
)
+
D
(
u
h
,
v
)
=
0
,
(
∂
u
h
∂
t
,
v
)
+
B
(
u
h
,
v
)
+
C
(
u
h
,
v
)
5
∀
v
∈
V
h
,
t
∈
J
,
(
u
h
(
x
,
0
)
,
v
)
=
(
u
0
h
(
x
)
,
v
)
,
∀
v
∈
V
h
,
x
∈
Ω
.The variation ofu over the element [
x
i
,
x
i
+
1
], i
=
0,1
,
…
,
N
-
1 is expressed as
(12)
u
e
=
∑
j
=
1
2
‍
L
j
(
x
)
u
j
(
t
)
.Forj
=
1,2, we choose v
=
L
j in problem (11) and substitute (12) into (11). Then an element’s contribution is obtained in the form of
(13)
∑
i
=
1
2
‍
{
(
∑
i
=
1
2
‍
u
i
L
i
)
2
(
∫
0
1
‍
L
i
·
L
j
d
x
+
μ
h
2
∫
0
1
‍
L
i
′
·
L
j
′
d
x
)
∂
u
i
∂
t
+
1
h
(
∫
0
1
‍
L
i
′
·
L
j
d
x
)
u
i
+
(
6
h
∫
0
1
‍
(
∑
i
=
1
2
‍
u
i
L
i
)
2
L
i
′
·
L
j
d
x
)
u
i
}
=
0
,
where the symbol (′) denotes differentiation with respect to η, which, in matrix form, is given by
(14)
(
A
e
+
B
e
)
∂
u
e
∂
t
+
C
e
u
e
+
D
e
(
u
e
)
u
e
=
0
,
where u
e
=
(
u
1
,
u
2
)
T are relevant nodal parameters. The element matrices are
(15)
A
j
k
e
=
∫
0
1
‍
L
j
L
k
d
η
,
B
j
k
e
=
μ
h
2
∫
0
1
‍
L
j
′
L
k
′
d
η
,
C
j
k
e
=
1
h
∫
0
1
‍
L
j
′
L
k
d
η
,
D
j
k
e
=
6
h
∫
0
1
‍
(
∑
m
=
1
2
‍
u
m
L
m
)
2
L
j
′
L
k
d
η
.For the element[
x
i
,
x
i
+
1
]
(
i
=
0,1
,
…
,
N
-
1
), the indices j and k take only the values 1 and 2 so that the matrices A
e, B
e, C
e, and D
e are 2 × 2:
(16)
A
e
=
(
1
3
1
6
1
6
1
3
)
,
B
j
k
e
=
μ
h
2
(
1
-
1
-
1
1
)
,
C
j
k
e
=
1
h
(
-
1
2
1
2
-
1
2
1
2
)
,
D
j
k
e
=
6
h
×
(
-
(
1
4
u
1
2
+
1
6
u
1
u
2
+
1
12
u
2
2
)
1
4
u
1
2
+
1
6
u
1
u
2
+
1
12
u
2
2
-
(
1
12
u
1
2
+
1
6
u
1
u
2
+
1
4
u
2
2
)
1
12
u
1
2
+
1
6
u
1
u
2
+
1
4
u
2
2
)
.Assembling contributions from all elements leads to the following matrix form of the coupled nonlinear ordinary differential equations:(17)
(
A
+
B
)
∂
u
∂
t
+
C
u
+
D
(
u
)
u
=
0
,
where u
=
(
u
0
,
u
1
,
…
,
u
N
)
T contains all the nodal parameters. The four assembled matrices are tridiagonal. The general row for each matrix has the following form:
(18)
A
:
(
1
6
,
2
3
,
1
6
)
,
B
:
μ
h
2
(
-
1,2
,
-
1
)
,
C
:
1
2
h
(
-
1,0
,
1
)
,
D
:
6
h
(
-
1
12
(
u
m
-
1
2
+
2
u
m
-
1
u
m
+
3
u
m
2
)
,
1
12
(
u
m
-
1
2
-
u
m
+
1
2
)
+
1
6
u
m
(
u
m
-
1
-
u
m
+
1
)
,
1
12
(
3
u
m
2
+
2
u
m
u
m
+
1
+
u
m
+
1
2
)
)
.


### 3.2. Full Discretization

We now consider a fully discrete scheme. Let0
=
t
0
<
t
1
<
⋯
<
t
M
=
T be a uniform partition of J into subintervals J
n
=
(
t
n
,
t
n
+
1
), n
=
0,1
,
…
,
M
-
1, with length Δ
t
=
t
n
+
1
-
t
n
=
T
/
M. For a generic function v of time, set v
n
=
v
(
t
n
). We use the Crank-Nicolson discretization method in (17) with
(19)
∂
u
∂
t
=
u
n
-
u
n
-
1
Δ
t
,
u
^
n
=
u
n
+
u
n
-
1
2
.
Then (17) can be written as
(20)
(
A
+
B
)
∂
u
∂
t
+
C
u
^
n
+
D
(
u
^
n
)
u
^
n
=
0
,
with u
0
=
u
0
h. This equation is symmetric about the point t
=
t
n
-
1
/
2, and one should, therefore, expect second order accuracy in time. However, the linearization decreases the order of the time discretization error to O
(
Δ
t
). This drawback can be overcome by using an extrapolation technique in the linearization of the nonlinear coefficient D. We choose u
¯
n
=
(
3
/
2
)
u
n
-
1
-
(
1
/
2
)
u
n
-
2, for n
≥
2, and define
(21)
(
A
+
B
)
∂
u
∂
t
+
C
u
^
n
+
D
(
u
¯
n
)
u
^
n
=
0
.
Now, the Crank-Nicolson method can be shown to produce an error of order O
(
(
Δ
t
)
2
). Combining with the error of finite element method, the Crank-Nicolson Galerkin finite element method of (6) produces an error of order O
(
(
Δ
t
)
2
+
h
2
) overall [17].We can obtain the following recurrence relation at pointx
m:
(22)
[
1
6
-
b
-
Δ
t
4
h
-
6
Δ
t
h
(
1
24
v
¯
m
-
1
2
+
1
12
v
¯
m
-
1
v
¯
m
+
1
8
v
¯
m
2
)
]
u
m
-
1
n
+
[
2
3
+
2
b
+
6
Δ
t
h
(
1
24
(
v
¯
m
-
1
2
-
v
¯
m
+
1
2
)
+
1
12
(
v
¯
m
-
1
-
v
¯
m
+
1
)
v
¯
m
)
]
u
m
n
+
[
1
6
-
b
+
Δ
t
4
h
+
6
Δ
t
h
(
1
8
v
¯
m
2
+
1
12
v
¯
m
v
¯
m
+
1
+
1
24
v
¯
m
+
1
2
)
]
u
m
+
1
n
=
[
1
6
-
b
+
Δ
t
4
h
+
6
Δ
t
h
(
1
24
v
¯
m
-
1
2
+
1
12
v
¯
m
-
1
v
¯
m
+
1
8
v
¯
m
2
)
]
u
m
-
1
n
-
1
+
[
2
3
+
2
b
-
6
Δ
t
h
(
1
24
(
v
¯
m
-
1
2
-
v
¯
m
+
1
2
)
+
1
12
(
v
¯
m
-
1
-
v
¯
m
+
1
)
v
¯
m
)
]
u
m
n
-
1
+
[
1
6
-
b
-
Δ
t
4
h
-
6
Δ
t
h
(
1
8
v
¯
m
2
+
1
12
v
¯
m
v
¯
m
+
1
+
1
24
v
¯
m
+
1
2
)
]
u
m
+
1
n
-
1
,
where v
¯
m and b are given by
(23)
v
¯
m
=
3
2
u
m
n
-
1
-
1
2
u
m
n
-
2
,
b
=
μ
h
2
.As (22) is used only for n
≥
2, we must compute u
1 by another method. For this, we choose a predictor-corrector method. Let n
=
1 in (21); we can get the first approximation u
1,0 with u
¯
1 replaced by u
0. Then we use (
1
/
2
)
(
u
1,0
+
u
0
) to substitute u
¯
1 in (21) and regard the result as the final approximation u
1. Or, more precisely, the concrete procedure is as follows. First, set
(24)
u
0
=
u
h
0
;
then we calculate
(25)
(
A
+
B
)
u
1,0
-
u
0
Δ
t
+
C
u
1,0
+
u
0
2
+
D
(
u
0
)
u
1,0
+
u
0
2
=
0
;
finally, we have
(26)
(
A
+
B
)
u
1
-
u
0
Δ
t
+
C
u
1
+
u
0
2
+
D
(
u
1,0
+
u
0
2
)
u
1
+
u
0
2
=
0
.Thus we can use the Thomas algorithm to solve the linear algebraic equations.In order to improve the accuracy in time, we also use Runge-Kutta discretization method. The fourth order Runge-Kutta discretization method in (17) can be described as follows:
(27)
u
n
+
1
=
u
n
+
K
1
+
2
K
2
+
2
K
3
+
K
4
6
,
-
(
A
+
B
)
K
1
=
[
C
u
n
+
D
(
u
n
)
u
n
]
Δ
t
,
-
(
A
+
B
)
K
2
=
[
C
(
u
n
+
1
2
K
1
)
+
D
(
u
n
+
1
2
K
1
)
(
u
n
+
1
2
K
1
)
]
Δ
t
,
-
(
A
+
B
)
K
3
=
[
C
(
u
n
+
1
2
K
2
)
+
D
(
u
n
+
1
2
K
2
)
(
u
n
+
1
2
K
2
)
]
Δ
t
,
-
(
A
+
B
)
K
4
=
[
C
(
u
n
+
K
3
)
+
D
(
u
n
+
K
3
)
(
u
n
+
K
3
)
]
Δ
t
.


## 3.1. Semidiscretization in Space

Let us consider a uniform1
D mesh with the mesh size h
=
x
i
+
1
-
x
i, i
=
0,1
,
…
,
N
-
1, which consists of N
+
1 points a
=
x
0
<
x
1
<
⋯
<
x
N
-
1
<
x
N
=
b. Then we obtain
(8)
∑
i
=
0
N
-
1
‍
(
∫
x
i
x
i
+
1
‍
∂
u
∂
t
·
v
d
x
+
∫
x
i
x
i
+
1
‍
μ
u
x
t
·
v
x
d
x
+
∫
x
i
x
i
+
1
‍
u
x
·
v
d
x
+
∫
x
i
x
i
+
1
‍
6
u
2
u
x
·
v
d
x
)
=
0
.
We make the transformation x
=
x
i
+
η
h, 0
≤
η
≤
1, i
=
0,1
,
…
,
N
-
1, in order to transfer an element into a standard interval. From (8), we derive the following equation:
(9)
∫
0
1
‍
∂
u
∂
t
·
v
d
η
+
μ
h
2
∫
0
1
‍
u
η
t
·
v
η
d
η
+
1
h
∫
0
1
‍
u
η
·
v
d
η
+
6
h
∫
0
1
‍
u
2
u
η
·
v
d
η
=
0
.Now, we define the finite dimensional subspaceV
h
⊂
V, V
h
=
span
⁡
{
L
1
,
L
2
}, where
(10)
L
1
=
η
,
L
2
=
1
-
η
are the linear basis functions on each element. Then the semidiscrete scheme for problem (3) is formulated as follows. Find u
h
:
J
→
V
h such that
(11)
(
∂
u
h
∂
t
,
v
)
+
B
(
u
h
,
v
)
+
C
(
u
h
,
v
)
+
D
(
u
h
,
v
)
=
0
,
(
∂
u
h
∂
t
,
v
)
+
B
(
u
h
,
v
)
+
C
(
u
h
,
v
)
5
∀
v
∈
V
h
,
t
∈
J
,
(
u
h
(
x
,
0
)
,
v
)
=
(
u
0
h
(
x
)
,
v
)
,
∀
v
∈
V
h
,
x
∈
Ω
.The variation ofu over the element [
x
i
,
x
i
+
1
], i
=
0,1
,
…
,
N
-
1 is expressed as
(12)
u
e
=
∑
j
=
1
2
‍
L
j
(
x
)
u
j
(
t
)
.Forj
=
1,2, we choose v
=
L
j in problem (11) and substitute (12) into (11). Then an element’s contribution is obtained in the form of
(13)
∑
i
=
1
2
‍
{
(
∑
i
=
1
2
‍
u
i
L
i
)
2
(
∫
0
1
‍
L
i
·
L
j
d
x
+
μ
h
2
∫
0
1
‍
L
i
′
·
L
j
′
d
x
)
∂
u
i
∂
t
+
1
h
(
∫
0
1
‍
L
i
′
·
L
j
d
x
)
u
i
+
(
6
h
∫
0
1
‍
(
∑
i
=
1
2
‍
u
i
L
i
)
2
L
i
′
·
L
j
d
x
)
u
i
}
=
0
,
where the symbol (′) denotes differentiation with respect to η, which, in matrix form, is given by
(14)
(
A
e
+
B
e
)
∂
u
e
∂
t
+
C
e
u
e
+
D
e
(
u
e
)
u
e
=
0
,
where u
e
=
(
u
1
,
u
2
)
T are relevant nodal parameters. The element matrices are
(15)
A
j
k
e
=
∫
0
1
‍
L
j
L
k
d
η
,
B
j
k
e
=
μ
h
2
∫
0
1
‍
L
j
′
L
k
′
d
η
,
C
j
k
e
=
1
h
∫
0
1
‍
L
j
′
L
k
d
η
,
D
j
k
e
=
6
h
∫
0
1
‍
(
∑
m
=
1
2
‍
u
m
L
m
)
2
L
j
′
L
k
d
η
.For the element[
x
i
,
x
i
+
1
]
(
i
=
0,1
,
…
,
N
-
1
), the indices j and k take only the values 1 and 2 so that the matrices A
e, B
e, C
e, and D
e are 2 × 2:
(16)
A
e
=
(
1
3
1
6
1
6
1
3
)
,
B
j
k
e
=
μ
h
2
(
1
-
1
-
1
1
)
,
C
j
k
e
=
1
h
(
-
1
2
1
2
-
1
2
1
2
)
,
D
j
k
e
=
6
h
×
(
-
(
1
4
u
1
2
+
1
6
u
1
u
2
+
1
12
u
2
2
)
1
4
u
1
2
+
1
6
u
1
u
2
+
1
12
u
2
2
-
(
1
12
u
1
2
+
1
6
u
1
u
2
+
1
4
u
2
2
)
1
12
u
1
2
+
1
6
u
1
u
2
+
1
4
u
2
2
)
.Assembling contributions from all elements leads to the following matrix form of the coupled nonlinear ordinary differential equations:(17)
(
A
+
B
)
∂
u
∂
t
+
C
u
+
D
(
u
)
u
=
0
,
where u
=
(
u
0
,
u
1
,
…
,
u
N
)
T contains all the nodal parameters. The four assembled matrices are tridiagonal. The general row for each matrix has the following form:
(18)
A
:
(
1
6
,
2
3
,
1
6
)
,
B
:
μ
h
2
(
-
1,2
,
-
1
)
,
C
:
1
2
h
(
-
1,0
,
1
)
,
D
:
6
h
(
-
1
12
(
u
m
-
1
2
+
2
u
m
-
1
u
m
+
3
u
m
2
)
,
1
12
(
u
m
-
1
2
-
u
m
+
1
2
)
+
1
6
u
m
(
u
m
-
1
-
u
m
+
1
)
,
1
12
(
3
u
m
2
+
2
u
m
u
m
+
1
+
u
m
+
1
2
)
)
.


## 3.2. Full Discretization

We now consider a fully discrete scheme. Let0
=
t
0
<
t
1
<
⋯
<
t
M
=
T be a uniform partition of J into subintervals J
n
=
(
t
n
,
t
n
+
1
), n
=
0,1
,
…
,
M
-
1, with length Δ
t
=
t
n
+
1
-
t
n
=
T
/
M. For a generic function v of time, set v
n
=
v
(
t
n
). We use the Crank-Nicolson discretization method in (17) with
(19)
∂
u
∂
t
=
u
n
-
u
n
-
1
Δ
t
,
u
^
n
=
u
n
+
u
n
-
1
2
.
Then (17) can be written as
(20)
(
A
+
B
)
∂
u
∂
t
+
C
u
^
n
+
D
(
u
^
n
)
u
^
n
=
0
,
with u
0
=
u
0
h. This equation is symmetric about the point t
=
t
n
-
1
/
2, and one should, therefore, expect second order accuracy in time. However, the linearization decreases the order of the time discretization error to O
(
Δ
t
). This drawback can be overcome by using an extrapolation technique in the linearization of the nonlinear coefficient D. We choose u
¯
n
=
(
3
/
2
)
u
n
-
1
-
(
1
/
2
)
u
n
-
2, for n
≥
2, and define
(21)
(
A
+
B
)
∂
u
∂
t
+
C
u
^
n
+
D
(
u
¯
n
)
u
^
n
=
0
.
Now, the Crank-Nicolson method can be shown to produce an error of order O
(
(
Δ
t
)
2
). Combining with the error of finite element method, the Crank-Nicolson Galerkin finite element method of (6) produces an error of order O
(
(
Δ
t
)
2
+
h
2
) overall [17].We can obtain the following recurrence relation at pointx
m:
(22)
[
1
6
-
b
-
Δ
t
4
h
-
6
Δ
t
h
(
1
24
v
¯
m
-
1
2
+
1
12
v
¯
m
-
1
v
¯
m
+
1
8
v
¯
m
2
)
]
u
m
-
1
n
+
[
2
3
+
2
b
+
6
Δ
t
h
(
1
24
(
v
¯
m
-
1
2
-
v
¯
m
+
1
2
)
+
1
12
(
v
¯
m
-
1
-
v
¯
m
+
1
)
v
¯
m
)
]
u
m
n
+
[
1
6
-
b
+
Δ
t
4
h
+
6
Δ
t
h
(
1
8
v
¯
m
2
+
1
12
v
¯
m
v
¯
m
+
1
+
1
24
v
¯
m
+
1
2
)
]
u
m
+
1
n
=
[
1
6
-
b
+
Δ
t
4
h
+
6
Δ
t
h
(
1
24
v
¯
m
-
1
2
+
1
12
v
¯
m
-
1
v
¯
m
+
1
8
v
¯
m
2
)
]
u
m
-
1
n
-
1
+
[
2
3
+
2
b
-
6
Δ
t
h
(
1
24
(
v
¯
m
-
1
2
-
v
¯
m
+
1
2
)
+
1
12
(
v
¯
m
-
1
-
v
¯
m
+
1
)
v
¯
m
)
]
u
m
n
-
1
+
[
1
6
-
b
-
Δ
t
4
h
-
6
Δ
t
h
(
1
8
v
¯
m
2
+
1
12
v
¯
m
v
¯
m
+
1
+
1
24
v
¯
m
+
1
2
)
]
u
m
+
1
n
-
1
,
where v
¯
m and b are given by
(23)
v
¯
m
=
3
2
u
m
n
-
1
-
1
2
u
m
n
-
2
,
b
=
μ
h
2
.As (22) is used only for n
≥
2, we must compute u
1 by another method. For this, we choose a predictor-corrector method. Let n
=
1 in (21); we can get the first approximation u
1,0 with u
¯
1 replaced by u
0. Then we use (
1
/
2
)
(
u
1,0
+
u
0
) to substitute u
¯
1 in (21) and regard the result as the final approximation u
1. Or, more precisely, the concrete procedure is as follows. First, set
(24)
u
0
=
u
h
0
;
then we calculate
(25)
(
A
+
B
)
u
1,0
-
u
0
Δ
t
+
C
u
1,0
+
u
0
2
+
D
(
u
0
)
u
1,0
+
u
0
2
=
0
;
finally, we have
(26)
(
A
+
B
)
u
1
-
u
0
Δ
t
+
C
u
1
+
u
0
2
+
D
(
u
1,0
+
u
0
2
)
u
1
+
u
0
2
=
0
.Thus we can use the Thomas algorithm to solve the linear algebraic equations.In order to improve the accuracy in time, we also use Runge-Kutta discretization method. The fourth order Runge-Kutta discretization method in (17) can be described as follows:
(27)
u
n
+
1
=
u
n
+
K
1
+
2
K
2
+
2
K
3
+
K
4
6
,
-
(
A
+
B
)
K
1
=
[
C
u
n
+
D
(
u
n
)
u
n
]
Δ
t
,
-
(
A
+
B
)
K
2
=
[
C
(
u
n
+
1
2
K
1
)
+
D
(
u
n
+
1
2
K
1
)
(
u
n
+
1
2
K
1
)
]
Δ
t
,
-
(
A
+
B
)
K
3
=
[
C
(
u
n
+
1
2
K
2
)
+
D
(
u
n
+
1
2
K
2
)
(
u
n
+
1
2
K
2
)
]
Δ
t
,
-
(
A
+
B
)
K
4
=
[
C
(
u
n
+
K
3
)
+
D
(
u
n
+
K
3
)
(
u
n
+
K
3
)
]
Δ
t
.


## 4. Stability Analysis

A linear stability analysis is made of the growth factorg of the error ε
j
n in a typical Fourier mode of amplitude ε
^
n:
(28)
ε
j
n
=
e
i
j
k
h
ε
^
n
,
where k is a mode number and h is the element size. The von Neumann stability method can only be applied to a linear scheme, so we linearize the nonlinear term by assuming that u in this term is locally constant. Substituting (28) into (22) gives
(29)
ε
^
n
+
1
=
g
ε
^
n
,
(30)
g
=
α
+
β
i
α
-
β
i
,
where
(31)
α
=
(
1
3
-
2
b
)
cos
⁡
(
k
h
)
+
2
3
+
2
b
,
β
=
Δ
t
2
h
(
1
+
ɛ
C
)
sin
(
k
h
)
.
Taking the modulus of (30) gives
(32)
|
g
|
=
g
g
¯
=
1
.
Therefore, our method is marginally stable.


## 5. Numerical Experiments

In this section, we present some numerical tests to check the efficiency and accuracy of our method, which are the propagation of single soliton and collision of two and three solitons at different time levels. Finally, we investigate the development of the Maxwellian initial condition into solitary waves.To illustrate the accuracy of the present method, we useL
2- and L
∞-error norms
(33)
L
2
=
∥
u
exact
,
n
-
u
n
∥
2
=
(
h
∑
i
=
1
N
‍
|
u
j
exact
,
n
-
u
j
n
|
2
)
1
/
2
,
L
∞
=
∥
u
exact
,
n
-
u
n
∥
∞
=
max
⁡
j
|
u
j
exact
,
n
-
u
j
n
|
to compare the numerical solution with the exact solution, which show the mean and maximum differences between the numerical and analytical solutions. The quantities I
1, I
2, and I
3 measure the conservation laws of our method during propagation.


### 5.1. Single Soliton

The analytical values of the three variants are(34)
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 purpose of comparison with the previous work [28], we choose the parameters x
0
=
40, c
=
0.05, h
=
k
=
0.2, T
=
10, and 0
≤
x
≤
100. Then the amplitude is 0.7071068 and the analytical values of three invariants are I
1
=
3.219174470, I
2
=
0.465531499, and I
3
=
0.008001323. The relevant numerical results which applied the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27) are presented in Tables 1 and 2, respectively. The profiles of solitary waves at time = 5 and time = 10 are given in Figure 1 and the error distribution in the wave profile in the three-dimensional view at different time level is illustrated in Figure 2. Table 3 presents the variants, L
∞ norm, and L
2 norm at time = 10 against the quadratic B-splines and the cubic B-splines in [28].Table 1
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, and 0
≤
x
≤
100 (using Crank-Nicolson time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
0
3.218839
0.4655250
8.007784E − 03
0
0
2
3.218839
0.4655233
8.007697E − 03
0.83581
E
-
04
0.45779
E
-
04
4
3.218854
0.4655215
8.007602E − 03
1.36403
E
-
04
0.45846
E
-
04
6
3.218842
0.4655196
8.007507E − 03
1.92290
E
-
04
0.68890
E
-
04
8
3.218819
0.4655177
8.007410E − 03
2.47879
E
-
04
0.91724
E
-
04
10
3.218794
0.4655159
8.007313E − 03
3.02940
E
-
04
1.13929
E
-
04Table 2
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, and 0
≤
x
≤
100 (using Runge-Kutta time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
0
3.218839
0.4655250
8.007778E − 03
0
0
2
3.218977
0.4655250
8.007779E − 03
2.613735
E
-
05
4.577868
E
-
05
4
3.219083
0.4655250
8.007780E − 03
3.432730
E
-
05
2.894973
E
-
05
6
3.219165
0.4655250
8.007780E − 03
4.744950
E
-
05
2.214867
E
-
05
8
3.219229
0.4655250
8.007780E − 03
6.052811
E
-
05
2.298029
E
-
05
10
3.219274
0.4655250
8.007780E − 03
7.200630
E
-
05
2.346961
E
-
05Table 3
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, 0
≤
x
≤
100, and time = 10.
Method
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
Our scheme in (22)
3.218794
0.4655160
8.007313E − 03
3.02940
E
-
04
1.13929
E
-
0
4
Our scheme in (27)
3.219274
0.4655250
8.007780E − 03
0.720063
E
-
0
4
0.234696
E
-
0
4
Quadratic [28]
3.215653
0.4655665
8.004883E − 03
1.99288
E
-
0
4
4.81156
E
-
0
4
Cubic [28]
3.215189
0.4655136
7.999173E − 03
4.53811
E
-
0
4
6.25887
E
-
0
4Figure 1
Solitary wave withc
=
0.05, h
=
0.2, k
=
0.2 and x
0
=
40, 0
≤
x
≤
100.Figure 2
Error graph forc
=
0.05, h
=
0.2, k
=
0.2 and x
0
=
40
,
0
≤
x
≤
100.We also choose the parametersc
=
0.3, h
=
0.1, k
=
0.01, and x
0
=
40, with the range [
0,100
]. Then the amplitude is 0.54772 and the analytical values of three invariants are I
1
=
3.581966678, I
2
=
1.345076492, and I
3
=
0.153723028. Tables 4 and 5 show the results about the three invariants, the L
2 norm, and the L
∞ norm at different times and compare these results with those by the collocation method with cubic B-splines in [27] using the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27), respectively. The changes of the variants are satisfactorily small, though our result for this case is not more accurate than the result in [27], where the simulation is done up to t
=
10. The profile of the approximate solution is given in Figures 3 and 4 shows the error profile of the single solitary waves.Table 4
Invariants and error norms for single solitary wave,c
=
0.3, h
=
0.1, k
=
0.01, x
0
=
40, and 0
≤
x
≤
100 (using Crank-Nicolson time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.581967
1.345076
0.1537230
0
3.581967
1.344973
0.1538264
0
0
2
3.581967
1.344973
0.1538263
1.42310
E
-
04
0.45779
E
-
04
4
3.581966
1.344973
0.1538262
2.65355
E
-
04
1.53988
E
-
04
6
3.581966
1.344973
0.1538261
3.74673
E
-
04
2.00152
E
-
04
8
3.581966
1.344972
0.1538260
4.77620
E
-
04
2.42916
E
-
04
10
3.581966
1.344972
0.1538260
5.77512
E
-
04
2.84373
E
-
04
[27]
3.58197
1.34508
0.153723
4.02927
E
-
04
2.06732
E
-
04Table 5
Invariants and error norms for single solitary wave,c
=
0.3, h
=
0.1, k
=
0.01, x
0
=
40, and 0
≤
x
≤
100 (using Runge-Kutta time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.581967
1.345076
0.1537230
0
3.581967
1.344973
0.1538264
0
0
2
3.581967
1.344973
0.1538263
1.396685
E
-
04
0.920780
E
-
04
4
3.581967
1.344973
0.1538262
2.592545
E
-
04
1.512530
E
-
04
6
3.581967
1.344973
0.1538262
3.647701
E
-
04
1.956752
E
-
04
8
3.581967
1.344973
0.1538262
4.638980
E
-
04
2.368004
E
-
04
10
3.581967
1.344973
0.1538262
5.600205
E
-
04
2.766940
E
-
04
[27]
3.58197
1.34508
0.153723
4.02927
E
-
04
2.06732
E
-
04Figure 3
Solitary wave withc
=
0.3, h
=
0.1, k
=
0.01 and x
0
=
40
,
0
≤
x
≤
100.Figure 4
Error graph forc
=
0.3, h
=
0.1, k
=
0.01 and x
0
=
40
,
0
≤
x
≤
100.The Runge-Kutta Galerkin method algorithm has been run for different space steps with a fixed time stepk
=
0.01. The error and convergence order of space for the L
2- and L
∞-norm are recorded in Tables 6 and 7, respectively, which verify that the L
2- and L
∞-error can achieve the order O
(
h
2
) in space. And the algorithm has also been run for different time steps with a fixed space step h
=
0.1. The error and convergence order of time for the L
2- and L
∞-norm are recorded in Tables 8 and 9, respectively.Table 6
The order of convergence in space forL
2 error norm with k
=
0.01 and c
=
0.3.
Space
Δh
=
0.1
Δh
=
0.2
Δh
=
0.4
Order
Order
(0.1/0.2)
(0.2/0.4)
2
1.396685 × 10−4
5.598770 × 10−4
2.258793 × 10−3
2.0031
2.0124
4
2.592545 × 10−4
1.038054 × 10−3
4.168527 × 10−3
2.0014
2.0057
6
3.647701 × 10−4
1.459612 × 10−3
5.846531 × 10−3
2.0005
2.0020
8
4.638980 × 10−4
1.855582 × 10−3
7.421648 × 10−3
2.0000
1.9999
10
5.600205 × 10−4
2.239531 × 10−3
8.948683 × 10−3
1.9996
1.9985Table 7
The order of convergence in space forL
∞ error norm with k
=
0.01 and c
=
0.3.
Space
Δh
=
0.1
Δh
=
0.2
Δt
=
0.4
Order
Order
(0.1/0.2)
(0.2/0.4)
2
0.920780 × 10−4
3.689931 × 10−4
1.480253 × 10−3
2.0027
2.0042
4
1.512530 × 10−4
6.055462 × 10−4
2.430139 × 10−3
2.0012
2.0047
6
1.956752 × 10−4
7.811765 × 10−4
3.126531 × 10−3
1.9972
2.0008
8
2.368004 × 10−4
9.454389 × 10−4
3.768796 × 10−3
1.9973
1.9950
10
2.766940 × 10−4
1.104758 × 10−3
4.416067 × 10−3
1.9974
1.9990Table 8
The order of convergence in time forL
2 error norm with h
=
0.1 and c
=
0.3.
Time
Δt
=
2.0
Δt
=
1.0
Δt
=
0.5
Order
Order
(2.0/1.0)
(1.0/0.5)
2
0.0511878
0.00418182
4.219760 × 10−4
3.613597
3.308899
4
0.0835505
0.00806197
8.133115 × 10−4
3.309253
3.309253
6
0.1027137
0.01176588
1.180928 × 10−3
3.125948
3.316617
8
0.1141693
0.01561415
1.547231 × 10−3
2.870249
3.335094
10
0.1220153
0.01982141
1.925846 × 10−3
2.621930
3.363496Table 9
The order of convergence in time forL
∞ error norm with h
=
0.1 and c
=
0.3.
Time
Δt
=
2.0
Δt
=
1.0
Δt
=
0.5
Order
Order
(2.0/1.0)
(1.0/0.5)
2
0.0318115
2.528470 × 10−3
2.505105 × 10−4
3.653213
3.335322
4
0.0498808
4.159059 × 10−3
4.456332 × 10−4
3.584155
3.222328
6
0.0594582
5.587773 × 10−3
5.905399 × 10−4
3.411531
3.242167
8
0.0641225
7.106547 × 10−3
7.270440 × 10−4
3.173611
3.289034
10
0.0661177
8.861734 × 10−3
8.671150 × 10−4
2.899376
3.353294


### 5.2.  Collision of Two Solitons

In this section, we study the interaction of two solitary waves with the initial conditions given by a linear sum of two separate solitary waves of various amplitudes:(35)
u
(
x
,
0
)
=
∑
i
=
1
2
‍
A
i
sech
⁡
(
p
i
(
x
-
x
i
)
)
,
where A
i
=
c
i, p
i
=
c
i
/
μ
(
c
i
+
1
), i
=
1
,
2, and x
i and c
i
(
i
=
1,2
) are arbitrary constants. The analytical values for the conservation laws in this case have the following form:
(36)
I
1
=
π
c
1
p
1
+
π
c
2
p
2
,
I
2
=
2
c
1
p
1
+
2
c
2
p
2
+
2
μ
c
1
p
1
3
+
2
μ
p
2
c
2
3
,
I
3
=
4
c
1
2
3
p
1
+
4
c
2
2
3
p
2
-
2
μ
p
1
c
1
3
-
2
μ
p
2
c
2
3
.
In our numerical scheme, we choose h
=
0.1, k
=
0.01, c
1
=
1, c
2
=
0.5, x
1
=
10, and x
2
=
40, with interval [
0,100
]. Then the amplitudes are in ratio 2 : 1, where 2
A
2
=
A
1 and the analytical values for the conservation are I
1
=
8.2905324294, I
2
=
5.224332548, and I
3
=
1.799113742. The changes of the invariants are satisfactorily small, since the changes of the invariants I
1, I
2, and I
3 are 1.3295 × 10−4, 2.2259 × 10−5, and 6.8989 × 10−5, respectively. The results are recorded in Table 10. A graph of the two solitary waves collision, plotted in a three-dimensional view at some discrete times, is shown in Figure 5.Table 10
Invariants and error norms for two solitary waves,c
1
=
1, c
2
=
0.5, x
1
=
10, x
2
=
40, h
=
0.1, k
=
0.01, and 0
≤
x
≤
100.
Time
I
1
I
2
I
3
C-N
R-K
C-N
R-K
C-N
R-K
0
8.288129
8.288129
5.222934
5.222933
1.800512
1.800513
2
8.289522
8.293010
5.222910
5.222951
1.800479
1.800507
4
8.289462
8.295348
5.222887
5.222954
1.800458
1.822507
6
8.289304
8.297822
5.222864
5.222956
1.800433
1.800507
8
8.289233
8.300322
5.222840
5.222958
1.800411
1.800508
10
8.289231
8.302630
5.222817
5.222960
1.800388
1.800509Figure 5
Interaction of two solitons at times fromt
=
2,4
,
6,8
,
10 for c
1
=
1, c
2
=
0.5, x
1
=
10, x
2
=
40, and 0
≤
x
≤
100.


### 5.3.  Collision of Three Solitons

In this section, we study the interaction of three solitary waves with the initial conditions given by a linear sum of three separate solitary waves of various amplitudes:(37)
u
(
x
,
0
)
=
∑
i
=
1
3
‍
A
i
sech
⁡
(
p
i
(
x
-
x
i
)
)
,
where A
i
=
c
i, p
i
=
c
i
/
μ
(
c
i
+
1
), i
=
1,2
,
3, and x
i and c
i
(
i
=
1,2
,
3
) are arbitrary constants. The analytical values for the conservation laws in this case have the following form:
(38)
I
1
=
π
c
1
p
1
+
π
c
2
p
2
+
π
c
3
p
3
,
I
2
=
2
c
1
p
1
+
2
c
2
p
2
+
2
c
3
p
3
+
2
μ
c
1
p
1
3
+
2
μ
p
2
c
2
3
+
2
μ
p
3
c
3
3
,
I
3
=
4
c
1
2
3
p
1
+
4
c
2
2
3
p
2
+
4
c
3
2
3
p
3
-
2
μ
p
1
c
1
3
-
2
μ
p
2
c
2
3
-
2
μ
p
3
c
3
3
.
In this case, we choose paraments as h
=
0.1, k
=
0.05, c
1
=
1, c
2
=
0.5, c
3
=
0.25, x
1
=
10, x
2
=
20, and x
3
=
40, with interval [
0,100
]. Then the amplitudes are in ratio 4 : 2 : 1, where A
1
=
2
A
2
=
4
A
3 and the analytical values for the conservation are I
1
=
11.802939794, I
2
=
6.416902131, and I
3
=
1.910917141. The changes of the invariants from the initial variants approach zero throughout and agree with the analytical values for the three invariants as presented in Table 11, which indicates that our scheme is satisfactorily conservative. The geometry of the initial state and the profiles at time t
=
10, t
=
20, and t
=
35 is shown graphically in Figure 6, respectively.Table 11
Invariants and error norms for three solitary waves,c
1
=
1, c
2
=
0.5, c
3
=
0.25, x
1
=
10, x
2
=
20, x
3
=
40, h
=
0.1, k
=
0.01, and 0
≤
x
≤
100.
Time
I
1
I
2
I
3
C-N
R-K
C-N
R-K
C-N
R-K
0
11.800513
11.800513
6.490318
6.490318
1.997082
1.997082
4
11.799364
11.807800
6.485747
6.490319
1.992012
1.996623
8
11.797850
11.812825
6.483120
6.490298
1.988853
1.996076
12
11.797570
11.817192
6.482677
6.490279
1.988441
1.996360
16
11.791837
11.820521
6.474939
6.490305
1.982852
1.996875
20
11.788577
11.823320
6.469884
6.490341
1.977926
1.996758Interaction of three solitary waves.
(a)
Time = 0
(b)
Time = 10
(c)
Time = 20
(d)
Time = 35


### 5.4.  The Maxwellian Initial Condition

Finally, we consider the development of the Maxwellian initial condition(39)
u
(
x
,
0
)
=
e
-
(
x
-
40
)
2
into a train of solitary waves. We choose different values μ
=
1.0, μ
=
0.5, and μ
=
0.1 in our numerical scheme. The comparison of the three variants I
1, I
2, and I
3 with the earlier result in [29] is presented in Table 12, here Crank-Nicolson method is used for time discretization. For μ
=
1.0, the changes of the variants I
1, I
2, and I
3 with respect to the initial values are 3.9492 × 10−6, 6.9575 × 10−4, and 3.1509 × 10−3, respectively, whereas they are 2.8210 × 10−6, 2.60000 × 10−4, and 7.0771 × 10−4 in [29]. For μ
=
0.5, the variants are changed by 7.8982 × 10−6, 6.8214 × 10−5, and 6.1410 × 10−4 in this case. The development of the Maxwellian initial condition is shown in Figure 7 with different parameter values, respectively. The smaller μ there is, the more the number of solitary waves will form. The simulations are done up to time = 10 in this case. From the previous work [25], the total number of solitary waves and the values μ have the approximate relation
(40)
N
≅
[
1
μ
5
]
.Table 12
Values ofI
1, I
2, and I
3 for Maxwellian initial condition when h
=
0.1 and k
=
0.01 and the space interval [
0,100
].
μ
Time
I
1
I
2
I
3
I
1 [29]
I
2 [29]
I
3 [29]
1
2
1.772528
2.499467
−0.3633570
1.772449
2.506352
−0.3668149
4
1.772544
2.499870
−0.3643380
1.772446
2.506235
−0.3666974
6
1.772532
2.500822
−0.3657774
1.772447
2.506171
−0.3666326
8
1.772524
2.500977
−0.3661331
1.772446
2.506123
−0.3665860
10
1.772521
2.501206
−0.3665079
1.772444
2.2060926
−0.3665553
0.5
2
1.772534
1.876446
0.2597450
1.772451
1.879888
0.2596494
4
1.772534
1.876534
0.2587968
1.772446
1.879855
0.2596828
6
1.772528
1.876600
0.2583802
1.772449
1.879841
0.2596973
8
1.772523
1.876581
0.2582386
1.772450
1.879834
0.2597050
10
1.772520
1.876574
0.2581499
1.772449
1.879828
0.2597092
0.1
2
1.770991
1.373320
0.7687222
1.772452
1.378607
0.7608777
4
1.769863
1.371236
0.7658462
1.772451
1.378577
0.7608364
6
1.768742
1.369139
0.7630269
1.772451
1.378546
0.7607937
8
1.767630
1.367058
0.7602389
1.772451
1.378515
0.7607529
10
1.766526
1.364993
0.7574795
1.772453
1.378483
0.7607117Maxwellian initial condition with different valuesμ.
(a)
μ
=
1.0
(b)
μ
=
0.5
(c)
μ
=
0.1


## 5.1. Single Soliton

The analytical values of the three variants are(34)
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 purpose of comparison with the previous work [28], we choose the parameters x
0
=
40, c
=
0.05, h
=
k
=
0.2, T
=
10, and 0
≤
x
≤
100. Then the amplitude is 0.7071068 and the analytical values of three invariants are I
1
=
3.219174470, I
2
=
0.465531499, and I
3
=
0.008001323. The relevant numerical results which applied the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27) are presented in Tables 1 and 2, respectively. The profiles of solitary waves at time = 5 and time = 10 are given in Figure 1 and the error distribution in the wave profile in the three-dimensional view at different time level is illustrated in Figure 2. Table 3 presents the variants, L
∞ norm, and L
2 norm at time = 10 against the quadratic B-splines and the cubic B-splines in [28].Table 1
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, and 0
≤
x
≤
100 (using Crank-Nicolson time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
0
3.218839
0.4655250
8.007784E − 03
0
0
2
3.218839
0.4655233
8.007697E − 03
0.83581
E
-
04
0.45779
E
-
04
4
3.218854
0.4655215
8.007602E − 03
1.36403
E
-
04
0.45846
E
-
04
6
3.218842
0.4655196
8.007507E − 03
1.92290
E
-
04
0.68890
E
-
04
8
3.218819
0.4655177
8.007410E − 03
2.47879
E
-
04
0.91724
E
-
04
10
3.218794
0.4655159
8.007313E − 03
3.02940
E
-
04
1.13929
E
-
04Table 2
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, and 0
≤
x
≤
100 (using Runge-Kutta time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
0
3.218839
0.4655250
8.007778E − 03
0
0
2
3.218977
0.4655250
8.007779E − 03
2.613735
E
-
05
4.577868
E
-
05
4
3.219083
0.4655250
8.007780E − 03
3.432730
E
-
05
2.894973
E
-
05
6
3.219165
0.4655250
8.007780E − 03
4.744950
E
-
05
2.214867
E
-
05
8
3.219229
0.4655250
8.007780E − 03
6.052811
E
-
05
2.298029
E
-
05
10
3.219274
0.4655250
8.007780E − 03
7.200630
E
-
05
2.346961
E
-
05Table 3
Invariants and error norms for single solitary wave,c
=
0.05, h
=
0.2, k
=
0.2, x
0
=
40, 0
≤
x
≤
100, and time = 10.
Method
I
1
I
2
I
3
L
2
L
∞
Analytical
3.219174
0.465531
8.001323E − 03
Our scheme in (22)
3.218794
0.4655160
8.007313E − 03
3.02940
E
-
04
1.13929
E
-
0
4
Our scheme in (27)
3.219274
0.4655250
8.007780E − 03
0.720063
E
-
0
4
0.234696
E
-
0
4
Quadratic [28]
3.215653
0.4655665
8.004883E − 03
1.99288
E
-
0
4
4.81156
E
-
0
4
Cubic [28]
3.215189
0.4655136
7.999173E − 03
4.53811
E
-
0
4
6.25887
E
-
0
4Figure 1
Solitary wave withc
=
0.05, h
=
0.2, k
=
0.2 and x
0
=
40, 0
≤
x
≤
100.Figure 2
Error graph forc
=
0.05, h
=
0.2, k
=
0.2 and x
0
=
40
,
0
≤
x
≤
100.We also choose the parametersc
=
0.3, h
=
0.1, k
=
0.01, and x
0
=
40, with the range [
0,100
]. Then the amplitude is 0.54772 and the analytical values of three invariants are I
1
=
3.581966678, I
2
=
1.345076492, and I
3
=
0.153723028. Tables 4 and 5 show the results about the three invariants, the L
2 norm, and the L
∞ norm at different times and compare these results with those by the collocation method with cubic B-splines in [27] using the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27), respectively. The changes of the variants are satisfactorily small, though our result for this case is not more accurate than the result in [27], where the simulation is done up to t
=
10. The profile of the approximate solution is given in Figures 3 and 4 shows the error profile of the single solitary waves.Table 4
Invariants and error norms for single solitary wave,c
=
0.3, h
=
0.1, k
=
0.01, x
0
=
40, and 0
≤
x
≤
100 (using Crank-Nicolson time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.581967
1.345076
0.1537230
0
3.581967
1.344973
0.1538264
0
0
2
3.581967
1.344973
0.1538263
1.42310
E
-
04
0.45779
E
-
04
4
3.581966
1.344973
0.1538262
2.65355
E
-
04
1.53988
E
-
04
6
3.581966
1.344973
0.1538261
3.74673
E
-
04
2.00152
E
-
04
8
3.581966
1.344972
0.1538260
4.77620
E
-
04
2.42916
E
-
04
10
3.581966
1.344972
0.1538260
5.77512
E
-
04
2.84373
E
-
04
[27]
3.58197
1.34508
0.153723
4.02927
E
-
04
2.06732
E
-
04Table 5
Invariants and error norms for single solitary wave,c
=
0.3, h
=
0.1, k
=
0.01, x
0
=
40, and 0
≤
x
≤
100 (using Runge-Kutta time discretization).
Time
I
1
I
2
I
3
L
2
L
∞
Analytical
3.581967
1.345076
0.1537230
0
3.581967
1.344973
0.1538264
0
0
2
3.581967
1.344973
0.1538263
1.396685
E
-
04
0.920780
E
-
04
4
3.581967
1.344973
0.1538262
2.592545
E
-
04
1.512530
E
-
04
6
3.581967
1.344973
0.1538262
3.647701
E
-
04
1.956752
E
-
04
8
3.581967
1.344973
0.1538262
4.638980
E
-
04
2.368004
E
-
04
10
3.581967
1.344973
0.1538262
5.600205
E
-
04
2.766940
E
-
04
[27]
3.58197
1.34508
0.153723
4.02927
E
-
04
2.06732
E
-
04Figure 3
Solitary wave withc
=
0.3, h
=
0.1, k
=
0.01 and x
0
=
40
,
0
≤
x
≤
100.Figure 4
Error graph forc
=
0.3, h
=
0.1, k
=
0.01 and x
0
=
40
,
0
≤
x
≤
100.The Runge-Kutta Galerkin method algorithm has been run for different space steps with a fixed time stepk
=
0.01. The error and convergence order of space for the L
2- and L
∞-norm are recorded in Tables 6 and 7, respectively, which verify that the L
2- and L
∞-error can achieve the order O
(
h
2
) in space. And the algorithm has also been run for different time steps with a fixed space step h
=
0.1. The error and convergence order of time for the L
2- and L
∞-norm are recorded in Tables 8 and 9, respectively.Table 6
The order of convergence in space forL
2 error norm with k
=
0.01 and c
=
0.3.
Space
Δh
=
0.1
Δh
=
0.2
Δh
=
0.4
Order
Order
(0.1/0.2)
(0.2/0.4)
2
1.396685 × 10−4
5.598770 × 10−4
2.258793 × 10−3
2.0031
2.0124
4
2.592545 × 10−4
1.038054 × 10−3
4.168527 × 10−3
2.0014
2.0057
6
3.647701 × 10−4
1.459612 × 10−3
5.846531 × 10−3
2.0005
2.0020
8
4.638980 × 10−4
1.855582 × 10−3
7.421648 × 10−3
2.0000
1.9999
10
5.600205 × 10−4
2.239531 × 10−3
8.948683 × 10−3
1.9996
1.9985Table 7
The order of convergence in space forL
∞ error norm with k
=
0.01 and c
=
0.3.
Space
Δh
=
0.1
Δh
=
0.2
Δt
=
0.4
Order
Order
(0.1/0.2)
(0.2/0.4)
2
0.920780 × 10−4
3.689931 × 10−4
1.480253 × 10−3
2.0027
2.0042
4
1.512530 × 10−4
6.055462 × 10−4
2.430139 × 10−3
2.0012
2.0047
6
1.956752 × 10−4
7.811765 × 10−4
3.126531 × 10−3
1.9972
2.0008
8
2.368004 × 10−4
9.454389 × 10−4
3.768796 × 10−3
1.9973
1.9950
10
2.766940 × 10−4
1.104758 × 10−3
4.416067 × 10−3
1.9974
1.9990Table 8
The order of convergence in time forL
2 error norm with h
=
0.1 and c
=
0.3.
Time
Δt
=
2.0
Δt
=
1.0
Δt
=
0.5
Order
Order
(2.0/1.0)
(1.0/0.5)
2
0.0511878
0.00418182
4.219760 × 10−4
3.613597
3.308899
4
0.0835505
0.00806197
8.133115 × 10−4
3.309253
3.309253
6
0.1027137
0.01176588
1.180928 × 10−3
3.125948
3.316617
8
0.1141693
0.01561415
1.547231 × 10−3
2.870249
3.335094
10
0.1220153
0.01982141
1.925846 × 10−3
2.621930
3.363496Table 9
The order of convergence in time forL
∞ error norm with h
=
0.1 and c
=
0.3.
Time
Δt
=
2.0
Δt
=
1.0
Δt
=
0.5
Order
Order
(2.0/1.0)
(1.0/0.5)
2
0.0318115
2.528470 × 10−3
2.505105 × 10−4
3.653213
3.335322
4
0.0498808
4.159059 × 10−3
4.456332 × 10−4
3.584155
3.222328
6
0.0594582
5.587773 × 10−3
5.905399 × 10−4
3.411531
3.242167
8
0.0641225
7.106547 × 10−3
7.270440 × 10−4
3.173611
3.289034
10
0.0661177
8.861734 × 10−3
8.671150 × 10−4
2.899376
3.353294


## 5.2.  Collision of Two Solitons

In this section, we study the interaction of two solitary waves with the initial conditions given by a linear sum of two separate solitary waves of various amplitudes:(35)
u
(
x
,
0
)
=
∑
i
=
1
2
‍
A
i
sech
⁡
(
p
i
(
x
-
x
i
)
)
,
where A
i
=
c
i, p
i
=
c
i
/
μ
(
c
i
+
1
), i
=
1
,
2, and x
i and c
i
(
i
=
1,2
) are arbitrary constants. The analytical values for the conservation laws in this case have the following form:
(36)
I
1
=
π
c
1
p
1
+
π
c
2
p
2
,
I
2
=
2
c
1
p
1
+
2
c
2
p
2
+
2
μ
c
1
p
1
3
+
2
μ
p
2
c
2
3
,
I
3
=
4
c
1
2
3
p
1
+
4
c
2
2
3
p
2
-
2
μ
p
1
c
1
3
-
2
μ
p
2
c
2
3
.
In our numerical scheme, we choose h
=
0.1, k
=
0.01, c
1
=
1, c
2
=
0.5, x
1
=
10, and x
2
=
40, with interval [
0,100
]. Then the amplitudes are in ratio 2 : 1, where 2
A
2
=
A
1 and the analytical values for the conservation are I
1
=
8.2905324294, I
2
=
5.224332548, and I
3
=
1.799113742. The changes of the invariants are satisfactorily small, since the changes of the invariants I
1, I
2, and I
3 are 1.3295 × 10−4, 2.2259 × 10−5, and 6.8989 × 10−5, respectively. The results are recorded in Table 10. A graph of the two solitary waves collision, plotted in a three-dimensional view at some discrete times, is shown in Figure 5.Table 10
Invariants and error norms for two solitary waves,c
1
=
1, c
2
=
0.5, x
1
=
10, x
2
=
40, h
=
0.1, k
=
0.01, and 0
≤
x
≤
100.
Time
I
1
I
2
I
3
C-N
R-K
C-N
R-K
C-N
R-K
0
8.288129
8.288129
5.222934
5.222933
1.800512
1.800513
2
8.289522
8.293010
5.222910
5.222951
1.800479
1.800507
4
8.289462
8.295348
5.222887
5.222954
1.800458
1.822507
6
8.289304
8.297822
5.222864
5.222956
1.800433
1.800507
8
8.289233
8.300322
5.222840
5.222958
1.800411
1.800508
10
8.289231
8.302630
5.222817
5.222960
1.800388
1.800509Figure 5
Interaction of two solitons at times fromt
=
2,4
,
6,8
,
10 for c
1
=
1, c
2
=
0.5, x
1
=
10, x
2
=
40, and 0
≤
x
≤
100.


## 5.3.  Collision of Three Solitons

In this section, we study the interaction of three solitary waves with the initial conditions given by a linear sum of three separate solitary waves of various amplitudes:(37)
u
(
x
,
0
)
=
∑
i
=
1
3
‍
A
i
sech
⁡
(
p
i
(
x
-
x
i
)
)
,
where A
i
=
c
i, p
i
=
c
i
/
μ
(
c
i
+
1
), i
=
1,2
,
3, and x
i and c
i
(
i
=
1,2
,
3
) are arbitrary constants. The analytical values for the conservation laws in this case have the following form:
(38)
I
1
=
π
c
1
p
1
+
π
c
2
p
2
+
π
c
3
p
3
,
I
2
=
2
c
1
p
1
+
2
c
2
p
2
+
2
c
3
p
3
+
2
μ
c
1
p
1
3
+
2
μ
p
2
c
2
3
+
2
μ
p
3
c
3
3
,
I
3
=
4
c
1
2
3
p
1
+
4
c
2
2
3
p
2
+
4
c
3
2
3
p
3
-
2
μ
p
1
c
1
3
-
2
μ
p
2
c
2
3
-
2
μ
p
3
c
3
3
.
In this case, we choose paraments as h
=
0.1, k
=
0.05, c
1
=
1, c
2
=
0.5, c
3
=
0.25, x
1
=
10, x
2
=
20, and x
3
=
40, with interval [
0,100
]. Then the amplitudes are in ratio 4 : 2 : 1, where A
1
=
2
A
2
=
4
A
3 and the analytical values for the conservation are I
1
=
11.802939794, I
2
=
6.416902131, and I
3
=
1.910917141. The changes of the invariants from the initial variants approach zero throughout and agree with the analytical values for the three invariants as presented in Table 11, which indicates that our scheme is satisfactorily conservative. The geometry of the initial state and the profiles at time t
=
10, t
=
20, and t
=
35 is shown graphically in Figure 6, respectively.Table 11
Invariants and error norms for three solitary waves,c
1
=
1, c
2
=
0.5, c
3
=
0.25, x
1
=
10, x
2
=
20, x
3
=
40, h
=
0.1, k
=
0.01, and 0
≤
x
≤
100.
Time
I
1
I
2
I
3
C-N
R-K
C-N
R-K
C-N
R-K
0
11.800513
11.800513
6.490318
6.490318
1.997082
1.997082
4
11.799364
11.807800
6.485747
6.490319
1.992012
1.996623
8
11.797850
11.812825
6.483120
6.490298
1.988853
1.996076
12
11.797570
11.817192
6.482677
6.490279
1.988441
1.996360
16
11.791837
11.820521
6.474939
6.490305
1.982852
1.996875
20
11.788577
11.823320
6.469884
6.490341
1.977926
1.996758Interaction of three solitary waves.
(a)
Time = 0
(b)
Time = 10
(c)
Time = 20
(d)
Time = 35


## 5.4.  The Maxwellian Initial Condition

Finally, we consider the development of the Maxwellian initial condition(39)
u
(
x
,
0
)
=
e
-
(
x
-
40
)
2
into a train of solitary waves. We choose different values μ
=
1.0, μ
=
0.5, and μ
=
0.1 in our numerical scheme. The comparison of the three variants I
1, I
2, and I
3 with the earlier result in [29] is presented in Table 12, here Crank-Nicolson method is used for time discretization. For μ
=
1.0, the changes of the variants I
1, I
2, and I
3 with respect to the initial values are 3.9492 × 10−6, 6.9575 × 10−4, and 3.1509 × 10−3, respectively, whereas they are 2.8210 × 10−6, 2.60000 × 10−4, and 7.0771 × 10−4 in [29]. For μ
=
0.5, the variants are changed by 7.8982 × 10−6, 6.8214 × 10−5, and 6.1410 × 10−4 in this case. The development of the Maxwellian initial condition is shown in Figure 7 with different parameter values, respectively. The smaller μ there is, the more the number of solitary waves will form. The simulations are done up to time = 10 in this case. From the previous work [25], the total number of solitary waves and the values μ have the approximate relation
(40)
N
≅
[
1
μ
5
]
.Table 12
Values ofI
1, I
2, and I
3 for Maxwellian initial condition when h
=
0.1 and k
=
0.01 and the space interval [
0,100
].
μ
Time
I
1
I
2
I
3
I
1 [29]
I
2 [29]
I
3 [29]
1
2
1.772528
2.499467
−0.3633570
1.772449
2.506352
−0.3668149
4
1.772544
2.499870
−0.3643380
1.772446
2.506235
−0.3666974
6
1.772532
2.500822
−0.3657774
1.772447
2.506171
−0.3666326
8
1.772524
2.500977
−0.3661331
1.772446
2.506123
−0.3665860
10
1.772521
2.501206
−0.3665079
1.772444
2.2060926
−0.3665553
0.5
2
1.772534
1.876446
0.2597450
1.772451
1.879888
0.2596494
4
1.772534
1.876534
0.2587968
1.772446
1.879855
0.2596828
6
1.772528
1.876600
0.2583802
1.772449
1.879841
0.2596973
8
1.772523
1.876581
0.2582386
1.772450
1.879834
0.2597050
10
1.772520
1.876574
0.2581499
1.772449
1.879828
0.2597092
0.1
2
1.770991
1.373320
0.7687222
1.772452
1.378607
0.7608777
4
1.769863
1.371236
0.7658462
1.772451
1.378577
0.7608364
6
1.768742
1.369139
0.7630269
1.772451
1.378546
0.7607937
8
1.767630
1.367058
0.7602389
1.772451
1.378515
0.7607529
10
1.766526
1.364993
0.7574795
1.772453
1.378483
0.7607117Maxwellian initial condition with different valuesμ.
(a)
μ
=
1.0
(b)
μ
=
0.5
(c)
μ
=
0.1


## 6. Conclusion

In this paper, we have developed a Galerkin linear finite element method to investigate the propagation of solitons and their interactions governed by the nonlinear MRLW equation. An extrapolation technique has been used to improve the accuracy in time for this numerical method. A linear stability analysis shows that this method is marginally stable. The high efficiency and accuracy of our method is tested by the numerical experiments: propagation of single soliton, collision of two and three solitons, and development of the Maxwellian initial condition into solitary waves. Moreover, this method numerically satisfies the conservation laws of mass, momentum, and energy.



---

*Source: 438289-2014-06-19.xml*