A FAST SECOND-ORDER DISCRETIZATION SCHEME FOR THE LINEARIZED GREEN-NAGHDI SYSTEM WITH ABSORBING BOUNDARY CONDITIONS

. In this paper, we present a fully discrete second-order finite-difference scheme with fast evaluation of the convolution involved in the absorbing boundary conditions to solve the one-dimensional linearized Green-Naghdi system. The Pad´e expansion of the square-root function in the complex plane is used to implement the fast convolution. By introducing a constant damping parameter into the governing equations, the convergence analysis is developed when the damping term fulfills some conditions. In addition, the scheme is stable and leads to a highly reduced computational cost and low memory storage. A numerical example is provided to support the theoretical analysis and to illustrate the performance of the fast numerical scheme


Introduction
Under the effect of the gravity, the motion of an irrotational and incompressible fluid is described by the free-surface Euler equations.Because of the complexity of this system, asymptotic models for the water wave problem [19] were developed over the years.In particular, the Green-Naghdi model [13] includes the dispersive effects and writes in the two-dimensional space as (1.1) In the above notations, we designate by  the fluid depth and by ⃗  the depth-averaged horizontal velocity, ()  stands for the derivation of  with respect to the time variable  and v is the material derivative.System (1.1) describes the bidirectional propagation of dispersive water waves in the shallow water regime [19].A one-dimensional simplification of (1.1), linearized around the steady state (, ) := ( 0 , 0) + ( 1 ,  2 ), with |( 1 ,  2 )| ≪ 1, can be derived as the one-dimensional linearized Green-Nagdhi (GN) system [18] ( 1 )  + ( 2 )  = 0, ( 2 )  + ( 1 )  = ( 2 )  ,  ∈ R,  > 0,  1 (, 0) =  1 (),  2 (, 0) =  2 (),  ∈ R, lim ||→+∞  1 (, ) = 0, lim ||→+∞  2 (, ) = 0,  > 0. (1.2) Since (1.2) is set in an unbounded domain, then suitable boundary conditions need to be introduced to get a finite spatial computational domain in view of a discretization.This problem is well-known in the literature and fits into the framework of designing Absorbing Boundary Conditions (ABCs) and Perfectly Matched Layers (when a fictitious layer is added) for systems of PDEs.We refer to [4,8,14,26] for overviews of the various approaches that can be used, their pros/cons and the related discretization aspects and computational difficulties that are met.For the GN system (1.2), a major contribution has been recently achieved by Kazakova and Noble in [18].In this work, exact ABCs for the fully discrete system of linearized Green-Naghdi equations (1.2) was proposed on staggered grids and the stability of the exact ABCs was proved.Nevertheless, and similarly to the one-dimensional linear Schrödinger equation, these absorbing boundary conditions require the expensive evaluation of nonlocal time convolution-type operators at the fictitious boundary points that lead to prohibitive computational costs and memory storage, in particular for long time computations.In addition, instabilities may arise if the evaluation of the ABCs through the -transform is not carefully implemented (see e.g.[1,9,17,23,27] for some examples).
Probably the most emblematic linear dispersive PDE analyzed in the literature for constructing ABCs is the Schrödinger equation.The first scheme with ABCs was introduced by Baskakov and Popov [11] for computational acoustics based on the parabolic equation.It is now well-known that the resulting scheme with ABCs suffered from stability issues [4].Since then, many developments allowed to solve numerous problems related to the one-dimensional case [2,4,5,8,9,12,16,17], and extensions of some of the methods were proposed for higherdimensional and nonlinear problems [3,6,7,10,15,20,23].In addition, some of these contributions also include the numerical analysis of the schemes and the derivation of fast and stable evaluation schemes of the nonlocal half-order time derivative operator arising in the definition of the ABCs.The case of the GN system remains much less studied [18] and therefore requires still more understanding.The aim of the present paper is to contribute to the development of ABCs for (1.2) by deriving alternative efficient formulations, and to carefully analyze the convergence of the fully discrete scheme.
To this end, we consider the recent approach introduced by Li et al. [21] for the one-dimensional Schrödinger equation and extend it to the GN system (1.2).More precisely, to overcome the numerical instability, a convergent numerical method, integrating a fast evaluation of the exact ABC, is proposed for solving the Cauchy problem (1.2).To this end, the GN system is first reformulated in Section 2 under an equivalent form by introducing a constant damping term, and then a modified Crank-Nicolson scheme is built in Section 3 to discretize the equivalent problem according to the time variable.More specifically, a semi-discrete ABC for the temporally discretized problem is derived for the Crank-Nicolson scheme based on the -transform.Then, a second-order finite difference-scheme for the full spatial discretization is considered.A fast algorithm is introduced in Section 4 to approximate the discrete convolution kernel involved in the exact semi-discrete ABC by using the Padé rational expansion of the square-root function [22].The damping parameter and the number of Padé terms are chosen to maintain the convergence order of the resulting discrete scheme [21].In Section 5, a convergence analysis for the proposed numerical method is developed, showing that the scheme is second-order both in space and time.A numerical example illustrates the properties of the scheme in Section 6. Section 7 finally concludes the paper.

Exact ABCs for the one-dimensional linearized Green-Naghdi system
We consider the initial-value problem for the linearized GN system set on the whole space (2.1) Let us introduce the new unknown functions:   (, ) =  −   (, ),  = 1, 2, where  > 0 is a parameter used to later control the stability of the fast algorithm.It is straightforward to check that the function   (, ) solves the following initial-value problem: To obtain exact ABCs for (2.2), we first consider the following exterior problem on the semi-infinite interval [ + , +∞): ) ) (2.3d) The Laplace transform in time (denoted by ̂︀ () for a function ()) of system (2.3) yields where C + stands for the right half-space of complex numbers with positive real part.After some manipulations, one has: ( + ) 2 ̂︀  2 (, ) = (( + ) 2 + 1)  ̂︀  2 (, ),  = 1, 2, which means that the general solution writes where √ • denotes the square-root with nonnegative real part and The behavior of the wave at infinity (2.6) implies that  2 () = 0.By differentiating (2.7), we obtain whose inverse Laplace transform yields an exact absorbing boundary condition for In (2.10), ℒ −1 denotes the inverse Laplace transform with respect to the variable  and we have A similar boundary condition can be derived at a left point  − <  + : These boundary conditions degenerate to the boundary conditions of the wave equation when  → 0. In view of (2.10) and (2.11), the solution of (2.2) is the same as the solution of the following problem in the bounded domain ( − ,  + ) where   denotes the outward normal derivative at the boundary points  ± .

Discretization of the one-dimensional GN system with exact semi-discrete ABC
In this section, we discretize the one-dimensional linearized Green-Naghdi system in time by using the Crank-Nicolson scheme and derive an associated exact semi-discrete ABC.Then, we propose a second-order finitedifference scheme for the spatial discretization.To this end, we first introduce the notations related to the -transform in the following subsection.

The 𝒵-transform of a sequence of functions
Let us consider a Hilbert space ℋ equipped with an inner product (•, •) ℋ and induced norm ‖ • ‖ ℋ .We introduce the semi-infinite sequence spaces: and with the inner product: (, ) The following Parseval's identity holds: In the above,  stands for the normalized Haar measure on the unit circle D of the complex plane, and (d) = 1 2 d through the change of variable  =   , with  ∈ [−, ).For a sequence  = {  } ∞ =0 ∈ ℓ 2 (ℋ), we define the operator  by:  = { +1 } ∞ =0 .The average operator  and the forward difference quotient operator   with step  are given by  = ( + )/2 and   = ( − )/ , respectively.We also introduce the following notations:   = ()  ,   = ()  and     = (  )  .

Exact ABCs for the semi-discretized one-dimensional linearized GN system
Let  > 0 be the uniform time step such that   =  =   , with  the maximal time of computation.Let us set the discrete times as:   =  , 0 ≤  ≤  .System (2.2) is semi-discretized in time following where    () ≈   (,   ), for  = 1, 2. We now assume that the initial data  1 and  2 are compactly supported in the interval [ − ,  + ].On [ + , +∞), the semi-discrete problem (3.4) reduces to Let us denote by ̃︀ (, ) the -transform of the sequence {  ()} ∞ =0 .Applying the -transform to (3.5), we obtain , setting, for all  > 0 and  > 0, The condition at infinity, i.e. lim →+∞ ̃︀  2 (, ) = 0, implies that ̃︀  2 = 0, leading to which corresponds to the semi-discretization of (2.9) at  =  + .Note that the function is analytic in the unit disk D. Thus, it admits a power series expansion  (3.4), originally defined on the whole space, can be reduced to the following semi-discrete problem on a bounded domain: (3.12)

Spatial discretization
Let  be a positive integer and ℎ = ( + −  − )/ the uniform mesh size.We define the mesh points: and  +1/2 =  0 + ( + 1/2)ℎ, for  = 0, 1, • • • ,  , where  0 and   +1 are two ghost points.In the time-stepping scheme (3.12), we use ( 2 )   to denote the numerical approximation of   2 (  ), with 0 ≤  ≤  + 1, and ( 1 )   to define that of )︂ , respectively.The linear operator which maps the ( + 2)-dimensional vector  = ( 0 , • • • ,   +1 ) to the  -dimensional vector ( 1 , • • • ,   ) will be denoted by .In addition, we introduce the Neumann and Dirichlet data associated with the ( + 2)-dimensional vector  as Let us define the inner product for two  -dimensional vectors and finally the inner product for two ( + 1)-dimensional vectors In the above expressions,  denotes the complex conjugate of a complex number .The induced norms are such that: Let us now introduce a second-order spatial discretization △ ℎ , which maps the ( + 2)-dimensional vector  to the  -dimensional vector space as )︂ • Thus, we have the discrete integration by parts formula Now, we define the vector ∇ ℎ ℋ  1 by where ℋ is any operator that only applies in the time direction.Similarly, the vector ∇ ℎ ℋ  2 is such that We can also introduce a vector △ ℎ ℋ  2 by Then, it is easy to see that the following identities hold: ) and changing the continuous operator   with its discrete analogue △ ℎ , we obtain the following fully discrete finite-difference scheme (3.14) In fact, at  ± , we have    2 ( ± , ) = ±   2 ( ± , ).Therefore, we can write that For (3.14), the expression of  1 at  1 can be written as Thus, the boundary condition for  1 can be supplied by  2 due to the staggered grid.

𝑛
In this section, we introduce a fast algorithm for approximating the discrete convolution product ( *  ±  2 )  arising in (3.14).The stability of the proposed fast algorithm will be analyzed in the next section.

Rational approximation of the convolution quadrature
In [22], for a nonnegative integer  > 0 and Re() ≥ −1, the Padé approximation of the function √ 1 +  can be expressed as where the coefficients are given by Based on this Padé approximation, a rational approximation of the square-root function √  on the closed right half complex plane can be written as ≡   (), Re() ≥ 0.
Thus, we deduce For all  > 0,  > 0, and for () defined by (3.6), we can introduce the rational approximation ̃︀  () () of the symbol ̃︀  () as We denote by  () * the convolution operator analogously defined as (3.11) by replacing the convolution coefficients with the series expansion coefficients of the function ̃︀  () ().By considering the rational approximation  () * of  * in (3.14), we obtain the following fully discrete scheme: In fact, equation (4.5) can be solved by the fast algorithm described in the next subsection.

Fast evaluation of
By applying (4.1) to (4.2), we obtain the sequence of equalities setting Thus, ̃︀  () () can be uniformly rewritten as can be implemented by a fast convolution in (4.5).For fixed , by defining with 1 ≤  ≤ 2 + 1, we derive that, Thus, the boundary term ( can be written in the form of fast convolution by (4.8), namely, From (4.9), we can see that the computational cost at the -th time step is (2 + 1), rather than () which is the computational cost of the convolution.Thus, the overall computational cost for all the  time steps is ((2+1)) instead of ( 2 ).Therefore, for large time steps , the total computational cost is greatly reduced.Such a cost is much lower than the one proposed in [18] for the linearized Green-Naghdi equation.It is similar to the techniques introduced for example in [21,27] to optimize both the memory storage and computational cost when evaluating convolutions.

Properties of the rational approximation ̃︀ 𝒯 (𝑚) (𝑧)
Let us now prove some properties of  ()  used to prove the error estimates.Proposition 4.1.Let us assume that the condition  ≥ 1 √ 2 is satisfied, the time step  is small enough and  is sufficiently large, i.e. it fulfills , for some  ∈ with Then, the following inequalities hold Re for any complex-valued sequence  = {  } ∞ =0 such that  0 = 0.
Before proving Proposition 4.1, we need to state a few lemmas.Let us introduce Then, by using (3.8), one can prove that the symbol ̃︀  () satisfies the following inequalities.
Recalling that () = with cos 2 for  ∈ (−∞, +∞).Using the above expression of √︀ (), we have where the last inequality is a consequence of: (1 − ) we see that the minimum value of leading to (4.18).Next, we have the sequence of equalities arg with For  ∈ [0, +∞), we obtain Similarly, we can write that arctan Let us recall that   () is defined by (4.1).Then, the following result was proved in [22].
Let us now consider the following lemma.

Error estimates of the scheme
Let us define the two error vectors where ( 1 ,  2 ) is solution to (2.12) and (  1 ,   2 ) is the solution to the discrete system (4.3)-(4.6).We first give the main result concerning the error estimate.
Theorem 5.1.Let us assume that the solutions  1 (, ) and  2 (, ) to system (2.1), or equivalently, the solutions  1 (, ) and  2 (, ) of (2.2), are sufficiently smooth.Let us suppose that  ≥ 1 √ 2 ,  is small enough and that  satisfies (4.10), with  and  given by (4.11).Then, we have the error bound max It is straightforward to check that the error vector (  1 ,   2 ) satisfies ) where and   ± are the interior/boundary truncation error vectors/numbers according to the time and space discretizations, i.e. ) The proof of Theorem 5.1 is presented in the next two subsections as a consequence of Propositions 5.2 and 5.3.

Estimate for truncation errors
Let us first prove the estimate for the truncation errors of the boundary and interior schemes.
Proposition 5.2.Under the conditions of Theorem 5.1, we have the following error estimate ,  being a strictly positive constant.
Proof of Proposition 5.2.The proof is separated into three estimates which are summed up at the end to prove the result.

Estimate of
Here, we prove (5.10) We divide the proof into two steps.

Error estimates
Let us state the error estimate for system (4.
where   is a constant depending on  .
Proof of Proposition 5.3.Due to taking the real part of the inner product between the left hand side of (5.2) and (5.30) using the inequality Summing up over , we obtain )︁ which leads to (5.33) The left hand side of (5.33) can be written as (5.34) By applying the discrete Green's formula (3.13), the boundary conditions (4.5), and (4.13), the right hand side of (5.33) can be rewritten following This leads to Again taking the inner product of (5.3) with (  + )△ ℎ   2 and then the real part, we obtain Re(( (5.37) The right hand side of (5.37) can be written as (5.38) Now by using the integration by part (3.13), the boundary conditions (4.5), and (4.14), the left term of (5.37) reads (5.39) Combining (5.38) and (5.39), we deduce from which we derive that Summing up the index  and using a summation by parts in time for (5.40) with  and  given by (4.11).Therefore for  fixed (with   =  ), the total computational cost to efficiently evaluate the convolution is ( ) = ( log( )) [21] by (4.9).
Example 6.1.To demonstrate the performance of our numerical scheme, we first consider a Gaussian initial distribution for the free-surface elevation and zero distribution for velocity, i.e.
We report on Figure 1a the amplitude of the numerical surface elevation  1 with ABCs on the computational domain.We set  =  = 640.The reference solution )︀ is computed for  = ℎ = 1 2560 with a CN scheme, in a very large computational domain to avoid any influence of the boundary condition.We also draw on Figure 1b the amplitude of the velocity  2 .Furthermore, on Figures 1c and 1d, we plot the error log 10 | ref  −   |,  = 1, 2, in the domain of computation.As observed, the numerical and reference solutions are very similar, the error being related to the second-order accuracy of the scheme.There is no reflection related to the absorbing boundaries.Finally, we report on Figures 1e and 1f the  ∞ -error max ∈[0,1] |  (, 1) −  ref  (,  )|,  = 1, 2, at final time  = 1, when recursively doubling the parameters of the discretization grid, i.e.  =  from 80 to 640.A second-order convergence rate in  ∞ -norm is observed.Example 6.2.To observe the dispersive behavior of the GN system, we consider now the following initial distribution for the free-surface elevation  1 (, 0) = exp(−400( − 0.5)) sin(20), ( and set to zero the initial velocity, i.e.  2 (, 0) = 0.The initial data  1 is small outside the spatial computation domain [0, 1] so that it can be considered as numerically compactly supported.We fix  = 10 −3 and the maximal evolution time  = 1.
We first plot on Figure 2a  )︀ is computed with  = ℎ = 1 6400 in a large enough domain with a CN scheme so that we do not see the effect of the boundary condition.Similarly, we draw on Figure 2b the amplitude of the velocity  2 .In addition, on Figures 2c and 2d, we report the error log 10 | ref  −   |,  = 1, 2, in the computational domain.As it can be observed, the numerical and reference solutions superpose, up to the second-order error of the scheme.No spurious reflection can be detected near the absorbing boundaries.For completeness, we plot on Figures 2e and 2f the  ∞ -error max ∈[0,1] |  (, 1) −  ref  (,  )|,  = 1, 2, at final time  = 1, when recursively doubling the parameters of the discretization grid, i.e.  =  from 160 to 1280.We observe a second-order convergence rate in  ∞ -norm.
To end, the CPU time (log scale, sec.) vs.  (log scale) is reported on Figure 3 for the fast evaluation of the convolution operator, fixing the number of spatial grid points to  = 160.The total number of time steps  increases from  = 1.2 × 10 5 to  = 7.2 × 10 5 , with step 1.2 × 10 5 .We observe a slope equal to 1, showing that the cost is linear according to log( ), i.e. as ( log  ) for the computational time.
We draw on the left of Figure 4 the amplitude of the velocity  2 until  = 2 to display the stability for a longer simulation time.In addition, on the right of Figure 4, we report the error log 10 | ref 2 −  2 |, in the computational domain.One can see the wave goes out of the computational domain.The error also increases due to the relation     (, ) =   (, ).

Conclusion
The one-dimensional linearized Green-Naghdi system in an unbounded domain was reformulated into an initial boundary-value problem in a bounded domain with transparent boundary conditions.A fully discrete    Crank-Nicolson finite-difference method was proposed to solve the reformulated initial boundary-value problem but with an exact semi-discrete ABC.A fast convolution algorithm is introduced to deal with the convolutions for the exact semi-discrete ABC by using the Padé rational expansion.A criterion determining the damping term was proposed to guarantee the convergence.In this case, it was proved theoretically that the corresponding numerical scheme can achieve a second-order accuracy both in space and time.A numerical example validates the accuracy and efficiency of the proposed numerical method.
The problem that still needs to be solved is that the damping term  − should satisfies the stability condition  ≥ 1 √ 2 .For a small dispersion parameter , the damping term  − which decays too fast will bring some numerical errors due to the relation     () =   ().We will deal with this problem in the forthcoming paper.Extensions to higher-dimensional problems still need further investigations.Finally, the variable coefficients and nonlinear cases of the Green-Naghdi system remain open problems as well as the case of the two-layer

Figure 3 .
Figure 3. Computational time for the evaluation of the convolution by the fast algorithm vs.  (for  = 160).