1. Introduction
Despite numerous studies devoted to capillary–gravity waves, this topic still fascinates researchers. The review by Dias & Kharif (Reference Dias and Kharif1999) briefly summarises clearly what was known on this subject at the end of the twentieth century. Monographs by Okamoto & Shoji (Reference Okamoto and Shoji2001) and Vanden-Broeck (Reference Vanden-Broeck2010) are other exhaustive sources of information on various types of capillary–gravity waves.
Travelling capillary–gravity waves of permanent form have been most deeply understood in the framework of weakly nonlinear and weakly dispersive equations, such as the Korteweg–de Vries (KdV) equation and the extended Korteweg–de Vries equation with fifth-order derivatives (KdV5). These equations model the unidirectional propagation of long waves in shallow water with some weak capillarity. Despite the apparent simplicity of the KdV5 model, it possesses a rich family of solutions. One of them consists of the so-called generalised solitary waves. These are solitary-wave pulses that are homoclinic to small-amplitude oscillatory waves. The formation of these waves is mathematically justified by the presence of a resonance in the dynamical system corresponding to the travelling waves (Lombardi Reference Lombardi2000). The existence of generalised solitary waves can be deduced from the existence of a phase shift when the action of the stationary problem changes sign (Bridges & Donaldson Reference Bridges and Donaldson2005). One can also use the existence of Smale’s horseshoe dynamics on the zero energy set (Buffoni, Champneys & Toland Reference Buffoni, Champneys and Toland1996a ). Then, it is relatively straightforward to construct a symbolic orbit that represents a generalised solitary wave. By using the first method one obtains a continuum family of solutions, while the latter gives ‘only’ a countable set of orbits. Benilov, Grimshaw & Kuznetsova (Reference Benilov, Grimshaw and Kuznetsova1993) showed that, for sufficiently weak surface tension, generalised solitary waves cannot exist. For non-vanishing values of the surface tension Grimshaw & Joshi (Reference Grimshaw and Joshi1995) constructed a one-parameter family of generalised solitary waves for the KdV5 equation using the methods of exponential asymptotics; see also Yang & Akylas (Reference Yang and Akylas1996). The existence of multi-pulse solutions was shown numerically and analytically (Calvo & Akylas Reference Calvo and Akylas1997; Champneys & Groves Reference Champneys and Groves1997). The main difficulties in computing numerically the generalised solitary waves for the KdV5 equation are well described by Boyd (Reference Boyd2007). The stability of multi-pulse solitary waves was studied by Chardard, Dias & Bridges (Reference Chardard, Dias and Bridges2009) using the Maslov theory.
Much less is known for the full water wave problem. However, the existence of generalised capillary–gravity solitary waves for the full Euler equations was shown by Sun (Reference Sun1991) and Sun & Shen (Reference Sun and Shen1993). As highlighted by Beale (Reference Beale1991), the generalised solitary waves stem from a resonance with periodic waves of the same speed. By using the method of boundary integral equations and Newton iterations to solve the resulting discrete system, Hunter & Vanden-Broeck (Reference Hunter and Vanden-Broeck1983) computed solitary and periodic capillary–gravity waves in finite depth (see also Vanden-Broeck (Reference Vanden-Broeck2010) and the references therein) while multi-modal solitary waves of depression are generated in Dias, Menasce & Vanden-Broeck (Reference Dias, Menasce and Vanden-Broeck1996) in both finite and infinite depth. Recently, generalised solitary waves were computed for the full water wave problem (Champneys, Vanden-Broeck & Lord Reference Champneys, Vanden-Broeck and Lord2002). The main purpose of the present paper is to show that there is a plethora (likely an infinite number) of generalised solitary waves for the full Euler equations.
In the present study, we consider a formulation for steady capillary–gravity solitary waves by following the pioneering work of Babenko (Reference Babenko1987). This formulation is based on the conformal mapping technique (Ovsyannikov Reference Ovsyannikov1974; Dyachenko et al. Reference Dyachenko, Kuznetsov, Spector and Zakharov1996a ; Dyachenko, Zakharov & Kuznetsov Reference Dyachenko, Zakharov and Kuznetsov1996b ) that has been successfully used to numerically compute periodic gravity waves in deep water (Choi & Camassa Reference Choi and Camassa1999) and in finite depth (Longuet-Higgins & Fox Reference Longuet-Higgins and Fox1996; Li, Hyman & Choi Reference Li, Hyman and Choi2004). More recently, this approach was adapted also to periodic capillary–gravity waves in deep water (Milewski, Vanden-Broeck & Wang Reference Milewski, Vanden-Broeck and Wang2010). The advantage of Babenko’s formulation is that it does not add nonlinearities to the Euler equations in the conformally mapped domain. For instance, the Babenko equation being quadratic in nonlinearity for pure gravity waves, it is easily solved numerically for solitary waves (Clamond & Dutykh Reference Clamond and Dutykh2013). This approach is extended here to compute generalised solitary waves.
In order to solve the Babenko equation numerically for capillary–gravity solitary waves, it is discretised using the Fourier-type pseudo-spectral method (Boyd Reference Boyd2000). The resulting system of nonlinear equations is solved using the well-known Levenberg–Marquardt (LM) method (Levenberg Reference Levenberg1944; Marquardt Reference Marquardt1963). This algorithm represents a mixture between the steepest descent far from the solution and the classical Newton method in the vicinity of the solution (Moré Reference Moré and Watson1978). It has been shown to be a robust nonlinear solver even in problems with millions of unknowns (Lourakis & Argyros Reference Lourakis and Argyros2005).
The paper is organised as follows. In § 2, we present the main constitutive assumptions of the mathematical model, as well as the conformal mapping technique and the derivation of a Babenko-like equation. The numerical resolution procedure is explained in § 3. Numerical results are presented in § 4. The new results on multi-hump generalised solitary waves are discussed in § 5.
2. Mathematical model
We consider a steady two-dimensional potential flow induced by a solitary wave in a horizontal channel of constant depth. The fluid is assumed to be inviscid and homogeneous. The pressure is equal to the surface tension at the impermeable free surface and the fixed horizontal seabed is impermeable as well. The flow is driven by the volumetric gravity force (directed downward) and by the capillary forces at the free surface.
Let
$(x,y)$
be a Cartesian coordinate system moving with the wave,
$x$
being the horizontal coordinate and
$y$
being the upward vertical one. The equations
$y=-d$
,
$y={\it\eta}(x)$
and
$y=0$
denote the positions of the bottom of the channel, of the free surface and of the mean water level, respectively. The latter implies that the Eulerian average
$\langle \bullet \rangle$
of the free surface is zero:
$$\begin{eqnarray}\langle {\it\eta}\rangle \equiv \lim _{L\rightarrow \infty }\frac{1}{2\,L}\int _{-L}^{L}{\it\eta}(x)\,\text{d}x=0.\end{eqnarray}$$
This definition of the average operator is not suitable for the computation of solitary waves. However, for a classical solitary wave, the mean condition (2.1) implies that
$$\begin{eqnarray}{\it\eta}(x)\rightarrow 0\quad \text{as }x\rightarrow \pm \infty ,\end{eqnarray}$$
the latter being more tractable for computations (see § 3). Also, for classical solitary waves, with
$|{\it\eta}|$
decaying faster than
$1/|x|$
in the far field, the condition (2.1) implies that the wave mass
$\int _{-\infty }^{\infty }{\it\eta}(x)\,\text{d}x$
is finite. Here, we consider finite-mass generalised solitary waves with a
$(2\ell )$
-periodic tail in the far field. For such waves, the condition (2.1) yields
$$\begin{eqnarray}\frac{1}{2\ell }\int _{x-\ell }^{x+\ell }{\it\eta}({\it\xi})\,\text{d}{\it\xi}\rightarrow 0\quad \text{as }x\rightarrow \pm \infty .\end{eqnarray}$$
The definition (2.3) of the mean level is more suitable than (2.1) for practical computations of generalised solitary waves (see § 3). Obviously, the condition (2.3) is equivalent to (2.2) when the amplitude of the tail tends to zero.
Let
${\it\phi}$
,
${\it\psi}$
,
$u$
and
$v$
be, respectively, the velocity potential, the stream function, and the horizontal and vertical velocities, such that
$u\equiv \partial _{x}{\it\phi}=\partial _{y}{\it\psi}$
and
$v\equiv \partial _{y}{\it\phi}=-\partial _{x}{\it\psi}$
. It is convenient to introduce the complex potential
$f\equiv {\it\phi}+\text{i}{\it\psi}$
(with
$\text{i}^{2}=-1$
) and the complex velocity
$w\equiv u-\text{i}v$
that are holomorphic functions of
$z\equiv x+\text{i}y$
(
$w=\text{d}f/\text{d}z$
). The complex conjugate is denoted with an asterisk (e.g.
$z^{\ast }=x-\text{i}y$
), while subscripts ‘
$b$
’ denote the quantities evaluated at the seabed – e.g.
$z_{b}(x)=x-\text{i}d$
,
${\it\phi}_{b}(x)={\it\phi}(x,y=-d)$
– and subscripts ‘
$s$
’ denote the quantities evaluated at the surface – e.g.
$z_{s}(x)=x+\text{i}{\it\eta}(x)$
,
${\it\phi}_{s}(x)={\it\phi}(x,y={\it\eta}(x))$
. Note that, e.g.,
$u_{s}=(\partial _{x}{\it\phi})_{s}\neq \partial _{x}{\it\phi}_{s}=u_{s}+{\it\eta}_{x}v_{s}$
. The traces of
${\it\psi}$
on the upper and lower boundaries,
${\it\psi}_{s}$
and
${\it\psi}_{b}$
, are constants because the surface and the bottom are streamlines. These constants are related to the mean flow by
$$\begin{eqnarray}-\!c\equiv \frac{1}{d}\left\langle \int _{-d}^{{\it\eta}}u\,\text{d}y\right\rangle =\frac{{\it\psi}_{s}-{\it\psi}_{b}}{d},\end{eqnarray}$$
so
$c$
is the wave phase velocity observed in the frame of reference without mean flow (
$c>0$
if the wave travels towards the increasing
$x$
-direction, so we consider
$c>0$
without loss of generality).
The dynamic condition can be expressed in term of the Bernoulli equation
$$\begin{eqnarray}2p+2gy+u^{2}+v^{2}=B,\end{eqnarray}$$
where
$p$
is the pressure divided by the density,
$g>0$
is the acceleration due to gravity and
$B$
is a Bernoulli constant. At the free surface
$y={\it\eta}(x)$
the pressure
$p$
reduces to the effect of the surface tension, that is
$p_{s}=-{\it\tau}\partial _{x}[{\it\eta}_{x}(1+{\it\eta}_{x}^{2})^{-1/2}]$
,
${\it\tau}$
being a (constant) surface tension coefficient (divided by the density). Averaging (2.5) written at the free surface, the definition of the mean level (2.1) yields an equation for the Bernoulli constant
$B$
,
$$\begin{eqnarray}B=\langle u_{s}^{2}+v_{s}^{2}\rangle .\end{eqnarray}$$
In the far field of a classical solitary wave, the free surface is horizontal and the flow is a uniform current. Thus,
${\it\eta}\rightarrow 0$
,
$u\rightarrow -c$
and
$v\rightarrow 0$
as
$x\rightarrow \pm \infty$
, where
$c$
is defined in (2.4). It follows that
$B=c^{2}$
for classical solitary waves. For the generalised solitary waves considered here, the far field involves a periodic wave of constant amplitude and it is therefore not a uniform current, implying that
$B\neq c^{2}$
in general. By analogy with the classical solitary wave, we define another phase velocity by
$$\begin{eqnarray}c^{\prime }\equiv \sqrt{B}.\end{eqnarray}$$
Many other phase velocities could of course be defined, such as the mean horizontal velocity at the bottom in the far field, but these considerations are secondary for the purpose of the present work. Note that the definition (2.7) implies that
$c^{\prime }>0$
without loss of generality.
2.1. Conformal mapping
Let us change to an independent variable:
$z\mapsto {\it\zeta}\equiv (\text{i}{\it\psi}_{s}-f)/c$
, which conformally maps the fluid domain
$$\begin{eqnarray}-\!\infty \leqslant x\leqslant \infty ,\quad -d\leqslant y\leqslant {\it\eta}(x)\end{eqnarray}$$
into the strip
$$\begin{eqnarray}-\!\infty \leqslant {\it\alpha}\leqslant \infty ,\quad -d\leqslant {\it\beta}\leqslant 0,\end{eqnarray}$$
where
${\it\alpha}\equiv \text{Re}({\it\zeta})$
and
${\it\beta}\equiv \text{Im}({\it\zeta})$
. The conformal mapping yields the Cauchy–Riemann relations
$x_{{\it\alpha}}=y_{{\it\beta}}$
,
$x_{{\it\beta}}=-y_{{\it\alpha}}$
, thence we have
$$\begin{eqnarray}\displaystyle \frac{c}{w}=-\frac{\text{d}z}{\text{d}{\it\zeta}},\quad \frac{u}{c}=\frac{-x_{{\it\alpha}}}{x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2}},\quad \frac{v}{c}=\frac{-y_{{\it\alpha}}}{x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2}},\quad \frac{u^{2}+v^{2}}{c^{2}}=\frac{1}{x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2}},\quad \quad & & \displaystyle\end{eqnarray}$$
while the pressure at the free surface can be conveniently written (see appendix A) as
$$\begin{eqnarray}\frac{p_{s}}{{\it\tau}}=-\left(\frac{\text{d}z_{s}^{\ast }}{\text{d}{\it\alpha}}\right)^{-1}\frac{\text{d}}{\text{d}{\it\alpha}}\left[\frac{\text{i}x_{{\it\alpha}}+y_{{\it\alpha}}}{\sqrt{x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2}}}\right]_{{\it\beta}=0}.\end{eqnarray}$$
With the change of dependent variables
$$\begin{eqnarray}\displaystyle x={\it\alpha}+X({\it\alpha},{\it\beta}),\quad y={\it\beta}+Y({\it\alpha},{\it\beta}), & & \displaystyle\end{eqnarray}$$
the Cauchy–Riemann relations
$X_{{\it\alpha}}=Y_{{\it\beta}}$
and
$X_{{\it\beta}}=-Y_{{\it\alpha}}$
hold, while the bottom (
${\it\beta}=-d$
) and the free surface (
${\it\beta}=0$
) impermeabilities yield
$$\begin{eqnarray}\displaystyle Y_{b}({\it\alpha})\equiv Y({\it\alpha},-d)=0,\quad Y_{s}({\it\alpha})\equiv Y({\it\alpha},0)={\it\eta}({\it\alpha}). & & \displaystyle\end{eqnarray}$$
The functions
$X$
and
$Y$
can be expressed in terms of
$X_{b}$
as (Clamond Reference Clamond1999, Reference Clamond2003)
$$\begin{eqnarray}\displaystyle & \displaystyle X({\it\alpha},{\it\beta})={\textstyle \frac{1}{2}}X_{b}({\it\zeta}+\text{i}d)+{\textstyle \frac{1}{2}}X_{b}({\it\zeta}^{\ast }-\text{i}d)=\text{cos}\,[({\it\beta}+d)\partial _{{\it\alpha}}]X_{b}({\it\alpha}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & Y({\it\alpha},{\it\beta})=\displaystyle {\displaystyle \frac{1}{2\text{i}}}X_{b}({\it\zeta}+\text{i}d)-\displaystyle {\displaystyle \frac{1}{2\text{i}}}X_{b}({\it\zeta}^{\ast }-\text{i}d)=\text{sin}\,[({\it\beta}+d)\partial _{{\it\alpha}}]X_{b}({\it\alpha}), & \displaystyle\end{eqnarray}$$
where an asterisk denotes the complex conjugate. Thus, the Cauchy–Riemann relations and the bottom impermeability are fulfilled identically. At the free surface
${\it\beta}=0$
, (2.14) yields
$$\begin{eqnarray}X_{s}({\it\alpha})=\text{cos}\,[d\partial _{{\it\alpha}}]X_{b}({\it\alpha})\quad \Longleftrightarrow \quad X_{b}({\it\alpha})=\text{sec}\,[d\partial _{{\it\alpha}}]X_{s}({\it\alpha}),\end{eqnarray}$$
and hence the relation (2.15) yields
$$\begin{eqnarray}Y_{s}({\it\alpha})=\text{tan}\,[d\partial _{{\it\alpha}}]X_{s}({\it\alpha})\equiv \mathscr{T}\{X_{s}\},\end{eqnarray}$$
which relates quantities written at the free surface only. The relation (2.17) can be trivially inverted giving, in particular,
$$\begin{eqnarray}(\partial _{{\it\alpha}}X)_{s}=\partial _{{\it\alpha}}\text{cot}\,[d\partial _{{\it\alpha}}]Y_{s}\equiv \mathscr{C}\{Y_{s}\}=\mathscr{C}\{{\it\eta}\},\end{eqnarray}$$
where
$\mathscr{T}$
and
$\mathscr{C}$
are pseudo-differential operators that, for a pure frequency, take the form
$$\begin{eqnarray}\displaystyle \mathscr{T}\{\text{e}^{\text{i}k{\it\alpha}}\}=\text{i}\tanh (kd)\,\text{e}^{\text{i}k{\it\alpha}},\quad \mathscr{C}\{\text{e}^{\text{i}k{\it\alpha}}\}=\left\{\begin{array}{@{}lr@{}}k\coth (kd)\text{e}^{\text{i}k{\it\alpha}} & (k\neq 0),\\ 1/d & (k=0).\end{array}\right. & & \displaystyle\end{eqnarray}$$
2.2. Babenko-like equation
Using (2.10a
), (2.11) and the relation
$u^{2}+v^{2}=ww^{\ast }$
, the Bernoulli equation (2.5) at the free surface can be written as
$$\begin{eqnarray}\displaystyle w_{s} & = & \displaystyle \displaystyle \frac{B-2g{\it\eta}-2p_{s}}{w_{s}^{\ast }}=\frac{2p_{s}+2g{\it\eta}-B}{c}\frac{\text{d}z_{s}^{\ast }}{\text{d}{\it\alpha}}\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{2g{\it\eta}-B}{c}\frac{\text{d}(x_{s}-\text{i}{\it\eta})}{\text{d}{\it\alpha}}-\frac{2{\it\tau}}{c}\frac{\text{d}}{\text{d}{\it\alpha}}\left[\frac{\text{i}x_{{\it\alpha}}+y_{{\it\alpha}}}{\sqrt{x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2}}}\right]_{{\it\beta}=0}.\end{eqnarray}$$
With
$w({\it\zeta})=u({\it\alpha},{\it\beta})-\text{i}v({\it\alpha},{\it\beta})$
being a holomorphic function such that
$\text{Im}(w)=0$
at the bottom, we have at the free surface – see above the derivation of (2.17) –
$$\begin{eqnarray}-\!v_{s}({\it\alpha})=\mathscr{T}\{u_{s}\}=\text{tan}\,[d\partial _{{\it\alpha}}]u_{s}({\it\alpha}).\end{eqnarray}$$
Applying the operator
$\text{cot}\,[d\partial _{{\it\alpha}}]=\partial _{{\it\alpha}}^{-1}\mathscr{C}$
to (2.21), then exploiting the relations (2.12a,b
), (2.18) and (2.20), after some algebra we obtain a pseudo-differential equation for
${\it\eta}$
:
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \mathscr{C}\left\{B{\it\eta}-\frac{g{\it\eta}^{2}}{2}+{\it\tau}-\frac{{\it\tau}(1+\mathscr{C}\{{\it\eta}\})}{\sqrt{(1+\mathscr{C}\{{\it\eta}\})^{2}+{\it\eta}_{{\it\alpha}}^{2}}}\right\}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad -\,g{\it\eta}(1+\mathscr{C}\{{\it\eta}\})+\frac{\text{d}}{\text{d}{\it\alpha}}\left\{\frac{{\it\tau}{\it\eta}_{{\it\alpha}}}{\sqrt{(1+\mathscr{C}\{{\it\eta}\})^{2}+{\it\eta}_{{\it\alpha}}^{2}}}\right\}=K,\end{eqnarray}$$
where
$K$
is an integration constant to be determined from the mean level condition (2.1).
$K=0$
for classical solitary waves (since
${\it\eta}(\infty )=0$
) but for generalised solitary waves
$K$
is not necessarily zero. However, via a change of definition of the mean water level (and therefore of
$d$
and
$B$
), it is always possible to take
$K=0$
without loss of generality. This is particularly convenient for a numerical resolution of (2.22), the mean level condition being enforced a posteriori as explained in the numerical procedure below.
The equation (2.22) is a Babenko (Reference Babenko1987) equation modified to incorporate capillarity (Buffoni, Dancer & Toland Reference Buffoni, Dancer and Toland2000). One advantage of Babenko-like formulations, compared to other techniques based on conformal mapping (Choi & Camassa Reference Choi and Camassa1999; Li et al. Reference Li, Hyman and Choi2004; Milewski et al. Reference Milewski, Vanden-Broeck and Wang2010), consists in reducing the degree of nonlinearity of the original Euler equations. This feature is advantageous, in particular, for numerical resolutions.
3. Numerical method
The generalised Babenko equation (2.22) is solved numerically using a Fourier collocation technique (Boyd Reference Boyd2000; Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2006) and a Levenberg–Marquardt–Gauß–Newton method for solving the nonlinear algebraic system of equations resulting from the discretisation (Nocedal & Wright Reference Nocedal and Wright2006). These numerical procedures being classical, they are only briefly described here.
3.1. Discretisation
The conformal domain
$-\infty <{\it\alpha}<\infty$
is truncated on a finite interval
$-{\it\Lambda}<{\it\alpha}\leqslant {\it\Lambda}$
and discretised at
$N$
equally spaced nodes
${\it\alpha}_{j}=-{\it\Lambda}+j{\rm\Delta}{\it\alpha}$
, with
${\rm\Delta}{\it\alpha}=2{\it\Lambda}/N$
and
$j=1,\ldots ,N$
, the values
${\it\eta}_{j}={\it\eta}({\it\alpha}_{j})$
being obtained from the discrete equation (2.22).
The generalised solitary waves obtained here can be seen as the nonlinear superposition of a classical solitary wave (possibly with multiple humps) and a periodic wave (with a single fundamental frequency). The former is thus called here the ‘solitary component’ of the wave, while the latter is called the ‘periodic component’.
${\it\Lambda}$
is then chosen large enough so that the ‘solitary component’ of the wave is zero to machine precision for
$|{\it\alpha}|>{\it\Lambda}_{0}$
with
$0<{\it\Lambda}_{0}<{\it\Lambda}$
.
${\it\Lambda}-{\it\Lambda}_{0}$
is also chosen large enough so that there are several periods of the wave ‘periodic component’.
The domain
${\it\alpha}\in \, ]-{\it\Lambda},{\it\Lambda}]$
is considered periodic in order to apply the fast Fourier transform (FFT).
$N$
is chosen large enough so that the Fourier spectrum is zero to machine precision for all the highest frequencies. Thus, the spatial and spectral truncation errors are extremely small (around machine precision). In practice, all the (pseudo) differential operations are performed in Fourier space, while all the nonlinear operations are performed in the (discrete)
${\it\alpha}$
-space. With this approach, we obtain a nonlinear system of algebraic equations that must be solved numerically with an iterative procedure.
Note that the
$(2{\it\Lambda})$
-periodic conformal
${\it\alpha}$
-domain corresponds to a
$(2L)$
-periodic physical
$x$
-domain with, in general,
$L\neq {\it\Lambda}$
. Note also that the computational nodes are not equally spaced in the physical domain. Note finally that in the far field of a generalised solitary wave, the free surface
${\it\eta}(x)$
is
$(2\ell )$
-periodic while
${\it\eta}({\it\alpha})$
is
$(2{\it\lambda})$
-periodic with, in general,
$\ell \neq {\it\lambda}$
.
3.2. Nonlinear equations solver
For pure gravity waves (
${\it\tau}=0$
), the equations can be solved efficiently with the simple Petviashvili scheme (Pelinovsky & Stepanyants Reference Pelinovsky and Stepanyants2004; Clamond & Dutykh Reference Clamond and Dutykh2013). When
${\it\tau}\neq 0$
, Petviashvili’s scheme fails due to the non-homogeneous character of the nonlinear terms involved in the capillary terms. Several variants of the Newton method are possible alternatives to the Petviashvili scheme. The reasons for the failure of some candidates (Boyd Reference Boyd2007) were taken into account in choosing a suitable method. The best results were obtained by employing the so-called Levenberg–Marquardt algorithm (Nocedal & Wright Reference Nocedal and Wright2006). We use the Matlab
® version of this algorithm implemented in the fsolve function of the Optimisation toolbox. Some computations were checked and confirmed with Nielsen (Reference Nielsen1999) and Lourakis (Reference Lourakis2004) implementations of this method.
For the computations herein below, the physical parameters are
$g=d=1$
, various
$B$
and
${\it\tau}$
being chosen for each computation. The Levenberg–Marquardt initial relaxation parameter is set to
$0.05$
and the tolerance is
$10^{-12}$
, the computations being performed in double precision (about 16 digits). Most of the time, the half of the computational domain
${\it\Lambda}/d_{0}$
varies between
$20$
and
$40$
. The number of Fourier modes is generally
$N=1024$
, but a larger number of nodes is sometimes used, for example with the solutions computed in figure 2.
3.3. Initial guess
In order to obtain multi-hump (classical and generalised) solitary waves, the initial guesses of the iterations contain the desired number of pulses. Depending on the initial guess, the algorithm may converge or not to a non-trivial solution. The solution thus obtained may not have the same number of pulses as the initial guess.
All the initial guesses we use satisfy the far-field condition
${\it\eta}=0$
to machine precision for
${\it\Lambda}_{0}<|{\it\alpha}|\leqslant {\it\Lambda}$
. Thus, the mean level condition (2.1) is numerically fulfilled according to (2.2) and (2.3), and
$K=0$
for the initial guesses. For all the computations, we choose unit gravity (
$g=1$
) and unit initial mean depth (
$d_{0}=1$
), while various surface tensions
${\it\tau}$
and initial Bernoulli constants
$B_{0}=c_{0}^{\prime 2}$
are considered.
In most of the cases, our implementation of the method makes use of the natural numerical continuation in the surface tension
${\it\tau}$
, as a way to improve the performance and to accelerate the convergence. Thus, for a given
$B_{0}$
, we start the algorithm with a small
${\it\tau}$
and subsequently increase
${\it\tau}$
little by little to reach the desired value.
During the computations, the mean level condition (2.1) is not enforced at each iteration, as it would greatly complicate the algorithm. It is simpler to release this constraint and carry out the computations until convergence is achieved. Eventually, the physical mean depth and Bernoulli constant are computed according to their definitions given in § 2. The procedure is detailed below.
3.4. Post-processing
An important aspect of the resolution method is the numerical definition of the averaging operations. In order to fulfil the condition (2.1) accurately, the computational domain has to be very large: this is beyond most computer capabilities. Such a large domain, however, is not necessary for computations. Indeed, if a classical solitary wave (and the ‘solitary part’ of a generalised solitary wave) is decaying rapidly, then there exists
${\it\Lambda}_{0}$
such that the solitary-wave contribution is smaller than the machine precision for all
$|{\it\alpha}|\geqslant {\it\Lambda}_{0}$
(the far field), so it does not contribute to the numerical solution in the far field. Therefore, the computational domain needs only to be large enough so that the ‘solitary part’ of the wave is zero to machine precision for
${\it\Lambda}_{0}\leqslant |{\it\alpha}|\leqslant {\it\Lambda}$
. Hence, in the far field, the flow is numerically
$2\ell$
-periodic to machine precision and can be expanded in Fourier polynomials as
$$\begin{eqnarray}\displaystyle {\it\eta}(x)\approx \mathop{\sum }_{m=-M}^{M}\mathfrak{h}_{m}\,\text{e}^{\text{i}m{\rm\pi}x/\ell },\quad |w_{s}(x)|^{2}\approx \mathop{\sum }_{m=-M}^{M}\mathfrak{q}_{m}\,\text{e}^{\text{i}m{\rm\pi}x/\ell }, & & \displaystyle\end{eqnarray}$$
where the unknown parameters
$\ell$
,
$\mathfrak{h}_{m}$
and
$\mathfrak{q}_{m}$
are all determined by nonlinear least-squares minimisation (Clamond & Barthélémy Reference Clamond and Barthélémy1995). According to (2.3), the mean water level
$\bar{{\it\eta}}$
is then
$$\begin{eqnarray}\bar{{\it\eta}}=\frac{1}{2\ell }\int _{-\ell }^{\ell }{\it\eta}(x)\,\text{d}x=\mathfrak{h}_{0}.\end{eqnarray}$$
In principle,
$\bar{{\it\eta}}$
should be zero to the accuracy involved in the numerical procedure. If this is not the case, it means that the depth
$d_{0}$
introduced a priori in the resolution procedure has ‘drifted’ and should then be redefined (renormalised) in order to compute the a posteriori (actual) depth
$d_{\infty }$
defined by
$$\begin{eqnarray}d_{\infty }\equiv d_{0}+\bar{{\it\eta}}=d_{0}+\mathfrak{h}_{0}.\end{eqnarray}$$
Similarly, the a posteriori Bernoulli constant
$B_{\infty }$
is then defined by
$$\begin{eqnarray}B_{\infty }\equiv \frac{1}{2\ell }\int _{-\ell }^{\ell }[2g{\it\eta}+|w_{s}(x)|^{2}]\,\text{d}x=2g\mathfrak{h}_{0}+\mathfrak{q}_{0},\end{eqnarray}$$
and the a posteriori phase velocity
$c_{\infty }^{\prime }$
is thus
$$\begin{eqnarray}c_{\infty }^{\prime }\equiv \sqrt{B_{\infty }-2g\bar{{\it\eta}}}=\sqrt{\mathfrak{q}_{0}}.\end{eqnarray}$$
In practice, the choice
$M=5$
is more than sufficient because the tails of the generalised solitary waves computed here are quasi-sinusoidal (see below).
3.5. Dimensionless parameters
In order to characterise the waves, we introduce the dimensionless parameters
$$\begin{eqnarray}\mathit{Fr}^{2}\equiv \frac{c^{\prime 2}}{gd},\quad \mathit{Bo}\equiv \frac{{\it\tau}}{gd^{2}},\end{eqnarray}$$
$\mathit{Fr}$
and
$\mathit{Bo}$
being, respectively, the Froude and Bond numbers. (Our definition of the Bond number is the most common in the theory of water waves, but it is the reciprocal of the Bond number usually defined in fluid mechanics.)
The Froude and Bond numbers are chosen a priori to compute the solutions, i.e. the procedure described above implies initial Froude and Bond numbers
$\mathit{Fr}_{0}$
and
$\mathit{Bo}_{0}$
, respectively. Once the computations have converged, we compute the a posteriori Froude and Bond numbers
$\mathit{Fr}_{\infty }$
and
$\mathit{Bo}_{\infty }$
, respectively.
4. Numerical results
We demonstrate here the efficiency of the algorithm by computing several known and new solutions.
Figure 1. Examples of solitary waves of depression for
$\mathit{Fr}=0.75$
and
$\mathit{Bo}=0.4$
.
4.1. Classical solitary waves
For supercritical Bond numbers
$\mathit{Bo}>1/3$
and suitable values of the Froude number
$\mathit{Fr}<1$
, the existence of isolated solitary waves of depression is known, theoretically and numerically, see Champneys et al. (Reference Champneys, Vanden-Broeck and Lord2002) and the references therein. A first check of the code is the computation of one of these waves corresponding to the values
$\mathit{Bo}=0.4$
and
$\mathit{Fr}=0.75$
(figure 1
a). The initial guess of the procedure was chosen simply as a negative localised bump. The algorithm’s convergence does not seem to be sensitive to the choice of the initial guess for such simple solutions.
Another type of solution, computed by the code and shown here in figure 1, consists of multi-modal solitary waves. The existence of these solutions in the KdV5 equation was shown theoretically and numerically (Champneys & Groves Reference Champneys and Groves1997). Investigating the neighbourhood of
$\mathit{Fr}=1$
and
$\mathit{Bo}=1/3$
, Buffoni, Groves & Toland (Reference Buffoni, Groves and Toland1996b
) proved, for the Euler equations, the existence of an infinite number of solitary waves. Using a boundary integral equation technique and a different numerical procedure (Hunter & Vanden-Broeck Reference Hunter and Vanden-Broeck1983), multi-modal solitary waves of depression were computed by Dias et al. (Reference Dias, Menasce and Vanden-Broeck1996).
In order to obtain multi-modal solitary waves of depression, the initial guess was chosen with several negative bumps. Three examples of multi-modal solitary waves of depression are given in figure 1 for
$\mathit{Bo}=0.4$
and
$\mathit{Fr}=0.75$
(identical to the uni-modal example). This shows numerically the non-unicity of the solution for a fixed pair of parameters (
$\mathit{Fr},\mathit{Bo}$
). Actually, solutions with any number of pulses could probably be obtained. For the same number of pulses, the distance between the pulses is also not unique (figure 1
b,c). Our numerical investigations suggest that the distance between the two pulses cannot be reduced below a certain value. However, the distance can be increased at will. In figure 1(b), the two pulses interact because
${\it\eta}(0)/d=-3.39\times 10^{-3}$
is larger than the numerical accuracy of the computations, this case being the shortest distance we were able to compute.
All the examples in figure 1 have identical Froude and Bond numbers, but also identical amplitudes
$\min ({\it\eta}/d)=-0.42396\pm 1.3\times 10^{-5}$
. For these examples, the algorithm converged to the desired precision (
$10^{-12}$
), the different solutions being obtained by just changing the initial guess. These solutions were first calculated with computational boxes of length
${\it\Lambda}/d=25$
and
$N=1024$
nodes, and subsequently verified with
${\it\Lambda}/d=50$
and
$N=2048$
.
Figure 2. Examples of multi-modal solitary waves with damped oscillations. (a)
$Fr=0.9$
,
$Bo=0.3$
; (b)
$Fr=0.99$
,
$Bo=0.3$
; (c)
$Fr=0.9945$
,
$Bo=0.3$
.
For
$\mathit{Fr}<1$
and
$\mathit{Bo}<1/3$
, solitary waves with tails oscillating around the rest level are obtained. A variety of multi-hump solitary waves with damped oscillatory tails is illustrated in figure 2 where some ‘exotic’ waves are formed as the Froude number approaches 1 (similar solutions have been computed by Dias et al. (Reference Dias, Menasce and Vanden-Broeck1996)). These examples illustrate the abilities of our numerical procedure. (For instance, the computation of figure 2(c) requires the limit of our hardware capabilities, but more powerful computers will permit the computation of more complicated solutions.)
4.2. Generalised solitary waves
The question of the existence of classical solitary waves of elevation has not been solved, to our knowledge, although some references suggest a negative answer (Champneys et al.
Reference Champneys, Vanden-Broeck and Lord2002). For
$\mathit{Bo}<1/3$
and
$\mathit{Fr}>1$
, what is known is the existence of generalised solitary waves (solitary waves that are homoclinic to oscillatory waves of exponentially small amplitude). Here, we focus on the numerical generation of waves of this kind.
Our experiments concern the range
$\mathit{Bo}<1/3$
and
$\mathit{Fr}>1$
, illustrating thus the influence of the capillarity (smaller or larger values of
$\mathit{Bo}$
) on the computations. The resulting numerical profiles are shown in figure 3 for the same a priori Froude number
$\mathit{Fr}_{0}=1.1$
and various Bond numbers. In all these cases, periodic tails emerged. The results presented here were obtained by incrementing the Bond number
$\mathit{Bo}$
. As a result, starting from small values of the Bond number close to
$0.1$
, the numerical results suggest the emergence of generalised solitary waves as a consequence of a resonance between solitary and periodic waves of the same speed (Beale Reference Beale1991).
Figure 3. Examples of generalised solitary waves with
$\mathit{Fr}_{0}=1.1$
. (a)
$\mathit{Bo}_{0}=0.1$
,
$\mathit{Bo}_{\infty }=0.0999998$
,
$\mathit{Fr}_{\infty }=1.099999$
; (b) zoom of panel (a); (c)
$\mathit{Bo}_{0}=0.15$
,
$\mathit{Bo}_{\infty }=0.1499$
,
$\mathit{Fr}_{\infty }=1.099$
; (d)
$\mathit{Bo}_{0}=0.25$
,
$\mathit{Bo}_{\infty }=0.248$
,
$\mathit{Fr}_{\infty }=1.095$
.
Figure 4. Fourier spectrum of figure 3(c).
Figure 5. Influence of the computational box for
$\mathit{Fr}_{0}=1.1$
and
$\mathit{Bo}_{0}=0.25$
. Thick blue line:
${\it\Lambda}/d_{0}=20$
,
$d_{\infty }/d_{0}=1.00914$
,
$\mathit{Fr}_{\infty }=1.0873$
,
$\mathit{Bo}_{\infty }=0.24549$
. (a) Thin red line:
${\it\Lambda}/d_{0}=30.11$
,
$d_{\infty }/d_{0}=1.00913$
,
$\mathit{Fr}_{\infty }=1.0874$
,
$\mathit{Bo}_{\infty }=0.2455$
; (b) thin red line:
${\it\Lambda}/d_{0}=30.5$
,
$d_{\infty }/d_{0}=1.0095$
,
$\mathit{Fr}_{\infty }=1.0868$
,
$\mathit{Bo}_{\infty }=0.24531$
.
For
$\mathit{Bo}=0.1$
, the tail is not visible in the wave of figure 3(a), but it is nevertheless present (see the zoom in figure 3
b). The amplitude of the ripples grows with the Bond number. These tails are quasi-sinusoidal, as can be seen from the Fourier transform of the free surface (figure 4), showing that the fundamental frequency is dominant with small and exponentially decaying harmonics.
With the emergence of an oscillating tail, the mean surface level is slightly moved, requiring a renormalisation of the physical parameters. However, this drift is small in all the solutions we computed, and at most a few per cent in the worst cases (figures 3 and 5).
The examples in figure 3 were obtained by continuation for the Bond number, increasing it progressively from
$\mathit{Bo}=0$
. If instead the solution is computed from the initial guess with a finite Bond number, a different solution may be obtained. This is illustrated for
$\mathit{Fr}_{0}=1.1$
and
$\mathit{Bo}_{0}=0.25$
, see figures 3(d) and 5.
For classical solitary waves, the length of the computational domain has no influence, provided that it is long enough. This is not the case for generalised solitary waves, where the tails are somehow quantised by the length of the domain. Nonlinear interactions between the solitary and periodic parts of the wave induce phase shifts, so the computational domain is generally not an exact multiple of tail wavelength (
$L/\ell$
and
${\it\Lambda}/{\it\lambda}$
are not integers, in general). Figure 5 illustrates the influence of the computational box on the results. The solution shown by the thick blue line was computed with
${\it\Lambda}=20d_{0}$
. Recomputing the solution in a longer domain with
${\it\Lambda}=30.11d_{0}$
, we obtained the same solution (see the thin red line in figure 5
a). This shows that, though the computational domain is periodic, the computed solution is an accurate approximation of an aperiodic wave. However, for domains of arbitrary lengths, different solutions are generally obtained, as shown for
${\it\Lambda}=30.5d_{0}$
in figure 5(b). This example clearly illustrates the influence of the computational domain on the tail characteristics that can be obtained with the algorithm described in this paper. Whether or not the wavelength and amplitude of the tail can vary continuously for exact generalised solitary waves is an open question that we cannot answer here.
Figure 6. Multi-hump generalised solitary waves with
$\mathit{Fr}_{0}=1.17$
and
$\mathit{Bo}_{0}=0.12$
. (a)
$\mathit{Fr}_{\infty }=1.2466$
,
$\mathit{Bo}_{\infty }=0.1181$
,
$d_{\infty }/d_{0}=1.008$
,
$\max ({\it\eta}/d_{\infty })=0.4018$
. (b)
$\mathit{Fr}_{\infty }=1.2396$
,
$\mathit{Bo}_{\infty }=0.1183$
,
$d_{\infty }/d_{0}=1.007$
,
$\max ({\it\eta}/d_{\infty })=0.4021$
.
Figure 7. Unequal multi-hump generalised solitary waves with
$\mathit{Fr}_{0}=1.1$
and
$\mathit{Bo}_{0}=0.25$
. (a)
$\mathit{Fr}_{\infty }=1.0876$
,
$\mathit{Bo}_{\infty }=0.24559$
; (b)
$\mathit{Fr}_{\infty }=1.0875$
,
$\mathit{Bo}_{\infty }=0.24557$
; (c)
$\mathit{Fr}_{\infty }=1.0871$
,
$\mathit{Bo}_{\infty }=0.24541$
; (d)
$\mathit{Fr}_{\infty }=1.0876$
,
$\mathit{Bo}_{\infty }=0.24558$
.
4.3. Multi-hump generalised solitary waves
A final type of waves computed and shown here consists of multi-hump generalised solitary waves. The existence of homoclinic connections with several loops near a resonance for a family of Hamiltonian systems has been recently analysed by Jézéquel, Bernard & Lombardi (Reference Jézéquel, Bernard and Lombardi2014). In the context of water waves, these solutions would correspond to multi-hump generalised solitary waves and we offer here some numerical evidence of their existence.
Taking a finite number of separated bumps as the initial iteration with
$\mathit{Fr}=1.17$
and
$\mathit{Bo}=0.12$
, a generalised two-pulse capillary–gravity solitary wave is shown in figure 6(a) while, for the same values of the parameters, a three-pulse solution is depicted in figure 6(b). For these two waves, the a posteriori parameters differ by less than one per cent, suggesting the non-unicity of the solution. This non-unicity is to be expected since it occurs already for classical solitary waves.
Our experiments suggest that the process of adding new pulses can be continued indefinitely, provided that the length and resolution of the computational domain are gradually increased as well. This experiment provides numerical evidence that the solutions should not be unique for a given set of Froude and Bond numbers, the number of solutions being likely infinite (countable infinity in the number of pulses, and a continuum of solutions with varying distance between pulses).
For a pair of classical solitary waves, we found that the two humps cannot be too close. This seems not to be the case with multi-hump generalised solitary waves because we succeeded in computed largely overlapping waves, such as in figure 7. These solutions were all obtained for
$\mathit{Fr}_{0}=1.1$
,
$\mathit{Bo}_{0}=0.25$
and
${\it\Lambda}=30d_{0}$
(though some results are plotted on restricted domains for better display). Except for figure 7(c), the other, very different solutions have the same (a posteriori) Froude and Bond numbers within some degree of accuracy (their discrepancies are smaller than
$0.01\,\%$
). This is strong numerical evidence of the non-unicity of generalised solitary waves.
5. Discussion
In this study, the problem of steady capillary–gravity solitary waves was reformulated on a flat domain and the corresponding Babenko-like equation was derived. This equation was discretised with a Fourier-type pseudo-spectral method. The resulting discrete nonlinear and non-local equation was solved using the Levenberg–Marquardt method. The accuracy and capabilities of the method were demonstrated via several examples.
Using this formulation we succeeded in computing, directly or by continuation in the Bond number
$\mathit{Bo}$
, classical and generalised solitary waves. Above the critical Bond number
$\mathit{Bo}=1/3$
, we found the classical localised solitary waves of depression that propagate with subcritical speeds
$\mathit{Fr}<1$
, in agreement with the predictions of the KdV5 model. We also showed that various (generalised) multi-pulse solitary waves of elevation, but also of depression (localised), can be successfully computed by our method. To our knowledge, these generalised multi-pulse solitary waves have never been computed in the context of capillary–gravity surface waves. The numerical simulations suggest the existence of an infinite number of such waves for equal Froude and Bond numbers.
The stability of these multi-hump generalised solitary waves is a question of physical and theoretical interest. Periodic (Stokes) surface waves are known to be unstable (McLean Reference McLean1982). The generalised solitary waves involve periodic waves in their far field, so they are likely unstable as well, but rigorous investigations have yet to be performed.
Numerical evidence is of course not rigorous mathematical proof, so mathematical investigations would be of great interest. Indeed, the lack of theoretical results for these multi-hump generalised solitary waves makes their physical analysis difficult. Such theoretical insights should also be beneficial for improving the computation of these waves; it would then be easier to compute other solutions, such as asymmetric generalised solitary waves (if they exist). We hope that the numerical evidence provided here will stimulate such theoretical investigations.
Acknowledgements
D.C. and D.D. acknowledge the support of the CNRS under the PEPS Inphyniti 2015 project FARA. A.D. has been supported by the project MTM2014-54710-P. The authors would like to thank Professor T. Lakoba (University of Vermont, USA) for stimulating discussions.
Appendix A. Details of calculation
With the conformal mapping of non-overturning waves (
$x_{{\it\alpha}}\geqslant 0$
), we have
$$\begin{eqnarray}{\it\eta}_{x}(1+{\it\eta}_{x}^{2})^{-1/2}=[y_{{\it\alpha}}(x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2})^{-1/2}]_{{\it\beta}=0},\end{eqnarray}$$
thence the horizontal derivative of (A 1) yields
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}x}\left[\frac{{\it\eta}_{x}}{\sqrt{1+{\it\eta}_{x}^{2}}}\right] & \hspace{-2.5pt}=\hspace{-2.5pt} & \displaystyle \displaystyle \left[\frac{x_{{\it\alpha}}y_{{\it\alpha}{\it\alpha}}-y_{{\it\alpha}}x_{{\it\alpha}{\it\alpha}}}{(x_{{\it\alpha}}^{2}+y_{{\it\alpha}}^{2})^{3/2}}\right]_{{\it\beta}=0}=\left[\frac{z_{{\it\alpha}}^{\ast }z_{{\it\alpha}{\it\alpha}}-z_{{\it\alpha}}z_{{\it\alpha}{\it\alpha}}^{\ast }}{2\text{i}(z_{{\it\alpha}}z_{{\it\alpha}}^{\ast })^{3/2}}\right]_{{\it\beta}=0}\nonumber\\ \displaystyle & \hspace{-2.5pt}=\hspace{-2.5pt} & \displaystyle \displaystyle \text{i}\left(\frac{\text{d}z_{s}^{\ast }}{\text{d}{\it\alpha}}\right)^{-1}\frac{\text{d}}{\text{d}{\it\alpha}}\left[\frac{z_{{\it\alpha}}^{\ast }}{z_{{\it\alpha}}}\right]_{{\it\beta}=0}^{1/2}=\left(\frac{\text{d}z_{s}^{\ast }}{\text{d}{\it\alpha}}\right)^{-1}\frac{\text{d}}{\text{d}{\it\alpha}}\left[\frac{\text{i}z_{{\it\alpha}}^{\ast }}{\sqrt{z_{{\it\alpha}}z_{{\it\alpha}}^{\ast }}}\right]_{{\it\beta}=0},\end{eqnarray}$$
from which the relation (2.11) is derived at once.
