T K Papathanasiou
email: theodosios.papathanasiou@brunel.ac.uk

Ch E Papoutsellis

G A Athanassoulis

Semi-explicit solutions to the water-wave dispersion relation and their role in the nonlinear Hamiltonian Coupled-Mode theory

Keywords: Dispersion relation, nonlinear water waves, Hamiltonian coupledmode theory, multimodal techniques, root approximation, Newton-Raphson iterations

, provides an efficient new approach for solving fully nonlinear water-wave problems over arbitrary bathymetry. This theory exactly transforms the free-boundary problem to a fixedboundary one, with space and time varying coefficients. In calculating these coefficients, heavy use is made of the roots of a local, water-wave dispersion relation with varying parameter, which have to be calculated at every horizontal position and every time instant. Thus, fast and accurate calculation of these roots, valid for all possible values of the varying parameter, are of fundamen-tal importance for the efficient implementation of HCMT. In this paper, new, semi-explicit and highly accurate root-finding formulae are derived, especially for the roots corresponding to evanescent modes. The derivation is based on the successive application of a Picard-type iteration and the Householder's root finding method. Explicit approximate formulae of very good accuracy are obtained, and machineaccurate determination of the required roots is easily achieved by no more than three iterations, using the explicit forms as initial values. Exploiting this procedure in the HCMT, results in an efficient, dimensionallyreduced, numerical solver able to treat fully non-linear water waves over arbitrary bathymetry. Applications to four demanding nonlinear problems demonstrate the efficiency and the robustness of the present approach. Specifically, we consider the classical tests of strongly nonlinear 1 steady wave propagation and the transformation of regular waves due to trapezoidal and sinusoidal bathymetry. Novel results are also given for the disintegration of a solitary wave due to an abrupt deepening. The derived root-finding formulae can be used with any other multimodal methods as well.

Introduction

The simulation of nonlinear water waves over variable bathymetry is a challenging task that involves various numerical techniques and extensive computations.

The complications due to the presence of the unknown free-surface elevation and the varying bathymetry are usually treated either by perturbation techniques or by direct numerical methods (DNM). The first approach, being in use for more than a century, has led to a plethora of approximate models, mainly classified as Boussinesq-type or Serre-Green-Naghdi-type models [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF], [START_REF] Madsen | A Boussinesq-type method for fully nonlinear waves interacting with a rapidly varying bathymetry[END_REF], [START_REF] Karambas | Boussinesq Model for Weakly Nonlinear Fully Dispersive Water Waves[END_REF], [START_REF] Bonneton | Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes[END_REF], [START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF], [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF], [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], [START_REF] Memos | On Higher-Order Boussinesq-Type Wave Models[END_REF]. DNM include finite-difference methods [START_REF] Bingham | On the accuracy of finite difference solutions for nonlinear water waves[END_REF], finite-element methods [START_REF] Gagarina | Variational space-time (dis) continuous Galerkin method for nonlinear free surface water waves[END_REF], [START_REF] Brink | Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves[END_REF] and boundary-element (BEM) methods [START_REF] Grilli | A fully non-linear model for threedimensional overturning waves over an arbitrary bottom[END_REF], [START_REF] Fructus | An efficient model for three-dimensional surface wave simulations. Part I: Free space problems[END_REF]. Perturbative approaches result in numerically efficient, dimensionally reduced models, limited to weakly nonlinear phenomena, while DNM are able to accurately treat strongly nonlinear/dispersive problems at the expense of efficiency, because of their high computational cost.

Recently, [START_REF] Athanassoulis | New form of the Hamiltonian equations for the nonlinear water-wave problem, based on a new representation of the DtN operator, and some applications[END_REF], see also [START_REF] Papoutsellis | A new efficient Hamiltonian approach to the nonlinear water-wave problem over arbitrary bathymetry[END_REF], proposed a new formulation of the exact waterwave problem over arbitrary (smooth) bathymetry, providing both dimensional reduction and high accuracy, comparable with that ensured by DNM. This formulation is based on the exact semi-separation of variables in the instantaneous (non-canonical) fluid domain, established in [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with nonplanar boundaries[END_REF], referred subsequently as AP17. The wave potential Φ = Φ(x, z, t), where x = (x 1 , x 2 ) is the horizontal position in the fluid domain, is expanded in a rapidly convergent series of the form Φ = ϕ n (x, t)Z n (z ; η, h), where the local vertical functions Z n are defined in the varying interval [-h(x), η(x, t) ], delimited by the local depth h(x) and the local, instantaneous free-surface elevation η(x, t). Introducing this expansion in Luke's variational principle, we obtain a convenient system of two nonlinear Hamiltonian evolution equations with respect to η(x, t) and the surface potential ψ(x, t) (see Eqs. [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF]), containing a non-local coefficient field, accounting for the substrate kinematics. The latter is determined by a linear, coupled-mode system of horizontal differential equations with variable coefficients, dependent also on the unknown free-surface elevation (see Eqs. [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF]). The last problem can be solved at each discrete time t using the free-surface elevation prediction from the previous time step, providing a versatile and numerically efficient substitute for the Dirichlet-to-Neumann (DtN) operator, introduced by Craig & Sulem [START_REF] Craig | Numerical Simulation of Gravity Waves[END_REF]. This approach is called subsequently Hamiltonian Coupled-Mode Theory or System (HCMT or HCMS), according to the context. In comparison with the nonlinear consistent coupled-mode system, derived previously in [START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF], the HCMS is more efficient and more accurate, since a double-series term appearing in the former has been summed analytically in the latter.

The approximations involved in the numerical implementation of the HCMS, apart from the inevitable discretization of the free-surface elevation η(x, t) and the depth function h(x), constitute of the truncation of the series expansion for the wave potential, and the calculation of the coefficients of the substrate system (see Eqs. ( 8) -( 9)). Since the series expansion is uniformly and rapidly convergent, even for strongly deformed boundaries (see AP17), the only critical issue is the accurate and fast calculation of the coefficients, which have to be calculated at any point of the spatial and temporal discretization. These coefficients are defined as integrals of the vertical basis functions, along local vertical intervals throughout the fluid.

The vertical basis functions are constructed by extending an L 2 -basis to an H 2 -basis (Sobolev space basis) as explained in AP17. The L 2 -basis may be constructed by means of the eigenfunctions of any regular Sturm-Liouville problems, defined in the local vertical intervals (-h, η) . Making the plausible choice to consider the Sturm-Liouville problem corresponding to linear waterwave problem in the strip (-h, η) , results in the usual eigenfunctions, with eigenvalues defined through the dispersion relation of linear water waves with a frequency parameter varying in space and time. Since the coefficients of the substrate kinematical problem, Eq. ( 8), are eventually expressed analytically in terms of these eigenvalues, the accurate and robust determination of the latter is of fundamental importance for the numerical implementation of the HCMS. Note that, in the computation of a demanding nonlinear problem we need to evaluate a number of 5 -8 eigenvalues for 10 8 (in 2D cases) up to 10 11 (in 3D cases) times. Thus, ideally, we would like to have explicit formulae or, at least, to ensure high accuracy with a few (1 -3) Newton-Raphson iterations. Even more essential is the robustness of the root finder, since in the case of non-convergence at any specific point-time the solution procedure will be stopped. Attempts to determine explicit approximate formulae for the solution of the linear water-wave dispersion relation goes back to 50s, aiming -at that timemainly to obtaining a convenient formula for the wavelength corresponding to a given frequency and depth. See, among others, [START_REF] Eckart | Propagation of gravity waves from deep to shallow water[END_REF], [START_REF] Olson | An explicit expression for the wavelength of a gravity wave[END_REF], [START_REF] Hunt | Direct solution of wave dispersion equation[END_REF], [START_REF] Nielsen | Explicit formulae for practical wave calculations[END_REF], [START_REF] Chen | Iterative and Padé solutions for the waterwave dispersion relation[END_REF], [START_REF] Wu | Wave numbers of linear progressive waves[END_REF], [START_REF] Mckee | Calculation of evanescent wave modes[END_REF], [START_REF] Fenton | On calculating the lengths of water waves[END_REF], [START_REF] Newman | Numerical solutions of the water-wave dispersion relation[END_REF], [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF]. More recently, several authors have focused on obtaining very accurate expressions for the real root, corresponding to the propagating mode. See e.g. [START_REF] Guo | Simple and explicit solution of wave dispersion equation[END_REF], [START_REF] Beji | Improved explicit approximation of linear dispersion relationship for gravity waves[END_REF], [START_REF] Simarro | Improved explicit approximation of linear dispersion relationship for gravity waves: Another discussion[END_REF], [START_REF] Vatankhah | Improved explicit approximation of linear dispersion relationship for gravity waves: a discussion[END_REF]. An interesting observation is that, among the aforementioned plethora of studies, only few refer in detail to the roots associated with evanescent modes. These include basically the analysis of [START_REF] Chen | Iterative and Padé solutions for the waterwave dispersion relation[END_REF], [START_REF] Mckee | Calculation of evanescent wave modes[END_REF] and [START_REF] Newman | Numerical solutions of the water-wave dispersion relation[END_REF], [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF]. This fact is partially justified by the physical significance of the propagating mode, linked to the real root. However, several recent advances regarding the modelling and solution of challenging ocean and coastal engineering problems necessitate the very accurate and effective (in terms of computational time) determination of a number of evanescent modes as well. Apart from our HCMT, which is a fully nonlinear approach needing a huge number of accurate evaluation of several eigenvalues, various other formulations also rely on the evaluation of several eigenvalues as, for example, the extended mild-slope equations [START_REF] Massel | Extended refraction-diffraction equation for surface waves[END_REF], [START_REF] Porter | Extensions of the mild-slope equation[END_REF], [START_REF] Chamberlain | Multi-mode approximations to wave scattering by an uneven bed[END_REF], and the consistent coupled-mode systems [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF], [START_REF] Belibassakis | A coupledmode model for the refraction-diffraction of linear waves over steep threedimensional bathymetry[END_REF], [START_REF] Belibassakis | A coupledmode model for water wave scattering by horizontal, non-homogeneous current in general bottom topography[END_REF], [START_REF] Papathanasiou | An Efficient Coupled-Mode / FEM Numerical Method for Linear Wave Propagation[END_REF].

In this study, new explicit and semi-explicit (requiring 1 to 3 iterations), highly accurate formulae for the imaginary roots of the water-wave dispersion relation are presented and analysed. These formulae occur from the iteration of recursive numerical root-finding schemes, especially fitted for the nonlinear equation under consideration. The explicit ones provide an accuracy of 10 -5 or better for all eigenvalues and for all values of the varying parameter, and can be used directly in the solvers of linear and non-strongly nonlinear water-wave problems. The semi-explicit ones provide machine accuracy of 10 -15 with two or three iterations, and are suitable for demanding long-time simulations using the HCMS. This suitability is clearly demonstrated by numerical simulations implemented by means of HCMS and the above root-finding methods, in three nonlinear problems; the propagation of highly nonlinear travelling waves over flat bottom, the transformation of a regular wave by a submerged trapezoidal bar [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF], [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF], and the reflection of nonlinear waves due to a sinusoidal bottom patch [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF]. After establishing its effectiveness and accuracy, the present method is utilized in studying the transformation (disintegration) of a solitary wave passing over an abrupt deepening. The phenomenon revealed in this case is studied herein for the first time, to the best of our knowledge.

The paper is organised as follows: in Section 2, the HCMT for the evolution of nonlinear water-waves over arbitrary bathymetry is briefly presented, and its efficient numerical implementation is discussed in Section 3. In Section 4, the root approximation strategy is developed and analysed. Two iterative procedures (one of second and one of third order) are formulated and error estimates are derived. The procedures for obtaining explicit approximate formulae for the roots, and a thorough analysis of the performance of these formulae either as explicit approximations or as improved initializations of iteration schemes, are the subject of Section 5. Section 6 is devoted to the investigation of the performance of the new formulae in the simulation of the three benchmark nonlinear water-wave problems mentioned above, and in the investigation of a new phenomenon occurring when a solitary wave pass over an abrupt deepening. In Section 7, a general discussion is presented and main conclusions are summarized. Finally, in a technical appendix, the proof of a convergence lemma is given.

Hamiltonian coupled-mode formalism for nonlinear water waves

In this section we briefly present a complete account of the new HCMT for fully nonlinear water waves over arbitrary bathymetry. To address the similarities and distinctions between the present theory and the well-developed, classical Hamiltonian theory by Zakharov [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] and Craig & Sulem [START_REF] Craig | Numerical Simulation of Gravity Waves[END_REF], we start by a quick review of the latter.

The classical Hamiltonian approach

The fully nonlinear water-wave problem in the two horizontal dimensions admits of a Hamiltonian formulation in terms of the free surface elevation η(x, t) and the trace of the wave potential on the free surface, ψ(x, t) = Φ(x, z = η(x, t), t), as canonical variables [START_REF] Craig | Numerical Simulation of Gravity Waves[END_REF], [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], [START_REF] Lannes | Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. In this context, the wave motion is governed by the two Hamiltonian evolution equations

∂ t η = G[η, h]ψ, (1a)
∂ t ψ = -gη - 1 2 |∇ x ψ| 2 + G[η, h]ψ + ∇ x ψ • ∇ x η 2 2 (1 + |∇ x η| 2 ) . ( 1b
)
where G [η, h] ψ denotes the Dirichlet-to-Neumann (DtN) operator defined by the formula

G[η, h]ψ = -∇ x η • [ ∇ x Φ] z = η + [ ∂ z Φ] z = η . (2)
In Eqs. ( 1)-( 2),

∇ x = (∂ x1 , ∂ x2 ) and ∇ = (∂ x1 , ∂ x2 , ∂ z )
denote the horizontal and three-dimensional (3D) gradients, respectively. Further, the wave potential Φ = Φ(x, z, t) satisfies the Laplace equation in the fluid domain together with the bottom impermeability condition on the seabed, and appropriate lateral conditions. Albeit a linear problem, the implementation of the non-local DtN operator, Eq. ( 2), which takes care of the substrate (interior and bottom) kinematics of the fluid, is the most crucial step for demanding numerical simulations of nonlinear water-waves, in terms of both accuracy and computational time. Craig & Sulem [START_REF] Craig | Numerical Simulation of Gravity Waves[END_REF] proposed a perturbative approach for computing the DtN operator in the case of a flat bottom, assuming periodic lateral conditions. This strategy, based on a functional Taylor-series expansion of the DtN operator around the state of zero free-surface elevation, η(x, t) = 0, has been further studied and advanced by many authors, see e.g. [START_REF] Nicholls | Traveling Water Waves: Spectral Continuation Methods with Parallel Implementation[END_REF], [START_REF] Craig | Water waves over a random bottom[END_REF], [START_REF] Xu | Numerical simulation of threedimensional nonlinear water waves[END_REF], leading to useful numerical schemes for solving various problems. This approach, being of perturbative character, is very efficient when slightly-to-moderately deformed fluid domains are considered. Extensions to large variations of the bathymetry profile become complicated, and of increased computational cost [START_REF] Craig | Hamiltonian long-wave expansions for water waves over a rough bottom[END_REF], [START_REF] Gouin | Development and validation of a non-linear spectral model for water waves over variable depth[END_REF], because of keeping more modes in the corresponding expansions. It should be noted here, that a perturbative method may fail as the number of modes increases. See, e.g. [START_REF] Nicholls | Stability of High-Order Perturbative Methods for the Computation of Dirichlet-Neumann Operators[END_REF].

A rapidly convergent "modal" series expansion of the wave potential

Recently, an alternative formulation has been proposed, based on an exact (nonperturbative), rapidly convergent series expansion of the unknown wave potential Φ(x, z, t) of the form

Φ(x, z, t) = ∞ n=-2 ϕ n (x, t)Z n (z ; η(x, t), h(x)) . (3)
As proved in detail in AP17, the vertical basis system {Z -2 , Z -1 , {Z n } n≥0 }, is composed of:

(i) The two specific functions

Z -2 (z; η, h) = µ 0 h 0 + 1 2h 0 (z + h) H 2 - µ 0 h 0 + 1 2h 0 H + 1, (4a)
Z -1 (z ; η, h) = µ 0 h 0 -1 2h 0 (z + h) H 2 + 1 h 0 (z + h) - µ 0 h 0 + 1 2h 0 H + 1 (4b)
where H = H(x, t) = η(x, t) + h(x) is the local depth of the fluid up to the instantaneous free surface, and µ 0 , h 0 are two auxiliary constants, which will be discussed in the sequel. These functions serve the purpose to free the expansion (3) of the boundary constraints imposed by the remaining part of the expansion (see Eqs. [START_REF] Bonneton | Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes[END_REF], below), resulting also in a significant acceleration of the convergence.

(ii) The set of local eigenfunctions {Z n = Z n (x, t)} n≥0 of a regular Sturm-Liouville problem, defined in the vertical interval [-h(x), η(x, t)]. In principle, any such problem can provide a theoretically acceptable system of local eigenfunctions. For physical reasons (see comments below), in constructing the HCMT, the choice has been made of the vertical Sturm-Liouville problem corresponding to linear water waves, with boundary conditions

∂ Z n ∂ z -µ 0 Z n = 0, at z = η(x, t), ∂ Z n ∂ z = 0, at z = -h(x). (5a, b)
Eigenfunctions {Z n } n≥0 , normalized to take the value 1 at z = η(x, t), are given by the equations:

Z 0 = cosh(k 0 (z + h) cosh(k 0 H) , Z n = cos(k n (z + h) cos(k n H) , n ≥ 1, (5c, d )
where k n = k n (x, t), n ≥ 0 are the roots of the following transcendental equations:

k 0 H tanh(k 0 H) = µ(x, t) (6a) k n H tan(k n H) = -µ(x, t), for n ≥ 1 (6b)
and µ(x, t) = µ 0 H(x, t). From the theory of Sturm-Liouville problems, it is known that the system {Z n } n≥0 is an L 2 -basis in each vertical interval (-h, η) [START_REF] Higgins | Completeness and basis propertries of sets of special functions[END_REF]. The addition of the two functions Z -2 , Z -1 makes the extended system {Z -2 , Z -1 , {Z n } n≥0 } a basis in the Sobolev space H 2 (-h, η), ensuring quick, point-wise convergence of the series itself to Φ(x, z, t), and of the term-wise differentiated series to the corresponding derivatives of Φ(x, z, t); see AP17 for a detailed proof. The additional basis functions, Z -2 , Z -1 , have been adopted by other researchers as well, who studied their role in the convergence of numerical schemes for solving the Helmholtz equation [START_REF] Hazard | An improved multimodal approach for non-uniform acoustic waveguides[END_REF], and ascertained their strong positive effects in solving problems in acoustical waveguides of irregular shape [START_REF] Mercier | Acoustic propagation in non-uniform waveguides: revisiting Webster equation using evanescent boundary modes[END_REF], [START_REF] Maurel | Improved multimodal admittance method in varying cross section waveguides[END_REF]. In our paper AP17, we have proved that, if η, h and Φ are sufficiently smooth functions, then the "modal" amplitudes ϕ n (x, t), and their derivatives ∇ x ϕ n (x, t), ∂ t ϕ n (x, t) decay as fast as O(n -4 ), while the decay of the first few modes (-2, -1, 0, 1, 2, 3) is even faster, namely exponential. This suggests that only a few modes are enough for accurate computations, as also confirmed by numerical experiments. The auxiliary constant h 0 , appearing in Eqs. (4), is introduced only for dimensional purposes, and its value is taken to be a characteristic depth of the studied configuration, e.g. the depth at the incident region, or the mean depth. As regards the auxiliary constant µ 0 , the following comments are in order. The essential role of this constant is to formulate the first boundary condition in [START_REF] Bonneton | Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes[END_REF] of the Sturm-Liouville problem defining the eigenfunctions {Z n } n≥0 . All theoretical statements made above remain valid for any value µ 0 > 0. Since, however, for a specific choice of µ 0 , the eigenfunctions {Z n } n≥0 become the physical modes of linear waves for angular frequency ω 0 = √ g µ 0 (at the local depth), it is preferable to select µ 0 in relation with a characteristic frequency of the problem, e.g. the frequency of the incident wave, or the central frequency of a wave packet or wave spectrum. In this way, the series expansion (3) encapsulates the physics of a "nearby" linear wave problem, even before using the nonlinear Hamiltonian equations (see Eqs. ( 7)), by means of which the solution will be determined taking fully into account all nonlinear features of the problem. For reasons explained above, this is not an approximation; it is, instead, an adaptation of the vertical expansion, making it converge even faster, since it takes a priori into account a part of the physical structure of the problem.

Comment on terminology. The two first terms of Eq. ( 3), ϕ -2 Z -2 and ϕ -1 Z -1 , are called, respectively, free-surface mode and sloping bottom mode, since they ensure the ability of the series to correctly satisfy the free-surface and bottom conditions. Also, they are referred to collectively as boundary modes. The term ϕ 0 Z 0 is referred to as the propagating mode, and the remaining terms, ϕ n Z n for n ≥ 1, are called evanescent modes, by analogy with the linear wave theory. It should be stressed, however, that the individual terms of the series (3) have not the usual meaning of modes (as in the linear wave propagation), since all terms together solve the nonlinear hydrodynamic problem. This is why we have used the word mode in quotes up to now, a convention which will be abandoned in the sequel, after these explanations.

The HCMT

Introducing the series expansion (3) of the wave potential into Luke's variational principle [START_REF] Luke | A variational principle for a fluid with a free surface[END_REF], and performing the variations with respect to the unknown fields η(x, t) and ϕ n (x, t), n ≥ -2, we obtain an infinite set of Euler-Lagrange equations. Elaborating further on this set of equations, we find the following two evolution equations with respect to η(x, t) and ψ(x, t) = n≥-2 ϕ n (x, t),

∂ t η = -(∇ x η) • (∇ x ψ) + ( |∇ x η | 2 + 1) h -1 0 ϕ -2 + µ 0 ψ , (7a)
∂ t ψ = -g η - 1 2 (∇ x ψ) 2 + 1 2 (|∇ x η | 2 + 1) h -1 0 ϕ -2 + µ 0 ψ 2 , (7b)
coupled, through the field ϕ -2 = ϕ -2 (x, t), with an infinite set of time-independent equations (with coefficients parametrically dependent on time) with respect to

ϕ n (x, t), n ≥ -2, given by ∞ n=-2 A mn ∇ 2 x + B mn • ∇ x + C m n ϕ n = 0, m ≥ -2, ( 8a
) ∞ n=-2 ϕ n = ψ. (8b)
The (x, t)-dependent matrix coefficients A mn = A mn (η, h), B mn = B 1 mn (η, h), B 2 mn (η, h) and C mn = C mn (η, h) are defined in terms of η(x, t) and h(x) through the vertical functions Z n = Z n (z; η(x, t), h(x)) , n ≥ -2, by the equations

A m n = η -h Z n Z m d z, (9a)
B m n = 2 η -h (∇ x Z n )Z m d z + (∇ x h) [Z m Z n ] z=-h , (9b)
C m n = η -h (∇ 2 x Z n + ∂ 2 z Z n ) Z m d z -N h • [(∇ x Z n , ∂ z Z n )Z m ] z=-h . (9c)
Equations ( 8) describe the kinematics of the fluid (substrate kinematics) at each time, for given bathymetry h(x), free-surface elevation η(x, t) and free-surface potential ψ(x, t). They form an elliptic problem which should be supplemented by appropriate lateral boundary conditions, dependent on the specific problem considered. A detailed description of the lateral boundary conditions for various specific problems, both in 2D and in 3D configurations, can be found in the Supplementary Material (Appendix C) of the paper [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with non-planar boundaries[END_REF], and in [START_REF] Papoutsellis | Implementation of a fully nonlinear Hamiltonian Coupled-Mode Theory, and application to solitary wave problems over bathymetry[END_REF], available online through the link https://arxiv.org/abs/1710.10847 . In the above formulation, Eqs. ( 7) bear in mind the classical Hamiltonian formulation, Eqs. ( 1) -( 2

), with ϕ -2 (x, t) being a nonlocal coefficient (different, yet) related with the DtN operator G [η, h] ψ through the equation G [η, h] ψ = -(∇ x η) • (∇ x ψ) + ( |∇ x η | 2 + 1) h -1 0 ϕ -2 + µ 0 ψ . (10)
Eq. ( 10) is obtained by comparing Eqs. ( 1) and ( 2) with Eq. ( 7); as has also been established independently in AP17. Further, as proved in the same paper, the approximation of ϕ -2 (x, t), and thus of the DtN operator through Eq. ( 10), resulting from solving the (truncated version of the) modal kinematical problem, Eq. ( 8), turns out to be amazingly effective for any (smooth) bathymetry. This is due to the fact that the truncated problem (8) exhibits a miraculous superconvergence with respect to the number of modes N tot kept in the truncated expansion; the error diminishes at a rate proportional to O(N -6.5 tot ), for any shape (however steep) of the free-surface and bottom boundaries. Accordingly, a small number of modes suffices to provide very accurate results for the modal amplitude ϕ -2 (x, t), making Eqs. [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF] essentially exact. This is also confirmed by the examples provided herein, in Section 6, as well as by other simulations concerning demanding problems with strong nonlinearity and dispersion; see e.g. [START_REF] Athanassoulis | New form of the Hamiltonian equations for the nonlinear water-wave problem, based on a new representation of the DtN operator, and some applications[END_REF], [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with nonplanar boundaries[END_REF], [START_REF] Papoutsellis | Implementation of a fully nonlinear Hamiltonian Coupled-Mode Theory, and application to solitary wave problems over bathymetry[END_REF].

3 Implementation of the HCMS

Scheme of numerical solution

In order to solve the evolution Eqs. [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF], we need to calculate the nonlocal coefficient ϕ -2 (x, t) at each time step. That is, we have to solve the elliptic problem (8) for given η(x, t) and ψ(x, t) (known from the previous iteration, or from the initial conditions). For this purpose, system ( 8) is truncated at a finite number of modes N tot = M + 3, where M is the order of the last kept evanescent mode, and discretised by using central finite differences of fourthorder accuracy. For example, in the case of one horizontal dimension (∂ x2 ≡ 0) and periodic lateral conditions, the discretised version of Eqs. (8) reads as follows:

M n=-2 - A i m n 12∆x 2 + B i m n 12∆x ϕ i-2 n + 4A i m n 3∆x 2 - 2B i m n 3∆x ϕ i-1 n - 5A i m n 2∆x 2 -C i m n ϕ i n + 4A i m n 3∆x 2 + 2B i m n 3∆x ϕ i+1 n - A i m n 12∆x 2 + B i m n 12∆x ϕ i+2 n = 0, i=1,...,N X , m=-2,...,Ntot-1, ( 11a
) M n=-2 ϕ i n = ψ i , i = 1, ..., N X (11b)
where ∆x is the mesh size, N X is the number of grid points, and the upper index ( • i ) denotes discrete local values; e.g.

ϕ i n = ϕ n (x i , t) and A i m n = A m n (η(x i ), h(x i )).
The convergence and accuracy of the above scheme is investigated in AP17; see also [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF].

Having solved the above linear system, the local values of the free surface mode ϕ i -2 = ϕ -2 (x i , t), i = 1, ..., N X , are extracted and used in the evolution Eqs. (7a,b). Introducing the vector field U = (η, ψ) T , system (7) is written as ∂ t U = N (t, U ) whereN (t, U ) is the vector function defined by the right-hand sides of Eqs. (7a,b). Using now a temporal discretization t n , n = 1, ... , N T , and the notation U n = U (x, t n ) = (η(x, t n ), ψ(x, t n )) T , system [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF] is marched in time by a straightforward adaptation of the classical Runge-Kutta method, based on the Butcher tableau RK41 given in [59, p. 91]. The whole scheme of numerical solution of the HCMS, Eqs. ( 7) and [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], is presented in Algorithm 1.

It should be noted that Algorithm 1 requires four evaluations of the matrix coefficients

A = (A m n ), B = B (1) m n , B (2) m n
and C = (C m n ) at each point of the spatial grid, for the implementation of one time step. Thus, a huge number of evaluations of these coefficients is required, which may count up to 10 11 times for a long-time simulation. Accordingly, the fast and accurate calculation of all elements of these coefficients is of fundamental importance for the efficient implementation of the HCMS. Modifications of the above numerical scheme for solving other variants of water-wave problems, e.g. problems with different boundary conditions (vertical impermeable walls or generating and absorbing layers) are possible, and have been also developed. More details can be found in [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF]Chapter 7] and [START_REF] Hunt | Direct solution of wave dispersion equation[END_REF].

Algorithm 1 Numerical solution of the HCMS, Eqs. ( 7), [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF] Given U 0 = (η 0 , ψ 0 ), Calculate k n (η 0 , h), Eqs. ( 6), and A(η 0 , h), B(η 0 , h), C(η 0 , h), Eqs. [START_REF] Memos | On Higher-Order Boussinesq-Type Wave Models[END_REF].

Solve CMS, Eqs. [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], with ψ = ψ 0 , to find ϕ -2 (x, t)

for m = 1 → N T do for l = 1 → 4 do Ū l ≡ (η l , ψl ) T = U m + δt l-1 j=1 a l j K m j
Calculate k n (η l , h), Eqs. [START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF], and A(η l , h), B(η l , h), C(η l , h), Eqs. [START_REF] Memos | On Higher-Order Boussinesq-Type Wave Models[END_REF] and Solve CMS, Eqs. [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], with ψ = ψ l , to find ϕ -2 (x, t)

K m l = N (t m + c l δt, Ū l )
end for

U m+1 ≡ ( η m + 1 , ψ m + 1 ) = U m + δt 4 j = 1 b j K m j end for
3.2 Analytic calculation of the (x, t)-varying matrix coefficients A,B and C

The best approach to ensure fast and accurate evaluations of all elements of the matrix coefficients A, B, C, is to derive closed-form analytic expressions of them. This requires extensive analytic manipulations, which have been performed and tested against numerical integration using Eqs. [START_REF] Memos | On Higher-Order Boussinesq-Type Wave Models[END_REF]. A systematic account of these calculations can be found in [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF]. The main findings are: (i) all coefficients can be expressed in closed form, in terms of h 0 , µ 0 , η, h, k n (η, h), and (ii) using these closed forms, accurate evaluations of the coefficients becomes about 150 times faster than by using numerical integration, at the same level of accuracy.

As an example, we give here the closed form expressions of the coefficients A mn and B mn for the case m , n ≥ 1:

A mn = H 2 + Hµ 2 0 -µ0 2k 2 n , m = n ≥ 1, 0, m = n, (12a)
B m n = 2 F (1) n I (1) mn + F (2) n I (2) mn + F (3) n I (3) mn + (∇ x h) [Z n Z m ] z=-h , (12b)
with

F (1) n = -µ 0 1 + H ∂ H k n k n ∇ x H, I (1) mn = A mn , F (2) n = -k n ∇ x h, I (2) mn =    µ 2 0 2k 3 m , m = n ≥ 1, - k -1 n µ 2 0 +kn-kn[Z n Zm] -h k 2 n -k 2 m , m = n, F (3) n = -∂ H k n ∇ x H, I (3) n =    -µ0+H(µ 2 0 -k 2 m ) 4k 3 m , m = n ≥ 1, - k -1 n µ0(Hµ 0 -1)+knH k 2 n -k 2 m , m = n,
where the derivative of k n with respect to H = η + h is obtained by applying the implicit function theorem to the defining relation, Eq. (6b), and reads as follows:

∂ H k n = k n ( k 2 n + µ 2 0 ) µ 0 -H(k 2 n + µ 2 0 ) , n ≥ 1. ( 13
)

Semi-analytic calculation of the roots k n (x, t) of the local dispersion relation

Having established explicit formulae for all coefficients A m n , B m n and C m n in terms of h 0 , µ 0 , η,h and k n (η, h) = k n (x, t), with h 0 ,µ 0 , η, h already known (η is known from the previous time step), the question of their fast, accurate and robust evaluation is reduced to the corresponding question for the roots k n (x, t) of the local dispersion relations, Eqs. [START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF]. Solving these equations by means of the Newton-Raphson method, or any other iterative procedure, is in principle a trivial task. However, since this evaluation is to be performed for 10 8 -10 11 times, it is imperative the convergence of the solution scheme to be a priori ensured, and the number of iterations to remain small. For that reasons, it is important to establish highly accurate initial values for all k n 's or, even better, to derive satisfactory closed-form approximations.

For the case of k 0 , explicit approximations, accurate up to the third or fourth decimal digit, for all values of the parameter µ = µ(x, t) = µ 0 H(x, t) are available in the literature; see e.g. [START_REF] Beji | Improved explicit approximation of linear dispersion relationship for gravity waves[END_REF], [START_REF] Simarro | Improved explicit approximation of linear dispersion relationship for gravity waves: Another discussion[END_REF], [START_REF] Vatankhah | Improved explicit approximation of linear dispersion relationship for gravity waves: a discussion[END_REF], [START_REF] Massel | Extended refraction-diffraction equation for surface waves[END_REF], [START_REF] Fenton | On calculating the lengths of water waves[END_REF], [START_REF] Householder | The numerical treatment of a single nonlinear equation[END_REF]. In this work, the approximation of [START_REF] Massel | Extended refraction-diffraction equation for surface waves[END_REF] is used, which is based on Newton-Raphson iterations of the form

j+1 κ 0 (µ) = j κ 2 0 (µ) + µ cosh 2 j κ 0 (µ) j κ 0 (µ) + 0.5 sinh (2 j κ 0 (µ)) , (14)
with initial guess

0 κ 0 (µ) = µ + µ 1.986 exp(-1.863 + 1.198µ 1.366 ) √ tanh µ . ( 15
)
The initiation formula [START_REF] Papoutsellis | A new efficient Hamiltonian approach to the nonlinear water-wave problem over arbitrary bathymetry[END_REF] has maximum relative error of approximately 10 -4 , and the methodology leads to relative errors of less than 10 -15 for any value of µ, within 2 iterations.

For the case of k n , n ≥ 1, a single, simple and sufficiently accurate approximation, valid for all values of µ = µ 0 H(x, t), is apparently not available. In addition, the application of the Newton-Raphson method to the solution of Eq. (6b) is more involved, especially for large values of µ and small values of n, where, the steep gradient of tan function may lead to a large value of iterations, or to non-convergence or, even worse, to convergence to an erroneous solution. Sections 4 and 5 of this paper are devoted to the derivation, analysis and application of new, highly accurate approximations for the roots of Eq. (6b), that optimise the computation of k n , n ≥ 1, for all values of µ and n, and, when used as initialization of an iterative solver, lead to machine precision tolerance 10 -15 with no more than 3 iterations.

Root Approximation Strategy

Eq. (6b) can be written in the form

g(κ ; µ) ≡ κ tan(κ) + µ = 0, ( 16
)
where κ n = k n (x, t)H(x, t) and µ = µ(x, t) = µ 0 H(x, t). The first three roots of Eq. ( 16) are depicted in Fig. 1, for three different values of the parameter µ = 0.1, 1.0, 5.0. The difficulties in obtaining efficient root-finding formulae for Eq. ( 16) are well demonstrated in this figure. First, the roots κ n are dependent on the parameter µ which, in the context of HCMT, is a variable quantity, µ = µ 0 H(x, t) spanning a wide range of values, since the method is applied to all depths and wave amplitudes. Recall also that µ 0 is an auxiliary numerical parameter, arbitrarily chosen, although some general principles for its appropriate (but nonunique) selection have been discussed in Section 2.2.

Second, the few first roots, which are the most important for applications, are associated to large values and steep gradients of the tan function, especially for large values of µ. Accordingly, the obvious initial choice 0 κ n = nπ, which is the asymptotic limit of κ n for n → ∞, is not appropriate in these cases, and may lead to a large number of Newton-Raphson iterations, or to nonconvergence or, even, to erroneous solutions.

Beginning with the observation that for any n ∈ N the roots of Eq. ( 16) satisfy

κ n ∈ (2n -1)π 2 , (2n + 1)π 2
and lim κ n → nπ, for all µ > 0, [START_REF] Craig | Numerical Simulation of Gravity Waves[END_REF] and following Newman [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF], we introduce the quantity

ε n (µ) = nπ -κ n (µ). (18)
Clearly, ε n (µ) satisfy the following relations:

ε n (µ) ∈ (0, π/2) , lim n→∞ ε n (µ) → 0, lim µ→∞ ε n (µ) → 0. ( 19
)
The values ε n (µ) are monotonically decreasing with respect to n for each µ > 0, and monotonically increasing with respect to µ for each n ∈ N. The relation between κ n (µ) and ε n (µ) is depicted in Fig. 2. Using Eq. ( 18), Eq. ( 16) is written as

(nπ -ε) tan(nπ -ε) = -µ, which, upon expanding tan(nπ -ε), becomes (nπ -ε) tan(ε) = µ. (20)
By rewriting Eq. ( 20) in the form tan(ε) = µ/κ, a Picard iteration is readily defined as follows

j+1 ε n = Arctan µ j κ n
, for j = 0, 1, 2, ... and an initial estimation 0 ε n , (21) where j κ n = nπ -j ε n . For later use we also write down the formula

cot j+1 ε n = j κ n /µ, (22)
which is, of course, equivalent with Eq. ( 21). On the other hand, simple algebraic manipulations show that the function g(κ; µ), Eq. ( 16), can be written in the form g(κ = nπ -ε ; µ) = -f (ε ; µ)/ tan(ε), where

f (ε ; µ) = nπ -ε -µ cot(ε). ( 23
)
Since tan(ε n ) = 0, ± ∞, the solutions ε n of Eq. ( 21) can also be found as the roots of the function f (ε ; µ). A Newton-Raphson procedure for the determination of ε n such that f (ε n ; µ) = 0, is given by the formula

j+1 ε n = j ε n - f ( j ε n ; µ) ∂ ε f ( j ε n ; µ)
, for j = 0, 1, 2, ... and an initial estimation 0 ε n , which specializes in the form

j+1 ε n = j ε n - nπ -j ε n -µ cot( j ε n ) µ -1 + µco t 2 ( j ε n ) , (25)
upon substituting f = f (ε ; µ) from Eq. ( 23) to Eq. ( 24). Further, a third-order formula can be obtained by employing the Householder method [63, Section 4.4], having the form

j+1 ε n = j ε n - f ( j ε n ; µ) ∂ ε f ( j ε n ; µ) 1 + f ( j ε n ; µ)∂ εε f ( j ε n ; µ) 2 (∂ ε f ( j ε n ; µ)) 2 . ( 26
)
The above iteration equation takes the specific form

j+1 ε n = j ε n - nπ -j ε n -µ cot( j ε n ) µ -1 + µ cot 2 ( j ε n ) × × 1 -2µ nπ -j ε n -µ cot( j ε n ) cot( j ε n ) + cot 3 ( j ε n ) 2 µ -1 + µ cot 2 ( j ε n ) 2 , (27)
upon substituting the function f from Eq. ( 23) into Eq. ( 26). In principle, using the general Householder method, formulae of order higher than 3 can be produced. The root approximation strategy proposed herein is based on successively applying Eq. ( 21) and Eq. ( 25) or Eq. ( 27). These two variants lead to a second and third order method respectively. The exact steps, involved in the general approximation algorithm, are summarised in Algorithm 2.

Algorithm 2 Root Approximation Strategy

Given 0 κ n , for j = 0, 1, 2, ...

(1) Calculate j+1 ε n = Arctan µ/ j κ n

(2) Calculate j+2 ε n by applying Eq. ( 24) or ( 26), replacing j by j + 1

(3) Set j+2 κ n = nπ -j+2 ε n (4) Replace j +2 by j +1 (in the left-hand side only) to derive the combined, two-step, iteration formula Applying now the above described procedure to the iteration formulae ( 25) and ( 27), we obtain the following two iteration schemes:

• 2 nd order method j+1 κ n = nπ + µ nπ -j κ n µ 2 + j κ 2 n -µ - µ 2 + j κ 2 n µ 2 + j κ 2 n -µ Arctan µ j κ n . (28)
• 3 rd order method

j+1 κ n = nπ + µ(nπ-j κn) µ 2 + j κ 2 n -µ - µ 2 + j κ 2 n µ 2 + j κ 2 n -µ Arctan µ j κn - j κnµ(µ 2 + j κ 2 n ) (µ 2 + j κ 2 n -µ) 3 nπ -j κ n -Arctan µ j κn 2 . ( 29
)
The following remarks are in order at this point:

Remark 4.1 Iteration formulae (25), ( 27) and ( 28), ( 29) are relatively simple. Each one of them involves only a few algebraic operations and the calculation of one transcendental function.

Remark 4.2 Setting 0 κ n = nπ in Eqs. ( 28), ( 29), explicit forms of reasonable accuracy are produced, which can in turn be used for the initiation of successive approximations through [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF] and (29). This direction is exploited further in the next section.

j+1 ε n -ε n ≤ r n (µ) j ε n -ε n 2 + O j ε n -ε n 3 , ( 30
) j+1 ε n -ε n ≤ R n (µ) j ε n -ε n 3 + O j ε n -ε n 4 , (31)
respectively. Then, the compound iterations, defined by Eqs. ( 28) and (29), are also of second and third order, respectively, but now satisfy the error estimates

j+1 κ n -κ n ≤ δ 2 n (µ)r n (µ) j κ n -κ n 2 + O j κ n -κ n 3 , ( 32
)
j+1 κ n -κ n ≤ δ 3 n (µ)R n (µ) j κ n -κ n 3 + O j κ n -κ n 4 , (33)
respectively, where δ n (µ) = 4µ/(2n -1) 2 π 2 + 4µ 2 < 0.32, for all µ and n. That is, the positive constants controlling the convergence rate are now considerably smaller, since they are multiplied by the small factors δ 2 n (µ) and δ 3 n (µ), respectively. The proof of the above proposition is given in Appendix A. The factor δ n (µ) takes its maximum value at µ = (n -1/2)π, and satisfies the inequality

δ n (µ) = 4µ (2n -1) 2 π 2 + 4µ 2 ≤ 1 (2n -1)π ≤ 1 π < 0.32.
Clearly, it decreases with increasing n, and its value tend to zero (for all n ∈ N) in both limiting cases µ → 0 and µ → ∞. A plot of the function δ n (µ), for n = 1, 2, 3, 4, is shown in Fig. 3.

Initiation, performance, and closed-form approximations

Improved initiation formulae for the iterations ( 28) and ( 29) can be derived by means of Eq. ( 21) and the Eqs. ( 28) and ( 29) themselves. Substituting 0 κ n = nπ in Eqs. ( 21) and [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF], and denoting the obtained we obtain 1 κ n by A n (µ) and B n (µ), respectively, we obtain the formulae

A n (µ) = nπ -Arctan µ nπ , (34)
B n (µ) = nπ - µ 2 + n 2 π 2 µ 2 + n 2 π 2 -µ Arctan µ nπ . ( 35
)
Figure 4: Relative % error of the computation of the first three roots of Eq. (6b) using A n (µ), Eq. ( 34), B n (µ), Eq. ( 35), and C n (µ), Eq. ( 36), plotted as a function of the parameter µ to machine precision after the second iteration. Especially, Eq. ( 29) initiated by E n (µ), leads to the determination of κ n (µ) with an error less than 10 -15 for the whole range of µ. Both methods appear to yield exact results for a particular value of µ, that depends on the number of the mode under consideration. This value of µ tends to be higher, as the number of mode increases.

The main advantage of the present schemes in comparison with those proposed by Newman [START_REF] Newman | Numerical solutions of the water-wave dispersion relation[END_REF], is the presence of powers of the small factor δ n (µ) in the error constant (see Proposition 1 and Remark 3 above). It is remarkable that the regions of increased errors in Fig. 5 and 6 coincide, more or less, with the region of high values δ n (µ), as shown in Figure 3. In any case, these errors are already small and disappear at the second iteration. The first six roots κ n , as obtained by two iterations of Eq. ( 29), initiated by E n (µ), are plotted in Figure 7, as functions of µ.

Before proceeding with the application of the above results to the simulation of demanding, nonlinear, water-wave problems, we shall discuss the enhancement of the Newton-Raphson (NR) method by the exploitation of Eq. ( 35) -( 37) for its initialization. For µ ∈ (0, 85), we compare the number of iterations required by the NR method in order to reach machine precision tolerance by starting from the following four different initial guesses: 0 κ n = nπ, B n (µ), C n (µ), D n (µ). Results for the first three roots, κ 1 , κ 2 , κ 3 , are shown in Figure 8. As one might expect, the choice 0 κ n = nπ is sufficient only when µ is relatively small. The number of iterations increases rapidly with increasing µ, leading eventually to overshooting or divergence of the NR method after some value µ * , which increases with the root index n. The initial guess 0 κ n = C n 28) initiated by D n (µ) for the first 6 roots of Eq. ( 16) leads to convergent NR iterations, but the number of iterations increases significantly for largeµ. On the other hand, the choices 0 κ n = B n , D n are the more robust, leading to convergence for the whole range of µ after 4-5 NR iterations.

Applications

In this section, we present a number of applications of the HCMT to specific, nonlinear, water-wave problems, and discuss in detail the role of the new approximation formulae for the roots κ n to the implementation of the method. Before proceeding to the specific applications, it is expedient to recall that the HCMS Eqs. [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF], closed by the substrate kinematical problem [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], is just an exact reformulation of the fully nonlinear water-wave problem, Eqs. ( 1), [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF], provided that the coefficients A mn , B mn and C mn are exactly calculated. Given that A mn , B mn and C mn have been analytically calculated in terms of κ n , the issue of the accuracy of the new Hamiltonian formulation is fully controlled by the accuracy of the determination of κ n . Having ensured the latter, by using the results of Sections 4 and 5, the new HCMS becomes very efficient because of the striking superconvergence of ϕ -2 with respect to the number of modes, N tot , retained in the truncated version of system (8), as already mentioned. Thus, only a few modes (up to five or six) are enough for accurate hydrodynamic calculations.

In problems studied in this section, the initialization of the evolutionary Hamiltonian equations ( 7) is performed by means of an accurately pre-calculated Figure 6: Logarithmic plot of the relative error of the approximation E n (µ) and the two first iterations of Eq. ( 29) initiated by E n (µ) for the first 6 roots of Eq. ( 16) Figure 7: Plot of the first six roots of Eq. ( 16), as approximated by Eq. ( 29), initiated by E n (µ), with two iterations Figure 8: Number of iterations of the Newton-Raphson method for the first three roots of Eq. ( 16), for initial guesses 0 κ n = nπ, B n (µ), C n (µ), D n (µ) steady travelling wave (in constant depth), of specific wavelength and height. Such a wave provides the initial conditions in the case studied in Subsection 6.1, while it plays the role of the excitation (incoming wave) in the cases studied in Sections 6.2 and 6.3. High-accuracy numerical simulations of steady travelling waves can be obtained by various methods/codes provided, for example, by [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF], [START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF], [START_REF] Clamond | Accurate fast computation of steady twodimensional surface gravity waves in arbitrary depth[END_REF] or an appropriate version of the present HCMT, given in [START_REF] Papoutsellis | Nonlinear water waves over varying bathymetry: Theoretical and numerical study using variational methods[END_REF]Chapter 6]. The latter method was used in all cases presented below. Note that all methods give practically identical results.

Evolution of steady travelling waves over flat bottom

The problems considered in this subsection are simple and permit an easy monitoring of the accuracy of the numerical solution. HCMS, Eqs. [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF], is solved under periodic lateral conditions, initialized by the fields (η(x, t 0 ), ψ(x, t 0 )) = (η 0 (x), ψ 0 (x)) corresponding to a travelling nonlinear waver with wave length equal to the horizontal extent of the domain. Thus, the evolution equations have to reproduce the initial state after one period, two periods etc. permitting a trivial accuracy check. Three different waves are considered, corresponding to deep, intermediate and shallow water conditions. Normalizing the depth to h 0 = 1, the wave lengths of the three cases are chosen to be λ = 1 m (deep water), λ = 5 m (intermediate-depth water), and λ = 18 m (shallow water). In all cases, the wave height has been selected to be the 80% of the maximum value for the specific choice of λ/h, as predicted by Williams [START_REF] Williams | Limiting gravity waves in water of finite depth[END_REF]. The parameter µ 0 is chosen to be µ 0 = (2π/λ) tanh(2πh/λ), that is, the value of µ 0 corresponding to the linear waves with the same wavelength. This is the most appropriate choice in this case, since it ensures that the representation of the flow field by the series (3) already encapsulates the physics of the corresponding linear wave, even before starting the procedure of the numerical solution. Of course, to obtain an accurate approximation, Eq. ( 16) is solved for the varying µ(x) = µ 0 (η(x)+ h), so that the local vertical basis, at each (x, t), to accurately represent the local wave potential field.

The free-surface wave profiles for the three cases studied are shown in the left panels of Figures 91011, along with the corresponding local values of µ = µ(x) = µ 0 (η(x) + h). The evolution problems are solved for a time span of three periods, by using N tot = 3, 4 , 5, 6 modes, that is, the modes -2, -1, supplemented by 1, 2, 3, 4 evanescent modes. The coefficients A mn , B mn , C mn are calculated by means of closed-form equations, like Eqs. [START_REF] Brink | Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves[END_REF], in terms of the roots κ n of Eq. ( 16). The whole numerical scheme is implemented as explained in Subection. 3.1.

For comparison purposes, Eq. ( 16) is solved by four different methods:

1. The standard Newton -Raphson (StNR) method applied to Eq. ( 16), initialized by 0 κ n = nπ, with tolerance 10 -15 (with as many iterations J iter (n) as needed),

2. The NR method applied to Eq. ( 16), now initialized by B n (µ), Eq. ( 35), with the same tolerance 10 -15 , called subsequently the improved NR (ImNR) method, 3. The 2 nd order semi-explicit (2SE) method, Eq. ( 28), initiated by D n (µ), Eq. ( 38), with J iter (1) = 3 and J iter (n) = 2, for n ≥ 2, 4. The 3 r d order semi-explicit (3SE) method, Eq. ( 29), initiated by E n (µ), Eq. ( 39), with J iter (1) = J iter (2) = 2, and J iter (n) = 1 for n ≥ 3.

A first important numerical finding from the applications presented below, is that the StNR method does not converge for κ 1 in the deep water case (λ/h = 1), since the values of µ = µ(x) become large (of the order of 10). The ImNR method, however, does converge in all studied cases. The two semi-explicit methods work with pre-specified number of iterations, being thus free of the convergence question.

In order to illustrate the nonlinear accuracy and efficiency of the present method we investigate the relative L 2 -error of the free surface elevation, computed by (cf. [START_REF] Bingham | On the accuracy of finite difference solutions for nonlinear water waves[END_REF])

Error = ||η 3 -η 0 || 2 3||η 0 || 2
where η 0 is the given initial wave profile, and η 3 is the computed wave profile after three periods. The horizontal grid consists of N X = 128 points in the cases λ/h = 1, 5, and N X = 256 points in the case λ/h = 18. The Courant-Friedrichs-Lewy (CFL) number, based on the group velocity of the nonlinear wave, is 0.7 for all cases. The roots κ n of Eq. ( 16) are computed by using the 2SE method. Results concerning the error are shown in Table 1. The convergence of the numerical scheme with respect to the number of modes N tot is clearly demonstrated in the same table. It is clearly seen that a total number of 5-6 modes is sufficient for an error of order 10 -4 -10 -5 , which is comparable with the corresponding results of Bingham & Zhang [10, Section 5], obtained by using 4 th order finite differences in the whole domain. The errors typically increase for larger simulation times and/or higher nonlinearity, which is expected in the context of the present explicit time-integration scheme. With the same spatio-temporal discretisation, a total number of modes N tot = 7 was needed for a stable simulation of duration 50T in the deep water case producing an error of 1.9 × 10 -3 . For the same duration, the errors in the intermediate and shallow water cases (with N tot = 6) are 1.0 × 10 -4 and 1.2 × 10 -4 . For higher nonlinearity (wave height 90% of the maximum value) the errors after 50T for the deep, intermediate and shallow case are 8.0×10 -3 , 8.8×10 -3 and 3.1×10 -2 respectively, by using N tot = 7. The above computations have been performed without the use of smoothing or filtering. The complete investigation of the numerical stability and limitations of the present scheme in conjunction with steady periodic waves will be the subject of another work. In the case of solitary waves, results can be found in [START_REF] Papoutsellis | Implementation of a fully nonlinear Hamiltonian Coupled-Mode Theory, and application to solitary wave problems over bathymetry[END_REF].

The relative computational time needed for the evaluation of κ n , n = 1, ... , 4 (N tot = 6) during the simulation of three periods, is documented for increasing spatial resolution N X = 32 , 64, 128, 256. Results for the above three cases are shown in the right panels of Figures 91011. The ImNR method and the 3SE method are the most expensive in all three cases. This is due to the relatively large number of iterations as concerns the ImNR, and to the more involved mathematical expressions as regards to 3SE. Thus, 2SE method turns out to be the best choice, being always the faster one, predominantly in deep water conditions. It should be stressed once again that the standard NR method, initialized by obvious choice 0 κ n = nπ, fails to converge in many cases, especially for deep water conditions. The improved NR method (that is, the NR method initiated by B n (µ) or D n (µ)) is proven to be convergent for all values of n and µ(x). These findings underline the significance of the proposed root approximation strategy and the explicit formulae given in Section 5. (*) In this case, the number of modes (N tot = 3) is not sufficient for the convergence of the algorithm.

Transformation of incident waves over a submerged bar

In order to further investigate the effect of semi-explicit methods derived in this paper, we consider the transformation of incident regular waves due to a submerged trapezoidal bar. This configuration has been investigated in the experiments of Beji & Battjes [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF] and Dingemans [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF], and is considered a standard benchmark test for the ability of numerical models to correctly simulate nonlinear and dispersive waves over variable bathymetry. A detailed investigation of the performance of HCMS in this experiment is presented in [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with nonplanar boundaries[END_REF]. The results obtained here correspond to N tot = 7 modes (the modes -2, -1, the propagating one, and four evanescent ones). The parameter µ 0 is chosen as µ 0 = ω 2 0 /g, where ω 0 is the circular frequency of the incident wave, ω 0 = 2π/T , with T = 2.02 sec. The spatio-temporal discretization is N X = 707 and N T = 3500. We consider the three different methods, ImNR, 2SE and 3SE, for calculating the four roots k n , n = 1, 2, 3, 4, at any horizontal position x and time t. In all three simulations, the matrix coefficients are evaluated at machine precision. The fluid domain at a specific simulation time is shown in Figure 12, together with the correspond- ing values of µ(x). The validity of the simulations is illustrated in Figure 13 by comparison of the time series of the free surface elevation at four indicative measuring stations (St. [START_REF] Karambas | Boussinesq Model for Weakly Nonlinear Fully Dispersive Water Waves[END_REF][START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF][START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF][START_REF] Gagarina | Variational space-time (dis) continuous Galerkin method for nonlinear free surface water waves[END_REF] used in [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF]. As expected, the computational results are identical and in very good agreement with experimental measurements. Differences are present only in the total simulation time. ImNR and 3SE methods result in practically the same computational time, while the 2SE method leads to 2% faster simulation. The small time enhancement observed is due to the fact that the time needed for the solution of the transcendental equation ( 16) is small in comparison with the time needed for the numerical solution of the HCMS, Eqs. ( 7), [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF]. Most important, however, is the robustness of the method.

Bragg scattering over sinusoidal bottom ripple patch

In this example, we study the transformation of an incident regular periodic wave propagating in a region having a sinusoidal patch in the seabed. We consider the configuration studied by Davies & Heathershaw [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], where this wave-bottom interaction was first investigated. The bottom sinusoidal patch, of wavelength l b , extends from x l to x r = x l + 10l b , and is defined by

z = -h m + d sin(k b x), x l < x < x r , -h m , elsewhere,
where k b = 2π/l b . The incident wave is chosen to have a wavenumber k = k b /2 (Bragg resonance condition), steepness kH/2 = 0.05 and period T = 1.28 sec. We consider HCMS with N tot = 6, N X = 1000, and µ 0 = k tanh(kh m ); see Figure 14. Simulations correspond to a duration of 40T , which is sufficient in order to estimate the reflection coefficient from the analysis of the local time series of the free surface elevation [START_REF] Goda | Estimation of incident and Reflected waves in random wave experiments[END_REF]. A comparison of experimental data and our computations, obtained by using the three methods for the computation of k n , n = 1, 2, 3, are shown in Figure 15. Results corresponding to the Higher Order Spectral (HOS) method [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF] have been digitized and are also shown in the same figure. Our three variants of computations are indistinguishable, as expected, and show excellent agreement with the measurements of [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF] over the rippled patch. HOS method slightly underestimates the reflection coefficient over this region. On the upwave side, all numerical computations are in agreement, differing from the measurements of [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF]. This is due to unwanted interference effects occurring in the experimental tank, as reported in [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF]. As in the previous case, the total computational time is approximately the same when using ImNR and 3SE methods, being 3% larger than the one corresponding to using 2SE method. Again, the small time enhancement is due to the fact that the time needed for the solution of Eq. ( 16) is small in comparison with the total time needed for the numerical solution of Eqs. ( 7), [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF].

Disintegration of a solitary wave due to an abrupt deepening

As a final example, we consider the transformation of a solitary wave propagating over an abrupt deepening, that is, a strong depth increase. The interesting feature of this numerical experiment is the transition from shallow to not-shallow water conditions which induces multiple space and time scales, possibly out of the reach of simplified models. To the best of our knowledge, this nonlinear phenomenon has not been previously examined. In the numerical simulation of the above described problem, the horizontal domain extents from x = 0 to x = 2000m. The bottom surface is defined by the equation

h(x) = -a -b tanh(0.2(x -625)),
where a > b > 0 are shape parameters defining the left and right (small and large, respectively, see Figure 16) depth of the seabed:

h l = a -b, h r = a + b,
This bathymetry represents a smooth but abrupt depth transition from h l = h 0 = 1m to a larger depth h r , which has been taken to be h r = 2h 0 , 4h 0 , 8h 0 . The initial solitary wave is centered at x = 550 m having an amplitude a = 0.3m (see Figure 16 (a)). Initial conditions (η(x, t 0 ), ψ(x, t 0 )) = (η 0 (x), ψ 0 (x)) are obtained by using the highly accurate solitary wave solutions of the complete nonlinear water wave problem, provided by [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. HCMS is implemented by using N tot = 8 modes and a spatio-temporal discretization δx = 0.2 m and δt = 0.03 s. The parameter µ 0 is chosen as µ 0 = k tanh(kh 0 ) with k = 2π/L, L being the support of the free surface elevation as computed by the code of [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF].

In the first few seconds of the simulation, as the solitary wave is moving over the shallow-flat part of the bottom, it travels steadily towards the right. The moment it encounters the depending, a violent change of shape takes place. The crest height of the solitary wave drops suddenly, and a small wave of depression emerges, back-propagating towards the left. As the leading solitary-like wave travels towards the right, moving over the deep part of the seabed (h r = 4m), its crest height steadily reduces and a highly oscillating wave trail, of smaller amplitude, is being developed and expanded horizontally. This trailing wave follows the leading wave, giving rise to a dispersive wave pattern that propagates over intermediate-to-deep water conditions. Besides, as the reflected wave propagates towards the left, also develops a small yet visible dispersive trail. Some snapshots of the above described wave-bottom interaction pattern are shown in Figure 17. The full simulation can be found in video 1. In order to study the influence of the intensity of the deepening on the phenomenon, two further simulations are performed, with h r = 2h 0 and h r = 8h 0 , keeping all other parameters the same. The time history of the maximum surface elevation along the evolution is shown in Figure 18. In all cases, the initial sudden decrease of the amplitude is followed by a small and brief bounce after which the amplitude continues to decrease. The larger h r the more rapid and significant the amplitude decrease. A long time after the interaction of the solitary wave with the abrupt deepening, the maximum elevation of the free surface continues to decrease (almost linearly) with time, at a slow rate. We believe that this phenomenon deserves a thorough experimental investigation, and we hope that it will be undertaken in the future. The novel Hamiltonian Coupled-Mode Theory, proposed by [START_REF] Athanassoulis | New form of the Hamiltonian equations for the nonlinear water-wave problem, based on a new representation of the DtN operator, and some applications[END_REF], [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with nonplanar boundaries[END_REF], is a nonperturbative, coupled-mode approach able to solve fully nonlinear water-wave problems in one or two horizontal dimensions, over varying bathymetry. In this approach, the evolution of the nonlocal Hamiltonian system requires the consecutive solution of a linear coupled-mode system with (x, t)-varying coefficients.

The efficient implementation of the whole numerical scheme heavily relies on the fast and accurate evaluation of these coefficients, which are expressed in terms of the roots of the water-wave dispersion relation with (x, t)-varying frequency parameter.

In this paper, we propose new, highly accurate, (semi)explicit formulae for the evaluation of wave numbers corresponding to evanescent modes. Their derivation is based on the successive application of a low-order iteration procedure (Picard iteration) and higher order methods, based on the general Householder formula for the solution of nonlinear algebraic equations. It is rigorously established that the compound iteration scheme, although retains the same order as that of the underlying high-order method, has a significantly lower error constant. The proposed 2 nd and 3 rd order semi-explicit methods are fairly simple and their first iteration provides quite accurate closed-form expressions for all eigenvalues k n and all values of the frequency parameter µ. Notwithstanding the usefulness of these expressions in all multimodal equations, their impact is really highlighted in connection with the HCMT, for solving fully nonlinear water wave problems over general bathymetry. In order to obtain stable longtime simulations in such demanding cases, computations at machine precision accuracy are necessary, and they are achieved by the present methods with not more than three iterations.

The main conclusions of the present paper are: (i) the 2 nd and 3 rd order methods provide highly accurate results for all shallowness conditions, in contrast with the Newton-Raphson method that diverges in the deep water case, if the initial guess in not judiciously chosen; (ii) the 2 nd order method outperforms the 3 rd order one, in terms of computational time; (iii) The HCMT is efficiently implemented by using the closed-form expressions for the coefficients of the kinematical problem, evaluated by the 2 nd order semi-explicit formulae for the local wavenumbers.

Figure 1 :Figure 2 :

 12 Figure 1: The first three roots of the transcendental Eq. (6b), for three different values of the parameter µ (µ = 0.1, 1.0, 5.0)

Remark 4 . 3 Proposition 4 . 1

 4341 The above-described, two-step (compound) procedure, Eqs.[START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF] and[START_REF] Guo | Simple and explicit solution of wave dispersion equation[END_REF], has improved convergence characteristics in comparison with the singlestep iterations[START_REF] Mckee | Calculation of evanescent wave modes[END_REF] and[START_REF] Newman | Numerical solutions of the water-wave dispersion relation[END_REF]. More precisely, as shown in the next proposition, although approximations (28) and (29) are of second and third order, respectively, the same as the parent single-step iterations, their error constants are significantly smaller. Assume that no error is introduced in all function calculations. Since, by their construction, the iteration formulae (25) and (27) are of second and third order, respectively, there exist positive constants r, R, depending only on n, µ, such that (for appropriate initial guess) the iteration methods (25) and (27) converge and satisfy the error estimates

Figure 3 :

 3 Figure 3: Plot of the factor δ n (µ), as a function of parameter µ, for different values of n ∈ N

Figure 5 :

 5 Figure 5: Logarithmic plot of the relative error of the approximation D n (µ) and the two first iterations of Eq. (28) initiated by D n (µ) for the first 6 roots of Eq.[START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with nonplanar boundaries[END_REF]

Table 1 : 4 . 1 ×

 141 L 2 -error on the free surface elevation for the simulation of strongly nonlinear travelling waves over flat bottomErroratt = 3T λ/h N tot = 3 N tot = 4 N tot = 5 N tot = 6 1(*) 6.0 × 10 -3 1.3 × 10 -3 1.9 × 10 -4 5 10 -3 3.6 × 10 -4 4.6 × 10 -5 9.1 × 10 -5 18 6.2 × 10 -3 3.3 × 10 -4 1.8 × 10 -4 2.6 × 10 -4

Figure 9 :Figure 10 : 5 0Figure 11 :

 910511 Figure 9: (a) Free surface elevation η(x) and log 10 µ(x) for the case λ/h = 1. (b) Relative computational time for the evaluation of κ n versus N X N T for the improved NR and the semi-explicit methods

Figure 13 :Figure 14 :

 1314 Figure13: Comparison of the computed free surface elevation with experimental data (grey bullets) of[START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF]

Figure 15 :

 15 Figure 15: Local Bragg reflection coefficient in the experiment of[START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF]

Figure 16 :

 16 Figure 16: (a) Initial free surface, and bottom surface for the simulation of a solitary wave propagating over an abrupt deepening, with h r = 4m. (b) Values of µ(x, t) =µ 0 (η(x, t) + h(x)) at the initial time, for the same case

Figure 17 :Figure 18 :7

 1718 Figure 17: Snapshots of the free surface elevation for the simulation of a solitary wave propagating over an abrupt deepening, with h r = 4m

Acknowledgement

This research has not been supported by any funding bodies. The authors would like to thank Mr. A. Charalampopoulos for his support in the numerical simulations.

Eqs. [START_REF] Porter | Extensions of the mild-slope equation[END_REF] and [START_REF] Chamberlain | Multi-mode approximations to wave scattering by an uneven bed[END_REF] can be considered either as improved initiations of the iterations [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF] and [START_REF] Guo | Simple and explicit solution of wave dispersion equation[END_REF], or as simple, closed-form approximations of the soughtfor roots κ n . The closed-form approximation

proposed by Chamberlain [START_REF] Guo | Simple and explicit solution of wave dispersion equation[END_REF], is also considered for comparison. Numerical experiments have been performed in order to test the accuracy of the approximations ( 34) - [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF]. "Exact" solutions are obtained using the fzero function of Matlab with initial values given by Eq. [START_REF] Chamberlain | Multi-mode approximations to wave scattering by an uneven bed[END_REF]. The percentage of relative error of A n (µ), B n (µ), C n (µ), as function of µ, is plotted in Figure 4, for the first three values of n. For A n (µ) and C n (µ) the maximum relative error is less than 10%, while for B n (µ) the maximum relative error is less than 1.5%. The error drops significantly as n increases. In all cases, the expression B n (µ) performs better than A n (µ) and C n (µ).

In accordance to the procedure described in Algorithm 2, the quantities nπ, A n (µ) and B n (µ) are three consecutive terms of a sequence converging to κ n . Motivated by this fact, we apply a Shanks' transformation [64, p. 369], also known as Aitken δ 2 -method or Steffensen's procedure [65, p. 103], to produce a possibly better approximation:

Substituting A n (µ) and B n (µ) from Eqs. [START_REF] Porter | Extensions of the mild-slope equation[END_REF], [START_REF] Chamberlain | Multi-mode approximations to wave scattering by an uneven bed[END_REF] into Eq. ( 37), we obtain

Formula ( 38) is "almost the same" as formula [START_REF] Chamberlain | Multi-mode approximations to wave scattering by an uneven bed[END_REF], and it is indeed more accurate than the latter. The maximum relative error in the computation of the first root is about 0.7%; see Figure 5. Again, the maximum relative error reduces significantly for increasing n. Finally, a different approximation is produced by setting 0 κ n = nπ in the third-order formula (29), obtaining:

The error of E n (µ) is plotted in Figure 6. Comparing Figures 5 and6, it is readily seen that E n (µ) behaves better that D n (µ) for large values of µ. Approximations D n (µ) and E n (µ) can be used for the initiation of the iteration formulae [START_REF] Chamberlain | On the solution of the dispersion relation for water waves[END_REF] and [START_REF] Guo | Simple and explicit solution of wave dispersion equation[END_REF], respectively. As shown in Figures 5 and6, they produce highly accurate results after the first iteration and, in general, results

Proof of proposition 1

The proof will be derived for the general two-step recursive formula (compound scheme)

where Eq. ( 40) is a first-order iteration scheme which, for an appropriate 0 x, converges to a simple root x = a and j+1 x = G( j x) is an iteration scheme of order p ≥ 2, converging to the same simple root. The last assumption implies, using a Taylor expansion with respect to x = a, that

Assume further that the functions ϕ, G are sufficiently smooth, and the distances 0 x -a , ϕ ( 0 x) -a are such that the employed schemes converge to a. Since ϕ(a) = G(a) = a and the derivatives G (k) (a) = 0 for all k < p, from Eqns. (A2), following the same procedure used for the derivation of (A3), results to

Thus, since G(ϕ(a)) = a, for any iteration performed via Eq. (A2) it is

where R G = G (p) (a) /p!. In the particular case of Picard iteration defined through Eq. ( 21), we have ϕ(x) = nπ-Arctan (µ/x), hence

From the definition of ε n it is ε n ∈ (0, π/2) and thus

which completes the proof.