 # 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.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

### 1.1 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 , an interval length and we consider the (continuous) problem:

 ⎧⎨⎩∂tu+a∂xu=0,t≥0,x∈(0,L),u(0,x)=f(x),x∈(0,L),u(t,0)=g(t),t≥0. (1.1)

The requirements on the initial and boundary data, namely and , will be made precise below. The solution to (1.1) is given by the method of characteristics, which yields the formula:

 ∀(t,x)∈R+×(0,L),u(t,x)={f(x−at),\rm if x≥at,g(t−xa),\rm if x≤at. (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, in (1.1)). The focus in [CL20] is on the outflow boundary ( here since 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 , that is meant to be large, and define the space step and the grid points by

 Δx:=LJ,xj:=jΔx(j∈Z).

The interval corresponds to the cells with , but considering the whole real line will be useful in some parts of the analysis. The time step is then defined as , where is a constant that is fixed so that assumption 1.1 below is satisfied. We use from now on the notation , ; the quantity will play the role of an approximation for the solution to (1.1) at the time on the cell .

We now examine an example where the exact solution to (1.1) is approximated by means of the Lax-Wendroff scheme. The approximation reads:

 un+1j=unj−λa2(unj+1−unj−1)+(λa)22(unj+1−2unj+unj−1),n∈N,j=1,…,J, (1.3)

where we recall that is a fixed constant and 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 in (1.1), namely:

 ∀j=1,…,J,u0j:=1Δx∫xjxj−1f(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 , the definition of the boundary (or ghost cell) values and . At the outflow boundary, we prescribe an extrapolation condition [Kre66, Gol77], the significance of which will be thoroughly justified in the next sections:

 unJ+1=2unJ−unJ−1,n∈N. (1.5)

Combining (1.3) with (1.5), the last interior cell value obeys the induction formula:

 un+1J=unJ−λa(unJ−unJ−1),n∈N,

which is nothing but the upwind scheme. It then only remains to determine the inflow numerical boundary condition . Since we wish to approximate the exact solution to (1.1) and is meant, at least, to approximate the trace , it seems reasonable at first sight to prescribe the Dirichlet boundary condition:

 un0=g(tn),n∈N. (1.6)

In the case of zero incoming boundary data (), and for any sufficiently smooth initial condition 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 in the maximum norm. Numerical simulations even predict that the rate of convergence should be , or at least close to , for smooth initial data. However, implementing the above numerical scheme111One can choose for instance , , , , and increase the integer geometrically. quickly shows that the rate of convergence falls down to when is nonzero and satisfies the compatibility conditions222The rate of convergence could be even smaller than when the compatibility conditions are not satisfied but that would just reflect the fact that the exact solution (1.2) is not smooth (for instance, not even continuous if ). described hereafter with the initial condition .

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 close to the incoming boundary and then using the advection equation (1.1) to substitute the normal derivatives for tangential derivatives , 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

in terms of the boundary datum but also of interior cell values with . 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].

### 1.2 The inverse Lax-Wendroff method

We first fix from now on some notations. In all this article, we are given some fixed integers and consider an explicit two time step approximation for the solution to (1.1):

 (1.7)

In (1.7), the numbers are defined in terms of the parameter and of the velocity (see, for instance, (1.3) for which ). These numbers are fixed, which means that (1.7) is linear with respect to . For simplicity, we follow [CL20] and choose as initial data for (1.7) the cell averages of the initial condition in (1.1

). This means that the vector

is defined by (1.4). For (1.7) to define inductively (with respect to ) the vector , we need to prescribe the ghost cell values and . 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 as:

 (D−v)j:=vj−vj−1,

and its iterates accordingly, we choose from now on an extrapolation order for the outflow boundary and prescribe:

 (Dkb−un)J+ℓ=0,n∈N,ℓ=1,…,p. (1.8)

The example (1.5) corresponds to (recall for the Lax-Wendroff scheme so there is only one ghost cell). It now remains to prescribe the inflow values . Unlike some previous works, we are going to prescribe Dirichlet boundary conditions, meaning for instance that the value will be determined in terms of the boundary datum only. Let us assume for a while that is a second order approximation of where is the exact solution (1.2) of the continuous problem (1.1). Then we formally have:

 un0≈u(tn,−Δx2)≈u(tn,0)−Δx2∂xu(tn,0),

where means “equal up to ”, and we then use (1.1) to get:

 un0≈u(tn,0)+Δx2a∂tu(tn,0)=g(tn)+Δx2ag′(tn).

The last term 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 should be “close” (whatever that means !) to some quantity of the form:

 K∑κ=0Δxκκ!aκακg(κ)(tn),

where is a truncation order and are numerical constants.

The general form of the Dirichlet boundary conditions that we consider below is:

 unℓ=K∑κ=0Δxκκ!aκακ,ℓg(κ)(tn),n∈N,ℓ=1−r,…,0,

where the ’s are numerical constants which will play a role (together with the truncation order ) in the consistency analysis. There are two main choices which we discuss in this article. The first one is given in [TS10, VS15]:

 ακ,ℓ:=(12−ℓ)κ,κ∈N,ℓ=1−r,…,0,

and is relevant if is eventually compared in the convergence analysis with , being the exact solution (1.2). The other possible choice we advocate is:

 ακ,ℓ:=(−1)κκ+1(ℓκ+1−(ℓ−1)κ+1),κ∈N,ℓ=1−r,…,0, (1.9)

and is relevant if is eventually compared (as in our main result, which is Theorem 1.2 below) in the convergence analysis with the average of on the cell . The truncation order is discussed with our main result in the following paragraph.

### 1.3 Results

We assume that the approximation (1.7) is consistent with the transport operator and that it defines a stable procedure on .

###### Assumption 1.1 (Consistency and stability without any boundary).

The coefficients in (1.7) satisfy (normalization), and for some integer , there holds:

 ∀m=0,…,k,p∑ℓ=−rℓmaℓ =(−λa)m,\rm(consistency of order k), (1.10) supθ∈[0,2π]∣∣ ∣∣p∑ℓ=−raℓeiℓθ∣∣ ∣∣ ≤1,\rm(ℓ2-stability on Z). (1.11)

For the Lax-Wendroff scheme (1.3), we have if , the integer equals , and (1.11) holds if and only if . In Theorem 1.2 below and all what follows, the velocity , the length , the parameter and the extrapolation order at the outflow boundary are given. Subsequent constants may depend on them. The integer 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:

 unℓ=k−1∑κ=0Δxκ(κ+1)!(−a)κ(ℓκ+1−(ℓ−1)κ+1)g(κ)(tn),n∈N,ℓ=1−r,…,0. (1.12)

Of course, prescribing (1.12) is meaningful only if is sufficiently smooth (say, ). One could push further the Taylor expansion in (1.12) and consider higher order correctors but it would require further smoothness on and it would eventually not improve our convergence result below, so fixing the truncation order 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.

###### Theorem 1.2 (Main convergence result).

Under assumption 1.1, there exists a constant such that for any final time , any integer , any data and satisfying the compatibility requirements at :

 ∀m=0,…,k,f(m)(0)=(−a)−mg(m)(0),

the solution to (1.4), (1.7), (1.8), (1.12) satisfies:

 sup0≤n≤T/Δtsup1≤j≤J∣∣∣unj−1Δx∫xjxj−1u(tn,x)dx∣∣∣≤CTeCT/LΔxmin(k,kb)−1/2(∥f∥Hk+1((0,L))+∥g∥Hk+1((0,T))), (1.13)

with the exact solution to (1.1), whose expression is given by (1.2).

Actually, the constant in (1.13) is independent of , which is consistent with the convergence result we shall prove below for the half-space problem on with inflow at . As in [CL20], the loss of in the rate of convergence of Theorem 1.2 looks somehow artificial and is mostly a matter of passing from the topology to . Our next result examines a situation where the optimal convergence rate 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 such that for any final time , any integer , any data , the solution to the scheme:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩u0j=1Δx∫xjxj−1f(x)dx,j≤J,(Dkb−un)J+ℓ=0,0≤n≤T/Δt,ℓ=1,…,p,un+1j=p∑ℓ=−raℓunj+ℓ,0≤n≤T/Δt−1,j≤J, (1.14)

satisfies the error estimate

 sup0≤n≤T/Δtsupj≤J∣∣∣unj−1Δx∫xjxj−1f(x−atn)dx∣∣∣≤CTΔxkb∥f∥Hk+1((0,L)),

as long as .

In other words, the technical assumption 3.2 hereafter, which is verified on many examples such as the Lax-Wendroff and schemes, allows to recover the optimal rate in the case . Of course, one would also like to improve the rate in the case , which is clearly the most natural choice. However, in that case, both the interior and boundary consistency errors scale like and, in the framework of assumption 1.1, stability in the interior domain is available only in the topology, so it is quite difficult to derive the convergence rate in the 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 in the future.

## 2 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 , and one for which there is only outflow at . The novelty here is the nonzero inflow source term so we first deal with that case.

### 2.1 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 such that for any final time , for any , for any initial condition and boundary source term satisfying the compatibility conditions:

 ∀m=0,…,k,f(m)(0)=(−a)−mg(m)(0), (2.1)

the solution to the numerical scheme:

 (2.2)

satisfies:

 sup0≤n≤T/Δt⎛⎝∑j≥1Δx(unj−1Δx∫xjxj−1u(tn,x)dx)2⎞⎠1/2≤CTΔxk(∥f∥Hk+1((0,+∞))+∥g∥Hk+1((0,T))),

where is the exact solution to the half-line transport problem:

 ⎧⎨⎩∂tu+a∂xu=0,t∈(0,T),x≥0,u(0,x)=f(x),x≥0,u(t,0)=g(t),t∈(0,T). (2.3)
###### Proof.

For convenience, we first extend into a function and then define:

 ∀x∈R,f♯(x):={f(x),% \rm if x>0,g♭(−x/a),\rm if x<0.

Since and satisfy the compatibility conditions (2.1), we have , and the exact solution to (2.3) is given by:

 ∀(t,x)∈[0,T]×(0,+∞),u(t,x)=f♯(x−at).

Let us now define:

 ∀j∈Z,∀n∈N,wnj:=1Δx∫xjxj−1f♯(x−atn)dx,

which corresponds to the cell average of the exact solution to (2.3). With the solution to the numerical scheme (2.2), we define the error , that is a solution to:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ε0j=0,j≥1,εnℓ=unℓ−wnℓ,0≤n≤T/Δt,ℓ=1−r,…,0,εn+1j=p∑ℓ=−raℓεnj+ℓ+Δten+1j,0≤n≤T/Δt−1,j≥1. (2.4)

The interior consistency error is easily estimated by means of the Cauchy-Schwarz inequality and Fourier analysis:

 ∑j≥1Δx(en+1j)2= ΔxΔt2∑j≥1(wn+1j−p∑ℓ=−raℓwnj+ℓ)2 = 1ΔxΔt2∑j≥1(∫xjxj−1(f♯(x−atn−aΔt)−p∑ℓ=−raℓf♯(x−atn+ℓΔx))dx)2 ≤ 1Δt2∫R(f♯(x−atn−aΔt)−p∑ℓ=−raℓf♯(x−atn+ℓΔx))2dx ≤ 12πΔt2∫R∣∣e−iaλΔxξ−p∑ℓ=−raℓeiℓΔxξ∣∣2∣∣ˆf♯(ξ)∣∣2dξ≤CΔx2k∫Rξ2(k+1)∣∣ˆf♯(ξ)∣∣2dξ,

where the final inequality comes from assumption 1.1 and the fact that the ratio is constant. Going back to the definition of , we have obtained the bound:

 sup0≤n≤T/Δt−1(∑j≥1Δx(en+1j)2)1/2≤CΔxk(∥f∥Hk+1((0,+∞))+∥g∥Hk+1((0,T))), (2.5)

for some constant that is independent of the final time and the data and .

We now turn to the boundary errors in (2.4), and wish to estimate the following quantities:

Let us consider an integer such that . From the definition of , , we have:

 unℓ−wnℓ= k−1∑κ=0Δxκ(κ+1)!(−a)κ(ℓκ+1−(ℓ−1)κ+1)g(κ)(tn)−1Δx∫xℓxℓ−1g♭(tn−x/a)dx = −(−a)−kΔx∫xℓxℓ−1xk∫10yk−1(k−1)!g(k)♭(tn−xya)dydx,

where we have used the Taylor formula333This is precisely at this point of the analysis that the definition of the coefficients in the inverse Lax-Wendroff method arises. Our choice in (1.12) is motivated by the fact that we compare the numerical solution with the cell average of the exact solution.. By the Cauchy-Schwarz inequality, we get:

 (unℓ−wnℓ)2≤CΔx∫xℓxℓ−1∫10x2ky2(k−1)g(k)♭(tn−xya)2dydx,

and we now apply the change of variables to get:

 (unℓ−wnℓ)2≤∫xℓxℓ−1(∫0v|v||u|2(k−1)g(k)♭(tn−ua)2du)dv.

Restricting to , we have:

 0∑ℓ=1−r(unℓ−wnℓ)2≤CΔx2k−1∫rΔx/a0g(k)♭(tn+τ)2dτ.

Summing now with respect to , we end up with the estimate:

 ∑0≤n≤T/Δt−1Δt0∑ℓ=1−r(εnℓ)2≤CΔx2k∫+∞0g(k)♭(t)2dt≤CΔx2k∥g∥2Hk((0,T)), (2.6)

for some constant that is independent of the final time and the data and .

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]:

 sup0≤n≤T/Δt(∑j≥1Δx(εnj)2)1/2≤C⎧⎨⎩Tsup1≤n≤T/Δt(∑j≥1Δx(enj)2)1/2+⎛⎝∑0≤n≤T/Δt−1Δt0∑ℓ=1−r(εnℓ)2⎞⎠1/2⎫⎪⎬⎪⎭.

The conclusion of Theorem 2.1 then comes from the combination of the estimates (2.5) and (2.6). ∎

### 2.2 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 such that for any final time , for any and for any initial condition , the solution to the numerical scheme:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩u0j=1Δx∫xjxj−1f(x)dx,j≤J,(Dkb−un)J+ℓ=0,0≤n≤T/Δt,ℓ=1,…,p,un+1j=p∑ℓ=−raℓunj+ℓ,0≤n≤T/Δt−1,j≤J,

satisfies:

 sup0≤n≤T/Δt⎛⎝∑j≤JΔx(unj−1Δx∫xjxj−1f(x−atn)dx)2⎞⎠1/2≤CTΔxmin(k,kb)∥f∥Hk+1((−∞,L)).

### 2.3 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 and some data , that satisfy the compatibility conditions stated in Theorem 1.2. We consider some function such that if and if . We then decompose the initial condition as:

 ∀x∈(0,L),f(x)=(1−χ(x/L))f(x)+χ(x/L)f(x).

Since vanishes on , we can extend it by zero to the interval and thus consider as an element of . Furthermore, the functions and satisfy the same compatibility conditions as and at . We can thus apply Theorem 2.1 to the sequence that is defined as the solution to the numerical scheme:

We obtain the estimate:

 sup0≤n≤T/Δt⎛⎝∑j≥1Δx(vnj−1Δx∫xjxj−1v(tn,x)dx)2⎞⎠1/2≤CTΔxk(∥f∥Hk+1((0,L))+∥g∥Hk+1((0,T))), (2.7)

where is the exact solution to the transport problem:

 ⎧⎨⎩∂tv+a∂xv=0,t∈(0,T),x≥0,v(0,x)=(1−χ(x/L))f(x),x≥0,v(t,0)=g(t),t∈(0,T).

Similarly, we can view as an element of that vanishes on . Theorem 2.2 then shows that the solution to the numerical scheme:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩w0j=1Δx∫xjxj−1χ(x/L)f(x)dx,j≤J,(Dkb−wn)J+ℓ=0,0≤n≤T/Δt,ℓ=1,…,p,wn+1j=p∑ℓ=−raℓwnj+ℓ,0≤n≤T/Δt−1,j≤J,

satisfies:

 sup0≤n≤T/Δt⎛⎝∑j≤JΔx(wnj−1Δx∫xjxj−1χ((x−atn)/L)f(x−atn)dx)2⎞⎠1/2≤CTΔxmin(k,kb)∥f∥Hk+1((0,L)). (2.8)

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 up to the threshold:

 N:=min(E(J/3−kbr),E(E(J/3)p)),

there holds:

 wn1−r=⋯=wn0=0,vnJ+1−kb=⋯=vnJ+p=0.

In particular, the solution to (1.7), (1.4), (1.8), (1.12) satisfies:

 ∀n=0,…,N,∀j=1−r,…,J+p,unj=vnj+wnj.

Combining then the error estimates (2.7) and (2.8), we obtain:

 (2.9)

where 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 between and , we split the solution to (1.7), (1.4), (1.8), (1.12) as the sum of the solution to the problem:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩~uNj=1Δx∫xjxj−1u(tN,x)dx,1≤j≤J,(Dkb−~uN+n)J+ℓ=0,0≤n≤N,ℓ=1,…,p,~uN+nℓ=k−1∑κ=0Δxκ(κ+1)!(−a)κ(ℓκ+1−(ℓ−1)κ+1)g(κ)(tN+n),0≤n≤N,ℓ=1−r,…,0,~uN+n+1j=p∑ℓ=−raℓ~uN+nj+ℓ,0≤n≤N−1,1≤j≤J,

and of the (presumably small) solution to the ‘error’ problem: