{'fraction_non_alphanumeric': 0.04193679092382496, 'fraction_numerical': 0.01819286871961102, 'mean_word_length': 4.61442220200182, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient computation of capillary-gravity generalized solitary waves

This paper is devoted to the computation of capillary-gravity solitary waves of the irrotational incompressible Euler equations with free surface. The numerical study is a continuation of a previous work in several points: an alternative formulation of the Babenko-type equation for the wave profiles, a detailed description of both the numerical resolution and the analysis of the internal flow structure under a solitary wave. The numerical code used in this study is provided in open source for those interested readers.


Introduction
The formation and dynamics of capillary-gravity waves on the surface of incompressible fluids constitute a research topic of permanent interest from experimental, theoretical and computational point of view [20,36,42,43]. These solutions have been also studied using analytical asymptotic methods in various model equations (such as KdV5), see [7][8][9] and also the contributions to the edited volume [40]. This paper presents in details an efficient numerical algorithm for computing various steady capillary-gravity solitary waves for the irrotational Euler equations with a free surface. In particular, this algorithm can be used to compute generalized multi-hump solitary waves [19]. In addition to the detailed algorithm, the present paper provides some physical characteristics of the generalized solitary waves that are not described in [19]. The highlights can be assembled in three points.
The first one concerns the mathematical formulation of the equations for the wave profiles. Different from other approaches used in the literature to reformulate the Euler system on the free surface, [24,38,42], the one employed in [19] is based on the conformal mapping technique [37], in order to transform the system into an integro-differential equation of Babenko type [3]. (This was previously developed for gravity solitary waves in [17,18,21].) As mentioned in [19], one of the advantages of this formulation is the preservation of the type of nonlinearity by the transformation [32]; it has implications on the convergence of the numerical method. As in [19], the numerical results on steady capillary-gravity waves presented here are obtained with a pseudo-spectral discretisation of the corresponding periodic problem and the application of the Levenberg-Marquardt (LM) algorithm, [26,28,30], to the resulting discrete systems for the Fourier components of the approximation.
A second highlight of the present paper is then a detailed description of the numerical procedure and its implementation to the model formulation. The performance of the method is shown through several numerical experiments computing different traveling waves.
Finally, the emergence of generalised solitary waves shown in [19], under gradually increasing values of the Bond number is used in a third highlighted point of the present study, where some properties of various physical fields under these waves are emphasised.
The manuscript is organised as follows. In Section 2, we summarise the application of the conformal mapping technique to the original Euler equations, explained in [19] in more detail. Section 3 is devoted to a detailed description of the numerical method while the main numerical experiments are presented in Section 4. Finally, the main conclusions and perspectives are outlined in Section 5. With the aim of involving a wider audience, the numerical code used in computations below is available to download as open source [22]. Thus, the claims made in this study can be easily verified by interested readers.

Mathematical model
Our starting point is the mathematical model for a potential flow induced by a solitary wave described in [19], which is summarised here. The physical assumptions involve an inviscid, homogeneous fluid in a horizontal channel of constant depth; the channel is modelled, above, by an impermeable free surface where the pressure is equal to the surface tension due to capillary forces, and bounded below by a fixed impermeable horizontal seabed. The mathematical description in Cartesian coordinates moving with the wave (with x as the horizontal coordinate and y the upward vertical one) is sketched in Figure 1: d stands for the mean depth of the channel, the bottom is posed at y = −d while y = η(x) denotes the free surface elevation from the mean water level at y = 0. The zero mean level condition of the free surface is redefined in order to cover the computations of generalised solitary waves (see [19] for details).
Let φ and ψ be the velocity potential and the stream function, respectively, and consider the complex potential f ≡ φ + iψ. They define a conformal mapping where ψ s and ψ b are the constant traces of ψ at the upper and lower boundaries, respectively, with −c standing for the mean flow velocity and satisfying The conformal mapping (2.1) transforms the fluid domain into the strip −∞ α ∞, −d β 0, where α ≡ Re(ζ) and β ≡ Im(ζ). On this domain, the free surface η can be modelled by a nonlocal Babenko-type equation of the form [19] In (2.3), B is the Bernoulli constant, τ is a surface tension parameter, K is an integration constant and C is a pseudo-differential operator with Fourier symbol The constant K is determined from the mean level condition on η. In the limit τ → 0 the Babenko formulation for pure gravity solitary waves derived in [18,21] is recovered. K = 0 for classical solitary waves because η → 0 as x → ±∞ and the wave mass´∞ −∞ η dx is finite. For generalised solitary waves K is not zero, in general, and it can be determined by the far field condition [19]. For the numerical resolution of (2.3), we start with initial guesses such that η → 0 as x → ±∞, so K = 0 initially. When oscillatory tails emerge during the iterations, K is no longer zero. In order to compute K dynamically at each iteration, one must know the nature (and location) of the tail, that is unknown a priori and thus only an a posteriori computation of K is practically reasonable after the algorithm has converged. The generalised solitary waves described in [19] have tails with one fundamental frequency, and an ad hoc procedure is then used to compute K. However, the possibility of more general tails has not been ruled out, so a generalised algorithm for computing K cannot be given until the nature of the tail is known. Therefore, we compute the numerical solutions of (2.3) with K = 0, that is always possible without loss of generality if one releases the mean level condition (this is equivalent to redefining the mean water depth).
With K = 0, the equation (2.3) is reformulated as follows. Applying the operator C −1 to both sides of (2.3) and, with T = C −1 d dα , we can write where we separated the linear and nonlinear (quadratic, rational) terms. The Fourier symbol of the linear operator (2.6) is In (2.8) the dispersion relation for infinitesimal capillary-gravity waves can be recognized [20,42].

The numerical method
Here, we describe a numerical procedure for discretising and solving the generalised Babenko equation (2.3).

Pseudospectral discretisation
For Λ > 0 sufficiently large, the periodic problem associated to (2.5) on (−Λ, Λ) is discretised with Fourier collocation techniques. For N > 1 and on an uniform grid α j = −Λ + jh, h = 2Λ/N , j = 1, . . . , N , the values η(α j ) of the solution of (2.5) are approximated by In (3.1)-(3.3), D α stands for the pseudo-spectral differentiation matrix [10,13,14,41], while C h and T h are discrete versions of the operators C and T , respectively. They are constructed as D α by using the discrete Fourier transform matrix with diagonal entries given by the coefficients displayed in (2.4). The operator T h is defined in the same way.) Thus, the discrete system (3.1)-(3.3) is implemented in the Fourier space. Finally, the nonlinear terms are computed by using the Hadamard product of vectors in C N . For these steady computations the use of the anti-aliasing rule was not necessary [6,10,43]. A major technical difficulty with nonlocal (generalised) solitary waves is that the phase of tail oscillations is sensitive to the position of the boundary. The possible tails are quantised by the computational box and since we realise the mean level condition during the computation, our method finds the corresponding solution for a given computational domain by making the adjustments in the mean water level.

The Levenberg-Maquardt algorithm
The system (3.1)-(3.3) was tentatively treated with several techniques. When τ = 0 (pure gravity waves) this was successfully solved using the Petviashvili scheme [18,21]. Note that when τ = 0, the nonlinear term (2.7) is homogeneous with degree two, while the Fourier symbol (2.8) associated to (2.6) is positive for all k when B > gd. These two properties are key arguments to explain this success [39] (see also [1,2]). However, they are not retained when τ = 0 and the various theories about the behaviour of the Petviashvili scheme cannot be applied here.
For the general case τ = 0, several variants of the Newton method can emerge as alternatives. The reasons for the failure of the classical implementations of this method [11] were taken into account. In our case, the translational invariance of the solitary waves implies that the corresponding Jacobian of F h in (3.1) is rank-deficient (in practice, it is nearly rankdeficient). This recommends to avoid not only the direct application of Newton's method but also the resolution in a least square sense with line search (the Gauß-Newton method). Instead of all this, the system (3.1)-(3.3) was numerically solved in a nonlinear, least square sense but with a trust region formulation. Among the standard techniques presented in the literature (see, e.g., [25] and the references therein), the most robust results were obtained employing the so-called Levenberg-Marquardt (LM) algorithm [26,30]. This is one of the most widely used methods for data-fitting nonlinear problems [28]. For the reader interest, a brief description of it in this context is given now.
· stands for the Euclidean norm in R N . As established in [33], the LM algorithm can be formulated as a damped Gauß-Newton method but combined with a trust region strategy. The Gauß-Newton method is a variant of the Newton method with line search, where the search direction p (ν) = p GN (ν), at each Newton step ν, is obtained by solving where r (ν) , J (ν) stand for the residual vector and its Jacobian at the ν-th iteration η h , respectively. (In our case, the Jacobian was approximated with finite differences, leading in fact to a quasi-Newton method.) Instead of (3.4), the LM algorithm computes p (ν) = p (ν) LM as a minimiser of the model function subject to the condition p ∆ (ν) , where ∆ (ν) > 0 denotes the corresponding radius of the trust region where the problem is constrained to ensure convergence of the algorithm (see [29,33,35] for details). Any solution p (ν) is characterised by the existence of a scalar λ (ν) 0 (the damping parameter) satisfying where I N is the N × N identity matrix. Condition (3.7) simply states that at least one of the nonnegative quantities, λ (ν) or ∆ (ν) − p (ν) , must vanish.) The literature contains several strategies for the computation of the damping parameter [30,34,35]; one of those proposed in the second reference is standard and was considered here by using the MATLAB Optimization Toolbox [12]: from an initial λ (0) , when the step is successful, λ (ν+1) is updated as λ (ν) /10; otherwise, λ (ν+1) is set to 10λ (ν) . Finally, the local convergence of the resulting method is linear [35], but close to quadratic [29], when the damping parameter is sufficiently small.

Numerical results
In numerical computations herein below, the parameters given in Table 1 will be used. The first two concern the model and the rest is associated to the numerical procedure described above. The definition of Froude (Fr) and Bond (Bo) numbers is the one given in [19]. In most of the cases, our implementation of the method will make use of the natural numerical continuation in the Bond number Bo, as a way to improve the performance and to accelerate the convergence. The numerical study is focused on classical, generalised and multi-pulse solitary waves, in several physical regimes.

Parameter
Value 20 Tolerance parameter on the increment norm, δ 10 −14 Table 1. Physical and numerical parameters used for the computation of generalised solitary waves by continuation in the Bond number Bo. The tolerance δ is the stopping criterium in the nonlinear solver.

Classical solitary wave of depression
The first group of numerical experiments concerns the generation of classical solitary waves of depression. With this name, we mean aperiodic solutions of (2.5) decaying to zero at infinity (in contrast with the generalised solitary waves, which will be considered later on) and with negative amplitude.
For supercritical Bond numbers Bo > 1/3 and suitable values of the Froude number Fr < 1, the existence of isolated solitary waves of depression is known, analytically and numerically (see [16] and references therein). A first check of the code is the computation of one of these waves. This is shown in Figure 2, corresponding to the values Bo = 0.45 and Fr = 0.6. The initial iteration of the procedure was chosen simply as a negative localised bump. There is no need to perform any post-processing, since the solution decays exponentially to zero, i.e., η(x) → 0 as |x| → ∞. So, the initial values of Fr and Bo parameters correspond perfectly to the fully converged solution.  The classical character of this computed wave is confirmed by Figure 2 (right), which shows that the associated orbit in the phase portrait (computed via the spectral approximation) is homoclinic at infinity to the origin.
The exponential decay has also been numerically checked. To this end and for several profiles, the rightest part of the computed wave was fitted by an exponential function of the form a e bx . By way of illustration, the fit corresponding to the solitary waves for Bo = 0.46, Fr = 0.6 and Bo = 0.46, Fr = 0.9 (computed on a longer interval) is shown in Figure 3. For the first case, the fitting coefficients with 95% confidence bound (computed with the MATLAB toolbox) are a = −0.9186, b = −1.414. The goodness of the fit was ensured by a SSE (sum of squares due to error) parameter of 2.024E−24 a R-squared parameter of 1 and a RMSE (root mean squared error) parameter of 9.526E−14.
The convergence of the method is illustrated by the following experiments. Figure 4 (left) displays the residual error as function of the iterations for the solitary wave of depression shown in Figure 2. This residual error is measured, at the n−th iteration η (n) h , as where the Euclidean norm is considered and F h {·} is given by (3.1). The vertical scale is logarithmic and the results confirm the convergence. On the other hand, Figure 4 (right) shows, in log-log scale, the relation between two consecutive residual errors. This serves as an estimate of the order of convergence. By fitting the logarithmic data to a straight line, the corresponding slope (with 95% confidence bound) is 1.179. When only errors up to 10 −7 are considered, the corresponding fitting line has a slope of approximately 2.156. Then, up to this level of errors, Figure 4 (right) suggests a quadratic order of convergence, becoming linear below this error tolerance. Moreover, the residual magnitude is known to be a rather pessimistic estimation of the accuracy in pseudo-spectral methods [10]. Consequently, the convergence of the residual (4.1) is a robust diagnostics of the algorithm convergence. As the Froude and Bond numbers approach the limiting values, the maximum negative excursion of the profiles of depression decreases and the code becomes less accurate. This was observed, for example, in the goodness of fit with an exponential when estimating the decay.
It should be noted that one may also study, by numerical means, the dependence of the asymptotic decay on the Froude and Bond numbers. This is illustrated in Figure 5  This subsection ends by visualising the internal hydrodynamics under classical solitary waves. This is illustrated by considering the computed wave of Figure 2. Thus, the velocity potential and the stream function are depicted in Figure 7, the total and dynamic pressures are represented in Figures 8, while horizontal and vertical speed and acceleration distribution are displayed in Figures 9 and 10, respectively. All are implemented by using the Cauchy integral formula, which allows to compute all physical fields of interest in the bulk of the fluid, [21]. The traces of the velocity potential, stream function, horizontal and vertical velocities are depicted in Figure 6.

Multi-modal classical solitary waves
Another type of solutions, which can be computed by our code, consists of multi-modal solitary waves. The existence of these solutions in the KdV5 equation was shown analytically and numerically in [15]. We refer to [19] for more details on these solutions in case of the full Euler equations.

Generalised solitary waves
The question of existence of classical solitary waves of elevation is not solved, to our knowledge, although some references suggest a negative answer [16]. For Bo < 1/3 and Fr > 1, what is known is the existence of generalised solitary waves, i.e. solitary waves that are homoclinic to exponentially small amplitude oscillatory waves. Furthermore, the amplitude of the corresponding asymptotic oscillations are of order less than for some s ∈ (0, π) and where ω = 0 satisfies the equation ω cosh(ω) = (1 + ω 2 Bo) sinh(ω) [27].
Here, we focus on the numerical generation of this kind of waves. Our first experiments concern the range Bo < 1/3 and Fr > 1, illustrating thus the influence of the capillarity (smaller or larger values of Bo) on the computations. Taking as initial guess a third-order asymptotic solution for solitary-gravity waves [23] (for Fr > 1 and 0 < Bo ≪ 1/3), the method has been run in the case of strong capillarity (Fr = 1.15 and Bo = 0.22). The resulting numerical profiles are shown in Figure 11. Once the solution is computed, the Froude and Bond numbers have to be actualised to take eventually into account for the periodic tail. The methodology of post-processing is explained in [19] and the effective values of the Froude Fr ∞ and Bond Bo ∞ numbers are provided in the corresponding captions in Figure 11. The corresponding phase plots are displayed in Figure 14. They show the characteristic behaviour of the orbit of a generalised solitary wave profile, homoclinic at infinity to a periodic wave.
From the computed profiles, some additional information concerning the oscillatory tails can be obtained. For instance, part of the profile can be taken and fitted to a sinusoidal wave in order to estimate the corresponding amplitude and wavelength. This fit is illustrated in Figure 13 which shows the data of part of the tail of Figure 11 (d) and two different sinusoidal curves. In Figure 13(a) the data are fitted by a Fourier sum (with two terms) while in Figure 13(b) this is made by a sum of sinusoidal functions with six different wavelengths.
x / d The corresponding coefficients and goodness of fit parameters are displayed in Table 2. In the first case, the best choice among different Fourier sums corresponds to taking two Fourier modes and leads to an approximate fundamental frequency of f = w/(2π) = 1.2634. On the other hand, in the second fitting curve we have β 1 = ω, β 6 = 2ω and the influence of the rest of the frequencies is limited by the amplitudes of the corresponding sinusoidal terms. Thus the two fits give the same essential information about the oscillatory tail: an approximate amplitude of 0.02517 and wavelength of λ = 1/f = 0.7915. The convergence process is illustrated as in the previous Section. Figure 12 displays the logarithm of the residual error (4.1) as function of the number of iterations. Note that the method performs better in the case of the nonlinear wave with weaker capillary effects, when the oscillations at infinity of the wave is of smaller amplitude, see (4.2). The order of convergence is suggested by the Figure 12 (right) (where the ratio between two consecutive errors is shown in log-log scale) in the same sense as in Figure 4 Figures 15-18 for the same profile, with the velocity potential (Fig. 15 left) and stream function (Fig. 15 right), total (Fig. 16 left) and dynamic (Fig. 16 right) pressures, horizontal (Fig. 17 left) and vertical (Fig. 17 right) velocities and accelerations (Fig. 18). The traces of the velocity potential, stream function, horizontal and vertical velocities on the free surface are depicted in Figure 19.
The dependence of the amplitude of the computed oscillatory tail on the Bond number can be observed in Figure 20. This shows the amplitude of the fitting sinusoidal wave as function of the Bond number for the profiles of Figure 11.

Multi-hump generalised solitary waves
As shown in [19], multi-hump generalised solitary waves (that is, homoclinic orbits with several loops near a resonance) can be computed using separated bumps as initial iteration. They are illustrated in Figure 21. We refer to [19] for other interesting examples.

Conclusions and perspectives
This study is a continuation of a previous work [19]. The goal here was to explain in details the numerical techniques used in [19] and to provide some additional insights.   Table 2. Coefficients and goodness of fit for the fitting curves of Figure 13: f (x) = a 0 + a 1 cos(ωx) + b 1 sin(ωx) + a 2 cos(2ωx) + b 2 sin(2ωx); f (x) = 6 j=1 α j sin(β j x + γ j ).
The Babenko equation was discretised with a Fourier-type pseudo-spectral method [10]. The resulting discrete nonlinear and nonlocal system was solved using the Levenberg-Marquardt (LM) method. The reason to use this technology is to overcome the drawbacks of the direct application of Newton's (and quasi-Newton's) methods in this and related problems [11]. The rank-deficient character of the Jacobian, due to the translational invariance of the solitary waves, was the main reason for this choice. Using this formulation, we succeeded to compute, by natural continuation in the Bond number Bo, the generalised solitary waves of elevation for Bo < 1/3 [16]. Above the critical Bond number Bo = 1/3, we found the classical localised solitary waves of depression which propagate with subcritical speeds Fr < 1, in agreement with the predictions of the KdV5 model. Here, we completed the numerical results with additional information provided by the computations: (i) on the asymptotic decay of the classical solitary waves; (ii) on the oscillatory tails of the computed generalised solitary wave profiles and (iii) on the internal hydrodynamics structure. In the latter case, we unveil what happens with various physical fields (velocities, pressure, accelerations) under a generalised solitary wave. This was made using Cauchy-type integral representations, available when we take the full advantage of the conformal mapping technique, [21].
The present work opens several research directions. First of all, the variability of solutions arising in several directions is a subject to be explored further. Then, the dynamics and   Table 2.
stability of these solutions have to be studied. This requires the corresponding unsteady formulation along with an efficient discretisation procedure. The generalised solitary waves are homoclinic to Stokes waves. We conjecture that most of the solutions computed in this study might be unstable, since the periodic (gravity) Stokes wave is unstable [31], mainly because of the Bespalov-Talanov-Benjamin-Feir instability [4,5]. The role of the surface tension has to be studied properly. Can it suppress some instabilities? Finally, the corresponding Babenko-type formulation for periodic capillary-gravity waves has  to be developed as well with the same questions about their dynamics and stability. It will be slightly different from the Babenko equation derived in [19] due to some extra parameters, which characterise periodic waves. The numerical procedure for the computation should also change in some suitable way. The code employed to obtain all the numerical results presented in this manuscript (but also the results of our previous work [19]) is freely available to download at the following URL address [22]: • https://github.com/dutykh/BabenkoCG/ We invite the readers to use this code to compute and to study generalised solitary waves.