New conformal mapping for adaptive resolving of the complex singularities of Stokes wave

TL;DR

This paper introduces a new conformal mapping technique that improves the numerical computation of Stokes waves by effectively managing complex singularities, enabling detailed analysis of waves approaching their limiting form.

Contribution

A novel conformal mapping method that accelerates convergence in computing Stokes waves and allows detailed exploration of their limiting behavior.

Findings

Enhanced numerical convergence for Stokes wave computations.

Successful resolution of waves near the limiting wave of greatest height.

Generalization of the mapping to handle multiple singularities.

Abstract

A new highly efficient method is developed for computation of traveling periodic waves (Stokes waves) on the free surface of deep water. A convergence of numerical approximation is determined by the complex singularites above the free surface for the analytical continuation of the travelling wave into the complex plane. An auxiliary conformal mapping is introduced which moves singularities away from the free surface thus dramatically speeding up numerical convergence by adapting the numerical grid for resolving singularities while being consistent with the fluid dynamics. The efficiency of that conformal mapping is demonstrated for Stokes wave approaching the limiting Stokes wave (the wave of the greatest height) which significantly expands the family of numerically accessible solutions. It allows to provide a detailed study of the oscillatory approach of these solutions to the limiting wave. Generalizations of the conformal mapping to resolve multiple singularities are also introduced.

Equations

$$z(w,t)=x(w,t)+\text{mbox}iy(w,t)$$

$$w:=u+\text{mbox}iv,-\infty <u<\infty ,$$

$${\hat{\tilde{z}}}_k\propto e^{-v_c|k|}\text{for}k\gg 1,$$

$$\tilde{z}(w,t)\equiv z(w,t)-w=\tilde{x}(w,t)+\text{mbox}iy(w,t)\text{with}\tilde{x}(w,t)\equiv x(w,t)-w$$

$$z(w)=x(w)+\text{mbox}iy(w)$$

$$q=2\text{arctan}\left[\frac{1}{L}\tan \frac{w}{2}\right],$$

$${\hat{\tilde{z}}}_k\propto e^{-(v_c/L)|k|}\text{for}k\gg 1$$

$$L_{optimal}\simeq {\left(\frac{v_c}{2}\right)}^{1/2}$$

$${\hat{\tilde{z}}}_k\propto e^{-(2v_c{)}^{1/2}|k|}\text{for}k\gg 1.$$

$$u=2\text{arctan}\left[L\tan \frac{q}{2}\right]$$

$$q=\frac{u}{L}$$

$$q_u=\frac{1}{u_q}=\frac{1}{L{\cos }^2\frac{u}{2}\left(1+\frac{1}{L^2}{\tan }^2\frac{u}{2}\right)}=\frac{{\cos }^2\frac{q}{2}\left(1+L^2{\tan }^2\frac{q}{2}\right)}{L}=\frac{1+L^2+(1-L^2)\cos q}{2L}.$$

$$\Delta q=q_u\Delta u+O(\Delta u^2)=\frac{1+L^2+(1-L^2)\cos q}{2L}\Delta u+O(\Delta u^2).$$

$$q_c=2\text{arctanh}\left[\frac{1}{L}\tanh \frac{v_c}{2}\right]$$

$$q_c=\frac{v_c}{L}+O\left(\frac{v_c^3}{L^3}\right)\simeq \frac{v_c}{L}\text{for}{v}_c\ll L$$

$$q=q_{\pm }:=\pm 2\arctan \frac{\text{mbox}i}{L}+2n\pi ,n=0,\pm 1,\pm 2,\dots .$$

$$q_{\pm }=\pm \pi \pm 2\text{mbox}iL+O(\text{mbox}i{L}^3).$$

$$\begin{array}{rl}& \left(c^2{\hat{k}}_u-1\right)y-\left(\frac{{\hat{k}}_u{y}^2}{2}+y{\hat{k}}_{u}y\right)=0,\end{array}$$

$${\hat{H}}_{u}f(u)=\frac{1}{\pi }\text{p.v.}\underset{-\infty }{\overset{+\infty }{\int }}\frac{f(u^{\prime })}{u^{\prime }-u}\text{mbox}d{u}^{\prime }$$

$$({\hat{H}}_{u}f{)}_k=\text{mbox}i\text{sign}(k)f_k,$$

$$f_k=\frac{1}{2\pi }\underset{-\pi }{\overset{\pi }{\int }}f(u)\exp \left(-\text{mbox}iku\right)\text{mbox}du,$$

$$f(u)=\sum _{k=-\infty }^{\infty }{f}_k\exp \left(\text{mbox}iku\right).$$

$$x=u-{\hat{H}}_{u}y,$$

$$\underset{-\pi }{\overset{\pi }{\int }}\eta (x)\text{mbox}dx=\underset{-\pi }{\overset{\pi }{\int }}y(u)x_u(u)\text{mbox}du=\underset{-\pi }{\overset{\pi }{\int }}y(u)[1+{\tilde{x}}_u(u)]\text{mbox}du=0,$$

$$\left(c^2{q}_u{\hat{k}}_q-1\right)y-\left(q_u\frac{{\hat{k}}_q{y}^2}{2}+q_u\hat{y}{\hat{k}}_{q}y\right)=0,$$

$${\hat{k}}_q:=-\frac{\partial }{\partial q}{\hat{H}}_q=-\frac{\partial }{\partial q}({\hat{H}}_q+const)$$

$${\hat{H}}_{q}f(q)=\frac{1}{\pi }\text{p.v.}\underset{-\infty }{\overset{+\infty }{\int }}\frac{f(q^{\prime })}{q^{\prime }-q}\text{mbox}d{q}^{\prime }$$

$$\left(c^2{\hat{k}}_q-\frac{2L}{1+L^2+(1-L^2)\cos q}\right)y-\left(\frac{{\hat{k}}_q{y}^2}{2}+y{\hat{k}}_{q}y\right)=0,$$

$$\underset{-\pi }{\overset{\pi }{\int }}y(q)[u_q+{\tilde{x}}_q(q)]dq=0$$

$${\hat{P}}_u^{-}=\frac{1}{2}(1+\text{mbox}i{\hat{H}}_u)\text{and}{\hat{P}}_u^{+}=\frac{1}{2}(1-\text{mbox}i{\hat{H}}_u)$$

Full text

hydrodynamics, gravity waves, conformal map

I Introduction

The potential flow of ideal fluid of infinite depth with free surface can be efficiently described through the time-dependent conformal mapping

[TABLE]

of the lower complex half-plane $\mathbb{C}^{-}:=\{w|v\leq 0\}$ of the auxiliary complex variable

[TABLE]

into the area $-\infty<x<\infty,$ $y\leq\eta(x,t)$ occupied by the fluid Ovsyannikov1973 ; MeisonOrzagIzraelyJCompPhys1981 ; TanveerProcRoySoc1991 ; TanveerProcRoySoc1993 ; DKSZ1996 ; ZakharovDyachenkoVasilievEuropJMechB2002 , where $y=\eta(x,t)$ is the coordinate of the free surface, $t$ is the time, $x$ and $y$ are the horizontal and vertical physical coordinates, respectively. Here the real line $v=0$ is mapped into the line $y=\eta(x,t)$ representing the free surface of the fluid (see Fig. 1 of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 for the schematic of the conformal mapping (1)). Refs. DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ; LushnikovStokesParIIJFM2016 used this conformal transformation extensively, both analytically and numerically, to reveal the structures of complex singularities of Stokes wave which is the fully nonlinear periodic gravity wave propagating with the constant velocity $c$  Stokes1847 ; Stokes1880 . Nonlinearity of Stokes wave increases with the increase of $H/\lambda,$ where $H$ is the Stokes wave height which is defined as the vertical distance from the crest to the trough of Stokes wave. We use scaled units at which without the loss of generality the spatial period is $\lambda=2\pi$ and $c=1$ for the linear gravity waves similar to Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 . In a Stokes wave $c>1$ and the limit $H\to 0,\ c\to 1$ corresponds to the linear gravity wave. The Stokes wave of the greatest height $H=H_{max}$ (also known as the limiting Stokes wave) has an angle of $2\pi/3$ radians at the crest, corresponding to a singularity $z\sim w^{2/3}$ Stokes1880sup . The non-limiting Stokes waves describe ocean swell and the slow time evolution of the Stokes wave toward its limiting form is one of the possible routes to wave-breaking and whitecapping in full wave dynamics. Wave-breaking and whitecapping carry away significant part of energy and momenta of gravity waves ZKPR2007 ; ZKP2009 . Here slow approach means the time scale which is much larger than the temporal period of the gravity wave of the same spatial period as for the given Stokes wave. Formation of limiting Stokes wave is also considered to be a probable final stage of evolution of a freak (or rogue) waves in the ocean resulting in formation of approximate limiting Stokes wave for a limited period of time with following wave breaking and disintegration of the wave or whitecapping and attenuation of the freak wave into wave of regular amplitude ZakharovDyachenkoProkofievEuropJMechB2006 ; RaineyLonguet-HigginsOceanEng2006 ; DyachenkoNewell2016 .

Thus the approach of non-limiting Stokes wave to the limiting Stokes wave has both significant theoretical and practical interests. It was studied in details in Refs. DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ; LushnikovStokesParIIJFM2016 how the complex singularity in $w$ plane approaches the real line (corresponds to the fluid’s free surface) from above during the transition from non-limiting Stokes wave to the limiting Stokes. Describing such transition is a numerically challenging task because in a Stokes wave the distance $v_{c}$ between the lowest branch points to the real line approaches zero, which implies slow decay of the Fourier coefficients:

[TABLE]

where $k$ is the Fourier wavenumber. Here, similar to Refs. DKSZ1996 ; DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 we separated $z(w,t)$ into a $2\pi$-periodic part $\tilde{z}(2\pi+w,t)=\tilde{z}(w,t)$ and a non-periodic part $w$ by introducing

[TABLE]

such that $\tilde{x}(-\pi,t)=\tilde{x}(\pi,t)=0$ while $x(-\pi,t)=-\pi$ and $x(\pi,t)=\pi$. Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 used up to $M=2^{27}\approx 134\times 10^{6}$ Fourier modes for $\hat{\tilde{z}}_{k}$ on the uniform grid which allowed to obtain the Stokes wave with $v_{c}=5.93824419892803271779\ldots\times 10^{-7}$ (the maximal Fourier mode $k_{max}$ resolved in these simulations corresponds to $k_{max}=M$).

Conformal mappings can be used for improving efficiency of simulations for the general periodic 1D system defined on the real line if such system allows analytic continuation to the strip containing the real axis (see e.g. BoydChebyshevFourierSpectralMethodsBook2001 for review). Assume that $v_{c}$ is the vertical distance from the real line to the complex singularity closest to the real line. Thus $v_{c}$ defines the thickness of the strip of analyticity in the direction where the singularity is nearest to the real line (Stokes wave is a special case because the thickness of strip is determined by the distance $v_{c}$ in the upper complex half-plane $w\in\mathbb{C}^{+}$ while the thickness is infinite below the real line). Then the FT for the system scales as in Eq. (3). The idea is to find a conformal transformation from $w$ to the new complex variable which makes the strip of analyticity thicker, i.e. to push all complex singularities of the system to the distance $\tilde{v}_{c}>v_{c}$ from the real line. Then FT in the new conformal variable scales as $\propto e^{-\tilde{v}_{c}|k|}$ for $k\gg 1,$ i.e. decays faster than in Eq. (3) speeding up numerical convergence. A similar idea can be applied to the nonperiodic systems holomorphic in a closed ellipse around the segment of the real line (with foci corresponding to the two ends of that segment) with e.g. rational spectral interpolants used instead of FT TeeTrefethenSIAMJSCICOMPUT2006 . In all such cases the spectral numerical methods including FT methods are highly efficient and typically having exponential convergence with the number of grid points $M$ used for the spectral collocation as exemplified by Eq. (3) if we use $k=k_{max}=M$ for the estimate of the numerical error. However exploiting such idea for the dynamics of the ideal fluid with free surface on infinite depth has previously met with obstacles because the water waves dynamics require to work with function holomorphic in the entire lower complex half plane instead of the strip. In other words, only singularities in the upper complex half-plane $w\in\mathbb{C}^{+}$ are allowed for the dynamics of the ideal fluid. In this paper we overcome that obstacle by the proper choice of the conformal map.

We focus on problems periodic in $x$ variable (with period $2\pi$) in the reference frame moving with the constant velocity $c$. The transformation (1) becomes independent of time:

[TABLE]

which parametrically defines the Stokes wave as $y(x)$ in physical coordinates with $x(u)$ and $y(u)$ being horizontal and vertical physical coordinates, respectively.

The main results of this paper is that we found a new conformal map

[TABLE]

which we demonstrate to be consistent with the fluid dynamics. We show below the extreme efficiency of equation (6) for simulations of Stokes wave. Here $q$ is the new complex coordinate and $L$ is the arbitrary positive constant which we adjust to optimize performance of simulations. For the general time-dependent problems $L$ would become time-dependent to account for the vertical motion of nearest singularities as well as both $w$ and $q$ can be translated in the horizontal direction with time to account for the horizontal motion of singularities. However such generalization to time-dependent problems is beyond the scope of this paper.

Figs. 1 and 2 show schematically that the new conformal map (6) zooms at the real line into the neighborhood of $w=0.$ Among all points on the real line $w=Re(w),$ that point $w=0$ is the closest to the lowest singularity (branch cut) of the strongly nonlinear Stokes wave which is located at $w=\mbox{i}v_{c},\ v_{c}\ll 1$. Then the uniform grid in the new variable $q$ corresponds to the highly nonuniform grid in the physical coordinates with the grid points concentrating at the neighborhood of the singularity as seen in Fig. 2. The substitution of $w=\mbox{i}v_{c}$   into equation (6) immediately reveals that the lowest singularity of Stokes wave is located at $q=\mbox{i}v_{c}/L+O(\mbox{i}v_{c}^{3}/L^{3})\simeq\mbox{i}v_{c}/L$ in $q$ plane. It means that the free parameter $L$ of the transformation (6) allows to change the position of the singularity in the complex $q$ plane. Here it is assumed that $v_{c}/L\ll 1$. Then FT in $q$ variable decay as

[TABLE]

which is much faster than (3) for $L\ \ll 1.$ It makes the new conformal map (6) highly efficient. Equation (6) has its own singularities at $q=\pm 2\arctan(\mbox{i}/L)+2n\pi=\pm 2\mbox{i}L+(2n+1)\pi+O(\mbox{i}L^{3}),\ n=0,\pm 1.,\pm 2,\ldots$ which approach the real line with the decrease of $L$ as schematically shown in Fig. 1. Balancing the contribution of singularities of Stokes wave and $q(w)$ (i.e. setting them to have the same distance to the real axis in $q$ plane) one obtains the optimal value

[TABLE]

which ensures the fastest possible convergence of Fourier modes as

[TABLE]

E.g., the simulation of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 with $M=2^{27}\simeq 1.3\cdot 10^{8}$ and $v_{c}=5.93824419892803271779\ldots\times 10^{-7}$ required running 64 cores computer cluster for $\sim 3$ months. In contrast, the simulations described in Section VI (they use the new conformal map (6)) allowed to achieve the same precision for the numerical grid with $M_{q}\simeq 4.2\cdot 10^{4}$ Fourier modes which takes a few minutes on the desktop computer. Respectively, by increasing $M_{q}$ (according to equation (9), one has to choose $M_{q}\sim M^{1/2}$ to reach the same precision as on the uniform grid) we were able to study Stokes waves with significantly smaller values of $v_{c}$ (down to $\sim 10^{-11}$) than in Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 .

The new conformal map (6) and its inverse provide the mapping between half-strips in $w$ and $q$ lower complex planes as shown in Fig. 3 by shaded areas. These areas extend all way down in the complex planes and correspond to the area occupied by fluid with the exception of the singularity points of the conformal map as detailed in Section II. These exceptional points result in the extra constant terms found in Section IV to ensure the exact solution of Euler equation through the conformal map.

The paper is organized as follows. In Section II we discuss the properties of the conformal map (6) in complex plane and corresponding discretization. In Secton III we introduce the closed equation for Stokes wave in new $q$ variable. In Section IV we analyze how to work with the projectors and Hilbert transformation in the new variable $q$. In Section V we describe the numerical algorithm used for obtaining the Stokes waves in the limit of the complex singularity approaching the real line. Section VI demonstrates the high efficiency of the new conformal map (6) and analyzes the results on the computed Stokes waves. Section VII provides a generalization of the conformal map (6) to adaptively resolve multiple singularities. In Section VIII the main results of the paper are discussed.

II New spatial coordinate for non-uniform grid

The conformal transformation (6) is $2\pi$-periodic in $w$. Inverting equation (6) at the real line $u=w,$ we obtain that

[TABLE]

which is $2\pi$-periodic in $q$. Also the real line $w=u$ maps into the real line $q=Re(q)$. Recalling that we assume $2\pi$ periodicity of Eq. (4) for Stokes wave, we conclude that it is sufficient to consider the conformal transformations (6) and (10) between half-strip $-\pi\leq u\leq\pi$, $-\infty<v\leq 0$ and $-\pi\leq Re(q)\leq\pi$, $-\infty<Im(q)\leq 0$ in $w$ and $q$, respectively. Here $2\pi$-periodicity is ensured by the limits $q\to\pm\pi$ for $u\to\pm\pi.$

If we assume that $|u|\ll L$ then equation (6) is reduced to

[TABLE]

which implies that taking numerical step $\Delta q\sim 1$ in $q$ space for $q$ near 0 is equivalent to taking the numerical step $\Delta u\sim L\ll 1$ in $u$ space. It ensures that the uniform grid in $q$ space is highly concentrated near $u=0$ in $u$ space, with a ”density” of grid points $q_{u}\sim 1/L$ near $u=0$. It allows to use much less grid points on the uniform grid in $q$ space in comparison with the uniform grid in $u$ space (the uniform grid was used previously in many simulations, see e.g. Refs. ZakharovDyachenkoVasilievEuropJMechB2002 ; ZakharovDyachenkoProkofievEuropJMechB2006 ; DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ) to archive the same precision of a numerical solution. To make these arguments precise we use the Jacobian $q_{u}$ of the transformation (6) is given by

[TABLE]

Then for a general value of $q$, the steps $\Delta q$ and $\Delta u$ are related as

[TABLE]

Fig. 2 shows $q(u)$ and $q_{u}$ for $L=0.0186\ldots$ with dots representing $M_{q}=1024$ points of discrete grids both in $u$ and $q$ spaces separated by $\Delta u$ and $\Delta q$, respectively.

The branch point singularity of Stokes wave, located at $w=\mbox{i}v_{c}$ in $w$ plane, corresponds to $q=\mbox{i}q_{c}$ in $q$ plane in accordance with the Eq. (6) where

[TABLE]

which implies that

[TABLE]

as schematically shown in Fig. 3. Thus $q_{c}$ is located significantly higher in $q$ plane compared to $v_{c}$ in $w$ plane providing the quantitative explantation of much quicker decay of Fourier spectrum (7) in $q$ variable compared to Eq. (3). However, the asymptotic (7) is valid provided $q=\mbox{i}q_{c}$ is closer to the real axis than the other parts of the mapping of the Stokes wave branch cut $w\in[\mbox{i}v_{c},\mbox{i}\infty)$ into $q$ plane. In particular, a one part $w\in[\mbox{i}v_{c},2\mbox{i}L+O(\mbox{i}L^{3}))$ of Stokes wave branch cut is mapped into $q\in[\mbox{i}q_{c},\mbox{i}\infty)$ and another part $w\in(2\mbox{i}L+O(\mbox{i}L^{3}),\mbox{i}\infty)$ is mapped into two branch cuts $q\in(\pm\pi+\mbox{i}\infty,\pm\pi+2\mbox{i}L+O(\mbox{i}L^{3}))$ by $2\pi$ periodicity in $q$ as sketched by vertical lines in Fig. 3. Here the branch points $\pm\pi+2\mbox{i}L+O(\mbox{i}L^{3})$ correspond to the singularities of the conformal map (10) in $q$ space.

These singularities are obtained from the Jacobian (12) which is nonsingular for any value $q\in\mathbb{C}$ but reaches zero (i.e. the singularity of $u_{q}$) at

[TABLE]

For $L\ll 1$ Eq. (16) reduces in the strip $-\pi\leq Re(q)\leq\pi$ to

[TABLE]

The locations $q_{+}$ and $\mbox{i}q_{c}$ of Stokes wave singularities in $q$ space are shown schematically in Fig. 3.

We note that the singularities $q_{-}=\pm\pi-2\mbox{i}L+O(\mbox{i}L^{3})$ are located in $q\in\mathbb{C}^{-}$ but they are invisible for any function $f(q)$ which is $2\pi$ periodic because the jumps at the corresponding branch cut $q\in(\pm\pi-2\mbox{i}L+O(\mbox{i}L^{3}),\pm\pi-\mbox{i}\infty)$ are identically zero, see also Refs. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ; LushnikovStokesParIIJFM2016 for somewhat similar discussion. Because of that we show points $q_{-}$ by filled circles for the transformation in right panel of Fig. 1 but do not show $q_{-}$ in right panel of Fig. 3. If instead of a Stokes wave one would consider a function $\tilde{z}(w)$ with the branch cut of finite extent then the end point of the mapping of that branch cut into $q$ plane would be not $q_{+}$ but other point located higher above the real axis. However we do not consider such functions in this paper.

The singularity $q_{+}$ dominates the asymptotic of Stokes wave FT in $q$ variable provided $|Im(q_{+})|<q_{c}$. Thus the best convergence of FT (faster decays of Fourier harmonics at large $k)$ occurs for $|Im(q_{+})|=q_{c}$ which together with Eqs. (15) and (17) give the optimal choice of the parameter $L=L_{optimal}$ given by Eq. (8) and valid for $v_{c}\ll 1.$ For $L=L_{optimal}$ both singularities of Stokes wave in $q$ space are located at a distance $\approx(2v_{c})^{1/2}$ from the real axis ensuring the FT asymptotic (9).

III Equation of Stokes wave

The closed equation for Stokes wave has the following form

[TABLE]

which is defined on the real line $w=u$ for the function $y(u)$ which is the imaginary part of Eq. (4). Eq. (18) was derived in Ref. BabenkoSovietMathDoklady1987 and later was independently obtained from results of Ref. DKSZ1996 in Ref. DyachenkoLushnikovKorotkevichJETPLett2014 from the exact Euler equations of free surface hydrodynamics. See also Ref. ZakharovDyachenkoPhysD1996 for somewhat similar equation. Here $\hat{k}$ is the positive-definite linear operator defined by $\hat{k}_{u}:=-\frac{\partial}{\partial u}\hat{H}=\sqrt{-\frac{\partial^{2}}{\partial u^{2}}}$ and $\hat{H}_{u}$ is the Hilbert transform,

[TABLE]

with p.v. designating a Cauchy principal value of integral and subscript in $u$ means that both the Hilbert transform and $\hat{k}_{u}$ are defined for the variable $u$. The Hilbert operator $\hat{H}_{u}$ is a multiplication operator on the Fourier coefficients:

[TABLE]

where $f_{k}$ are the Fourier coefficients (harmonics):

[TABLE]

of the periodic function $f(u)=f(u+2\pi)$ represented through the Fourier series

[TABLE]

Here $\text{sign}(k)=-1,0,1$ for $k<0,\ k=0$ and $k>0$, respectively.

After solving Eq. (18) numerically as described in Refs. DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 , we recover the real part $x(u)$ of Eq. (5) from $y(u)$ as

[TABLE]

which follows from the analyticity of $\tilde{z}(w)$ (4) in $\mathbb{C}^{-}$ (see a derivation of Eq. (23) e.g. in Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ). Then Stokes wave solution is represented in the parametric form $(x(u),y(u)).$

Eq. (18) was derived in Ref. DyachenkoLushnikovKorotkevichJETPLett2014 under the assumption that

[TABLE]

meaning that the mean elevation of the free surface is set to zero. Equation (24) reflects a conservation of the total mass of fluid.

In this paper instead of solving Eq. (18) in $u$-variable, we transform it into $q$-variable using (10). Then we solve the resulting equation numerically in a more efficient way using the procedure described in Section V. We express $u$ through $q$ as given by Eq. (10), and obtain from Eqs. (18) and (12) that

[TABLE]

where the operators

[TABLE]

and

[TABLE]

now act in $q$ space with $q_{u}$ given by Eq. (12). Here and below we abuse notation and use the same symbol $y$ for both functions of $u$ and $q$ (in other words, we assume that $\tilde{y}(q)=y(u)$ and remove $\tilde{~}{}$ sign). The comparison of Eqs. (18) and (25) together with Eqs. (19) and (27) reveals that we simply replaced $\hat{H}_{u}$ by $\hat{H}_{q}+const$, where the explicit expression for a constant is not important for solving Eq. (25) because it includes derivatives over $q$ thus removing this constant. The justification of the validity of this nontrivial replacement is provided in Section IV. We also note by comparison of the definitions of $\hat{k}_{u}$ and $\hat{k}_{q}$ above in this Section that FT of $\hat{k}_{q}$ has the same meaning of the multiplication on $|k|$ but this time in Fourier space of $q$.

Because the Jacobian $q_{u}$ is nonzero for any real values $q\in[-\pi,\pi],$ one immediately obtains from equation (25) a more compact expression

[TABLE]

which we use for simulations. The mean level zero condition (24) is transformed to

[TABLE]

in the $q$ variable.

IV Projectors and Hilbert transformation in $q$ variable

In this Section we justify the use of the operator $\hat{k}_{q}$ in Eqs. (25) and (28). It is convenient to introduce the operators

[TABLE]

which are the projector operators of a general periodic function into a functions analytic in $w\in\mathbb{C}^{-}$ and $w\in\mathbb{C}^{+}$ correspondingly. To understand the action of these projector operators, we introduce the splitting of a general $2\pi$ periodic function $f(u)$ with the Fourier series (22) as

[TABLE]

where

[TABLE]

is the analytical (holomorphic) function in $\mathbb{C}^{+}$ and

[TABLE]

is the analytical function in $\mathbb{C}^{-}$ as well as $f_{0,u}$ is the zero Fourier harmonic defined through Eq. (21) as $f_{0,u}=f_{k}|_{k=0}=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}f(u)du$. Together with the property

[TABLE]

which follows from Eq. (20) we obtain that

[TABLE]

i.e. the functions which are holomorphic in $\mathbb{C}^{-}$ and $\mathbb{C}^{+}$, respectively.

Here we use the notation $\hat{P}_{u}^{-}$, $\hat{P}_{u}^{+}$ and $\hat{H}_{u}$ for the projectors and Hilbert transform in $u$ space. Similarly, we introduce projector operators $\hat{P}_{q}^{-}$, $\hat{P}_{q}^{+}$ and Hilbert transform $\hat{H}_{q}$ in the $q$-variable. We can use two approaches to determine the form of projectors $\hat{P}_{u}^{-}$, $\hat{P}_{u}^{+}$ in the $q$ space. The first approach is to analyze how Fourier series transforms as we make a change of variables from $u$ to $q$. The second approach is to use definition of these operators through complex contour integrals and see how these integrals transform as we make a change of variables from $u$ to $q$. In this paper we focus on the second approach as well as we provide the expressions only for $\hat{P}_{q}^{-}$ and $\hat{H}_{q}$. The expression for $\hat{P}_{q}^{+}$ can be derived in a similar way but it is not need for the computation of Stokes wave.

Using the Sokhotskii-Plemelj theorem (see e.g. Gakhov1966 ; PolyaninManzhirov2008 )

[TABLE]

where $\mbox{i}0$ means $\mbox{i}\epsilon,\ \epsilon\to 0^{+}$, we rewrite Eq. (30) as follows

[TABLE]

We now use $2\pi$ periodicity of $f(u)$ to reduce Eq. (37) into the integral over one period

[TABLE]

Using equations (10) and (12) we transform Eq. (38) to $q$-variable as follows

[TABLE]

Similar to Eq. (31) we write $f(q)$ in $q$ space as follows

[TABLE]

where

[TABLE]

is analytic function in $\mathbb{C}^{+}$ and

[TABLE]

is the analytical function in $\mathbb{C}^{-}$ and $f_{0,q}$ is the zero Fourier harmonic $f_{0,q}=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}f(q)dq.$

We evaluate integrals in Eq. (IV) using (40) by closing complex contours in $q\in\mathbb{C}^{+}$ for $f^{+}(q)$ and $q\in\mathbb{C}^{-}$ for $f^{-}(q),$ respectively as shown in Fig. 4. For $f_{0,q}$ it can be done in both ways giving the same result. The zeros of the denominator are located at $q^{\prime}=q-\mbox{i}0$ and $q^{\prime}=\pm 2\arctan{\frac{\mbox{i}}{L}}$. We calculate the residues to obtain:

[TABLE]

Using equations (40) and (43) yields:

[TABLE]

Defining the projector $\hat{P}^{-}_{q}$ in $q$ space similar to Eq. (35) as

[TABLE]

we obtain from Eq. (44) that

[TABLE]

where

[TABLE]

is the constant. We define $\hat{H}_{q}$ through a relation similar to Eq. (30) in $q$ space

[TABLE]

and we obtain from Eq. (46) and (48) that

[TABLE]

Thus the operators $\hat{P}^{-}$ and $\hat{H}$ in $u$ and $q$ spaces are the same except for the shift by a constant $c_{shift}$ and $-2\mbox{i}c_{shift}$ respectively. These constants result from the singularities (16) of the conformal map (10). The explicit expression for $c_{shift}$ is calculated from the values of $f^{-,q}(q)$ as follows. We notice that for Stokes wave $y(u)$ is an even real function, $y(u)=y(-u)\in\mathbb{R},$ which implies that $y(q)=y(-q)\in\mathbb{R}$ in $q$ variable. By taking $f(q)=y(q)$ we obtain that $f^{+,q}(q)=f^{-,q}(-q)$. The analytical continuation of $f^{-,q}(q)$ from the real line $q=Re(q)$ into the complex value $q=-2\arctan{\frac{\mbox{i}}{L}}$ is trivially done by plugging the complex value of $q$ into the series (42) reducing Eq. (47) to

[TABLE]

for the even real function $y(q).$ For more general non-even solution of Eq. (28) (corresponds to higher order progressive waves, which have more than one different peaks per $2\pi$ spatial period ChenSaffmanStudApplMath1980 ) we generally obtain nonzero value of $c_{shift}$ by a similar procedure as follows. We recover $f^{+,q}(q)$ from $f^{-,q}(q)$ using the identity

[TABLE]

which follows from the condition that $f(u)$ is the real-valued function. Here $\bar{f}(q)$ means the complex conjugation of the function $f(q)$ for real values of $q$, i.e. $\bar{f}(q)\equiv\overline{f(\bar{q})}$ for complex values of $q.$ For Eq. (42) it implies that $\bar{f}^{-,q}(q)=\sum\limits_{k=-\infty}^{-1}\bar{f}_{k}\exp\left(-\mbox{i}kq\right).$ Then $c_{shift}$ results from the analytical continuation of $f^{-,q}(q)$ and $\bar{f}^{-,q}(-q)$ from the real line into $q=-2\arctan{\frac{\mbox{i}}{L}}$ together with Eqs. (47) and (51).

Note that we do not need the explicit value for $c_{shift}$ to solve Eqs. (25) and (28) because they both include derivatives over $q$ which removes $c_{shift}$. However, to obtain $x(u)$ one generally needs the value of $c_{shift}$ (which produces only a trivial shift in the horizontal direction). Using Eqs. (10), (49), (47), we transform Eq. (23) into the variable $q$ as follows

[TABLE]

We conclude that this section has justified the derivation of Eqs. (25) and (28) from Eq. (18).

V Numerical algorithm for computing Stokes wave

We solve Eq. (28) numerically using the generalized Petviashvili method (GPM) LY2007 ; PelinovskyStepanyantsSIAMNumerAnal2004 and the Newton Conjugate Gradient method proposed in Refs. JiankeYang2009 ; YangBook2010 . These numerical methods are similar to the numerical solution of Eq. (18) in Refs. DyachenkoLushnikovKorotkevichJETPLett2014 ; DyachenkoLushnikovKorotkevichPartIStudApplMath2016 . For both methods $y(q)$ is expanded in cosine Fourier series and the operator $\hat{k}_{q}$ (26) is evaluated numerically using Fast Fourier Transform (FFT) on the uniform grid with $M_{q}$ points discretization of the interval $-\pi\leq q<\pi$.

As alternative to solving Eq. (28), we also numerically solved the equivalent equation

[TABLE]

which is the analog of equation

[TABLE]

derived in Ref. LushnikovStokesParIIJFM2016 . Eq. (54) is equivalent to Eq. (18) and is obtained by applying the projector operator $\hat{P}_{u}^{-}$ (30) to equation (18) together with the condition (24). In a similar way. Eq. (53) is obtained by applying the projector operator $\hat{P}_{q}^{-}$ (48) to equation (28) together with the condition (29).

Solving Eq. (53) numerically instead of Eq. (28) typically provides 1-2 extra digits of accuracy in Stokes wave height $H$ as well as in the accuracy of the solution spectrum and $v_{c}$. The extra cost is however that we have to solve Eq. (53) for the complex-valued function $\tilde{z}(q)$ instead of the real valued function $y(q)$ in Eq. (28) which doubles memory requirements and the number of numerical operations.

After we obtain a numerical solution for $z(q),$ we use it to determine the value of $v_{c}$ via one of three numerical methods:

(i) The first method uses a least squares fit of Fourier spectrum of a solution $\tilde{z}$ to the asymptotic series described in Eq. (41) of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 . Working in $u$ variable this method allows one to obtain $v_{c}$ with the absolute accuracy about $10^{-10}-10^{-11}$ in double precision (DP) using 7-12 terms in the series of Eq. (41) of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 . While working in $q$ variable, the second singularity (16) located at $q=q_{+}=\pm\pi+2\mbox{i}L+O(\mbox{i}L^{3})$ introduces a contribution to the Fourier spectrum of the same order as the main singularity (14) (located at $q=\mbox{i}q_{c}\simeq\mbox{i}v_{c}/L$) if the parameter $L$ is chosen close to $L_{optimal}\simeq(v_{c}/2)^{1/2}$ (8), so we typically can get only 1-2 digits of precision in $v_{c}$. In order to obtain $v_{c}$ with higher accuracy one needs to remap the solution via Fourier interpolation to a uniform grid for the new variable $\tilde{q}$ with the larger value of the parameter $\tilde{L}$ (we found that a factor 8 or 16 is typically enough to obtain $v_{c}$ with maximum possible accuracy in DP). This pushes the second singularity much further away from the real line compared to the first one so that the main contribution to the tail of Fourier spectrum in $\tilde{q}$ space comes from the first singularity (at a distance $\simeq v_{c}/\tilde{L}$ from the real line) which will allow us to find $v_{c}/\tilde{L}$ (and consequently $v_{c}$) using the same fitting procedure as in $u$ space. Using this approach for solutions in $q$ space we were able to recover $v_{c}$ with absolute accuracy about $10^{-9}-10^{-10}$ in DP. We typically used it for solution with $M_{q}<10^{5}$ Fourier harmonics since Fourier interpolation procedure uses $O(M_{q}^{2})$ operations and becomes slow for larger $M_{q}$.

(ii) The second method is described in Section 6.1 a of Ref. LushnikovStokesParIIJFM2016 and based on the compatibility of the series expansions at points $\zeta=\pm\mbox{i}\chi_{c}$ in the axillary space $\zeta=\tan[\frac{w}{2}]$ (which implies that $\chi_{c}=\tanh[\frac{v_{c}}{2}])$ with the equation (18) of Stokes wave. Current realization of this algorithm also requires Padé approximation of the solution in the axillary space $\zeta$ (described in Section 4 of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ) for calculation of coefficients of series expansion at the point $\zeta=-\mbox{i}\chi_{c}=-\tan[\frac{\mbox{i}v_{c}}{2}]$. This method is so far the most accurate but requires $O(M_{q}\times N_{d})$ operations ($N_{d}$ is the number of poles) for finding Padé approximation of a solution thus slow for large $M_{q}$. We typically used that method for $M_{q}>10^{5},$ where the the small value of $v_{c}$ required us to use quadruple (quad) precision with 32 digits accuracy both to obtain $\tilde{z}$ and recover $v_{c}.$ The absolute accuracy for $v_{c}$ in this method was $\sim 10^{-26}.$

(iii) The third method is described in Section 4.3 of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 and uses nonlinear fit of the crest of a solution to a series (4.13) in Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 . One can work in either $q$ space to find $v_{c}/L$ or $u$ space to find directly $v_{c}$. The method was used as the substitute of the method (ii) for the smallest values $v_{c}\lesssim 10^{-10}$ we achieved, where Padé approximation become computationally challenging taking more computer time than the calculation of $\tilde{z}$ itself. The absolute accuracy of that method for $v_{c}$ is $\sim 10^{-20}.$

To summarize, we used the first two methods for finding $v_{c}$ for solution with $v_{c}>10^{-8}$ in DP, the second method for solutions with with $v_{c}\gtrsim 10^{-10}$ in quad precision and the third method for solution with $v_{c}\lesssim 10^{-10}$ (with $M_{q}\gtrsim 10^{6}$). We also performed a multiprecision simulations with a variable precision arithmetics with $\sim 200$ digits for selected values of parameters as described in the next Section.

VI Results of numerical simulations

Previous results summarized in Fig. 2 of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 (they are also reproduced in the left part of the curve of Fig. 5) showed a nontrivial dependence of the Stokes wave speed $c$ on the height $H$ with $H$ monotonically approaching the maximum value $H_{max}$ and $c$ approaching a finite value non-monotonically while oscillating with an amplitude that decreases approximately two orders in magnitude every half of such oscillation. Computing Stokes wave solutions in $u$ space as in Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 allowed to resolve about 1.5 of such oscillations (see Fig. 5) while implementing the approach described in this paper (solving in $q$ space) allowed to resolve about 3.5 of such oscillations.

One example of the numerical solution (corresponds to the most extreme wave of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 ) is given in the Introduction. Another example of a less steep wave solution for $c=1.0924$ computed using double precision resulting in $H/\lambda=0.1404429731116977$ and $v_{c}=0.0006925714\ldots$ is given in Fig. 6 in variables $u$ (left panel) and $q$ (right panel) with the corresponding spectra of $\tilde{z}(u)$ and $\tilde{z}(q)$ showed in Fig. 7 (both spectra have only negative components of $k$ since both $\tilde{z}(w)$ and $\tilde{z}(q)$ are holomorphic in $\mathbb{C}^{-}$). Here $M=64536$ on a uniform grid and $M_{q}=1024$ on a nonuniform grid with $L=0.018608751114420542$. It demonstrates that for this particular case one needs 64 times less Fourier harmonics in $q$ space compared to the $u$ space in order to resolve the solution up to DP round-off error. The speed up factor could be roughly estimated as $1/L=(v_{c}/2)^{1/2}$ that becomes significant as we go to lower values of $v_{c}$.

High precision and range of our simulation parameters allow to reveal the asymptotic behavior of Stokes wave as it approaches the limiting form as well as make a comparison with the theory of Stokes wave. We start by analyzing the dependencies of wave speed $c$ and height $H$ on the parameter $\chi_{c}$ for the obtained family of Stokes waves, where

[TABLE]

is the distance to the singularity of a Stokes wave solution to the real line in the axillary space $\zeta=\tan[\frac{w}{2}]$. Notice, that for the highly nonlinear Stokes waves $v_{c}\rightarrow 0$ and $\chi_{c}\simeq v_{c}/2$, while for the almost linear Stokes waves $v_{c}\rightarrow\infty$ and $\chi_{c}\rightarrow 1$. Fig. 8 shows $|c_{lim}-c|$ and $(H_{max}-H)/\lambda$ vs. $\chi_{c}$ for computed Stokes waves in the log-log scale together with the corresponding fitting curves. Here $c_{lim}$ and $H_{max}/\lambda$ are the speed and the scaled height of the limiting Stokes waves, respectively. We use the numerical values $c_{lim}^{GL}=1.0922850485861$ and $H^{GL}_{max}/\lambda=0.1410634839798$ found by I.S. Gandzha and V. P. Lukomsky in Ref. GandzhaLukomskyProcRoySocLond2007 with the claimed accuracy in 11 digits. Comparable accuracy was also achieved in Ref. MaklakovEuroJnlAppliedMath2002 . It is seen from Fig. 8 that $|c_{lim}-c|$ experiences oscillations with their envelope being the excellent fit to the linear law

[TABLE]

Fig. 8. shows that the dependence of $H_{max}-H$ on $\chi_{c}$ at the leading order fits well to the the power law

[TABLE]

while experiencing small oscillations with the vanishing amplitude as $\chi_{c}\rightarrow 0.$ The scaling (57) was proposed in Ref. DyachenkoLushnikovKorotkevichJETPLett2014 from simulations and can be extracted at the leading order from the analytical Stokes wave solution of Section 8 of Ref. LushnikovStokesParIIJFM2016 . Using the scaling (57) we fit the simulation data into the model

[TABLE]

where the constant $\Delta H_{max}$ accounts for the accuracy in the numerical value of $H^{GL}_{max}$ and $A$ is the another fitting parameter. Using the smallest values of $\chi_{c}\lesssim 10^{-10}$ achieved in simulations, we obtained from that fit the estimate $H^{fit}_{max}/\lambda=(H^{GL}_{max}+\Delta H_{max})/\lambda=0.141063483977\pm 10^{-11},$ i.e. $\Delta H_{max}/\lambda=-2.8\cdot 10^{-12}$ which is consistent with 11 digits accuracy of $H^{GL}_{max}$. The highest wave that we computed in QP has $H_{max}^{lowerbound}/\lambda=0.1410634805062790\ldots$ (for $c=1.09228504858750000$) which provides the best lower bound $H_{max}/\lambda$ from our simulations. That lower bound is within $\simeq 3.5\times 10^{-9}$ from $H^{fit}_{max}/\lambda$ which is more than 3 orders in magnitude of improvement compare with the simulations of Ref. DyachenkoLushnikovKorotkevichPartIStudApplMath2016 .

To focus on the corrections beyond the leading order scalings (56) and (57), we plot $|c_{lim}-c|/\chi_{c}$ and $(H_{max}-H)/(\lambda\chi_{c}^{2/3})$ vs. $\chi_{c}$ in Fig. 9. It seen on left panel that the simulation data for $(c_{lim}-c)/\chi_{c}$ are well fit onto the sin-log model

[TABLE]

Here we used the data points from the left-most 1.5 oscillations to find the fitting values $\alpha=0.395,\ \omega_{1}=0.716$ and $\varphi_{1}=2.01$. The points from the steepest waves are the most sensitive to the value of $c_{lim}$ on this plot. Adjusting the value $c_{lim}$ and observing the changes in the plot while assuming that in the limit $\chi_{c}\rightarrow 0$ the proposed sin-log model is valid we estimated that $c_{lim}=1.0922850485861\pm 5\times 10^{-13}$ which is again consistent with 11 digits accuracy of $c_{lim}^{GL}$.

One can compare Eq. (59) with the expression

[TABLE]

which was obtained by M.S. Longuet-Higgins and M.J.H. Fox in Ref. Longuet-HigginsFoxJFM1978 by matched asymptotic expansions. Here $\epsilon:=2^{-1/2}q$ and $q$ is the particle speed at the wave crest in a frame of reference moving with the phase speed $c$. To find $q$ we notice that the complex velocity $V:=v_{x}-\mbox{i}v_{y}$ is given by $V=\Pi_{u}/z_{u}$, where $v_{x}$ and $v_{y}$ are the horizonal and vertical velocities in physical coordinates in the rest frame and $\Pi$ is complex potential which for Stokes wave is given by $\Pi=c(z-w)$ (see e.g. the Appendix B of Ref.LushnikovStokesParIIJFM2016 ). It implies using the analytic solution of Section 8 of Ref. LushnikovStokesParIIJFM2016 that $\epsilon=2^{-1/2}q=2^{-1/2}|V|_{w=0}-c|=2^{-1/2}c|(z_{u}-1)/z_{u}-1||_{w=0}=cP\chi_{c}^{1/3}+O(\chi_{c}^{2/3}),$ where $P\sim 1$ is the constant. Then it is seen that Eqs. (59) and (60) are consistent if we additionally notice that $2.143\ln\epsilon\simeq 0.714\ln\chi_{c}+const$ which is within the accuracy of the numerical value $\omega_{1}=0.716$ in the parameter fit of Eq. (59). In addition, the coefficient $1.1931$ in right-hand side of Eq. (60) is the numerical approximation of Ref. Longuet-HigginsFoxJFM1978 for $c_{lim}^{2}$. Thus Fig. 9 reproduces 3 oscillations of Eq. (60).

The definition of the Stokes wave height results in $H=y(0)-y_{b}$, where $y_{b}:=y(\pm\pi)$. Then the analytical Stokes wave solutuon of Section 8 of Ref. LushnikovStokesParIIJFM2016 implies that the next correction beyond the leading order model (58) is given by

[TABLE]

where $y_{b,lim}$ is the value of $y_{b}$ for the limiting Stokes wave (we approximate it by the most extreme wave in our simulations). To check the accuracy of Eq. (61) we divided it by $\lambda\chi_{c}$ and compared the right-hand side with the left-hand side on right panel of Fig. 10 showing excellent agreement. Oscillations both in $c^{2}(\chi_{c})$ (as seen on left panel of Fig. 9) and in $y_{b}(\chi_{c})$ are comparable in amplitude both contributing to that agreement. Inspired by the model (59), we fit the data $H(\chi_{c})$ to the following model $H^{GL}_{max}+\Delta H_{max}-H(\chi_{c})-A\chi_{c}^{2/3}\simeq B\chi_{c}\cos[\omega_{2}\ln(\chi_{c})+\varphi_{2}]$, where $\Delta H_{max},A,B,\omega_{2},\varphi_{2}$ are unknown constants. Using the data only from the large oscillation at the smallest $\chi_{c}$ on right panel of Fig. 10, we obtained that $\Delta H_{max}/\lambda=1.3\times 10^{-12},A=0.1289,B=0.120,\omega_{2}=0.702,\varphi_{2}=-2.00$. This model fit is shown on right panel of Fig. 10 by the solid curve. Notice, that $|\omega_{1}|\approx|\omega_{2}|$. We expect that if the suggested models are correct then $|\omega_{1}|$ should be equal to $|\omega_{2}|$ as $\chi_{c}\rightarrow 0$.

It is also instructive to relate the physical variables $H_{max}-H$ and $c_{lim}-c$ directly bypassing the use of $\chi_{c}$. Via Eq. (58) we approximate $\chi_{c}$ at the leading order through $H^{GL}_{max}-H.$ We plug in that approximation into Eq. (59) to obtain at the leading order that

[TABLE]

Left panel of Fig. 10 shows $(c_{lim}-c)\lambda^{3/2}/(H_{max}^{GL}-H)^{3/2}$ vs. $(H_{max}^{GL}-H)/\lambda$ together with the model $(c_{lim}-c)\lambda^{3/2}/(H_{max}-H)^{3/2}\simeq D\cos[\omega_{3}\ln[(H^{GL}_{max}-H)/\lambda]+\varphi_{3}].$ The fitting constants of the model are $D=8.55,\ \omega_{3}=1.076$ and $\varphi_{3}=-2.04$ being consistent with the leading order expression (VI).

VII Generalization of the conformal map to resolve multiple singularities

Assume that we aim to approximate a general $2\pi$ periodic function $f(u),\ u\in\mathbb{R}$ which has multiple complex singularities in its analytical continuation into the complex plane $w=u+\mbox{i}v.$ Example of such function is higher order progressive waves, which have more than one different peaks per $2\pi$ spatial period ChenSaffmanStudApplMath1980 . We would like to efficiently approximate $f(u)$ thought the Fourier series in the new variable $q(u).$ A generalization of the conformal transformation (6) to take into account these multiple complex singularities of the function $f(u)$ is given by

[TABLE]

where $L_{j}>0$, $\beta_{j}>0$ and $-\pi<u_{j}<\pi$ are real constants. A condition $\sum\limits_{j=1}^{N}\beta_{j}=1$ ensures that $-\pi\leq q\leq\pi.$ The Jacobian of Eq. (63) given by

[TABLE]

is positive-definite. It ensures that Eq. (63) is one-to-one map between $-\pi\leq q<\pi$ and $-\pi\leq u<\pi$.

Consider an analytical continuation of $f(u)$ into the complex plane $w=u+\mbox{i}v$ and choose $N$ complex singularities of $f(w)$ which are closest to the segment $-\pi\leq u<\pi$ of the real line $w=u$. Then $u_{j}$ can be chosen either close or equal to the projections of the positions of complex singularities of $f(w)$ into the real line $w=u.$ Similar to the case $N=1$ described in previous sections, the constants $L_{j}$ can be chosen to move the complex singularities of $f(q)$ further away from the real line of $q$ compared with locations of singularities of $f(w)$. Then the uniform grid in $q$ corresponds to the nonuniform grid in $u$ concentrating at neighborhoods of each $N$ singularities of $f(w)$ thus greatly improving the approximation of $f(u)$ by FT in $q$. The constants $\beta_{j}>0$ can be either chosen equals, $\beta_{j}=1/N,j=1,2,\ldots,N,$ or nonequal (to provide stronger weights for most dangerous complex singularities). Note that $u_{j}$ can be chosen the same for several values of $j$ (but with different values of $L_{j})$ thus simultaneously resolving several singularities (or several segments of the vertical branch cut) located along the same vertical line in $w$ plane.

Another transformation $q(u)$ to resolve multiple singularities is given by a generalization of Eq. (10) as follows

[TABLE]

which is explicit expression for $u(q)$ but implicit for the inverse $q(u).$ Here $L_{j}>0$, $\alpha_{j}>0$ and $-\pi<q_{j}<\pi$ are real constants. A condition $\sum\limits_{j=1}^{N}\alpha_{j}=1$ ensures that $-\pi\leq q\leq\pi.$ Similar to Eq. (63), Eq. (65) has a strictly positive Jacobian which ensures one-to-one map between $-\pi\leq q<\pi$ and $-\pi\leq u<\pi$.

VIII Conclusion

In conclusion, we found the new transformation (6) which allows to move the lowest complex singularity $w=\mbox{i}v_{c}$ of the function $f(u)$ away from the real line thus greatly improving the efficiency of FT of that function in the new variable $q$. Number of Fourier modes needed to reach the same precision of approximation of $f$ in $q$ variable is reduced by the factor $\sim v_{c}^{1/2}$ for $v_{c}\ll 1$ compare with FT in $u$. We showed that the new transformation (6) is consistent with the dynamics of two dimensional Euler equation with free surface. We demonstrated the effciency of Eq. (6) for simulations of Stokes wave by improving the numerical performance in many orders of magnitude. It allows to reveal the details of the oscillatory behaviour of the parameters of Stokes wave as it approaches the limiting wave. We suggested the generalizations (63) and (65) of Eq. (6) to resolve multiple singularities.

The work of D.S. and P.L. was partially supported by the National Science Foundation Grant DMS-1412140.