https://doi.org/10.1007/s42286-020-00035-9 O R I G I N A L A R T I C L E
Hyperbolicity of the Modulation Equations for the Serre–Green–Naghdi Model
Sergey Tkachenko1· Sergey Gavrilyuk1· Keh-Ming Shyue2
Received: 4 November 2019 / Accepted: 2 April 2020
© Springer Nature Switzerland AG 2020
Abstract
The Serre–Green–Naghdi (SGN) system is one of the most useful dispersive models for the description of long water waves having good mathematical and physical properties.
First, the model is a mathematically justified approximation of the exact water-wave problem. Second, the SGN equations are the Euler–Lagrange equations derived from Hamilton’s principle of surface water waves with an approximate Lagrangian func- tional which appears naturally from the full Lagrangian. Finally, the equations are Galilean invariant which is necessary for physically relevant mathematical models.
In this work, we derive the Whitham modulation equations for the SGN model and prove that these equations are strictly hyperbolic for any wave amplitude. This prop- erty can be interpreted by concluding that the periodic wave solutions of the SGN equations are modulationally stable. Numerical simulations of the full SGN equations are also shown, and these results confirm the theoretical stability results found using the Whitham modulation theory.
Keywords Dispersive shallow water equations· Whitham equations · Hyperbolicity
1 Introduction
One usually uses two methods to obtain the modulation equations for a reversible dispersive system: the Whitham averaged Lagrangian method [44] or averaging of the corresponding conservation laws (cf. [5,29]). Both of them give the same system for
B
Sergey Gavrilyuksergey.gavrilyuk@univ-amu.fr Sergey Tkachenko
sergey.tkachenko@univ-amu.fr Keh-Ming Shyue
shyue@ntu.edu.tw
1 Aix Marseille Univ, CNRS, IUSTI, UMR 7343, Marseille, France
2 Institute of Applied Mathematical Sciences, National Taiwan University, Taipei 106, Taiwan
slowly varying wave train characteristics (wave length and amplitude, for example).
If the corresponding system of modulation equations is hyperbolic (elliptic), one says that the corresponding wave trains are modulationally stable (unstable). The relation between the modulational instability and, for example, the classical spectral instability is not completely understood. An important result for a class of Hamiltonian systems has been obtained in [4]: the hyperbolicity of modulation equations is necessary for the spectral stability of periodic traveling waves.
The modulational instability (wavetrain instability) is often generically called
‘Benjamin–Feir instability’ even if, formally, this last instability concerns only about the surface gravity water waves (Benjamin and Feir [2] and Zakharov [45] in the case of deep-water waves, and Benjamin [3] in the case of finite-depth waves). We refer interested readers to the article [46] where the history of the modulational instabil- ity theory is presented. In the case of small amplitudes, the hyperbolicity (ellipticity) condition can easily be formulated in terms of the non-linear amplitude dependent dispersion relation (see [44], chapter 15). A recent application of such an approach can be found in [36] for a ‘conduit’ equation. However, the study of hyperbolicity of the modulation equations in the case of large-amplitude solutions is a more diffi- cult problem. For integrable systems, the hyperbolicity of modulation equations and existence of the Riemann invariants were established, for example, for Korteweg-de Vries equation (KdV equation) [44], for nonlinear Schrödinger equation (NLS equa- tion) [39], for sine-Gordon equation [15,26], and for Benjamin–Ono equation [7] (see also the book [29] for further references). For non-integrable systems, one can cite [27] where the Whitham equation was studied and modulational instability for short enough waves was shown, or [35] and [28] where the regions of modulational and spectral stability for roll waves to the Saint-Venant equations were determined. In general, a ‘right choice’ of unknowns in which the modulation equations are written is necessary to have explicit (or almost explicit) expressions of the corresponding characteristic values. Such a choice is not at all obvious.
The aim of this work is to study the modulational stability of periodic waves to Serre–Green–Naghdi equations (SGN equations). One-dimensional SGN equations can be written in the Eulerian coordinates in the form [22,23,42,43]:
ht + (hu)x = 0,
(hu)t+ (hu2+ p)x= 0, (he)t+ (hue + pu)x = 0,
(1)
with
p= gh2 2 +1
3h2D2h
Dt2, e =u2 2 +gh
2 +1 6
Dh Dt
2
, D
Dt = ∂
∂t + u ∂
∂x. (2) Here, h is the fluid depth, u is the averaged over the fluid depth velocity, and p is the integrated over the fluid depth pressure. If L0is a characteristic wave length, and H0 is the characteristic water depth; we define the dimensionless small parameter β = H02/L20. The SGN equations are obtained by depth-averaging the Euler system
and keeping in the resulting set of equations only first-order terms inβ without mak- ing any assumptions on the amplitude of the waves. The third equation in (1) (the energy equation) is a consequence of the mass and momentum equations (the first two equations). The fourth conservation law (generalized Bernoulli conservation law) can also be written here (cf. [17,19]).
Mathematical justification of this model and some related systems can be found in [8,30,33,37,41]. A variational formulation of the SGN equations is given in [17,32,38].
The linear stability of solitary waves of small amplitude to the SGN equations was established in [32]. Also, it has been mentioned there that numerically, the solitary waves are stable for any wave amplitude. Recent years have seen increased activity in both the study of qualitative properties of the solutions to the SGN system and in the development of numerical discretization techniques [6,9,14,18,21,31,34].
In [11], the Riemann problem for the SGN equations was examined. Earlier, the Riemann problem was mainly studied for integrable systems in [24] for the KdV equa- tion, and in [10,25] for the NLS equation. Recently, this problem has received much attention for non-integrable systems of equations mainly because of the dispersive shocks commonly present in physics [12]. In [11], the wave number, amplitude, aver- age fluid depth, and average velocity have been chosen as primary variables to study the dispersive shocks of the SGN equations. This choice is quite natural, because, for example, the leading edge of the dispersive shock corresponds to the limit of small wave numbers, while the trailing edge is the limit of small amplitudes. Therefore, the asymptotic study in the limit of small wave numbers or small amplitudes is important to predict the solution behaviour.
To capture better the case of moderate and large amplitude waves, one can try to use other variables. Even if a priori, they may be not necessarily physically tractable, they could be useful to parameterize globally the generic solution to the modulation equations. In particular, it could help to determine the regions of modulational stability and instability for waves of arbitrary amplitude.
The structure of the article is as follows. The system of four modulation equations is derived in Sects.2,3, and4. In Sect.5, the averaged quantities are expressed as func- tions of the roots of the third-order polynomial determining the fluid depth behaviour, and the phase velocity. The non-conservative form of the modulation equations and their hyperbolicity analysis are given in Sects.6,7. Numerical tests showing the wave- train stability for the full SGN equations are presented in Sect.8. Technical details are described in Appendix.
2 Averaged Conservation Laws for the SGN Model
A formal derivation of the modulation equations to the SGN equations (even in a more general formulation which contained, in particular, equations of bubbly fluids) can be found in [16,21]. However, the analysis of the hyperbolicity for such a general formulation was not performed there. Here, we will concentrate on SGN equations and will use the approach based on the averaging of conservation laws.
Suppose that the unknowns h, u (and also p and e which are functions of these variables and their derivatives) depend on the rapid travelling coordinateξ = x − Dt
and slow variables X = μx, T = μt (see AppendixAfor the details). Here, D is the travelling wave velocity. Let us introduce the followingμ–expansion ansatz for h, u,
p, and e:
h(ξ, X, T ) = h0(ξ, X, T ) + μh1(ξ, X, T ) + O(μ2), u(ξ, X, T ) = u0(ξ, X, T ) + μu1(ξ, X, T ) + O(μ2), p(ξ, X, T ) = p0(ξ, X, T ) + μp1(ξ, X, T ) + O(μ2), e(ξ, X, T ) = e0(ξ, X, T ) + μe1(ξ, X, T ) + O(μ2).
Here, all the terms are supposed to be L-periodic with respect toξ, where L is also a slowly varying function of X and T . Taking into account the following transformations of the partial derivatives with respect to time and space:
∂
∂t = −D ∂
∂ξ + μ ∂
∂T, ∂
∂x = ∂
∂ξ + μ ∂
∂ X; (3)
in zero-order approximation with respect toμ, we obtain:
− Dh0ξ+ (h0u0)ξ = 0,
− D(h0u0)ξ+ (h0u20+ p0)ξ = 0,
− D(h0e0)ξ+ (h0u0e0+ p0u0)ξ = 0.
(4)
After averaging the first order equations over the period L, one gets the following system:
(h0)T + (h0u0)X = 0, (h0u0)T + (h0u20+ p0)X = 0, (h0e0)T +
h0u0e0+ p0u0
X = 0.
(5)
Notice that here we used the fact that the averaging procedure and the derivation with respect to slow variables commute (cf. [5,29,44]). System (5) will be written in closed form. For this, explicit periodic solutions to (4) will be presented.
3 Stationary Periodic Solutions
We will now show that the equations of zero -order approximation (4) admit periodic solutions. We rewrite (4) in the form:
− Dh+ (hu)= 0,
− D(hu)+ (hu2+ p)= 0,
− D(he)+ (hue + pu)= 0.
Here and further, the zero index is omitted and primes stand for ∂
∂ξ. The first equation reads:
−Dh+ (hu)= 0.
The integration gives:
h(u − D) = m = const. (6)
When we write here and further, m = const, for example, we mean that m does not depend on rapid variableξ: it is a function of just two variables T and X. The second equation can be integrated as:
p = i −m2
h , i = const. (7)
The second-derivative term in (2) can be transformed. Keeping only the zero powers ofμ, we rewrite the second derivative of h as:
D2h
Dt2 = (u − D)
(u − D)h
= (u − D)h h
(u − D)h h h
= m h
mh h
=m2 h
h h
.
Thus, we have:
D2h Dt2 = m2
h
h h
. Hence, the pressure expression in (2) reads:
p= gh2 2 +1
3m2h
h h
.
Replacing the pressure expression into (7), one obtains:
m2 h + gh2
2 +1 3m2h
h h
= i.
Multiplying both side of the equation by h/m2h2, we have:
h h3 + gh
2m2 +1 3
h h
h h
= i h m2h2, which can be rewritten in the form:
1 6
h h
2 + h
h3+ gh
2m2 = i h m2h2. Integrating the equation once leads to:
1 6
h h
2
− 1
2h2 + gh
2m2 = − i
m2h + , = const.
Thus, the equation for h is given as:
h2
= 3 − 6i
m2h+ 6h2− 3g
m2h3. (8)
Denote the polynomial in the right-hand side as F3(h):
F3(h) = 3 − 6i
m2h+ 6h2− 3g m2h3, or, equivalently:
F3(h) = 3g m2
m2 g −2i
g h+2m2 g h2− h3
= 3g
m2(h − h0)(h − h1)(h2− h), where h0 < h1 < h2are the roots of F3(h). Using Vieta’s formulas one can write F3(h) in the following way:
F3(h) = 3
I3(I3− I2h+ I1h2− h3), where
I1= h0+ h1+ h2, I2= h0h1+ h1h2+ h0h2, I3= h0h1h2.
(9)
Identifying the coefficients of F3(h), one obtains:
I1= 2m2
g , I2=2i
g, I3= m2
g . (10)
We are searching for the periodic solutions of (8), oscillating between two real positive roots h1and h2. Since they are real, the third root h0is real too. Moreover, the last formula in (10) implies that h0is necessarily positive.
The periodic solution that oscillates between h1and h2is given by the formula:
h= h1+ (h2− h1)cn2(α ξ; k), α2=3 4
(h2− h0)
h0h1h2 , k2= h2− h1
h2− h0. (11) Here, the Jacobi elliptic function cn(u; k) is defined as:
cn(u; k) = cos(ϕ(u, k)), whereϕ(u, k) is obtained implicitly from the relation:
ϕ
0
dθ
1− k2sin2(θ) = u.
The wavelength L can explicitly be given as:
L =
ξ2
ξ1
dξ = 2
h2 h1
√dh
F3(h) = 2
h0h1h2
3
h2 h1
√dh P3(h), where the interval[ξ1, ξ2] has the length L and
P3(h) = (h − h0)(h − h1)(h2− h). (12) The wavelength is thus completely defined by the roots h0, h1, and h2. The averaging of any arbitrary function of f(h) reads:
f(h) = 1 L
ξ2
ξ1
f(h)dξ = 2 L
h2
h1
f(h)dh
√F3(h) =
h2
h1
f(h)dh
√P3(h) h2
h1
√dh
P3(h). (13)
4 Averaged Equations
Consider the first-order part of (5) (‘zero’ index is omitted):
(h)T + (hu)X = 0, (hu)T + (hu2+ p)X = 0, (he)T +
hue+ pu
X = 0.
(14)
In the following, we will express all averaged quantities in (14) in terms of four unknowns: h0, h1, h2, and D.
The flux in the first equation of (14) is:
hu= h(u − D + D) = m + Dh = m + Dh = hU, U = hu
h . (15) We introduced here the depth averaged velocity U . In terms of this velocity the mass equation can be rewritten in standard form:
(h)T + (hU)X = 0. (16)
Since
hu2= h(u − D + D)2
= h(u − D)2+ 2h(u − D)D + D2h
=
h2(u − D)2 h
+ 2h(u − D)D + D2h
= m2h−1+ 2Dm + D2h,
and
p= i − m2h−1, the flux in the second equation of (14) is:
hu2+ p = i + 2Dm + D2h = hU2+ i − m2 h .
The last two terms represent a combination of the pressure (defined up to multiplicative constant which is the fluid density) first integrated over the water depth and then averaged over the wave period, and the corresponding quadratic velocity correlation.
Together, the two terms form an ‘effective pressure’. One can prove that such an effective pressure is always positive:
i−m2 h = g
2
⎛
⎝I2− 2I3
h2 h1
√dh P3(h)
h2
h1
√hdh P3(h)
⎞
⎠ > 0.
In the third equation, we need to calculate he, hue, and pu. Let us remark that for the travelling wave solutions, one has:
Dh Dt
2
=
(u − D)h2
= m2h2 h2. Also, it follows from (8) that
h2 h
= −6i
m2 + 3h−1+ 6h − 3g m2h2,
and
h h
2
= 6 + 3h−2− 6i
m2h−1− 3g m2h.
The averaged energy is:
e= u2 2 +gh
2 +1 6
Dh Dt
2
= u2 2 +gh
2 +1 6
(u − D)h2
= 1 2
(u − D)2+ 2D(u − D) + D2 +g
2h+1 6
h2(u − D)2h2 h2
= 1 2
h2(u − D)2 h2
+ D
h(u − D) h
+1
2D2+g 2h+1
6m2
h2 h2
= 1
2m2h−2+ Dmh−1+1 2D2+g
2h+m2 6
h h
2
.
Thus, the averaged energy reads:
e= 1
2D2+ m2 + m2h−2+ (Dm − i)h−1.
It can be also expressed as a function of U (instead of D) and h0, h1, h2. The volume average energy is:
he=hu2 2 + gh2
2 +h 6
Dh Dt
2
=1 2
2Dm+ m2h−1+ D2h
+gh2 2 +h
6
(u − D)h2
=1 2
2Dm+ m2h−1+ D2h
+g 2h2+1
6
h2(u − D)2h2 h
=1 2
2Dm+ m2h−1+ D2h
+g
2h2+m2 6
h2 h
.
Thus:
he= Dm − i + m2h−1+
1
2D2+ m2
h,
and
hue= h(u − D + D)e = me + Dhe = me + Dhe
= 3
2D2m− i D + m3 + m3h−2+ (2Dm2− mi)h−1+
1
2D3+ m2D
h.
Also, one has:
pu= iu − m2uh−1
= ih
h(u − D + D) − m2h
h2(u − D + D)
= i
hh(u − D) + i D −m2
h2h(u − D) − m2Dh−1
= i D − m3h−2− (Dm2− mi)h−1.
Now we are able to write the modulation equations because all quantities we need are given explicitly:
hu= hD + m, hu2+ p = hD2+ 2m D + i,
he=1
2h D2+ m D − i + m2h + m2h−1, hue+ pu = 1
2h D3+ m2hD + 3
2m D2+ m3 + m2h−1D.
(17)
The last step would be to replace the integration constants m, i , and by their expres- sions in terms of invariants Ii, i= 1, 2, 3, using (10):
m2= gI3(m = sgn(m) g I3), i =1
2g I2,
= 1 2
I1
I3.
(18)
The negative (positive) sign of m corresponds to the right (left) facing periodic waves.
We introduce the synthesis of both notations in order to obtain the simplest form of the equations. Basically, we will describe everything in terms of m, i , and I1. Using (18), one can eliminate the dependence on in (17):
hu= hD + m, hu2+ p = hD2+ 2m D + i,
he= 1
2h D2+ m D − i +1
2g I1h+ m2h−1, hue+ pu = 1
2h D3+1
2g I1h D+3
2m D2+1
2g I1m+ m2h−1D.
(19)
Complemented by equation (1/L)t+ (D/L)x = 0 (conservation of waves, see AppendixA), equations (5) will finally be written as:
LT − L DX+ DLX = 0, hT +
m+ hD
X = 0,
m+ hD
T +
h D2+1
2g I2+ 2m D
X
= 0,
1
2h D2+1
2g I1h−1
2g I2+ gI3h−1+ m D
T
+
1
2h D3+1
2g I1h D+ gI3h−1D+3
2m D2+1 2mg I1
X
= 0.
(20)
The depth averaged velocity U = uh/h is linked to the phase speed D by (15):
D= U −m h.
If one uses U instead of D, Eq. (20) complemented by Eq. (30) for the wavelength (see AppendixA) will be written as:
(1/L)T +
U/L − m/(Lh)
X = 0, hT +
hU
X = 0,
hU
T +
hU2+ P
X = 0, P = 1
2g I2− m2/h,
1
2hU2+ hE
T
+
1
2hU3+ U
P + hE
− m3h−1/h + m3/h2
X
= 0, hE = m2h−1−1
2m2/h + 1 2gh
I1− I2/h .
(21)
Let us remark that (21) are Galilean invariant, i.e., the equations will not change after the change of variables:
X → X + V T , U → U + V , V = const.
The other variables do not change. The Galilean invariance implies, in particular, that the corresponding characteristic polynomial is also Galilean invariant. This fact will be used later. It appears that the quasilinear form of (20) in variables D, h0, h1, and h2
is simpler than that in variables U , h0, h1, and h2. We now rewrite (20) in quasilinear form.
5 Expressions for the Main Averaged Variables
The expressions of h, h−1, and L in terms of h0, h1, and h2 are (for proof, see AppendixB):
h= h0+ (h2− h0)E(k)
K(k), h−1= (n, k)
h2K(k), L = 4
h0h1h2
3
K(k)
√h2− h0
.
Then, one can write the following differentials : dh= 0dh0+ 1dh1+ 2dh2, dh−1= 0dh0+ 1dh1+ 2dh2, d L= 0dh0+ 1dh1+ 2dh2, d I1= dh0+ dh1+ dh2,
d I2= (h1+ h2)dh0+ (h0+ h2)dh1+ (h0+ h1)dh2, d I3= h1h2dh0+ h0h2dh1+ h0h1dh2,
dm= m 2
dh0
h0 +dh1
h1 +dh2
h2
.
Herei, i, and i (i = 0, 1, 2) read (for proof, see AppendixC):
0= 1
2− h2− h0
2(h1− h0) E2(k) K2(k), 1= h2− h0
2(h2− h1)−h2− h0
h2− h1
E(k)
K(k)+ (h2− h0)2 2(h2− h1)(h1− h0)
E2(k) K2(k), 2= − h1− h0
2(h2− h1)+h2− h0
h2− h1
E(k)
K(k)− h2− h0
2(h2− h1) E2(k) K2(k),
0= 1
2h0(h1− h0) E(k) K(k)− 1
2h0h2
(n, k)
K(k) − 1
2h2(h1− h0)
(n, k)E(k) K2(k) ,
1= 1
2h1(h2− h1)− h2− h0
2h1(h2− h1)(h1− h0) E(k)
K(k)− 1
2h1(h2− h1) (n, k)
K(k)
+ h2− h0
2h2(h2− h1)(h1− h0)
(n, k)E(k) K2(k) ,
2= − 1
2h2(h2− h1)+ 1 2h2(h2− h1)
E(k)
K(k)+ h1
2h22(h2− h1) (n, k)
K(k)
− 1
2h2(h2− h1)
(n, k)E(k) K2(k) , 0= 2
√3
√
h0h1h2
(h1− h0)√ h2− h0
E(k) + h1h2
√h2− h0
√h0h1h2
K(k)
,
1= 2
√3
−
√h2− h0
√h0h1h2
(h2− h1)(h1− h0)E(k) + h0h22 (h2− h1)√
h2− h0
√h0h1h2
K(k)
,
2= 2
√3
√
h0h1h2
(h2− h1)√ h2− h0
E(k) − h0h21 (h2− h1)√
h2− h0
√h0h1h2
K(k)
.
The formulas fork, kand k, k = 0, 1, 2 were verified by hand calculations and with Wolfram Mathematica. One must pay attention to the fact that the complete elliptic integrals which we use depend on elliptic modulus k, while Wolfram Mathematica uses
the definition from Abramovitz and Stegun [1] where the complete elliptic integrals depend on parameter m= k2(do not confound the notations m with m coming from the mass conservation equation).
6 Nonconservative Modulation Equations
The modulation equations (20) can be written in the following developed form:
0h0T+ 1h1T + 2h2T − L DX + D 0h0X+ D 1h1X+ D 2h2X= 0, 0h0T+ 1h1T + 2h2T + hDX +
D0+ m 2h0
h0X +
D1+ m 2h1
h1X+
D2+ m 2h2
h2X= 0, h DT+
D0+ m 2h0
h0T +
D1+ m 2h1
h1T +
D2+ m 2h2
h2T +
2h D+ 2m DX +
D20+1
2g(h1+ h2) + m h0
D
h0X
+
D21+1
2g(h0+ h2) + m h1
D
h1X+
D22+1
2g(h0+ h1) + m h2
D
h2X= 0,
h D+ m DT +
1 2
D2+ gI1
0+ m2 0+1
2g(h − h1− h2) + gh1h2h−1+ m 2h0D
h0T
+
1 2
D2+ gI1
1+ m2 1+1
2g(h − h0− h2) + gh0h2h−1+ m 2h1
D
h1T
+
1 2
D2+ gI1
2+ m2 2+1
2g(h − h0− h1) + gh0h1h−1+ m 2h2
D
h2T +
3
2h D2+1
2g I1h+ m2h−1+ 3m D
DX
+
1
2(D2+ gI1)D0+ m2D 0+1
2gh D+ gh1h2h−1D +3
4 m h0
D2+1 4gm I1
h0 +1 2gm
h0X +
1
2(D2+ gI1)D1+ m2D 1+1
2gh D+ gh0h2h−1D +3
4 m h1
D2+1 4gm I1
h1 +1 2gm
h1X
+
1
2(D2+ gI1)D2+ m2D 2+1
2gh D+ gh0h1h−1D +3
4 m h2D2+1
4gm I1
h2 +1 2gm
h2X = 0.
(22)
Or, in matrix form:
AUT + BUX = 0, where
U=
⎡
⎢⎢
⎣ D h0
h1
h2
⎤
⎥⎥
⎦ , A=
⎡
⎢⎢
⎣
0 a12 a13 a14
0 a22 a23 a24
a31 a32 a33 a34
a41 a42 a43 a44
⎤
⎥⎥
⎦ , B=
⎡
⎢⎢
⎣
b11 b12 b13 b14
b21 b22 b23 b24
b31 b32 b33 b34
b41 b42 b43 b44
⎤
⎥⎥
⎦ .
The coefficients of A are given by:
a11= 0, a12= 0, a13= 1, a14= 2, a21= 0, a22= 0, a23= 1, a24= 2, a31= h, a32= D0+ m
2h0, a33= D1+ m
2h1, a34= D2+ m 2h2, a41= hD + m,
a42= 1 2
D2+ gI1
0+ m2 0+1
2g(h − h1− h2) + gh1h2h−1+ m
2h0D, a43= 1
2
D2+ gI1
1+ m2 1+1
2g(h − h0− h2) + gh0h2h−1+ m
2h1D, a44= 1
2
D2+ gI1
2+ m2 2+1
2g(h − h0− h1) + gh0h1h−1+ m
2h2D.
The coefficients of B are given by:
b11= −L, b12 = D 0, b13= D 1, b14 = D 2, b21= h, b22= D0+ m
2h0, b23= D1+ m
2h1, b24= D2+ m 2h2, b31= 2hD + 2m,
b32= D20+1
2g(h1+ h2) + m h0D, b33= D21+1
2g(h0+ h2) + m h1D, b34= D22+1
2g(h0+ h1) + m h2D, b41= 3
2h D2+1
2g I1h+ m2h−1+ 3m D,
b42= 1
2(D2+ gI1)D0+ m2D 0+1
2gh D+ gh1h2h−1D+3 4
m h0
D2 +1
4gm I1
h0 +1 2gm, b43= 1
2(D2+ gI1)D1+ m2D 1+1 2gh D +gh0h2h−1D+3
4 m h1
D2+1 4gm I1
h1 +1 2gm, b44= 1
2(D2+ gI1)D2+ m2D 2+1
2gh D+ gh0h1h−1D +3
4 m h2
D2+1 4gm I1
h2 +1 2gm.
The eigenvalues are the roots of the fourth-order polynomial:
det(B − λA) = 0. (23)
To simplify the computations, let us recall that the introduction of the depth averaged velocity U = uh/h allows us to rewrite the mass conservation equation in standard form (16) and to establish the Galilean invariance of the modulation equations. There- fore, to check the hyperbolicity for any U is equivalent to check the hyperbolicity for U = 0. Since:
U =m h + D,
we will put into the coefficients of the matrices A and B the value of D corresponding to U = 0:
D= −m h.
Let us also remark that for the hyperbolicity study, one can always take h0= 1 in the coefficients of the polynomial (23) (if h0 = α > 0, all the eigenvalues will only be multiplied by√α). The corresponding symmetry relations are the consequences of the fact that the periodic solution is determined in terms of the third degree polynomial.
7 Hyperbolicity Region
Since the roots of polynomial (12) satisfy the inequality : 1 = h0 < h1 < h2, one can parameterize h1and h2as : h1= s, h2 = s + τ, where s > 1 and τ > 0. Such a parameterization allows us to use a standard ‘Cartesian’ frame for the computation of the eigenvalues of (23). The eigenvalues thus are given explicitly as functions of s andτ. We used Wolfram Mathematica for such a computation. The numerical results show that the eigenvalues are all real in a large region = {(s, τ)|1 < s < 100, 0 <
τ < 100}. First, for each pair (s, τ) from , we computed the numerical values of matrices A and B. Then, the corresponding eigenvalues were computed as the roots of the fourth-order polynomial (23). Moreover, one can find that the resultant of the
Fig. 1 The resultant R of the polynomial (23) and its derivative does not change sign in the region 1< s <
100 and 0< τ < 100 where all roots are real. The resultant is plotted here in a smaller region of s and τ.
Thus, the polynomial (23) has no multiple real roots. It means that the system of modulation equations is strictly hyperbolic
Fig. 2 The region s> 1, τ > 0 is divided by a smooth curve into two sub-regions, ‘grey’ and ‘white’.
If m< 0, in the grey sub-region, one has three positive and one negative eigenvalues, while in the white sub-region, one has two positive and two negative eigenvalues. If m> 0, the signs of the roots will only change, since we have taken vanishing U
Fig. 3 The case of right facing waves is considered (m< 0). The behavior of the eigenvalues is shown as functions of the parameter md = h2− h1
h2− h0 = τ
τ + s − 1along straight lines. On the left figure, positive eigenvaluesλ1> λ2> λ3> 0 are shown along the line s = 1.2−τ, τ ∈ (0, 0.2) which belongs entirely to the ‘grey’ region. Two of them become multiple (shown by black dots) only at the boundary of the domain:
in the limit of small amplitude (linear waves limit) or infinite length (solitary waves limit). On the right figure, the behavior of all eigenvalues is shown along the line s= 12 − τ, τ ∈ (0, 11). This line crosses the boundary between the ‘white’ and ‘grey’ regions at the point s≈ 3.6697 and τ ≈ 8.3303 corresponding to md≈ 0.757 (shown by ‘cross’) where the lowest positive eigenvalue changes the sign from positive to negative when mdincreases
corresponding polynomial (23) and its derivative (denoted further as R) does not change sign (see Fig.1). Thus, the polynomial (23) has no multiple real roots.
Since all the roots are real and different, the system of modulation equations for the SGN equations is thus strictly hyperbolic. The fact that periodic waves of all lengths are modulationally stable corroborates the results of [32] where the spectral stability of solitary waves (the limit s→ 1) has been proven for small amplitudes and numerically confirmed for large amplitudes. The whole hyperbolicity region s > 1, τ > 0 is divided by a smooth curve corresponding toλ = 0 into two sub-regions, ‘grey’ and
‘white’ (see Fig.2). If m< 0, in the grey sub-region, one has three positive and one negative eigenvalues, while in the white sub-region, one has two positive and two negative eigenvalues. If m> 0, the signs of the roots will only change, since we have taken U vanishing. Figure3shows the behavior of eigenvalues for m< 0. The roots become multiple only at the boundary of the domain corresponding to the linear waves and solitary waves.
8 Numerical Results
To study the modulational stability of periodic waves to SGN model numerically, we use the following set of parameters: h0= 1 m, h1= 1.5 m, h2= 2 m (i.e. s = 1.5 m andτ = 0.5 m), and g = 10 m/s2, as an example, for the solution of the water height described by (11) over a single wave length L which corresponds to 7.4163 m approximately. To setup the problem, a wave train is formed initially that consists of N aforementioned single stationary wave solutions (see Sect.3) in a domain of
Fig. 4 The initial conditions for the modulational test of periodic solutions of SGN model; water height h is shown on the left, and the phase portrait graphed in the(h, h ˙h) plane is shown on the right. Three different perturbation amplitudes,i.e., a= 10−3(first row), a= 10−2(second row), and a= 10−1(third row), are considered here with N= 50
size L1 = N × L, where the travelling wave speed of this wave train is taken be D= −m/h ≈ 3.1688 m/s (this allows us to take vanishing the averaged over period the mass weighted velocity U = hu/h). With that, we then introduce perturbations to the height of the wave train h(x) as:
Fig. 5 Numerical results for a modulational stability test of periodic waves shown for four time instants : t= 200, 400, 800, 1200 s. The graphs are displayed in the same manner as in Fig.4
˜h(x) := h(x)
1+ a cos
2πx L1
, (24a)
and define the perturbed velocity of the wave by
˜u(x) := m
˜h(x)+ D, (24b)
Fig. 5 continued
where a is a small parameter. For this problem, periodic boundary conditions are assumed and used on the left and right of the interval[0, L1].
In the numerical simulations of the SGN model performed below, we take N = 50, and perturbation amplitudes a = 10− j for j = 1, 2, 3 in the runs. We show the initial condition of the test in Fig.4, and the computed solutions at four different times t = 200, 400, 800, 1200 s in Fig.5, where both the water height and the phase portrait in the(h, h ˙h)-plane are present. The choice of h ˙h variable is natural, because for the
Fig. 5 continued
travelling wave solutions h ˙h= mdh
dξ, so up to a multiplicative constant,(h, h ˙h)-plane is nothing than the classical phase space
h,dhdξ
. The periodic wave remains stable when the smaller values a= 10−2and 10−3are used. To see the limit of linear stability, we have also taken a large-amplitude perturbation (a= 10−1). The periodic wave train becomes unstable: we are too far from the classical ‘small perturbation analysis’. The numerical results are obtained using a hyperbolic–elliptic splitting method proposed by the authors [21] with 400 meshes for each wave length.
Fig. 5 continued
9 Conclusions and Perspectives
We have derived Whitham equations to the SGN model and show that they are strictly hyperbolic for arbitrary wave amplitudes, i.e., the corresponding periodic wave trains are modulationally stable. This corroborates the results [32] where the linear stability of solitary waves (which can be considered as the limit of periodic waves of large length) has been proven. The existence of the Riemann invariants and study of characteristic fields (are they genuinely degenerate or genuinely nonlinear in the sense of Lax)
is a hopeless problem due to quite complex expressions of coefficients of Whitham equations. Numerical solutions of Riemann problem for modulation equations can probably clarify this structure.
Acknowledgements Part of this work was performed during the visit of SG to the Department of Mathe- matics of National Taiwan University. SG was partially supported by l’Agence Nationale de la Recherche, France (grant numbers ANR-11-LABEX-0092, and ANR-11-IDEX-0001-02). KMS was partially sup- ported by MOST 105-2115-M-002-013-MY2 and 107-2115-M-002-014.
A Multiscale Decomposition
The classical Whitham method [44] consists in decomposing the scales in the following way (for simplicity, we consider just the velocity variable u):
u(x, t) = u
θ(x, t), μx, μt
, θ(x, t) = (μx, μt)
μ . (25)
Here,θ is a fast phase variable, is a slow phase variable, and μ is a small parameter.
The solution is supposed to be 2π-periodic with respect to θ. The definitions of the local wave frequencyω and the local wave number κ:
∂θ
∂t = −ω, ∂θ
∂x = κ, (26)
automatically imply the evolution equation forκ:
κt+ ωx = 0. (27)
Written in slow variables X = μx, T = μt, Eq. (26) is equivalent to:
∂
∂T = −ω, ∂
∂ X = κ, and (27) reads as:
κT + ωX = 0. (28)
One can also define the travelling wave coordinateξ = x − Dt, and the phase velocity D= ω/κ. The solution is decomposed as:
u(x, t) = u(x − Dt, μx, μt) = uθ(x,t)
κ , μx, μt
= u
(X, T ) μκ , X, T
. (29)
The wavelength L is defined as:
L =2π κ .
Thus, sinceθ ∈ [0, 2π], then ξ = μκ ∈ 0,2κπ
= [0, L]. Therefore, u is an L- periodic function with respect to the travelling coordinateξ:
u(ξ + L, X, T ) = u(ξ, X, T ).
Finally, (28) reads:
(1/L)T + (D/L)X = 0. (30)
The approach employing the travelling coordinate is equivalent to the one using the phase variable. The consistency equation can always be written in any of the two forms presented above: (28) or (30).
B Computation of Elliptic Integrals Let:
h2> h1> h0, P3(h) = (h − h0)(h − h1)(h2− h). (31) Then, one has:
1 2
h2 h1
√h dh
P3(h) =
h2− h0E(k) +√h0K(k) h2− h0
, 1
2
h2
h1
√dh
P3(h) = K(k)
√h2− h0
, 1
2
h2
h1
h−1dh
√P3(h) = (n, k) h2
√h2− h0
.
Here, K(k), E(k), and (n, k) are the complete elliptic integrals of the first, second, and third type, respectively:
K(k) =
π 2
0
dθ 1− k2sin2θ,
E(k) =
π
2
0
1− k2sin2θdθ, (n, k) =
π
2
0
dθ (1 − n sin2θ)
1− k2sin2θ. The characteristic n and elliptic modulus k are defined as:
n= h2− h1
h2 , k2=h2− h1
h2− h0.
