High order numerical schemes for transport equations on bounded domains

This article is an account of the NABUCO project achieved during the summer camp CEMRACS 2019 devoted to geophysical fluids and gravity flows. The goal is to construct finite difference approximations of the transport equation with nonzero incoming boundary data that achieve the best possible convergence rate in the maximum norm. We construct, implement and analyze the so-called inverse Lax-Wendroff procedure at the incoming boundary. Optimal convergence rates are obtained by combining sharp stability estimates for extrapolation boundary conditions with numerical boundary layer expansions. We illustrate the results with the Lax-Wendroff and O3 schemes.


Context
The goal of this article is to propose a high order numerical treatment of nonzero incoming boundary data for the advection equation. The methodology is developed here for the one-dimensional problem but it is our hope that the tools used below will be useful for higher dimensional problems. We are thus given a fixed constant velocity a > 0, an interval length L > 0 and we consider the (continuous) problem: x ∈ (0, L) , u(t, 0) = g(t) , t ≥ 0 . (1.1) The requirements on the initial and boundary data, namely f and g, will be made precise below. The solution to (1.1) is given by the method of characteristics, which yields the formula: (1. 2) The question we address is how to construct high order numerical approximations of the solution (1.2) to (1.1) by means of (explicit) finite difference approximations. This problem has been addressed in [CL20] in the case of zero incoming boundary data (that is, g = 0 in (1.1)). The focus in [CL20] is on the outflow boundary (x = L here since a is positive), for which extrapolation numerical boundary conditions are analyzed. Fortunately for us, a large part of the analysis in [CL20] can be used here as a black box and we therefore focus on the incoming boundary. To motivate the analysis of this paper, let us present a very simple -though illuminating-example for which we just need to introduce the basic notations that will be used throughout this article. In all what follows, we consider a positive integer J, that is meant to be large, and define the space step ∆x and the grid points (x j ) j∈Z by ∆x := L J , x j := j ∆x (j ∈ Z) .
The interval (0, L) corresponds to the cells (x j−1 , x j ) with j = 1, . . . , J, but considering the whole real line {j ∈ Z} will be useful in some parts of the analysis. The time step ∆t is then defined as ∆t := λ ∆x, where λ > 0 is a constant that is fixed so that assumption 1.1 below is satisfied. We use from now on the notation t n := n ∆t, n ∈ N; the quantity u n j will play the role of an approximation for the solution u to (1.1) at the time t n on the cell (x j−1 , x j ).
We now examine an example where the exact solution to (1.1) is approximated by means of the Lax-Wendroff scheme. The approximation reads: λ a 2 (u n j+1 − u n j−1 ) + (λ a) 2 2 (u n j+1 − 2 u n j + u n j−1 ) , n ∈ N , j = 1, . . . , J , (1.3) where we recall that λ = ∆t/∆x is a fixed constant and a > 0 is the transport velocity in (1.1). The initial condition for (1.3) is defined, for instance, by computing the cell averages of the initial condition f in (1.1), namely: ∀ j = 1, . . . , J , u 0 j := 1 ∆x x j x j−1 f (x) dx . (1.4) Without any boundary, the Lax-Wendroff scheme is a second order approximation to the transport equation [GKO95]. We would like, of course, to maintain the second order accuracy property when implementing (1.3) on an interval. This implementation, however, requires, at each time iteration n, the definition of the boundary (or ghost cell) values u n 0 and u n J+1 . At the outflow boundary, we prescribe an extrapolation condition [Kre66,Gol77], the significance of which will be thoroughly justified in the next sections: approximate the trace u(t n , 0), it seems reasonable at first sight to prescribe the Dirichlet boundary condition: u n 0 = g(t n ) , n ∈ N . (1.6) In the case of zero incoming boundary data (g = 0), and for any sufficiently smooth initial condition f that is "flat" at the incoming boundary, the main result of [CL20] shows that the above numerical scheme (1.3), (1.4), (1.5), (1.6) converges towards the exact solution to (1.1) with a rate of convergence 3/2 in the maximum norm. Numerical simulations even predict that the rate of convergence should be 2, or at least close to 2, for smooth initial data. However, implementing the above numerical scheme 1 quickly shows that the rate of convergence falls down to 1 when g is nonzero and satisfies the compatibility conditions 2 described hereafter with the initial condition f . Our goal is to provide with a thorough treatment of nonzero incoming boundary data and to design numerical boundary conditions that recover the optimal rate of convergence in the maximum norm (at least, the same rate of convergence as the one in [CL20] for zero boundary data). The strategy is not new and is now referred to as the inverse Lax-Wendroff method. It consists, as detailed below, in writing Taylor expansions with respect to the space variable x close to the incoming boundary and then using the advection equation (1.1) to substitute the normal derivatives ∂ m x u(t, 0) for tangential derivatives ∂ m t u(t, 0), the latter being computed thanks to the boundary condition in (1.1). This strategy is available when the boundary is non-characteristic [BGS07].
The inverse Lax-Wendroff method is a general strategy that has been followed in various directions. We refer for instance to [TS10,ST17,FY13,VS15,DDJ18] for various implementations related to either hyperbolic or kinetic partial differential equations. In these works, most of the time, the incoming numerical boundary condition prescribes the ghost cell value u n 0 in terms of the boundary datum g but also of interior cell values u n j with j ≥ 1. This is the reason why stability is a real issue in these works, see for instance the discussion in [VS15, Section 4], and many rigorous justifications are still open. We develop here a simplified version of some of those previously proposed boundary treatments, but we rigorously justify the convergence with an (almost) optimal rate of convergence. As in [CL20], the key ingredient in our analysis is an unconditional stability result for the Dirichlet boundary conditions which dates back to [GT78,GT81], see an alternative proof in [CG11].

The inverse Lax-Wendroff method
We first fix from now on some notations. In all this article, we are given some fixed integers p, r ∈ N and consider an explicit two time step approximation for the solution to (1.1): (1.7) In (1.7), the numbers a −r , . . . , a p are defined in terms of the parameter λ and of the velocity a (see, for instance, (1.3) for which p = r = 1). These numbers are fixed, which means that (1.7) is linear with respect to (u n j ). For simplicity, we follow [CL20] and choose as initial data for (1.7) the cell averages of the initial condition f in (1.1). This means that the vector (u 0 1 , . . . , u 0 J ) is defined by (1.4). For (1.7) to define inductively (with respect to n) the vector (u n 1 , . . . , u n J ), we need to prescribe the ghost cell values u n 1−r , . . . , u n 0 and u n J+1 , . . . , u n J+p . As explained above, we focus here on the inflow boundary and we therefore follow the extrapolation boundary treatment of [CL20] for the outflow boundary. Namely, if we define the finite difference operator D − as: and its iterates D m − accordingly, we choose from now on an extrapolation order k b ∈ N for the outflow boundary and prescribe: (1.8) The example (1.5) corresponds to k b = 2 (recall p = 1 for the Lax-Wendroff scheme so there is only one ghost cell). It now remains to prescribe the inflow values u n 1−r , . . . , u n 0 . Unlike some previous works, we are going to prescribe Dirichlet boundary conditions, meaning for instance that the value u n 0 will be determined in terms of the boundary datum g only. Let us assume for a while that u n j is a second order approximation of u(t n , (x j−1 + x j )/2) where u is the exact solution (1.2) of the continuous problem (1.1). Then we formally have: where ≈ means "equal up to O(∆x 2 )", and we then use (1.1) to get: The last term ∆x/(2 a) g ′ (t n ) in the previous equality is precisely the correction that is required to recover the second order accuracy when dealing with the Lax-Wendroff scheme (compare with (1.6)). More generally speaking, we could have pushed further the above Taylor expansion and obtained as a final (formal !) result that u n 0 should be "close" (whatever that means !) to some quantity of the form: where K is a truncation order and α 0 , . . . , α K are numerical constants. The general form of the Dirichlet boundary conditions that we consider below is: where the α κ,ℓ 's are numerical constants which will play a role (together with the truncation order K) in the consistency analysis. There are two main choices which we discuss in this article. The first one is given in [TS10,VS15]: and is relevant if u n j is eventually compared in the convergence analysis with u(t n , (x j−1 + x j )/2), u being the exact solution (1.2). The other possible choice we advocate is: and is relevant if u n j is eventually compared (as in our main result, which is Theorem 1.2 below) in the convergence analysis with the average of u(t n , ·) on the cell (x j−1 , x j ). The truncation order K is discussed with our main result in the following paragraph.

Results
We assume that the approximation (1.7) is consistent with the transport operator and that it defines a stable procedure on ℓ 2 (Z).
Assumption 1.1 (Consistency and stability without any boundary). The coefficients a −r , ..., a p in (1.7) satisfy a −r a p = 0 (normalization), and for some integer k ≥ 1, there holds: For the Lax-Wendroff scheme (1.3), we have p = r = 1 if λ a = 1, the integer k equals 2, and (1.11) holds if and only if λ a ≤ 1. In Theorem 1.2 below and all what follows, the velocity a > 0, the length L > 0, the parameter λ = ∆t/∆x and the extrapolation order k b ∈ N at the outflow boundary are given. Subsequent constants may depend on them. The integer k ≥ 1 is also fixed such that assumption 1.1 holds. We consider the initial condition (1.4) and its evolution by the numerical scheme (1.7), (1.8), the inflow ghost cell values being given by: (1.12) Of course, prescribing (1.12) is meaningful only if g is sufficiently smooth (say, g ∈ C k−1 ). One could push further the Taylor expansion in (1.12) and consider higher order correctors but it would require further smoothness on g and it would eventually not improve our convergence result below, so fixing the truncation order K = k − 1 seems to be the most convenient choice. Our main convergence result is the extension of the main result in [CL20] to the case of nonzero boundary data.
Actually, the constant C in (1.13) is independent of L ≥ 1, which is consistent with the convergence result we shall prove below for the half-space problem on R + with inflow at x = 0. As in [CL20], the loss of 1/2 in the rate of convergence of Theorem 1.2 looks somehow artificial and is mostly a matter of passing from the ℓ ∞ n ℓ 2 j topology to ℓ ∞ n,j . Our next result examines a situation where the optimal convergence rate min(k, k b ) can be obtained. In order to simplify (and shorten) the proof of Theorem 1.3, we only examine here the case of a half-space with extrapolation outflow conditions. The extension of the techniques to the case of an interval is left to the interested reader. Theorem 1.3 (Optimal rate of convergence for the outflow problem). Under assumption 1.1 and under the additional assumption 3.2 stated hereafter, there exists a constant C > 0 such that for any final time T ≥ 1, any integer J ∈ N * , any data f ∈ H k+1 ((−∞, L)), the solution to the scheme: (1.14) satisfies the error estimate as long as k b < k.
In other words, the technical assumption 3.2 hereafter, which is verified on many examples such as the Lax-Wendroff and O3 schemes, allows to recover the optimal rate k b = min(k b , k) in the case k b < k. Of course, one would also like to improve the rate min(k b , k) − 1/2 in the case k b = k, which is clearly the most natural choice. However, in that case, both the interior and boundary consistency errors scale like ∆x k and, in the framework of assumption 1.1, stability in the interior domain is available only in the ℓ 2 j topology, so it is quite difficult to derive the convergence rate k in the ℓ ∞ j topology. Theorem 1.3 already indicates that combining the approach of [CL20] with other techniques (here, boundary layer expansions) may improve some results. We hope to deal with the case k b = k in the future.

Convergence analysis for the inverse Lax-Wendroff method
This Section is devoted to the proof of Theorem 1.2. In order to shorten the exposition, we shall use some results of [CL20] as a black box and we refer the interested reader to [CL20] for more details. Following [CL20], we shall prove Theorem 1.2 by using a stability estimate for (1.7), (1.8), (1.12) and a superposition argument, which amounts to considering separately two half-space problems: one in which there is only inflow at x = 0, and one for which there is only outflow at x = L. The novelty here is the nonzero inflow source term so we first deal with that case.

Convergence analysis on a half-line for the inflow problem
We focus here on the inflow source term, and therefore start by proving the main convergence estimate that is the new ingredient for the proof of Theorem 1.2.
Theorem 2.1 (Convergence estimate for the inflow problem). Under assumption 1.1, there exists a constant C > 0 such that for any final time T ≥ 1, for any J ∈ N * , for any initial condition f ∈ H k+1 ((0, +∞)) and boundary source term g ∈ H k+1 ((0, T )) satisfying the compatibility conditions: the solution (u n j ) j≥1−r,n∈N to the numerical scheme: where u is the exact solution to the half-line transport problem: Proof. For convenience, we first extend g into a function g ♭ ∈ H k+1 ((0, +∞)) and then define: Since f and g satisfy the compatibility conditions (2.1), we have f ♯ ∈ H k+1 (R), and the exact solution u to (2.3) is given by: Let us now define: which corresponds to the cell average of the exact solution to (2.3). With (u n j ) j≥1−r,0≤n≤T /∆t the solution to the numerical scheme (2.2), we define the error ε n j := u n j − w n j , that is a solution to: (2.4) The interior consistency error (e n+1 j ) j≥1,0≤n≤T /∆t−1 is easily estimated by means of the Cauchy-Schwarz inequality and Fourier analysis: where the final inequality comes from assumption 1.1 and the fact that the ratio ∆t/∆x is constant. Going back to the definition of f ♯ , we have obtained the bound: for some constant C that is independent of the final time T ≥ 1 and the data f and g.
We now turn to the boundary errors in (2.4), and wish to estimate the following quantities: Let us consider an integer n such that 0 ≤ n ≤ T /∆t − 1. From the definition of w n ℓ , ℓ ≤ 0, we have: where we have used the Taylor formula 3 . By the Cauchy-Schwarz inequality, we get: and we now apply the change of variables (x, y) → (x y, x) to get: Restricting to ℓ = 1 − r, . . . , 0, we have: Summing now with respect to n, we end up with the estimate: for some constant C that is independent of the final time T ≥ 1 and the data f and g.
We now apply the main stability estimate for the error problem (2.4), for which we refer to the seminal papers [GT78,GT81] and to the more recent works [CG11,CL20]: The conclusion of Theorem 2.1 then comes from the combination of the estimates (2.5) and (2.6).

Convergence estimate for the outflow problem
We recall the convergence estimate obtained in [CL20] for the complementary half-space problem where extrapolation conditions are prescribed at the outflow boundary. We refer the reader to [CL20] for the details.
Theorem 2.2 (Convergence estimate for the outflow problem [CL20]). Under assumption 1.1, there exists a constant C > 0 such that for any final time T ≥ 1, for any J ∈ N * and for any initial condition f ∈ H k+1 ((−∞, L)), the solution (u n j ) j≤J+p,0≤n≤T /∆t to the numerical scheme:

Proof of Theorem 1.2
It remains to combine the convergence estimates of Theorems 2.1 and 2.2 to prove Theorem 1.2. We use a slight modification of the superposition argument in [CL20] in order to cope with the nonzero incoming condition, but we basically follow the same lines. Let us consider a final time T ≥ 1 and some data f ∈ H k+1 ((0, L)), g ∈ H k+1 ((0, T )) that satisfy the compatibility conditions stated in Theorem 1.2. We consider some function χ ∈ C ∞ (R) such that χ(x) = 0 if x ≤ 1/3 and χ(x) = 1 if x ≥ 2/3. We then decompose the initial condition f as: Since (1 − χ(·/L)) f vanishes on (2 L/3, L), we can extend it by zero to the interval (L, +∞) and thus consider (1 − χ(·/L)) f as an element of H k+1 ((0, +∞)). Furthermore, the functions (1 − χ(·/L)) f and g satisfy the same compatibility conditions as f and g at t = x = 0. We can thus apply Theorem 2.1 to the sequence (v n j ) j≥1−r,0≤n≤T /∆t that is defined as the solution to the numerical scheme: We obtain the estimate: (2.7) where v is the exact solution to the transport problem: Similarly, we can view χ(·/L) f as an element of H k+1 ((−∞, L)) that vanishes on (−∞, L/3). Theorem 2.2 then shows that the solution (w n j ) j≤J+p,0≤n≤T /∆t to the numerical scheme: Using the support property of the function χ and the fact that the scheme (1.7) is explicit with a finite stencil, we find that for all time iteration n up to the threshold: there holds: Combining then the error estimates (2.7) and (2.8), we obtain: where u is the exact solution to (1.1). It remains, as in [CL20], to iterate in time the error estimate (2.9). We follow again the argument in [CL20]. For any time iteration n between N and 2 N , we split the solution (u n j ) 1−r≤j≤J+p,0≤n≤T /∆t to (1.7), (1.4), (1.8), (1.12) as the sum of the solution to the problem: and of the (presumably small) solution to the 'error' problem: Since the initial condition u(· − a t N ) and the boundary source term g(t N + ·) satisfy the compatibility conditions at the corner t = x = 0, we can apply the first step of the proof (leading to the error estimate (2.9)) for the (ũ N +n j ) part, and we apply the stability estimate of [CL20, Proposition 4.1] for the (ε N +n j ) part. This leads to the second error estimate: and, more generally, to: The end of the proof is the same as in [CL20] and we refer the interested reader to that reference for the details.

High order outflow boundary layer analysis
In the present section, we explain how the analysis of [BC17], which dealt with the case of the Dirichlet boundary condition at the outflow boundary, can be extended to the case of high order extrapolation (1.8).
The goal is to obtain an accurate description of the numerical solution close to the outflow boundary by means of a boundary layer expansion. The leading order term in the expansion corresponds to the exact solution to the transport equation. However, this leading order term does not satisfy the extrapolation condition (1.8), leading to a consistency error of magnitude O(∆x k b ) on the boundary. Under some mild structural assumption on the numerical scheme (1.7), we show below that this O(∆x k b ) error on the boundary gives rise to a boundary layer term which scales as O(∆x k b +1/2 ) in the ℓ 2 j norm. This gain of a factor ∆x 1/2 enables us to recover the optimal convergence rate k b in the maximum norm on the whole spatial domain for k b < k.

An introductive example
Let us go back for a while to the case of the Lax-Wendroff scheme (1.3), which we consider here on the left half space: At the outflow boundary, we impose the first order extrapolation condition (which corresponds to k b = 1 while k = 2 for the Lax-Wendroff scheme): We start with some smooth initial condition f defined on (−∞, L) which we project as a piecewise constant function: The exact solution to the transport equation on (−∞, L) with initial condition f is u(t, x) = f (x − a t) (recall a > 0). Hence the consistency analysis of the Lax-Wendroff scheme indicates that u n j reads: where the first term in the expansion on the right hand side yields an O(∆x 2 ) consistency error in the interior domain, but also an O(∆x) consistency error on the boundary. If we wish to push forward the above expansion, we need to take into account the boundary consistency error and introduce a corrector which will hopefully not alter the interior consistency error. This can be achieved by observing that the sequence: is kept unchanged by the Lax-Wendroff scheme on Z, and belongs to ℓ 2 (−∞, J) (we assume 0 < λ a < 1 so |κ| > 1). Hence, to remove the boundary consistency error, we can add a corrector on the right hand side of (3.1) in the following way: where w n is defined in such a way that the two first terms on the right hand side satisfy the first order extrapolation condition, that is: If f is sufficiently smooth, then w n is O(1) and the first corrector on the right hand side of (3.2) is O(∆x) in ℓ ∞ j . Note however that the ℓ 2 norm (in space) of this boundary layer corrector scales as ∆x 3/2 since the sequence (κ j ) is square integrable on (−∞, 0). Another important observation at this point is that defining w n requires the real number κ not to equal 1. This fact follows here from a mere verification but it is a general consequence of the analysis in [Gol77] of the Lopatinskii determinant associated with the boundary condition (1.8) (see also the proof of Lemma 3.5 below).
At this stage, the error analysis amounts to studying the system satisfied by the sequence: the main point being that there is no boundary forcing term, and since the added boundary layer corrector is O(∆x 3/2 ) in ℓ 2 , the initial condition and interior consistency errors will be O(∆x 3/2 ). Overall, the stability estimate for the Lax-Wendroff scheme with first order extrapolation at the boundary yields the convergence estimate: By the triangle inequality, we thus obtain: and this immediately gives the uniform convergence estimate: The above brief sketch is made complete and rigorous below in the general framework of Theorem 1.3.

Discrete steady states
Formalizing somehow the previous example in a more general framework, let us now introduce the following definition.
Definition 3.1 (Steady state for the numerical scheme). A sequence (v j ) j∈Z is called a (discrete) steady state for the scheme (1.7) if it is kept unchanged by the time iteration process on Z, that is, if it satisfies: In order to characterize the discrete steady states, it is natural to introduce the characteristic polynomial: (3.4) From the consistency property (1.10), any constant sequence is a discrete steady state for (1.7), the same property being available for the continuous model (namely, the transport operator). However, the discrete nature of the differentiation operator involved in the numerical scheme (1.7) allows the existence of many other discrete steady states. The latter play an important role when considering the half-space problem with some discrete boundary conditions. From the non-characteristic assumption a > 0, it follows that among the roots of A, X = 1 is always a simple root. Let us now introduce the whole set of (pairwise distinct) roots of A together with their multiplicities through the factorization of A in C[X]: (3.5) Clearly, looking at the degree of the polynomial A, one has the equality τ σ=1 µ σ = r + p .
For convenience, we order the roots of A with decreasing modulus: To make the analysis more intelligible, we will work under the following assumption, which was already present in [BC17].
Assumption 3.2. The characteristic polynomial A defined in (3.4) has a unique root (equal to 1) on the unit circle S 1 = {z ∈ C , |z| = 1}. In other words, we assume: As observed on the above example of the Lax-Wendroff scheme, the steady states we are looking at should decrease rapidly as j tends to −∞, so that they provide with a localized correction (near the boundary) to the usual convergence analysis and belong to ℓ 2 (−∞, J). We are therefore only concerned with those roots of A that have modulus larger than 1. Lemma 3.3 below gives the precise number of such roots (counted with their multiplicities). We refer to [BC17, Lemma 2.1] for the proof.  assumptions 1.1 and 3.2, letting κ 1 , . . . , κ τ + be the roots of A that belong to U := {z ∈ C , |z| > 1} with their corresponding multiplicities µ 1 , . . . , µ τ + , then one has τ + σ=1 µ σ = p . (3.7) A direct consequence of Lemma 3.3 is the following description of steady states for (1.7) that belong to ℓ 2 (−∞, J). The proof follows from the standard description of the set of solutions to the recurrence relation (3.3).
Lemma 3.4. The set of discrete steady states for the scheme (1.7) that belong to ℓ 2 (−∞, J) is the finite dimensional linear subspace spanned by the p linearly independent sequences ρ (σ,ν) : Equivalently, such discrete steady states in ℓ 2 (−∞, J) read: Let us detail the parametrization of the set of (stable) discrete steady states on the two main examples we are concerned with. For the Lax-Wendroff scheme (1.3), one has: The (two simple) roots of A are 1 and: with κ ∈ U assuming, as usual, 0 < λ a < 1. For the half space problem on (−∞, J), κ is therefore the unique stable root, and 1 counts as an unstable root (see [BC17]). In particular, assumption 3.2 is satisfied. The set of solutions to (3.3) that belong to ℓ 2 (−∞, J) is the one-dimensional subspace spanned by the sequence (κ j−J ) j∈Z .
Let us now consider the so-called O3 scheme, which is a convex combination of the Lax-Wendroff and Beam-Warming schemes, see [Str62,Des08]. We now have p = 1 and r = 2, and the scheme reads: with, again, 0 < λ a < 1. Assumption 1.1 is then satisfied (with k = 3). The roots of the corresponding characteristic polynomial A are: each of them being simple. The root κ − is the only one in U and κ + belongs to the open unit disk D, which is consistent with Lemma 3.4 (p = 1). In particular, assumption 3.2 is satisfied.

The boundary layer expansion. Proof of Theorem 1.3
We now start proving Theorem 1.3, and for that, we consider some initial condition f ∈ H k+1 ((−∞, L)) which, for convenience, we extend to the whole real line R as an element of H k+1 (R). Our aim is to compare the solution to the scheme (1.14) (which is set on a half line) with the piecewise constant projection of the exact solution to the transport equation. We thus introduce the notation: The consistency analysis in [CL20] of the scheme (1.14) amounts to considering the numerical scheme satisfied by the error (u n j − ω n j ). It is proved in [CL20] that the resulting boundary consistency errors have size O(∆x k b ), while the interior consistency errors have size O(∆x k ). Here we have k b < k so the worst term is on the boundary. Following the arguments in [BC17], we are therefore going to introduce a boundary layer corrector in order to remove the boundary consistency error, up to introducing new initial and interior consistency errors, whose size will be proven to be O(∆x k b +1/2 ) hence the final result of Theorem 1.3. Let us make this argument precise.
The consistent expansion of the numerical solution (u n j ) takes the form of a corrected version of (ω n j ), involving now a boundary layer term (v n j ) ∈ ℓ 2 (−∞, J) as for the above introductive example. The aim is to reduce the magnitude, at the boundary, of the following error: The definition of (v n j ) j≤J+p,n∈N is chosen so as to correct the consistency error at the boundary. The simplest way to do so consists in chosing (v n j ) j≤J+p,n∈N so as to get precisely in the ghost cells the relations (D k b − ε n ) J+ℓ = 0, ℓ = 1, . . . , p. From now on, we formulate the problem in such a way to normalize the generating sequences according to the value of J. In view of Lemma 3.4, the problem to be solved writes: where the sequences ρ (σ,ν) are defined in (3.8). Equivalently to (3.12), we can look for the boundary layer corrector (v n j ) j≤J+p,n∈N under the form (3.14) where p n,σ ∈ C µσ−1 [X] for all index 1 ≤ σ ≤ τ + . The existence of the corrector (v n j ) is given by the following result. We recall that in the framework of Theorem 1.3, there holds k b < k.
Lemma 3.5. Consider the initial condition f ∈ H k+1 ((−∞, L)) extended to the whole real line R. Then the boundary layer problem (3.12)-(3.13) admits a unique solution (v n j ) j≤J+p,n∈N , and this solution satisfies the estimate: where the constant C > 0 is independent of ∆x > 0, J, L and f .
Proof. Let us fix some integer n ∈ N. The solution (u n j ) j≤J+p to (1.14) solves the homogeneous boundary condition (1.8), thus equivalently to (3.13) one has to find the vector of coordinates z ∈ C p solution to the linear system , and the p × p matrix A k b is defined as follows: where we have relabeled the sequences ρ (σ,ν) , σ = 1, . . . , τ + , ν = 0, . . . , µ σ − 1 as ρ (1) , . . . , ρ (p) in order to make the definition of A k b easier to read. The latter matrix is somehow the k b th-order discrete derivative of the so-called confluent Vandermonde matrix. It seems possible to compute the determinant of A 0 , see [HJ94], and then to extend this result to higher values of k b but we prefer to avoid such complicated computations. From the identity of dimensions, we shall just prove that the matrix A k b is one-to-one, in other words we shall prove that the problem (3.12)-(3.13), or equivalently (3.14)-(3.13), admits a trivial kernel.
Dealing with discrete derivatives of products of polynomial and/or geometric sequences, the divided difference algebra appears as a suitable tool in our analysis. For more details we refer the interested reader to [Ste39,Pop40,dB05]. For consistency in the notation, we recall hereafter the recursive definition of divided differences, but specified for the case of consecutive integer abscissae. Being given a sequence of complex numbers (w j ) j∈Z , one defines: The quantity (D k b − w) j is directly related to the divided difference w[j − k b , . . . , j] by the equality: (3.17) Importantly, we may also use the Leibniz formula for divided differences of a product of two sequences: In terms of the D − operator, using the relation (3.17), the Leibniz formula (3.18) rewrites under the more recognizable form: Let us continue with the representation formula (3.9) of the solution to the boundary layer problem. Looking at the kernel of the linear problem (3.13), we have to find polynomials (p σ ) 1≤σ≤τ + with respective degrees less than or equal to (µ σ − 1) 1≤σ≤τ + , satisfying the set of equations: where we denote, with a slight abuse in the notation, κ σ for the corresponding geometric sequence (κ m ) m∈Z , for any σ = 1, . . . , τ + . Actually, from the identity (3.17) and by induction on the integer m (or using (3.16)), it is easy to prove that the m-th order divided difference of κ σ is given by: Let us introduce, for any integer σ and any polynomial p σ with degree less than or equal to µ σ − 1, the following polynomial Q σ also with degree less than or equal to µ σ − 1: With these notations, the equations to solve now equivalently read: (3.20) Actually, the above set of equations (3.20) exactly corresponds to the generalized Lagrange-Hermite interpolation problem, which is known to be invertible. Thus one necessarily has Q σ = 0 for any σ = 1, . . . , τ + . It then only remains to deduce that any of the polynomials p σ is also zero.
Observe that for any integer k with 0 ≤ k ≤ k b , from the divided difference algebra, the polynomial p σ [X − k b , . . . , X − k] has degree less than µ σ − (k b − k), see (3.17). Thus the highest degree polynomial involved in the sum (3.19) is p σ [X − k b ] (for k = k b ). Since we know that Q σ is zero, then p σ is also necessarily zero (consider the highest degree coefficient). The injectivity of the boundary layer problem (3.14)-(3.13) is proved, and the matrix A k b is therefore invertible.
As a consequence, there exist some uniquely determined coefficients (β σ,ν,ℓ ) that depend only on the considered scheme and on the extrapolation order k b (but neither on the initial condition f nor on the time index n), such that the solution to (3.12)-(3.13) has the form: (3.21) Using now triangular inequalities, we obtain, for some constant C > 0, the upper bound: On the one side, we recall the definition (3.8) of the sequences ρ (σ,ν) in Lemma 3.4, hence the estimate: where the constant C > 0 is independent of J and ∆x. On the other side, from [CL20, Lemma 3.6] and the continuity of the reflection operator from H k b +1 ((−∞, L)) to H k b +1 (R), we have the upper bound: and thus the required estimate (3.15) follows.
The interested reader will find in [Gol77] a similar argument to the one developed in the proof of Lemma 3.5. In [Gol77], the analysis of the determinant of the matrix A k b arises from the verification of the so-called Uniform Kreiss-Lopatinskii Condition (a condition whose significance is based on the work [GKS72]). Let us now prove Theorem 1.3. The error (ε n j ) j≤J+p,n∈N introduced in (3.11), and fully defined through Lemma 3.5, satisfies the following set of equations 4 : Considering the scheme (3.23), the error (ε n j ) j≤J+p,0≤n≤T /∆t obeys the stability estimate applicable in the case of the homogeneous extrapolation boundary condition, see [CL20, Proposition 3.4]: (3.24) It therefore remains to estimate the initial and interior consistency errors in (3.23): • The initial consistency error. Estimating the initial condition (ε 0 j ) j≤J directly follows from the estimate (3.15) in Lemma 3.5: • The interior consistency error. I. Estimating the interior consistency error (e n j ) related to the projected exact solution (ω n j ) has already been achieved in [CL20] so we just report the result: • The interior consistency error. II. Estimation of the new error term related to (δ n j ). Observe that, first due to the steady states decomposition from Lemma 3.4, and then using successively (3.12) and (3.21), the interior consistency error arising from the boundary layer corrector rewrites as: Thus, from Cauchy-Schwarz inequalities, there exists a constant C such that In the above formula, the discrete in time derivative of ω n j rewrites, for any j ≤ J + p as dx .
Since f ∈ H k+1 (R) with k > k b , we have at least f ∈ H k b +2 (R) and therefore F ∈ H k b +1 (R) with from which we deduce, using again [CL20, Lemma 3.6]: Thus, using again the upper bound (3.22), the above estimate and the H k b +2 -continuity of the extension operator, we end up with: Let us now come back to the stability estimate (3.24) and use the three above consistency estimates to get (recall T ≥ 1 and k b < k): From the constructive formula for the boundary layer corrector (v n j ), we have derived the bound (3.15) which, by the triangle inequality, yields the convergence estimate (recall ε n j = ω n j − u n j + ∆x k b v n j ): Using now the (crude) estimate: we complete the proof of Theorem 1.3.

Numerical experiments 4.1 The Lax-Wendroff scheme
We report in this paragraph on various numerical experiments with the Lax-Wendroff scheme (1.3) (which corresponds to p = r = 1). Assumption 1.1 is satisfied provided that λ a ≤ 1, and the order k equals 2. In all what follows, we choose a = 1 and λ = 5/6. The interval length is L = 6 and the final time T equals 8. The initial condition is f (x) = sin x and the boundary source term is g(t) = − sin t so that the exact solution to (1.1) is u(t, x) = sin(x − t). With the values of J reported in Table 1 below, we implement the Lax-Wendroff scheme (1.3) with the following numerical boundary conditions: u n J+1 = u n J , (first order outflow extrapolation condition), u n 0 = − sin t n , (Dirichlet inflow condition (1.6)), − sin t n − (∆x/2) cos t n , (inverse Lax-Wendroff inflow condition (1.12)).
The errors, as measured in the statement of Theorem 1.2, are reported in Table 1 below for each of the two cases (either the Dirichlet inflow condition (1.6) or the inverse Lax-Wendroff inflow condition (1.12)). In either case, the observed convergence rate is 1 since increasing J by a factor 2 decreases the error of the same factor 2. This behavior is fully justified by Theorem 1.3 since we have k b < k here.
The errors, as measured in the statement of Theorem 1.2, are reported in Table 2 below for each of the two cases (either the Dirichlet inflow condition (1.6) or the inverse Lax-Wendroff inflow condition (1.12)). For the Dirichlet inflow condition, the observed convergence rate is 1 again (despite the more accurate outflow treatment), but one recovers the convergence rate 2 with the inverse Lax-Wendroff inflow condition (1.12). However, proving rigorously that this numerical scheme converges with the rate 2 in the maximum norm might be very difficult (it might actually even be wrong !), even for smooth data, since the Lax-Wendroff scheme is known to be unstable in ℓ ∞ (Z). Improving the convergence rate 3/2 of Theorem 1.2 in the case of the Lax-Wendroff scheme with second order extrapolation outflow condition is left to a future work.
The measured errors are reported in Table 3 below. They correspond to a rate of convergence 3. Let us observe that the O3 scheme is known to be stable in ℓ ∞ (Z), see [Tho65,Des08], hence there is a genuine hope of proving rigorously that this rate of convergence does indeed hold (for smooth compatible data). Such a justification is also left to a future work.