# Propagation of Water Waves over Uneven Bottom under the Effect of Surface Tension

Authors: Juan Carlos Muñoz Grajales
Journal: International Journal of Differential Equations
 (2015)

---

## Abstract

We establish existence and uniqueness of solutions to the Cauchy problem associated with a new one-dimensional weakly-nonlinear, weakly-dispersive system which arises as an asymptotical approximation of the full potential theory equations for modelling propagation of small amplitude water waves on the surface of a shallow channel with variable depth, taking into account the effect of surface tension. Furthermore, numerical schemes of spectral type are introduced for approximating the evolution in time of solutions of this system and its travelling wave solutions, in both the periodic and nonperiodic case.


---

## Body

## 1. Introduction

In this paper we study the propagation of water waves on the surface of a shallow channel with variable depth, considering the effect of surface tension. To describe this phenomenon we will derive a new water wave model from Euler’s equations (in dimensionless variables) for an inviscid, incompressible liquid bounded above by a free surface and bounded below by an impermeable bottom topography [1]:(1)
β
ϕ
x
x
+
ϕ
y
y
=
0
f
o
r
-
H
x
γ
<
y
<
α
η
x
,
t
,
-
∞
<
x
<
∞
,with the nonlinear free surface conditions(2)
η
t
+
α
ϕ
x
η
x
-
1
β
ϕ
y
=
0
,
η
+
ϕ
t
+
α
2
ϕ
x
2
+
1
β
ϕ
y
2
-
β
σ
∂
x
η
x
1
+
β
α
2
η
x
2
=
0
,at y
=
α
η
(
x
,
t
). Here ϕ
(
x
,
y
,
t
) denotes the potential velocity and η
(
x
,
t
) the wave elevation measured with respect to the undisturbed free surface y
=
0. The dimensionless parameters α and β are small positive real numbers which measure the strength of nonlinear and dispersive effects, respectively. The parameter γ measures the ratio inhomogeneities/wavelength and the parameter σ is associated with the surface tension. The Neumann condition at the impermeable bottom is(3)
ϕ
y
+
β
γ
H
′
x
γ
ϕ
x
=
0
.The bottom topography is described by y
=
-
H
(
x
/
γ
), where(4)
H
x
γ
=
1
+
n
x
γ
,
w
h
e
n
0
<
x
<
L
1
,
w
h
e
n
x
≤
0
o
r
x
≥
L
. We will see in Section 3 that system (1)–(3) is equivalent to leading order provided that α
≪
1, β
≪
1 to the weakly-nonlinear, weakly-dispersive system(5)
I
-
β
2
∂
ξ
2
v
t
=
2
β
3
∂
ξ
4
Φ
-
∂
ξ
1
+
α
v
M
2
Φ
ξ
,
I
-
β
2
∂
ξ
2
Φ
t
=
β
σ
v
ξ
ξ
-
3
M
′
M
4
v
ξ
+
3
M
′
2
M
5
v
-
M
′
′
M
4
v
-
1
M
v
-
α
2
M
2
Φ
ξ
2
,
v
ξ
,
0
=
v
0
ξ
,
Φ
ξ
,
0
=
Φ
0
ξ
.The variables ξ and t are the space and time coordinates, respectively, v
=
M
(
ξ
)
η
(
ξ
,
t
) with η
(
ξ
,
t
) denoting the elevation of the free surface, and Φ
(
ξ
,
t
) represents the potential fluid velocity measured at the channel’s bottom. The metric coefficient M
(
ξ
) is related to the channel’s bottom y
=
-
H
(
x
/
γ
) and it is defined as(6)
M
ξ
=
π
4
β
∫
-
∞
∞
H
x
ξ
0
,
-
β
/
γ
c
o
s
h
2
π
/
2
β
ξ
0
-
ξ
d
ξ
0
,where (
x
,
y
)
→
(
ξ
,
ζ
~
) is the coordinate transformation used in the derivation of system (5) to map the original physical channel with variable depth onto a strip in the complex plane. This strategy of change of variable was introduced by Hamilton [2] and later used successfully by Nachbin in [3] to study wave propagation over a channel with a highly variable topography. Observe from (6) that the coefficient M
(
ξ
) is infinitely differentiable although the function describing the channel’s bottom is not smooth. In case of constant depth, the coefficient M
(
ξ
) is identically one and system (5) reduces to a system derived by Quintero and Montes in [4].Wave-topography interaction has been the subject of considerable mathematical research [5–18]. The physical applications range from coastal surface waves [19] to atmospheric flows over mountain ranges [20, 21]. In particular, the interaction of waves with fine features of the topography is of great interest. As pointed out in the introduction to theOrography proceedings [21] of the European Centre for Medium-Range Weather Forecasts (ECMWF), “the representation … of subgrid-scale orographic processes is recognized as crucial to numerical weather prediction at all time ranges.” In the atmospheric literature,orography implies mountain ranges [20].In previous works, some weakly-nonlinear, weakly-dispersive models have been developed to describe the interaction of a long pulse with small amplitude that propagates on the surface of a channel with a variable bottom [2, 3, 5, 7–9, 22–24]. However, these models either neglect the effect of surface tension on the free surface where the wave propagates or are not applicable to bottoms described by a discontinuous/nondifferentiable function or the wave elevation is removed in their physical derivation. For instance, Milewski [23] derived a bidirectional scalar Benney-Luke type equation in terms of the potential velocity which includes the surface tension effect and the influence of a variable bottom. However, the asymptotical derivation of this model implies eliminating the wave elevation and consequently neglecting several second-order terms in the parameters α and β.Some of the features of new formulation (5) are the following:(i)
It can be applied to study wave propagation over a shallow channel with a discontinuous or nondifferentiable bottom, provided that the bottom’s fluctuations satisfy|
n
|
<
1. This is a consequence of introducing the conformal mapping (
x
,
y
)
→
(
ξ
,
ζ
~
) mentioned above. Observe that all coefficients of the reduced equations (5) in the new coordinate system result in being infinitely differentiable since they depend on the smooth function M
(
ξ
).
(ii)
One-dimensional system (5) models bidirectional waves and it incorporates the simultaneous effects of surface tension and variable depth upon the shape of a water wave that propagates on the surface of an irregular shallow channel.
(iii)
Furthermore, in the derivation of system (5), we do not eliminate the wave elevation (which prevents us from neglecting additional terms of order O
(
α
2
,
β
2
)), as done, for example, in [23]. Therefore, the new formulation (5) is expected to be a more accurate approximation of the full potential theory equations.In the present paper, we establish existence and uniqueness of a solution to Cauchy problem (5) using classical semigroup theory and Banach’s fixed point principle. In [4] the well-posedness of system (5) was analyzed but only in the case that M
(
ξ
)
≡
1. In second place, we formulate a Galerkin-spectral numerical scheme of spectral accuracy in space to approximate the solutions of system (5) based on the Fourier basis and using an implicit-explicit (IMEX) second-order strategy for time stepping. This type of temporal discretization is described, for instance, in [25], and it has been used in conjunction with spectral methods [26, 27] for the time integration of spatially discretized PDEs of diffusion-convection type. IMEX schemes have also been successfully applied to the incompressible Navier-Stokes equations [28] and in environmental modelling studies [29]. On the other hand, we develop a numerical solver to compute travelling wave solutions of system (5) by using a Fourier-collocation strategy combined with a Newton-type iterative procedure. We also indicate how to determine appropriate starting points for this iterative process in order to achieve convergence. Existence results on travelling wave solutions (in the periodic and nonperiodic case) of system (5) with M
≡
1 have also been established in [4]. Travelling wave solutions exist as a consequence of a balance between nonlinear and dispersive effects present in a system; these waves travel with a constant speed without any temporal evolution in shape or size when the frame of reference moves with the same speed of the wave. In the last decades, the study of travelling waves has grown enormously because they appear in several and varied fields of application, such as fluid mechanics, optics, acoustics, oceanography, and astronomy. Thus, to determine existence and properties of such type of solutions is a fundamental problem in the theory of ordinary and partial differential equations of great interest for both pure and applied mathematicians.The rest of this paper is organized as follows. In Section2, we introduce the functional spaces and notation to be employed in the paper. In Section 3, we present in detail the derivation of model (5) starting from the potential theory equations including the surface tension effect. In Section 4, we discuss existence and uniqueness of a solution of the Cauchy problem associated with system (5) by using Banach’s fixed point principle and semigroup theory. In Section 5, we introduce the numerical schemes for approximating the evolution of a solution of system (5) and computing their travelling wave solutions (periodic and nonperiodic cases). Section 6 presents a set of numerical simulations to check the accuracy of the numerical schemes developed in the paper. Furthermore, some numerical experiments are included to illustrate the interaction between the wave-topography effects and surface tension. Finally in Section 7 are the conclusions of the paper.


## 2. Preliminaries

To analyze existence of solutions of problem (5), we will use the standard notation. For 1
≤
p
≤
∞, we will denote by L
p
(
R
) (or simply L
p) the Banach space of measurable functions in R such that ∫
R
|
f
(
x
)
|
p
d
x
<
∞ if 1
≤
p
<
∞ and e
s
s
s
u
p
R
|
f
|
<
∞, if p
=
∞. We define the norm in L
p
(
R
) for 1
≤
p
<
∞ by(7)
f
L
p
=
∫
R
f
x
p
1
/
pand in L
∞
(
R
) by f
∞
=
e
s
s
s
u
p
R
|
f
|. L
2
(
R
) is a Hilbert space for the scalar product(8)
f
,
g
=
∫
R
f
x
g
x
¯
d
x
.We set f
=
f
L
2. For function f
∈
L
1
(
R
), the Fourier transform is defined as(9)
F
f
y
=
f
^
y
=
∫
R
f
x
e
-
i
x
y
d
x
,
y
∈
R
,and the inverse Fourier transform is defined by(10)
F
-
1
f
y
=
f
ˇ
y
=
1
2
π
∫
R
f
x
e
i
x
y
d
x
,
y
∈
R
.We will also denote by F
(
f
) (or f
^) and F
-
1
(
f
) (or f
ˇ) the extensions of these operators to L
2
(
R
). The convolution of two functions, f
,
g
∈
L
2
(
R
), is defined as(11)
f
∗
g
x
=
∫
R
f
x
-
y
g
y
d
y
.We recall that f
∗
g
^
(
y
)
=
f
^
(
y
)
g
^
(
y
). For s
∈
R, we define the Sobolev space H
s
(
R
) (sometimes written for simplicity as H
s) as the completion of the Schwartz space (rapidly decaying functions) defined as (12)
S
R
≔
f
∈
C
∞
R
:
x
ν
∂
x
β
f
∞
<
∞
,
f
o
r
a
n
y
ν
,
β
∈
Z
+
,with respect to the norm(13)
f
H
s
=
∫
R
1
+
y
2
s
f
^
y
2
d
y
1
/
2
.For simplicity, we also denote this norm by f
s. The inner product in H
s
(
R
) is defined as(14)
f
,
g
s
=
∫
R
1
+
y
2
s
f
^
y
g
^
y
¯
d
y
.The product norm in the space H
s
(
R
)
×
H
s
(
R
) is defined by (15)
ζ
,
u
H
s
×
H
s
=
ζ
H
s
+
u
H
s
,for (
ζ
,
u
)
∈
H
s
(
R
)
×
H
s
(
R
). Sometimes, we will also use the equivalent product norm(16)
ζ
,
u
H
s
×
H
s
=
ζ
H
s
2
+
u
H
s
2
1
/
2
.For T
>
0, we will denote by C
(
[
0
,
T
]
,
H
s
(
R
)
) the space of continuous functions f
:
[
0
,
T
]
→
H
s
(
R
), that is, the space of continuous functions t
→
f
(
t
,
·
)
∈
H
s
(
R
), t
∈
[
0
,
T
], with the supremum norm (17)
f
C
0
,
T
=
s
u
p
t
∈
0
,
T
f
t
,
·
H
sand the product norm(18)
ζ
,
u
C
0
,
T
2
=
ζ
C
0
,
T
+
u
C
0
,
T
,for (
ζ
,
u
)
∈
C
(
[
0
,
T
]
,
H
s
(
R
)
×
H
s
(
R
)
).


## 3. Governing Equations

We start by presenting the potential theory formulation for Euler’s equations (in dimensionless variables) for an inviscid, incompressible liquid bounded above by a free surface and bounded below by an impermeable bottom topography and including the effect of surface tension [1]:(19)
β
ϕ
x
x
+
ϕ
y
y
=
0
f
o
r
-
H
x
γ
<
y
<
α
η
x
,
t
,
-
∞
<
x
<
∞
,with the nonlinear free surface conditions(20)
η
t
+
α
ϕ
x
η
x
-
1
β
ϕ
y
=
0
,
η
+
ϕ
t
+
α
2
ϕ
x
2
+
1
β
ϕ
y
2
-
β
σ
∂
x
η
x
1
+
β
α
2
η
x
2
=
0
,at y
=
α
η
(
x
,
t
). Here ϕ
(
x
,
y
,
t
) denotes the potential velocity and η
(
x
,
t
) the wave elevation measured with respect to the undisturbed free surface y
=
0. The dimensionless parameters α and β measure the strength of nonlinear and dispersive effects, respectively. The parameter γ measures the ratio inhomogeneities/wavelength and σ is related to the surface tension effects. The Neumann condition at the impermeable bottom is(21)
ϕ
y
+
β
γ
H
′
x
γ
ϕ
x
=
0
.The bottom topography is described by y
=
-
H
(
x
/
γ
), where(22)
H
x
γ
=
1
+
n
x
γ
,
w
h
e
n
0
<
x
<
L
1
,
w
h
e
n
x
≤
0
o
r
x
≥
L
.The bottom profile is described by the (possibly rapidly varying) function -
n
(
x
/
γ
). The topography is rapidly varying when γ
≪
1. The scale L represents the total length of the irregular section of the coast. The undisturbed depth is given by y
=
-
1 and the topography can be of large amplitude provided that |
n
|
<
1. The fluctuations n are not assumed to be small, nor continuous, nor slowly varying.Let us consider a symmetric flow domain by reflecting the original one about the undisturbed free surface (cf. Figure1). This domain is denoted by Ω
z where z
=
x
+
i
β
y and can be considered as the conformal image of the strip Ω
w where w
=
ξ
+
i
ζ
~ with |
ζ
~
|
≤
β. Then z
=
x
(
ξ
,
ζ
~
)
+
i
β
y
(
ξ
,
ζ
~
)
=
x
(
ξ
,
ζ
~
)
+
i
y
~
(
ξ
,
ζ
~
) with x and y
~, a pair of harmonic functions on Ω
w. Following the strategy suggested by Hamilton in [2] and Nachbin [3] within the weakly-nonlinear, weakly-dispersive regime (α
≪
1, β
≪
1) and using the relationships (23)
ϕ
x
=
1
J
y
~
ζ
~
ϕ
ξ
-
y
~
ξ
ϕ
ζ
~
,
ϕ
y
~
=
1
J
-
x
ζ
~
ϕ
ξ
+
x
ξ
ϕ
ζ
~
,
ϕ
x
2
+
ϕ
y
~
2
=
1
J
ϕ
ξ
2
+
ϕ
ζ
~
2
,
∂
x
=
1
J
y
~
ζ
~
∂
ξ
-
y
~
ξ
∂
ζ
~
,where(24)
J
=
x
ξ
y
~
ζ
~
-
y
~
ξ
x
ζ
~
=
y
~
ζ
~
2
+
y
~
ξ
2
,and the variable free surface coefficient M
(
ξ
) defined as(25)
M
ξ
≔
y
~
ζ
~
ξ
,
0
=
1
+
π
4
β
∫
-
∞
∞
n
x
ξ
0
,
-
β
/
γ
c
o
s
h
2
π
/
2
β
ξ
0
-
ξ
d
ξ
0
,the potential theory equations can be approximated to order O
(
α
,
β
) in the orthogonal curvilinear coordinates (
ξ
,
ζ
) with ζ
=
ζ
~
/
β
+
1 by the equation(26)
β
ϕ
ξ
ξ
+
ϕ
ζ
ζ
=
0
,
a
t
0
<
ζ
<
1
+
α
N
ξ
,
t
,with conditions at the free surface ζ
=
1
+
α
N
(
ξ
,
t
),(27)
J
N
t
+
α
ϕ
ξ
N
ξ
-
1
β
ϕ
ζ
=
0
,
η
+
ϕ
t
+
α
2
J
ϕ
ξ
2
+
1
β
ϕ
ζ
2
-
β
σ
1
J
M
∂
ξ
1
M
η
ξ
=
O
α
2
,
α
β
,
β
2and condition at the channel’s bottom(28)
ϕ
ζ
=
0
,
a
t
ζ
=
0
.Observe that the change of variables ζ
~
→
ζ lets the origin of the curvilinear coordinate system at the bottom. The Jacobian for the (
ξ
,
ζ
~
)
→
(
x
,
y
~
) coordinate transformation is represented by |
J
|, and ζ
=
1
+
α
N
(
ξ
,
t
) corresponds to the position of the free surface in the curvilinear coordinate system. By performing an asymptotic simplification as in [1] (page 464) through a power series expansion in terms of the dispersion parameter β near the bottom of the channel in the form(29)
ϕ
ξ
,
ζ
,
t
=
∑
j
=
0
k
-
1
j
β
j
ζ
2
j
2
j
!
∂
ξ
2
j
Φ
ξ
,
t
+
O
β
k
+
1
,into (26)–(28) for α
,
β being small, we find that free surface conditions (27) can be approximated to order O
(
α
), O
(
β
) by the equations(30)
M
ξ
η
t
+
1
+
α
M
ξ
η
Φ
ξ
ξ
-
β
6
Φ
ξ
ξ
ξ
ξ
=
O
α
2
,
α
β
,
β
2
,
(31)
η
+
Φ
t
-
β
2
Φ
ξ
ξ
t
+
α
2
M
2
ξ
Φ
ξ
2
-
β
σ
1
M
∂
ξ
1
M
η
ξ
=
O
α
β
,
β
2
.Here Φ
(
ξ
,
t
)
=
ϕ
(
ξ
,
0
,
t
) denotes the potential fluid velocity at the channel’s bottom ζ
=
0. We point out that (30) and (31) with σ
=
0 (no surface tension) correspond to those derived in [3] (bottom of page 915 and top of page 916). In the derivation of (30)-(31), we used the relationship(32)
J
ξ
,
t
=
M
ξ
2
+
O
α
2
,which means that, at leading order, the Jacobian of the conformal coordinate transformation is time independent.Figure 1
The symmetric domain in the complexz-plane, where z
=
x
(
ξ
,
ζ
~
)
+
i
y
~
(
ξ
,
ζ
~
). The lower half (x
∈
[
-
5,5
]
,
y
∈
[
-
3,0
]) is the physical channel with y
=
ζ
~
=
0 indicating the undisturbed free surface. Superimposed in this complex z-plane domain are the (curvilinear) coordinate level curves from the w-plane system ξ
ζ
~. The polygonal line at the bottom of the figure is a schematic representation of the topography (where ζ
~
=
±
β). This figure was generated using SC Toolbox [30].Let us introduce the new variablev
(
ξ
,
t
)
:
=
M
(
ξ
)
η
(
ξ
,
t
). Thus observe that from (30) (33)
v
t
=
M
η
t
=
-
∂
ξ
2
Φ
+
O
α
,
β
.This relationship allows us to change the form of dispersive terms in the equations above. In particular, by using the decomposition (34)
-
β
6
∂
ξ
4
Φ
=
β
2
∂
ξ
4
Φ
-
2
β
3
∂
ξ
4
Φ
=
-
β
2
∂
ξ
2
v
t
-
2
β
3
∂
ξ
4
Φ
+
O
α
β
,
β
2
,in (30), we obtain that system (30)-(31) can be approximated to order O
(
α
,
β
) by the following equations:(35)
I
-
β
2
∂
ξ
2
v
t
=
2
β
3
∂
ξ
4
Φ
-
∂
ξ
1
+
α
v
M
2
Φ
ξ
,
I
-
β
2
∂
ξ
2
Φ
t
=
β
σ
1
M
3
v
ξ
ξ
-
3
M
′
M
4
v
ξ
+
3
M
′
2
M
5
v
-
M
′
′
M
4
v
-
1
M
v
-
α
2
M
2
Φ
ξ
2
.We point out that the system above is applicable to modelling of propagation of water waves over an arbitrary rapidly varying depth. In case of a slowly varying channel’s topography and n
(
x
)
=
β
n
0
(
x
) with n
0 being a stationary random process with standard deviation and correlation length of order one, the metric coefficient M
(
ξ
) can be expanded as (see [9])(36)
M
ξ
=
1
+
β
n
0
ξ
+
O
β
.Thus(37)
β
σ
1
M
3
v
ξ
ξ
=
β
σ
1
-
3
β
n
0
ξ
+
O
β
v
ξ
ξ
,and system (35) can be approximated to order O
(
α
,
β
) by(38)
I
-
β
2
∂
ξ
2
v
t
=
2
β
3
∂
ξ
4
Φ
-
∂
ξ
1
+
α
v
M
2
Φ
ξ
,
(39)
I
-
β
2
∂
ξ
2
Φ
t
=
β
σ
v
ξ
ξ
-
3
M
′
M
4
v
ξ
+
3
M
′
2
M
5
v
-
M
′
′
M
4
v
-
1
M
v
-
α
2
M
2
Φ
ξ
2
.Note that the coefficient M
(
ξ
) is smooth even when the function describing the bottom y
=
-
H
(
x
/
γ
) is discontinuous or nondifferentiable. The function M
(
ξ
) is time independent and becomes identically one for a channel with constant depth. In this case, system (38)-(39) reduces to that studied in [4]. Moreover, in applications, this coefficient is bounded and i
n
f
R
M
(
ξ
)
>
0. These properties will be important to obtain existence and uniqueness results for system (38)-(39). We also point out that the function M
(
ξ
) actually depends on the dispersion parameter β. For this reason, it will be denoted by M
(
ξ
,
β
) whenever we need to emphasize this dependence.


## 4. Existence and Uniqueness

System (38)-(39) can be written as(40)
v
Φ
t
=
A
v
Φ
+
G
v
Φ
,where(41)
A
v
Φ
≔
0
I
-
β
2
∂
ξ
2
-
1
I
-
2
β
3
∂
ξ
2
-
∂
ξ
2
I
-
β
2
∂
ξ
2
-
1
β
σ
∂
ξ
2
0
v
Φ
,
G
v
Φ
≔
-
I
-
β
2
∂
ξ
2
-
1
α
M
2
v
Φ
ξ
ξ
I
-
β
2
∂
ξ
2
-
1
-
3
β
σ
M
′
M
4
v
ξ
+
3
β
σ
M
′
2
M
5
v
-
β
σ
M
′
′
M
4
v
-
1
M
v
-
α
2
M
2
Φ
ξ
2
.System (40) is supplemented with the initial conditions(42)
v
ξ
,
0
=
v
0
ξ
,
Φ
ξ
,
0
=
Φ
0
ξ
.Taking Fourier transform in the spatial variable ξ in system (40), we have(43)
U
^
t
=
A
y
U
^
+
G
U
^
,
U
^
y
,
0
=
U
^
0
y
,where(44)
U
≔
v
Φ
,
U
0
=
v
0
Φ
0
,
A
y
=
0
y
2
1
+
2
β
/
3
y
2
1
+
β
/
2
y
2
-
β
σ
y
2
1
+
β
/
2
y
2
0
.


### 4.1. Analysis of the Linear Semigroup

Consider system (40) with G
≡
0,(45)
U
^
t
=
A
y
U
^
,
U
^
y
,
0
=
U
^
0
y
,which has unique solution in the form(46)
U
^
y
,
t
=
e
A
y
t
U
^
0
y
,where(47)
e
A
y
t
=
a
11
y
,
t
a
12
y
,
t
a
21
y
,
t
a
22
y
,
t
=
cos
⁡
λ
y
t
Γ
y
sin
⁡
λ
y
t
Ψ
y
sin
⁡
λ
y
t
cos
⁡
λ
y
t
,
λ
y
=
2
y
2
β
σ
3
+
2
β
y
2
3
2
+
β
y
2
,
Γ
y
=
3
+
2
β
y
2
3
β
σ
,
Ψ
y
=
-
3
β
σ
3
+
2
β
y
2
.By using the inverse Fourier transform in (46), we have that(48)
U
x
,
t
=
S
t
U
0
x
,where(49)
S
t
f
=
F
-
1
e
A
·
t
f
^
·
,with(50)
f
=
f
1
f
2
.Theorem 1.
The family of linear operators(
S
(
t
)
)
t
≥
0 is a C
0-semigroup in H
s
-
2
×
H
s
-
1. Furthermore A
:
H
s
-
1
×
H
s
→
H
s
-
2
×
H
s
-
1 is its infinitesimal generator in H
s
-
2
×
H
s
-
1.Proof.
Letf
=
(
f
1
,
f
2
)
T. Then(51)
S
t
f
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
f
^
1
+
a
12
f
^
2
F
-
1
a
21
f
^
1
+
a
22
f
^
2
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
f
^
1
+
F
-
1
a
12
f
^
2
H
s
-
2
2
+
F
-
1
a
21
f
^
1
+
F
-
1
a
22
f
^
2
H
s
-
1
2
≤
C
∫
R
1
+
y
2
s
-
2
a
11
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
a
12
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
21
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
22
2
f
^
2
2
d
y
≤
C
∫
R
1
+
y
2
s
-
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
3
+
2
β
y
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
1
3
+
2
β
y
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
f
^
2
2
d
y
≤
C
f
1
H
s
-
2
2
+
f
2
H
s
-
1
2
.Hereafter C will denote a generic constant and we recall that 0
<
β
<
1. Therefore S
(
t
), t
≥
0, is a family of continuous linear operators in H
s
-
2
×
H
s
-
1. On the other hand, it is easy to see that S
(
0
)
=
I, S
(
t
+
s
)
=
S
(
t
)
S
(
s
)
=
S
(
s
)
S
(
t
), t
,
s
≥
0. Finally (52)
S
t
f
-
f
H
s
-
2
×
H
s
-
1
2
≤
F
-
1
a
11
f
^
1
+
a
12
f
^
2
-
f
^
1
F
-
1
a
21
f
^
1
+
a
22
f
^
2
-
f
^
2
H
s
-
2
×
H
s
-
1
2
≤
C
∫
R
1
+
y
2
s
-
2
a
11
-
1
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
a
12
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
21
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
22
-
1
2
f
^
2
2
d
y
.Observe that a
12
(
y
,
t
)
→
0, a
21
(
y
,
t
)
→
0 as t
→
0
+ and a
11
(
y
,
t
)
→
1, a
22
(
y
,
t
)
→
1, as t
→
0
+. By using Lebesgue’s dominated convergence theorem, we can conclude that(53)
S
t
f
-
f
H
s
-
2
×
H
s
-
1
⟶
0
,as t
→
0
+. We conclude that S
(
t
), t
≥
0, is a strongly continuous semigroup in H
s
-
2
×
H
s
-
1.
On the other hand,(54)
A
f
H
s
-
2
×
H
s
-
1
2
≤
∫
R
1
+
y
2
s
-
2
F
I
-
β
2
∂
ξ
2
-
1
I
-
2
β
3
∂
ξ
2
∂
ξ
2
f
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
F
I
-
β
2
∂
ξ
2
-
1
β
σ
∂
ξ
2
f
1
2
d
y
≤
∫
R
1
+
y
2
s
-
2
1
+
2
β
/
3
y
2
2
y
4
1
+
β
/
2
y
2
2
f
^
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
β
2
σ
2
y
4
1
+
β
/
2
y
2
2
f
^
1
2
d
y
≤
C
f
1
s
-
1
2
+
f
2
s
2
=
C
f
H
s
-
1
×
H
s
2
.To see that A is the infinitesimal generator of the semigroup S
(
t
), t
≥
0, observe that, for f
∈
H
s
-
1
×
H
s, (55)
S
h
f
-
f
h
-
A
f
H
s
-
2
×
H
s
-
1
2
=
F
-
1
e
A
h
-
I
h
f
^
-
A
f
^
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
·
,
h
f
^
1
+
a
12
·
,
h
f
^
2
-
f
^
1
h
-
I
-
β
2
∂
ξ
2
-
1
I
-
2
β
3
∂
ξ
2
-
∂
ξ
2
f
2
H
s
-
2
2
+
F
-
1
a
21
·
,
h
f
^
1
+
a
22
·
,
h
f
^
2
-
f
^
2
h
-
I
-
β
2
∂
ξ
2
-
1
β
σ
∂
ξ
2
f
1
H
s
-
1
2
=
∫
R
1
+
y
2
s
-
2
cos
⁡
λ
h
-
1
f
^
1
+
Γ
sin
⁡
λ
h
f
^
2
h
-
y
2
1
+
2
β
/
3
y
2
f
^
2
1
+
β
/
2
y
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
Ψ
sin
⁡
λ
h
f
^
1
+
cos
⁡
λ
h
-
1
f
^
2
h
-
-
β
σ
y
2
f
^
1
1
+
β
/
2
y
2
2
d
y
.But by virtue of (56)
l
i
m
h
→
0
cos
⁡
λ
h
-
1
h
=
0
,
l
i
m
h
→
0
Γ
s
i
n
λ
h
h
=
Γ
λ
=
y
2
1
+
2
β
/
3
y
2
1
+
β
/
2
y
2
,
l
i
m
h
→
0
Ψ
s
i
n
λ
h
h
=
Ψ
λ
=
-
β
σ
y
2
1
+
β
/
2
y
2and using again Lebesgue’s dominated convergence theorem, we get that(57)
l
i
m
h
→
0
S
h
f
-
f
h
-
A
f
H
s
-
2
×
H
s
-
1
2
=
0
.


### 4.2. Analysis of the Nonlinear Term

Theorem 2.
Lets
>
3
/
2. The application G maps H
s
-
1
×
H
s on itself and(58)
G
U
H
s
-
1
×
H
s
≤
C
U
H
s
-
1
×
H
s
.Furthermore,(59)
G
U
1
-
G
U
2
H
s
-
1
×
H
s
2
≤
C
1
+
Φ
1
H
s
2
v
1
-
v
2
H
s
-
1
2
+
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
,where C
>
0 is a constant and U
i
=
v
i
Φ
i, i
=
1,2.Proof.
In the first place, letU
=
(
v
,
Φ
)
T: (60)
G
U
H
s
-
1
×
H
s
2
≤
C
∫
R
1
+
y
2
s
-
1
1
1
+
β
/
2
y
2
2
v
M
2
Φ
ξ
ξ
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
M
4
v
ξ
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
2
M
5
v
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
′
M
4
v
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
v
M
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
Φ
ξ
2
M
2
^
2
d
y
.Taking into account the inequalities (61)
v
M
2
Φ
ξ
H
s
-
2
≤
C
v
M
2
H
s
-
2
Φ
ξ
H
s
-
1
≤
C
1
M
2
H
s
-
1
v
H
s
-
1
Φ
H
s
,
M
′
M
4
v
ξ
H
s
-
2
≤
C
M
′
M
4
H
s
-
1
v
ξ
H
s
-
2
≤
C
M
′
M
4
H
s
-
1
v
H
s
-
1
,
M
′
2
M
5
v
H
s
-
2
≤
C
M
′
2
M
5
H
s
-
1
v
H
s
-
2
≤
C
M
′
2
M
5
H
s
-
1
v
H
s
-
1
,
M
′
′
M
4
v
H
s
-
2
≤
C
M
′
′
M
4
H
s
-
1
v
H
s
-
1
,
v
M
H
s
-
2
≤
C
1
M
H
s
-
1
v
H
s
-
2
≤
C
1
M
H
s
-
1
v
H
s
-
1
,
Φ
ξ
2
M
2
H
s
-
2
≤
C
Φ
ξ
M
2
Φ
ξ
H
s
-
2
≤
C
Φ
ξ
H
s
-
1
Φ
ξ
M
2
H
s
-
2
≤
C
Φ
s
1
M
2
H
s
-
1
Φ
ξ
H
s
-
2
≤
C
Φ
s
1
M
2
H
s
-
1
Φ
H
s
, valid for s
>
3
/
2 (applying Corollary  3.16 in [31]), and using the fact that the function M is bounded, we get that (62)
G
U
H
s
-
1
×
H
s
2
≤
C
v
M
2
Φ
ξ
H
s
-
2
2
+
M
′
M
4
v
ξ
H
s
-
2
2
+
M
′
2
M
5
v
H
s
-
2
2
+
M
′
′
M
4
v
H
s
-
2
2
+
v
M
H
s
-
2
2
+
Φ
ξ
2
M
2
H
s
-
2
2
≤
C
v
H
s
-
1
2
+
Φ
H
s
2
.On the other hand, for U
i
=
(
v
i
,
Φ
i
)
T
∈
H
s
-
1
×
H
s, we have that (63)
G
v
1
Φ
1
-
G
v
2
Φ
2
H
s
-
1
×
H
s
2
≤
C
v
1
M
2
Φ
1
ξ
-
v
2
M
2
Φ
2
ξ
H
s
-
2
2
+
M
′
M
4
v
1
-
v
2
ξ
H
s
-
2
2
+
M
′
2
M
5
v
1
-
v
2
H
s
-
2
2
+
M
′
′
M
4
v
1
-
v
2
H
s
-
2
2
+
1
M
v
1
-
v
2
H
s
-
2
2
+
1
M
2
Φ
1
ξ
2
-
Φ
2
ξ
2
H
s
-
2
2
≤
C
v
1
Φ
1
ξ
-
v
2
Φ
2
ξ
H
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
ξ
2
-
Φ
2
ξ
2
H
s
-
2
2
≤
C
v
1
-
v
2
Φ
1
ξ
H
s
-
2
+
v
2
Φ
1
ξ
-
Φ
2
ξ
H
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
ξ
+
Φ
2
ξ
H
s
-
1
2
Φ
1
ξ
-
Φ
2
ξ
H
s
-
1
2
≤
C
v
1
-
v
2
H
s
-
1
2
Φ
1
ξ
H
s
-
2
2
+
v
2
H
s
-
1
2
Φ
1
ξ
-
Φ
2
ξ
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
≤
C
1
+
Φ
1
H
s
2
v
1
-
v
2
H
s
-
1
2
+
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
.


### 4.3. Local Existence and Uniqueness

Theorem 3.
Lets
>
3
/
2 and ϕ
=
(
v
0
,
Φ
0
)
T
∈
H
s
-
1
×
H
s. Then there exists T
(
s
,
ϕ
)
>
0 and unique U
=
(
v
,
Φ
)
T
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
), which satisfies the integral equation(64)
U
t
=
S
t
ϕ
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
.Proof.
LetM
,
T
>
0 be fixed constants and consider the nonlinear operator(65)
Ψ
U
t
=
S
t
ϕ
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
,defined in the complete metric space(66)
X
M
T
=
U
∈
C
0
,
T
;
H
s
-
1
×
H
s
:
s
u
p
t
∈
0
,
T
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
≤
M
.Let us prove that if U
∈
X
M
(
T
), then Ψ
U
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
). Indeed if U
∈
X
M
(
T
), then (67)
Ψ
U
t
-
Ψ
U
τ
H
s
-
1
×
H
s
≤
S
t
-
S
τ
ϕ
H
s
-
1
×
H
s
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
-
S
τ
-
t
′
G
U
t
′
d
t
′
H
s
-
1
×
H
s
.We remark that if U
=
(
v
,
Φ
)
T
∈
X
M
(
T
), then(68)
U
t
H
s
-
1
×
H
s
≤
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
+
S
t
ϕ
H
s
-
1
×
H
s
≤
M
+
ϕ
H
s
-
1
×
H
s
.Given that (
S
(
t
)
)
t
≥
0 is a C
0-semigroup, then (
S
(
t
)
-
S
(
τ
)
)
ϕ
H
s
-
1
×
H
s
→
0 as τ
→
t. To analyze the second expression, let us suppose that τ
>
t
>
0. We obtain that (69)
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
-
∫
0
τ
S
τ
-
t
′
G
U
t
′
d
t
′
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
-
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
+
∫
t
τ
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≔
I
1
+
I
2
.Observe that since a
i
j
(
y
,
τ
-
t
′
)
→
a
i
j
(
y
,
t
-
t
′
) as τ
→
t,  i
,
j
=
1,2, and using Lebesgue’s dominated convergence theorem, we get that I
1
→
0 as τ
→
t.
On the other hand, due to Theorems1 and 2, we arrive at (70)
I
2
=
∫
t
τ
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
t
τ
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
t
τ
U
H
s
-
1
×
H
s
d
t
′
≤
C
τ
-
t
M
+
ϕ
H
s
-
1
×
H
s
.Therefore we get that I
2
→
0 as τ
→
t. This means that Ψ
U
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
).
Let us prove existence of small timeT
>
0 such that Ψ
(
X
M
(
T
)
)
⊂
X
M
(
T
). Let U
=
(
v
,
Φ
)
T
∈
X
M
(
T
). For any t
∈
[
0
,
T
], (71)
Ψ
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
G
U
t
′
H
s
-
1
×
H
s
≤
C
∫
0
t
G
U
t
′
H
s
-
1
×
H
s
≤
C
T
s
u
p
t
∈
0
,
T
U
t
H
s
-
1
×
H
s
≤
C
T
M
+
ϕ
H
s
-
1
×
H
s
≤
M
,for T being small enough. Therefore Ψ
(
U
)
∈
X
M
(
T
).
Finally, let us show that timeT
^
∈
(
0
,
T
] exists so that the operator Ψ is a contraction on X
M
(
T
^
). Let U
1
=
(
v
1
,
Φ
1
)
T, U
2
=
(
v
2
,
Φ
2
)
T
∈
X
M
(
T
). Then using Theorem 2, we arrive at (72)
Ψ
U
1
t
-
Ψ
U
2
t
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
G
U
1
t
′
-
G
U
2
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
0
t
G
U
1
-
G
U
2
H
s
-
1
×
H
s
d
t
′
≤
∫
0
t
1
+
Φ
1
H
s
v
1
-
v
2
H
s
-
1
+
v
2
H
s
-
1
+
Φ
1
+
Φ
2
H
s
Φ
1
-
Φ
2
H
s
d
t
′
≤
C
T
1
+
M
+
ϕ
H
s
-
1
×
H
s
v
1
-
v
2
H
s
-
1
+
3
M
+
3
ϕ
H
s
-
1
×
H
s
Φ
1
-
Φ
2
H
s
≤
C
T
3
M
+
3
ϕ
H
s
-
1
×
H
s
s
u
p
t
∈
0
,
T
U
1
-
U
2
H
s
-
1
×
H
s
.We see that it is possible to select time T
^
∈
(
0
,
T
] such that(73)
C
T
^
3
M
+
3
ϕ
H
s
-
1
×
H
s
<
M
,so that the operator Ψ is a contraction on the closed ball X
M
(
T
^
)
⊂
C
(
[
0
,
T
^
]
;
H
s
-
1
×
H
s
). The fixed point principle guarantees the existence of unique solution (74)
U
∈
C
0
,
T
^
;
H
s
-
1
×
H
sof integral equation (64). Finally, uniqueness of this solution in the metric space C
(
[
0
,
T
^
]
;
H
s
-
1
×
H
s
) can be established via Gronwall’s lemma.


## 4.1. Analysis of the Linear Semigroup

Consider system (40) with G
≡
0,(45)
U
^
t
=
A
y
U
^
,
U
^
y
,
0
=
U
^
0
y
,which has unique solution in the form(46)
U
^
y
,
t
=
e
A
y
t
U
^
0
y
,where(47)
e
A
y
t
=
a
11
y
,
t
a
12
y
,
t
a
21
y
,
t
a
22
y
,
t
=
cos
⁡
λ
y
t
Γ
y
sin
⁡
λ
y
t
Ψ
y
sin
⁡
λ
y
t
cos
⁡
λ
y
t
,
λ
y
=
2
y
2
β
σ
3
+
2
β
y
2
3
2
+
β
y
2
,
Γ
y
=
3
+
2
β
y
2
3
β
σ
,
Ψ
y
=
-
3
β
σ
3
+
2
β
y
2
.By using the inverse Fourier transform in (46), we have that(48)
U
x
,
t
=
S
t
U
0
x
,where(49)
S
t
f
=
F
-
1
e
A
·
t
f
^
·
,with(50)
f
=
f
1
f
2
.Theorem 1.
The family of linear operators(
S
(
t
)
)
t
≥
0 is a C
0-semigroup in H
s
-
2
×
H
s
-
1. Furthermore A
:
H
s
-
1
×
H
s
→
H
s
-
2
×
H
s
-
1 is its infinitesimal generator in H
s
-
2
×
H
s
-
1.Proof.
Letf
=
(
f
1
,
f
2
)
T. Then(51)
S
t
f
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
f
^
1
+
a
12
f
^
2
F
-
1
a
21
f
^
1
+
a
22
f
^
2
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
f
^
1
+
F
-
1
a
12
f
^
2
H
s
-
2
2
+
F
-
1
a
21
f
^
1
+
F
-
1
a
22
f
^
2
H
s
-
1
2
≤
C
∫
R
1
+
y
2
s
-
2
a
11
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
a
12
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
21
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
22
2
f
^
2
2
d
y
≤
C
∫
R
1
+
y
2
s
-
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
3
+
2
β
y
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
1
3
+
2
β
y
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
f
^
2
2
d
y
≤
C
f
1
H
s
-
2
2
+
f
2
H
s
-
1
2
.Hereafter C will denote a generic constant and we recall that 0
<
β
<
1. Therefore S
(
t
), t
≥
0, is a family of continuous linear operators in H
s
-
2
×
H
s
-
1. On the other hand, it is easy to see that S
(
0
)
=
I, S
(
t
+
s
)
=
S
(
t
)
S
(
s
)
=
S
(
s
)
S
(
t
), t
,
s
≥
0. Finally (52)
S
t
f
-
f
H
s
-
2
×
H
s
-
1
2
≤
F
-
1
a
11
f
^
1
+
a
12
f
^
2
-
f
^
1
F
-
1
a
21
f
^
1
+
a
22
f
^
2
-
f
^
2
H
s
-
2
×
H
s
-
1
2
≤
C
∫
R
1
+
y
2
s
-
2
a
11
-
1
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
2
a
12
2
f
^
2
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
21
2
f
^
1
2
d
y
+
C
∫
R
1
+
y
2
s
-
1
a
22
-
1
2
f
^
2
2
d
y
.Observe that a
12
(
y
,
t
)
→
0, a
21
(
y
,
t
)
→
0 as t
→
0
+ and a
11
(
y
,
t
)
→
1, a
22
(
y
,
t
)
→
1, as t
→
0
+. By using Lebesgue’s dominated convergence theorem, we can conclude that(53)
S
t
f
-
f
H
s
-
2
×
H
s
-
1
⟶
0
,as t
→
0
+. We conclude that S
(
t
), t
≥
0, is a strongly continuous semigroup in H
s
-
2
×
H
s
-
1.
On the other hand,(54)
A
f
H
s
-
2
×
H
s
-
1
2
≤
∫
R
1
+
y
2
s
-
2
F
I
-
β
2
∂
ξ
2
-
1
I
-
2
β
3
∂
ξ
2
∂
ξ
2
f
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
F
I
-
β
2
∂
ξ
2
-
1
β
σ
∂
ξ
2
f
1
2
d
y
≤
∫
R
1
+
y
2
s
-
2
1
+
2
β
/
3
y
2
2
y
4
1
+
β
/
2
y
2
2
f
^
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
β
2
σ
2
y
4
1
+
β
/
2
y
2
2
f
^
1
2
d
y
≤
C
f
1
s
-
1
2
+
f
2
s
2
=
C
f
H
s
-
1
×
H
s
2
.To see that A is the infinitesimal generator of the semigroup S
(
t
), t
≥
0, observe that, for f
∈
H
s
-
1
×
H
s, (55)
S
h
f
-
f
h
-
A
f
H
s
-
2
×
H
s
-
1
2
=
F
-
1
e
A
h
-
I
h
f
^
-
A
f
^
H
s
-
2
×
H
s
-
1
2
=
F
-
1
a
11
·
,
h
f
^
1
+
a
12
·
,
h
f
^
2
-
f
^
1
h
-
I
-
β
2
∂
ξ
2
-
1
I
-
2
β
3
∂
ξ
2
-
∂
ξ
2
f
2
H
s
-
2
2
+
F
-
1
a
21
·
,
h
f
^
1
+
a
22
·
,
h
f
^
2
-
f
^
2
h
-
I
-
β
2
∂
ξ
2
-
1
β
σ
∂
ξ
2
f
1
H
s
-
1
2
=
∫
R
1
+
y
2
s
-
2
cos
⁡
λ
h
-
1
f
^
1
+
Γ
sin
⁡
λ
h
f
^
2
h
-
y
2
1
+
2
β
/
3
y
2
f
^
2
1
+
β
/
2
y
2
2
d
y
+
∫
R
1
+
y
2
s
-
1
Ψ
sin
⁡
λ
h
f
^
1
+
cos
⁡
λ
h
-
1
f
^
2
h
-
-
β
σ
y
2
f
^
1
1
+
β
/
2
y
2
2
d
y
.But by virtue of (56)
l
i
m
h
→
0
cos
⁡
λ
h
-
1
h
=
0
,
l
i
m
h
→
0
Γ
s
i
n
λ
h
h
=
Γ
λ
=
y
2
1
+
2
β
/
3
y
2
1
+
β
/
2
y
2
,
l
i
m
h
→
0
Ψ
s
i
n
λ
h
h
=
Ψ
λ
=
-
β
σ
y
2
1
+
β
/
2
y
2and using again Lebesgue’s dominated convergence theorem, we get that(57)
l
i
m
h
→
0
S
h
f
-
f
h
-
A
f
H
s
-
2
×
H
s
-
1
2
=
0
.


## 4.2. Analysis of the Nonlinear Term

Theorem 2.
Lets
>
3
/
2. The application G maps H
s
-
1
×
H
s on itself and(58)
G
U
H
s
-
1
×
H
s
≤
C
U
H
s
-
1
×
H
s
.Furthermore,(59)
G
U
1
-
G
U
2
H
s
-
1
×
H
s
2
≤
C
1
+
Φ
1
H
s
2
v
1
-
v
2
H
s
-
1
2
+
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
,where C
>
0 is a constant and U
i
=
v
i
Φ
i, i
=
1,2.Proof.
In the first place, letU
=
(
v
,
Φ
)
T: (60)
G
U
H
s
-
1
×
H
s
2
≤
C
∫
R
1
+
y
2
s
-
1
1
1
+
β
/
2
y
2
2
v
M
2
Φ
ξ
ξ
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
M
4
v
ξ
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
2
M
5
v
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
M
′
′
M
4
v
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
v
M
^
2
d
y
+
C
∫
R
1
+
y
2
s
1
1
+
β
/
2
y
2
2
Φ
ξ
2
M
2
^
2
d
y
.Taking into account the inequalities (61)
v
M
2
Φ
ξ
H
s
-
2
≤
C
v
M
2
H
s
-
2
Φ
ξ
H
s
-
1
≤
C
1
M
2
H
s
-
1
v
H
s
-
1
Φ
H
s
,
M
′
M
4
v
ξ
H
s
-
2
≤
C
M
′
M
4
H
s
-
1
v
ξ
H
s
-
2
≤
C
M
′
M
4
H
s
-
1
v
H
s
-
1
,
M
′
2
M
5
v
H
s
-
2
≤
C
M
′
2
M
5
H
s
-
1
v
H
s
-
2
≤
C
M
′
2
M
5
H
s
-
1
v
H
s
-
1
,
M
′
′
M
4
v
H
s
-
2
≤
C
M
′
′
M
4
H
s
-
1
v
H
s
-
1
,
v
M
H
s
-
2
≤
C
1
M
H
s
-
1
v
H
s
-
2
≤
C
1
M
H
s
-
1
v
H
s
-
1
,
Φ
ξ
2
M
2
H
s
-
2
≤
C
Φ
ξ
M
2
Φ
ξ
H
s
-
2
≤
C
Φ
ξ
H
s
-
1
Φ
ξ
M
2
H
s
-
2
≤
C
Φ
s
1
M
2
H
s
-
1
Φ
ξ
H
s
-
2
≤
C
Φ
s
1
M
2
H
s
-
1
Φ
H
s
, valid for s
>
3
/
2 (applying Corollary  3.16 in [31]), and using the fact that the function M is bounded, we get that (62)
G
U
H
s
-
1
×
H
s
2
≤
C
v
M
2
Φ
ξ
H
s
-
2
2
+
M
′
M
4
v
ξ
H
s
-
2
2
+
M
′
2
M
5
v
H
s
-
2
2
+
M
′
′
M
4
v
H
s
-
2
2
+
v
M
H
s
-
2
2
+
Φ
ξ
2
M
2
H
s
-
2
2
≤
C
v
H
s
-
1
2
+
Φ
H
s
2
.On the other hand, for U
i
=
(
v
i
,
Φ
i
)
T
∈
H
s
-
1
×
H
s, we have that (63)
G
v
1
Φ
1
-
G
v
2
Φ
2
H
s
-
1
×
H
s
2
≤
C
v
1
M
2
Φ
1
ξ
-
v
2
M
2
Φ
2
ξ
H
s
-
2
2
+
M
′
M
4
v
1
-
v
2
ξ
H
s
-
2
2
+
M
′
2
M
5
v
1
-
v
2
H
s
-
2
2
+
M
′
′
M
4
v
1
-
v
2
H
s
-
2
2
+
1
M
v
1
-
v
2
H
s
-
2
2
+
1
M
2
Φ
1
ξ
2
-
Φ
2
ξ
2
H
s
-
2
2
≤
C
v
1
Φ
1
ξ
-
v
2
Φ
2
ξ
H
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
ξ
2
-
Φ
2
ξ
2
H
s
-
2
2
≤
C
v
1
-
v
2
Φ
1
ξ
H
s
-
2
+
v
2
Φ
1
ξ
-
Φ
2
ξ
H
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
ξ
+
Φ
2
ξ
H
s
-
1
2
Φ
1
ξ
-
Φ
2
ξ
H
s
-
1
2
≤
C
v
1
-
v
2
H
s
-
1
2
Φ
1
ξ
H
s
-
2
2
+
v
2
H
s
-
1
2
Φ
1
ξ
-
Φ
2
ξ
s
-
2
2
+
v
1
-
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
≤
C
1
+
Φ
1
H
s
2
v
1
-
v
2
H
s
-
1
2
+
v
2
H
s
-
1
2
+
Φ
1
+
Φ
2
H
s
2
Φ
1
-
Φ
2
H
s
2
.


## 4.3. Local Existence and Uniqueness

Theorem 3.
Lets
>
3
/
2 and ϕ
=
(
v
0
,
Φ
0
)
T
∈
H
s
-
1
×
H
s. Then there exists T
(
s
,
ϕ
)
>
0 and unique U
=
(
v
,
Φ
)
T
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
), which satisfies the integral equation(64)
U
t
=
S
t
ϕ
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
.Proof.
LetM
,
T
>
0 be fixed constants and consider the nonlinear operator(65)
Ψ
U
t
=
S
t
ϕ
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
,defined in the complete metric space(66)
X
M
T
=
U
∈
C
0
,
T
;
H
s
-
1
×
H
s
:
s
u
p
t
∈
0
,
T
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
≤
M
.Let us prove that if U
∈
X
M
(
T
), then Ψ
U
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
). Indeed if U
∈
X
M
(
T
), then (67)
Ψ
U
t
-
Ψ
U
τ
H
s
-
1
×
H
s
≤
S
t
-
S
τ
ϕ
H
s
-
1
×
H
s
+
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
-
S
τ
-
t
′
G
U
t
′
d
t
′
H
s
-
1
×
H
s
.We remark that if U
=
(
v
,
Φ
)
T
∈
X
M
(
T
), then(68)
U
t
H
s
-
1
×
H
s
≤
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
+
S
t
ϕ
H
s
-
1
×
H
s
≤
M
+
ϕ
H
s
-
1
×
H
s
.Given that (
S
(
t
)
)
t
≥
0 is a C
0-semigroup, then (
S
(
t
)
-
S
(
τ
)
)
ϕ
H
s
-
1
×
H
s
→
0 as τ
→
t. To analyze the second expression, let us suppose that τ
>
t
>
0. We obtain that (69)
∫
0
t
S
t
-
t
′
G
U
t
′
d
t
′
-
∫
0
τ
S
τ
-
t
′
G
U
t
′
d
t
′
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
-
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
+
∫
t
τ
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≔
I
1
+
I
2
.Observe that since a
i
j
(
y
,
τ
-
t
′
)
→
a
i
j
(
y
,
t
-
t
′
) as τ
→
t,  i
,
j
=
1,2, and using Lebesgue’s dominated convergence theorem, we get that I
1
→
0 as τ
→
t.
On the other hand, due to Theorems1 and 2, we arrive at (70)
I
2
=
∫
t
τ
S
τ
-
t
′
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
t
τ
G
U
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
t
τ
U
H
s
-
1
×
H
s
d
t
′
≤
C
τ
-
t
M
+
ϕ
H
s
-
1
×
H
s
.Therefore we get that I
2
→
0 as τ
→
t. This means that Ψ
U
∈
C
(
[
0
,
T
]
;
H
s
-
1
×
H
s
).
Let us prove existence of small timeT
>
0 such that Ψ
(
X
M
(
T
)
)
⊂
X
M
(
T
). Let U
=
(
v
,
Φ
)
T
∈
X
M
(
T
). For any t
∈
[
0
,
T
], (71)
Ψ
U
t
-
S
t
ϕ
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
G
U
t
′
H
s
-
1
×
H
s
≤
C
∫
0
t
G
U
t
′
H
s
-
1
×
H
s
≤
C
T
s
u
p
t
∈
0
,
T
U
t
H
s
-
1
×
H
s
≤
C
T
M
+
ϕ
H
s
-
1
×
H
s
≤
M
,for T being small enough. Therefore Ψ
(
U
)
∈
X
M
(
T
).
Finally, let us show that timeT
^
∈
(
0
,
T
] exists so that the operator Ψ is a contraction on X
M
(
T
^
). Let U
1
=
(
v
1
,
Φ
1
)
T, U
2
=
(
v
2
,
Φ
2
)
T
∈
X
M
(
T
). Then using Theorem 2, we arrive at (72)
Ψ
U
1
t
-
Ψ
U
2
t
H
s
-
1
×
H
s
≤
∫
0
t
S
t
-
t
′
G
U
1
t
′
-
G
U
2
t
′
H
s
-
1
×
H
s
d
t
′
≤
C
∫
0
t
G
U
1
-
G
U
2
H
s
-
1
×
H
s
d
t
′
≤
∫
0
t
1
+
Φ
1
H
s
v
1
-
v
2
H
s
-
1
+
v
2
H
s
-
1
+
Φ
1
+
Φ
2
H
s
Φ
1
-
Φ
2
H
s
d
t
′
≤
C
T
1
+
M
+
ϕ
H
s
-
1
×
H
s
v
1
-
v
2
H
s
-
1
+
3
M
+
3
ϕ
H
s
-
1
×
H
s
Φ
1
-
Φ
2
H
s
≤
C
T
3
M
+
3
ϕ
H
s
-
1
×
H
s
s
u
p
t
∈
0
,
T
U
1
-
U
2
H
s
-
1
×
H
s
.We see that it is possible to select time T
^
∈
(
0
,
T
] such that(73)
C
T
^
3
M
+
3
ϕ
H
s
-
1
×
H
s
<
M
,so that the operator Ψ is a contraction on the closed ball X
M
(
T
^
)
⊂
C
(
[
0
,
T
^
]
;
H
s
-
1
×
H
s
). The fixed point principle guarantees the existence of unique solution (74)
U
∈
C
0
,
T
^
;
H
s
-
1
×
H
sof integral equation (64). Finally, uniqueness of this solution in the metric space C
(
[
0
,
T
^
]
;
H
s
-
1
×
H
s
) can be established via Gronwall’s lemma.


## 5. Numerical Schemes

In this section we describe the numerical scheme we propose to approximate solutions of system (38)-(39). In the first place, it is convenient to rewrite it as(75)
I
-
β
2
∂
ξ
2
v
t
=
2
β
3
∂
ξ
3
u
-
∂
ξ
1
+
α
v
M
2
u
,
I
-
β
2
∂
ξ
2
u
t
=
β
σ
∂
ξ
3
v
+
β
σ
-
3
M
′
M
4
v
ξ
+
3
M
′
2
M
5
v
-
M
′
′
M
4
v
-
∂
ξ
1
M
v
-
α
2
∂
ξ
1
M
2
u
2
,subject to the initial conditions v
(
ξ
,
0
)
=
v
0
(
ξ
), u
(
ξ
,
0
)
=
u
0
(
ξ
), where we introduced the new variable u
≔
Φ
ξ. In the numerical solver to be introduced, the computational domain [
0
,
L
] is discretized by N equidistant points, with spacing Δ
ξ
=
L
/
N, and the unknowns v and u are expanded as truncated Fourier series in space with time-dependent coefficients:(76)
v
ξ
,
t
=
∑
j
v
^
j
t
e
i
w
j
ξ
u
ξ
,
t
=
∑
j
u
^
j
t
e
i
w
j
ξ
,with (77)
w
j
=
2
π
j
L
,
j
=
-
N
2
+
1
,
…
,
0
,
…
,
N
2
.The time-dependent coefficients v
^
j
(
t
), j
=
-
N
/
2
+
1
,
…
,
0
,
…
,
N
/
2, are calculated by means of the equation(78)
v
^
j
t
=
1
L
∫
0
L
v
ξ
,
t
e
-
i
w
j
ξ
d
ξand analogously for u
^
j
(
t
). Substituting these expressions into (75) and projecting the resulting equations with respect to the L
2-orthonormal basis ϕ
j
=
L
-
1
/
2
e
i
w
j
ξ and the inner product(79)
f
,
g
=
∫
0
L
f
ξ
g
ξ
¯
d
ξ
,it follows that(80)
v
^
j
′
t
=
γ
1
u
^
j
+
N
1
,
u
^
j
′
t
=
γ
2
v
^
j
+
N
2
,where(81)
γ
1
=
-
2
β
/
3
i
w
j
3
1
+
β
/
2
w
j
2
,
γ
2
=
-
β
σ
i
w
j
3
1
+
β
/
2
w
j
2
,
N
1
=
-
i
w
j
P
j
1
+
α
v
/
M
2
u
1
+
β
/
2
w
j
2
,
N
2
=
i
w
j
P
j
β
σ
-
3
M
′
/
M
4
v
ξ
+
3
M
′
2
/
M
5
v
-
M
′
′
/
M
4
v
-
1
/
M
v
-
α
/
2
M
2
u
2
1
+
β
/
2
w
j
2
,and P
j
[
·
] denotes the operator(82)
P
j
g
=
1
L
∫
0
L
g
ξ
e
-
i
w
j
ξ
d
ξ
.When the period L is taken large enough, this numerical scheme can be applied to approximate the rapidly decaying solutions of system (75) on the entire real line R. This technique was used successfully in [9].


### 5.1. Temporal Discretization

Note that (80) can be seen as a system of ordinary differential equations where the unknowns are the time-dependent Fourier coefficients of the solutions. To solve it, we use an implicit-explicit scheme (IMEX) in the form(83)
v
^
j
n
+
1
-
v
^
j
n
Δ
t
=
γ
1
u
^
j
n
+
1
+
u
^
j
n
2
+
3
2
N
1
n
-
1
2
N
1
n
-
1
,
u
^
j
n
+
1
-
u
^
j
n
Δ
t
=
γ
2
v
^
j
n
+
1
+
v
^
j
n
2
+
3
2
N
2
n
-
1
2
N
2
n
-
1
.Here v
n, u
n denote the approximations of the unknowns v
(
ξ
,
t
), u
(
ξ
,
t
), respectively, at time t
=
n
Δ
t, where Δ
t is the time step of the method, and N
1
n
,
N
2
n denote the approximations of the functions N
1
,
N
2 evaluated at time n
Δ
t. Similarly, v
^
j
n
,
u
^
j
n denote the approximations to the Fourier transforms of the functions v and u, respectively, with respect to the variable ξ, evaluated at time n
Δ
t. The numerical approach adopted for solving system (75) ensures the scheme results to be linearly unconditionally stable which can be easily verified. Further, observe that the linear dispersive terms are approximated by using an implicit strategy, in contrast to the nonlinear terms and terms where the variable coefficient M is present which are treated in explicit form. Implicit-explicit schemes (IMEX) were already applied in [25] for scalar dispersive evolution equations. The main advantage of the numerical scheme described here is that at each time step we can solve explicitly the approximations v
^
j
n
+
1, u
^
j
n
+
1, j
=
-
N
/
2
+
1
,
…
,
N
/
2, of the Fourier coefficients of the unknowns v
(
ξ
,
t
) and u
(
ξ
,
t
) from (83), without using implicit Newton-type iterations. Thus, the scheme results in being cheap and its computer implementation is easier.In the scheme, the spatial derivativesu
ξ and v
ξ are computed by exact differentiation of the truncated Fourier series. For instance,(84)
v
ξ
ξ
,
t
=
∑
j
i
w
j
P
j
v
·
,
t
e
i
w
j
ξ
.The numerical calculations presented in this paper were carried out in double precision by using MATLAB R2012b on a Mac platform. The Fourier-type integral appearing in the operator P
j
[
·
] (see (82)) is approximated through the well-known Fast Fourier Transform (FFT) routine.


### 5.2. Approximating Travelling Wave Solutions

In this section we are interested in computing solutions of system (75) over a channel with flat bottom (i.e., when the metric coefficient M is equivalent to 1) of the form (85)
v
ξ
,
t
=
v
~
ξ
-
c
t
,
u
ξ
,
t
=
u
~
ξ
-
c
t
,which are named as travelling wave solutions. The parameter c is called the wave velocity.After dropping the tildes, travelling wave solution(
v
,
u
) with speed c of system (75) must satisfy the following equations:(86)
v
′
′
u
′
′
=
1
β
/
4
c
2
-
2
β
σ
/
3
c
2
2
3
σ
c
2
c
v
-
1
+
α
v
u
c
u
-
v
-
α
2
u
2
.In first place, we are interested in finding approximations to even periodic travelling wave solutions (
v
,
u
) with period 2
l, l
>
0 of system (75). Thus let us introduce truncated cosine expansions for v and u, (87)
v
x
≈
v
0
+
∑
n
=
1
N
/
2
v
n
cos
⁡
n
π
l
x
,
u
x
≈
u
0
+
∑
n
=
1
N
/
2
u
n
cos
⁡
n
π
l
x
,where(88)
v
0
=
1
l
∫
0
l
v
x
d
x
=
1
2
l
∫
0
2
l
v
x
d
x
,
v
n
=
2
l
∫
0
l
v
x
cos
⁡
n
π
x
l
d
x
=
2
2
l
∫
0
2
l
v
x
c
o
s
n
π
x
l
d
x
,and analogous expressions for u. By substituting expressions (87) into (86), evaluating them at the N
/
2
+
1 collocation points(89)
x
j
=
2
l
j
-
1
N
,
j
=
1
,
…
,
N
2
+
1
,we obtain a system of N
+
2 nonlinear equations in the form(90)
F
v
0
,
v
1
,
…
,
v
N
/
2
,
u
0
,
u
1
,
…
,
u
N
/
2
=
0
,where the N
+
2 coefficients v
n
,
u
n are the unknowns. Nonlinear system (90) can be solved by Newton’s iteration. Computation of the cosine series in (87) and the integrals in (88) is performed using the FFT (Fast Fourier Transform) algorithm. The Jacobian of the vector field F
:
R
N
+
2
→
R
N
+
2 is approximated by the second-order accurate formula(91)
J
i
,
j
F
x
≈
F
i
x
+
h
e
j
-
F
i
x
-
h
e
j
2
h
,
j
=
1
,
…
,
N
+
2
,where e
j
=
(
0
,
…
,
1
,
…
,
0
) and h
=
0.01. We stop Newton’s iteration when the relative error between two successive approximations and the value of the vector field F are smaller than 1
0
-
12.The starting point for Newton’s procedure in the periodic frame is taken as(92)
u
0
x
=
c
o
s
π
x
l
,
v
0
x
=
c
o
s
π
x
l
.In second place, we are also interested in approximating solitary wave solutions of system (75), that is, when the functions v
,
u and their derivatives decay to zero at ±
∞. In this case, we can also approximate v
,
u by a truncated Fourier series as in (87) using a large enough length L. To derive an appropriate initial point for Newton’s iteration, we compute an approximate solitary wave solution of system (75) with M
≡
1 and α
,
β
,
σ being small. To accomplish this, we will use original system (38)-(39) in the variables (
v
,
Φ
). Let us observe that from (39)(93)
v
=
-
I
-
β
2
∂
ξ
2
Φ
t
+
β
σ
v
ξ
ξ
-
α
2
Φ
ξ
2
=
-
Φ
t
+
O
α
,
β
,
σ
.Substituting this into (38) and neglecting second-order terms in α
,
β
,
σ, we derive the Benney-Luke type equation for the potential Φ:(94)
-
Φ
t
t
+
Φ
ξ
ξ
+
β
1
-
σ
Φ
ξ
ξ
t
t
-
2
β
3
Φ
ξ
ξ
ξ
ξ
-
2
α
Φ
ξ
Φ
ξ
t
-
α
Φ
t
Φ
ξ
ξ
=
0
.We look for a solitary wave solution of (94) in the form(95)
Φ
ξ
,
t
=
Φ
~
ξ
-
c
t
.Therefore, abandoning the tildes, the function Φ must satisfy the following equation:(96)
1
-
c
2
Φ
′
′
+
β
1
-
σ
c
2
-
2
β
3
Φ
′
′
′
′
+
3
α
c
2
Φ
′
2
′
=
0
.Integrating the equation above and using the fact that Φ and their derivatives decay to zero at ±
∞, we arrive at(97)
1
-
c
2
Φ
′
+
β
1
-
σ
c
2
-
2
3
Φ
′
′
′
+
3
α
c
2
Φ
′
2
=
0
.Multiplying the previous equation by 2
Φ
′
′ and integrating again, we obtain that(98)
1
-
c
2
Φ
′
2
+
β
1
-
σ
c
2
-
2
3
Φ
′
′
2
+
α
c
Φ
′
3
=
0
.Furthermore, letting u
=
Φ
′, we find from (93)(99)
v
=
c
u
-
β
c
2
u
′
′
+
β
σ
c
u
′
′
-
α
2
u
2
.Assume the solution form of Φ
′ to be(100)
u
=
Φ
′
x
=
A
s
e
c
h
2
B
x
,where A
,
B are constants. Then, by replacing this into (98), we get that the constants A
,
B are given by(101)
A
=
c
2
-
1
α
c
,
B
=
c
2
-
1
4
β
1
-
σ
c
2
-
2
/
3
.Now substituting into (99), we arrive at the following approximation of the wave elevation:(102)
v
=
c
A
sech
2
⁡
B
x
-
α
2
A
2
sech
4
⁡
B
x
+
β
c
σ
-
1
2
4
A
B
2
sech
2
⁡
B
x
tanh
2
⁡
B
x
-
2
A
B
2
sech
4
⁡
B
x
.


## 5.1. Temporal Discretization

Note that (80) can be seen as a system of ordinary differential equations where the unknowns are the time-dependent Fourier coefficients of the solutions. To solve it, we use an implicit-explicit scheme (IMEX) in the form(83)
v
^
j
n
+
1
-
v
^
j
n
Δ
t
=
γ
1
u
^
j
n
+
1
+
u
^
j
n
2
+
3
2
N
1
n
-
1
2
N
1
n
-
1
,
u
^
j
n
+
1
-
u
^
j
n
Δ
t
=
γ
2
v
^
j
n
+
1
+
v
^
j
n
2
+
3
2
N
2
n
-
1
2
N
2
n
-
1
.Here v
n, u
n denote the approximations of the unknowns v
(
ξ
,
t
), u
(
ξ
,
t
), respectively, at time t
=
n
Δ
t, where Δ
t is the time step of the method, and N
1
n
,
N
2
n denote the approximations of the functions N
1
,
N
2 evaluated at time n
Δ
t. Similarly, v
^
j
n
,
u
^
j
n denote the approximations to the Fourier transforms of the functions v and u, respectively, with respect to the variable ξ, evaluated at time n
Δ
t. The numerical approach adopted for solving system (75) ensures the scheme results to be linearly unconditionally stable which can be easily verified. Further, observe that the linear dispersive terms are approximated by using an implicit strategy, in contrast to the nonlinear terms and terms where the variable coefficient M is present which are treated in explicit form. Implicit-explicit schemes (IMEX) were already applied in [25] for scalar dispersive evolution equations. The main advantage of the numerical scheme described here is that at each time step we can solve explicitly the approximations v
^
j
n
+
1, u
^
j
n
+
1, j
=
-
N
/
2
+
1
,
…
,
N
/
2, of the Fourier coefficients of the unknowns v
(
ξ
,
t
) and u
(
ξ
,
t
) from (83), without using implicit Newton-type iterations. Thus, the scheme results in being cheap and its computer implementation is easier.In the scheme, the spatial derivativesu
ξ and v
ξ are computed by exact differentiation of the truncated Fourier series. For instance,(84)
v
ξ
ξ
,
t
=
∑
j
i
w
j
P
j
v
·
,
t
e
i
w
j
ξ
.The numerical calculations presented in this paper were carried out in double precision by using MATLAB R2012b on a Mac platform. The Fourier-type integral appearing in the operator P
j
[
·
] (see (82)) is approximated through the well-known Fast Fourier Transform (FFT) routine.


## 5.2. Approximating Travelling Wave Solutions

In this section we are interested in computing solutions of system (75) over a channel with flat bottom (i.e., when the metric coefficient M is equivalent to 1) of the form (85)
v
ξ
,
t
=
v
~
ξ
-
c
t
,
u
ξ
,
t
=
u
~
ξ
-
c
t
,which are named as travelling wave solutions. The parameter c is called the wave velocity.After dropping the tildes, travelling wave solution(
v
,
u
) with speed c of system (75) must satisfy the following equations:(86)
v
′
′
u
′
′
=
1
β
/
4
c
2
-
2
β
σ
/
3
c
2
2
3
σ
c
2
c
v
-
1
+
α
v
u
c
u
-
v
-
α
2
u
2
.In first place, we are interested in finding approximations to even periodic travelling wave solutions (
v
,
u
) with period 2
l, l
>
0 of system (75). Thus let us introduce truncated cosine expansions for v and u, (87)
v
x
≈
v
0
+
∑
n
=
1
N
/
2
v
n
cos
⁡
n
π
l
x
,
u
x
≈
u
0
+
∑
n
=
1
N
/
2
u
n
cos
⁡
n
π
l
x
,where(88)
v
0
=
1
l
∫
0
l
v
x
d
x
=
1
2
l
∫
0
2
l
v
x
d
x
,
v
n
=
2
l
∫
0
l
v
x
cos
⁡
n
π
x
l
d
x
=
2
2
l
∫
0
2
l
v
x
c
o
s
n
π
x
l
d
x
,and analogous expressions for u. By substituting expressions (87) into (86), evaluating them at the N
/
2
+
1 collocation points(89)
x
j
=
2
l
j
-
1
N
,
j
=
1
,
…
,
N
2
+
1
,we obtain a system of N
+
2 nonlinear equations in the form(90)
F
v
0
,
v
1
,
…
,
v
N
/
2
,
u
0
,
u
1
,
…
,
u
N
/
2
=
0
,where the N
+
2 coefficients v
n
,
u
n are the unknowns. Nonlinear system (90) can be solved by Newton’s iteration. Computation of the cosine series in (87) and the integrals in (88) is performed using the FFT (Fast Fourier Transform) algorithm. The Jacobian of the vector field F
:
R
N
+
2
→
R
N
+
2 is approximated by the second-order accurate formula(91)
J
i
,
j
F
x
≈
F
i
x
+
h
e
j
-
F
i
x
-
h
e
j
2
h
,
j
=
1
,
…
,
N
+
2
,where e
j
=
(
0
,
…
,
1
,
…
,
0
) and h
=
0.01. We stop Newton’s iteration when the relative error between two successive approximations and the value of the vector field F are smaller than 1
0
-
12.The starting point for Newton’s procedure in the periodic frame is taken as(92)
u
0
x
=
c
o
s
π
x
l
,
v
0
x
=
c
o
s
π
x
l
.In second place, we are also interested in approximating solitary wave solutions of system (75), that is, when the functions v
,
u and their derivatives decay to zero at ±
∞. In this case, we can also approximate v
,
u by a truncated Fourier series as in (87) using a large enough length L. To derive an appropriate initial point for Newton’s iteration, we compute an approximate solitary wave solution of system (75) with M
≡
1 and α
,
β
,
σ being small. To accomplish this, we will use original system (38)-(39) in the variables (
v
,
Φ
). Let us observe that from (39)(93)
v
=
-
I
-
β
2
∂
ξ
2
Φ
t
+
β
σ
v
ξ
ξ
-
α
2
Φ
ξ
2
=
-
Φ
t
+
O
α
,
β
,
σ
.Substituting this into (38) and neglecting second-order terms in α
,
β
,
σ, we derive the Benney-Luke type equation for the potential Φ:(94)
-
Φ
t
t
+
Φ
ξ
ξ
+
β
1
-
σ
Φ
ξ
ξ
t
t
-
2
β
3
Φ
ξ
ξ
ξ
ξ
-
2
α
Φ
ξ
Φ
ξ
t
-
α
Φ
t
Φ
ξ
ξ
=
0
.We look for a solitary wave solution of (94) in the form(95)
Φ
ξ
,
t
=
Φ
~
ξ
-
c
t
.Therefore, abandoning the tildes, the function Φ must satisfy the following equation:(96)
1
-
c
2
Φ
′
′
+
β
1
-
σ
c
2
-
2
β
3
Φ
′
′
′
′
+
3
α
c
2
Φ
′
2
′
=
0
.Integrating the equation above and using the fact that Φ and their derivatives decay to zero at ±
∞, we arrive at(97)
1
-
c
2
Φ
′
+
β
1
-
σ
c
2
-
2
3
Φ
′
′
′
+
3
α
c
2
Φ
′
2
=
0
.Multiplying the previous equation by 2
Φ
′
′ and integrating again, we obtain that(98)
1
-
c
2
Φ
′
2
+
β
1
-
σ
c
2
-
2
3
Φ
′
′
2
+
α
c
Φ
′
3
=
0
.Furthermore, letting u
=
Φ
′, we find from (93)(99)
v
=
c
u
-
β
c
2
u
′
′
+
β
σ
c
u
′
′
-
α
2
u
2
.Assume the solution form of Φ
′ to be(100)
u
=
Φ
′
x
=
A
s
e
c
h
2
B
x
,where A
,
B are constants. Then, by replacing this into (98), we get that the constants A
,
B are given by(101)
A
=
c
2
-
1
α
c
,
B
=
c
2
-
1
4
β
1
-
σ
c
2
-
2
/
3
.Now substituting into (99), we arrive at the following approximation of the wave elevation:(102)
v
=
c
A
sech
2
⁡
B
x
-
α
2
A
2
sech
4
⁡
B
x
+
β
c
σ
-
1
2
4
A
B
2
sech
2
⁡
B
x
tanh
2
⁡
B
x
-
2
A
B
2
sech
4
⁡
B
x
.


## 6. Description of the Numerical Experiments

In this section we compute some approximations to solutions of system (75) using the numerical schemes described in the previous section. In first place, some travelling wave solutions are approximated by using the numerical scheme discussed in Section 5.2 and then these solutions are checked through computer simulations conducted with numerical solver (83). Finally, we apply it to compute the evolution of a gaussian pulse over a channel with variable depth taking into account the surface tension effect.


### 6.1. Periodic Travelling Wave Solutions

Using Newton’s iteration as explained in Section5.2, we develop some numerical simulations to compute periodic travelling wave solutions of system (75). The results are shown in Figures 2, 3, 4, and 5 for different values of the modelling parameters α
,
β
,
σ and wave speed c. The initial step for Newton’s scheme is taken as in (92). The number of FFT point is 2
7 and we recall that here the channel has constant depth; that is, M
≡
1. In Figure 6 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of numerical scheme (83) at time t
=
10 (using N
=
2
10 and Δ
t
=
0.01), with the profile in Figure 2 shifted at a distance of c
t
=
10
×
1.2
=
12 to the right. Observe that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order 1
0
-
3. Similar results were obtained for the other periodic travelling wave solutions shown in Figures 3, 4, and 5.Figure 2
Periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0, wave speed c
=
1.2, and period L
=
10, obtained after 7 Newton’s iterations.Figure 3
Periodic travelling wave solution of system (75) for α
=
β
=
1, σ
=
0.01, wave speed c
=
0.5, and period L
=
10, obtained after 21 Newton’s iterations.Figure 4
Periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, wave speed c
=
1.2, and period L
=
10, obtained after 7 Newton’s iterations.Figure 5
Periodic travelling wave solution of system (75) for α
=
β
=
0.5, σ
=
0.5, wave speed c
=
0.5, and period L
=
10, obtained after 12 Newton’s iterations.Figure 6
Evolution of a periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0, wave speed c
=
1.2, and period L
=
10, using the numerical scheme (83).


### 6.2. Solitary Wave Solutions

Numerical experiments in the case of solitary wave solutions are shown in Figures7, 8, and 9 for different values of the modelling parameters α
,
β
,
σ and wave speed c. The initial step for Newton’s scheme is taken as in (100)-(102) centered at the position ξ
=
50. The number of FFT point is 2
10 and we recall that here the channel has constant depth; that is, M
≡
1. In this nonperiodic scenario, the spatial computational domain is the interval [
0,100
] which is large enough so that the pulses do not reach the computational boundaries within the time interval modeled. In Figure 10 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of the numerical scheme (83) at time t
=
10 (using N
=
2
10 and Δ
t
=
0.01), with the profile in Figure 2 shifted at a distance of c
t
=
10
×
1.1
=
11 to the right. As in the experiments in the previous section, we see that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order 1
0
-
3. Similar results were obtained for the other periodic travelling wave solutions shown in Figures 8 and 9.Figure 7
Solitary wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, and wave speed c
=
1.1, obtained after 5 Newton’s iterations.Figure 8
Solitary wave solution of system (75) for α
=
β
=
0.5, σ
=
0, and wave speed c
=
1.1, obtained after 5 Newton’s iterations.Figure 9
Solitary wave solution of system (75) for α
=
β
=
0.5, σ
=
0.5, and wave speed c
=
1.5, obtained after 5 Newton’s iterations.Figure 10
Evolution of a solitary wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, and wave speed c
=
1.1, using numerical scheme (83).


### 6.3. Variable Depth

Finally in this section we include some numerical simulations to illustrate the dynamics of incoming Gaussian pulses in the form(103)
u
ξ
,
0
=
v
ξ
,
0
=
e
-
20
x
-
28
2
,governed by (75), under the simultaneous effects of surface tension (σ), dispersion (β), and channel’s topography (coefficient M
(
ξ
)). We further remark that the incoming pulses are located at two units to the left of the irregular part of the channel which covers the interval [
30,40
]. In order to focus on these three phenomena, we only consider the linear case; that is, α
=
0.In Figure11 is displayed the output of numerical scheme (83) for the value of the dispersion parameter β
=
0.01 and without surface tension effect, that is, σ
=
0, for a constant depth channel. The numerical parameters are N
=
2
11 and Δ
t
=
0.01 and the computational domain is the interval [
0,60
] which again is large enough so that the pulses do not reach the computational boundaries. We point out that, as a consequence of dispersion (which produces a difference between the phase velocity of the Fourier modes contained in the pulse), the incident pulses progressively develop an oscillatory tail that propagates to the left, while the leading wave propagates to the right. We contrast this experiment to the dynamics of a pulse under the influence of surface tension which is presented in Figure 12. Note that, in this case, both the oscillatory tail and the wave front propagate to the right side of the computational domain. We conclude that the surface tension effect alters the dispersive characteristics of model (75).Figure 11
Evolution of a Gaussian pulse over a flat channel without surface tension effect. Parameters:α
=
0, β
=
0.01, and σ
=
0.Figure 12
Evolution of a Gaussian pulse over a flat channel with surface tension effect. Parameters:α
=
0, β
=
0.01, and σ
=
0.6.We can anticipate the influence of the additional third-order termβ
σ
∂
ξ
2
v (due to the surface tension effects) on the solutions of system (75) (with α
=
0 and flat channel M
≡
1) by analyzing its linear dispersion relation given by(104)
C
σ
2
=
w
2
k
2
=
1
+
2
β
/
3
k
2
1
+
β
/
2
k
2
2
+
β
σ
k
2
1
+
2
β
/
3
k
2
1
+
β
/
2
k
2
2
>
C
0
2
.We recall that this expression arises by studying the propagation of solutions of system (75) in the form u
=
u
0
e
i
(
k
ξ
-
w
t
), v
=
v
0
e
i
(
k
ξ
-
w
t
), where u
0
,
v
0 are constants.In Figure13, the phase velocities C
σ (with σ
=
0.6) and C
0 (no surface tension) are compared as functions of β
k (k is called the wavenumber). Observe that, in the presence of surface tension, the phase speed C
σ is larger than 1, and it approximates to the value 8
σ
/
3 as β
k tends to infinity. In contrast, without surface tension effects, the phase speed C
0 is lower than 1 and Fourier components with high wavenumber travel more slowly than the Fourier modes with small wavenumber. These observations are in accordance with the results presented in Figures 11 and 12. We further point out that the fact that the phase velocity in system (75) is bounded by 8
σ
/
3 for all β
k has an advantage in improving the stability of numerical solvers formulated to approximate solutions of this system. In particular, the properties of numerical solver (83) employed in the present paper, such as the convergence, stability, and error estimates, will be studied in a future research.Figure 13
The phase velocity for system (75).After analyzing the flat bottom case, we consider a numerical simulation in Figure15 where the channel’s topography is variable. The topography considered is shown in Figure 14 together with the metric coefficient M
(
ξ
,
β
) computed using the Schwarz-Christoffel Toolbox designed for MATLAB by Driscoll [30]. Note that the irregular part covers the interval [
30,40
] (10 units long). A nonnegligible value of surface tension σ
=
0.6 is incorporated into model equations (75). Again the level of dispersion is taken as β
=
0.01. The numerical parameters are the same as in the previous computer simulation. In addition to the effect of the surface tension on the direction of propagation of the oscillatory tail, we now observe irregular fluctuations (coda) developing behind the wave front and propagating to the left end of the computational domain. This tail signal is due to the multiple reflections of the incoming pulses produced by every variation of the channel’s depth. These phenomena have been studied in previous works using other water wave models that neglect the surface tension effect [9–11]. In the present paper, our purpose was to initiate the study of these wave phenomena and numerical scheme (83) and system (75) have proved to be valuable tools in developing this research in a future work.Figure 14
The metric coefficientM
(
ξ
,
β
=
0.01
) (in solid line) compared with a (nondifferentiable) piecewise linear topography (in dashed line). Observe that the coefficient M is smooth even though the physical topography is not differentiable. The irregular part of the bottom is located over the interval [
30,40
] and the channel’s bottom is flat outside this interval.Figure 15
Evolution of a Gaussian pulse over a channel with the variable depth shown in Figure14 and taking into account surface tension effect. Parameters: α
=
0, β
=
0.01, and σ
=
0.6.To see in more detail the effect of surface tension on the wave-topography interaction, in Figure16, we show the evolution of the Gaussian pulse (105)
u
ξ
,
0
=
v
ξ
,
0
=
e
-
20
x
-
198
2on the surface of a channel with a large irregular bottom covering the interval [
200,300
] (10 times larger than in the previous numerical simulation) in the case of σ
=
0 (no surface tension) and σ
=
0.6. As in the previous experiments, the physical depth is described by a linear piecewise function and the coefficient M
(
ξ
,
β
) is computed again using Driscoll’s package for MATLAB. The numerical parameters for this computer simulation are N
=
2
14, Δ
t
=
0.0025, α
=
0
,
β
=
0.01, and the computational domain is the interval [
0,400
]. Observe that the pulse without the presence of surface tension appears to be delayed with respect to the pulse propagating under the influence of a nonnegligible value of surface tension σ
=
0.6 and the long channel’s topography introduced. These phenomena and further research using full model equations (35) for the case of a rapidly varying bottom (i.e., γ
≪
1) will be studied in detail in a future work. We remark that numerical scheme (83) could be adapted to solve full system (35).Figure 16
Propagation of a Gaussian pulse over a channel with variable depth. Observe the effect of surface tension on the dynamics of the pulse. Parameters:α
=
0, β
=
0.01. Solid line: σ
=
0.6. Dashed line: σ
=
0.


## 6.1. Periodic Travelling Wave Solutions

Using Newton’s iteration as explained in Section5.2, we develop some numerical simulations to compute periodic travelling wave solutions of system (75). The results are shown in Figures 2, 3, 4, and 5 for different values of the modelling parameters α
,
β
,
σ and wave speed c. The initial step for Newton’s scheme is taken as in (92). The number of FFT point is 2
7 and we recall that here the channel has constant depth; that is, M
≡
1. In Figure 6 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of numerical scheme (83) at time t
=
10 (using N
=
2
10 and Δ
t
=
0.01), with the profile in Figure 2 shifted at a distance of c
t
=
10
×
1.2
=
12 to the right. Observe that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order 1
0
-
3. Similar results were obtained for the other periodic travelling wave solutions shown in Figures 3, 4, and 5.Figure 2
Periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0, wave speed c
=
1.2, and period L
=
10, obtained after 7 Newton’s iterations.Figure 3
Periodic travelling wave solution of system (75) for α
=
β
=
1, σ
=
0.01, wave speed c
=
0.5, and period L
=
10, obtained after 21 Newton’s iterations.Figure 4
Periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, wave speed c
=
1.2, and period L
=
10, obtained after 7 Newton’s iterations.Figure 5
Periodic travelling wave solution of system (75) for α
=
β
=
0.5, σ
=
0.5, wave speed c
=
0.5, and period L
=
10, obtained after 12 Newton’s iterations.Figure 6
Evolution of a periodic travelling wave solution of system (75) for α
=
β
=
0.3, σ
=
0, wave speed c
=
1.2, and period L
=
10, using the numerical scheme (83).


## 6.2. Solitary Wave Solutions

Numerical experiments in the case of solitary wave solutions are shown in Figures7, 8, and 9 for different values of the modelling parameters α
,
β
,
σ and wave speed c. The initial step for Newton’s scheme is taken as in (100)-(102) centered at the position ξ
=
50. The number of FFT point is 2
10 and we recall that here the channel has constant depth; that is, M
≡
1. In this nonperiodic scenario, the spatial computational domain is the interval [
0,100
] which is large enough so that the pulses do not reach the computational boundaries within the time interval modeled. In Figure 10 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of the numerical scheme (83) at time t
=
10 (using N
=
2
10 and Δ
t
=
0.01), with the profile in Figure 2 shifted at a distance of c
t
=
10
×
1.1
=
11 to the right. As in the experiments in the previous section, we see that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order 1
0
-
3. Similar results were obtained for the other periodic travelling wave solutions shown in Figures 8 and 9.Figure 7
Solitary wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, and wave speed c
=
1.1, obtained after 5 Newton’s iterations.Figure 8
Solitary wave solution of system (75) for α
=
β
=
0.5, σ
=
0, and wave speed c
=
1.1, obtained after 5 Newton’s iterations.Figure 9
Solitary wave solution of system (75) for α
=
β
=
0.5, σ
=
0.5, and wave speed c
=
1.5, obtained after 5 Newton’s iterations.Figure 10
Evolution of a solitary wave solution of system (75) for α
=
β
=
0.3, σ
=
0.3, and wave speed c
=
1.1, using numerical scheme (83).


## 6.3. Variable Depth

Finally in this section we include some numerical simulations to illustrate the dynamics of incoming Gaussian pulses in the form(103)
u
ξ
,
0
=
v
ξ
,
0
=
e
-
20
x
-
28
2
,governed by (75), under the simultaneous effects of surface tension (σ), dispersion (β), and channel’s topography (coefficient M
(
ξ
)). We further remark that the incoming pulses are located at two units to the left of the irregular part of the channel which covers the interval [
30,40
]. In order to focus on these three phenomena, we only consider the linear case; that is, α
=
0.In Figure11 is displayed the output of numerical scheme (83) for the value of the dispersion parameter β
=
0.01 and without surface tension effect, that is, σ
=
0, for a constant depth channel. The numerical parameters are N
=
2
11 and Δ
t
=
0.01 and the computational domain is the interval [
0,60
] which again is large enough so that the pulses do not reach the computational boundaries. We point out that, as a consequence of dispersion (which produces a difference between the phase velocity of the Fourier modes contained in the pulse), the incident pulses progressively develop an oscillatory tail that propagates to the left, while the leading wave propagates to the right. We contrast this experiment to the dynamics of a pulse under the influence of surface tension which is presented in Figure 12. Note that, in this case, both the oscillatory tail and the wave front propagate to the right side of the computational domain. We conclude that the surface tension effect alters the dispersive characteristics of model (75).Figure 11
Evolution of a Gaussian pulse over a flat channel without surface tension effect. Parameters:α
=
0, β
=
0.01, and σ
=
0.Figure 12
Evolution of a Gaussian pulse over a flat channel with surface tension effect. Parameters:α
=
0, β
=
0.01, and σ
=
0.6.We can anticipate the influence of the additional third-order termβ
σ
∂
ξ
2
v (due to the surface tension effects) on the solutions of system (75) (with α
=
0 and flat channel M
≡
1) by analyzing its linear dispersion relation given by(104)
C
σ
2
=
w
2
k
2
=
1
+
2
β
/
3
k
2
1
+
β
/
2
k
2
2
+
β
σ
k
2
1
+
2
β
/
3
k
2
1
+
β
/
2
k
2
2
>
C
0
2
.We recall that this expression arises by studying the propagation of solutions of system (75) in the form u
=
u
0
e
i
(
k
ξ
-
w
t
), v
=
v
0
e
i
(
k
ξ
-
w
t
), where u
0
,
v
0 are constants.In Figure13, the phase velocities C
σ (with σ
=
0.6) and C
0 (no surface tension) are compared as functions of β
k (k is called the wavenumber). Observe that, in the presence of surface tension, the phase speed C
σ is larger than 1, and it approximates to the value 8
σ
/
3 as β
k tends to infinity. In contrast, without surface tension effects, the phase speed C
0 is lower than 1 and Fourier components with high wavenumber travel more slowly than the Fourier modes with small wavenumber. These observations are in accordance with the results presented in Figures 11 and 12. We further point out that the fact that the phase velocity in system (75) is bounded by 8
σ
/
3 for all β
k has an advantage in improving the stability of numerical solvers formulated to approximate solutions of this system. In particular, the properties of numerical solver (83) employed in the present paper, such as the convergence, stability, and error estimates, will be studied in a future research.Figure 13
The phase velocity for system (75).After analyzing the flat bottom case, we consider a numerical simulation in Figure15 where the channel’s topography is variable. The topography considered is shown in Figure 14 together with the metric coefficient M
(
ξ
,
β
) computed using the Schwarz-Christoffel Toolbox designed for MATLAB by Driscoll [30]. Note that the irregular part covers the interval [
30,40
] (10 units long). A nonnegligible value of surface tension σ
=
0.6 is incorporated into model equations (75). Again the level of dispersion is taken as β
=
0.01. The numerical parameters are the same as in the previous computer simulation. In addition to the effect of the surface tension on the direction of propagation of the oscillatory tail, we now observe irregular fluctuations (coda) developing behind the wave front and propagating to the left end of the computational domain. This tail signal is due to the multiple reflections of the incoming pulses produced by every variation of the channel’s depth. These phenomena have been studied in previous works using other water wave models that neglect the surface tension effect [9–11]. In the present paper, our purpose was to initiate the study of these wave phenomena and numerical scheme (83) and system (75) have proved to be valuable tools in developing this research in a future work.Figure 14
The metric coefficientM
(
ξ
,
β
=
0.01
) (in solid line) compared with a (nondifferentiable) piecewise linear topography (in dashed line). Observe that the coefficient M is smooth even though the physical topography is not differentiable. The irregular part of the bottom is located over the interval [
30,40
] and the channel’s bottom is flat outside this interval.Figure 15
Evolution of a Gaussian pulse over a channel with the variable depth shown in Figure14 and taking into account surface tension effect. Parameters: α
=
0, β
=
0.01, and σ
=
0.6.To see in more detail the effect of surface tension on the wave-topography interaction, in Figure16, we show the evolution of the Gaussian pulse (105)
u
ξ
,
0
=
v
ξ
,
0
=
e
-
20
x
-
198
2on the surface of a channel with a large irregular bottom covering the interval [
200,300
] (10 times larger than in the previous numerical simulation) in the case of σ
=
0 (no surface tension) and σ
=
0.6. As in the previous experiments, the physical depth is described by a linear piecewise function and the coefficient M
(
ξ
,
β
) is computed again using Driscoll’s package for MATLAB. The numerical parameters for this computer simulation are N
=
2
14, Δ
t
=
0.0025, α
=
0
,
β
=
0.01, and the computational domain is the interval [
0,400
]. Observe that the pulse without the presence of surface tension appears to be delayed with respect to the pulse propagating under the influence of a nonnegligible value of surface tension σ
=
0.6 and the long channel’s topography introduced. These phenomena and further research using full model equations (35) for the case of a rapidly varying bottom (i.e., γ
≪
1) will be studied in detail in a future work. We remark that numerical scheme (83) could be adapted to solve full system (35).Figure 16
Propagation of a Gaussian pulse over a channel with variable depth. Observe the effect of surface tension on the dynamics of the pulse. Parameters:α
=
0, β
=
0.01. Solid line: σ
=
0.6. Dashed line: σ
=
0.


## 7. Conclusions

In this paper we introduced a new model for describing the propagation of a small amplitude water wave over the surface of a shallow channel with a possibly varying depth. It was derived formally as an asymptotical approximation of the full potential theory equations in the weakly-nonlinear, weakly-dispersive regime. Existence and uniqueness of a solution to the initial value problem associated with (5) were also studied using Banach’s fixed point argument and the classical semigroup theory in appropriate Sobolev spaces. The main features of model (5) considered in this work are the following:(i)
It can be applied to study wave propagation over a shallow channel with a discontinuous or nondifferentiable bottom provided that the bottom’s fluctuations satisfy|
n
|
<
1. Observe that in the derivation of system (5) the function n was not assumed to be small nor continuous since the conformal mapping has a smoothing effect on the topography. Note that the effective variable coefficient M
(
ξ
) that appears in final system (5) is defined as a convolution between the function describing the bottom’s fluctuations n and the s
e
c
h
2 function (see (25)). This fact implies that M
(
ξ
) and all variable coefficients in system (5) are infinitely differentiable, even when the function n is discontinuous or nondifferentiable. In Figure 14, the coefficient M
(
ξ
) is displayed for a (nondifferentiable) piecewise linear topography. We see that the coefficient M
(
ξ
) is a regularized version of the physical topography.
(ii)
One-dimensional system (5) models bidirectional waves and it incorporates the simultaneous effects of surface tension and variable depth upon the shape of a water wave that propagates on the surface of an irregular shallow channel.
(iii)
Furthermore, in the derivation of system (5), we do not eliminate the wave elevation η (which prevents us from neglecting additional terms of order O
(
α
2
,
β
2
)), as done, for example, in the model derived by Milewski [23]. This means that we did not reduce (30)-(31) to only one equation for the potential velocity Φ where η does not appear. Therefore, new formulation (5) is expected to be a more accurate approximation of the original potential theory equations.In second place, we introduced a numerical Galerkin-spectral discretization of second order in time and spectral accuracy in space to approximate solutions of system (5) and a Fourier-collocation strategy combined with a Newton-type iterative procedure to compute its travelling wave solutions. In all computer simulations, these solutions were observed to propagate approximately without changing their shape within a large time interval.



---

*Source: 805625-2015-09-30.xml*