1 Introduction
Many physical phenomena and mathematical problems related to surface gravity water waves remain unknown, not well understood or unsolved, even in the ‘simple’ case of travelling waves of permanent form in two-dimensional irrotational motion (Clamond Reference Clamond2012; Constantin Reference Constantin2012). Travelling waves are of special interest because complex sea states are often described as the superposition and interaction of such waves. Surveys of analytical and numerical models, their limitations and open questions can be found in dedicated articles (Fenton Reference Fenton1988, Reference Fenton1999; Dias & Kharif Reference Dias and Kharif1999; Groves Reference Groves2004; Strauss Reference Strauss2010) and books (Okamoto & Shōji Reference Okamoto and Shōji2001; Vanden-Broeck Reference Vanden-Broeck2010; Constantin Reference Constantin2011).
Since exact analytic solutions for irrotational steady surface gravity waves are still unknown, and probably will never be known, only analytic or numerical approximations are accessible. Simple analytic approximations are interesting for physical insights, but they are of limited accuracy. Even formal analytic solutions in terms of small parameter expansions have limited accuracy because they generally converge slowly (Schwartz Reference Schwartz1974) when they converge (shallow water expansions are divergent for all amplitudes (Germain Reference Germain1967)). Thus, their numerical calculation suffers large truncation errors and is prone to important accumulation of round-off and cancellation errors. Even when a simple analytic approximation would be sufficient for a given application, its computation may be practically intractable. An example is the Korteweg–de Vries (KdV) cnoidal wave analytic solution which cannot be easily computed for very long waves, as demonstrated in the present paper. Therefore, only numerical approximations of the original equations can provide highly accurate solutions which are necessary for many applications. For instance, for stability analysis and interactions of travelling waves using accurate numerical models, too crude (approximately six digits accuracy, say) an initial condition for the travelling waves may lead to incorrect behaviour, especially for long-time simulations. A lack of initial accuracy may then lead to erroneous physical interpretations. Another example is the numerical investigation of mathematical conjectures. With arbitrary precision computations, one can check whether a conjecture is likely to be true or not, or can formulate new conjectures worthy of investigation. Indeed, some questions of mathematical interest (e.g. rigorous proofs of unicity and monotonicity regarding velocity, acceleration and pressure fields) remain open for all waves and not only for very large amplitudes (see the references in Constantin (Reference Constantin2012) for reviews).
Several algorithms have been proposed in the literature for the computation of steady surface wave solutions of the irrotational Euler equations (Okamoto & Shōji Reference Okamoto and Shōji2001). The focus has been more on the computation of the almost highest gravity waves (see, e.g., Williams Reference Williams1985; Byatt-Smith Reference Byatt-Smith2001; Maklakov Reference Maklakov2002) or ‘exotic’ capillary–gravity waves (e.g. Vanden-Broeck Reference Vanden-Broeck2010; Clamond, Dutykh & Durán Reference Clamond, Dutykh and Durán2015) than the computation of arbitrary wavelengths. Actually, none of these algorithms are capable of computing long waves in shallow water, not even the ones of small amplitude, because they are too demanding. Indeed, all of these algorithms lead to the resolution of a discrete system of nonlinear equations, this system being large for long waves. The resolution of this system is generally performed with (the like of) Newton or, better, Levenberg–Marquardt iterations (Lourakis & Argyros Reference Lourakis and Argyros2005). Although robust and effective, these methods are computationally very demanding because each iteration requires
$O(N^{3})$
operations,
$N$
being the number of unknowns (e.g. the Fourier coefficient for periodic water waves). Indeed, although their theoretical minimum complexity is
$O(N^{2})$
, usually Levenberg–Marquardt iterations use direct solvers that have
$O(N^{3})$
complexity. In the simplest case, they rely on Cholesky decomposition which needs
$O(N^{3})$
operations. In the most numerically robust version of the Levenberg–Marquardt method, its implementation relies on a rank-revealing QR algorithm, which is also
$O(N^{3})$
(with a large constant in front of
$N^{3}$
). Thus, when
$N$
is very large, the computational time may be prohibitive and the accumulation of round-off errors significant, even if an algorithm of complexity
$O(N^{2})$
can be used. It is well known that for steep waves and cnoidal waves in shallow water, the number of Fourier modes needed for accurate resolution is very large, especially when using conformal mapping. Therefore, an algorithm with complexity lower than
$O(N^{2})$
is desirable to achieve these computations. It is the purpose of the present paper to describe such an algorithm with overall complexity
$O(N\log N)$
, which is thus suitable for arbitrary precision computations of all waves of practical interest in arbitrary depth.
For a simple wave equation, Petviashvili (Reference Petviashvili1976) proposed an algorithm based on stabilised fixed-point iterations for computing solitary waves. This algorithm is interesting because it is very simple to implement and has complexity
$O(N)$
. However, Petviashvili’s method generally works only for equations with special nonlinear terms, typically autonomous equations with homogeneous nonlinearities (Yang Reference Yang2010). Several variants have then been proposed to somehow extend the scope of these modified Petviashvili methods (Ablowitz & Musslimani Reference Ablowitz and Musslimani2005; Lakoba & Yang Reference Lakoba and Yang2007; Álvarez & Duràn Reference Álvarez and Duràn2014a
). Instead of tweaking the algorithm, our approach here is to rewrite the Euler equations in a form suitable for their numerical resolution via the classical Petviashvili method. By doing so, we obtain an algorithm suitable for the computation of waves in arbitrary depth, in particular long waves with wavelengths of tenths to millions times the mean water depth. To the best of our knowledge, it is the first algorithm capable of computing such long waves.
The Petviashvili method works for the Euler equations if they are rewritten in the form of a Babenko equation (Babenko Reference Babenko1987). This was successfully implemented to compute solitary gravity waves (Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014). Unfortunately, this algorithm does not work for periodic waves. In infinite depth, periodic waves were successfully computed with a modified Petviashvili method by Dyachenko, Lushnikov & Korotkevich (Reference Dyachenko, Lushnikov and Korotkevich2013). Their algorithm does not work in finite depth, however. In order to overcome these drawbacks, we propose here a simple change of variable that transforms the Babenko equation into a form tractable with the classical Petviashvili method. Our algorithm is not designed for the most extreme waves: for a given wavelength, it works only for all wave heights of less than approximately
$99\,\%$
of the maximum one. Thus, we are able to rapidly compute waves of practical interest in arbitrary depth (infinite, finite and shallow) and to arbitrary precision.
This paper is organised as follows. In § 2, we present the physical assumptions and the mathematical definitions and notations. In § 3, we introduce a conformal mapping and we give precise definitions of all the variables and parameters in the conformal space. In § 4, we derive a Babenko equation written in a suitable form for a fast numerical resolution. The numerical algorithm is subsequently described in § 5 and relevant numerical examples are provided in § 6. A summary and perspectives are outlined in § 8.
2 Definitions and notations
We consider steady two-dimensional potential flows due to surface gravity waves in constant depth
$d$
. The fluid is of (positive) constant density
$\unicode[STIX]{x1D70C}$
, the pressure is zero at the impermeable free surface, while the seabed is fixed, horizontal and impermeable.
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 wave is (
$2\unicode[STIX]{x03C0}/k$
)-periodic and
$x=0$
is the abscissa of a wave crest. (The fundamental wavenumber
$k$
is zero for solitary and aperiodic waves.) By
$y=-d$
,
$y=\unicode[STIX]{x1D702}(x)$
and
$y=0$
we denote respectively the equations of the bottom, the free surface and the mean water level (see figure 1
a). The latter implies that
$\langle \unicode[STIX]{x1D702}\rangle =0$
(where
$\langle \bullet \rangle$
represents the Eulerian average operator over one spatial period (wavelength)), i.e.
$$\begin{eqnarray}\langle \unicode[STIX]{x1D702}\rangle \stackrel{def}{=}\frac{k}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}/k}^{\unicode[STIX]{x03C0}/k}\unicode[STIX]{x1D702}(x)\,\text{d}x=0.\end{eqnarray}$$
Here,
$a\stackrel{def}{=}\unicode[STIX]{x1D702}(0)$
denotes the wave crest amplitude and
$b\stackrel{def}{=}-\unicode[STIX]{x1D702}(\unicode[STIX]{x03C0}/k)$
is the wave trough amplitude, so that
$H\stackrel{def}{=}a+b$
is the total wave height. A wave steepness
$\unicode[STIX]{x1D700}$
is then classically defined as
$\unicode[STIX]{x1D700}\stackrel{def}{=}kH/2$
. With
$\unicode[STIX]{x1D719}$
denoting the velocity potential, the classical equations of motion are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{xx}+\unicode[STIX]{x1D719}_{yy}=0\quad \text{for }-d\leqslant y\leqslant \unicode[STIX]{x1D702}(x), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{y}=0\quad \text{at }y=-d, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{y}-\unicode[STIX]{x1D702}_{x}\unicode[STIX]{x1D719}_{x}=0\quad \text{at }y=\unicode[STIX]{x1D702}(x), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 2g\unicode[STIX]{x1D702}+\unicode[STIX]{x1D719}_{x}^{2}+\unicode[STIX]{x1D719}_{y}^{2}=B\quad \text{at }y=\unicode[STIX]{x1D702}(x), & \displaystyle\end{eqnarray}$$
where
$g>0$
is the (constant) acceleration due to gravity and
$B$
is the Bernoulli constant.
Definition sketch: (a)
$z$
-plane; (b)
$\unicode[STIX]{x1D701}$
-plane.
Let
$\unicode[STIX]{x1D713}$
,
$u$
and
$v$
be the stream function and the horizontal and vertical velocities respectively, such that
$u=\unicode[STIX]{x1D719}_{x}=\unicode[STIX]{x1D713}_{y}$
and
$v=\unicode[STIX]{x1D719}_{y}=-\unicode[STIX]{x1D713}_{x}$
. It is convenient to introduce the complex potential
$f\stackrel{def}{=}\unicode[STIX]{x1D719}+\text{i}\unicode[STIX]{x1D713}$
(with
$\text{i}^{2}=-1$
) and the complex velocity
$w\stackrel{def}{=}u-\text{i}v$
which are holomorphic functions of
$z\stackrel{def}{=}x+\text{i}y$
(i.e.
$f=f(z)$
and
$w=\text{d}f/\text{d}z$
). The complex conjugate is denoted with a star (e.g.
$z^{\ast }=x-\text{i}y$
), while overbars denote the quantities written at the seabed, e.g.
$\bar{z}(x)=x-\text{i}d$
,
$\bar{\unicode[STIX]{x1D719}}(x)=\unicode[STIX]{x1D719}(x,y=-d)$
, and overtildes denote the quantities written at the surface, e.g.
$\tilde{z}(x)=x+\text{i}\unicode[STIX]{x1D702}(x)$
,
$\tilde{\unicode[STIX]{x1D719}}(x)=\unicode[STIX]{x1D719}(x,y=\unicode[STIX]{x1D702}(x))$
. (It should be noted that, e.g.,
$\tilde{u} =\widetilde{~\unicode[STIX]{x1D719}_{x}\,}\neq \tilde{\unicode[STIX]{x1D719}}_{x}=\tilde{u} +\unicode[STIX]{x1D702}_{x}\tilde{v}$
.) As the free surface and the bottom are streamlines,
$\tilde{\unicode[STIX]{x1D713}}$
and
$\bar{\unicode[STIX]{x1D713}}$
are constants. One can then take
$\tilde{\unicode[STIX]{x1D713}}=0$
or
$\bar{\unicode[STIX]{x1D713}}=0$
without loss of generality (gauge condition for the stream function).
The pressure field can be obtained from the Bernoulli equation,
$$\begin{eqnarray}2p+2gy+u^{2}+v^{2}=B,\end{eqnarray}$$
where
$p$
is the pressure divided by the density. At the free surface, as the pressure is zero (i.e.
$\tilde{p}=0$
), the Bernoulli constant
$B$
is defined by averaging (2.6) applied at the free surface and using the condition (2.1), i.e.
$$\begin{eqnarray}B=\langle \tilde{u} ^{2}+\tilde{v}^{2}\rangle .\end{eqnarray}$$
From the incompressibility and the irrotationality, it follows that
$B$
can also be obtained from the expression at the bottom (Clamond & Constantin Reference Clamond and Constantin2013),
$$\begin{eqnarray}B=\langle \bar{u}^{2}\rangle .\end{eqnarray}$$
Then, from the Bernoulli equation averaged at the bed,
$\langle \bar{p}\rangle =gd$
. More generally,
$B$
equals
$u^{2}+v^{2}$
averaged along any streamline (in the frame of reference moving with the wave).
Let
$-c_{S}$
be the mean flow velocity, defined as
$$\begin{eqnarray}c_{S}\stackrel{def}{=}-\left\langle \frac{1}{d}\int _{-d}^{\unicode[STIX]{x1D702}}u(x,y)\,\text{d}y\right\rangle =\frac{\bar{\unicode[STIX]{x1D713}}-\tilde{\unicode[STIX]{x1D713}}}{d}.\end{eqnarray}$$
Thus,
$c_{S}$
is the phase velocity of the wave observed in the frame of reference without mean flow (Stokes’ second definition of phase celerity), which is also the frame where the wave impulse is zero (appendix B). Another important quantity is the phase velocity
$c_{E}$
observed in the frame of reference without mean velocity at the seabed (which is also the one where the circulation is zero; cf. appendix B),
$$\begin{eqnarray}c_{E}\stackrel{def}{=}-\langle \bar{u}\rangle =-\langle u(x,y=-d)\rangle .\end{eqnarray}$$
This is Stokes’ first definition of phase celerity. Since the motion is irrotational,
$c_{E}$
can be obtained by averaging
$u$
along any horizontal line
$y=\text{const.}$
(but not along the wavy surface, i.e.
$c_{E}\neq -\langle \tilde{u} \rangle$
). Many other phase velocities can of course be defined, but
$c_{S}$
and
$c_{E}$
are two velocities of special interest here. It should be noticed that
$B=c_{S}^{2}=c_{E}^{2}$
in deep water and for solitary waves (see below), and that neither
$c_{S}$
nor
$c_{E}$
is the linear phase velocity
$c_{0}\stackrel{def}{=}\sqrt{(g/k)\tanh (kd)}$
if the wave amplitude is not zero. A discussion on the Bernoulli constant and phase velocities can be found in Clamond (Reference Clamond2017a
). Further information regarding these phase velocities and integral quantities is given in appendices A–C.
3 Conformal mapping
Let us make the change of independent complex variable
$z\mapsto \unicode[STIX]{x1D701}\stackrel{def}{=}(\text{i}\tilde{\unicode[STIX]{x1D713}}-f)/c_{R}$
,
$c_{R}\neq 0$
being a velocity of reference. In practice, one could take
$c_{R}=c_{S}$
,
$c_{R}=c_{E}$
,
$c_{R}=c_{0}$
,
$c_{R}=B^{1/2}$
or
$c_{R}=(c_{E}c_{S})^{1/2}$
, for example, but another convenient choice can be made depending on the problem at hand. One should take
$c_{R}>0$
if the wave travels towards the increasing
$x$
direction in a ‘fixed’ frame of reference and
$c_{R}<0$
is the wave travels towards the decreasing
$x$
direction. Here, without loss of generality, we consider only waves travelling towards the increasing
$x$
direction.
This change of variable conformally maps the fluid fundamental domain
$$\begin{eqnarray}0\leqslant x\leqslant 2\unicode[STIX]{x03C0}/k,\quad -d\leqslant y\leqslant \unicode[STIX]{x1D702}(x)\end{eqnarray}$$
into the rectangle (see figure 1 b)
$$\begin{eqnarray}0\leqslant \unicode[STIX]{x1D6FC}\leqslant 2\unicode[STIX]{x03C0}c_{E}/kc_{R},\quad -dc_{S}/c_{R}\leqslant \unicode[STIX]{x1D6FD}\leqslant 0,\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}\stackrel{def}{=}\operatorname{Re}(\unicode[STIX]{x1D701})$
and
$\unicode[STIX]{x1D6FD}\stackrel{def}{=}\operatorname{Im}(\unicode[STIX]{x1D701})$
. For convenience, we introduce the apparent wavenumber 
and apparent depth 
in the conformal plane,
which are generally different from the corresponding quantities in the physical plane. It should be noted that
$c_{R}$
does not appear in the expression of 
, so no choice of
$c_{R}$
can enforce the equality 
Conversely, the choice
$c_{R}=\sqrt{c_{E}c_{S}}$
yields 
, so, with this peculiar choice of
$c_{R}$
, the areas of the fundamental periods are identical in physical and conformal planes. For our numerical resolution, we found it convenient to take
$c_{R}=c_{E}$
(see § 5 below).
Since
$\text{d}z/\text{d}\unicode[STIX]{x1D701}=z_{\unicode[STIX]{x1D6FC}}=-\text{i}z_{\unicode[STIX]{x1D6FD}}$
(with subscripts denoting partial derivatives), we have the Cauchy–Riemann relations
$x_{\unicode[STIX]{x1D6FC}}=y_{\unicode[STIX]{x1D6FD}}$
and
$x_{\unicode[STIX]{x1D6FD}}=-y_{\unicode[STIX]{x1D6FC}}$
, while the complex velocity and the velocity components are
$$\begin{eqnarray}\frac{w}{c_{R}}=-\left(\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right)^{-1},\quad \frac{u}{c_{R}}=\frac{-x_{\unicode[STIX]{x1D6FC}}}{x_{\unicode[STIX]{x1D6FC}}^{2}+y_{\unicode[STIX]{x1D6FC}}^{2}},\quad \frac{v}{c_{R}}=\frac{-y_{\unicode[STIX]{x1D6FC}}}{x_{\unicode[STIX]{x1D6FC}}^{2}+y_{\unicode[STIX]{x1D6FC}}^{2}},\quad \frac{u^{2}+v^{2}}{c_{R}^{2}}=\frac{1}{x_{\unicode[STIX]{x1D6FC}}^{2}+y_{\unicode[STIX]{x1D6FC}}^{2}};\end{eqnarray}$$
thence,
$$\begin{eqnarray}x_{\unicode[STIX]{x1D6FC}}=y_{\unicode[STIX]{x1D6FD}}=\frac{-c_{R}u}{u^{2}+v^{2}},\quad y_{\unicode[STIX]{x1D6FC}}=-x_{\unicode[STIX]{x1D6FD}}=\frac{-c_{R}v}{u^{2}+v^{2}}.\end{eqnarray}$$
With these relations, one can compute all of the physical quantities of interest.
3.1 Conformal averaging operator
The Bernoulli constant, from the relation (2.7), is defined in the conformal plane by
the integral being computed keeping
$\unicode[STIX]{x1D6FD}$
constant, in particular at the free surface
$\unicode[STIX]{x1D6FD}=0$
. The relation (3.6) shows that it is convenient to introduce the average operator over one period in the
$\unicode[STIX]{x1D6FC}$
variable (
$\unicode[STIX]{x1D6FD}$
being kept constant),
for any quantity
$(\cdots \,)$
. Thus, at the free surface and at the bottom, we have respectively
and conversely
As averaged physical quantities are defined in the physical plane and not in the conformal plane, the connections between the averaging operators are useful to express these physical quantities in the conformal plane. Conversely, these connections are also useful to express averaged quantities in the conformal plane (needed for a numerical resolution) in their physical plane counterparts. For easy reference, we give several such relations in appendix C.
3.2 Resolution of the conformal mapping
With the change of dependent variables
the Cauchy–Riemann relations
$X_{\unicode[STIX]{x1D6FC}}=Y_{\unicode[STIX]{x1D6FD}}$
and
$X_{\unicode[STIX]{x1D6FD}}=-Y_{\unicode[STIX]{x1D6FC}}$
hold, while the bottom (
) and the free surface (
$\unicode[STIX]{x1D6FD}=0$
) impermeabilities yield
At the boundaries of the fundamental period (i.e.
$\unicode[STIX]{x1D6FC}=0$
and 
), we have from (3.12a
)
and more generally 
. Therefore, the function
$X$
is 
-periodic. The functions
$X$
and
$Y$
can be expressed in terms of
$\bar{X}$
– i.e. the function
$X$
written at the bottom – as (see Clamond (Reference Clamond1999, Reference Clamond2003) for detailed derivations)
where a star denotes the complex conjugate. Thus, the Cauchy–Riemann relations and the bottom impermeability are fulfilled identically. At the free surface
$\unicode[STIX]{x1D6FD}=0$
, (3.15) yields 
, which can be inverted as 
. The relation (3.16) can then be rewritten with quantities expressed at the free surface only, i.e.
where
$\mathscr{H}$
is an anti-adjoint pseudo-differential operator acting on a pure frequency as
Indeed, as the differential operator
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}$
corresponds to a simple multiplication by
$\text{i}\unicode[STIX]{x1D705}$
in Fourier space (where
$\unicode[STIX]{x1D705}$
is the frequency), the pseudo-differential operators are defined in Fourier space by substituting
$\text{i}\unicode[STIX]{x1D705}$
for
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}$
. (It should be noted that
$\unicode[STIX]{x1D705}=nk$
for the
$n$
th Fourier mode of a (
$2\unicode[STIX]{x03C0}/k$
)-periodic function.) Thus, for instance, the operator 
is in Fourier space a multiplication by 
. It should be noted that the formal Taylor expansions, such as in (3.15) and (3.16), of the pseudo-differential operators are practically useless, especially in deep water (
$d\rightarrow \infty$
) where the operator
$\mathscr{H}$
becomes the classical Hilbert transform.
3.3 Averaging of dependent functions
At the bottom, by averaging
$x$
over one wavelength, one obtains easily
$$\begin{eqnarray}\langle \langle \bar{x}\rangle \rangle =c_{R}c_{E}^{-1}\langle \langle \unicode[STIX]{x1D6FC}\rangle \rangle +\langle \langle \bar{X}\rangle \rangle =\unicode[STIX]{x03C0}k^{-1}+\langle \langle \bar{X}\rangle \rangle .\end{eqnarray}$$
From the definition of the average operator (3.7), we have also
$$\begin{eqnarray}\langle \langle \bar{x}\rangle \rangle =\langle -c_{E}^{-1}x\bar{\unicode[STIX]{x1D719}}_{x}\rangle =\unicode[STIX]{x03C0}k^{-1}.\end{eqnarray}$$
By comparing (3.19) and (3.20), we obtain
$$\begin{eqnarray}\langle \langle \bar{X}\rangle \rangle =0,\end{eqnarray}$$
meaning that
$\bar{X}$
is a periodic function that averages zero. Thus, with the relations (3.17) to (3.21), the boundedness and periodicity of
$\bar{X}$
imply that
$$\begin{eqnarray}\langle \langle X(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})\rangle \rangle =0,\quad \langle \langle Y(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})\rangle \rangle =0.\end{eqnarray}$$
Hence, the functions
$X$
and
$Y$
have zero average in the
$\unicode[STIX]{x1D6FC}$
variable. The relation (3.17) can thus be inverted without ambiguities giving, in particular,
where
$\mathscr{C}$
is a self-adjoint positive-definite pseudo-differential operator such that
In deep water
$d\rightarrow \infty$
, we have
$c_{E}=c_{S}\stackrel{def}{=}c$
,
$B=c^{2}$
, and we thus take
$c_{R}=c$
because it is an obvious choice. Then, as
$d\rightarrow \infty$
, the operator
$\mathscr{C}$
tends to
$\mathscr{C}_{\infty }$
, with
$$\begin{eqnarray}\mathscr{C}_{\infty }\{\text{e}^{\text{i}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FC}}\}=|\unicode[STIX]{x1D705}|\text{e}^{\text{i}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FC}},\quad \mathscr{C}_{\infty }^{-1}\{\text{e}^{\text{i}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FC}}\}=\left\{\begin{array}{@{}ll@{}}|\unicode[STIX]{x1D705}|^{-1}\text{e}^{\text{i}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FC}} & (\unicode[STIX]{x1D705}\neq 0),\\ 0 & (\unicode[STIX]{x1D705}=0).\end{array}\right.\end{eqnarray}$$
It should be noted that
$\mathscr{C}_{\infty }^{-1}$
is singular for
$\unicode[STIX]{x1D705}=0$
, but the special choice
$\mathscr{C}_{\infty }^{-1}\{1\}\stackrel{def}{=}0$
does not matter as long as
$\mathscr{C}_{\infty }^{-1}$
is applied to a function averaging to zero (or to any known value that can be enforced).
In summary, we have obtained the special relations
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{x}_{\unicode[STIX]{x1D6FC}}=c_{R}c_{E}^{-1}+\mathscr{C}\{{\tilde{Y}}\}=c_{R}c_{S}^{-1}+\mathscr{C}\{{\tilde{y}}\}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{C}\{{\tilde{Y}}\}=\mathscr{C}\{{\tilde{y}}\}+c_{R}c_{S}^{-1}-c_{R}c_{E}^{-1} & \displaystyle\end{eqnarray}$$
and the averaged quantities
$$\begin{eqnarray}\langle \langle \mathscr{C}\{{\tilde{Y}}\}\rangle \rangle =0,\quad \langle \langle {\tilde{y}}\rangle \rangle =(c_{S}c_{E}^{-1}-1)d,\quad \langle \langle \mathscr{C}\{{\tilde{y}}\}\rangle \rangle =c_{R}c_{E}^{-1}-c_{R}c_{S}^{-1}.\end{eqnarray}$$
3.4 Mean level condition
The definition of the mean level in the transformed domain remains to be considered. Using the relations (3.12) and (3.28), the mean level condition (2.1) becomes
$$\begin{eqnarray}0=\langle \langle {\tilde{y}}\tilde{x}_{\unicode[STIX]{x1D6FC}}\rangle \rangle =\langle \langle ({\tilde{Y}}-d(1-c_{S}/c_{E}))(\tilde{X}_{\unicode[STIX]{x1D6FC}}+c_{R}/c_{E})\rangle \rangle =\langle \langle {\tilde{Y}}\tilde{X}_{\unicode[STIX]{x1D6FC}}\rangle \rangle -d(1-c_{S}/c_{E})c_{R}/c_{E};\end{eqnarray}$$
thence,
$$\begin{eqnarray}\langle \langle {\tilde{Y}}\mathscr{C}\{{\tilde{Y}}\}\rangle \rangle =d(c_{E}-c_{S})c_{R}c_{E}^{-2}=-c_{R}c_{E}^{-1}\langle \langle {\tilde{y}}\rangle \rangle .\end{eqnarray}$$
This relation is fundamental, in particular for a numerical resolution. With the previous results, several useful averaged relations can be expressed in terms of the important parameter
$\langle \langle {\tilde{y}}\rangle \rangle$
(cf. appendix C).
3.5 Celerities
The definition (3.6) of the Bernoulli constant written at the free surface and at the bottom yields
$$\begin{eqnarray}\frac{B}{c_{R}c_{E}}=\left\langle \left\langle \frac{c_{R}c_{S}^{-1}+\mathscr{C}\{{\tilde{y}}\}}{(c_{R}c_{S}^{-1}+\mathscr{C}\{{\tilde{y}}\})^{2}+{\tilde{y}}_{\unicode[STIX]{x1D6FC}}^{2}}\right\rangle \right\rangle =\left\langle \left\langle \frac{1}{c_{R}c_{S}^{-1}+\mathscr{S}\{{\tilde{y}}\}}\right\rangle \right\rangle ,\end{eqnarray}$$
the second equality deriving from the relations
where 
is a pseudo-differential operator acting on a pure frequency as
It should be noted that, as
$|\unicode[STIX]{x1D705}|\rightarrow \infty$
,
$\mathscr{S}$
decays exponentially fast in Fourier space, unlike
$\mathscr{C}$
which grows linearly. It should be noted also that
$\mathscr{S}\{{\tilde{y}}\}\rightarrow 0$
as
$d\rightarrow \infty$
; hence,
$B=c_{E}c_{S}$
in deep water (together with
$c_{E}=c_{S}=c$
, as mentioned above) and thus
$B=c^{2}$
.
4 Babenko equations
Using (3.4a
) and the relation
$u^{2}+v^{2}=ww^{\ast }$
, the Bernoulli equation (2.6) at the free surface can be written as
$$\begin{eqnarray}\tilde{w}=\frac{B-2g\unicode[STIX]{x1D702}}{\tilde{w}^{\ast }}=\frac{2g{\tilde{y}}-B}{c_{R}}\frac{\text{d}\tilde{z}^{\ast }}{\text{d}\unicode[STIX]{x1D6FC}}=\frac{(2g{\tilde{y}}-B)\tilde{x}_{\unicode[STIX]{x1D6FC}}}{c_{R}}-\text{i}\frac{(2g{\tilde{y}}-B){\tilde{y}}_{\unicode[STIX]{x1D6FC}}}{c_{R}}.\end{eqnarray}$$
As
$w=u-\text{i}v$
is a holomorphic function such that
$\operatorname{Im}(w)=0$
at the bottom, we have at the free surface (see the derivation of (3.17) in the previous section)
thence, with (3.26) and (4.1),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}(B{\tilde{y}}-g{\tilde{y}}^{2}) & = & \displaystyle \mathscr{H}\{(2g{\tilde{y}}-B)(c_{R}c_{S}^{-1}+\mathscr{C}\{{\tilde{y}}\})\}\nonumber\\ \displaystyle & = & \displaystyle 2gc_{R}c_{S}^{-1}\mathscr{H}\{{\tilde{y}}\}+2g\mathscr{H}\{{\tilde{y}}\mathscr{C}\{{\tilde{y}}\}\}-B{\tilde{y}}_{\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$
After simplifications and application of the antiderivative operator
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}^{-1}$
, one obtains at once the Babenko equation
$$\begin{eqnarray}Bg^{-1}{\tilde{y}}-{\textstyle \frac{1}{2}}{\tilde{y}}^{2}+K_{1}=c_{R}c_{S}^{-1}\mathscr{C}^{-1}\{{\tilde{y}}\}+\mathscr{C}^{-1}\{{\tilde{y}}\mathscr{C}\{{\tilde{y}}\}\},\end{eqnarray}$$
where
$K_{1}$
is an integration constant. By averaging the equation over one wavelength, then using the relations (C 3) and (C 5), one obtains an expression for the constant
$K_{1}$
,
$$\begin{eqnarray}K_{1}={\textstyle \frac{1}{2}}\langle \langle {\tilde{y}}^{2}\rangle \rangle -Bg^{-1}\langle \langle {\tilde{y}}\rangle \rangle \geqslant 0.\end{eqnarray}$$
For numerical resolution, it is convenient to make the change of dependent variable
${\tilde{y}}(\unicode[STIX]{x1D6FC})\stackrel{def}{=}\unicode[STIX]{x1D6F6}(\unicode[STIX]{x1D6FC})+\unicode[STIX]{x1D6FF}$
, where
$\unicode[STIX]{x1D6FF}$
is a constant at our disposal. Thus, (4.4) becomes
$$\begin{eqnarray}(Bg^{-1}-2\unicode[STIX]{x1D6FF})\unicode[STIX]{x1D6F6}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6F6}^{2}+K_{2}=(1+\unicode[STIX]{x1D6FF}d^{-1})c_{R}c_{S}^{-1}\mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\}+\mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\mathscr{C}\{\unicode[STIX]{x1D6F6}\}\},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle K_{2} & = & \displaystyle K_{1}-(d-B/g)\unicode[STIX]{x1D6FF}-{\textstyle \frac{3}{2}}\unicode[STIX]{x1D6FF}^{2}\nonumber\\ \displaystyle & = & \displaystyle (\unicode[STIX]{x1D6FF}-Bg^{-1})\langle \langle \unicode[STIX]{x1D6F6}\rangle \rangle +{\textstyle \frac{1}{2}}\langle \langle \unicode[STIX]{x1D6F6}^{2}\rangle \rangle -(d+\unicode[STIX]{x1D6FF})\unicode[STIX]{x1D6FF}\end{eqnarray}$$
and we have (see appendix C for details)
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \unicode[STIX]{x1D6F6}\rangle \rangle =(c_{S}c_{E}^{-1}-1)d-\unicode[STIX]{x1D6FF}=\langle \langle {\tilde{y}}\rangle \rangle -\unicode[STIX]{x1D6FF}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \unicode[STIX]{x1D6F6}\mathscr{C}\{\unicode[STIX]{x1D6F6}\}\rangle \rangle =-c_{R}c_{S}^{-1}(1+2\unicode[STIX]{x1D6FF}d^{-1})\langle \langle \unicode[STIX]{x1D6F6}\rangle \rangle -c_{R}c_{S}^{-1}(1+\unicode[STIX]{x1D6FF}d^{-1})\unicode[STIX]{x1D6FF}, & \displaystyle\end{eqnarray}$$
giving two equations for
$\unicode[STIX]{x1D6FF}$
and
$c_{S}/c_{E}$
. Alternatively, a related equation is obtained by applying the operator
$\mathscr{C}$
to (4.6),
In practice, the parameter
$\unicode[STIX]{x1D6FF}$
may be conveniently chosen such that:
(i)
$\unicode[STIX]{x1D6FF}=0$
, so we are dealing with the free surface
${\tilde{y}}$
;(ii)
$\unicode[STIX]{x1D6FF}=\langle \langle {\tilde{y}}\rangle \rangle$
, so we are dealing with the zero-mean dependent variable
${\tilde{Y}}$
;(iii)
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
, so
$\unicode[STIX]{x1D6F6}\geqslant 0$
and the wave appears to somehow look like a ‘solitary wave’ on
$[-\unicode[STIX]{x03C0}/k,\unicode[STIX]{x03C0}/k]$
;(iv)
$K_{2}=0$
, so (4.6) is homogeneous.
Properties (iii) and (iv) are both desirable features for the use of Petviashvili iterations, leading to a simple and fast numerical scheme. However, no choice of
$\unicode[STIX]{x1D6FF}$
can enforce (iii) and (iv) simultaneously, so we proceed as follows.
It should be noted first that in deep water
$d\rightarrow \infty$
,
$\mathscr{C}^{-1}\rightarrow \mathscr{C}_{\infty }^{-1}$
is singular and (4.6) is problematic because
$\mathscr{C}_{\infty }^{-1}$
is applied to the quantity
$\unicode[STIX]{x1D6F6}\mathscr{C}_{\infty }\{\unicode[STIX]{x1D6F6}\}$
which has a non-zero (unknown a priori) mean value. This problem is overcome by applying the operator
$\mathscr{C}_{\infty }$
to (4.6) (the constant
$K_{2}$
then vanishes), thus yielding the alternative homogeneous equation
$$\begin{eqnarray}\left(\frac{B}{g}-2\unicode[STIX]{x1D6FF}\right)\mathscr{C}_{\infty }\{\unicode[STIX]{x1D6F6}\}-\mathscr{C}_{\infty }\left\{\frac{\unicode[STIX]{x1D6F6}^{2}}{2}\right\}=\left(1+\frac{\unicode[STIX]{x1D6FF}}{d}\right)\frac{c_{R}}{c_{S}}\mathscr{C}_{\infty }\circ \mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\}+\mathscr{C}_{\infty }\circ \mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\mathscr{C}\{\unicode[STIX]{x1D6F6}\}\},\end{eqnarray}$$
where
The equation (4.11) for deep water with
$\unicode[STIX]{x1D6FF}=0$
was solved by Dyachenko et al. (Reference Dyachenko, Lushnikov and Korotkevich2013) using the generalised Petviashvili method (GPM) because the classical Petviashvili method (CPM) does not converge (with
$\unicode[STIX]{x1D6FF}=0$
). Here, we solve (4.11) with
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
using classical Petviashvili iterations that converge in deep water and in finite depth.
5 Numerical resolution
We detail here how the numerical algorithm is implemented. We mostly focus on periodic waves (
$k>0$
), the case of solitary waves (
$k=0$
) being subsequently briefly described. Details of alternative computations of solitary waves are given in Clamond & Dutykh (Reference Clamond and Dutykh2013) and Dutykh & Clamond (Reference Dutykh and Clamond2014).
5.1 User provided parameters
For practical applications, the user generally wants to define the wave by the parameters
$g$
,
$d$
,
$k$
and
$H$
. Two of these parameters can be freely chosen to compute the solution. In deep water and finite depth (i.e. not shallow), a wave is most often defined only by the two dimensionless parameters
$kd$
and
$\unicode[STIX]{x1D700}=kH/2$
, from which the four physical parameters are defined as follows in the numerical resolution. In shallow water, the dimensionless height
$H/d$
is often used instead of
$\unicode[STIX]{x1D700}$
. It should be noticed that
$\unicode[STIX]{x1D700}$
is not a good parameter to characterise the steepness in shallow water (
$\unicode[STIX]{x1D700}=0$
for all solitary waves), while
$H/d$
is not suitable in deep water (
$H/d=0$
for all waves in infinite depth). Alternative definitions of the wave steepness, valid for all depths and all wavelengths, could be introduced, such as
$kH/\tanh (kd)$
. However, as this alternative steepness parameter is generally not used by practitioners, it is not considered here.
The water depth is considered to be infinite if
$$\begin{eqnarray}1-\tanh (kd)\leqslant \text{tol},\end{eqnarray}$$
where
$\text{tol}$
is the numerical tolerance used for the computations. Working with
$n$
-digit arithmetic, we generally take
$\text{tol}=10^{2-n}$
, yielding
$\text{tol}=10^{-14}$
in double precision,
$\text{tol}=10^{-32}$
in quadruple precision and
$\text{tol}=10^{-69}$
in octuple precision (IEEE Standard 754-2008). If the inequality (5.1) is fulfilled, we consider that
$d=\infty$
and we choose (dimensionless) units such that
$g=k=1$
(this choice is always possible via a suitable scaling). Otherwise, we choose
$g=d=1$
(which is also always possible without loss of generality). Many other scalings could be used, but the ones above are the most common in both deep water and finite depth.
The parameters
$g$
,
$d$
and
$k$
being now defined, and the parameter
$\unicode[STIX]{x1D700}$
being chosen, the total wave height is obviously
$H=2\unicode[STIX]{x1D700}/k$
. It should be noted that the wave can be defined by other pairs of dimensionless parameters than
$\{kd,\unicode[STIX]{x1D700}\}$
or
$\{kd,H/d\}$
. For instance, a parameter like
$|\tilde{u} _{0}/c_{R}|$
(where
$\tilde{u} _{0}$
is the horizontal velocity at the crest) is often used to investigate the (almost) highest waves (Maklakov Reference Maklakov2002; Maklakov & Petrov Reference Maklakov and Petrov2015). However, this parameter involves quantities generally unknown from, e.g., experimental measurements, so it is not convenient for practical use. Thus, the definition a wave by such parameters would require important pre- and postprocessing for practical applications, which is not always trivial.
5.2 Computational parameters
In addition to the physical parameters defined above, it is convenient (but unessential) to take
$c_{R}=c_{E}$
and to introduce the dimensionless parameter
$\unicode[STIX]{x1D70E}\stackrel{def}{=}c_{S}/c_{E}$
(hence, 
, 
). In infinite depth,
$\unicode[STIX]{x1D70E}=1$
, but in finite depth,
$\unicode[STIX]{x1D70E}$
is unknown and therefore must be computed (see below). We also chose
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
because we shall use Petviashvili iterations.
The fundamental period is discretised with
$2N$
equally spaced points 
for
$j=1,2,\ldots ,2N$
. The number of nodes should be large enough to ensure the desired accuracy. For infinitesimal short waves (i.e.
$\unicode[STIX]{x1D700}\ll 1$
and
$kd$
not small),
$N=32$
(say) may be sufficient to ensure machine precision when working in double precision. However, the number of required nodes increases rapidly as the wave height and the wavelength increase. Thus, for steep and cnoidal waves, the number of nodes must be substantially larger. Since we are solving the equation with an algorithm requiring
$O(N\log N)$
operations,
$N$
can be large without significant speed burden. Therefore, one can take
$N=1024$
as default when working in double precision because it is often sufficient. If more nodes are needed, the cost of increasing
$N$
is low because the algorithm is very fast. However,
$N$
cannot be too large in order to avoid large accumulation of round-off errors. (When working in double precision, it is wise to take
$N\leqslant 2^{19}=524\,288$
, say.) If an extremely large number of nodes are required, then quadruple or higher precision is necessary to avoid a significant accumulation of round-off errors.
5.3 Classical Petviashvili method
Petviashvili iterations represent a simple and efficient algorithm for computing solitary wave solutions of the Babenko equation (Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014). This algorithm is very easy to implement and runs rapidly since each iteration requires
$O(N)$
operations. It should be noted that Petviashvili iterations are used here together with fast Fourier transforms (FFTs) that require
$O(N\log N)$
operations, so the overall complexity of our algorithm is
$N\log N$
.
In order to apply Petviashvili’s method, we separate the linear and nonlinear terms and rewrite (4.11) as
$\mathscr{L}\{\unicode[STIX]{x1D6F6}\}=\mathscr{N}\{\unicode[STIX]{x1D6F6}\}$
, with
$$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{L}\{\unicode[STIX]{x1D6F6}\}\stackrel{def}{=}(B/g-2\unicode[STIX]{x1D6FF})\mathscr{C}_{\infty }\{\unicode[STIX]{x1D6F6}\}-(1+\unicode[STIX]{x1D6FF}/d)\unicode[STIX]{x1D70E}^{-1}\mathscr{C}_{\infty }\circ \mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{N}\{\unicode[STIX]{x1D6F6}\}\stackrel{def}{=}{\textstyle \frac{1}{2}}\mathscr{C}_{\infty }\{\unicode[STIX]{x1D6F6}^{2}\}+\mathscr{C}_{\infty }\circ \mathscr{C}^{-1}\{\unicode[STIX]{x1D6F6}\mathscr{C}\{\unicode[STIX]{x1D6F6}\}\}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}=c_{S}/c_{E}$
,
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
and
$\unicode[STIX]{x1D6F6}=y-\unicode[STIX]{x1D6FF}$
. From the approximation
$\unicode[STIX]{x1D6F6}^{(i)}$
obtained at the
$i$
th iteration, the next approximation
$\unicode[STIX]{x1D6F6}^{(i+1)}$
is obtained via the Petviashvili weighted fixed-point iteration
$$\begin{eqnarray}\unicode[STIX]{x1D6F6}^{(i+1)}=S_{i}^{2}\times \mathscr{L}^{-1}\circ \mathscr{N}\{\unicode[STIX]{x1D6F6}^{(i)}\},\quad S_{i}\stackrel{def}{=}\langle \langle \unicode[STIX]{x1D6F6}^{(i)}\mathscr{L}\{\unicode[STIX]{x1D6F6}^{(i)}\}\rangle \rangle /\langle \langle \unicode[STIX]{x1D6F6}^{(i)}\mathscr{N}\{\unicode[STIX]{x1D6F6}^{(i)}\}\rangle \rangle .\end{eqnarray}$$
At each iteration, zero values of
$\unicode[STIX]{x1D6F6}^{(i)}$
at the trough and the wave height are enforced via the renormalisation
This renormalisation improves the convergence, especially for steep waves.
In finite depth, the operator
$\mathscr{L}^{-1}$
is singular but
$\mathscr{L}^{-1}\circ \mathscr{N}$
is regular. This is because both operators involve the factor
$\mathscr{C}_{\infty }$
which is zero for zero frequency (this factor is introduced to kill the constant
$K_{2}$
). In doing this, the Babenko equation has been singularised at zero frequency, but this singularity is only apparent (i.e. it is movable in the sense that
$\mathscr{C}_{\infty }^{-1}\circ \mathscr{C}_{\infty }$
is identity). However, it is not necessary to define
$\mathscr{L}^{-1}\circ \mathscr{N}$
explicitly at zero frequency, as the mean value of
$\unicode[STIX]{x1D6F6}^{(i+1)}$
is enforced by the renormalisation (5.5). Indeed, by setting
$\mathscr{L}^{-1}\{1\}\equiv 0$
arbitrarily,
$\unicode[STIX]{x1D6F6}^{(i+1)}$
computed with (5.4) is obtained modulo an unknown (generally incorrect) mean value, the right value of
$\langle \langle \unicode[STIX]{x1D6F6}\rangle \rangle$
being subsequently enforced via (5.5).
As initial guess
$\unicode[STIX]{x1D6F6}^{(0)}$
, we take the linear approximation
unless a better guess is provided by the user, for instance from another calculation with slightly different parameters
$kd$
and
$\unicode[STIX]{x1D700}$
(useful for analytic continuations). We found that, with the initial guess (5.6), the Petviashvili iterations always converge, even for large waves in shallow water, so we did not try out other guesses. The convergence from the initial guess (5.6) illustrates the robustness of the method.
The Petviashvili iterations (5.4) involve the unknown parameters
$B$
,
$\unicode[STIX]{x1D70E}$
and
$\unicode[STIX]{x1D6FF}$
via the definition of the operators
$\mathscr{C}$
,
$\mathscr{L}$
and
$\mathscr{N}$
. Therefore,
$B$
,
$\unicode[STIX]{x1D70E}$
and
$\unicode[STIX]{x1D6FF}$
must be computed from
$\unicode[STIX]{x1D6F6}^{(i)}$
before (5.4) can be used. (It would be the same with any other iterations, such as Newton or Levenberg–Marquardt methods.) These parameters are obtained as follows.
5.4 Computation of the unknown parameters
We first compute
${\tilde{Y}}=\unicode[STIX]{x1D6F6}^{(i)}-\langle \langle \unicode[STIX]{x1D6F6}^{(i)}\rangle \rangle$
. In deep water,
$\unicode[STIX]{x1D70E}=1$
and, according to (3.30), we have
$$\begin{eqnarray}\langle \langle {\tilde{y}}\rangle \rangle =-\langle \langle {\tilde{Y}}\mathscr{C}_{\infty }\{{\tilde{Y}}\}\rangle \rangle ;\end{eqnarray}$$
thence,
${\tilde{y}}={\tilde{Y}}+\langle \langle {\tilde{y}}\rangle \rangle$
by definition of
${\tilde{Y}}$
.
In finite depth, in general,
$\unicode[STIX]{x1D70E}\neq 1$
is unknown (
$\unicode[STIX]{x1D70E}=1$
only for solitary waves) and must be computed. To do so, the relation (3.30) is rewritten as the equation
$$\begin{eqnarray}E(\unicode[STIX]{x1D70E})\stackrel{def}{=}\langle \langle {\tilde{Y}}\mathscr{C}\{{\tilde{Y}}\}\rangle \rangle +(\unicode[STIX]{x1D70E}-1)d=0.\end{eqnarray}$$
It should be recalled here that
$\mathscr{C}$
depends on
$\unicode[STIX]{x1D70E}$
, so (5.8) is a nonlinear equation for
$\unicode[STIX]{x1D70E}$
. Equation (5.8) is thus solved with Newton iterations,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{j+1}=\unicode[STIX]{x1D70E}_{j}-\frac{E(\unicode[STIX]{x1D70E}_{j})}{E^{\prime }(\unicode[STIX]{x1D70E}_{j})},\quad E^{\prime }(\unicode[STIX]{x1D70E})\stackrel{def}{=}\frac{\text{d}E(\unicode[STIX]{x1D70E})}{\text{d}\unicode[STIX]{x1D70E}}=d-d\langle \langle {\tilde{Y}}\mathscr{S}^{2}\{{\tilde{Y}}\}\rangle \rangle ,\end{eqnarray}$$
with
$\mathscr{S}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}\csc [\unicode[STIX]{x1D70E}d\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FC}}]$
. In practice, one Newton iteration is sufficient because the initial guess
$\unicode[STIX]{x1D70E}_{0}$
is given by the approximation of
$\unicode[STIX]{x1D70E}$
obtained at the previous iteration from
$\unicode[STIX]{x1D6F6}^{(i-1)}$
which, if the CPM converges, is close to the exact solution. Once
$\unicode[STIX]{x1D70E}$
has been obtained, we compute
$\langle \langle {\tilde{y}}\rangle \rangle =(\unicode[STIX]{x1D70E}-1)d$
and
$\unicode[STIX]{x1D6FF}=\langle \langle {\tilde{y}}\rangle \rangle -\langle \langle \unicode[STIX]{x1D6F6}\rangle \rangle$
, so these parameters are now known for the
$i$
th Petviashvili iteration. It should be emphasised that accurate computations of
$\unicode[STIX]{x1D70E}$
and
$\langle \langle {\tilde{y}}\rangle \rangle$
are absolutely crucial to ensure that the mapping
$z\mapsto \unicode[STIX]{x1D701}$
is conformal and that the still water level is where it should be.
Finally, the Bernoulli constant
$B$
is obtained from (4.10) applied at the crest (
$\unicode[STIX]{x1D6FC}=0$
) and at the trough (
), i.e.
thence,
with the notations
$[f]_{a}\stackrel{def}{=}f(a)$
and
$[f]_{a}^{b}\stackrel{def}{=}f(b)-f(a)$
. Relation (5.12) is obtained by subtracting (5.10) and (5.11); thus,
$K_{2}$
vanishes and does not need to be computed.
All of the parameters involved in the Babenko equation are now defined and the Petviashvili iterations (5.4) can be applied until the desired accuracy is reached.
5.5 Postprocessing
After convergence of the Petviashvili iterations, all of the parameters of interest can be computed. The celerity
$c_{E}$
is given by
$$\begin{eqnarray}c_{E}=\sqrt{B}\left\langle \left\langle \frac{1+\mathscr{C}\{{\tilde{Y}}\}}{(1+\mathscr{C}\{{\tilde{Y}}\})^{2}+{\tilde{Y}}_{\unicode[STIX]{x1D6FC}}^{2}}\right\rangle \right\rangle ^{-1/2};\end{eqnarray}$$
thence,
$c_{S}=\unicode[STIX]{x1D70E}c_{E}$
. These parameters being defined, all of the integral quantities in appendix B are easily computed.
Often, users want to know the velocity and other fields inside the bulk at a given location
$z$
. This can be obtained from the integrals provided in appendix D and discretised according to the trapezoidal rule (Dutykh & Clamond Reference Dutykh and Clamond2014). This is very simple to implement and also very accurate, provided that
$z$
is not too close to the free surface. (Typically, the distance between
$z$
and the free surface should be larger than
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}$
.)
5.6 Remark
Since we use a spectral method to solve a nonlinear pseudo-differential equation, some aliasing errors may occur. With a number of Fourier modes large enough to achieve spectral accuracy, the magnitude of all of the neglected high-frequency modes is smaller than the machine precision (approximately
$10^{-16}$
in double precision). As the equation is quadratic in nonlinearity, in theory, a two-thirds rule (Boyd Reference Boyd2000) should be applied in order to avoid aliasing. However, even if the two-thirds rule is not applied, the aliasing error is approximately machine precision, provided that the aliased frequencies (i.e. the upper third of the spectra) are small (say, of magnitudes of at most
$10^{-8}$
in double precision) so that their products are numerically zero. This means that an anti-aliasing filter is not necessary if
$N$
is large enough, as can be seen in the spectra below. Of course, one can easily introduce an anti-aliasing (i.e. low-pass) filter if this turns out to be necessary for some computations.
6 Numerical examples
The algorithm described above has been implemented in Matlab™ and it is freely downloadable (Clamond Reference Clamond2017b ). This program was written with clarity in mind, so it can be easily understood, modified and translated into any programming language. In particular, this program can be easily modified to run in arbitrary precision, provided that this feature is available to the user (we use the Advanpix Multiprecision Toolbox (Advanpix 2017)).
Influence of the number of Fourier modes for
$kd=\infty$
,
$\unicode[STIX]{x1D700}=0.4401$
. (a) Normal solution obtained with
$N=1024$
. (b) Ghost solution obtained with
$N=2048$
.
6.1 Deep water
In deep water (
$d\rightarrow \infty$
), periodic waves with identical crests are obtained with our algorithm provided that
$\unicode[STIX]{x1D700}\lesssim 0.44$
. Thus, the algorithm converges for rather steep waves (up to approximately
$99.3\,\%$
of the highest waves), the maximum steepness being
$\unicode[STIX]{x1D700}\approx 0.443164$
(Maklakov Reference Maklakov2002).
It should be noted that with
$\unicode[STIX]{x1D6FF}=0$
, the CPM diverges for all steepnesses. This shows that the choice
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
improves the convergence of the CPM significantly, rending the use of the GPM unnecessary. It should also be noted that Fenton’s algorithm (Fenton Reference Fenton1988) converges for
$\unicode[STIX]{x1D700}\lesssim 0.36$
and that it is much slower than the method described here. For
$\unicode[STIX]{x1D700}\lesssim 0.36$
, Fenton’s algorithm and our one match up to approximately six digits (for non-infinitesimal waves), which corresponds to the accuracy of Fenton’s algorithm.
Actually, the present algorithm can converge also for
$\unicode[STIX]{x1D700}>0.44$
, but then ‘ghost’ solutions (Domokos & Holmes Reference Domokos and Holmes2003) are obtained (similar to the one in figure 2
b). For
$\unicode[STIX]{x1D700}=0.44$
, variation of
$N$
leads to the same solution in double and quadruple precision. However, for
$\unicode[STIX]{x1D700}=0.4401$
, with
$N=1024$
, a ‘normal’ solution is obtained (figure 2
a), but a ‘ghost’ (spurious) solution is obtained with
$N=2048$
(figure 2
b). The algorithm behaviour for the highest computable waves is discussed in § 7.
Decay of the Fourier coefficients for
$\unicode[STIX]{x1D700}=0.4$
and
$kd=\infty$
: blue, 16 digits; red, 34 digits; orange, 71 digits.
For
$\unicode[STIX]{x1D700}\leqslant 0.44$
, the algorithm converges rapidly to the solution. (On a 2012 MacBook Pro laptop computer, with
$N=2048$
,
$\text{tol}=10^{-15}$
,
$\unicode[STIX]{x1D700}=0.4$
and in double precision, the solution is obtained in approximately half a second.) Actually, any arbitrary accuracy can be achieved provided that
$N$
is large enough (figure 3). For instance, for
$\unicode[STIX]{x1D700}=0.4$
, the solution is obtained to machine double precision with
$N=512$
(figure 3
a). A value of
$N=1024$
is not sufficient to achieve machine quadruple precision (figure 3
b), the latter being obtained for
$N=2048$
(figure 3
c). However,
$N=2048$
is not sufficient to achieve full octuple precision, which can be obtained with larger
$N$
. Similarly, any accuracy can be obtained provided that
$N$
is large enough. It should be noted that figure 3 clearly shows that aliasing errors are not significant although no special anti-aliasing techniques were applied.
This test illustrates the accuracy and the robustness of the algorithm. Indeed, some algorithms diverge when
$N$
is too large although they converge for smaller
$N$
(for the same steepness); this is a problem not faced by the algorithm described here.
Examples of steep waves in finite depth for
$kd=1$
: (a,c)
$\unicode[STIX]{x1D700}=0.3$
; (b,d)
$\unicode[STIX]{x1D700}=0.31$
; (a,b) free surfaces; (c,d) spectra.
6.2 Finite depth
We found that the CPM and GPM are both divergent when applied to the Babenko equation (4.11) with
$\unicode[STIX]{x1D6FF}=0$
. Conversely, the CPM converges well if one takes
$\unicode[STIX]{x1D6FF}=\min ({\tilde{y}})$
, so we did not try the GPM. As for the deep water case, the algorithm converges for all but the highest waves, i.e. the algorithm converges up to the first maximum of
$B(\unicode[STIX]{x1D700})$
. For higher waves, the algorithm diverges or converges to a ghost solution, the maximum steepness computable depending on the depth (see § 7).
For example, we consider the case
$kd=1$
with steepnesses
$\unicode[STIX]{x1D700}=0.3$
and
$\unicode[STIX]{x1D700}=0.31$
. Although these two large steepnesses are close, they correspond to quite different free surfaces (figure 4
a,b), as can be seen in their spectra (figure 4
c,d). With
$N=2048$
, the case
$\unicode[STIX]{x1D700}=0.3$
is resolved to machine double precision (figure 4
c), while the case
$\unicode[STIX]{x1D700}=0.31$
is resolved only to a mild accuracy (figure 4
d). Machine double precision is achieved for
$\unicode[STIX]{x1D700}=0.31$
with
$N=8192$
, however. This shows that as the wave approaches the highest one, the number
$N$
of Fourier modes has to be increased dramatically (exponentially fast) in order to reach a full spectral precision. This is inherent to the solution formulated in the conformal plane and it has nothing to do with the numerical algorithm described in the present paper. Indeed, for the highest wave with an angular crest, the Fourier spectrum decays algebraically while smaller waves have spectra decaying exponentially rapidly at high frequencies.
This example illustrates, like the one in deep water, the need for a large number of Fourier modes in order to achieve the full precision of a given floating point format. This example also illustrates the rapid increases of
$N$
necessary for an accurate resolution as the steepness increases, and, therefore, the need for a fast algorithm. Since the algorithm described here has an overall complexity
$O(N\log N)$
, the necessity of large
$N$
is not problematic and accurate computations are rapidly achieved.
For
$kd=1$
, the highest computable wave has steepness
$\unicode[STIX]{x1D700}\approx 0.3146$
, which is approximately
$99.6\,\%$
of the highest wave. Limiting values for other
$kd$
values are given in table 3, showing that the highest computable wave, for any given depth and wavelength, is approximately
$99\,\%$
of the maximum one.
6.3 Shallow water
When the depth
$d$
over wavelength
$L=2\unicode[STIX]{x03C0}/k$
ratio is very small, i.e.
$kd\ll 1$
, we are dealing with the so-called shallow water situation. It is well known that Stokes’ expansion fails to approximate such solutions of finite amplitude. Some shallow water approximations have then been proposed to approximate periodic waves, the so-called cnoidal waves. Many numerical algorithms devoted to resolution of the full equations also fail in shallow water. For instance, the algorithm of Fenton (Reference Fenton1988) works only for
$L/d\lesssim 30$
. (For
$L/d>30$
, Fenton’s algorithm does not converge or converges to ghost solutions with spurious oscillations; see figure 3-2 in Fenton (Reference Fenton2015).)
The present algorithm works in shallow water without difficulties. For instance, with
$L/d=1000$
and
$H/d=0.4$
(i.e.
$\unicode[STIX]{x1D700}\approx 0.00126$
and
$kd\approx 0.00628$
), the solution is obtained (in approximately
$0.4~\text{s}$
on a MacBook Pro laptop computer from 2012) using
$N=8192$
Fourier modes which are necessary to achieve machine double precision (figure 5). Another more extreme example is the case
$L/d=10\,000$
and
$H/d=0.7$
, which is computed in double precision with
$N=2^{19}$
and
$\text{tol}=10^{-14}$
after
$237$
iterations in approximately
$1~\text{min}$
. It should be noted that these solutions are obtained from the initial guess (5.6), which is not at all a decent approximation of the solution, thus illustrating the robustness of the algorithm.
Cnoidal wave in very shallow water (
$L/d=1000$
,
$H/d=0.4$
).
Steep cnoidal wave (
$L/d=71$
,
$H/d=0.802$
).
As the steepness or the wavelength increases, the number of Fourier modes
$N$
required to achieve spectral accuracy increases rapidly. For example, with
$L/d=71$
and
$H/d=0.802$
(i.e.
$\unicode[STIX]{x1D700}\approx 0.0355$
and
$kd\approx 0.0885$
), a steep cnoidal wave is obtained (figure 6
a), its computation to full spectral accuracy requiring
$N=2^{17}=13\,1072$
Fourier modes (figure 6
b). (This result was obtained in less than three minutes. A rough estimate suggests that the same computation with Newton or Levenberg–Marquardt iterations, instead of Petviashvili ones, would take several days (possibly weeks) on the same computer.)
For steep cnoidal waves in very shallow water, the number of necessary Fourier modes can be prohibitively large. An alternative is to compute a shorter cnoidal wave and to eventually increase the length of the trough. Indeed, as a cnoidal wave surface decays rapidly from the crest, it rapidly reaches its minimum to machine precision (figure 6 c). The example of figure 6 shows that longer cnoidal waves can be obtained to machine precision by increasing the length of the trough up to the desired wavelength, then redefining the mean water level and the mean depth, as well as all of the related parameters (renormalisation). A similar procedure can be used to compute solitary waves, as shown below.
The KdV parameter
$m$
for
$H/d=0.1$
.
The KdV parameter
$m$
for
$H/d=0.5$
.
6.4 Comparison with a KdV cnoidal wave
Korteweg & de Vries (Reference Korteweg and de Vries1895) proposed an analytic approximation for small-amplitude long periodic waves in shallow water. They coined the term ‘cnoidal’ wave because this approximation can be expressed in terms of the Jacobian cn function. In our notation, the KdV analytic solution can be conveniently written as
$$\begin{eqnarray}\unicode[STIX]{x1D702}=a\frac{Kdn^{2}(\unicode[STIX]{x1D705}x\mid m)-E}{K-E},\quad k=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}}{K},\quad H=\frac{mKa}{K-E},\quad (\unicode[STIX]{x1D705}d)^{2}=\frac{3H}{4md},\end{eqnarray}$$
where
$dn$
is the elliptic function of Jacobi of parameter
$m$
(
$0\leqslant m\leqslant 1$
), and
$K\stackrel{def}{=}K(m)$
and
$E\stackrel{def}{=}E(m)$
are the complete elliptic integrals of the first and second kind respectively (Abramowitz & Stegun Reference Abramowitz and Stegun1965).
Although it is an analytic approximation, the KdV cnoidal wave requires significant computations. Indeed, as a wave is generally defined for given height
$H$
and wavenumber
$k$
, the parameter
$m$
must be determined by numerically solving the equations in (6.1) relating the parameters. For very long waves,
$m$
is very close to one, to such an extent that it cannot be practically computed (see tables 1 and 2). For instance, for a very long small-amplitude cnoidal wave with
$L/d=10^{6}$
and
$H/d=10^{-1}$
, we have
$1-m\approx 5.78\times 10^{-118936}$
, a value that cannot be easily computed. This problem becomes more severe as the amplitude increases (table 2). Thus, for very long waves, the KdV analytic cnoidal solution is useless for practical applications, even if only crude approximations are sufficient. It is then more efficient to solve the KdV equation numerically, for instance with FFT and Petviashvili’s method, as illustrated here for the Babenko equation. However, in doing so, solving the KdV equation is not much less demanding than solving the Babenko equation, so the latter should be preferred. For the extreme example
$H/d=10^{-1}$
and
$L/d=10^{6}$
, the Babenko solution is computed to spectral double precision with
$N=2^{21}$
in approximately
$100~\text{s}$
using our algorithm, while the KdV analytic solution cannot be computed in double precision. (This example is given to illustrate the efficiency of the method, not for its practical interest. On the same hardware, for the same wave, rough estimates indicate that algorithms of complexity
$O(N^{2})$
would take months, and ones of complexity
$O(N^{3})$
would take thousands of years. For a direct determination of
$m$
from the relations (6.1), one would have to use something like two hundred thousand digit computation. An alternative is to derive better conditioned relations via some non-trivial mathematical manipulations, thus losing the analytical simplicity of KdV.) Of course, this drawback is not limited to the KdV analytic solution, it is a fortiori present in all cnoidal-like approximations, such as the solutions of the Boussinesq-like equations. These considerations demonstrate that our algorithm for solving the irrotational Euler equations is also a suitable alternative to simple analytic models.
6.5 Solitary waves
As solitary waves decay exponentially rapidly, their surface elevation reaches zero to machine precision close to the crest. Thus, solitary waves can be efficiently computed in a periodic box, provided that the box is long enough so that the periodisation does not affect the solution. This numerical trick is well known and has been used by many authors.
Various parameters for the highest computable waves.
The steep cnoidal example of figure 6 reaches its minimal elevation to machine precision before its trough at
$x=L/2$
(see figure 6
c). Thus, this cnoidal wave can be considered as a solitary wave computed in a periodic box, but with a different still water level. The actual depth for the solitary wave is
$d_{\infty }\stackrel{def}{=}d+\unicode[STIX]{x1D702}(L/2)=d-b$
, the surface elevation is
$\unicode[STIX]{x1D702}_{\infty }=\unicode[STIX]{x1D702}+b$
and the dimensionless amplitude is
$H/d_{\infty }$
. For the example of figure 6, we obtain
$H/d_{\infty }=0.8236847804878956$
. This result is surprising because the direct computation of solitary waves via the CPM converges only for
$H/d\lesssim 0.79$
(Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014). The solitary wave thus obtained has been compared with an approximation obtained with Tanaka’s method using
$1024$
nodes (Tanaka’s algorithm (Tanaka Reference Tanaka1986) is far too slow to use it with
$2^{18}$
nodes). We found that the two solutions match up to approximately five digits, which is consistent with the accuracy of Tanaka’s method (Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014), confirming the solution obtained by the CPM after renormalisation.
7 Remarks on the highest computable waves
In order to characterise the highest computable waves, we consider several dimensionless parameters. For a given
$kd$
, we denote by
$\hat{\unicode[STIX]{x1D700}}$
,
$\unicode[STIX]{x1D700}_{B}$
and
$\unicode[STIX]{x1D700}_{C}$
the steepnesses corresponding to respectively the highest wave (with a
$120^{\circ }$
inner angle at the crest), the first maximum of the Bernoulli constant
$B$
and the highest computable wave. We also consider the maximum inclination of the free surface
$\unicode[STIX]{x1D703}_{max}$
and the parameters
$$\begin{eqnarray}\unicode[STIX]{x1D717}\stackrel{def}{=}-{\textstyle \frac{1}{2}}\log (1-\unicode[STIX]{x1D700}/\hat{\unicode[STIX]{x1D700}}),\quad A\stackrel{def}{=}\log (\tilde{u} _{1}/\tilde{u} _{0}),\end{eqnarray}$$
where
$\tilde{u} _{0}$
and
$\tilde{u} _{1}$
are the horizontal velocities at the crest and the trough respectively. It should be noted that
$\unicode[STIX]{x1D717}\rightarrow \infty$
and
$A\rightarrow \infty$
as
$\unicode[STIX]{x1D700}\rightarrow \hat{\unicode[STIX]{x1D700}}$
; a wave is then considered steep if these parameters are larger than 2 (Maklakov Reference Maklakov2002).
The Bernoulli constant in deep water as a function of the steepness.
The maximum computable steepness
$\unicode[STIX]{x1D700}_{C}$
is determined (by dichotomy) for each
$kd$
only up to the fourth decimal place, which is sufficient here for the discussion. All of the computations of
$\unicode[STIX]{x1D700}_{C}$
in table 3 were performed in double precision, with
$N=13\,1072$
Fourier positive modes and with tolerance
$\text{tol}=10^{-12}$
for the iterations. The corresponding parameters
$A$
,
$\unicode[STIX]{x1D717}$
and
$\unicode[STIX]{x1D703}_{max}$
are given by truncated decimal expansions (i.e. not rounded to the nearest decimal approximation). In table 3, the steepnesses of the highest waves
$\hat{\unicode[STIX]{x1D700}}$
and the steepnesses
$\unicode[STIX]{x1D700}_{B}$
corresponding to the first maximum of the Bernoulli constant were kindly provided by Professor Dmitri Maklakov, who guaranteed that the six decimals are correct. The values of
$\unicode[STIX]{x1D700}_{B}$
were confirmed by our algorithm, thus providing another validation of the method.
As shown in table 3, the highest computable waves are rather steep, with
$\unicode[STIX]{x1D703}_{max}\approx 30^{\circ }$
,
$A$
and
$\unicode[STIX]{x1D717}$
being significantly larger than 2. More interestingly,
$\unicode[STIX]{x1D700}_{C}$
always exceeds
$\unicode[STIX]{x1D700}_{B}$
,
$\unicode[STIX]{x1D700}_{C}$
corresponding at least to
$99\,\%$
of the maximum steepness
$\hat{\unicode[STIX]{x1D700}}$
. This is somewhat surprising because similar algorithms for periodic waves in deep water (Dyachenko et al.
Reference Dyachenko, Lushnikov and Korotkevich2013) and solitary waves (Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014) have maximum computable steepnesses
$\unicode[STIX]{x1D700}_{C}\approx \unicode[STIX]{x1D700}_{B}$
.
Indeed, in their implementation of the GPM for the Babenko equation (4.11) with
$\unicode[STIX]{x1D6FF}=0$
, Dyachenko et al. (Reference Dyachenko, Lushnikov and Korotkevich2013) found convergence for
$\unicode[STIX]{x1D700}\lesssim 0.436$
, and, for solitary waves, the present authors (Clamond & Dutykh Reference Clamond and Dutykh2013; Dutykh & Clamond Reference Dutykh and Clamond2014) noticed that the CPM converges for
$H/d\lesssim 0.79$
. These maximum computable steepnesses correspond to the first maximum of the Bernoulli constant (figure 7, table 3). The limiting value
$\unicode[STIX]{x1D700}_{C}\approx \unicode[STIX]{x1D700}_{B}$
was then explained by the fact that, for steeper waves, the lack of one-to-one correspondence between the parameters prevents the algorithm from ‘deciding’ which solution should be retained.
The situation is a bit different here, where we found that systematically
$\unicode[STIX]{x1D700}_{C}>\unicode[STIX]{x1D700}_{B}$
(see table 3), even for solitary waves (see § 6.5). Clearly, the modifications and rescaling, introduced to make the CPM work for arbitrary wavelengths, have a beneficial effect on the steepest computable waves. A possible explanation is that the rescaling somehow enlarges the region of one-to-one correspondence between the computed parameters. This is a conjecture requiring rigorous mathematical investigations that are far beyond the purpose of the present paper.
8 Discussion
We have described an efficient algorithm for computing steady surface gravity waves for an ideal homogeneous fluid in irrotational motion. After analytic transformations (conformal mapping, rewriting of the conditions at the free surface, change of dependent variables) we ended up with a modified Babenko-like equation that can be solved numerically via the CPM. The algorithm thus obtained is very rapid – with complexity
$O(N\log N)$
, where
$N$
is the number of Fourier modes – and any accuracy can be reached for all depths, provided that
$N$
is large enough and that the steepness is not too close to the limiting one. All waves of practical interest can therefore be computed. To the best of our knowledge, it is the first algorithm that is uniformly valid for all wavelength-over-depth ratios, in practice and not only in principle, and that is moreover accurate in the sense that arbitrary precision numerical solutions can be obtained.
With this algorithm, the computation of steady gravity waves for the irrotational Euler equation is not more demanding than the numerical resolution of simplified models such as KdV. It should be noted that some simplified water wave models, e.g. some variants of the Boussinesq equations for long waves in shallow water, have inhomogeneous nonlinear terms. Therefore, the Petviashvili (or similar) method does not work and Newton or Levenberg–Marquardt methods should be used instead. As the latter has complexity
$O(N^{3})$
, the numerical resolution of these simplified models is more demanding than our algorithm for the irrotational Euler equations.
It is often believed that the CPM works only for localised (solitary) waves. Here, we disproved this belief with strong numerical evidence. Rigorous mathematical results are scarce (Pelinovsky & Stepanyants Reference Pelinovsky and Stepanyants2004; Álvarez & Duràn Reference Álvarez and Duràn2014b ) for the Petviashvili (or similar) method. Deeper mathematical understanding would be beneficial for improving the method, in particular for the computation of almost highest waves. We hope that the numerical evidence presented here will stimulate such investigations.
There exist steady surface gravity waves with different crests and asymmetric waves (Okamoto & Shōji Reference Okamoto and Shōji2001). However, these solutions exist close to the limiting ones where the CPM does not work, at least not as formulated in this paper. Improvement of this fixed-point iteration method may increase the range of computable steepnesses. However, this will not be sufficient for the computation of extreme waves using Fourier decomposition together with the conformal mapping, because a huge number of Fourier modes are needed to achieve high accuracy. This can be understood by considering the limiting wave with a
$120^{\circ }$
inner angular crest. Such solutions have a power
$2/3$
singularity at the crest in the conformal plane (Stokes Reference Stokes1880), and, therefore, the
$n$
th Fourier coefficient decays like
$n^{-5/3}$
as
$n\rightarrow \infty$
. Thus, if the Fourier expansion is truncated after the
$N$
th term, the error decays like
$N^{-2/3}$
as
$N\rightarrow \infty$
. In practice, this means that in order to increase the accuracy by two digits, the number of computed Fourier modes must be multiplied by (roughly) one thousand! Clearly, the ‘brute force’ approach consisting of massively increasing
$N$
is inefficient for extreme waves, even on a powerful computer. To overcome this difficulty, several authors have proposed a change of independent variable such that the corresponding Fourier spectrum decays more rapidly. Another possibility for accurate computations of extreme waves consists of a totally different mathematical reformulation of the problem, so that the numerical resolution is less demanding, as suggested in Clamond (Reference Clamond2018). Needless to say, such approaches are suitable if the number of modes
$N$
required for high accuracy is relatively small and the algorithmic complexity is
$O(N^{2})$
or
$O(N^{3})$
, or if the algorithmic complexity remains
$O(N\log N)$
with
$N$
not too large. To the best of our knowledge, an efficient algorithm for arbitrary precision calculation of extreme waves has yet to be discovered.
Here, we introduced simple tricks in order to successfully apply the CPM for the computation of steady water waves of arbitrary wavelength, in arbitrary depth and to arbitrary precision. These tricks can certainly be used for other equations involved in fluid mechanics.
Acknowledgements
This work was supported by the Spanish Ministerio de Economía y Competitividad under the Research Grant MTM2014-54710-P. The authors would like to thank Professor A. Durán for useful discussions. The authors are also grateful to Professor D. Maklakov for providing some valuable data in table 3 and for comments that helped to improve the paper.
Appendix A. Potential and stream function in a ‘fixed’ frame of reference
The definition (2.10) of
$c_{E}$
implies that the function
$\unicode[STIX]{x1D6F7}\stackrel{def}{=}\unicode[STIX]{x1D719}+c_{E}x$
averages zero along any horizontal line
$y=\text{const.}$
, in particular at the bed
$y=-d$
. This interesting property suggests the introduction of a stream function
$\unicode[STIX]{x1D6F9}$
in this frame of reference where the mean horizontal velocity is zero at the bed and such that
$$\begin{eqnarray}\unicode[STIX]{x1D6F9}\stackrel{def}{=}\unicode[STIX]{x1D713}-\bar{\unicode[STIX]{x1D713}}+c_{E}(y+d),\quad \bar{\unicode[STIX]{x1D6F9}}=0,\quad \tilde{\unicode[STIX]{x1D6F9}}=c_{E}(\unicode[STIX]{x1D702}+d)-c_{S}d;\end{eqnarray}$$
thence,
$$\begin{eqnarray}\langle \unicode[STIX]{x1D6F9}(x,y=\text{const.})\rangle =0,\quad \langle \tilde{\unicode[STIX]{x1D6F9}}(x)\rangle =\langle \unicode[STIX]{x1D6F9}(x,y=\unicode[STIX]{x1D702}(x))\rangle =(c_{E}-c_{S})d.\end{eqnarray}$$
The definitions of
$c_{E}$
and
$c_{S}$
also yield (by integrating by parts and exploiting the irrotationality)
$$\begin{eqnarray}\displaystyle \langle \tilde{\unicode[STIX]{x1D6F9}}\rangle & = & \displaystyle \left\langle \int _{-d}^{\unicode[STIX]{x1D702}}u\,\text{d}y+c_{E}d\right\rangle =\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}u\,\text{d}y-d\bar{u}\right\rangle =\left\langle \unicode[STIX]{x1D702}\tilde{u} -\int _{-d}^{\unicode[STIX]{x1D702}}yu_{y}\,\text{d}y\right\rangle \nonumber\\ \displaystyle & = & \displaystyle \left\langle \unicode[STIX]{x1D702}\tilde{u} -\int _{-d}^{\unicode[STIX]{x1D702}}yv_{x}\,\text{d}y\right\rangle =\left\langle \unicode[STIX]{x1D702}\tilde{u} +\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}_{x}\tilde{v}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\int _{-d}^{\unicode[STIX]{x1D702}}yv\,\text{d}y\right\rangle =\langle \unicode[STIX]{x1D702}\tilde{\unicode[STIX]{x1D719}}_{x}\rangle ;\end{eqnarray}$$
thence,
$$\begin{eqnarray}\langle \unicode[STIX]{x1D702}\tilde{\unicode[STIX]{x1D719}}_{x}\rangle =(c_{E}-c_{S})d,\end{eqnarray}$$
which is exploited in the relation (C 1). It should be noticed that the quantity
$(c_{S}-c_{E})d$
(related to the wave impulse; see appendix B) must remain bounded as
$d\rightarrow \infty$
, implying that
$c_{S}\rightarrow c_{E}$
as
$d\rightarrow \infty$
. It should be noted also that, as the velocity magnitude varies monotonically from the bottom to the surface, we necessarily have
$c_{E}\geqslant c_{S}$
(Constantin Reference Constantin2013) and therefore
$\langle \langle {\tilde{y}}\rangle \rangle \leqslant 0$
(see (C 1)).
Appendix B. Integral quantities
The wave can be characterised by several integral parameters (McCowan Reference McCowan1891; Starr Reference Starr1947; Longuet-Higgins Reference Longuet-Higgins1975, Reference Longuet-Higgins1984). In the frame of reference moving with the wave, there are three physically important constants: the fluid flow
${\mathcal{Q}}_{0}$
, the momentum flux
${\mathcal{S}}_{0}$
and the energy flux
${\mathcal{F}}_{0}$
, defined by
$$\begin{eqnarray}{\mathcal{Q}}_{0}\stackrel{def}{=}\int _{-d}^{\unicode[STIX]{x1D702}}u\,\text{d}y=\tilde{\unicode[STIX]{x1D713}}-\bar{\unicode[STIX]{x1D713}}=-c_{S}d,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{S}}_{0} & \stackrel{def}{=} & \displaystyle \int _{-d}^{\unicode[STIX]{x1D702}}(p+u^{2})\,\text{d}y=[(p+u^{2})(y+d)]_{-d}^{\unicode[STIX]{x1D702}}-\int _{-d}^{\unicode[STIX]{x1D702}}(p_{y}+2uu_{y})(y+d)\,\text{d}y\nonumber\\ \displaystyle & = & \displaystyle \tilde{u} ^{2}(\unicode[STIX]{x1D702}+d)+\int _{-d}^{\unicode[STIX]{x1D702}}(g-uv_{x}-vu_{x})(y+d)\,\text{d}y\nonumber\\ \displaystyle & = & \displaystyle (\tilde{u} ^{2}+\tilde{v}^{2})(\unicode[STIX]{x1D702}+d)+\frac{1}{2}g(\unicode[STIX]{x1D702}+d)^{2}-\unicode[STIX]{x2202}_{x}\int _{-d}^{\unicode[STIX]{x1D702}}uv(y+d)\,\text{d}y\nonumber\\ \displaystyle & = & \displaystyle B(\unicode[STIX]{x1D702}+d)+2gd(\unicode[STIX]{x1D702}+d)-\frac{3}{2}g(\unicode[STIX]{x1D702}+d)^{2}-\unicode[STIX]{x2202}_{x}\int _{-d}^{\unicode[STIX]{x1D702}}uv(y+d)\,\text{d}y,\end{eqnarray}$$
$$\begin{eqnarray}{\mathcal{F}}_{0}\stackrel{def}{=}\int _{-d}^{\unicode[STIX]{x1D702}}\left[p+\frac{1}{2}u^{2}+\frac{1}{2}v^{2}+gy\right]u\,\text{d}y=\int _{-d}^{\unicode[STIX]{x1D702}}\frac{1}{2}Bu\,\text{d}y=\frac{1}{2}B{\mathcal{Q}}_{0}.\end{eqnarray}$$
These quantities are related to other averaged quantities of physical interest (see below). In particular, by averaging
${\mathcal{S}}_{0}$
over one wavelength and exploiting the relation (2.8) and the impermeability of the free surface, one obtains at once
$$\begin{eqnarray}{\mathcal{S}}_{0}=Bd+{\textstyle \frac{1}{2}}gd^{2}-3{\mathcal{V}},\end{eqnarray}$$
where
${\mathcal{V}}\stackrel{def}{=}\langle (g\unicode[STIX]{x1D702}^{2})/2\rangle$
is the potential energy of the gravity force.
Other integral quantities can be defined relative to the uniform flow of speed
$-c_{R}$
, i.e. in the fixed frame of reference where the phase velocity is
$c_{R}$
. The integral quantities of interest here are the following.
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \text{Circulation:}\quad {\mathcal{C}}\stackrel{def}{=}\langle \tilde{u} +c_{R}+\tilde{v}\unicode[STIX]{x1D702}_{x}\rangle =c_{R}-c_{E}.\\ \displaystyle \text{Impulse:}\quad {\mathcal{I}}\stackrel{def}{=}\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}(u+c_{R})\,\text{d}y\right\rangle =(c_{R}-c_{S})d.\\ \displaystyle \text{Kinetic energy:}\quad {\mathcal{K}}\stackrel{def}{=}\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}\frac{1}{2}[(u+c_{R})^{2}+v^{2}]\,\text{d}y\right\rangle =\frac{1}{2}c_{R}{\mathcal{I}}-\frac{1}{2}dc_{S}{\mathcal{C}}.\\ \displaystyle \text{Radiation stress:}\quad {\mathcal{S}}_{xx}\stackrel{def}{=}\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}[p+(u+c_{R})^{2}+gy]\,\text{d}y\right\rangle =2c_{R}{\mathcal{I}}-2{\mathcal{V}}+(B-c_{R}^{2})d.\\ \displaystyle \text{Momentum flux:}\quad {\mathcal{S}}\stackrel{def}{=}\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}[p+(u+c_{R})^{2}]\,\text{d}y\right\rangle ={\mathcal{S}}_{xx}-{\mathcal{V}}+\frac{1}{2}gd^{2}.\\ \displaystyle \text{Energy flux:}\quad {\mathcal{F}}\stackrel{def}{=}\left\langle \int _{-d}^{\unicode[STIX]{x1D702}}[p+\frac{1}{2}(u+c_{R})^{2}+\frac{1}{2}v^{2}+gy](u+c_{R})\,\text{d}y\right\rangle \\ \hspace{62.00006pt}={\displaystyle \frac{1}{2}}(B-c_{R}^{2})c_{R}d+{\displaystyle \frac{1}{2}}(B+c_{R}^{2}){\mathcal{I}}+({\mathcal{K}}-2{\mathcal{V}})c_{R}.\\ \displaystyle \text{Group celerity:}\quad c_{g}\stackrel{def}{=}{\mathcal{F}}/({\mathcal{K}}+{\mathcal{V}}).\end{array}\right\}\end{eqnarray}$$
The equalities in these integral relations are easily obtained via some trivial derivations. It should be noted that the radiation stress defined here differs from the definition of Longuet-Higgins (Reference Longuet-Higgins1975); that is,
${\mathcal{S}}_{xx}^{LH}\stackrel{def}{=}{\mathcal{S}}_{xx}-{\mathcal{V}}$
. It should be noted also that the group celerity defined above is not the linear one, i.e.
$c_{g}\neq \unicode[STIX]{x2202}(kc_{0})/\unicode[STIX]{x2202}k$
if
$H\neq 0$
.
Appendix C. Relations between averaged quantities in physical and conformal planes
Averaging in the physical plane is different from averaging in the conformal plane. For practical applications, many averaged quantities may have to be computed. For an easy reference, we give below various connections between averaged quantities in the physical and conformal planes. In particular, we have the special ‘conformal averaged’ relations at the free surface
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle {\tilde{y}}\rangle \rangle =-c_{E}^{-1}\langle \unicode[STIX]{x1D702}\tilde{\unicode[STIX]{x1D719}}_{x}\rangle =(c_{S}c_{E}^{-1}-1)d, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \mathscr{C}\{{\tilde{y}}\}\rangle \rangle =c_{R}c_{E}^{-1}-c_{R}c_{S}^{-1}=c_{R}c_{S}^{-1}d^{-1}\langle \langle {\tilde{y}}\rangle \rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \mathscr{C}^{-1}\{{\tilde{y}}\}\rangle \rangle =(c_{S}c_{E}^{-1}-1)c_{S}c_{R}^{-1}d^{2}=c_{S}c_{R}^{-1}d\langle \langle {\tilde{y}}\rangle \rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle {\tilde{y}}\mathscr{C}\{{\tilde{y}}\}\rangle \rangle =(c_{E}-c_{S})c_{R}c_{E}^{-1}c_{S}^{-1}d=-c_{R}c_{S}^{-1}\langle \langle {\tilde{y}}\rangle \rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \mathscr{C}^{-1}\{{\tilde{y}}\mathscr{C}\{{\tilde{y}}\}\}\rangle \rangle =(1-c_{S}c_{E}^{-1})d^{2}=-d\langle \langle {\tilde{y}}\rangle \rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \tilde{u} \rangle \rangle =-c_{E}^{-1}\langle \tilde{u} (\tilde{u} +\tilde{v}\unicode[STIX]{x1D702}_{x})\rangle =-,c_{E}^{-1}\langle \tilde{u} ^{2}+\tilde{v}^{2}\rangle -c_{E}^{-1}B, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \tilde{v}\rangle \rangle =-c_{E}^{-1}\langle \tilde{v}(\tilde{u} +\tilde{v}\unicode[STIX]{x1D702}_{x})\rangle =-c_{E}^{-1}\langle \unicode[STIX]{x1D702}_{x}(\tilde{u} ^{2}+\tilde{v}^{2})\rangle =0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \tilde{u} ^{2}+\tilde{v}^{2}\rangle \rangle =B-2g\langle \langle {\tilde{y}}\rangle \rangle =B-2g(c_{S}c_{E}^{-1}-1)d, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \tilde{u} ^{-1}\rangle \rangle =-c_{E}^{-1}\langle \tilde{u} ^{-1}(\tilde{u} +\tilde{v}\unicode[STIX]{x1D702}_{x})\rangle =-c_{E}^{-1}\langle 1+\unicode[STIX]{x1D702}_{x}^{2}\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \tilde{w}^{-1}\rangle \rangle =-\frac{1}{c_{E}}\left\langle \frac{\tilde{u} +\tilde{v}\unicode[STIX]{x1D702}_{x}}{\tilde{u} -\text{i}\tilde{v}}\right\rangle =-\frac{1}{c_{E}}\left\langle \frac{1+\unicode[STIX]{x1D702}_{x}^{2}}{1-\text{i}\unicode[STIX]{x1D702}_{x}}\right\rangle =-\frac{\langle 1+\text{i}\unicode[STIX]{x1D702}_{x}\rangle }{c_{E}}=-\frac{1}{c_{E}}, & \displaystyle\end{eqnarray}$$
with
${\tilde{y}}(\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D702}(\tilde{x}(\unicode[STIX]{x1D6FC}))$
, and at the bottom
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \bar{y}\rangle \rangle =c_{E}^{-1}\langle d\bar{\unicode[STIX]{x1D719}}_{x}\rangle =-d, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \bar{u}\rangle \rangle =-c_{E}^{-1}\langle \bar{u}^{2}\rangle =-c_{E}^{-1}B, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \bar{u}^{2}\rangle \rangle =B-2g\langle \langle \bar{y}\rangle \rangle -2\langle \langle \bar{p}\rangle \rangle =B+2gd+2c_{E}\langle \bar{p}\bar{u}\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \langle \bar{u}^{-1}\rangle \rangle =-c_{E}^{-1}, & \displaystyle\end{eqnarray}$$
with
$\bar{y}=-d$
. We also have the special ‘physical averaged’ relations
$$\begin{eqnarray}\displaystyle & \displaystyle \langle x\rangle =\langle \langle \widetilde{xx_{\unicode[STIX]{x1D6FC}}}\rangle \rangle =\langle \langle \overline{xx_{\unicode[STIX]{x1D6FC}}}\rangle \rangle =\unicode[STIX]{x03C0}/k, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D702}\rangle =\langle \langle \widetilde{yx_{\unicode[STIX]{x1D6FC}}}\rangle \rangle =-c_{E}\langle \langle {\tilde{y}}\tilde{u} /(\tilde{u} ^{2}+\tilde{v}^{2})\rangle \rangle =0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \tilde{u} \rangle =-c_{E}\langle \langle \tilde{u} ^{2}/(\tilde{u} ^{2}+\tilde{v}^{2})\rangle \rangle =\langle \tilde{\unicode[STIX]{x1D719}}_{x}/(1+\unicode[STIX]{x1D702}_{x}^{2})\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \tilde{v}\rangle =-c_{E}\langle \langle \tilde{u} \tilde{v}/(\tilde{u} ^{2}+\tilde{v}^{2})\rangle \rangle =\langle \tilde{\unicode[STIX]{x1D719}}_{x}\unicode[STIX]{x1D702}_{x}/(1+\unicode[STIX]{x1D702}_{x}^{2})\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \tilde{u} ^{2}+\tilde{v}^{2}\rangle =\langle \tilde{\unicode[STIX]{x1D719}}_{x}^{2}/(1+\unicode[STIX]{x1D702}_{x}^{2})\rangle . & \displaystyle\end{eqnarray}$$
All of these relations can be easily obtained from their definition and using the transformations (3.8)–(3.11) between the averaging operators. Other relations can be similarly obtained.
Appendix D. Velocity and pressure fields in the fluid
In the numerical procedure described below, we use conformal mapping and a Fourier pseudo-spectral method to solve the equations. This means that we obtain a discrete approximation equally spaced along each streamline. However, for practical applications, it is often necessary to determine the fields (velocity, pressure, etc.) at various positions that are not necessarily the nodes used for the computation. This information can be obtained as follows.
Let
$W(z)\stackrel{def}{=}c_{E}+w(z)$
be the complex velocity observed in the frame of reference where the fluid velocity averages to zero at the bottom. In this frame of reference, it follows that the complex potential
$F(z)\stackrel{def}{=}c_{E}z+f(z)$
(i.e.
$W=\text{d}F/\text{d}z$
) is a periodic function, bounded in the whole fluid domain (
$F$
is unbounded in any other frame of reference).
As the complex velocity is known at the free surface from our approximation procedure,
$W$
at any complex abscissa
$z$
can be obtained from the Cauchy integral
where
$\unicode[STIX]{x1D703}=2\unicode[STIX]{x03C0}$
if
$z$
is strictly inside the fluid domain (i.e.
$\operatorname{Im}(z)<\unicode[STIX]{x1D702}$
),
$\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}$
if
$z$
is at the free surface (i.e.
$\operatorname{Im}(z)=\unicode[STIX]{x1D702}$
) and
$\unicode[STIX]{x1D703}=0$
if
$z$
is strictly above the free surface (i.e.
$\operatorname{Im}(z)>\unicode[STIX]{x1D702}$
). As the bottom impermeability is taken into account via the method of images, the Cauchy integral (D 1) yields for any
$z$
below the free surface
$$\begin{eqnarray}W(z)=\frac{\text{i}}{2\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }\left[\frac{c_{E}\tilde{z}^{\prime }(\unicode[STIX]{x1D6FC})-c_{R}}{\tilde{z}(\unicode[STIX]{x1D6FC})-z}-\frac{c_{E}\tilde{z}^{\prime \ast }(\unicode[STIX]{x1D6FC})-c_{R}}{\tilde{z}^{\ast }(\unicode[STIX]{x1D6FC})-2\text{i}d-z}\right]\,\text{d}\unicode[STIX]{x1D6FC},\end{eqnarray}$$
where
$\tilde{z}(\unicode[STIX]{x1D6FC})=(c_{R}/c_{E})\unicode[STIX]{x1D6FC}+\tilde{X}(\unicode[STIX]{x1D6FC})+\text{i}\unicode[STIX]{x1D702}(\unicode[STIX]{x1D6FC})$
and
$\tilde{z}^{\prime }(\unicode[STIX]{x1D6FC})=\text{d}\tilde{z}/\text{d}\unicode[STIX]{x1D6FC}=(c_{R}/c_{S})+\mathscr{C}\{\unicode[STIX]{x1D702}\}(\unicode[STIX]{x1D6FC})+\text{i}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D6FC})$
,
$\tilde{X}$
and
$\unicode[STIX]{x1D702}$
being known from the numerical resolution of the Babenko equation.
The integral relation (D 2) is not suitable for periodic domains. For the latter, the infinite integral is replaced by one over one period involving the Hilbert kernel,
From this relation, we obtain the derivative of
$W$
(required to compute the acceleration field),
and the complex potential,
such that
$W=\text{d}F/\text{d}z$
and
$\operatorname{Im}(F)=0$
at the bed.
