Afaf Bouharguane
email: afaf.bouharguane@math.u-bordeaux1.fr

Benjamin Melinand
email: benjamin.melinand@math.u-bordeaux1.fr

A SPLITTING METHOD FOR DEEP WATER WITH BATHYMETRY

In this paper we derive and prove the wellposedness of a deep water model that generalizes the Saut-Xu system for nonflat bottoms. Then, we present a new numerical method based on a splitting approach for studying this system. The advantage of this method is that it does not require any low pass filter to avoid spurious oscillations. We prove a local error estimate and we show that our scheme represents a good approximation of order one in time. Then, we perform some numerical experiments which confirm our theoretical result and we study three physical phenomena : the evolution of water waves over a rough bottom; the evolution of a KdV soliton when the shallowness parameter increases; the homogenization effect of rapidly varying topographies on water waves.

INTRODUCTION

1.1. Presentation of the problem. Understanding the influence of the topography on water waves is an important issue in oceanography. Many physical phenomena are linked to the variation of the topography : shoaling, rip currents, diffraction, Bragg reflection. Since the direct study on the Euler equations is quite involved, several authors derived and justified asymptotic models according to different small parameters. A usual way to derive asymptotic models is to start from the Zakharov/Craig-Sulem-Sulem formulation [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF][START_REF] Craig | Numerical simulation of gravity waves[END_REF][START_REF] Craig | Nonlinear modulation of gravity waves: a rigorous approach[END_REF], which is a good formulation for irrotational water waves, and to expand the Dirichlet-Neumann operator. Then, in the shallow water regime for example, several models were obtained like the Saint-Venant equations or the Green-Naghdi or Boussinesq equations, see [START_REF] Alvarez-Samaniego | Large time existence for 3D water-waves and asymptotics[END_REF][START_REF] Lannes | The water waves problem[END_REF], [START_REF] Iguchi | A shallow water approximation for water waves[END_REF] for instance. The present paper addresses the influence of the bathymetry in deep water, in the sense explained below.

In this paper, a denotes the typical amplitude of the water waves, L the typical length, H the typical height and a bott the typical amplitude of the bathymetry. Then, we introduce three parameters : ε = a H the nonlinearity parameter, µ = H 2 L 2 the shallowness parameter and β = a bott H the bathymetric parameter. We recall that assuming µ small leads to shallow water models. In deep water, which is typically the case when µ is of order 1, it is quite common to assume that the steepness parameter ε √ µ = a L is small. The first asymptotic model with a small steepness assumption was derived by Matsuno in 2D for a flat and non-flat bottom and weakly transverse 3D water waves [START_REF] Matsuno | Nonlinear evolutions of surface gravity waves on fluid of finite depth[END_REF][START_REF] Matsuno | Nonlinear evolution of surface gravity waves over an uneven bottom[END_REF]. Then, Choi extended this result in 3D for flat bottom [START_REF] Choi | Nonlinear evolution equations for two-dimensional surface waves in a fluid of finite depth[END_REF]. Finally, Bonneton and Lannes gave a 3D version in the case of a non-flat bottom [START_REF] Bonneton | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. It is important to notice that these models are only formally derived. It is proven in [START_REF] Alvarez-Samaniego | Large time existence for 3D water-waves and asymptotics[END_REF] that smooth enough solutions to theses models are close to the solutions of the water waves equations but, to the best of our knowledge, the wellposedness of the Mastuno equations, even in the case of a flat bottom, is still an open problem. This system could be illposed (see Ambrose, Bona and Nicholls [START_REF] Ambrose | On ill-posedness of truncated series models for water waves[END_REF]). To avoid this difficulty, Saut and Xu ( [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF]) developed an equivalent system to the Matsuno system which is consistent with the water waves problem and with the same accuracy. Then, they proved that this new system is wellposed. However, this model is for a flat bottom. In this paper, we derive (see Section 2), use, and prove the wellposedness of a generalization of the Saut-Xu system with a non-flat bottom which is the following system (1.1)

       ∂ t ζ-H µ v+ε √ µ 1 2 v∂ x ζ+ 1 2 H µ (v∂ x H µ ζ)+H µ (ζ∂ x H µ v)+ζ∂ x v = β √ µ∂ x (B µ v) ∂ t v + ∂ x ζ + 3ε √ µ 2 v∂ x v - ε √ µ 2 ∂ x ζH µ ∂ x ζ - ε √ µ 2 vH 2 µ ∂ x v = 0,
where (see Subsection 1.2 for the notations)

H µ = - tanh( √ µD) D ∂ x and B µ = sech( √ µD) (b sech( √ µD) • ) .
Many authors developed numerical approaches to study the impact of the bottom on water waves, see for instance [START_REF] Mei | Resonant reflection of surface water waves by periodic sandbars[END_REF], [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF], [START_REF] Smith | An operator expansion formalism for nonlinear surface waves over variable depth[END_REF], [START_REF] Guyenne | Numerical simulation of solitary waves on plane slopes[END_REF], [START_REF] Cathala | Asymptotic shallow water models with non smooth topographies[END_REF], [START_REF] Nachbin | A terrain-following Boussinesq system[END_REF], [START_REF] Hamilton | Differential equations for long-period gravity waves on fluid of rapidly varying depth[END_REF], [START_REF] Bonneton | Recent advances in serre-green naghdi modelling for wave transformation, breaking and runup processes[END_REF], [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive green-naghdi model[END_REF]). However to the best of our knowledge, when one works with deep water, there is no convergence result in the literature. After the original work of Craig and Sulem ( [START_REF] Craig | Numerical simulation of gravity waves[END_REF]) and the paper of Craig et al. ( [START_REF] Craig | Hamiltonian long-wave expansions for water waves over a rough bottom[END_REF]), Guyenne and Nicholls ([19]) developed a numerical method based on a pseudospectral method and a fourth-order Runge-Kutta scheme for the time integration. The linear terms are solved exactly whereas the nonlinear terms are viewed as source terms. Their approach has been developed for the whole water waves equations but we could easily adapt it to our system. However with their scheme, we observe spurious oscillations in the wave profile that lead to instabilities. These errors seem to appear when the nonlinear part is evaluated via the Fourier transform. This is the aliasing phenomenon. Guyenne and Nicholls also observe these oscillations and, to fix it, they apply at every time step a low-pass filter. The scheme that we propose in this paper avoids this low-pass filter.

We present a new numerical method based on a splitting approach for studying nonlinear water waves in the presence of a bottom. We remark that the Saut-Xu system contains a dispersive part and a nonlinear transport part. Thus, the splitting method becomes an interesting alternative to solve the system since this approach is commonly used to split different physical terms, see for instance [START_REF] Ropp | Stability of operator splitting methods for systems with indefinite operators: advectiondiffusion-reaction systems[END_REF]. We also motivate our decomposition by the fact that, due to the pseudodifferential operator, some terms in the dispersive part may be computed efficiently using the fast Fourier transform. The transport part is computed by a Lax-Wendroff method. Various versions of the splitting method have been developed for instance for the nonlinear Schrodinger, the viscous Burgers equation, Kortewegde-Vries equations [START_REF] Carles | On Fourier time-splitting methods for nonlinear Schrödinger equations in the semiclassical limit[END_REF][START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF][START_REF] Sacchetti | Spectral splitting method for nonlinear Schrödinger equations with singular potential[END_REF][START_REF] Taha | Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrödinger equation[END_REF]. Thanks to this splitting, we only use a pseudospectral method for the nonlocal terms (contrary to [START_REF] Craig | Numerical simulation of gravity waves[END_REF][START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]), which limits the aliasing phenomenon and allows us to avoid a low-pass filter.

We denote by Φ t the nonlinear flow associated to the Saut-Xu system (1.1), Φ t A and Φ t D , respectively, the evolution operator associated with the transport part (see equation (3.1)) and with the dispersive part (see (3.2)). We consider the Lie formula defined by

(1.2) Y t = Φ t A Φ t D .
The Saut-Xu system (1.1) is a quasilinear system. This implies derivatives losses in the proof of the convergence. In Theorem 4.6, we show that the numerical solution converges to the solution of the Saut-Xu system (1.1) in the

H N + 1 2 × H N -norm for initial data in H N + 1 2 -7 × H N -7 (R), where N ≥ 7.
Notice that it is not hard to generalize the present work to the Lie formula Φ t D Φ t A . We also make the choice to prove a convergence result for a Lie splitting but our proof can be adapted to a Strang splitting or a more complex one. Finally, notice that our scheme can be used for other equations.

The paper is organised as follows. In the next section, we extend the Saut-Xu system by adding a topography term and we prove a local wellposedness result. We also show that the flow map Φ t is uniformly Lipschitzean. In section 3, we split the problem and we give some estimates on Φ t A and Φ t D . In Section 4, we prove a local error estimate and we show that the Lie method represents a good approximation of order one in time (Theorem 4.6). Finally, in Section 5, we perform some numerical experiments which confirm our theoretical result and we illustrate three physical phenomena : the evolution of water waves over a rough bottom; the evolution of a KdV soliton when the shallowness parameter increases and the homogenization effect of rapidly varying topographies on water waves.

Notations and assumptions.

• x denotes the horizontal variable and z the vertical variable. In this paper, we only study the 2D case (x ∈ R). • We assume that

(1.3) 0 ≤ ε, β ≤ 1, ∃µ max > µ min > 0, µ max ≥ µ ≥ µ min .
• We denote δ = max(ε, β).

• We denote Λ = (1-∂ 2 x ) 1/2 and H s (R) = u ∈ L 2 (R), ||u|| H 2 = ||Λ s u|| L 2 < ∞ the usual Sobolev space for s ≥ 0. • Let f ∈ C 0 (R) and m ∈ N such that f 1+|x| m ∈ L ∞ (R). We define the Fourier multiplier f (D) : H m (R) L 2 (R) as ∀u ∈ H m (R) , f (D)u(ξ) = f (ξ) u(ξ).
• D denotes the Fourier multiplier corresponding to ∂x i . • We denote by C(c 1 , c 2 , ...) a generic positive constant, strictly positive, which depends on parameters c 1 , c 2 , • • • .

THE SAUT-XU SYSTEM

In this part, we extend the Saut-Xu system ( [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF]) for a non-flat bottom. Then, we give a wellposedness result that generalizes the one of Saut and Xu.

The Matsuno system, which is a full dispersion model for deep waters, is an asymptotic model of the water waves equations with an accuracy of order O δ 2 . Bonneton and Lannes [START_REF] Bonneton | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF] formulated it in the following way in the presence of a non flat topography (2.1)

∂ t ζ -1 √ µν H µ v + ε ν (H µ (ζ∂ x H µ v) + ∂ x (ζv)) = β ν ∂ x (B µ v) ∂ t v + ∂ x ζ + ε ν v∂ x v -ε √ µ∂ x ζH µ ∂ x ζ = 0, where ζ = ζ(t, x) is the free surface, v = v(t, x) is the horizontal velocity at the surface, ν = tanh( √ µ) √ µ
, H µ and B µ are Fourier multipliers,

H µ = - tanh( √ µD) D ∂ x and B µ = sech( √ µD) (b sech( √ µD) • ) ,
and -1 + βb is the topography. (We erased the fluid part)

In [START_REF] Alvarez-Samaniego | Large time existence for 3D water-waves and asymptotics[END_REF], Alvarez-Samaniego and Lannes show that this model is consistent with the Zakharov/Craig-Sulem-Sulem formulation when β = 0 and it is not painful to generalize their result to the case when β = 0. In [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF], Saut and Xu obtained a new model with the same accuracy with the Matsuno system thanks to a nonlinear change of variables. Notice that this change of variables is inspiblack by [START_REF] Bona | Long wave approximations for water waves[END_REF].

The advantage of this model is that they proved a local wellposedness on large time for this new model. We follow their approach. We define new variables

(2.2) v = v + ε √ µ 2 vH µ ∂ x ζ and ζ = ζ - ε √ µ 4 v 2 .
Then, up to terms of order O δ 2 , ζ and v satisfy (we omit the tildes for the sake of simplicity)

(2.3)        ∂ t ζ+ ε ν - ε √ µ 2 v∂ x ζ- 1 √ µν H µ v+ ε ν 1 2 H µ (v∂ x H µ ζ)+H µ (ζ∂ x H µ v)+ζ∂ x v = β ν ∂ x (B µ v) ∂ t v + ε ν + ε √ µ 2 v∂ x v + ∂ x ζ - ε √ µ 2 ∂ x ζH µ ∂ x ζ - ε 2ν vH 2 µ ∂ x v = 0.
Since our motivation is the study of water waves in deep water (µ close to 1), we assume that ν = 1 µ . Hence we study the following system which is a variable bottom analog of the system of Saut and Xu

       ∂ t ζ-H µ v+ε √ µ 1 2 v∂ x ζ+ 1 2 H µ (v∂ x H µ ζ)+H µ (ζ∂ x H µ v)+ζ∂ x v = β √ µ∂ x (B µ v) ∂ t v + ∂ x ζ + 3ε √ µ 2 v∂ x v - ε √ µ 2 ∂ x ζH µ ∂ x ζ - ε √ µ 2 vH 2 µ ∂ x v = 0,
In the following, we denote U = (ζ, v) t and we define the energy of the system for N ∈ N by (2.4)

E N (U) = 1 √ µ Λ N ζ 2 2 + |D| 1 2 Λ N ζ 2 2 + |v| 2 H N ,
where Λ = 1 + |D| 2 and D = -i∇. We also denote by E N µ the energy space related to this norm.

Remark 2.1. Notice that if µ satisfies condition (1.3), the energy E N is equivalent to the

H N + 1 2 × H N -norm.
The main result of this section is the following local wellposedness result.

Theorem 2.2. Let N ≥ 2, U 0 ∈ H N + 1 2 (R) × H N (R) and b ∈ L ∞ (R).
We assume that ε, β, µ satisfy Condition (1.3) and

|U 0 | H N + 1 2 ×H N + |b| L ∞ ≤ M.
Then, there exists a time T 0 = T 0 M, 1 µ min , µ max independent of ε, µ and β and a unique solution U ∈ C 0, T 0 δ , E N µ of the system (1.1) with initial data U 0 . Furthermore, we have the following energy estimate, for all t ∈ 0, T 0 δ ,

E N (U(t, •)) ≤ e δC 0 t E N (U 0 ) , where C 0 = C M, 1 µ min , µ max .
Proof. We refer to Paragraph IV in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF] for a complete proof and we focus only on the bottom contribution. For 0 ≤ α ≤ N , we denote

U (α) = (∂ α x ζ, ∂ α x v). Then, applying ∂ α x to System (1.1), we get ∂ t U (α) + LU (α) + ε √ µ 2 1 {α =0} B[U]U (α) = β √ µ (∂ x ∂ α x (B µ v) , 0) t + ε √ µG α ,
where

L = 0 -1 √ µν H µ ∂ x 0 , B[U] = H µ (vH µ ∂ x • ) + v∂ x H µ ( • H µ ∂ x ζ) -∂ x ζH 2 µ -∂ x ζH µ ∂ x -H µ ∂ x ζ∂ x 3v∂ x -vH 2 µ ∂ x ,
and

G α = (G α 1 , G α 2 ) t with G α 1 = ∂ α x g(ζ, v) - 1 2 1≤γ≤α-1 C γ α H µ (∂ γ x vH µ ∂ 1+α-γ x ζ) + ∂ γ x v∂ 1+α-γ x ζ - 1 2 ∂ x ζ(H 2 µ + 1)∂ α x v G α 2 = 1 2 1≤γ≤α-1 C γ α ∂ 1+γ x ζH µ ∂ 1+α-γ x ζ + 1≤γ≤α C γ α - 3 2 ∂ γ x v∂ 1+α-γ x v + 1 2 ∂ γ x vH 2 µ ∂ 1+α-γ x v
where

g(ζ, v) = -[H µ , ζ]H µ ∂ x v -ζ(H 2 µ + 1)∂ x v.
Then we can show, as in Paragraph IV. B in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF] (see the paragraph called Estimate on G α ) that

(2.5) |G α | 2 + |D| 1 2 G α 2 ≤ C 1 µ min E N (U) .
Like Saut and Xu we define a symmetrizer for

L + B[U] (2.6) S = D tanh( √ µD) 0 0 1 . Notice that (S •, •) is a norm equivalent to √ E 0 .
Then, as in Paragraph IV. B in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF] (see the paragraph called Estimate on II), we get

(2.7) L + ε √ µ 2 1 {α =0} B[U] U (α) , SU (α) ≤ εC 1 µ min , µ max E 2 (U)E N (U) .
Furthermore, for the bottom contribution, we easily get

(2.8) D tanh( √ µD) ∂ α x ζ, ∂ x ∂ α x sech ( √ µD) (b sech ( √ µD) v) ≤ C 1 µ min |b| ∞ E N (U) .
Finally, from Equations (2.5), (2.7), (2.8), we obtain,

E N (U) ≤ E N (U 0 ) + δC 1 µ min , µ max t 0 E N (U) 3 2 + E N (U) (s)ds,
and there exists a time T > 0, such that, for all t ∈ 0, T δ ,

E N (U(t, •)) ≤ C 1 µ min , µ max , E N (U 0 ) .
The energy estimate follows from the Gronwall Lemma.

In order to use a Lady Windermere's fan argument (a well-known telescopic identity used to relate global and the local error), to prove the convergence of the numerical scheme, we need a Lipschitz property for the flow of the Saut-Xu system (1.1). We first give a control of the differential of the flow with respect to the initial datum.

Proposition 2.3. Let N ≥ 2, V 0 ∈ H N + 1 2 (R) × H N (R) , U 0 ∈ H N +1+ 1 2 (R) × H N +1 (R)
, and b ∈ L ∞ (R). We assume that ε, β, µ satisfy Condition (1.3) and

|V 0 | H N + 1 2 ×H N + |U 0 | H N +1+ 1 2 ×H N +1 + |b| L ∞ ≤ M.
Then, there exists a time T = T M, 1 µ min , µ max independent of the parameters ε, µ and β such that Φ t (U 0 ) • (V 0 ) exists on 0, T δ . Furthermore, we have, for all 0 ≤ t ≤ T δ ,

Φ t (U 0 ) • (V 0 ) H N + 1 2 ×H N ≤ C 1 µ min , µ max , M |V 0 | H N + 1 2 ×H N .
Proof. We denote by U(t) = (ζ(t), v(t)) the solution of the Saut-Xu system (1.1) with initial data U 0 . We denote also (η(t),

w(t)) = Φ t (U 0 ) • (V 0 ). Then, (η, w) satisfy the following system (2.9) ∂ t η w +L η w +ε √ µN [(ζ, v)]∂ x η w +ε √ µN [(η, w)]∂ x ζ v = β √ µ (∂ x (B µ w) , 0) t ,
where

L = 0 -1 √ µν H µ ∂ x 0 and N [(ζ, v)] = 1 2 H µ (vH µ •) + 1 2 v H µ (ζH µ •) + ζ -1 2 ∂ x ζH µ 3 2 v -1 2 vH 2 µ . For 0 ≤ α ≤ N , we denote V (α) = (∂ α x η, ∂ α x w).
Then, applying ∂ α x to System (2.9), we get

∂ t V (α) +LV (α) + ε √ µ 2 1 {α =0} B[U]V (α) + B[V]∂ α x U = √ µβ ∂ x ∂ α x (B µ w) 0 + ε √ µJ α ,
where

B[U] = H µ (vH µ ∂ x • ) + v∂ x H µ ( • H µ ∂ x ζ) -∂ x ζH 2 µ -∂ x ζH µ ∂ x -H µ ∂ x ζ∂ x 3v∂ x -vH 2 µ ∂ x , J α = -∂ α x N [(ζ, v)]∂ x η w + N [(η, w)]∂ x ζ v + 1 2 B[U]V (α) + B[V]∂ α x U .
Then, we can show, as in Paragraph IV. B in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF], that (2.10)

|J α | 2 + |D| 1 2 J α 2 ≤ ε √ µC 1 µ min E N (V) .
We recall that we can symmetrize L thanks to

S = D tanh( √ µD) 0 0 1 .
We define the energy associated to this symmetrizer

F α (V) = D tanh( √ µD) ∂ α x η 2 2 + |∂ α x w| 2 2 , and F N (V) = 0≤α≤N F α (V) .
We have, for α = 0,

d dt F α (V) = ε √ µ J α , SV (α) - ε √ µ 2 B[U]V (α) , SV (α) + B[V]∂ α x U, SV (α) + β √ µ ∂ x ∂ α x (B µ v) , SV (α) = I + II + III + IIII.
We can estimate I thanks to estimate (2.10) and II as in Paragraph IV. B in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF]. For IIII, we can proceed as in the previous theorem. For III, we get, thanks to Proposition A.1,

|III| ≤ ε √ µ |(ζ, v)| H N +1+ 1 2 ×H N +1 |(η, w)| H N + 1 2 ×H N
Then, we obtain

d dt F N (V) ≤ δ √ µC(M ) F N (V) + F N (V) ,
and the result follows.

Proposition 2.4. Let N ≥ 2, U 0 , V 0 ∈ H N +1+ 1 2 × H N +1 (R) and b ∈ L ∞ (R).
We assume that ε, β, µ satisfy Condition (1.3) and

|V 0 | H N +1+ 1 2 ×H N +1 + |U 0 | H N +1+ 1 2 ×H N +1 + |b| L ∞ ≤ M.
Then, there exists a time T independent of ε, µ and β and two unique solutions U, V of the system (1.1) on 0, T δ with initial data U 0 and V 0 . Furthermore, we have the following Lipschitz estimate, for all 0 ≤ t ≤ T δ ,

(2.11) |U(t, •) -V(t, •)| H N + 1 2 ×H N ≤ K |U 0 -V 0 | H N + 1 2 ×H N , where K = C 1 µ min , µ max , M .
Proof. The existence of U, V and T follow from the previous theorem. Furthermore, we have

U(t) -V(t) = 1 s=0 Φ t (V 0 + s (U 0 -V 0 )) • (U 0 -V 0 ) .
The result follows from Proposition 2.3.

A SPLITTING SCHEME

In this section, we split the Saut-Xu system (1.1) and we give some estimates for the sub-problems. We consider, separately, the transport part

(3.1) ∂ t ζ + ε √ µ 2 H 2 µ + 1 v ∂ x ζ = 0 ∂ t v + 3ε √ µ 2 v∂ x v = 0, and the dispersive part (3.2) ∂ t ζ -H µ v +ε √ µ 1 2 H µ (v∂ x H µ ζ)+H µ (ζ∂ x H µ v)+ζ∂ x v -1 2 ∂ x ζH 2 µ v = β √ µ∂ x (B µ v) ∂ t v +∂ x ζ - ε √ µ 2 ∂ x ζH µ ∂ x ζ - ε √ µ 2 vH 2 µ ∂ x v = 0.
We denote by Φ t A the flow of System (3.1) and by Φ t D the flow of System (3.2). Remark 3.1. Notice that we keep the term ζ∂ x v in the first equation and we decompose v∂ x ζ as

v∂ x ζ = ∂ x ζ H 2 µ + 1 v -∂ x ζH 2 µ v .
This will be useful for the local wellposedness of the dispersive part.

In the following, we prove the local existence on large time for Systems (3.1) and (3.2).

3.1. The transport equation. The system (3.1) is a transport equation. Then, it is easy to get the following result.

Proposition 3.2. Let s 1 ≥ 0, s 2 > 3
2 and M > 0. We assume that ε, µ satisfies Condition (1.3). Then, there exists a time

T 1 = T 1 (M, µ max ) > 0, such that if |ζ 0 | H s 1 + |v 0 | H s 2 ≤ M, we have a unique solution (ζ, v) ∈ C 0, T 1 ε , H s 1 (R) × H s 2 (R) , to System (3.1)
with initial data (ζ 0 , v 0 ). Furthermore, we have, for all t ≤ T 1 ε ,

(3.3) |ζ(t, •)| H s 1 + |v(t, •)| H s 2 ≤ C(M, µ max ). Finally, if s 2 ≥ 4 and M 1 = max 0≤t≤ T 1 ε |v(t, •)| H s 2 -2
, then for all t ≤ T 1 ε , we have

(3.4) |ζ(t, •)| H s 1 + |v(t, •)| H s 2 ≤ e εC 1 t |U 0 | H s 1 ×H s 2 ,
where C 1 > 0 depends on M 1 and µ max .

Proof. The proof follows from the fact that the quasilinear system (3.1) is symmetric. Thanks to the Coifman-Meyer estimate (see Proposition A.3), we get

d dt |ζ(t, •)| 2 H s 1 + |v(t, •)| 2 H s 2 ≤ Cε √ µ |ζ(t, •)| 2 H s 1 + |v(t, •)| 2 H s 2 3 2 .
Then, we see that the energy is bounded uniformly with respect to ε and µ and we get Estimate (3.3).

For the second estimate, using the same trick that in Lemma 3.1 in [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF], we notice that, if

s 1 ≥ 4, d dt |ζ(t, •)| 2 H s 1 + |v(t, •)| 2 H s 2 ≤ ε √ µ |v(t, •)| H s 2 -2 |ζ(t, •)| 2 H s 1 + |v(t, •)| 2 H s 2 .
By applying the Gronwall lemma, we get the result.

3.2. The dispersive equation. The system (3.2) contains all the dispersive terms of the Saut-Xu system. We have the following estimate for the flow. Proposition 3.3. Let N ≥ 2, and b ∈ L ∞ (R). We assume that ε, β, µ satisfy Condition (1.3). Then, there exists a time T 2 = T 2 M, 1 µ min , µ max such that if

|ζ 0 | H N + 1 2 + |v 0 | H N + |b| L ∞ ≤ M, we have a unique solution (ζ, v) ∈ C 0, T 2 δ , H N + 1 2 (R) × H N (R) to the system (3.
2) with initial data (ζ 0 , v 0 ). Furthermore, we have, for all t ≤ T 2 δ ,

(3.5) |ζ(t, •)| H N + 1 2 + |v(t, •)| H N ≤ C M, µ max , 1 µ min .
Finally, if N ≥ 7, and

M 1 = max 0≤t≤ T 2 δ |ζ(t, •)| H N + 1 2 -2 + |v(t, •)| H N -2 ,
then for all t ≤ T 2 δ , we have

(3.6) |ζ(t, •)| H N + 1 2 + |v(t, •)| H N ≤ e δC 2 t |U 0 | H N +1/2 ×H N ,
where C 2 is a positive constant which depends on µ max , 1 µ min , M 1 . Proof. The proof is an adaptation of the proof of Theorem 2.2 and part IV in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF]. We notice that, in the proof of Saut and Xu, the transport part can be treated separately and does not influence the control of the other terms. Hence, we can use the same symmetrizer S that in Theorem 2.2 (see (2.6)) and we get

d dt E N (ζ, v) ≤ C 1 µ min ε ν E N (ζ, v) 3 2 + β ν E N (ζ, v) .
Then, by Remark 2.1, we get Estimate (3.5). Furthermore, we notice that we use the same trick as in Lemma 3.1 in [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF]. By keeping the same notations as in Theorem 2.2, we get from Equations (2.7) and (2.8) that

d dt E N (U) ≤ δC 1 µ min , µ max GN H 1/2 ×L 2 + E 2 (U) + |b| ∞ E N (U) where GN = ( GN 1 , GN 2 ) t with GN 1 = -∂ N x [H µ , ζ]H µ ∂ x v + ζ(H 2 µ + 1)∂ x v - 1 2 1≤γ≤N -1 C γ N H µ (∂ γ x vH µ ∂ 1+N -γ x ζ) - 1 2 ∂ x ζ(H 2 µ + 1)∂ N x v GN 2 = 1 2 1≤γ≤N -1 C γ N ∂ 1+γ x ζH µ ∂ 1+N -γ x ζ + 1 2 1≤γ≤N C γ N ∂ γ x vH 2 µ ∂ 1+N -γ x v.
To explain how we can adapt the trick used in Lemma 3.1 in [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF], we focus our attention to one term. For

1 ≤ γ ≤ N -1, we have to control ∂ 1+γ x ζH µ ∂ 1+N -γ x ζ L 2 . If γ ≤ N 2 ,
we get from Propositions A.4 and A.1 that

∂ 1+γ x ζH µ ∂ 1+N -γ x ζ L 2 ≤ ∂ 1+γ x ζ H 1 H µ ∂ 1+N -γ x ζ L 2 ≤ C (µ max ) |ζ| H 2+ N 2 |ζ| H N , whereas if γ > N 2 , we have ∂ 1+γ x ζH µ ∂ 1+N -γ x ζ L 2 ≤ ∂ 1+γ x ζ L 2 H µ ∂ 1+N -γ x ζ H 1 ≤ C (µ max ) |ζ| H N |ζ| H 2+ N 2 .
We can mimic this method to control the other terms of GN and, thanks to Propositions A.1, A.2 and A.4, we obtain if N ≥ 7 that

d dt E N (ζ, v) ≤ δC 1 µ min , µ max |ζ(t, •)| H N + 1 2 -2 + |v(t, •)| H N -2 + |b| L ∞ E N (ζ, v) .
Then, Estimate (3.6) follows.

Remark 3.4. Under the assumption of Proposition 3.3 and if N ≥ 7, we get from relations (3.4) and (3.6) that, there exists a time T 3 > 0, such that for all t ∈ 0, T 3 δ ,

|Y t U 0 | H N +1/2 ×H N ≤ e C 3 δt |U 0 | H N +1/2 ×H N ,
where

C 3 = C |U 0 | H N +1/2-2 ×H N -2 , µ max , 1 µ min and T 3 = C |U 0 | H N +1/2 ×H N , µ max , 1 µ min .

ERROR ESTIMATES

The goal of this part is to prove the main result of this paper (Theorem 4.6). Our analysis is based on energy estimates.

4.1.

The local error estimate. The local error is the following quantity

(4.1) e (t, U 0 ) = Φ t U 0 -Y t U 0 .
Our approach is similar to the one developed in [START_REF] Chartier | Uniformly accurate time-splitting methods for the semiclassical schrödinger equation part 2 : numerical analysis of the linear case[END_REF]. We use the fact that Φ t U 0 satisfies a symmetrizable system. Therefore, e satisfies this system up to a remainder and then, we can control e thanks to energy estimates. In the following we give different technical lemmas in order to control the local error. We recall that the transport operator is the operator

A A (ζ, v) = - ε √ µ 2 H 2 µ + 1 v ∂ x ζ 3v∂ x v .
The following proposition gives an estimate of the differential of the transport operator.

A (ζ, v).(η, w) H s 1 ×H s 2 ≤ εC(µ max ) |(ζ, v)| H s 1 +1 ×H s 2 +1 |(η, w)| H s 1 +1 ×H s 2 +1 .
Proof. We have

A (ζ, v).(η, w) = - ε √ µ 2 H 2 µ + 1 v ∂ x η + H 2 µ + 1 w ∂ x ζ 3w∂ x v + 3v∂ x w ,
and the estimate follows from Propositions A.1 and A.4.

We can do the same for the dispersive part (using also Proposition A.1). We recall that the dispersive operator is the operator

D D(ζ, v) = H µ v + ε √ µ 1 2 H µ (v∂ x H µ ζ) + H µ (ζ∂ x H µ v) + ζ∂ x v -1 2 ∂ x ζH 2 µ v -β √ µ∂ x (B µ v) -∂ x ζ + ε √ µ 2 ∂ x ζH µ ∂ x ζ + ε √ µ 2 vH 2 µ ∂ x v Lemma 4.2. Let s > 0, ε, β, µ satisfying Condition (1.3) and b ∈ L ∞ (R). Then, D (ζ, v).(η, w) H s ×H s ≤ C(µ max ) 1+β |b| L ∞+ε |(ζ, v)| H s+1 ×H s+1 |(η, w)| H s+1 ×H s+1 .
Furthermore, we have to control the derivative of the flow Φ t A with respect to the initial data. We denote it by

Φ t A . Lemma 4.3. Let s 1 , s 2 ≥ 0, M > 0, ε, β, µ satisfying Condition (1.3) and b ∈ L ∞ (R). Let (ζ 0 , v 0 ) ∈ H s 1 +1 × H s 2 +1 (R d ) such that, |(ζ 0 , v 0 )| H s 1 +1 ×H s 2 +1 ≤ M.
Then, there exists a time T = T (M, µ max ), such that Φ t A (ζ 0 , v 0 )•(η 0 , w 0 ) exists for all t ∈ 0, T δ and if we denote

η w = Φ t A (ζ 0 , v 0 ) • (η 0 , w 0 ) , for all 0 ≤ t ≤ T δ , |(η, w) (t, •)| H s 1 ×H s 2 ≤ |(η 0 , w 0 )| H s 1 ×H s 2 C (µ max , M ) .
Proof. The quantity (η, w) satisfies the following linear system

∂ t η + ε √ µ 2 H 2 µ + 1 v∂ x η + ε √ µ 2 H 2 µ + 1 w∂ x ζ = 0, ∂ t w + 3ε √ µ 2 v∂ x w + 3ε √ µ 2 w∂ x v = 0, where (ζ, v) = Φ t A (ζ 0 , v 0 ).
The result follows from energy estimates, the Gronwall lemma and Proposition 3.2.

In the following, we use the fact Φ t A • Φ t D satisfies the Saut-Xu system (1.1) up to a remainder. The following lemma is the key point for the control of this remainder.

Lemma 4.4. Let N ≥ 2, M > 0, ε, β, µ satisfying Condition (1.3) and b ∈ L ∞ (R). Let U = (ζ, v) ∈ H N + 1 2 × H N (R d ) such that, |b| L ∞ + |U| H N + 1 2 ×H N (R) ≤ M.
Then, there exists a time T = T M, µ max , 1 µ min > 0, such that Φ t A (U) exists for all 0 ≤ t ≤ T δ , and furthermore,

Φ t A (U) • D (U) -D Φ t A (U) H N -2 ×H N -2 ≤ εC M, µ max , 1 µ min t.
Proof. The existence of T follows from Proposition 3.2. Then, we notice that

Φ t A (U)• D(U) -D Φ t A (U) = t 0 A (Φ s A (U)) • (Φ s A ) (U)•D(U) -D (Φ s A (U)) • A (Φ s A (U)).
Using Lemmas 4.1, 4.2 and Proposition 3.2, we get,

Φ t A (U)• D(U) -D Φ t A (U) H N -2 ×H N -2 ≤ C (µ max , M ) t 0 ε (Φ s A ) (U)• D (U) H N -1 ×H N -1 + |A (Φ s A (U))| H N -1 ×H N -1 .
Then, using Lemma 4.3, the product estimate A.4 and the expression of A, we obtain

Φ t A (U)•D(U)-D Φ t A (U) H N -2 ×H N -2 ≤ εC(µ max , M ) t 0 |D (U)| H N -1 ×H N -1 + |Φ s A (U)| 2 H N ×H N.
Finally, the result follows from the expression of D, the product estimate A.4 and Proposition A.1.

We can now give the main result of this part, the local error estimate.

Proposition 4.5. Let N ≥ 4, M > 0, ε, β, µ satisfying Condition (1.3) and b ∈ L ∞ (R). Let U 0 = (ζ 0 , v 0 ) such that, |b| L ∞ + |U 0 | H N + 1
2 ×H N ≤ M. Then, there exists a time T 4 = T 4 M, 1 µ min , µ max > 0, such that the local error e (t, U) defined in (4.1) exists for all 0 ≤ t ≤ T 4 δ , and furthermore,

|e (t, U 0 )| H N -4+ 1 2 ×H N -4 ≤ δC 4 t 2 , where C 4 = C 1 µ min , µ max , M .
Proof. From Propositions 3.2 and 3.3, we obtain the existence of T . We denote

U(t) = ζ(t) v(t) = Φ t (U 0 ) and V(t) = η(t) w(t) = Φ t A Φ t D (U 0 ) .
Then, from Theorem 2.2 and Propositions 3.2 and 3.3 we also have, for all

0 ≤ t ≤ T δ , (4.2) |U(t, •)| H N + 1 2 ×H N + |V(t, •)| H N + 1 2 ×H N ≤ C 1 µ min , µ max , M .
We know that (ζ, v) satisfy the Saut-Xu system (1.1). Furthermore, (η, w) also satisfy the Saut-Xu system (1.1) up to a remainder

∂ t η w = A (η, w) + D (η, w) + R(t),
where

R(t) = Φ t A Φ t D (U 0 ) • D Φ t D (U 0 ) -D Φ t A Φ t D (U 0 )
. Therefore, the local error e satisfies the following system

(4.3) ∂ t e = 0 H µ -∂ x 0 e + 0 β √ µB µ 0 0 e + T µ ((ζ, v) , (η, w)) -R(t),
where the operator T µ (U, V) is quadratic and satisfies the following estimate, for 0 ≤ s ≤ N -1,

(4.4) |T µ ((ζ, v) , (η, w))| H s ×H s ≤ εC 1 µ min , µ max , M |e| H s+1 ×H s+1 .
Then, since e |t=0 = 0,

e(t, •) = t 0 ∂ t e(s, •)ds,
and since e satisfies (4.3), we obtain, thanks to Estimates (4.2), (4.4) and Lemma 4.4,

(4.5) |e (t, •)| H N -2 ×H N -2 ≤ C 1 µ min , µ max , M t.
Furthermore, we recall that the Saut-Xu system (1.1) is symmetrizable thanks to the symmetrizer (see Theorem 2.2)

S = D tanh( √ µD) 0 0 1 .
Therefore, applying S to the system (4.3), and using the fact that (S•, •) is a norm equivalent to the H (S∂ α x e, ∂ α x e). Then, we get

F(e)(t) ≤ δC 1 µ min , µ max , M t 0 F(e)(s) + s F(e)(s)ds. Denoting M(t) = max [0,t]
F(e)(t), we have

M(t) ≤ δC 1 µ min , µ max , M t 0 M(s) + sds,
and the result follows from the Grönwall's lemma.

Global error estimate.

In this part, we prove our main result. We denote by

U k = Y ∆t k U 0
the approximate solution and by U(t k ) := Φ k∆t U 0 the exact solution at the time t k = k∆t.

Theorem 4.6.

Let N ≥ 7, M > 0, ε, β, µ satisfying Condition (1.3) and b ∈ L ∞ (R). Let U 0 = (ζ 0 , v 0 ) such that, |b| L ∞ + |U 0 | H N + 1 2 ×H N ≤ M.
Let U 0 (t, •) the solution of the Saut-Xu equations (1.1) with initial data U 0 defined on 0, T δ . Then, there exist constants γ, ν, ∆t 0 , C 0 > 0 such that for all ∆t ∈]0, ∆t 0 ] and for all n ∈ N such that 0 ≤ n∆t ≤ T δ ,

|U n | H N + 1 2 ×H N ≤ ν and Φ n∆t (U 0 ) -Y ∆t n (U 0 ) H N -4+ 1 2 ×H N -4 ≤ γ∆t. Proof.
The proof is based on a Lady's Windermere's fan argument and is similar to the one in [START_REF] Carles | On Fourier time-splitting methods for nonlinear Schrödinger equations in the semiclassical limit[END_REF] (see also [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF]). In order to simplify the notations, we forget the dependence on 1 µ min and µ max in all the We denote by X N the following space

X N = H N + 1
2 × H N . Thanks to Theorem 2.2, there exists ρ such that, for all t k = k∆t ∈ 0, T δ ,

|U(t k )| X N = Φ t k U 0 X N ≤ ρ.
We prove by induction that there exists ∆t 0 , γ, ν such that if 0 < ∆t ≤ ∆t 0 , for all n ∈ N with n∆t ≤ T δ ,

(i) |U n -U(t n )| X N -4 ≤ γ∆t, (ii) |U n | X N ≤ e C 3 (M 1 )δn∆t |U 0 | X N ≤ M 0 , (iii) |U n | X N -2 ≤ M 1 , (iv) |U n | X N -4 ≤ 2ρ,
with

M 1 = e C 3 (2ρ)T M, M 0 = e C 3 (M 1 )T M, γ = T max(K, 1)C 4 e C 0 (M 0 )T M 0 , ∆t 0 = min T, T 0 (M 0 ), T 3 (M 0 ), T 4 (M 0 ), ρ γ , K = K M 0 e T C 0 (M 0 )
where C 0 , T 0 , C 3 , T 3 , C 4 , T 4 , K are constants from Theorem 2.2, Remark 3.4, Proposition 4.5 and Inequality (2.11). The above properties are satisfied for n = 0. Let n ≥ 1, and suppose that the induction assumptions are true for 0 ≤ k ≤ n -1. First, we have the following telescopic series (see [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF] or [START_REF] Carles | On Fourier time-splitting methods for nonlinear Schrödinger equations in the semiclassical limit[END_REF])

(4.6) U n -U(t n ) = 0≤k≤n-1 Φ (n-k-1)∆t Y ∆t U k -Φ (n-k-1)∆t Φ ∆t U k .
For k ≤ n -2, since Y ∆t U k = U k+1 , using the induction assumption (ii), we have

|Y ∆t U k | X N -3 ≤ M 0 ,
and from Theorem 2.2, we get

|Φ ∆t U k | X N -3 ≤ e C 0 (M 0 )T M 0 .
Therefore, from Proposition 2.4 and up to replacing K = K M 0 e T C 0 (M 0 ) with max(K, 1), we obtain, for k ≤ n -1 and n∆t ≤ T δ ,

Φ (n-k-1)∆t Y ∆t U k -Φ (n-k-1)∆t Φ ∆t U k X N -4 ≤ K Y ∆t U k -Φ ∆t U(t k ) X N -4 .
Then, using Proposition 4.5 and Inequality (ii), we infer 2 . Therefore, using the telescopic series (4.6), we get

Φ (n-k-1)∆t Y ∆t U k -Φ (n-k-1)∆t Φ ∆t U k X N -4 ≤ δC 4 e C 0 (M 0 )T M 0 K(∆t)
|U n -U(t n )| X N -4 ≤ nC 4 e C 0 (M 0 )T M 0 Kδ(∆t) 2 ≤ C 4 e C (M 0 )T M 0 KT ∆t.
For Estimate (ii), using Remark 3.4 and the induction assumptions (iii) and (ii), we have,

|U n | X N = Y ∆t (U n-1 ) X N ≤ e δC 3 (M 1 )∆t |U n-1 | X N ≤ e C 3 (M 1 )δn∆t |U 0 | X N ≤ M 0 .
We get Estimate (iii) in the same way, using the induction assumptions (iv) and (iii). Finally, for Estimate (iv), using (i), we have

|U n | X N -4 ≤ |U n -U(t n )| X N -4 + |U(t n )| X N -4 ≤ γ∆t + ρ ≤ 2ρ.

NUMERICAL EXPERIMENTS

The aim of this section is to numerically verify the Lie method convergence rate in O(∆t) for the Saut-Xu system (1.1) and to illustrate some physical phenomena.

In other works and particularly on the whole water waves problem (see for example [START_REF] Craig | Numerical simulation of gravity waves[END_REF], [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF], [START_REF] Nicholls | Stability of high-order perturbative methods for the computation of Dirichlet-Neumann operators[END_REF] and references therein), several authors use a discrete Fourier transform even for the transport part. They observe spurious oscillations in the wave profile that lead to instabilities. These errors seem to appear when they evaluate the nonlinear part via Fourier transform because additional terms appear in the approximation, this is the aliasing phenomenon. To fix this problem, they apply at every time step a low-pass filter. The main interest of our scheme is that we do not need one because we use a finite difference method to approximate the nonlinear part.

For the dispersive equation (3.2), we use the forward Euler discretization in time and for the spatial discretization we consider the Fast Fourier Transform (FFT) implemented in Matlab. In this scheme, the interval [0, 1] is discretized by N equidistant points, with spacing ∆x = 1/N . The spatial grid points are then given by x j = j/N , j = 0, ..., N . Therefore, if u j (t) denotes the approximate solution to u(t, x j ), the discrete Fourier transform of the sequence {u j } N -1 j=0 is defined by

û(k) = F d k (u j ) = N -1 j=0
u j e -2iπjk/N , for k = 0, • • • , N -1, and the inverse discrete Fourier transform is given by

u j = F -d j (û k ) = 1 N N -1 k=0 ûk e 2iπkx j , for j = 0, • • • , N -1.
Here F d denotes the discrete Fourier transform and F -d its inverse.

Then, in what follows the numerical scheme to solve (3.2) is given by (5.1)

ζ n+1 j v n+1 j = ζ n j v n j -∆t F n j + S n j G n j where S n j = β √ µF -d j (ikF d k (B µ v n j )
) and F n j = I 1 + I 2 with

I 1 = F -d j i tanh( √ µk) -1 + ε √ µ 2 F d k (v n j F -d j (k tanh( µk) ζn k )) + ε √ µF d k (ζ n j F -d j (k tanh( µk)v n k )) I 2 = ζ n j F -d j (ikv n k ) + 1 2 F -d ik ζn k F -d j tanh( √ µk) 2 vn k .
To approximate the equation (3.1), we use the following finite difference scheme (5.2)

ζ n+1 j v n+1 j = ζ n j v n j -∆t ε √ µ 2 G n 1 3G n 2
where

G 1 = w n j ζ n j+1 -ζ n j-1 2∆x - ∆t 2∆x 2 (w n j ) 2 ζ n j-1 -2ζ n j + ζ n j+1 with w n j = -F -d j tanh( √ µk) 2 vn k + v n j and G n 2 = (v n j+1 ) 2 -(v n j-1 ) 2 2∆x - ∆t 2∆x 2 v n j+1/2 (v n j+1 ) 2 -(v n j ) 2 -v n j-1/2 (v n j ) 2 -(v n j-1 ) 2 with v n j±1/2 = v n j +v n j±1 2 .
We remarked that for our numerical simulations, it is not necessary to decompose the term v∂ x ζ (see Remark 3.1) to get the numerical convergence. Indeed, it seems that since the time step is chosen very small, we obtain a solution of the dispersive equation for each iteration. In this case, we do not need to evaluate the term ∂ x ζH 2 µ v. To ensure the validity of our numerical simulations, we have to be careful of the numerical instability, that why the time and the space steps are chosen in a way that the following CFL condition is satisfied:

(5.3) |v| ∆t ∆x < 1.
5.1. Example 1: Convergence curve. In this example, we consider the following initial data:

ζ 0 (x) = sech √ 3 2 x , v 0 = ζ 0 .
with two different bathymetries: a bump ( b(x) = cos(x) ) and a ripple bottom

b(x) = 0.5 -1 18 (x -8) 2 if 5 ≤ x ≤ 11 0 otherwise .
Note that in order to avoid numerical reflections due to the boundaries and justify of the use of the Fast Fourier Transform, we decide to take rapidly decreasing initial data. Figures 1 and2 display the evolution for different times of the free surface ζ for these two test cases. We decided to take ε = 0.1, µ = 1, β = 1 2 , N = 2 8 , ∆x = 2L/N, T = 10 where N is the mesh modes number, L = 30 the length of the domain and T the final time. Note that the time step ∆t is chosen iteratively in a way that the CFL condition (5.3) is satisfied. Figures 3 displays the convergence curve for this example. We plot the logarithm of the error (in norm H 1 × L 2 ) in function of the logarithm of the time step ∆t. The convergence numerical order is then given by the slope of this curve. For reference, a small line (the dashed line) of slope one is added in this figure. We see that the numerical rate of convergence is greater than 1.

Example 2: Non smooth topographies.

In this example we study the evolution of water waves over a rough bottom. This problem is still a mathematical issue. Many models derived from the Euler equations suppose that the bathymetry is smooth. Even worse, a non smooth bathymetry introduces singular terms in these models. This issue is particularly easy to see for shallow waters models. To handle this, Hamilton ([20]) and Nachbin ([30]) used a coformal mapping to derive long waves models. Notice also the work of Cathala ([10]) who derived alternatives Saint-Venant equations and Boussinesq systems with non smooth topographies which do not involve any singular terms. We notice that the Saut-Xu equations (1.1) can handle a non smooth topography (see Theorem (2.2)) and our numerical scheme too (see Theorem (4.6)). In the following, we give an example with a non smooth bathymetry. We consider the following initial conditions and bathymetry

ζ 0 (x) = v 0 (x) = e -x 2 and b(x) = β 4 (1 + tanh(100(x -2)))(1 -tanh(100(x -8))).
We decided to take µ = 1, β = 0.5, ε = 0.1, N = 2 8 , ∆x = 2L/N where N is the mesh modes number, L = 20 the length of the domain. The time step ∆t is chosen iteratively in a way that the CFL condition (5.3) is satisfied. Figure 4 displays the evolution of the surface ζ for different times over the bottom.

5.3. Example 3: Boussinesq regime. In Section 3, we crucially use the fact that µ is bounded from below. In this example, we test our scheme for small values of µ (also called the shallow water regime). We show that our scheme is still valid even if we do not have a proof of the convergence of our scheme in this regime. In the shallow water regime, there is a huge literature for asymptotic models (see for instance [START_REF] Lannes | The water waves problem[END_REF]). Among all these asymptotic models, we have the KdV equation. It is a model obtained under the Boussinesq regime, i.e. when ε = µ, β = 0 and µ small. In the following, we formally derive a KdV equation from the Saut-Xu equations and we give numerical simulations in this setting.

We recall that, without the assumption ν = 1 µ , the Saut-Xu equations are given by the system (2.3). Notice also that (5.4)

H µ = - √ µ∂ x - 1 3 µ 3 2 ∂ 3 x + O(µ 2 ).
Then if we assume that µ = ε, ν = 1 (since ν ∼ 1 if µ is small) and we drop all the terms of order O(µ 2 ) in System (2.3), we obtain the following equations (5.5)

     ∂ t ζ + ∂ x v + µv∂ x ζ - 1 2 µ 3 2 v∂ x ζ + 1 3 µ∂ 3 x v + µζ∂ x v = 0, ∂ t v + ∂ x ζ + µv∂ x v + µ 3 2 1 2 v∂ x v = 0.
Formally, the solutions of this system are close to the solutions of (2.3) with an accuracy of order O(µ 2 ). Notice that this system is not a standard Boussinesq system (in the sense of [START_REF] Bona | Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I. Derivation and linear theory[END_REF] or [START_REF] Lannes | The water waves problem[END_REF]) because of our nonlinear change of variables (2.2). Using the approach developed in [START_REF] Schneider | Corrigendum: The long-wave limit for the water wave problem I. The case of zero surface tension [mr1780702[END_REF], [START_REF] Bona | Long wave approximations for water waves[END_REF], [START_REF] Alvarez-Samaniego | Large time existence for 3D water-waves and asymptotics[END_REF] (see also Part 7.1.1 in [START_REF] Lannes | The water waves problem[END_REF]) we can check that, formally, the following KdV equation is an asymptotic model of the system (5.5)

(5.6)

∂ t f + 3 2 f ∂ x f + 1 6 ∂ 3 x f = 0.
This means that if we solve (5.5) with the initial data (f 0 , f 0 ) and (5.6) with the initial datum f 0 , the solution (ζ, v) (t, x) of (5.5) is close to (f, f ) (µt, x -t). Furthermore, if we take f 0 (x) = αsech 2 3 4 αx , the solution f of the KdV equation with this initial datum is the soliton f (t, x) = f 0 (x -ct) with c = α 2 . Hence, in this case, the solution of (5.5) and (2.3) are close to a soliton.

In the following we check that the solution to (1.1) is indeed close to the KdV solution when µ is small. We simulate one soliton. We took v 0 (x) = ζ 0 (x) = sech 2 √ 3 2 x , ε = µ = 0.01, α = 1 and the final time is T = 10. We decided to take N = 2 9 , ∆x = 2L/N where N is the mesh modes number, L = 30 the length of the domain. The time step ∆t is chosen iteratively in a way that the CFL condition (5.3) is satisfied. Figure 5 represents the evolution of this soliton at different times. Hence, our sheme is still valid when µ is small. In deep water (µ not small), the KdV approximation ceases to be a good approximation. In order to get some insight on the range of validity of the KdV approximation, we compare in Figure 6 the solution of (1.1) to the exact soliton after a time T = 10 for various values of µ. We took the same numerical parameters that before. We notice that even for ε = µ = 0.1 and a final time T = 1 µ , the KdV approximation remains a good approximation of the Saut-Xu system. In this example we study the evolution of water waves over a rapidly varying periodic bottom. We assume that µ = 1. This problem is linked to the Bragg reflection phenomenon (see for instance [START_REF] Mei | Resonant reflection of surface water waves by periodic sandbars[END_REF], [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF], [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]). We take (5.7)

ζ 0 = v 0 = sech 2 √ 3 2
x and b(x) = cos(αx).

We decided to take N = 2 9 , ∆x = 2L/N where N is the mesh modes number, L = 30 the length of the domain. The time step ∆t is chosen iteratively in a way that the CFL condition (5.3) is satisfied. Figure 7 compares the evolution of water waves when we take α = 10 (blue line) and when we take b(x) = 0 (blue line). Figure 8 displays the difference between the case of a flat bottom and the case of a bottom of the form b(x) = cos(αx) for different values of α. We observe an homogenization effect when α is large. It seems that when α goes to infinity, a solution of the Saut-Xu equations converges to a solution of the Saut-Xu equations with a flat bottom (corresponding to the mean of b). Notice that this result is different from what we could see in the literature ( for instance [START_REF] Chupin | Roughness effect on Neumann boundary condition[END_REF] or [START_REF] Craig | Water waves over a rough bottom in the shallow water regime[END_REF]), since we take a bottom of the form b(x) = cos(αx) and not of the form b(x) = 1 α cos(αx). Our numerical simulations suggest therefore a homogenization effect for large amplitude bottom variations that has not been investigated so far.

Lemma 4 . 1 .

 41 Let s 1 , s 2 ≥ 0 and ε, µ satisfying Condition (1.3). Then,

1 2 ×

 2 L 2 -norm, we obtain, thanks to estimates (4.2), (4.4) and (4.5) and Lemma 4.4, d dt F(e) ≤ C 1 µ min , µ max , M βF(e) + εt F(e) , where F(e) = |α|≤N -4

FIGURE 1 .

 1 FIGURE 1. Upper: Evolution of the free surface for different times t=2.5,5,7.5 and 10. Lower: bottom topography and initial condition.

FIGURE 2 .

 2 FIGURE 2. Upper: Evolution of the free surface for different times t=2.5,5,7,5 and 10 . Lower: bottom topography and initial condition.

FIGURE 3 .

 3 FIGURE 3. Convergence curve (for the H 1 × L 2 -norm) for the Lie method for two bottoms : bump (black line) and ripple bottom (blue line) for T = 10.

FIGURE 4 .

 4 FIGURE 4. Evolution of the free surface (blue lines) for different times t=0,3,6,9 and 12 over a rough bottom (black line).

FIGURE 5 .

 5 FIGURE 5. Evolution of the soliton at different times t = 0, 3, 6, 9 (ε = 0.01).

FIGURE 6 .

 6 FIGURE 6. Difference after a time T = 10 between a real soliton and a soliton generated by our scheme with the same initial data and for different values of ε = µ. Abscissa : value of ε; Ordinate : quotient of the difference after a final time T = 10 by the maximum of the soliton.

FIGURE 7 .

 7 FIGURE 7. Comparison between the evolution of a water wave (blue line) over a bottom of the form b(x) = cos(10x) (dashed line) and the evolution of a water wave over a flat bottom (black line). ε = 0.05, β = 0.5.

FIGURE 8 .

 8 FIGURE 8. Difference between a water wave over a rapidly varying topography b(x) = cos(αx) and a water wave over a flat bottom. Abscissa : value of α; Ordinate : quotient of the difference after a final time T = 10 by the maximum of ζ 0 .

ACKNOWLEDGMENTS

The second author has been partially funded by the ANR project Dyficolti ANR-13-BS01-0003.

APPENDIX A

In this part, we give some estimate for the operator H µ and some standard product and commutator estimates. For the estimates for H µ , we refer to part III in [START_REF] Saut | Well-posedness on large time for a modified full dispersion system of surface waves[END_REF]. For the other estimates we refer to [START_REF] Alinhac | Opérateurs pseudo-différentiels et théorème de Nash-Moser[END_REF] and [START_REF] Lannes | Sharp estimates for pseudo-differential operators with symbols of limited smoothness and commutators[END_REF]. We recall that H µ is defined by

First, we show that H µ is a zero-order operator.

Proposition A.1. Let s ≥ 0 and µ satisfying Condition (1.3). Then,

Furthermore, for all s ≥ r ≥ 0,

Then, we give a commutator estimate for H µ .

Proposition A.2. Let s ≥ 0, t 0 > 1 2 , r ≥ 0, and µ satisfying Condition (1.3). Then,

and

We recall the well-known Coifman-Meyer estimate. We recall also that Λ is the Fourier multiplier

Then we have the following commutator estimate

We recall also the following product estimate.

Proposition A.4. Let s 1 , s 2 , s such that s 1 + s 2 ≥ 0, s ≤ min (s 1 , s 2 ) and s < s 1 + s 2 -1 2 . Let u ∈ H s 1 (R) and v ∈ H s 2 (R). Then,


Olivier Sans

Remi Coletta

Gilles Trombettoni

An Interval Filtering Operator for Upper and Lower Bounding in Constrained Global Optimization

This paper presents a new interval-based operator for continuous constrained global optimization. It is built upon a new filtering operator, named TEC, which constructs a bounded subtree using a Branch and Contract process and returns the parallel-to-axes hull of the leaf domains/boxes. Two extensions of TEC use the information contained in the leaf boxes of the TEC subtree to improve the two bounds of the objective function value: (i) for the lower bounding, a polyhedral hull of the (projected) leaf boxes replaces the parallel-to-axes hull, (ii) for the upper bounding, a good feasible point is searched for inside a leaf box of the TEC subtree, following a look-ahead principle. The algorithm proposed is an auto-adaptive version of TEC that plays with both extensions according to the outcomes of the upper and lower bounding phases. Experimental results show a significant speed-up on several state-of-the-art instances.

I. INTRODUCTION

Interval Branch & Bound algorithms are used to solve constrained global optimization problems 1 in a reliable way [START_REF] Van Hentenryck | Numerica : A Modeling Language for Global Optimization[END_REF], [START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF], [START_REF] Lebbah | An Efficient and Safe Framework for Solving Optimization Problems[END_REF], [START_REF] Ninin | A Reliable Affine Relaxation Method for Global Optimization[END_REF], [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF]. They interleave branching, filtering and bounding steps for calculating an optimal solution and its cost with a bounded error or a proof of infeasibility. The branching step selects a variable and splits its domain into two sub-domains. To limit the combinatorial explosion, the filtering step can contract the domains without loss of solutions. In several interval Branch & Bounds, the objective function value is handled as an additional variable, so that the filtering step is also used for finding a good lower bound, i.e. a point such that no other feasible point exists with a cost lower than it. The lower bounding step is used to prove the optimality of the solution.

In this paper, we first introduce a new filtering operator, called TEC (Tree for Enforcing Consistency), for contracting the domains and improving the lower bound of the objective function value. TEC is a generalization to several dimensions of the constructive disjunction used by the operator CID [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF]. It constructs a subtree, bounded in size by a parameter, by using a Branch & Contract process, and returns a parallelto-axes hull of the leaf boxes. Two extensions of TEC exploit the information contained in the leaf boxes of the TEC subtree. In a first extension, the basic parallel-to-axes hull is 1 In this paper, we consider minimization, without loss of generality.

replaced by a polyhedral hull of the (projected) leaf boxes for improving the filtering capacity of TEC and the lower bounding. This polyhedral hull is achieved by a rigorous variant of the Graham computational geometry algorithm. A second extension of TEC focuses on the upper bounding step of global optimization, which searches for a feasible point with a cost better than the cost of the best solution found so far. This extension searches for a good feasible point inside a selected leaf box of the TEC subtree, thus following a look-ahead principle.

Finally, we introduce an advanced version of TEC endowed with both extensions. This LGT (Lazy Graham TEC) algorithm calls the standard hull or the Graham polyhedral hull algorithm in an adaptive way during the optimization process. When the upper bound has not been improved for long, LGT intensifies its effort for proving optimality: it tries to improve the lower bound of the objective function value by switching to the polyhedral hull of the TEC subtree during the next nodes. If a new upper bound is found,

LGT switches back to the less time-consuming box hull, preferring investing in the upper bounding phase.

After presenting the background in Section II, TEC is described in Section III, and LGT is detailed in Section IV. Section V carries out an empirical study of these operators implemented in the Ibex interval library [START_REF] Chabert | Contractor Programming[END_REF], and compares them to state-of-the-art operators on constrained optimization problems belonging to the Coconut benchmark [START_REF] Shcherbina | Benchmarking Global Optimization and Constraint Satisfaction Codes[END_REF]. Finally, Section VI presents the related work.

II. BACKGROUND An interval [x

i ] = [x i , x i ] defines the set of reals x i s.t. x i ≤ x i ≤ x i , where x i and x i are floating-point numbers. IR denotes the set of all intervals. A box [X] is a cartesian product of intervals [x 1 ] × ... × [x i ] × ... × [x n ]. ω([x i ]) denotes the width x i -x i of an interval [x i ].
The width of a box is given by the width of its largest dimension (i.e., max i=1..n ω([x i ])). Since the union of several boxes is generally not a box, the box hull operator returns the minimal box enclosing all of them.

Interval arithmetic [START_REF] Moore | Interval Analysis[END_REF] has been defined to extend to IR elementary functions over R. For instance, the interval sum is defined by

[x 1 ]+[x 2 ] = [x 1 +x 2 , x 1 +x 2 ]. When a function
f is a composition of elementary functions, an extension of f to intervals must be defined to ensure a conservative image computation.

Definition 1: Extension of a function f : R n → R to IR.

[f ] : IR n → IR is said to be an extension of f to intervals if:

∀[X] ∈ IR n [f ]([X]) ⊇ {f (X), X ∈ [X]} ∀X ∈ R n f (X) = [f ](X)
In our context, the expression of a function f is always a composition of elementary functions (e.g., sin, cos, +, ×, exp, log, ...). The natural extension [f ] N is then simply a composition of the corresponding interval operators.

This paper deals with continuous global optimization under inequality constraints defined by:

min X∈[X] f (X) s.t. g(X) ≤ 0 (1)
where f : R n → R is the (non convex) real-valued objective function and g : R n → R m is a (non convex) vector-valued function. X = (x 1 , ..., x i , ..., x n ) is a vector of variables varying in a box [X]. A point X is said to be feasible if it satisfies the constraints. An interval Branch & Bound (B&B) scheme for continuous constrained global optimization is described below. Algorithms 1 and 2 correspond to a simplified version of the algorithms implemented in IbexOpt [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF].

The algorithm is launched with the vector of constraints g, the objective function f and with the input domain initializing a list Boxes of boxes to be handled. obj is the absolute or relative precision required on the objective function value and is used in the stopping criterion. The algorithm outputs a feasible floating point X f of cost f and a proof that no feasible points can be found with a cost lower than f min .

Let X f * and f * be the optimal real point and its cost. The algorithm returns

f min ≤ f * ≤ f with an absolute ( f -f min ≤ obj ) or relative ( f -fmin | f | ≤ obj ) precision.
The variable x obj (corresponding to the objective function value) is added to the problem (to the vector X of variables) with the constraint f (X) = x obj . The algorithm maintains two main types of information during the iterations:

• f : the value of the best feasible point X f found so far, • f min : the minimal value of the lower bounds x obj of the boxes to explore. The lower bound f min guarantees that, in every box, no feasible point exists with an objective function value lower than it. The procedure SelectBox selects the next node to handle. The selected box [X] is then split into two subboxes [X] 1 and [X] 2 along one dimension (selected by a branching strategy). Both sub-boxes are then handled by the Contract&Bound procedure (see Algorithm 2).

A constraint x obj ≤ fobj is first added to the problem for decreasing the upper bound of the objective function in Algorithm 1: IntervalBranch&Bound (f, g, X, box, obj , sol )

1 f min ← -∞; f ← +∞ 2 Boxes ← {box} 3 while Boxes = ∅ and f -f min > obj and f -f min | f | > obj do 4 [X] ← SelectBox (Boxes, criterion) 5 Boxes ← Boxes\{[X]} 6 ([X] 1 , [X] 2 ) ← Bisect ([X]) 7 ([X] 1 , X f , f , Boxes) ← Contract&Bound ([X] 1 , f, g, X, obj , X f , f , Boxes) 8 ([X] 2 , X f , f , Boxes) ← Contract&Bound ([X] 2 , f, g, X, obj , X f , f , Boxes) 9 PushIfNonEmpty([X] 1 , Boxes) 10 PushIfNonEmpty([X] 2 , Boxes) 11 f min ← min [X]∈Boxes x obj
the box. Algorithm 2 then contracts the handled box without loss of feasible part [START_REF] Chabert | Contractor Programming[END_REF]. In other words, some infeasible parts at the bounds of the domain are discarded using filtering operators. First, CPContract uses a filtering operator which comes from constraint programming [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF], [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF], [START_REF] Van Hentenryck | Numerica : A Modeling Language for Global Optimization[END_REF]. Then, ConvexContract uses a convexification algorithm [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF] which relies on a polyhedral convex relaxation of the system to contract the box. These contractions work on the extended box including the objective function variable x obj and on the associated constraint, so improving x obj amounts to improving the lower bound of the objective function image (lower bounding).

Algorithm 2: Contract&Bound ([X], f, g, X, obj , X f , f , Boxes)

1 g ← g ∪ {x obj ≤ f -obj } 2 [X] ← CPContract ([X], g ∪ {f (X) = x obj }) 3 [X] ← ConvexContract ([X], g ∪ {f (X) = x obj }) 4 if [X] = ∅ then 5 (X f , cost) ← FeasibleSearch ([X], f, g , obj ) 6 if cost < f then 7 f ← cost 8 Boxes ←FilterOpenNodes (Boxes, f -obj ) 9 return ([X], X f , f , Boxes)
The last part of the procedure carries out upper bounding. FeasibleSearch calls one or several heuristics searching for a feasible point X f that improves f , the best cost found so far. If the upper bound of the objective value is improved, the FilterOpenNodes procedure performs a type of garbage collector on all the open nodes by removing from Boxes all the nodes having x obj > fobj . After the calls to Contract&Bound, Algorithm 1 pushes the two sub-boxes in the set Boxes of open nodes.

Several filtering operators have been proposed by the CP community. The HC4 algorithm [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] is a well-known constraint propagation algorithm. The 3B algorithm [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF], like SAC [START_REF] Debruyne | Some Practicable Filtering Techniques for the Constraint Satisfaction Problem[END_REF] for finite CSP, relies on the shaving principle to prune a sub-interval on a bound of an interval. The shaving consists in first restricting an interval to a sub-interval on a bound. If the application of a constraint propagator, like HC4, on the corresponding subproblem leads to a wipeout, the sub-interval is removed from the domain. Otherwise, the sub-interval is kept and the pruning effort performed by the constraint propagator is lost. The constructive disjunction permits to exploit this pruning effort by performing the box hull of all the boxes calculated by the subproblems. CID uses this principle to improve the contraction. First, CID selects a variable x i and splits its interval [x i ] in s cid sub-intervals. Then, each corresponding subproblem is contracted by a constraint propagator like HC4, and the box hull of the contracted boxes is performed. This process is achieved on each variable while a contraction is obtained. 3BCID [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] is an efficient algorithm that hybridizes 3B and CID. This state-of-the-art variant of CID is used in the experiments shown in Section V.

III. TEC

Contrarily to 3B, CID is able to prune values from all variable domains when it performs shaving tests on one of them. TEC generalizes the CID principle. TEC steps up the shaving process by shaving a subset of variables simultaneously and uses constructive disjunction to be able to contract all variable domains. Indeed, TEC constructs a bounded subtree, called here TEC tree, and performs the box hull of all generated leaves (boxes). The TEC tree is constructed in breadth-first search by a Branch & Contract process where each node is handled by constraint propagation. Algorithm 3 describes this process.

A. TEC Algorithm

Algorithm 3 is called at line 2 of Algorithm 2. It is launched with the vector of constraints, and an initial box to contract. The combinatorial process is achieved by the loop at line 4. First, the initial box is contracted by a constraint propagator Φ and is pushed into the list Boxes. At each iteration of the loop, the first box [X] of Boxes is extracted, following a "First In, First Out" principle. If the size of [X] is inferior to a given precision sol , [X] is stored in AtomicBoxes (line 7). Otherwise, [X] is split into two subboxes along one dimension. Both are contracted by Φ and are pushed at the end of Boxes. This process is repeated until the number of contractions reaches the value given by the T ECnodes parameter or until Boxes is empty. Finally, TEC returns the box hull of the leaf boxes of the TEC tree.

B. Implementation Choices for TEC

Several implementation choices have been validated by preliminary experiments. They have been obtained on a subset of the benchmark used in this paper (see Section V).

[X] ← Pop(Boxes) 6 if ω([X]) ≤ sol then 7 Push([X], AtomicBoxes) 8 else 9 ([X] 1 , [X] 2 ) ← Bisect([X]) 10 PushIfNonEmpty(Φ([X] 1 , g , sol ), Boxes) 11 PushIfNonEmpty(Φ([X] 2 , g , sol ), Boxes) 12 #nodes ← #nodes + 2 13 return BoxHull(Boxes ∪ AtomicBoxes)
Variable Selection Heuristic TEC selects a subset of variables to contract the box. Indeed, selecting a good subset of variables is very important to apply a strong contraction. The importance of a variable depends on constraints, but also on the box studied. Consequently, the importance of a variable changes during the search. Therefore we compared different dynamical variable selection heuristics available in Ibex [START_REF] Chabert | Contractor Programming[END_REF]. The tests led us to choose the SmearSumRel2 heuristic. This was not a surprise since SmearSumRel is state-of-the-art both in constraint satisfaction and global optimization. Although it is not always the best, this strategy appears to be more robust than its competitors on the tested benchmark.

Subfiltering Operator

To choose the subfiltering operator Φ, we compared the different filtering operators available in Ibex. Experimental studies showed that using HC4 in TEC provides the best results, due to its low computational time.

TECnodes Value

Experiments showed that the most significant impact on performance is brought by the size of the TEC tree. We compared several values of TECnodes (see Algorithm 3) for each instance of our benchmark. Depending on the instance considered, the optimal value of T ECnodes varied from 10 to 100. Further work could investigate how to automatically adapt TECnodes according to the instance. Lacking such an adaptive strategy, we used 25 as a fixed value for all benchmarks in the following experiments, since we identify it as a robust value.

Overall, the different implementation choices led to a current version of TEC that is comparable in performance to 3BCID.

IV. LGT

Although TEC is comparable in performance to the stateof-the-art filtering operator 3BCID, its tree (i.e., the TEC tree) contains additional information.

LGT uses this additional information in an auto-adaptive way to better improve the two bounds of the objective function value. During the optimization process, LGT can switch between two versions of TEC, namely TEC-UB and Graham-TEC.

For improving the lower bounding and the contraction, Graham-TEC adds new linear constraints in the polytope generated by the polyhedral convexification algorithm, e.g. X-Newton (see the ConvexContract procedure in Algorithm 2). For the upper bounding, TEC-UB makes better use of the In-X-Newton procedure.

Let us first present the convexification algorithms X-Newton and In-X-Newton, available in our optimization code IbexOpt, for supporting the observation that led us to propose TEC-UB (see Section IV-B) and Graham-TEC (see Section IV-C).

A. Brief Description of X-Newton and In-X-Newton X-Newton [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF] is a convexification algorithm which relies on an outer interval linear form of the constraints to contract the box. It can implement the ConvexContract procedure called in Algorithm 2.

A polyhedral convexification is built upon an adaptation of a specific first order interval Taylor form of a nonlinear function [START_REF] Lin | LP Strategy for the Interval-Newton Method in Deterministic Global Optimization[END_REF]. Consider a function g j : R n → R defined on a domain [X]. For any variable x i ∈ X, let

[a i ] = ∂gj ∂xi N ([X]
) be the natural interval evaluation of the partial derivative of g j w.r.t. x i . The idea is to (lower) tighten g j (X), for every point X in [X], using an interval Taylor form g l j (X) expanded at a vertex of the box, e.g. X.

∀X ∈ [X] g l j (X) = g j (X)

+ i a i (x i -x i ) ≤ g j (X) (2)
Thus, considering an inequality constraint g j (X) ≤ 0, Expression (2) defines a hyper-plane g l j (X) ≤ 0 bounding the solution set from below: g l j (X) ≤ g j (X) ≤ 0. Note that, contrarily to the usual case where the Taylor expansion point is taken in the middle of the box, this interval Taylor form is polyhedral because for all x i , (x i -x i ) ≥ 0.

X-Newton defines a hyper-plane for the objective function, and for each constraint to obtain a polytope enclosing the solution set. Then, X-Newton contracts the interval of each variable x i by calling a Simplex algorithm twice, improving potentially both bounds of [x i ]. Improving x obj performs lower bounding.

In-X-Newton [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF] is a heuristic searching for a feasible point (see Algorithm 2) which relies on an inner interval linear form computing an inner region. An inner region r in is a feasible subset of the box [X], i.e., r in ⊂ [X] and all points X ∈ r in satisfy all the constraints g(X) ≤ 0. The inner regions computed by In-X-Newton are polytopes. Symmetrically to Eq. ( 2):

∀X ∈ [X] g l j (X) = g j (X)+ i a i (x i -x i ) ≥ g j (X) (3)
In-X-Newton defines a hyper-plane for each constraint to obtain an inner region. Then, In-X-Newton defines a hyper-plane for the objective function and calls a Simplex algorithm minimizing x obj to find the best solution inside the inner region. This solution is also a solution of the original system. Even if it is generally not the optimal one, it can be used to update the upper bound.

Both algorithms use an outer (resp. inner) interval linear form. In fact, the generated polytope can greatly overestimate (resp. underestimate) the set of solutions. In some cases, the overestimation may produce a polytope larger than the box, and the underestimation may produce an empty inner region even if the box contains solutions. In practice, these cases occur more often in large boxes. This is mainly explained by the first order interval Taylor form. Close to the expansion point (e.g., X in Eq. ( 2)), the approximation is tight, even exact at the expansion point. Thus, the smaller the box is, the sharper will be the approximation of the solution set.

B. TEC-UB

To improve the upper bounding, TEC-UB simply calls In-X-Newton to search for a good feasible point in a small leaf box of the TEC tree (see Fig. 1). Thus, In-X-Newton is called twice:

• once on the contracted hull box (like in the standard interval B&B), and • another time on the selected leaf box of the TEC tree. As explained in Section IV-A, the smaller the box is, the smaller is the underestimation. Therefore, TEC-UB selects the smallest box in the leaf boxes of the TEC tree (i.e., min [X] ω([X])). To break ties, the selection is made using the objective function, i.e. minimizing x obj , for finding a better solution. Difference between TEC and TEC-UB on a system with one variable x and the objective variable x obj . The region removed by contraction appears light greyed out and the set of solutions is dark greyed out (or red). [X] 1 , [X] 2 and [X] 3 are the leaf boxes of the TEC tree. Left: Just one call to In-X-Newton is achieved on the contracted box (i.e., the box hull of the leaf boxes). Right: Two calls to In-X-Newton are performed, the second one on the selected leaf box [X] 3 .

C. Graham-TEC

Another possible way to improve the outer approximation of the solution set is to add hyper-planes to the polytope generated in ConvexContract. Graham-TEC uses the additional information contained in the leaf boxes of the TEC tree to learn implicit binary constraints. In other words, Graham-TEC:

1) constructs the TEC tree in the same way as TEC, 2) computes, for each pair of variables, a polyhedral convex hull of the leaf boxes of the TEC tree (see Algorithm 4), 3) transmits to the ConvexContract operator the polytope made of these binary constraints and the box hull of the TEC tree leaves. Finally, ConvexContract (e.g., X-Newton, Quad [START_REF] Lebbah | Efficient and safe global constraints for handling numerical constraint systems[END_REF]) enriches its polytope with the additional constraints given in input, thus computing generally a sharper contracted box than TEC would do.

Convex Hull of the TEC Tree Leaves

Let us detail hereafter how the polyhedral convex hull of boxes is computed. Algorithm 4 is called at the end of TEC (see Algorithm 3). The algorithm is launched with the set L of the TEC tree leaf boxes and with the box hull [H] of these leaf boxes. γ is a parameter used to decide whether a constraint learned is kept or not. γ measures a gain relative to each interval size of [H].

Algorithm 4: ConvexHull(L, [H], γ)

1 S ← ∅ / * a set of linear constraints * / 2 foreach (x i , x j ) ∈ X 2 with i < j do

3 N E ← ∅; N W ← ∅; SW ← ∅; SE ← ∅ / * quadrant * / 4 foreach [X] ∈ L do 5 N W ← N W ∪ (x i , x j ); N E ← N E ∪ (x i , x j ) 6 SW ← SW ∪ (x i , x j ); SE ← SE ∪ (x i , x j ) 7 foreach quadrant ∈ {N E, N W, SW, SE} do 8 if ¬Cover(quadrant, [H], γ) then 9 F ront ← GrahamFront(quadrant) 10 S ← S ∪ GetConstraints(F ront, x i , x j , quadrant) 11 return S
For each pair of variables (x i , x j ), Algorithm 4 constructs the convex hull of the leaf boxes projected on (x i , x j ) and computes linear constraints from this (2-dimensional) convex hull. The projection of a box on (x i , x j ) is a rectangle, in which each vertex belongs to a quadrant (see Fig. 2). Thus, at lines 5 and 6, the 4 points of each projected leaf box are pushed into 4 lists corresponding to the 4 quadrants. Then, the 2-dimensional convex hull is constructed quadrant by quadrant. Some quadrants are rapidly treated. Let [H] ij be the rectangle projected from [H] onto the variables x i and x j . When a vertex of [H] ij and a point of a projected leaf box in L coincide, no constraint can be found on this (covered) quadrant, as shown in Fig. 2. It is therefore useless to handle the points of this quadrant. Moreover, we consider that a quadrant is also covered if a point of a projected leaf box is too close (on each dimension) to the corresponding vertex of [H] ij (see Fig. 2). The idea is to learn only impacting constraints.

Then, the GrahamFront function, based on the Graham Scan algorithm, is called on each non covered quadrant. Graham Scan is a classical computational geometry algorithm that constructs a convex hull of given 2-dimensional points [START_REF] Graham | An Efficient Algorithm for Determining the Convex Hull of a Finite Planar Set[END_REF]. GrahamFront just computes the front of the convex hull relative to the quadrant studied and returns the corresponding list of points, in counterclockwise order. We call this list of points the front of a quadrant. Let r be the number of leaf boxes, the complexity of GrahamFront is dominated by the sort (counterclockwise order) in O(r.log(r)) [START_REF] Graham | An Efficient Algorithm for Determining the Convex Hull of a Finite Planar Set[END_REF].

Finally, GetConstraints is called to compute the constraints corresponding to the front of a quadrant. Still with the objective of only keeping impacting constraints, when a constraint does not give a significant gain (see Fig. 2), it is removed. The impacting constraints are added to the ConvexContract polytope.

Computation of Rigorous Constraints

The algorithm GetConstraints (see Algorithm 5) generates rigorous constraints from the front of a quadrant. Each pair of consecutive points ((x i , y i ), (x i+1 , y i+1 )) of the F ront list is iteratively handled. The goal is to rigorously approximate the line y = ax + b that passes through the 2 points. Indeed, a = yi+1-yi xi+1-xi and b = y i -(x i × a) are generally not floating-point numbers. The round-off error may generate a constraint (halfspace) which does not include one of both points. To avoid this drawback, a and b are rigorously calculated using interval arithmetic [START_REF] Moore | Interval Analysis[END_REF]. From the two points (x i , y i ) and (x i+1 , y i+1 ), interval arithmetic calculates the intervals [a] and [b] enclosing the respective real values. Then, to generate a rigorous constraint, the adequate bound of [a] and [b] is selected depending on the quadrant (see line 7 to line 13). The reader will be able to easily check the choice of bounds. Fig. 3 shows the principle in the quadrant SE.

Algorithm 5: GetConstraints(F ront, x, y, quadrant) The line 14 adds a condition to guarantee that the previous points satisfy the constraint. Theoretically, the round-off of a may produce this pathological case when the 3 points are quasi-aligned. This case never occurred in our tests.

1 P ← ∅ / * a set of rigorous hyperplanes * / 2 foreach (x i , y i ) ∈ F ront do 3 [a] ← ([y i+1 ] -[y i ])/([x i+1 ] -[x i ]) 4 [b] ← [y i ] -([x i ] × [a])

D. The Lazy-Graham-TEC Auto-adaptive Algorithm

Lazy-Graham-TEC (LGT) is an advanced version of TEC switching between TEC-UB and Graham-TEC in an autoadaptive way during the optimization process. The main parameter of LGT is a number p of nodes in the search tree.

LGT starts by using TEC-UB in order to better improve the upper bounding. When LGT does not find a new upper bound after a given number p of nodes, LGT is betting that the best feasible point in a local minimum (or an optimal solution) has been found. Therefore, the stronger filtering operator Graham-TEC is called to more quickly reach another region of the tree search, or prove optimality of the solution by improving the lower bound.

When a new upper bound is found, LGT switches back to TEC-UB for better intensifying the search.

LGT performs a new phase of at least p nodes with TEC-UB, before switching again to Graham-TEC, and so on.

V. EXPERIMENTS

We have run experiments on constrained global optimization instances. Section V-A reports the implementation and the protocol used. Section V-B compares the different versions of TEC and their settings. Section V-C presents the results obtained by LGT, compared to HC4 and 3BCID.

A. Experimental Protocol and Implementation

All the versions of TEC have been implemented in the Interval-Based EXplorer Ibex [START_REF] Chabert | Contractor Programming[END_REF]. This free C++ library has facilitated the implementation of our filtering operators by providing us a direct access to the interval Branch & Bound IbexOpt [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF] presented in Section II, the X-Newton, In-X-Newton algorithms and filtering operators as HC4 and 3BCID. All the experiments were run on the same computer (Intel X64, 2.6GHz, 32Go RAM).

All the experiments were performed on constrained global optimization instances solved by IbexOpt. All the parameters of IbexOpt have been fixed to a given set of adequate values common to all the tested instances. Thus, SmearSumRel was used as variable selection heuristic, X-Newton as the algorithm used by the procedure ConvexContract, and the algorithms In-HC4 [START_REF] Araya | Upper Bounding in Inner Regions for Global Optimization under Inequality Constraints[END_REF] and In-X-Newton as FeasibleSearch heuristic searching for a good feasible point (upper bounding). The precision settings were fixed to obj = 1.e-8 and sol = obj 100 . A timeout of 1800 seconds was chosen. The parameters used to construct the TEC tree were assigned with the values specified in Section III-B, for all versions of TEC.

Tests have been performed on a benchmark of 82 satisfiable instances of constrained optimization problems belonging to the Coconut benchmark [START_REF] Shcherbina | Benchmarking Global Optimization and Constraint Satisfaction Codes[END_REF]. These instances correspond to all the instances from the series 1 and 2 of the Coconut benchmark that have a reasonable size 3 (between 6 and 50 variables) and are significant (solved in more than 1 second by at least one approach).

Note that the equations h(X) = 0 in these instances have been relaxed by inequalities -1.e-8 ≤ h(X) ≤ +1.e-8 because IbexOpt needs to achieve upper bounding in non empty feasible regions [START_REF] Araya | Upper Bounding in Inner Regions for Global Optimization under Inequality Constraints[END_REF] (see Section IV-A). Note that other interval Branch & Bound solvers like Icos [START_REF] Lebbah | An Efficient and Safe Framework for Solving Optimization Problems[END_REF] and GlobSol [START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF] can provide a guaranteed optimum also on problems having non relaxed equations. A comparison between Ibex and these solvers can be found in [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF].

B. Comparing the Different Versions of TEC

Table I shows a synthetic comparison between TEC and the advanced versions of TEC, i.e. TEC-UB, Graham-TEC and LGT. We report the time and number of nodes gain ratios, comparing each strategy s with the reference strategy TEC. We take into account the CPU time and the number of nodes obtained by a strategy s on an instance i. These experimental results are averaged over 10 trials obtained with different seeds, due to the random choices made by IbexOpt. The following formulas are used to display the time comparison:

• Total Gain(s):

i T ime T EC (i) -i T ime s (i) i T ime T EC (i) (4)
• Max Loss(s):

max i T ime s (i) -T ime T EC (i) T ime s (i) (5)
• Max Gain(s):

max i T ime T EC (i) -T ime s (i) T ime T EC (i) (6)
Same formulas are used for comparing the number of nodes. Graham-TEC shows the lowest gain on the total time and the second best gain on the total number of nodes. This confirms that the convex hull computed by Graham-TEC produces a significant gain of contraction compared to the box hull achieved by TEC, but at a high cost. The improvement of the upper bounding brought by TEC-UB also produces a gain on the total number of nodes, at lower computational time than Graham-TEC. Finally, LGT exploits these two approaches and obtains the best total gain both in time and number of nodes. These results seem to confirm the relevance of the adaptation policy proposed and show a complementarity between the improvement of the upper bounding (TEC-UB) and the improvement of the contraction and the lower bounding (Graham-TEC) leading to a drop in the number of nodes.

To reach a good value of the parameter p used by LGT (number of nodes visited with no improvement of the upper bound), we tested several values. Below p = 500, the results did not show good performances and above p = 1500 too. A robust tested value was 1000 with the γ parameter set to 0.05. Note that first attempts to tune p in an auto-adaptive way during the optimization process were not successful.

C. Results Obtained by LGT

Compared to the HC4 constraint propagation algorithm,

LGT obtains better results (see Fig. 4). Three additional instances (haldmads, ex6_2_10 and disc2) are solved by LGT. When we extend the timeout to 7200 s, only ex6_2_10 is solved in 2530 s by HC4.

Compared to the 3BCID state-of-the-art domain shaving algorithm, LGT also obtains good results as shown in Fig. 5. Two instances are solved by LGT whereas 3BCID solves them with an extended timeout (disc2 in 3885 s and ex_6_2_10 in 3463 s). However, one instance (ex8_4_5) is solved by 3BCID but is not solved by LGT, even when the timeout is extended to 7200 s. Among all the filtering operators and all versions of TEC tested, only 3BCID reaches the optimal solution of this instance in the extended timeout 7200 s.

VI. RELATED WORK

Polyhedral convexification of the solution set consists in computing an outer approximation of the solution space of a system by convexifying each inequality constraint individually. A polyhedral enclosure of the solution set leads to a reduction of the search space and/or an improvement of the lower bound by using a call to a linear programming solver. The idea behind Graham-TEC is significantly different since the linear constraints learned depend on the filtering operation. They do not correspond to the initial constraints.

In Pelleau et al. [START_REF] Pelleau | Octagonal Domains for Continuous Constraints[END_REF], the notion of box encountered in continuous constraint programming is called into question. Each pair (i, j) of variables is associated with an octagon. This octagon can be viewed as the intersection of two boxes: the initial box and a box rotated by 45 degrees. For each pair (i, j), 4 constraints ±x i ± x j ≤ c k are added. To translate the constraints in the rotated box, the variables x i and x j are replaced by their expressions in the new base. Then, domains are contracted with these additional constraints. The principles are significantly different from LGT and TEC, but a common point resides in the additional binary constraints that are not defined by the user.

An arborescent filtering operator, called Box-k, was proposed in 2009 [START_REF] Araya | Filtering Numerical CSPs Using Well-Constrained Subsystems[END_REF]. Several well-constrained n × n subsystems of equations were initially selected in the system. At each node of the search tree, a contraction, based on an arborescent exploration, is performed on each subsystem. Unlike the TEC tree, an interval Newton is also called at each node of the tree. Under some conditions, Newton can use the property that the subsystem is well-constrained to converge more quickly towards the solutions (Newton behaves like a global constraint on the subsystem). The main difference between TEC and Box-k is that TEC is called at each node visited on all the system whereas the subsystem handled by Box-k is structured. Box-k provides good results only on very sparse systems.

VII. CONCLUSION

We have proposed a new filtering operator TEC and several ways to exploit the TEC tree in global optimization: TEC-UB and Graham-TEC which improve the upper and the lower bounding, respectively. The adaptive version LGT can switch between TEC-UB and Graham-TEC during search. The experimental results underline that LGT is better than each of them. Finally, results show that LGT brings a significant speed-up on several state-of-the-art instances compared to the filtering operators HC4 and 3BCID.

LGT is implemented in IbexOpt which is rigorous only for inequalities. It would be interesting to adapt the current algorithm for making it work in rigorous global optimization codes (also for equations) like Icos [START_REF] Lebbah | An Efficient and Safe Framework for Solving Optimization Problems[END_REF] and GlobSol [START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF].

Algorithm 3 :

 3 TEC (box, g , sol , T ECnodes, Φ) 1 Boxes ← {Φ(box, g , sol )} / * queue (FIFO) of boxes * / 2 AtomicBoxes ← ∅ / * set of atomic boxes * / 3 #nodes ← 1 4 while Boxes = ∅ and #nodes ≤ T ECnodes do 5

Figure 1 .

 1 Figure 1.Difference between TEC and TEC-UB on a system with one variable x and the objective variable x obj . The region removed by contraction appears light greyed out and the set of solutions is dark greyed out (or red). [X] 1 , [X] 2 and [X] 3 are the leaf boxes of the TEC tree. Left: Just one call to In-X-Newton is achieved on the contracted box (i.e., the box hull of the leaf boxes). Right: Two calls to In-X-Newton are performed, the second one on the selected leaf box [X] 3 .

Figure 2 .

 2 Figure 2. Illustration of ConvexHull on the projection of the leaf boxes [X] 1 ,...,[X] 4 and the hull box [H] on (x i , x j ). N W , SW , SE and N E are the 4 quadrants. Left: The quadrants N W and SW are covered. SE is also covered since on the x j -axis we have d < γ × ω([x j ]), and on the x i -axis we have 0 < γ × ω([x i ]). Only the 4 points of N E are handled to learn constraints, the others are ignored. Right: The two lines represent the two constraints learned. The dashed one does not give a sufficient gain on one of the two dimensions (d < γ × ω([x i ])). It is not learned.

Figure 3 .

 3 Figure 3. Rigorous calculation on the quadrant SE. The dashed line is the one with a, b ∈ R. The second line corresponds to the rigorous approximation of the first line with a and b rounded to the floating-point numbers a and b.

7 p 8 case NE 9 pNW 11 pSW 13 p

 7891113 ← a × x -y ≤ -b ← a × x -y ≥ -b 10 case ← a × x -y ≥ -b 12 case ← a × x -y ≤ -b 14 if the previous points satisfy p then P ← P ∪ {p} 15 return P ⊆ P

Figure 4 .

 4 Figure 4. Pairwise comparison between LGT and HC4. The coordinates of each point represents the CPU time (in second) required by both competitors to solve one instance.

Figure 5 .

 5 Figure 5. Pairwise comparison between LGT and 3BCID (parameter s=10 splits). The coordinates of each point represents the CPU time (in second) required by competitors to solve an instance.

Table I COMPARISON

 I OF TEC-UB, GR A H A M-TEC ANDLGT VERSUS TEC.

			Times			Nodes
		Total Gain Max Loss Max Gain Total Gain Max Loss Max Gain
	TEC-UB	26.4 %	40.7 %	71.8 %	28.1 %	53.1 %	80.7 %
	Graham-TEC 15.7 %	64.0 %	41.4 %	39.3 %	26.1 %	69.6 %
	LGT	35.4 %	27.3 %	88.3 %	41.2 %	51.5 %	91.3 %

SmearSumRel is a variant of the well-known smear function[START_REF] Kearfott | INTBIS, a portable interval Newton/Bisection package[END_REF]. The smear(x i , c j ) function measures an impact of the variable x i on the function of c j , depending on the width of x i and on the partial derivative of the c j function w.r.t. x i evaluated on the box [X].

Instances with more than 50 variables are generally not solved in reasonable time by Interval Branch & Bound algorithms without a tuning of the algorithm specific to each instance.