# Solitary Wave Formation from a Generalized Rosenau Equation

Authors: J. I. Ramos; C. M. García-López
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A generalized viscous Rosenau equation containing linear and nonlinear advective terms and mixed third- and fifth-order derivatives is studied numerically by means of an implicit second-order accurate method in time that treats the first-, second-, and fourth-order spatial derivatives as unknown and discretizes them by means of three-point, fourth-order accurate, compact finite differences. It is shown that the effect of the viscosity is to decrease the amplitude, curve the wave trajectory, and increase the number and width of the waves that emerge from an initial Gaussian condition, whereas the linear convective term pushes the wave front towards the downstream boundary. It is also shown that the effect of the nonlinear convective term is to increase the steepness of the leading wave front and the number of sawtooth waves that are generated behind it, while that of the first dispersive term is to increase the number of waves that break up from the initial condition as the coefficient that characterizes this term is decreased. It is also shown that, for reasons of stability, the second dispersion coefficient must be much smaller than the first one and its effects on wave propagation are relatively small.


---

## Body

## 1. Introduction

In his 1986 and 1988 papers, Rosenau [1, 2] developed a formalism to treat the dynamics of discrete dense systems that can deal with wave-wave and wave-wall interactions that cannot be treated with the Korteweg-de Vries (KdV) equation. Such a formalism leads to a quasi-continuum which is endowed with the leading effects due to discreteness.Rosenau’s original equation may be written as(1)
u
t
±
u
+
u
2
x
+
u
x
x
x
x
t
=
0
, where the subscripts t and x denote partial differentiation with respect to the time and the spatial coordinate, respectively. In the absence of the fifth-order derivative term, (1) reduces to a first-order nonlinear wave equation which has analytical solutions and may have shock wave solutions if u
x
(
0
,
x
)
<
0.The Rosenau equation has been the subject of several analytical and numerical studies. Barreto et al. [3] prove the existence of solutions of (1) with the plus sign in the advection-like term in moving domains by making use of the Galerkin method, multiplier techniques, and energy estimates. They also proved analogous results for the Benjamin-Bona-Mahony or regularized long-wave (RLW) equation with a linear advective term; that is,(2)
u
t
+
u
+
u
2
x
-
u
x
x
t
=
0
, which models long waves in a nonlinear dispersive medium.A one-dimensional generalization of the Rosenau equation may be written as(3)
u
t
+
ϕ
x
u
=
f
t
,
x
, where f
(
t
,
x
) denotes a forcing function and(4)
ϕ
u
=
α
u
+
ϵ
u
p
-
μ
u
x
-
δ
u
x
t
+
θ
u
m
x
x
-
β
u
u
x
x
-
ν
u
x
2
-
ψ
u
x
x
x
t
-
ω
u
x
x
x
x
.Equation (4) with ϕ
(
u
)
=
±
(
u
+
u
2
)
+
u
x
x
x
t and f
(
t
,
x
)
=
0 corresponds to the original Rosenau equation and includes the linear (first-order) wave equation, the inviscid and viscous Burgers equations, the RLW, modified RLW, and generalized RLW equations, the Korteweg-de Vries (KdV) equation, the K
(
m
,
n
) Rosenau-Hyman equation, the Camassa-Holm equation, the Olver-Rosenau equation, the Kawahara equation, the Cooper-Shepherd-Sodano equation, and combinations thereof; these equations are presented in Table 1.Table 1
Some typical equations included in (3) and (4).
Name
Equation
Linear wave equation
u
t
+
c
u
x
=
0
Inviscid Burgers equation
u
t
+
u
u
x
=
0
Inviscid generalized Burgers equation
u
t
+
u
p
x
=
0
Diffusion equation
u
t
=
α
u
x
x
Equal-width (EW) equation
u
t
+
u
2
x
=
δ
u
x
x
t
Viscous EW equation
u
t
+
u
2
x
=
δ
u
x
x
t
+
μ
u
x
x
RLW equation
u
t
+
u
x
+
ϵ
u
2
x
=
δ
u
x
x
t
Viscous RLW equation
u
t
+
u
x
+
ϵ
u
2
x
=
δ
u
x
x
t
+
μ
u
x
x
Generalized RLW (GRLW) equation
u
t
+
u
x
+
ϵ
u
p
x
=
δ
u
x
x
t, p
≠
2
Viscous GRLW equation
u
t
+
u
x
+
ϵ
u
p
x
=
δ
u
x
x
t
+
μ
u
x
x, p
≠
2
Rosenau's original equation
u
t
±
u
+
u
2
x
+
u
x
x
x
x
t
=
0
Rosenau-Hyman'sK
(
m
,
n
) equation
u
t
±
u
m
x
+
u
n
x
x
x
=
0
Camassa-Holm equation
u
t
-
u
x
x
t
+
3
u
u
x
=
2
u
x
u
x
x
+
u
u
x
x
x
Cooper-Shepherd-Sodano equation
u
t
+
u
l
-
2
u
x
-
p
u
p
-
1
u
x
2
x
+
2
α
u
p
u
x
x
x
=
0
Olver-Rosenau equation
u
±
u
x
x
t
=
a
u
x
+
1
2
u
2
±
u
2
x
u
±
u
x
x
x
Korteweg-de Vries equation
u
+
u
u
x
+
δ
u
x
x
x
=
0
Kawahara (K) equation
u
t
+
α
u
u
x
+
β
u
x
x
x
-
γ
u
x
x
x
x
x
=
0
Modified Kawahara equation
u
t
+
α
u
2
u
x
+
β
u
x
x
x
-
γ
u
x
x
x
x
x
=
0
Generalized Kawahara equation
u
t
+
α
u
p
u
x
+
β
u
x
x
x
-
γ
u
x
x
x
x
x
=
0, p
≠
1,2Park [4] studied the global existence and uniqueness of solutions of the following multidimensional generalized Rosenau equation:(5)
u
t
+
c
Δ
u
t
=
∇
·
f
u
,where Δ
≡
∇
2 is the Laplacian or harmonic operator, while the same author [5] provided pointwise decay estimates for (3) with ϕ
(
u
)
=
u
p
+
1
/
p
+
1
+
u
x
x
x
x
t, that is, a Rosenau-Burgers equation. Wang and Xu [6] studied the existence and global solution of the second-order in time Rosenau equation u
t
t
+
u
x
x
x
x
t
t
-
γ
u
x
x
=
(
f
(
u
)
)
x
x and proved the existence of finite-time blowup for f
(
u
)
=
-
β
|
u
|
p
u for β
>
0 and p
>
0, by means of a potential well technique, while H. Wang and S. Wang [7] studied the long-time behavior of small solutions of the Cauchy problem for a (second-order in time) Rosenau equation, obtained its global small solution, and analyzed the decay and scattering of such a solution. A similar study for the (first-order in time) Rosenau-Burgers (R-B) equation was performed by Liu and Mei [8] who proved that the solution of a nonlinear parabolic equation is a better asymptotic profile of the R-B equation. H. Wang and S. Wang [9] proved the global existence of the solution to the Cauchy problem for the nth dimensional Rosenau equation with a damping term when the initial data are small, while Kim and Lee [10] analyzed the convergence of a semidiscretization of u
t
+
(
a
(
t
,
x
)
u
x
x
t
)
x
x
=
ψ
x
(
u
).Esfahani [11, 12] obtained solitary wave solutions to the generalized Rosenau-KdV and Rosenau-RLW equations, respectively, by means of the sech and trigonometric function methods, respectively, while Razborova et al. [13] studied the solitons, shock waves, and conservation properties of Rosenau-KdV-RLW equations with power nonlinearities by means of a semi-inverse variational method. Razborova et al. [14] used a perturbation method to study shallow water waves governed by the same equation. Choo et al. [15] obtained a posteriori error estimates of the Rosenau equation by means of a discontinuous Galerkin method and analyzed the stability of the dual problem.In Table2, some numerical methods that have been used to study one-dimensional Rosenau equations are summarized. These methods include finite difference and finite element techniques for the space discretization, linear and iterative implicit time discretizations, and time integration with explicit Runge-Kutta procedures of third- and fifth-order accuracy.Table 2
Examples of generalized Rosenau equationsu
t
+
ψ
(
u
,
u
x
,
u
x
x
)
+
θ
u
x
x
x
x
t = 0, their acronyms, and numerical methods for their solution.
E.A.
ψ
(
u
)
Numerical method
R
3
2
u
+
12
u
5
-
u
3
dGM (θ
=
1
/
2) [15]
R
u
+
1
2
u
2
CN, iterative SOR FDM [22]
R
u
+
1
2
u
2
Conservative IL3TL FDM [23]
R
u
+
u
2
EB Galerkin FEM [24]
R
u
+
u
p
Conservative IL3TL FDM [25]
R
u
+
u
p
+
1
p
+
1
Conservative IL3TL FDM [26]
R-B
u
+
1
2
u
2
-
u
x
Conservative IL3TL FDM [27]
R-B
u
+
1
2
u
2
-
u
x
Conservative IL3TL FDM [28]
R-B
u
+
1
2
u
2
-
u
x
CN FDM [29]
R-B
u
+
1
2
u
2
-
u
x
RK3 QBSPCM [30]
R-B
β
u
+
1
2
u
2
-
α
u
x
Conservative IL3TL FDM [31]
R-B
u
+
u
p
+
1
p
+
1
-
u
x
Implicit IL2TL FDM [32]
R-B
u
+
u
p
+
1
p
+
1
-
u
x
Conservative IL3TL FDM [33]
R-KdV
u
+
1
2
u
2
+
u
x
x
Conservative IL3TL FDM [34]
R-KdV
u
+
1
2
u
2
+
u
x
x
RK3 QBSPCM [35]
R-KdV
u
+
α
u
2
+
β
u
x
x
Conservative IL3TL FDM [36]
R-KdV
u
+
u
p
+
u
x
x
Conservative CN FDM [37]
R-KdV
u
+
u
p
+
u
x
x
EB FEM [38]
R-KdV
u
+
u
p
+
u
x
x
CN FEM [38]
R-KdV
u
+
u
p
+
u
x
x
IL3TL FEM [38]
R-RLW
u
+
u
p
-
u
x
t
Conservative IL2TL FDM [39]
R-RLW
u
+
u
p
-
u
x
t
Conservative IL3TL FDM [40]
R-RLW
u
+
1
2
u
2
-
u
x
t
Galerkin FEM with CN [41]
R-RLW
u
+
1
2
u
2
-
u
x
t
Conservative CN FDM [42]
R-RLW
u
+
1
2
u
2
-
u
x
t
IL3TL FDM [43]
R-RLW
u
+
1
2
u
2
-
u
x
t
RK54 QBSPCM [44]
R-RLW
u
+
α
u
2
-
γ
u
x
t
Conservative IL3TL FDM [36]
R-KdV-RLW
u
+
α
u
2
+
β
u
x
x
-
γ
u
x
t
Conservative IL3TL FDM [36]
Equation acronym (E.A.): R  =  Rosenau, RLW  =  regularized long-wave, B  =  Burgers, and KdV =  Korteweg-de Vries. Numerical method: dGM  =  discontinuous Galerkin method, FDM  =  finite difference method, FEM  =  finite element method, QBSPCM  =  quintic B-splines collocation method, SOR  =  successive overrelaxation, CN  =  Crank-Nicolson, RK3 =  third-order accurate Runge-Kutta method, IL2TL  =  implicit linear two-time level, IL3TL  =  implicit linear three-time level, and EB  =  Euler's backward time discretization. Unless stated otherwiseθ = 1.Most of the methods presented in Table2 are second-order accurate in both space and time, while the one presented in this paper is fourth-order accurate in space and second-order accurate in time. Moreover, the accuracy of most of the methods presented in Table 2 has been assessed by what is referred to as the method of manufactured solutions; that is, the coefficients of (3) and (4) are found so that the solution of (3) is of a specified travelling-wave, for example, hyperbolic secant, type and the initial condition corresponds to the exact solution at t
=
0. For such an initial condition, u
(
t
,
x
) is very smooth and does not exhibit steep gradients; as a consequence, a relative small number of grid points are required to obtain very accurate numerical solutions. However, if the initial profile does not correspond to an exact solution of (3) and u
x
(
0
,
x
) is negative, the leading part of the initial condition steepens due to the nonlinear advective term and, in the absence of dispersion and diffusion, would result in the formation of a shock wave [16, 17]. For the same type of initial conditions but with dissipation and no dispersion, the initial steepening is eventually balanced by viscous dissipation and a travelling wave which may be referred to as Taylor’s wave is formed [17–19]; such a wave has a finite thickness. On the other hand, for the same initial conditions as the ones discussed above but with dispersion and no diffusion, the initial steepening of the leading part of the initial profile caused by the nonlinear advective term is eventually somewhat balanced by dispersion and a dispersive shock wave or conservative undular bore may form [17, 20, 21].In this paper, a generalized viscous Rosenau equation that includes linear and nonlinear advective terms is studied numerically by means of a second-order accurate, linearized Crank-Nicolson method and three-point, fourth-order accurate, compact operator discretizations for the first-, second-, and fourth-order spatial derivatives. The breakup of the initial condition is studied as a function of the linear and nonlinear convective terms and the viscous and two dispersive terms that appear in the equation.The paper has been arranged as follows. In the next section, the generalized Rosenau equation considered in the study reported here is presented and the linear stability of its linear counterpart is analyzed. This is followed by a section where the numerical method employed to solve the equation and its linear stability are considered. The fourth section presents an exhaustive numerical study of the effects of the linear and nonlinear advective, viscous, and dispersive terms and initial conditions on wave breakup, generation, and propagation. Finally, a short concluding section summarizes the most important findings reported in the paper.


## 2. Governing Equation

In this paper, the following generalized Rosenau equation is considered:(6)
u
t
+
α
u
x
+
ϵ
u
4
x
=
μ
u
x
x
+
δ
u
x
x
t
+
γ
u
x
x
x
x
t
,
t
>
0
,
-
∞
<
x
<
∞
,where α, ϵ, μ, δ, and γ are constants that are associated with the linear and nonlinear advective, dissipative, and third- and fifth-order dispersive terms, respectively. Equation (6) is a particular case of (3) and (4), includes both linear and nonlinear terms, and will be solved subject to the following initial(7)
u
0
,
x
=
f
x and boundary conditions(8)
u
t
,
±
∞
=
u
x
x
t
,
±
∞
=
0
.Before proceeding with the discretization of (6), it is convenient to analyze its linear counterpart.


### 2.1. Linear Analysis

The linear counterpart of (6) may be written as(9)
u
t
+
α
u
x
=
μ
u
x
x
+
δ
u
x
x
t
+
γ
u
x
x
x
x
t
,
t
>
0
,
-
∞
<
x
<
∞
, which, upon introducing u
t
,
x
=
A
exp
⁡
(
I
(
k
x
-
ω
t
)
), where I
2
=
-
1, ω and k denote the (angular) frequency and wavenumber, respectively, and A is the amplitude, becomes(10)
I
ω
1
+
δ
k
2
-
γ
k
4
=
μ
k
2
+
I
k
α
, which, for real k and complex ω, that is, ω
=
ω
R
+
I
ω
I, where ω
R and ω
I denote the real and imaginary parts, respectively, of ω, becomes, upon separating the real and imaginary parts,(11)
ω
I
1
+
δ
k
2
-
γ
k
4
=
-
μ
k
2
,
(12)
ω
R
1
+
δ
k
2
-
γ
k
4
=
α
k
.The above relations indicate that instabilities arise wheneverω
I
>
0 and this condition requires that (
1
+
δ
k
2
-
γ
k
4
)
<
0 for μ
>
0 which is not fulfilled if δ
≥
0 and γ
≤
0. Furthermore, (11) and (12) indicate that neither ω
R nor ω
I is defined if (
1
+
δ
k
2
-
γ
k
4
)
=
0 provided that neither α nor μ is zero, respectively.For small wavenumber, that is,k
≪
1, (11) and (12) indicate that ω
R
=
α
k
+
O
(
k
3
) and ω
I
=
-
μ
k
2
+
O
(
k
4
), thus showing that long waves are stable (the viscosity coefficient μ
≥
0). On the other hand, for very large wavenumbers, that is, very small wavelengths, these equations show that ω
R
=
-
α
/
γ
k
-
3
+
O
(
k
-
5
) and ω
I
=
μ
/
γ
k
-
2
+
O
(
k
-
4
) which indicate that short waves are unstable if γ
>
0. From this analysis, it may be concluded that no linear instabilities occur if δ and γ are positive and negative, respectively, and, for these values, S
(
k
;
δ
,
γ
)
≡
1
+
δ
k
2
-
γ
k
4 is a positive monotonously increasing function of the wavenumber and achieves its minimum value of one for k
=
0. Moreover, these conditions on δ and γ also imply that the phase velocity, that is, v
p
h
≡
ω
R
/
k
=
α
/
S
(
k
;
δ
,
γ
), is a monotonously decreasing function of k and is equal to α, the linear wave speed in (6), for k
=
0; on the other hand, ω
I
=
0 for k
=
0; that is, there is no damping for k
=
0.Forγ
<
0, the minimum value of ω
I occurs at k
=
(
-
γ
)
-
1
/
4 and is equal to -
μ
/
δ
+
2
-
γ, while the group velocity, that is, v
g
≡
d
ω
R
/
d
k
=
α
k
1
-
δ
k
2
+
3
γ
k
4
/
S
2, is not defined whenever S
=
0 and v
g
=
α
+
O
(
k
2
) and 3
α
γ
k
-
4
+
O
(
k
-
6
) for k
≪
1 and k
≫
1, respectively. The group velocity is positive for k
∈
[
0
,
k
c
), nil for k
=
k
c, and negative for k
∈
(
k
c
,
∞
), where k
c
2
=
δ
-
δ
2
-
12
γ
/
12
γ and γ
<
0.


### 2.2. Invariants

Integration of (6) from -
∞ to ∞ yields(13)
I
1
t
=
∫
-
∞
∞
u
t
,
x
d
x
=
C
1
, where C
1 is a constant equal to ∫
-
∞
∞
‍
u
(
0
,
x
)
d
x
=
∫
-
∞
∞
‍
f
(
x
)
d
x.Multiplication of (6) by u
(
t
,
x
) and integration of the resulting equation from -
∞ to ∞ yield(14)
d
I
2
d
t
t
=
d
d
t
∫
-
∞
∞
E
t
,
x
d
x
=
-
μ
∫
-
∞
∞
u
x
2
d
x
, where the total energy density is(15)
E
t
,
x
=
1
2
u
2
+
δ
u
x
2
-
γ
u
x
x
2 and contains three semipositive terms if δ
>
0 and γ
<
0 and, for these values, I
2
(
t
) decreases with time if μ is not nil. Note that the third term on the right-hand side of (15) is positive if γ
<
0.


## 2.1. Linear Analysis

The linear counterpart of (6) may be written as(9)
u
t
+
α
u
x
=
μ
u
x
x
+
δ
u
x
x
t
+
γ
u
x
x
x
x
t
,
t
>
0
,
-
∞
<
x
<
∞
, which, upon introducing u
t
,
x
=
A
exp
⁡
(
I
(
k
x
-
ω
t
)
), where I
2
=
-
1, ω and k denote the (angular) frequency and wavenumber, respectively, and A is the amplitude, becomes(10)
I
ω
1
+
δ
k
2
-
γ
k
4
=
μ
k
2
+
I
k
α
, which, for real k and complex ω, that is, ω
=
ω
R
+
I
ω
I, where ω
R and ω
I denote the real and imaginary parts, respectively, of ω, becomes, upon separating the real and imaginary parts,(11)
ω
I
1
+
δ
k
2
-
γ
k
4
=
-
μ
k
2
,
(12)
ω
R
1
+
δ
k
2
-
γ
k
4
=
α
k
.The above relations indicate that instabilities arise wheneverω
I
>
0 and this condition requires that (
1
+
δ
k
2
-
γ
k
4
)
<
0 for μ
>
0 which is not fulfilled if δ
≥
0 and γ
≤
0. Furthermore, (11) and (12) indicate that neither ω
R nor ω
I is defined if (
1
+
δ
k
2
-
γ
k
4
)
=
0 provided that neither α nor μ is zero, respectively.For small wavenumber, that is,k
≪
1, (11) and (12) indicate that ω
R
=
α
k
+
O
(
k
3
) and ω
I
=
-
μ
k
2
+
O
(
k
4
), thus showing that long waves are stable (the viscosity coefficient μ
≥
0). On the other hand, for very large wavenumbers, that is, very small wavelengths, these equations show that ω
R
=
-
α
/
γ
k
-
3
+
O
(
k
-
5
) and ω
I
=
μ
/
γ
k
-
2
+
O
(
k
-
4
) which indicate that short waves are unstable if γ
>
0. From this analysis, it may be concluded that no linear instabilities occur if δ and γ are positive and negative, respectively, and, for these values, S
(
k
;
δ
,
γ
)
≡
1
+
δ
k
2
-
γ
k
4 is a positive monotonously increasing function of the wavenumber and achieves its minimum value of one for k
=
0. Moreover, these conditions on δ and γ also imply that the phase velocity, that is, v
p
h
≡
ω
R
/
k
=
α
/
S
(
k
;
δ
,
γ
), is a monotonously decreasing function of k and is equal to α, the linear wave speed in (6), for k
=
0; on the other hand, ω
I
=
0 for k
=
0; that is, there is no damping for k
=
0.Forγ
<
0, the minimum value of ω
I occurs at k
=
(
-
γ
)
-
1
/
4 and is equal to -
μ
/
δ
+
2
-
γ, while the group velocity, that is, v
g
≡
d
ω
R
/
d
k
=
α
k
1
-
δ
k
2
+
3
γ
k
4
/
S
2, is not defined whenever S
=
0 and v
g
=
α
+
O
(
k
2
) and 3
α
γ
k
-
4
+
O
(
k
-
6
) for k
≪
1 and k
≫
1, respectively. The group velocity is positive for k
∈
[
0
,
k
c
), nil for k
=
k
c, and negative for k
∈
(
k
c
,
∞
), where k
c
2
=
δ
-
δ
2
-
12
γ
/
12
γ and γ
<
0.


## 2.2. Invariants

Integration of (6) from -
∞ to ∞ yields(13)
I
1
t
=
∫
-
∞
∞
u
t
,
x
d
x
=
C
1
, where C
1 is a constant equal to ∫
-
∞
∞
‍
u
(
0
,
x
)
d
x
=
∫
-
∞
∞
‍
f
(
x
)
d
x.Multiplication of (6) by u
(
t
,
x
) and integration of the resulting equation from -
∞ to ∞ yield(14)
d
I
2
d
t
t
=
d
d
t
∫
-
∞
∞
E
t
,
x
d
x
=
-
μ
∫
-
∞
∞
u
x
2
d
x
, where the total energy density is(15)
E
t
,
x
=
1
2
u
2
+
δ
u
x
2
-
γ
u
x
x
2 and contains three semipositive terms if δ
>
0 and γ
<
0 and, for these values, I
2
(
t
) decreases with time if μ is not nil. Note that the third term on the right-hand side of (15) is positive if γ
<
0.


## 3. Finite Difference Discretization

Equation (6) may be written as the following system of equations:(16)
u
t
+
α
F
+
4
ϵ
u
3
F
=
μ
G
+
δ
G
t
+
γ
H
t
,
(17)
F
=
u
x
,
(18)
G
=
u
x
x
,
(19)
H
=
G
x
x
.Using a second-order accurate time discretization and linearizing the nonlinear terms with respect to the previous time level, (16) becomes(20)
Δ
u
=
Δ
t
2
μ
2
G
n
+
Δ
G
+
δ
Δ
G
+
γ
Δ
H
-
α
2
F
n
+
Δ
F
-
4
ϵ
2
T
n
+
P
n
Δ
u
+
Q
n
Δ
F
, where(21)
Δ
u
=
u
n
+
1
-
u
n
,
T
=
u
3
F
,
P
=
∂
T
∂
u
=
3
u
2
F
,
Q
=
∂
T
∂
F
=
u
3
,
Δ
t is the time step, and the superscript n denotes the nth time level.Equation (20) together with (17)–(19) represents a linear fourth-order ordinary differential equation for Δ
u at each time level and is subject to the boundary conditions specified in (8); note that, for example, Δ
G
=
Δ
u
x
x. However, (20) may not be solved analytically due to the dependence of G
n, F
n, H
n, T
n, P
n, and Q
n on x, but it becomes a discrete equation upon approximating these variables and F, G, and H by finite differences.The infinite domain(
-
∞
,
+
∞
) was truncated to a finite (
a
,
b
) one, where |
a
|
<
∞ and |
b
|
<
∞ were chosen so that the locations of these boundaries do not affect the numerical results, and the finite domain was discretized into N equally spaced nonoverlapping intervals, so that x
i
=
a
+
(
i
-
1
)
h,  i
=
1,2
,
…
,
N
+
1, x
1
=
a, x
N
+
1
=
b, and h
=
b
-
a
/
N is the grid spacing.In this paper, the following three-point, compact, fourth-order accurate, finite difference approximations for (17)–(19) have been used:(22)
1
6
F
i
+
1
+
4
F
i
+
F
i
-
1
=
1
2
h
u
i
+
1
-
u
i
-
1
,
(23)
1
12
G
i
+
1
+
10
G
i
+
G
i
-
1
=
1
h
2
u
i
+
1
-
2
u
i
+
u
i
-
1
,
(24)
1
12
H
i
+
1
+
10
H
i
+
H
i
-
1
=
1
h
2
G
i
+
1
-
2
G
i
+
G
i
-
1
.From (8), u
1
=
u
N
+
1
=
0 and G
1
=
G
N
+
1
=
0. Furthermore, by assuming that (6) is valid at the boundaries and making use of (8), it is easy to obtain that(25)
α
F
t
,
a
=
γ
H
t
t
,
a
,
α
F
t
,
b
=
γ
H
t
t
,
b
, where the boundary conditions have been applied at the boundaries of the truncated domain. This approximation is valid provided that a and b are sufficiently far away.Equation (25) clearly indicates that if H
(
t
,
a
)
=
H
(
t
,
b
)
=
0, that is, H
1
=
H
N
+
1
=
0, then F
(
t
,
a
)
=
F
(
t
,
b
)
=
0, that is, F
1
=
F
N
+
1
=
0; these conditions have been used in (22) and (24) to determine F
i and H
i, with i
=
2,3
,
…
,
N, imply that u
x
(
t
,
a
)
=
u
x
(
t
,
b
)
=
u
x
x
x
x
(
t
,
a
)
=
u
x
x
x
x
(
t
,
b
)
=
0 and, together with (8), result in eight boundary conditions for a fourth-order partial differential equation.It has been found that provided that the solution of (6) is sufficiently smooth near the boundaries of the truncated domain, this overspecification of the boundary conditions does not affect the numerical solution of that equation. A similar behavior was also observed in numerical experiments where F
1 and F
N
+
1 were determined by means of fourth-order accurate, one-sided finite difference approximations that make use of five grid points, one located at the boundary and four interior ones, expanding the values of, say, u
i with i
=
2,3
,
4,5, in Taylor’s series expansions about x
1
=
a and eliminating u
x
x, u
x
x
x, and u
x
x
x
x at x
1 from these expansions; an analogous procedure was also used to determine a fourth-order accurate approximation to H at the boundaries, except that, in this case, the Taylor series expansions of H
i with i
=
2,3
,
4,5
,
6,7 were performed with respect to u
x
x, say, at x
1 and u
x, u
x
x
x, u
x
x
x
x
x, u
x
x
x
x
x
x, and u
x
x
x
x
x
x
x at x
1 were eliminated from these expansions. These one-sided approximations, however, result in nonsymmetric, nontridiagonal matrices for F
i and H
i (cf. (22) and (24)) and were found to be less robust than those obtained with (22)–(24) and F
1
=
F
N
+
1
=
H
1
=
H
N
+
1
=
0, especially when the solution of (6) approached the boundaries.The finite difference method presented in this section is linear at each time step and second- and fourth-order accurate in time and space, respectively, and its stencil consists of only three grid points. Moreover, (20) at x
i with (22)–(24) may be written as a block tridiagonal matrix for (
u
i
,
F
i
,
G
i
,
H
i
)
T with i
=
2,3
,
…
,
N, where the superscript T denotes transpose that may be solved by means of the block tridiagonal matrix method. Alternatively, one may easily determine F
=
A
u, G
=
B
u, and H
=
C
G
=
C
B
u from (22)–(24), respectively, where the matrices A, B
, and C can be easily determined from those equations and, for example, G
=
(
G
2
,
G
3
,
…
,
G
N
)
T; the values of F, G, and H thus obtained can then be substituted into (20) applied at all the interior grid points to obtain a tridiagonal matrix for u.Although the three-point compact operator method presented here is formally fourth-order accurate in space, it must be noted that this scheme as well as many higher-order ones may result in oscillations in regions where steep gradients exist or the solution is not properly resolved. These oscillations are a consequence of the fact that compact and higher-order methods result in finite difference equations which have more solutions than those of the continuous problem. The spurious solutions of these methods may not cause numerical problems if their magnitude is smaller than those associated with the main roots of the characteristic polynomial of the finite difference equation.


### 3.1. Linear Stability of the Finite Difference Discretization

Consider the linear counterpart of (6), that is, (9), which upon time discretization may be written as (cf. (20))(26)
Δ
u
i
=
Δ
t
2
μ
2
G
i
n
+
Δ
G
i
+
δ
Δ
G
i
+
γ
Δ
H
i
-
α
2
F
i
n
+
Δ
F
i
, and assume that(27)
u
i
n
=
U
n
exp
⁡
I
i
k
h
,
F
i
n
=
F
n
exp
⁡
I
i
k
h
,
(28)
G
i
n
=
G
n
exp
⁡
I
i
k
h
,
H
i
n
=
H
n
exp
⁡
I
i
k
h
.Substitution of the above expressions in (22)–(24) yields(29)
F
n
=
I
P
U
n
,
P
=
3
h
sin
⁡
2
ψ
3
-
2
s
i
n
2
ψ
,
G
n
=
-
Q
U
n
,
Q
=
12
h
2
s
i
n
2
ψ
3
-
s
i
n
2
ψ
,
H
n
=
-
Q
G
n
, where ψ
=
1
/
2
k
h.Substitution of (27) and (29) into (26) yields(30)
U
n
+
1
U
n
=
M
1
-
I
N
M
2
+
I
N
, where(31)
M
1
=
M
-
1
2
μ
Δ
t
Q
,
M
2
=
M
+
1
2
μ
Δ
t
Q
,
M
=
1
+
δ
Q
-
γ
Q
2
,
N
=
1
2
α
Δ
t
P
, and the condition of linear stability, that is, |
U
n
+
1
|
≤
|
U
n
|, is satisfied provided that μ
Δ
t
Q
≥
0 which is indeed the case because Q
>
0 and μ
≥
0. The method is, however, dispersive and the phase of U
n
+
1
/
U
n is given by θ
n
u
m
=
arctan
⁡
-
2
M
N
/
M
2
-
μ
Δ
t
/
2
Q
2
-
N
2 which for k
h
≪
1, μ
k
2
Δ
t
≪
1, and |
α
|
k
Δ
t
≪
1 is tan
⁡
θ
≈
-
α
k
Δ
t.


## 3.1. Linear Stability of the Finite Difference Discretization

Consider the linear counterpart of (6), that is, (9), which upon time discretization may be written as (cf. (20))(26)
Δ
u
i
=
Δ
t
2
μ
2
G
i
n
+
Δ
G
i
+
δ
Δ
G
i
+
γ
Δ
H
i
-
α
2
F
i
n
+
Δ
F
i
, and assume that(27)
u
i
n
=
U
n
exp
⁡
I
i
k
h
,
F
i
n
=
F
n
exp
⁡
I
i
k
h
,
(28)
G
i
n
=
G
n
exp
⁡
I
i
k
h
,
H
i
n
=
H
n
exp
⁡
I
i
k
h
.Substitution of the above expressions in (22)–(24) yields(29)
F
n
=
I
P
U
n
,
P
=
3
h
sin
⁡
2
ψ
3
-
2
s
i
n
2
ψ
,
G
n
=
-
Q
U
n
,
Q
=
12
h
2
s
i
n
2
ψ
3
-
s
i
n
2
ψ
,
H
n
=
-
Q
G
n
, where ψ
=
1
/
2
k
h.Substitution of (27) and (29) into (26) yields(30)
U
n
+
1
U
n
=
M
1
-
I
N
M
2
+
I
N
, where(31)
M
1
=
M
-
1
2
μ
Δ
t
Q
,
M
2
=
M
+
1
2
μ
Δ
t
Q
,
M
=
1
+
δ
Q
-
γ
Q
2
,
N
=
1
2
α
Δ
t
P
, and the condition of linear stability, that is, |
U
n
+
1
|
≤
|
U
n
|, is satisfied provided that μ
Δ
t
Q
≥
0 which is indeed the case because Q
>
0 and μ
≥
0. The method is, however, dispersive and the phase of U
n
+
1
/
U
n is given by θ
n
u
m
=
arctan
⁡
-
2
M
N
/
M
2
-
μ
Δ
t
/
2
Q
2
-
N
2 which for k
h
≪
1, μ
k
2
Δ
t
≪
1, and |
α
|
k
Δ
t
≪
1 is tan
⁡
θ
≈
-
α
k
Δ
t.


## 4. Results

In this section, some sample results illustrating the numerical solution to (6) are presented for several values of the parameters that appear in that equation. Unless otherwise stated, the initial condition used in the study is(32)
u
0
,
x
=
A
exp
⁡
-
x
-
30
σ
2
, which corresponds to a Gaussian function of amplitude A
=
1 and standard deviation σ
2
=
20.The numerical experiments reported here were performed witha
=
0 and b
=
100, and the numerical method was stopped whenever the influence of the boundaries on the solution was noticeable. This explains why, in some presented graphs, the time axis is smaller than in other ones.Most of the calculations were performed withN
=
1000, that is, h
=
0.1, and Δ
t
=
0.0001 and correspond to nearly grid-independent results. When either steep gradients or radiation tails were observed, calculations were performed with smaller time steps and grid spacings in order to ensure almost grid-independent results.The accuracy of the results was assessed in terms of the discreteL
1 and L
2 norms of the differences between two solutions corresponding to different time steps and/or grid sizes; they were also assessed in terms of the invariants reported before, but the preservation of invariants was found not to be a good indicator of the accuracy of the method because the invariants represent global quantities that are obtained from a spatial integration of the solution at each time level, and numerical quadrature is a smoothing operation.In many of the results presented here, it was observed that|
I
1
(
30
)
-
I
1
(
0
)
|
/
I
1
(
0
) was less than 1
0
-
4 for h
=
0.1 and Δ
t
=
0.0001, when the invariants were evaluated by means of a second-order accurate trapezoidal rule (cf. (13)). Further comments on the first invariant are made at the end of this section when assessing the effects of the initial conditions on wave formation.Before proceeding with the presentation and discussion of results, it is convenient to analyze the different terms that appear in (6) even though the superposition principle is not applicable to such a nonlinear equation. By considering only the first and second terms of (6), that is, u
t
+
α
u
x
=
0, one obtains the wave solution u
(
t
,
x
)
=
f
(
x
-
α
t
) that indicates that u
(
t
,
x
) remains constant along the characteristic lines x
-
α
t
=
C, where C is a constant and f denotes a function. If only the first and third terms of (6) are considered, that is, u
t
+
4
ϵ
u
3
u
x
=
0, the solution of the resulting equation may be written as u
(
t
,
x
)
=
f
(
x
-
4
ϵ
u
3
t
) provided that shock waves are not formed; since the speed of the waves is equal to 4
ϵ
u
3 and, therefore, increases with u, wave steepening and shock formation may occur if u
x
(
0
,
x
)
<
0.If only the first and fourth terms of (6) are considered, that is, u
t
=
μ
u
x
x, the solution of the resulting equation may be written as u
(
t
,
x
)
=
∫
-
∞
∞
‍
Θ
(
x
-
z
)
f
(
z
)
d
z, where Θ
(
t
,
x
)
=
1
/
4
π
μ
t
e
x
p
-
x
2
/
4
μ
t and u
(
0
,
x
)
=
f
(
x
), and decays in both space and time. On the other hand, if only the first and fifth terms of (6) are accounted for, that is, u
t
=
δ
u
x
x
t, an integration of this equation yields u
=
δ
u
x
x
+
g
(
x
), the integration of which yields (for δ
>
0)(33)
u
t
,
x
=
C
cosh
⁡
x
δ
+
B
sinh
⁡
x
δ
+
δ
∫
g
s
sinh
⁡
x
-
s
δ
d
s
, where C and B are integration constants; if δ
<
0, the solution to u
t
=
δ
u
x
x
t contains trigonometric terms and is not reported here. Note that the homogeneous part of (33) increases as |
x
| increases unless C
=
B
=
0.Finally, if only the first and sixth terms of (6) are considered, that is, u
t
=
γ
u
x
x
x
x
t, an integration of this equation yields u
=
γ
u
x
x
x
x
+
g
(
x
), the integration of which (for γ
>
0) may be written as(34)
u
t
,
x
=
E
cosh
⁡
θ
x
+
B
sinh
⁡
θ
x
+
C
cos
⁡
θ
x
+
D
sin
⁡
θ
x
+
u
p
x
, where E, B, C, and D are integration constants and θ
4
=
γ
-
1, and indicates that u
(
t
,
x
) increases as |
x
| increases unless E
=
B
=
0. The particular solution u
p
(
x
) may be found without much difficulty by means of Lagrange’s variation of parameters, but it is not reported here.Forγ
<
0, the homogeneous solution of u
=
γ
u
x
x
x
x
+
g
(
x
) is(35)
u
h
t
,
x
=
cosh
⁡
θ
x
E
cos
⁡
θ
x
+
B
sin
⁡
θ
x
+
sinh
⁡
θ
x
C
cos
⁡
θ
x
+
D
sin
⁡
θ
x
, where θ
4
=
(
-
4
γ
)
-
1.A comparison between (33) and (34) indicates that the former grows faster than the latter as |
x
| increases if 1
/
δ
>
1
/
γ
1
/
4 (with δ
>
0 and γ
>
0), that is, if γ
≫
δ
2, whereas the opposite holds if γ
≪
δ
2. On the other hand, a comparison between (33) and (35) shows that the former grows faster than the latter if 1
/
δ
>
1
/
2
(
-
γ
)
1
/
4 (with δ
>
0 and γ
<
0), that is, if -
γ
≫
δ
2
/
4, whereas the opposite holds if -
γ
≪
δ
2
/
4.Since, according to the analysis performed in Section2, γ
<
0 for (linear) stability, the contribution of the third-order derivative term of (9) is much more important than that of the fifth-order derivative one if δ
≫
-
γ
k
2; moreover, since, the shortest wave that can be represented in a grid is 2
h and the largest wavenumber that may be resolved is, therefore, π
/
h, the condition δ
≫
-
γ
k
2 requires that δ
≫
-
1000
γ. This implies that δ must be much larger than |
γ
|. For γ
>
0, a similar condition holds.Effect of
μ. Figures 1 and 2 show the solution of (6) for two different values of the viscosity coefficient μ. Figure 1 indicates that the initial Gaussian profile undergoes a large change initially, and its leading edge steepens whereas its amplitude decreases as time increases. The wave profile also widens as time increases; this widening is a consequence of the value of μ considered in Figure 1, while the steepening of the leading front is a consequence of the nonlinear advective terms, as discussed at the beginning of this section. Figure 1 also shows that, for μ
=
0.5, u
(
t
,
x
)
≥
0.Figure 1
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0.5, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure 2
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0.01, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.It should be pointed out that the combs/peaks observed in the amplitude of the leading wave in some figures presented in this paper are due to the fact that only a limited number of grid points and time levels are represented.As the value ofμ is decreased, the effects of viscosity become less important than the advective and dispersive ones, the wave’s leading front steepening increases, and sawtooth waves may form behind the leading front; these waves have a typical N-structure that has been observed, for example, in nonlinear acoustics [18], and correspond to dispersive shock waves [17, 20, 21]. The number and amplitude of these sawtooth waves increase as μ is decreased, and a radiation tail may form as indicated in Figure 2 which corresponds to μ
=
0.01 and the values of the parameters shown in Table 3. Figure 2 also shows the presence of four waves whose amplitude decreases slowly with time.Table 3
Values of the parameters used in the numerical simulations (var. = variable).
Figure
α
ϵ
μ
δ
γ
Figure1
1
1
0.5
0.1
0.00001
Figure2
1
1
0.01
0.1
0.00001
Figure3
1
1
var.
0.1
0.00001
Figure4
0
1
0.005
0.1
0.00001
Figure5
1
1
0.005
0.1
0.00001
Figure6
var.
1
0.005
0.1
0.00001
Figure7
1
0.1
0.005
0.1
0.00001
Figure8
1
0.5
0.005
0.1
0.00001
Figure9
1
var.
0.005
0.1
0.00001
Figure10
1
1
0.005
0.5
-
0.00001
Figure11
1
1
0.005
0.05
-
0.00001
Figure12
1
1
0
0.1
0.00001
Figure13
1
1
0
0.01
0.00001
Figure14
1
1
0
var.
0.00001
Figure15
1
1
0
0.1
var.
Figure16
1
1
0
0.01
0.00001
Figure17
1
0
0
0.1
0.00001
Figure18
0
1
0
0.1
0.00001A similar structure to that of Figure2 has also been observed for μ
=
0, thus indicating that the number of waves increases from 1 to 4 as μ is decreased from unity. The width of these waves increases whereas their height decreases as the viscosity coefficient is increased. Figure 2 also shows that the waves formed from the breakup of the initial Gaussian condition used in this study exhibit curved trajectories. Solitary waves of the EW, RLW, and GRLW have also been found to follow curved trajectories in viscous media [45–47].In Figure3, the u profile is shown at t
=
2 and t
=
20 in order to illustrate the initial adjustment of the wave profile from its initial Gaussian shape, the formation of N- or dispersive shock waves behind the leading front, and the number of waves that are formed as functions of the viscosity coefficient. In particular, Figure 3 shows that the amplitude of the leading wave that emerges first from the initial Gaussian profile has a larger amplitude than those that emerge later.Figure 3
Wave profile att
=
2 (a) and 20 (b) as a function of μ for initial Gaussian conditions with A
=
1, σ
2
=
20, α
=
1, ϵ
=
1, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols ×, ⋄, and ∘ correspond to μ
=
0, 0.1, and 0.05, resp.).
(a)
(b)Effect of
α. Figures 4 and 5 illustrate the effects of the linear advection term on the solution of (6). These figures show that there is an initial adjustment of the u profile that generates a steep leading front that becomes a solitary wave at later times. The u profile behind the trailing edge of the leading wave also steepens and results in the formation of another wave that travels at a smaller speed than the leading one. Behind the second leading wave, a third wave is formed and this process continues as indicated in Figures 4 and 5 until the nonlinear advective terms are smaller than the linear ones. The amplitude of the leading wave is larger than those of the trailing ones; the amplitude of the latter decreases as the distance to the upstream boundary decreases.Figure 4
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
0, ϵ
=
1, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure 5
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure4 also shows (cf. Figures 1 and 2) that the trailing waves are initially of the N- or dispersive shock wave type and are caused by the steepening of the u profile associated with the nonlinear advective term and dispersion. This is more clearly visible in Figure 5 that corresponds to a linear convective field with α
=
1.Asα is increased, the speed of the leading wave front increases but the number of waves that emerge from the breakup of the initial Gaussian profile decreases; this effect is due to the pushing effect associated with the linear advective term. Although not shown here, similar results to those presented in Figures 4 and 5 have been observed for γ
=
0 and the same values of the parameters as those of Figures 4 and 5.Although not shown here, numerical experiments performed with the same values of the parameters as those of Table3 but with α
<
0 show similar behavior to that observed in Figures 4 and 5, except that the leading wave speed increases and the leading wave’s front becomes steeper as α is increased. This is illustrated in Figure 6 that shows the u profile at t
=
2 and 20 for positive and negative values of α and indicates that the amplitude of the leading wave is greater than those of the trailing ones. Figure 6 also shows that as the magnitude of α for α
<
0 is increased, it takes a longer time to break up the initial Gaussian profile due to the opposing linear convective term.Figure 6
Wave profile att
=
2 (a) and 20 (b) as a function of α for initial Gaussian conditions with A
=
1, σ
2
=
20, ϵ
=
1, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols ×, ∘, ⋄, +, and ∗ correspond to α
=
1, 0.5, 0, −1, and −0.5, resp.).
(a)
(b)Effect of
ϵ. The effects of the nonlinear advective term on wave breakup and propagation are illustrated in Figures 7 and 8. For ϵ
=
0, (6) is linear and may be solved analytically by means of the Fourier transform in x. This solution is analogous to the initial Gaussian condition and moves at a speed approximately equal to α, and its amplitude decreases while its width increases for μ
≠
0. The numerical solution obtained with the compact operator method presented in this paper may, however, exhibit a small oscillation at the wave’s trailing edge if the number of grid points is not enough to ensure an adequate spatial resolution there; such an oscillation does not occur if the trailing edge is properly resolved.Figure 7
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
0.1, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure 8
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
0.5, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Asϵ is increased, the magnitude of the nonlinear convective terms also increases; this results in a steepening of the leading front and the formation of an N- or dispersive shock wave as illustrated in Figure 7 which corresponds to ϵ
=
0.1. This figure also shows a small radiation tail behind the N-wave.Forϵ
=
0.5, the results presented in Figure 8 show that the initial Gaussian profile generates three solitary waves, an N- or dispersive shock wave, and some radiation which are clearly visible at t
=
30. This figure also indicates that the amplitude of the waves decreases from the leading to the trailing one, and the amplitude of each wave also decreases slowly with time due to the value of μ
=
0.005 employed in the simulations.In Figure9, u
(
2
,
x
) and u
(
20
,
x
) are presented as functions of x for several values of ϵ. This figure shows that, for ϵ
=
0, the profiles of u are almost identical to the initial one and that the wave propagates at a speed approximately equal to unity, as discussed previously. The u
(
2
,
x
) profiles corresponding to ϵ
=
0.5 and 1 exhibit similar trends but the amplitude of the secondary wave is smaller for ϵ
=
1 than for ϵ
=
0.5; this is due to the fact that the wave steepening and wave breakup increase as ϵ is increased as discussed at the beginning of this section. The u profiles at t
=
20 show that the number of dispersive shock waves that result from the breakup of the initial Gaussian condition increases as ϵ is increased.Figure 9
Wave profile att
=
2 (a) and 20 (b) as a function of ϵ for initial Gaussian conditions with A
=
1, σ
2
=
20, α
=
1, μ
=
0.005, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols ×, ∘, ⋄, +, and ∗ correspond to ϵ
=
1, 0.5, 0, −1, and −0.5, resp.).
(a)
(b)Effect of
δ. Figures 10 and 11 illustrate the effects of δ here referred to as the first dispersion coefficient on wave breakup. For δ
=
0.5, the initial Gaussian profile splits into two right-propagating waves; the leading one propagates at a higher speed than the trailing one, and some radiation is observed between the upstream boundary, that is, x
=
a, and the trailing wave. Figure 10 also shows the initial steepening of the front of the initial Gaussian condition and the smaller slope of the leading edge of the trailing wave.Figure 10
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0.005, δ
=
0.5, γ
=
-
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure 11
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0.005, δ
=
0.05, γ
=
-
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Forδ
=
0.05, the results presented in Figure 11 indicate that six waves are formed; the leading wave propagates at a faster speed than the trailing ones, and the speed of the latter decreases from the downstream to the upstream boundary. Figure 11 also shows that the leading wave follows a curved trajectory, and the curvature of the trailing waves decreases from the downstream to the upstream boundary. The separation between waves also decreases from the downstream to the upstream boundaries.Although not shown here, it has been observed that five waves are present att
=
30 for δ
=
0.1, thus indicating that the number of waves generated from the breakup of the initial Gaussian distribution increases as δ is decreased from unity. Furthermore, for the same values of the parameters as those for Figures 10 and 11 but δ
=
0.005, it has been observed that the leading part of the initial Gaussian condition steepens as a consequence of the nonlinear convective terms and results in large gradients of u at the wave front. The interaction among the steepening associated with the nonlinear convective terms and the third- and fifth-order derivative ones results in the formation of sawtooth or dispersive shock waves in the trailing part of the leading wave. The amplitude of the sawtooth waves decreases as the gradients of u decrease and no sawtooth waves are observed behind the last trailing wave.Although not shown here, similar results to those reported in Figures10 and 11 have been observed for the same parameters as those of these figures except that γ
=
0.00001, thus indicating that the fifth-order derivative term that appears in (6) does not play an important role in wave propagation for α
=
ϵ
=
1, μ
=
0.005, -
0.00001
≤
γ
≤
0.00001, and δ
≥
0.01.For the same parameters as those for Figures10 and 11 but smaller values of δ, the contribution of the third-order derivative or first dispersive term decreases, and the number of sawtooth waves formed behind the leading wave’s trailing edge decreases and their amplitude increases as δ is decreased from 0.01.As indicated at the beginning of this section, for the initial Gaussian conditions considered in this study, the nonlinear convective terms result in the formation of a shock wave forμ
=
δ
=
γ
=
0 and a very steep front for small μ
≠
0 and δ
=
γ
=
0; therefore, it may be concluded that the sawtooth waves observed for small value of δ are caused by the fifth-order derivative term that appears in (6). Note that, for the EW, RLW, and GRLW equations [45–47] which do not contain u
x
x
x
x
t in (6), no sawtooth waves are observed.The effects of bothμ and δ on wave propagation are shown in Figures 12 and 13 that should be compared with Figures 1 and 2 and Figures 10 and 11, respectively, in order to observe the effects of these parameters. For μ
=
0 and δ
=
0.1, four waves can be seen at t
=
25; the fourth one is followed by an N-wave and some radiation. The amplitude of succesive waves increases and the separation between them decreases from the downstream to the upstream boundary, and their trajectories are straight lines compared with the curved ones observed in Figure 2; this indicates once again that not only does the viscosity decrease the wave amplitude, it also increases the curvature of the waves’ trajectories. Furthermore, for δ
=
0.01 and μ
=
0, the results presented in Figure 13 indicate that twelve waves followed by an N one and some radiation are observed at t
=
30; these waves move along straight lines and their speeds decrease from the downstream to the upstream boundary.Figure 12
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0, δ
=
0.1, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.Figure 13
Wave propagation for initial Gaussian conditions withA
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0, δ
=
0.01, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.A summary of the results discussed in this subsection is presented in Figure14 that shows the u
(
t
,
x
) profiles at t
=
2 and 20 for several values of δ and indicates that the number of waves that break up from the initial Gaussian one increases as δ is decreased from unity. For δ
=
1 and 0.5, only two right-propagating waves and a small oscillatory tail are observed. Figure 14 also shows the initial steepening and wave breakup at t
=
2.Figure 14
Wave profile att
=
2 (a) and 20 (b) as a function of δ for initial Gaussian conditions with A
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols ×, ∘, ⋄, +, and ∗ correspond to δ
=
1, 0.5, 0.1, 0.05, and 0.01, resp.).
(a)
(b)Effect of
γ. As discussed previously in this paper, the linear stability of (6) places some limitations on the sign of γ which is here referred to as the second dispersion coefficient. Since the numerical calculations are performed in finite domains with a finite number of grid points and the largest wavenumber that can be resolved is π
/
h, small wave instabilities have been observed for values of γ that do not satisfy the condition described at the beginning of this section. This condition implies that |
γ
| should be smaller than approximately one-thousandth of δ. For such small values of γ, very few differences in wave propagation were observed as indicated in Figure 15 that shows u
(
2
,
x
) and u
(
20
,
x
) for several values of γ. This figure indicates that the breakup of the initial Gaussian condition and the waves generated in its breakup are nearly independent of γ.Figure 15
Wave profile att
=
2 (a) and 20 (b) as a function of γ for initial Gaussian conditions with A
=
1, σ
2
=
20, α
=
1, ϵ
=
1, μ
=
0, δ
=
0.1, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols ×, ∘, ⋄, +, and ∗ correspond to γ
=
0, 1
0
-
4, 1
0
-
5, -
1
0
-
5, and -
1
0
-
4, resp.).
(a)
(b)Effect of the Initial Conditions. For the same initial Gaussian condition as (32) and A
=
1, but σ
=
1
/
0.1, Figure 16 shows that nine solitary waves and a radiation tail are present at t
=
20; the amplitude of these waves remains constant after an initial transient, and their speeds decrease while the separation between successive waves increases from the downstream to the upstream boundary. As σ is decreased (for A
=
1), that is, as the Gaussian initial condition becomes narrower and, therefore, its slope increases, the number of waves also increases.Figure 16
Wave propagation for initial Gaussian conditions withA
=
1, σ
=
1
/
0.1, α
=
1, ϵ
=
1, μ
=
0, δ
=
0.01, γ
=
0.00001, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001.For the initial Gaussian condition of (7), the first invariant is I
1
=
π
A
σ which clearly increases with the amplitude and standard deviation; furthermore, the largest slope of the Gaussian condition occurs at x
-
x
0
=
±
σ
/
2 and its absolute value is A
/
σ
2
/
e. This implies that the slope of the initial condition increases as A is increased and as σ is decreased, and an increase in slope results in an increase of the nonlinear advective terms and a consequent larger steepening of the propagating wave initially, as discussed at the beginning of this section. Such a steepening is eventually balanced by dispersion and results in the formation of a leading wave as illustrated, for example, in Figure 2 that propagates towards the downstream boundary; however, steepening may still occur behind the leading wave due to the nonlinear advective terms until they are again balanced by dispersion and a second solitary wave emerges and propagates towards the downstream boundary. This process continues until the slope of the u profile is so small that the nonlinearities become smaller than dispersion; when this occurs, dispersion dominates and oscillatory tails such as the ones shown, for example, in Figure 2, may be observed.In order to further assess the effects of the initial conditions on both the first invariant and the wave breakup and propagation, numerical experiments were performed with the values ofA and σ shown in Table 4 and the initial condition of (32) and the results are illustrated in Figures 17 and 18. Figure 17 indicates that the first invariant is very well preserved, except for (
A
,
σ
-
2
)
=
(
1,1
) which corresponds to the largest slope of the initial conditions shown in Table 4. As discussed above, the initial steepness of the wave propagation and the number of solitary waves that emerge from the initial Gaussian profile increase as the slope of the initial condition increases as indicated in Table 4 and Figure 18; therefore, it may be concluded that the preservation of I
1 depends strongly on the numerical resolution because as stated above, the number of solitary waves generated from the breakup of the initial conditions increases as the slope of the initial condition is increased (cf. Figure 18(a)) and an adequate number of grid points is needed to accurately resolve these solitary waves; the number of grid points must be increased as A is increased and/or σ is decreased.Table 4
Values of the parameters and number of solitary waves (NS) att
=
25 for initial Gaussian conditions.
Figure18
σ
-
2
A
NS
Left
0.1
1
4
Left
0.05
1
5
Left
0.0125
1
8
Right
0.01
0.1
1
Right
0.04
0.2
1
Right
1
1
2Figure 17
First invariant as a function of time for initial Gaussian conditions withα
=
1, ϵ
=
0, μ
=
0, δ
=
0.1, γ
=
1
0
-
5, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols + and 1, ∗ and 2, and × and 3 correspond to (
A
,
σ
-
2
)
=
(
1,0.1
), (
1,0.05
), and (
1,0.0125
), resp., for (a) and (
0.1,0.01
), (
0.2,0.04
), and (
1,1
), resp., for (b)).
(a)
(b)Figure 18
Amplitude profile att
=
25 as a function of space for initial Gaussian conditions with α
=
0, ϵ
=
1, μ
=
0, δ
=
0.1, γ
=
1
0
-
5, x
0
=
30, a
=
0, b
=
100, N
=
1000, and Δ
t
=
0.0001 (the symbols +, ∗, and × correspond to (
A
,
σ
-
2
)
=
(
1,0.1
), (
1,0.05
), and (
1,0.0125
), resp., for (a) and (
0.1,0.01
), (
0.2,0.04
), and (
1,1
), resp., for (b)).
(a)
(b)Figure18(b) also indicates that, for the same value of I
1
(
0
), the number of solitary waves increases as A is increased and, therefore, as the slope of the initial profile increases; on the other hand, for a fixed value of A, an increase of σ results in an increase of I
1
(
0
) and a decrease of the initial slope. Figure 18(a) indicates that the number of solitary waves that are formed from the breakup of the initial Gaussian profile increases as σ is increased and, therefore, as I
1
(
0
) is increased, even though the slope of the initial profile decreases as σ is increased.Figure18(b) clearly shows that the initial steepening of the leading front increases as A is increased and as σ is decreased; the results shown in Figure 18(a) correspond to different values of I
1. When A and σ are varied so that I
1 is the same, the results presented in Figure 18(b) show that wave breakup only occurs for (
A
,
σ
)
=
(
1,1
) for which the maximum slope of the initial condition is unity. For the other values of A and σ considered in Figure 18(b), I
1
=
π, but both the initial amplitude and the initial slope of the u profile are so small that the nonlinear advective terms are much smaller than the linear ones, dispersion is much more important than the nonlinearities, and, therefore, no wave breakup takes place.The results described heretofore correspond to the initial Gaussian condition of (32). Numerical experiments have also been performed with other initial conditions such as, for example,(36)
u
0
,
x
=
A
c
o
s
h
-
2
ν
x
-
30 and show similar trends to the ones described above, that is, an initial steepening of the leading front of the initial condition, the formation of a leading solitary wave, and so forth, provided that the slope of the initial condition is larger than a critical value that depends on A and ν. In addition, the amplitude of the solitary waves formed after the breakup of the initial Gaussian condition decreases whereas the separation between successive waves and the wave’s width increase from the downstream to the upstream boundary.


## 5. Conclusions

Wave generation from a generalized RLW-Rosenau equation subject to initial Gaussian conditions has been studied numerically by means of a second-order accurate finite difference method in time that employs fourth-order accurate finite difference discretizations for the first-, second-, and fourth-order spatial derivatives and is linearly stable, as a function of the linear and nonlinear advective terms and the viscosity and the two dispersion coefficients. It has been shown that the value of the viscosity does not only affect the amplitude and width of the waves; it also affects their curvature and the number of waves that are generated.The linear advective term was found to push the waves towards either the downstream or upstream boundaries, depending on its direction, whereas the nonlinear advective terms cause a steepening of the leading front and the formation of sawtooth waves and a radiation tail. The number of waves generated by nonlinearities depends not only on their magnitude but also on those of the viscosity and the two dispersion terms.The magnitude of the coefficient of the second dispersion term, that is, that associated with the fourth-order spatial derivative, was found to be limited by stability considerations and provided exponential solutions whose growth may be larger than those associated with the dispersion term that contains a second-order spatial derivative. Such a behavior is not observed if the second dispersion coefficient is much smaller than the first one, but, in this case, it has been found that the effects of the second dispersion term are small.It has also been found that the slope of the initial conditions plays a paramount role in determining the wave formation/breakup and number of solitary waves and that narrower Gaussian conditions result in more solitary waves than wider ones provided that the largest slope of the initial condition is sufficiently large to cause wave steepening and the first dispersion term is large enough so that a balance between nonlinearities and dispersion is reached and a solitary wave may form. In the absence of dispersion, the nonlinearities result in shock wave formation, whereas, in the absence of nonlinearities, neither solitary waves nor shock waves are formed.For the same value of the first invariant, it has been observed that solitary waves may only form when the slope of the initial conditions exceeds a threshold value. Below this value, the nonlinear advective terms are smaller than the linear ones and no solitary wave emerges from the initial conditions.It has also been found that the numerical method used in this study preserves very well the first invariant and predicts with a great accuracy the decay of the second one for negative values of the second dispersion coefficient in viscous media.



---

*Source: 4618364-2016-01-06.xml*