Christophe Besse
email: christophe.besse@math.univ-toulouse.fr

Benoît Mésognon-Gireau
email: benoit.mesognon-gireau@insa-toulouse.fr

Pascal Noble
email: pascal.noble@math.univ-toulouse.fr

Artificial boundary conditions for the linearized Benjamin-Bona-Mahony equation

We consider various approximations of artificial boundary conditions for linearized Benjamin-Bona-Mahoney equation. Continuous (respectively discrete ) artificial boundary conditions involve non local operators in time which in turn requires to compute time convolutions and invert the Laplace transform of an analytic function (respectively the Z-transform of an holomorphic function). In this paper, we derive explicit transparent boundary conditions both continuous and discrete for the linearized BBM equation. The equation is discretized with the Crank Nicolson time discretization scheme and we focus on the difference between the upwind and the centered discretization of the convection term. We use these boundary conditions to compute solutions with compact support in the computational domain and also in the case of an incoming plane wave which is an exact solution of the linearized BBM equation. We prove consistency, stability and convergence of the numerical scheme and provide many numerical experiments to show the efficiency of our tranparent boundary conditions.

Introduction

The Benjamin-Bona-Mahony (BBM) equation is a classical nonlinear, dispersive equation which model the unidirectional propagation of weakly nonlinear, long waves in the presence of dispersion. It is written as

(1) 1 - µ 6 ∂ 2 x ∂ t u + ∂ x u + 3ε 2 u∂ x u = 0, ∀t > 0, ∀x ∈ R.
where ε > 0 is the nonlinearity parameter, µ the shallowness/dispersion parameter and ε, µ have the same order (see [L] for more details on the derivation of this particular equation). It is usually proposed as an analytically advantageous alternative to the well-known Korteweg-de Vries equation

(2)

∂ t u + ∂ x u + 3ε 2 u∂ x u + µ 6 ∂ xxx u = 0, ∀t > 0, ∀x ∈ R.
Note that these two equations are asymptotically equivalent in the limit ε = µ → 0 since, ∂ 3 xxx u = ∂ 3 txx u + O(µ) but enjoy different properties. Moreover, they both possess solitary waves and cnoidal waves solutions. The most striking difference lies in the dispersive properties of each equation. Indeed, the dispersion relation for plane wave solutions of the linearized (ε = 0) (KdV) and (BBM) equations of the form u(x, t) = e i(kx-ωt) are given by

ω kdv (k) = k - µ 6 k 3 , ω bbm (k) = k 1 + µ 6 k 2 .
Hence, the sign of the phase speed depends on the wave number and is unbounded in the case of the Korteweg-de Vries equation whereas the phase speed is always non negative and bounded for the Benjamin-Bona-Mahony equation which is more consistent with the full water wave problem. The difference between the two equations is also seen regarding the regularizing effect of the dispersive terms on the solutions of the Burgers equation

∂ t u + ∂ x u + u∂ x u = 0.
Dispersive regularization of hyperbolic conservation laws is known to generate so-called dispersive shock waves. In contrast to their diffusive counterparts, dispersive shocks have an oscillatory structure and expand with time so that the Rankine-Hugoniot jump conditions are not satisfied. There is a huge literature on these particular patterns for the Korteweg-de Vries equation whereas less is known for the Benjamin-Bona-Mahony equation. In the case of the BBM equation, shock type solutions without oscillations can be constructed whereas they do not exist for the KdV equation: see [EHS] for more details.

The numerical simulation of such patterns is a hard task: usually, such equations are solved by using spectral techniques which are particularly suitable to describe oscillatory phenomena but suppose that periodic boundary conditions are imposed to the edge of the computational domains. Moreover, due to the fact that the oscillatory part of the dispersive shock wave expands in time, one has to take larger and larger computational domain which, in turn, imply high computational costs [K]. In addition, one should mention that the dynamic of dispersive equations is dramatically changed depending they are set on the whole space or in a periodic domain: in the latter case, small amplitude waves cannot scatter to infinity and stay in the computational domain forever.

Instead, one can imagine a more appropriate strategy based on the transparent boundary conditions (TBC): this consists in deriving suitable boundary conditions so that the solution calculated in the computational domain is an approximation of the exact solution restricted to the computational domain. These artificial boundary conditions are called absorbing boundary conditions (ABC) if they lead to a well-posed initial boundary value problem where some energy is absorbed at the boundary. See [AABES] for a review on the techniques used to construct such transparent or artificial boundary conditions for the Schrödinger equation.

In this paper, we focus on the linearized BBM equation

(3) ∂ t (u -ε∂ xx u) + c∂ x u = 0 ∀t > 0, ∀x ∈ R,
where ε is the dispersion parameter and c > 0 is a velocity. We complement this equation with an initial datum u(0, x) which has compact support in [x , x r ].

The aim of this paper is to derive transparent boundary conditions for (3),to give consistency error estimates and to prove convergence of the numerical approximation and to implement these boundary conditions numerically. At this stage, there are two distinct strategies. The first strategy consists in computing the continuous transparent boundary conditions and then implement these conditions numerically. Though, there is no clear strategy about how to discretize these boundary conditions which may lead to ill-posed numerical problem (see [H] for more details in the case of the wave equation). A more robust strategy consists in discretizing the equations both in time and space and then deriving the suitable artificial boundary conditions for the fully discrete problem by using the Z-transformation: see [BEL-V] for an application in the case of the linearized KdV equation. This is the strategy we will follow in this paper. We shall mention here that in contrast to [BEL-V], we can compute these discrete transparent boundary conditions explicitly.

The paper is organized as follow. In section 2, we apply the ideas found in [ZWH] to obtain the continuous transparent boundary conditions for the linearized BBM equation (3). In section 3, we present two appropriate space and time discretization of (3) and explain the procedure to derive discrete transparent boundary conditions. We focused on and compare two spatial discretizations of the convection term which are the upwind scheme and the centered scheme. The first scheme is of particular importance when we will later have to deal with the nonlinear hyperbolic term. The centered scheme is higher order but one has to be careful in a nonlinear setting where it can generate numerical dispersion and instability. In both cases, we prove the stability, the consistence of the discrete transparent boundary conditions with the continuous ones and the convergence of the numerical scheme. Note that this latter point was already carried out (even at the non linear level) in [AKO] but with simpler Dirichlet boundary conditions which are not suitable to deal with travelling wave solutions (wave packets, solitons, dispersive shocks). Finally, in section 4, we present some numerical simulation to illustrate our findings: we performed three types of simulation. The two first numerical exemples are directly inspired by the theory: we consider a Gaussian and a wave packet initial data. The last exemple is inspired by problem of wave generation in a computational domain. Here, we show how to inject a travelling wave solution of the linear BBM in the computational domain.

Transparent boundary conditions in the continuous case

In this section, we show how to derive transparent boundary conditions in the continuous case and show their stability: this will be a guideline in the construction of discrete transparent boundary conditions. The construction of continuous artificial boundary conditions associated equation ( 3) is established by considering the problem on the complementary of [x , x r ],

∂ t (u -ε∂ xx u) + c∂ x u = 0, ∀t > 0, ∀x < x or ∀x > x r , u(0, x) = 0, ∀x < x or ∀x > x r , u(•, x) → 0, x → ±∞.
The Laplace transform of u with respect to time, defined for all s ∈ C such that (s) > 0 by

L(u(., x))(s) = ∞ 0 exp(-s t) u(t, x)dt satisfies s L (u(., x)) -ε∂ 2 x L (u(., x)) + c ∂ x L (u(x, .)) = 0, ∀x < x or ∀x > x r , L (u(., x)) → 0, x → ±∞.
Then L(u(., x)) is given by

L (u(., x)) (s) = A + r exp (ξ + (s)(x -x r )) + A - r exp (ξ -(s)(x -x r )) , ∀x > x r , L (u(., x)) (s) = A + exp (ξ + (s)(x -x )) + A -exp (ξ -(s)(x -x )) , ∀x < x , with ξ ± (s) = 1 2ε c s ± c s 1 + 4ε c 2 s 2
, where the branch of the complex square root is taken such that its real part is positive. We assume that c > 0: then , for all s ∈ C such that (s) > 0, one has (ξ + (s)) > 0 and (ξ -(s)) < 0. As a consequence, the condition L (u(., x)) as x → ±∞ implies A - l = 0 and A + r = 0 so that one has the boundary conditions:

(4) ∂ x L (u(., x r )) (s) = ξ -(s)L(u(., x r ))(s), ∂ x L (u(., x )) (s) = ξ + (s)L (u(., x )) (s).
Back to the original coordinates, one finds

(5)

∂ 2 tx u(t, x r ) = 1 2ε c u(t, x r ) - 1 2 √ ε t 0 J 0 ( c s 2 √ ε )(c 2 + 4ε∂ 2 t )u(t -s, x r )ds , ∂ 2 tx u(t, x ) = 1 2ε c u(t, x ) + 1 2 √ ε t 0 J 0 ( s c 2 √ ε )(c 2 + 4ε∂ 2 t )u(t -s, x )ds ,
where J 0 is the Bessel function of the first kind defined for all t ∈ R by J 0

(t) = 1 2π π -π exp(-i t sin(θ))dθ.
We prove that these boundary conditions satisfy the following property:

Proposition 2.1. The linearized BBM equation ( 3) with the boundary conditions ( 5) is

L ∞ t (L 2 x ([x ; x r ])) stable: xr x u 2 + ε(∂ x u) 2 (t, y)dy - xr x u 2 + ε(∂ x u) 2 (0, y)dy + c 4π R 1 - 4ε c 2 ξ 2 |L(U (., x ))(ξ)| 2 + |L(U (., x r ))(ξ)| 2 dξ = 0, with U (t, x) = u(t, x)1 [0,T ] (t)
for T > 0 sufficiently large.

Proof. We multiply equation ( 3) by u and integrate on the interval [x , x r ]:

d dt xr x u 2 (t, y) 2 dy -ε xr x ∂ 3 txx u u(t, y)dy + c u 2 2 xr x = 0.
Then, by integrating by parts, one finds

(6) d dt xr x u 2 (t, y) 2 + ε (∂ x u) 2 (t, y) 2 dy + c u 2 2 -ε∂ 2 xt u u xr x = 0. Let us denote E(t) = xr x u 2 (t, y) 2 + ε (∂ x u) 2 (t, y) 2 dy.
Then, by integrating (6) with respect to time, one finds

E(t) -E(0) + t 0 c u 2 2 -ε∂ 2 xt u u xr x dt = 0. Set J = t 0 c u(x , s) 2 2 -ε∂ 2 xt u(x , s) u(x , s)ds and U (t) = u(t, x )1 [0,T ] , V (t) = ∂ x u(t, x )1 [0,T ] with T
sufficiently large: one has

J = t 0 c U (t) 2 2 -εU (t)V (t)dt = R c U (t) 2 2 -εU (t)V (t)dt = 1 2 π R c |L(U )(ξ)| 2 2 -iξεL(U )(ξ)L(V )(ξ)dξ = 1 2 π R c |L(U )(ξ)| 2 2 -iξεL(U )(ξ) 1 2ε c iξ + c iξ 1 - 4ε c 2 ξ 2 L(U )(ξ)dξ = - c 4π R 1 - 4ε c 2 ξ 2 |L(U )(ξ)| 2 dξ.
Similarly, if one denotes

U (t) = u(t, x r )1 [0,T ] , V (t) = ∂ x u(t, x r )1 [0,T ]
with T sufficiently large: one has

J r = t 0 c u 2 2 -ε∂ 2 xt u(x r , s) u(x r , s)ds = c 4π R 1 - 4ε c 2 ξ 2 |L(U )(ξ)| 2 dξ,
which completes the proof of the proposition.

Though being completely explicit, the transparent boundary conditions (5) are hardly implemented numerically. A major issue is to make the discretization of those boundary conditions with the numerical scheme used to discretize the BBM equation. An alternative approach consists in deriving directly discrete transparent boundary conditions. For that purpose, we mimic the previous analysis at the discretized level: the analogous of the Laplace transform is the Z-transform whereas the ODE problem is transformed into a linear recurrence problem.

Discrete Artificial boundary conditions

In this section, we present two different Crank Nicolson type space-time discretizations of the linear BBM equation (3) on the bounded space domain [x ;

x r ] and we implement discrete boundary conditions for the two schemes. The first discretization leads to an upwind transport and centered dispersion scheme and the second one to a centered transport and dispersion scheme. These two schemes preserve the L 2 stability of the equation (see later Proposition 3.6 and 3.11) which is essential for this kind of hyperbolic structure.

As one shall see in Section 4, the numerical implementation of the second scheme gives a better approximation of the exact solution than the first one. In particular, the convergence rate with respect to the space step is only of order 1 in the first case compared to an order 2 in the second case. However, the centered dispersion/upwind transport scheme is more adapted to deal with a non-linearity of the form u∂ x u and thus in view of working with the fully non-linear BBM equation (1), we choose to present both schemes in this paper.

Let us now introduce some notations that will be used throughout the paper. We choose a time step δt > 0 and a space step δx > 0. We set:

J + 1 = x r -x δx .
The goal is to compute an approximation u n j of the exact solution u of the linear BBM equation (3) on the finite domain [x ; x r ] at points jδx and instants nδt with 0 ≤ j ≤ J + 1 and n ∈ N:

∀n ∈ N, ∀j ∈ [0; J + 1],
u n j ≈ u(nδt, x + jδx).

Note that the choice of J + 1 instead of J is made to have J "interior points" and two "exterior points" corresponding to x = x and x = x r . This latter two points play an important role in the derivation of the discrete transparent boundary conditions.

Centered dispersion/Upwind transport

In this section, we consider the following Crank Nicolson discretization of the linear BBM equation (3):

u n+1 j - ε δx 2 (u n+1 j+1 -2u n+1 j + u n+1 j-1 ) -u n j - ε δx 2 (u n j+1 -2u n j + u n j-1 + cδt 2δx (u n+1 j -u n+1 j-1 ) + (u n j -u n j-1 ) = 0, ∀ j ∈ N. (7)
For the sake of clarity, one sets:

α = ε δx 2 , λ = cδt 2δx .
The numerical scheme then reads for j ∈ N:

-(α + λ)u n+1 j-1 + (1 + 2α + λ)u n+1 j -αu n+1 j+1 = -(α -λ)u n j-1 + (1 + 2α -λ)u n j -αu n j+1 .
In view of a numerical approximation of the solutions of the linear BBM equation on the bounded space domain [x ;

x r ], one looks for the computation of u j n for all n ∈ N and all j ∈ [0; J + 1]. The equations for j = 1 and j = J involve the values u k 0 and u k J+1 for k = n and k = n + 1 which corresponds to discrete boundary conditions. The goal of this section is to derive transparent boundary conditions for the upwind Crank Nicolson scheme (7).

We first compute these boundary conditions in terms of Z-transforms of the sequences (u n j ) n∈N for j = 0 and j = J + 1. We then give an explicit formulation of these boundary conditions. Finally, we prove that this scheme coupled with the derived discrete boundary conditions is convergent to the solution of the continuous boundary problem with continuous transparent boundary conditions (5).

Discrete transparent boundary conditions for the upwind Crank Nicolson scheme

The discrete equivalent of the Laplace transform of the continuous case (see Section 2) is the Z-transform. We denote ûj (z) the Z-transform of (u n j ) n∈N :

ûj (z) = Z{(u n j ) n∈N }(z) = ∞ k=0 u k j z -k , |z| ≥ R > 0
(see also Definition A.1). In order to derive the discrete boundary conditions, we apply the Z-transform to the numerical scheme (7) and we obtain, using the shift result of Theorem A.2:

(8) 2αû j+1 -2(1 + 2α + λδ(z))û j + 2(α + λδ(z))û j-1 = 0, with (9) δ(z) = z + 1 z -1 .
Note that the initial value does not enter into account here since u 0 j = 0 for all j ≤ 0 and j ≥ J + 1 (recall that the initial data is compactly supported in [x r , x ]). From now on we choose R > 1 in the definition of the Z-transform u given by Definition A.1. The solutions of the second order linear recursion (8) are written as ( 10)

ûj (z) = A + r + (z) j + A -r -(z) j , ∀j ≤ 0, ∀j ≥ J + 1,
where the roots r ± (z) are given by

(11) r ± (z) = 1 + 2α + λδ(z) ± 1 + 4α + 2λδ(z) + λ 2 δ(z) 2 2α .
In order to describe the transparent boundary conditions, we first need to determine the modulus of r ± (z) for all z such that |z| > ∆ > 1. Note that r ± (z) are smooth functions of λδ(z). In what follows, we set

p = λδ(z) ∈ C + = {z ∈ C, | (z) > 0} (notice that the function δ defined by (3.1.1) maps the disc D(0, R) into C + ).
Proposition 3.1. For all p ∈ C + , the functions r + (p) and r -(p) defined as

r + (p) = 1 + 2α + p + 1 + 4α + 2p + p 2 2α , r -(p) = 1 + 2α + p -1 + 4α + 2p + p 2 2α satisfy |r -(p)| < 1 < |r + (p)|.
Proof. In order to simplify notations, we set a = 2α > 0. First, let us show that |r ± (p)| = 1 for all p ∈ C + . Assume that r ± (p) = e iθ . Then, the following equation is satisfied

ae 2iθ -2(1 + a + p)e iθ + (a + 2p) = 0.
Necessarily, one has p = a + 1 -a cos(θ) (cos(θ) -1) -i sin(θ)

. By introducing g(θ) = cos(θ) -1, we obtain

p = g(θ) 1 -ag(θ) g 2 (θ) + sin 2 (θ)
.

We have

1 ≤ 1 -ag(θ) ≤ 1 + 2a and -2 ≤ g(θ) ≤ 0. This implies that (p) ≤ 0 which is in contradiction with (p) > 0. As p → 1, one has |r + (p)| → ∞ and |r -(p)| → 1 2 < 1. Then,
(12) û1 (z) = r + (z)û 0 (z), ûJ+1 (z) = r -(z)û J (z).
The product of the two roots is 2αr

-r + = 2(α + λδ(z)) which gives r -(z) = 2(α + λδ(z)) 1 + 2α + p + 1 + 4α + 2λδ(x) + λ 2 δ(z) 2 .
Next, defining

α = 1 + 4α, ξ 2 = λ 2 + 2λ + α, η 2 = λ 2 -2λ + α, Λ = ξ η , µ = α -λ 2 ξη ,
we have

1 + 4α + 2λδ(x) + λ 2 δ(z) 2 = η 2 (z -1) 2 (Λ 2 z 2 -2Λµz + 1).
The expression for µ is given by

µ = α -λ 2 α2 + 2( α -2)λ 2 + λ 4 . Since α = 1 + 4α > 1, one has |µ| < 1.
The new expressions of the boundary conditions are ( 13)

2α(z -1)û 1 = (1 + 2α + λ)z û0 + (λ -1 -2α)û 0 + η Λ 2 z 2 -2Λµz + 1û 0 and (14) (2(α + λ)z + 2(λ -α))û J = (1 + 2α + λ)z ûJ+1 + (λ -1 -2α)û J+1 + η Λ 2 z 2 -2Λµz + 1û J+1 .

Explicit transparent boundary conditions

In this section, we write the boundary conditions ( 13) and ( 14) with original coordinates. This is done by inverting explicitly the expression ( 13) and ( 14) via the Z-transform. Let us write:

Λ 2 z 2 -2Λµz + 1 = (Λz -2µz + Λ -1 z -1 ) Λz Λ 2 z 2 -2Λµz + 1
and using the fact that (see Appendix A for more details)

Z -1 Λz Λ 2 z 2 -2Λµz + 1 = Λ -n P n (µ), |z| > R = max{|Λ -1 z 1 |, |Λ -1 z 2 |}
where z 1,2 are the roots of Λ 2 z 2 -2Λµz + 1, we have

(15) (1 + 2α + λ + ηΛ)u n+1 0 -2αu n+1 1 = (1 + 2α -λ)u n 0 -2αu n 1 -η n p=0 s p u n-p 0 and (16) -2(α + λ)u n+1 J + (1 + 2α + λ + ηΛ)u n+1 J+2 = -2(α -λ)u n J + (1 + 2α -λ)u n J+1 -η n p=0 s p u n-p J+1
where:

(17) ∀p ∈ N, s p (µ) = Λ -p (P p+1 (µ) -2µ P p (u) + P p-1 (µ))
with the convention P -1 (µ) = 0. For a numerical implementation of the discrete boundary conditions ( 15) and ( 16), one must check that the coefficients s p in the convolution products are bounded uniformly with respect to p ∈ N. One can use the following decay property of the Legendre polynomials [ET]:

Lemma 3.2. Set µ = cos(θ) ∈] -1, 1[ with θ ∈]0, π[, then P n (µ) = √ 2 √ π sin(θ) cos ((n + 1/2)θ -π/4) √ n + O(n -3/2 ).
The bound for the error terms holds uniformly in the interval η ≤ θ ≤ π -η with η > 0.

As a result, one obtains that (s n (µ)) n∈N satisfy

s n (µ) = O(Λ -n n -3/2 ) with Λ > 1
which is a really fast decay. However, one can derive a recursion formula for the coefficients s n . For that purpose, one recalls the recursion formula satisfied by the Legendre's polynomials:

Proposition 3.3. For all n ∈ N, let P n be the n-th Legendre's polynomial. Then, the following recursion formula stands:

∀n ∈ N * , (n + 1)P n+1 (x) = (2n + 1)xP n (x) -nP n-1 (x).
Using the recursion formula given by Proposition 3.3, it is possible to compute the terms of the sequence (s n ) n∈N by recursion. More precisely, one gets the following result: Proposition 3.4. One has, for all n ≥ 3:

s n = µ Λ 2n -1 n + 1 s n-1 - 1 Λ 2 n -2 n + 1 s n-2 .
Proof. Using the recursion formula given by Proposition 3.3 and the definition of s n given by ( 17), one gets, for all n ≥ 3:

s n = Λ -n ( P n+1 (µ) -P n-1 (µ) 2n + 1 ).
One can check using Proposition 3.3 that the following recursion formula stands:

∀n ≥ 3, P n+1 (µ) -P n-1 (µ) 2n + 1 = µ P n (µ) -P n-2 (µ) n + 1 - n -2 (n + 1)(2n -3) (P n-1 (µ) -P n-3 (µ)
and the result follows.

The recursion formula given by Proposition 3.4 avoid the computation of Λ -n for large values of n. Thus, the numerical computation of the s n is simplified by using this recursion, since it is necessary to use at least one recursion to compute the Legendre polynomials. Now let us discuss the numerical implementation of the upwind Crank Nicolson scheme with transparent boundary conditions ( 15) and ( 16). Defining

a 0 = 1 + 2α + λ, a -= α + λ, b 0 = 1 + 2α -λ and b -= α -λ, the linear system we have to solve is (18)        a 0 + ηΛ -2α -a -a 0 -α . . . . . . . . . -a -a 0 -α -2a -a 0 + ηΛ               u n+1 0 u n+1 1 . . . u n+1 J u n+1 J+1        =        b 0 -2α -b -b 0 -α . . . . . . . . . -b -b 0 -α -2b -b 0               u n 0 u n 1 . . . u n J u n J+1        -        η n p=0 σ p u n-p 0 0 . . . 0 η n p=0 σ p u n-p J+1       
.

Convergence results

In this section, we prove stability and convergence result for the upwind Crank Nicolson scheme (7) with boundary conditions ( 13) and ( 14). It is very important to note that throughout this section, we do not need the explicit expression of the discrete boundary conditions, but only their expression in term of Z-transform. Thus, the methods used in this section can be adapted to cases where the explicit inversion via the Z-transform of the transparent boundary conditions is not possible, as for instance in the Korteweg-de-Vries case (see [BEL-V]).

We first prove a stability result for the upwind Crank Nicolson scheme (7) with boundary conditions ( 13) and ( 14). For the linear BBM equation set on R t × [x ; x r ], multiplying the equation by u and integrating by parts on [x ; x r ] gives

∂ t ( 1 2 |u| 2 L 2 x ([x ;xr]) + ε 2 |∂ x u| 2 L 2 x ([x ;xr]) ) = - c 2 [u(x) 2 ] xr x + ε[∂ 2 tx uu] xr x .
This suggests that a natural energy for the linear BBM equation is given by 1 2 |u| 2

L 2 x ([x ;xr]) + ε 2 |∂ x u| 2 L 2
x ([x ;xr]) . One can derive a similar discrete energy estimate for the fully discrete equation, by using analogous derivatives operators. To this purpose, we define the following discrete operators: Definition 3.5. For all sequences (u n j ) (j,n)∈[0;J+1]×N we set, for all n ∈ N, all j ∈ [1, J]:

D + t u n j = u n+1 j -u n j δt , D + x u n j = u n j+1 -u n j δx , D - x u n j = u n j -u n j-1 δx .
Inspired by the continuous case, one can define an "energy" for the solutions of the fully discrete problem:

(19) ∀n ∈ N, E n (u) = J j=1 (u n j ) 2 + ε J+1 j=1 (D - x u j ) 2 .
We also define the following "discrete Sobolev space"

H 1 (0; J + 1) = {(w 0 , ..., w J+1 ) ∈ R [0;J+1] } endowed with the norm | • | H 1 (0;J+1) defined by |(w 0 , ..., w J+1 )| H 1 (0;J+1) = (E N ((w 2 j ) j )) 1/2 .
We now prove the following stability result:

Theorem 3.6. The upwind Crank-Nicolson scheme for the linearized BBM equation (7) with transparent boundary conditions ( 12) is H 1 stable in the following sense:

∀N ≥ 1, E N (u) ≤ E 0 (u).
Proof. Inspired by [AES], and as in the continuous case, we compute a discrete "time derivative" of the energy for all n ∈ N:

D + t (E n )(u) = D + t ( J j=1 (u n j ) 2 + ε J+1 j=1 (D - x u n j ) 2 ) = J j=1 (u n+1 j ) 2 -(u n j ) 2 δt + ε J+1 j=1 (D - x u n+1 j ) 2 -(D - x u n j ) 2 δt = J j=1 u n+1 j -u n j δt (u n+1 j + u n j ) + ε J+1 j=1 D - x (u n+1 j -u n j ) δt (D - x (u n+1 j + u n j )). (20)
We set, for all j ∈ [0; J + 1], all n ∈ N :

u n+1/2 j = u n+1 j + u n j .
We now use a discrete integration by part type result, which can be easily proved, for all sequences

(u j ), (v j ) ∈ R [0;J+1] : J j=1 u j D - x v j = - J j=1 (D + x u j )v j + u J+1 v J -u 1 v 0 δx .
By using this result in equation ( 20), one gets:

(21)

D + t (E n )(u) = J j=1 (D + t u n j )u n+1/2 j -ε J j=1 (D + t D + x D - x u n j )u n+1/2 j + ε δx (D + t D - x u n J+1 ))u n+1/2 J+1 - ε δx (D + t D - x (u n 1 ))u n+1/2 0 .
Recall that (u n j ) (j,n) is a solution of the upwind Crank Nicolson scheme (7) which can be written as:

D + t (u n j -εD + x D - x u n j ) + c 2 D - x u n+1/2 j = 0.
We then obtain:

D + t (E n )(u) = - c 2 J j=1 (D - x u n+1/2 j )u n+1/2 j + ε δx (D + t D - x u n J+1 ))u n+1/2 J+1 - ε δx (D + t D - x (u n 1 ))u n+1/2 0 . The term c 2 J j=1 (D - x u n+1/2 j )u n+1/2 j
is the discrete equivalent of the term c xr x u∂ x udx in the continuous case of the energy estimate. This latter term is c 2 [u 2 (x)] xr x . In the discrete case, we express the discrete space derivative operator as follow, for all sequence (v j ) ∈ R [0;J+1] :

D - x (v j ) = - (v j-1 + 2v j -v j-1 ) δx + D + x (v j ) = -δxD - x D + x (v j ) + D + x (v j ).
Note that we wrote the discrete upwind derivative as a viscosity term plus the discrete downwind derivative. The aim is to obtain a symmetric expression by discrete integration by parts. Thus, by using as before a discrete integration by parts, one can prove by an easy computation that:

c 2 J j=1 (D - x u n+1/2 j )u n+1/2 j = cδx 4 J-1 j=0 (D + x u n+1/2 j ) 2 + c 4δx ((u n+1/2 J ) 2 + (u n+1/2 0 ) 2 ).
We now take the summation of D + t (E n )(u) for n from 0 to N . One obtains:

E N -E 0 δt = A + B + C with A = - N n=0 cδx 4 J-1 j=0 (D + x u n+1/2 j ) 2 , B = N n=0 ε δx (D + t D - x u n J+1 ))u n+1/2 J+1 - c 4 (u n+1/2 J ) 2 , C = N n=0 - ε δx (D + t D - x (u n 1 ))u n+1/2 0 - c 4 (u n+1/2 0 ) 2 .
We now prove that E N -E 0 ≤ 0. One has obviously that A ≤ 0. To prove that B ≤ 0, we need to use precisely the relation between u J+1 and u J , given by the transparent boundary condition. For this purpose, we recall the Plancherel formula (see also Theorem A.4):

∞ n=0 u n v n = 1 2π 2π 0 ( u v)(e iϕ )dϕ.
In order to apply this formula to the finite sum, we define

u N j = N n=0 u j z n .
We now write:

N n=0 ε δ (D + t D - x u n J+1 ))u n+1/2 J - c 4 (u n+1/2 J ) 2 = 1 2π π -π ε δx 2 e iϕ -1 δt ( u N J+1 -u N J )(e iϕ ) u N J (e iϕ )e iϕ + 1dϕ - c 4δx 1 2π π -π u N
J (e iϕ ) u N J (e iϕ )(e iϕ + 1)e iϕ + 1)dϕ.

We now use the discrete boundary condition (12) satisfied by u J , to rewrite this latter term under the form:

N n=0 ε δx (D + t D - x u n J+1 ))u n+1/2 J - c 4 (u n+1/2 J ) 2 = - α 2πδt π -π | u N J+1 (e iϕ )| 2 2i sin(ϕ) 1 - 1 r -(e iϕ ) - λ 2α δ(e iϕ ) 1 |r -(e iϕ )| 2 dϕ
(just notice that replacing u N J+1 by u J+1 in the integral does not change anything, since the integral of the extra terms vanishes). This integral is real, so one only has to check that the real part is negative. One recalls that:

∀z = 1, r -(z) = 1 + 1 2α + λ 2α δ(z) - 1 2α 1 + 4α + 2λδ(z) + λ 2 δ(z) 2 .
By noticing that δ(z) = 1 i tan(ϕ/2) for z = e iϕ , one gets that B is of the sign of ( 22)

sin(ϕ)Im( 1 2α 1 + 4α + 2 λ i tan(ϕ/2) - λ 2 tan(ϕ/2) 2 ).

Now recall that

√

• denotes here the square root with positive real part. Therefore, the imaginary part of the square root of ( 22) is of the sign of -1 tan(ϕ/2) . Finally, B is of the sign ofsin(ϕ) tan(ϕ/2) which is negative on ] -π; π[. The same technique applies to prove that C ≤ 0, by using this time the relation between u 1 and u 0 given by ( 12).

Let us now state a consistency result. For the sake of clarity, in the following result, [x ; x r ] = [0; 1]. Proposition 3.7. Let δt, δx > 0. Let u be a smooth function such that u satisfies the transparent boundary condition (4). We define the Z-transform of u(•, x), for all x ∈ [0; 1] by:

∀z = 0, u(z, x) = ∞ n=0 u(nδt, x) z n .
One has, for all compact K ⊂ C + , all s ∈ K:

(e sδ -1) u(e sδt , δx) -(e sδt -1)r + (e sδt ) u(e sδt , 0) = δx(O(δx 2 + δt 2 )) (e sδt -1) u(e sδt , 1) -(e sδt -1)r -(e sδt ) u(e sδt , Jδx) = δx(O(δx 2 + δt 2 ))
where r ± is defined by (11), and O(δx 2 ) + O(δt 2 ) is a term bounded in absolute value by a constant of the form M (δt 2 + δx 2 ), where M only depends on u, K.

Remark 3.8.

-Proposition 3.7 states that the discrete boundary conditions are consistent with the continuous boundary conditions, up to a O(δt 2 ) + O(δx 2 ) term. Moreover, the error made by replacing u(nδt, 0) and u(nδt, (J + 1)δx) by their expressions given by the discrete boundary condition

(12) is of size δx 2 + δt 2 .
-The consistency result of the proposition is expressed with the (z -1) factor as in (13) and ( 14) since these two expressions are the ones which are actually implemented numerically.

-A remarkable result is that the upwind transport/centered dispersion Crank Nicolson scheme (7) is only consistent at order O(δt 2 + δx) with the linear BBM equation whereas the discrete transparent boundary conditions for this scheme are consistent at a higher order with the continuous transparent boundary conditions for the BBM equation (5).

Proof. First of all, we prove that the Z-transform is an approximation of the Laplace transform. More precisely, one has, for all smooth function f such that f (0) = 0, and all δt > 0:

(23) L(f )(s) = δt f (e sδt ) + O(δt 2 )
Indeed, by setting g(t, s) = e -st f (t), one finds

L(f )(s) = R e -st f (t)dt = ∞ n=0 (n+1)δt nδt e -st f (t)dt = δt ∞ n=0 g(nδt) + g((n + 1)δt) 2 - ∞ n=0 (n+1)δt nδt ((n + 1)δt -t)(t -nδt) 2 g (t)dt = δt ∞ n=1 g(nδt) -δt 2 1 0 u(1 -u) 2 ∞ n=0 δtg ((n + u)δt)) du = δt f (e sδt ) + O(δt 2 ),
where for all n ∈ N, t n ∈]nδt, (n + 1)δt[ and O(δt 2 ) denotes any term bounded in absolute value by M δt 2 where M only depends on f . Now, one computes, recalling the expression of r + given by (11):

u(z, δx) -r + (z) u(z, 0) = u(z, δx) - 1 2α + 1 + λδ(z) 2α + 1 2α (1 + λδ(z)) 2 + 4α u(z, 0) = u(z, δx) -u(z, 0) - 1 2α 1 + λδ(z) + (1 + λδ(z)) 2 + 4α u(z, 0).
We now replace z by e sδt . We notice that:

δ(z) = 2 + O(δt) sδt , α = ε δx 2 , λ = cδt 2δx
and therefore, one gets by using a Taylor expansion of u(z, δx) with respect to δx → 0:

(24) u(z, δx) -r + (z) u(z, 0) = δx δt δt(∂ x u(e sδt , 0) -( c 2εs + c 2εs 1 + 4εs 2 c 2 u(e sδt , 0)) - δx 2 δt δt( u(e sδt , 0) -ε∂ 2 x u(e sδt , 0) + c s ∂ x u(e sδt , 0)) + δx δt O(δx 2 + δt 2 )
We now use the relation ( 23) to get that the first line of ( 24) is the Laplace transform of the continuous transparent boundary condition satisfied by u up to a O(δt 2 ) term, and the second line is the Laplace transform of the linear BBM equation satisfied by u up to the same precision. Therefore, one gets the desired result by noticing that z -1 = δt + O(δt 2 ). The analysis about the right end points j = J, J + 1 is similar.

We can now prove the following convergence theorem:

Theorem 3.9. Let u be the solution of the linearized BBM equation (3) with boundary conditions (5).

We define the sequence (v n 0 , ..., v n J+1 ) n∈N by:

∀k ∈ [0; J + 1], ∀n ∈ N, v n k = u(nδt, kδx).
We define the sequence (u n 0 , ..., u n J+1 ) n∈N built by induction by the Crank-Nicolson scheme (7) with transparent boundary conditions (12) initialized with (u(0, 0), ...u(0, (J + 1)δx)). Then, one has:

∀T > 0, ∃C T , √ δx sup n,nδt<T |(u n 0 , ..., u n J+1 ) -(v n 0 , ...v n J+1 )| H 1 (0;J+1) ≤ C T (δx + δt 2 )
where

|(w 0 , ..., w J+1 )| H 1 (0;J+1) = (E N ((w 2 j ) j )) 1/2 and E N is defined by (19).
Proof. The proof is an easy adaptation of the proof of Theorem 3.6 with a use of the consistency result of Proposition 3.7. We fix a T > 0. We define w n = u n -v n for all n such that nδt ≤ T . It is easy to check, by using a Taylor expansion of the solution u of (3) on t n , x j that:

∀1 ≤ j ≤ J, ∀0 ≤ nδt ≤ T, D + t (w n j -εD + x D - x w n j ) + c 2 D - x w n+1/2 j = ε n j with (25) ∀1 ≤ j ≤ J ∀0 ≤ nδt ≤ T, |ε n j | ≤ C T (δt 2 + δx 2 )
where C T only depends on T and u. Then, using the same computation as in the proof of Theorem 3.6, one can obtain:

D + t (E n )(w) = - c 2 J j=1 (D - x w n+1/2 j )w n+1/2 j + ε δx (D + t D - x w n J+1 ))w n+1/2 J+1 - ε δx (D + t D - x (w n 1 ))w n+1/2 0 + J j=1 ε n j w n+1/2 j .
Still adapting the proof of Theorem 3.6, one finds, taking the summation for n = 0, N with (N +1)δt ≤ T :

E N +1 -E 0 δt = A + B + C + D with A = - N n=0 cδx 4 J-1 j=0 (D + x w n+1/2 j ) 2 , B = N n=0 ε δx (D + t D - x w n J+1 ))w n+1/2 J+1 - c 4 (w n+1/2 J ) 2 , C = N n=0 - ε δx (D + t D - x (w n 1 ))w n+1/2 0 - c 4 (w n+1/2 0 ) 2 and D = N n=0 J j=1
ε n j w n j .

One has obviously A ≤ 0. We control B and C using the same technique as in Theorem 3.6. For instance, for C, one writes, with the same notations as in the proof of Theorem 3.6:

- N n=0 ε δx (D + t D - x w n 1 ))w n+1/2 0 - c 4 (w n+1/2 0 ) 2 = - 1 2π π -π ε δx 2 e iϕ -1 δt ( w N 1 -w N 0 )(e iϕ ) w N 0 (e iϕ )e iϕ + 1dϕ - c 4δx 1 2π π -π w N
0 (e iϕ ) w N 0 (e iϕ )(e iϕ + 1)e iϕ + 1dϕ.

Now, one writes

w 1 (z) -w 0 (z) = (r + (z) -1) w 0 (z) + v 1 (z) -r + (z) v 0 (z)
and thus C can be written as:

C = - α 2πδt π -π | u N J+1 (e iϕ )| 2 2i sin(ϕ) r + (e iϕ ) -1 + λ 2α δ(e iϕ ) | u N 0 (e iϕ )| 2 dϕ - 1 2π π -π ε δx 2 δt (e iϕ -1)( v 1 -r + v 0 )(e iϕ ) w N 0 (e iϕ )e iϕ + 1dϕ.
Therefore, the first term of C is negative as in the proof of Theorem 3.6, and the second one is, according to Proposition 3.7, bounded by C δxδt (δx 2 + δt 2 ) √ E N (w). The same estimates goes for B. The term D is bounded, using a convexity identity by C T N δx (δt + δx 2 ) 2 + 1 2 N n=0 E n (w) + E n+1 (w) (the 1 δx factor comes from the summation for j = 1 to J of the estimate (25), with J + 1 = 1 δx ). Finally, one obtained (recall that N δt ≤ T ):

δxE N +1 ≤ δxE 0 + C T (δt + δx 2 ) 2 + δxδt 1 2 N n=0 E n (w) + E n+1 (w)
and thus one gets the Theorem by finite induction for N δt ≤ T , noticing that E 0 = 0.

Centered dispersion and advection

In this section, we consider the case of a centered space discretization and Crank Nicolson time discretization and derive exact discrete transparent boundary conditions. First, the Crank-Nicolson time discretization of (3) reads ( 26)

u n+1 -u n δt -ε ∂ 2 x u n+1 -∂ 2 x u n δt + c∂ x u n+1 + u n 2 = 0.

Discrete transparent boundary conditions for the centered Crank Nicolson scheme

We perform the Z-transform of (26), using she shift result of Theorem A.2:

(27) û -ε∂ xx û + cδt 2 z + 1 z -1 ∂ x û = 0.
Now we discretize in space by a centered finite difference scheme

(28) ûi -ε ûi+1 -2û i + ûi-1 δx 2 + cδt δx z + 1 4(z -1) (û i+1 -ûi-1 ) = 0.
In order to simplify notations, we set

α = ε δx 2 , λ = cδt 2δx , δ(z) = z + 1 z -1 . Then, equation (28) reads (29) (2α -λδ(z)) ûi+1 -2(1 + 2α)û i + (2α + λδ(z)) ûi-1 = 0.
The solutions of (29) are written as

(30) ûj (z) = A + r + (z) j + A -r -(z) j , ∀j ≤ 1, ∀j ≥ J, r ± (z) = 1 + 2α ± 1 + 4α + λ 2 δ(z) 2 2α -λδ(z) .
In order to describe the transparent boundary conditions, we first need to determine the modulus of r ± (z) for all z such that |z| > ∆ > 1. Note that r ± (z) are smooth functions of λδ(z). In what follows, we set

p = λδ(z) ∈ C + = {z ∈ C, | (z) > 0}.
Proposition 3.10. For all p ∈ C + , the functions r + (p) and r -(p) defined as

r + (p) = 1 + 2α + 1 + 4α + p 2 2α -p , r -(p) = 1 + 2α -1 + 4α + p 2 2α -p . satisfy |r -(p)| < 1 < |r + (p)|.
Proof. In order to simplify notations, we set a = 2α > 0. First, let us show that |r ± (p)| = 1 for all p ∈ C + . Assume that r ± (p) = e iθ . Then, the following equation is satisfied 31)

û1 (z) = r + (z)û 0 (z), ûJ+1 (z) = r -(z)û J (z).

Explicit transparent boundary conditions

In order to write discrete transparent boundary conditions for the numerical scheme, we have to invert the conditions (31) and boundary conditions are written with the help of discrete convolution that involve boundary terms. For that purpose, we write the boundary conditions ( 31) as

(32) (2α -p(z))û 1 (z) = (1 + 2α)û 0 (z) + 1 + 4α + p(z) 2 û0 (z), (2α + p(z))û J (z) = (1 + 2α)û J+1 (z) + 1 + 4α + p(z) 2 ûJ+1 (z).
We multiply both equations by z -1 so that (32) now read ( 33)

((2α -λ)z -(λ + 2α)) û1 (z) = (1 + 2α)(z -1)û 0 (z) + Λ z z 2 -2µz + 1 (z -2µ + z -1 )û 0 (z), ((λ + 2α)z + λ -2α) ûJ (z) = (1 + 2α)(z -1)û J+1 (z) + Λ z z 2 -2µz + 1 (z -2µ + z -1 )û J+1 (z),
where, in order to simplify notations, we set

Λ = 1 + 4α + λ 2 , µ = 1 + 4α -λ 2 1 + 4α + λ 2 . Now, since |µ| < 1, one can use the definition of Legendre polynomials (see Appendix A) to get z z 2 -2µz + 1 = ∞ n=0 P n (µ) z -n .
As a result, by using the convolution result of Theorem A.3, the boundary conditions ( 33) are written as:

(34)

(1 + 2α + Λ)u n+1 0 + (λ -2α)u n+1 1 = (1 + 2α) u n 0 -(λ + 2α)u n 1 -Λ n p=0 s p (µ) u n-p 0 , (1 + 2α + Λ)u n+1 J+1 -(λ + 2α)u n+1 J = (1 + 2α) u n J+1 + (λ -2α)u n J -Λ n p=0 s p (µ) u n-p J+1 ,
where s p (µ) = P p+1 (µ) -2µ P p (u) + P p-1 (µ) for all p ∈ N with the convention P -1 (µ) = 0. Note that we can also write the boundary conditions (34) as a directe discretization of the continuous transparent boundary conditions ( 5):

(35)

(u n+1 1 -u n 1 ) -(u n+1 0 -u n 0 ) δxδt = 1 2ε   c u n+1 1 + u n 1 2 + δx δt   (u n+1 0 -u n 0 ) + Λ(u n+1 0 + n p=0 s p (µ)u n-p 0 )     , (u n+1 J+1 -u n J+1 ) -(u n+1 J -u n J ) δxδt = 1 2ε   c u n+1 J + u n J 2 + δx δt   (u n+1 J+1 -u n J+1 ) -Λ(u n+1 J+1 + n p=0 s p (µ)u n-p J+1 )     .
By using the recurrence relation satisfied by the Legendre polynomials, one finds that s p (µ) is written as

s p (µ) = P p+1 (µ) -µ P p (µ) p , ∀p ∈ N * .
As a result, one obtains that (s p (µ)) p∈N satisfy s p (µ) = O(p -3/2 ). For a numerical implementation of such boundary conditions, we define a 0 = 1 + 2α, a -= λ -2α, a + = λ + 2α, so that the linear system we have to solve is

(36)        a 0 + Λ a + -a -2a 0 a + . . . . . . . . . -a -2a 0 a + -a -a 0 + Λ               u n+1 0 u n+1 1 . . . u n+1 J u n+1 J+1        =        a 0 -a - a + 2a 0 -a - . . . . . . . . . a + 2a 0 -a - a + a 0               u n 0 u n 1 . . . u n J u n J+1        -        Λ n p=0 σ p u n-p 0 0 . . . 0 Λ n p=0 σ p u n-p J+1       
.

Convergence results

We prove in this section some convergence results for the centered Crank Nicolson scheme (26) with transparent boundary conditions (31). All these results do not need the explicit expression of these transparent boundary conditions, but only their expressions in term of Z-transform. Let us check the stability of the discrete transparent boundary conditions. In order to simplify notations, we set v n j = (u n j + u n+1 j )/2. We prove the following H 1 -stability result:

Proposition 3.11. Let us set

E n = δx (u n 0 ) 2 + (u n J+1 ) 2 4 + δx J j=1 (u n j ) 2 2 + δx ε 2 J j=0 u n j+1 -u n j δx 2 ∀n ∈ N.
Then any sequence (u n j ) j∈{0,...,J+1},n∈N satisfying the Crank Nicolson scheme (37) and together with the boundary conditions ( 35) satisfies the energy estimate E N ≤ E 0 for all N ∈ N * .

Proof. Recall that the Crank Nicolson scheme reads

(37) u n+1 j - ε δx 2 (u n+1 j+1 -2u n+1 j + u n+1 j-1 ) -u n j - ε δx 2 (u n j+1 -2u n j + u n j-1 ) + cδt 2δx v n j+1 -v n j-1 = 0,
for all j = 1, . . . , J. We multiply equation ( 37) by v n j and sum over all j = 1, . . . , J. By using the discrete transparent boundary conditions ( 35), one finds

J j=1 (u n+1 j ) 2 2 - J j=1 (u n j ) 2 2 + ε 2   J j=0 u n+1 j+1 -u n+1 j δx 2 - J j=0 u n j+1 -u n j δx 2   + v n 0   u n+1 0 -u n 0 2 + Λ 2   u n+1 0 + n p=0 s p (µ)u n-p 0     (38) + v n J+1   u n+1 J+1 -u n J+1 2 + Λ 2   u n+1 J+1 + n p=0 s p (µ)u n-p J+1     = 0.
In order to simplify notations, we set

E n = δx (u n 0 ) 2 + (u n J+1 ) 2 4 + δx J j=1 (u n j ) 2 2 + δx ε 2 J j=0 u n j+1 -u n j δx 2 .
Note that E n is a consistent approximation of the energy E(n δt) introduced in section (2) up to order 1 with respect to the space variable. The equation ( 38) reads

(39) E n+1 -E n + Λδx 2   v n 0   u n+1 0 + n p=0 s p (µ)u n-p 0   + v n J+1   u n+1 J+1 + n p=0 s p (µ)u n-p J+1     = 0.
We sum equation ( 39) over all n = 0, . . . , N -1 with N ∈ N * : one finds (40)

E N -E 0 + Λδx 2   N -1 n=0 v n 0 (u n+1 0 + n p=0 s p (µ)u n-p 0 ) + N -1 n=0 v n J+1 (u n+1 J+1 + n p=0 s p (µ)u n-p J+1 )   = 0.
To finish the proof, it remains to determine the sign of the last term in (40). To this end, we define (for fixed N ) the two sequences

ũn := u n 0 , ∀n = 0, . . . , N, 0, n > N,
and

ṽn := u n * sn (µ) := u n+1 + n p=0 s p (µ) u n-p , ∀n ∈ N. The Z-transform Z[ũ n ](z) = û(z) is analytic for |z| > 0, since it is a finite sum. The Z-transform Z[ṽ n ](z) = v(z) then satisfies v(z) = (z + ŝ(z))û(z) with ŝ(z) = ∞ p=0 s p (µ) z p and is analytic for |z| ≥ 1.
Using Plancherel's Theorem A.4 for Z-transforms, we have

N -1 n=0 v n 0 (u n+1 0 + n p=0 s p (µ)u n-p 0 ) = 1 2π π -π e iθ + 1 2 e iθ + ŝ(e iθ )|û(e iθ )| 2 dθ.
We deduce that

N -1 n=0 v n 0 (u n+1 0 + n p=0 s p (µ)u n-p 0 ) = 1 2π π -π cos θ 2 cos θ 2 + ( 2(cos(θ) -µ)) |û(θ)| 2 dθ > 0.
Using a similar trick for the other boundary term, one finds

N -1 n=0 v n J+1 (u n+1 J+1 + n p=0 s p (µ)u n-p J+1 ) > 0.
From ( 40), we deduce that E N ≤ E 0 for all N ∈ N * which completes the proof of the proposition.

We now state the following consistency result:

Proposition 3.12. Let δt, δx > 0. Let u satisfies the linearized BBM equation ( 3) and the transparent boundary conditions (5). We define the Z-transform of (u(nδt, x)) n∈N for all x ∈ [0; 1] by:

û(z, x) = ∞ n=0 u(nδt, x) z n , ∀z = 0.
One has, for all compact set K ⊂ C + , all s ∈ K:

δx(z -1) û(e s δt , δx) -r + (e s δt )û(e s δt , 0) = O(δt 2 + δx 2 ), δx(z -1) û(e s δt , 1) -r -(e s δt )û(e s δt , 1 -δx) = O(δt 2 + δx 2 ).
where r ± is defined by (30) and O(δx 2 ) + O(δt 2 ) is a term bounded in absolute value by a constant of the form M (δt 2 + δx 2 ) where M only depends on u, K.

Proof. We focus on the left boundary condition. Recall that u, solution of ( 3) and ( 5), satisfies the boundary condition

L ∂u ∂x (. , x ) (s) = 1 2εs c + c 2 + 4ε s 2 L(u(. , x ))(s),
where

L(f )(s) = ∞ 0 e -s t f (t)dt. Since u(0, x) = ∂u ∂x (0, x) = 0 for all x ∈ {0, δx, 1 -δx, 1}, one has δtû(e s δt , x) = δt ∞ n=0 e -s nδt u(nδt, x) = L(u(. , x))(s) + O(δt 2 ), ∀x ∈ {0, δx, 1 -δx, 1}, δt ∂ û ∂x (e s δt , x) = δt ∞ n=0 e -s nδt ∂u ∂x (nδt, x) = L ∂u ∂x (. , x) (s) + O(δt 2 ), ∀x ∈ {0, δx, 1 -δx, 1},
Now, let us consider the error of consistency E defined as

δx -1 E = (2α(z -1) -λ(z + 1))û(z, δx) -(1 + 2α)(z -1) + (1 + 4α)(z -1) 2 + λ 2 (z + 1) 2 û(z, 0), = δt δx z + 1 2 2ε( 2(z -1) δt(z + 1) ) û(z, δx) -û(z, 0) δx -δx 2(z -1) δt(z + 1) û(z, 0) + c û(z, δx) -û(z, 0) δx - δt δx z + 1 2   c + c 2 + (4ε + δx 2 ) 2(z -1) δt(z + 1) 2   û(z, 0).
Now, let us set z = e s δt , one has

2E = δt δx sû(e sδt , 0) -sε ∂ 2 û ∂x 2 (e sδt , 0) + cδx ∂ û ∂x (e sδt , 0) + O(δx 2 + s 2 δt 2 ) , δt 2ε s ∂ û ∂x (e sδt , 0) -c + c 2 + 4ε s 2 û(e sδt , 0) + O(δx 2 + s 2 δt 2 ) = δxL ∂ ∂t (u -ε ∂ 2 u ∂x 2 )(., 0) + c ∂u ∂x (., 0) (s) + O(δx 2 + s 2 δt 2 ) -2ε sL ∂u ∂x (., 0) (s) -c + c 2 + 4ε s 2 L(u(., 0))(s) + O(δx 2 + s 2 δt 2 )
Since u is a solution of ( 3) and ( 5), one deduces that E = O(δt 2 + δx 2 ). We deal similarly with the right boundary and this completes the proof of the proposition.

We now state a convergence Theorem. Note that the convergence rate is better than one of Theorem 3.9 since the centered dispersion and transport Crank Nicolson scheme has a better consistency with the linear BBM equation than the centered dispersion and upwind transport one.

Theorem 3.13. Let u be the solution of the linearized BBM equation (3) with boundary conditions (5). Let δt, δx > 0. We define the sequence (v n 0 , ..., v n J+1 ) n∈N by:

∀k ∈ [0; J + 1], ∀n ∈ N, v n k = u(nδt, kδx).
We define the sequence (u n 0 , ..., u n J+1 ) n∈N built by induction by the Crank-Nicolson scheme (26) with transparent boundary conditions (33) initialized with (u(0, 0), ...u(0, (J + 1)δx)). Then, one has:

∀T > 0, ∃C T , √ δx sup n,nδt<T |(u n 0 , ..., u n J+1 ) -(v n 0 , ...v n J+1 )| H 1 (0;J+1) ≤ C T (δx 2 + δt 2 )
where |(w 0 , ..., w J+1 )| H 1 (0;J+1) = (E N ((w 2 j ) j )) 1/2 and E N is defined in Proposition 3.11.

Proof. The proof is the same as for Theorem 3.9. The only difference is that the centered dispersion and transport Crank Nicolson scheme ( 26) is consistent at order 2 in space and time with the linear BBM equation.

Numerical Results

In this section, we present various numerical results to test the accuracy of our transparent boundary conditions. The first numerical test is made with a gaussian initial datum with a small standard deviation.

The spectrum associated to this datum is therefore wide. We analyze the dispersive effects related to ε.

The second numerical scheme concerns a wave packet as initial datum. The spectrum is more narrow and the dispersive effects are different. Finally, we show how to take into account an incoming wave. We denote by UCN and CCN respectively the Centered dispersion/Upwind transport numerical scheme introduced in section 3.1 and the Centered dispersion and advection numerical scheme introduced in section 3.2.

Numerical example 1: Gaussian initial condition

In order to show the efficiency of the approximation of the artificial boundary conditions, we compute numerical solutions and compare them to reference solutions. We consider here the solution to (3) with the initial datum u 0 (x) = exp(-400(x -1 2 ) 2 ).

The computational domain is (t, x) ∈ [0, 1] × [0, 1] meshed with (N + 1) × (J + 2) nodes. The evolution of the solution depends on ε and c. In order to check the order of the numerical scheme, we define e (n) the relative 2 -error at time t = nδt given by:

e (n) = u ref (t n , •) -u n (•) 2 / u ref (t n , •) 2 ,
where u n is the solution to the numerical scheme and where we use trapezoidal rule to compute the 2norm. Since the equations are linear, we compute reference solutions by using fast Fourier transform on a larger domain [-40, 80]. Thanks to the definition of e (n) , we consider two error functions: the maximum of e (n) with respect to 0 < n < N and the 2 -error given by

rel.ErrT m = max 0<n≤N e (n) , rel.ErrL2 = δt N n=1 (e (n) ) 2 1/2
, which corresponds respectively to the discrete version of L ∞ t L 2 x and L 2 t L 2 x errors functions. For the UCN numerical scheme, the error functions rel.ErrT m and rel.ErrL2 is bounded by {rel.ErrT m, rel.ErrL2} ≤ C t δt 2 + C x δx. and for the CCN numerical scheme, we have {rel.ErrT m, rel.ErrL2} ≤ C t δt 2 + C x δx 2 .

We consider here a positive velocity c = 2 and various values of ε, namely ε = 10 -p , p = 1, 2, 3, 4. The smallest values of ε lead to higher oscillations in the propagative right part. We plot the behavior of the reference solution for ε = 10 -4 on Figure 1. Note that the characteristics in the (x, t) plane have all a slope close to 1/2 which is close to 1/c. It is noticeable that, as ε > 0 is increased, a part of the energy carried by the solution travel to the left even if phase velocity are always positive for the BBM equation. This is easily explained by the fact that the BBM equation is dispersive and the dispersion relation is not monotone with respect to k > 0 the spatial wavenumber. More precisely, if we consider plane wave solutions ae i(kx-ωt) , the dispersion relation is

ω(k) = ck 1 + εk 2
which leads to the phase and group velocities respectively

v ϕ = ω k = c 1 + εk 2 , v g = ω (k) = c(1 -εk 2 ) (1 + εk 2 ) 2 .
If c > 0, the phase velocity is always positive. However, if |k| > 1/ √ ε, the group velocity is negative. We can see the effect of this process on the Figure 2 when ε = 10 -3 .

In order to look at the behavior of the error functions with respect to δt, we take J = 2 19 which leads to δx = 2 -19 ≈ 2. 10 -6 . We are forced to take a very small value for δx to be sure that the dominating error term is linked to C t and more specifically when we study the UCN scheme. The errors are plotted on Fig. 3. The second order accuracy with respect to time step is preserved. We clearly identify the saturation process linked to the dominating error due to the diffusive behaviour with respect to δx. Since the CCN scheme has a second order accuracy with respect to δx, the error level are lower than that of UCN.

We plot on Fig. 4 the behavior of the errors with respect to δx for N = 10 4 . The accuracy is in

Numerical example 2: Small wave packet

In order to observe more clearly the dispersive behavior of the BBM equation (in particular negative group velocities), we consider here the solution to (3) with the initial datum

u 0 (x) = exp(-400(x - 1 2 ) 2 ) sin(20πx).
The evolution of the reference solution is plotted in Fig. 5 for c = 2 and ε = 10 -3 . The computational domain is (t, x) ∈ [0, 1] × [0, 1] meshed with (N + 1) × (J + 2) nodes. We clearly observe that a large part of the energy carried out by the wave packet travel to the left, illustrating again the fact that for wavenumbers larger than k c = 1/ √ ε, the group velocity is negative. On the other hand, we also see in front of the wave packet small amplitude oscillations/large scale travelling wave with a group velocity v g ≈ v φ (0) = c. We present the evolution of the errors with respect to δt and δx for ε = 10 -3 on Fig. 6. As in the first numerical experiment, we get the accuracy O(δx + δt 2 ) for UCN scheme and O(δx 2 + δt 2 ) for CCN scheme.

Numerical example 3: Incoming traveling waves

In this section, we show how to inject an incoming travelling wave in a computational domain through transparent boundary conditions. This is of particular interest for the applications we have in mind: in particular, we wish to inject particular solutions of the full BBM equations (likes solutions, cnoidal waves or trains of solitons) in the computational domain in order to describe more realistic situations. Here, we carry out such an analysis in the linear case. We follow here the method presented in [AMP] for the Schrodinger-Poisson system.

Let us denote u in (x, t) = α+β cos(kx-ω(k)t) a plane wave solution of the linearized BBM equation, α being the mean and β the amplitude of the wave. We search transparent boundary conditions for the linear BBM condition with an initial data u 0 satisfying u 0 (x) = u in (x), ∀x ≤ x and u 0 (x) = 0, ∀x ≥ x r . For that purpose, we introduce a cut-off function χ such that χ(x) = 1 ∀x ≤ x and χ(x) = 0 ∀x ≥ x r and decompose u as u(x, t) = χ(x)u in (x, t) + v(x, t). The new unknown v is compactly supported in [x , x r ] at time t = 0 and satisfies the BBM equation with a source term:

∂ t (v -ε∂ xx v) + c∂ x v = G ε (x, t), ∀x ∈ [x , x r ], with G(x, t) = ε χ (x)∂ t u in (x, t) + 2χ (x)∂ xt u in (x, t) -cχ (x)u in (x, t). The function G is compactly supported in [x , x r ].
The derivation of the transparent boundary conditions for v is similar to the homogeneous u in = 0. Back to the original coordinates, the transparent boundary conditions for u are written as (41)

∂ 2 ∂x∂t u(t, x r ) = 1 2ε c u(t, x r ) - 1 2 √ ε t 0 J 0 ( c s 2 √ ε )(c 2 + 4ε∂ 2 t )u(t -s, x r )ds , ∂ 2 ∂x∂t (u -u in )(t, x ) = 1 2ε c (u -u in )(t, x ) + 1 2 √ ε t 0 J 0 ( s c 2 √ ε )(c 2 + 4ε∂ 2 t )(u -uin)(t -s, x )ds .
The derivation of discrete transparent boundary conditions is completely similar. The continuous plane wave solution is replaced by the discrete plane wave solution u in n,j = α + β cos(j kδx -nω(k)δt), ω(k) = 2 δt arctan λ sin(kδx) 1 + εk 2 sinc 2 (kδx/2) and discrete transparent boundary conditions are applied to the sequence u n j -u in n,j both at j = 0 and j = J + 1 (with the notation u in n,J = u in n,J+1 = 0). We present on Figures 7 and 8 the results of simulation when one sets an incoming wave u in with α = 0 and β = 0.5 on the left boundary with k = 2πp for various p, with c = 4 and ε = 10 -3 . The computational domain is (t, x) ∈ [0, 10]×[-2, 20] meshed with N = 20000 and J = 2 15 . The initial datum is u 0 (x) = 0.

When the wave number is small compared (case p = 1) to the critical wave number k c = 1/ √ ε, the group velocity is close to the wave velocity and the incoming wave propagate into the computational domain. At the end of the computation, we then obtain the profile of the incoming wave which is spatially periodic. The method also work with incoming wave with a non zero spatial mean. When the wave number is increased, both the group velocity and phase velocity are smaller and it takes much more time for the incoming wave to propagate in the computational domain (cases p = 2, 3). For wavenumbers larger than k c , the "incoming" wave does not invade the computational domain: indeed, we only observe the part of the solution that describes the transition between the incoming wave at x = -∞ and the trivial solution u = 0 at x = +∞.

We also show on Figure 9 the propagation of an incoming front u in with α = 2, β = 0.5, p = 5 and ε = 10 -3 . We plot the evolution of the wave on Figure 10 at times t = 2.5, t = 5, t = 7.5 and t = 10. We do not see any reflection at boundaries which show the effectivness of our TBCs.

Conclusion

In this paper, we derived discrete transparent boundary conditions for the linear Benjamin-Bona-Mahony equation for a upwind and centered discretization of the convection term, the dispersive term being modeled with finite difference centered scheme. Both schemes are proved to be stable, consistant and convergent. Note that in our proof, we do not need an explicit form of the boundary conditions (rather their expression with the help of the Z-transform) and that this result could be extended to linear Korteweg-de Vries equations [BEL-V].

In addition, we showed how to use those transparent boundary conditions to deal with the problem of wave generation in water wave problems. In practice, we will have to deal with non-linear equations. In order to derive transparent boundary conditions in the nonlinear case, we will adapt our strategy to linear equations with variable coefficients and then adopt a fixed point strategy: see [AABES] for more details in the case of nonlinear Schrodinger equations. We shall use this strategy to study accurately the interaction of solitons in BBM equations like [E-McG, DP] (where non physical boundary conditions were used). We also expect that this strategy can also be used in the numerical simulation of Boussinesq systems of PDEs which are in some cases nothing but a system of linearized BBM equations only coupled by nonlinear convection terms [AABCW]. An other question of interest is of course to derive discrete transparent boundary conditions in the case of the Serre-Green-Naghdi equations [L] which are physically more relevant for the water wave problem.

A The Z-transform

The discrete analogous of the Laplace transform is the so called Z-transform. In this section, we list some results used in this paper concerning this transformation. More precise details on the Z-transform and its properties are given in [JEI].

  by using the continuity of p → |r + (p)| and p → |r -(p)| on C + , one finds |r -(p)| < 1 < |r + (p)|. This completes the proof of the proposition. Since |r + (z)| > 1 and |r -(z)| < 1, the discrete transparent boundary conditions are written as

  (a -p)e 2iθ -2(1 + a)e iθ + (a + p) = 0. Then, necessarily, one has p = i a + 1 -a cos(θ) sin(θ) which is in contradiction with (p) > 0. As p → a, one has |r + (p)| → ∞ and |r -(p)| → a a + 1 < 1. Then, by using the continuity of p → |r + (p)| and p → |r -(p)| on C + \ {a}, one finds |r -(p)| < 1 < |r + (p)|. This completes the proof of the proposition. Since |r + (z)| > 1 and |r -(z)| < 1, the discrete transparent boundary conditions are written as (

Figure 1 :"Figure 2 :

 12 Figure 1: Evolution of the reference solution for c = 2 and ε = 10 -4

Figure 3 :

 3 Figure 3: Evolution of the errors with respect to δt for c = 2 and various ε

Figure 4 :

 4 Figure 4: Evolution of the errors with respect to δx for c = 2 and various ε

Figure 5 :

 5 Figure 5: Evolution of the reference solution for ε = 10 -3

Figure 6 :

 6 Figure 6: Evolution of the errors with respect to δt and δx for ε = 10 -3

Figure 7 :

 7 Figure 7: Evolution of an incoming wave for p = 1 (left) and p = 2 (right)

Figure 8 :Figure 9 :

 89 Figure 8: Evolution of an incoming wave for p = 3 (left) and p = 5 (right)

Figure 10 :

 10 Figure 10: Evolution of a propagative incoming front at times t = 2.5, t = 5, t = 7.5 and t = 10

Research of C.B. was partially supported by the ANR project BoND ANR-13-BS01-0009-01. Research of P.N. was partially supported by the ANR project BoND ANR-13-BS01-0009-01.

where R is the radius of convergence of the series ∞ n=0 u n z n . Note that in the previous definition, the radius R depends a priori on the sequence (u n ). We now state a couple of analogue results of the Laplace transform properties in the context of the Z-transform. One has the following discrete differentiation Theorem: Theorem A.2. When these quantities are well defined, one has:

We also have the following behavior of the Z-transform regarding the convolution:

) and one has:

By a classical series-integral commutation, one gets the following Plancherel Theorem:

We now state the following inversion result (which is a consequence of Cauchy Integral Theorem):

Theorem A.5. Let (u n ) be a complex sequence and u its Z-transform with radius of convergence R u .

Then, one gets the following relation:

where R > R u and C(0, R) denotes the circle of center 0 and radius R.

A useful result is the following inverse Z-transform:

where z 1 , z 2 are the square roots of z 2 -2µz + 1 and P n is the n-th Legendre polynomial. We deduce by scaling the variables the following inverse formula:

for all |z| > max(|Λ -1 z 1 |, |Λ -1 z 2 |) where z 1 , z 2 are the square roots of z 2 -2µz + 1 and P n is the n-th Legendre polynomial.