Delphine Bresch-Pietri

Miroslav Krstic

Robustness to Diffusion of Prediction-based Control for Convection Processes

In this paper, we consider an Ordinary Differential Equation (ODE) with convection and diffusion in the actuation path. We prove that a prediction-based controller, designed to compensate for the sole convective PDE, actually achieves exponential stabilization of the complete plant, provided that diffusion is small enough. Our result is obtained in L p norm and covers two cases, full-state feedback and boundary feedback. Simulation results emphasize the validity of this approach.

I. INTRODUCTION

For linear systems subject to input-delay, prediction-based control strategies, more commonly known as Smith Predictor (see [START_REF] Artstein | Linear systems with delayed controls: a reduction[END_REF], [START_REF] Manitius | Finite spectrum assignment problem for systems with delays[END_REF], [START_REF] Smith | A controller to overcome dead time[END_REF],) are state-of-the-art for systems with constant input time-delays (see for instance [START_REF] Bekiaris-Liberis | Nonlinear Control Under Nonconstant Delays[END_REF], [START_REF] Bresch-Pietri | Adaptive control scheme for uncertain time-delay systems[END_REF], [START_REF] Gu | Survey on recent results in the stability and control of time-delay systems[END_REF], [START_REF] Jankovic | Recursive predictor design for linear systems with time delay[END_REF], [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems[END_REF], or [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF] and the references therein). Grounding on the use of a prediction of the system state on a time horizon equal to the delay, this technique aims at compensating it, which notably improves transient performances.

When the actuation path is subject to both convection and diffusion, this technique has been recently extended in [START_REF] Krstic | Compensating actuator and sensor dynamics governed by diffusion PDEs[END_REF], [START_REF] Susto | Control of PDE-ODE cascades with neumann interconnections[END_REF] (see also [START_REF] Bresch-Pietri | Adaptive compensation of diffusionadvection actuator dynamics using boundary measurements[END_REF] for an adaptive version) to compensate for both effects. However, while delay dynamics can be finitely stabilized, diffusion compensation should then be understood differently as it introduces an infinite relative degree. Correspondingly, while delay compensation requires to predict future values of the ODE state, in the diffusion case, the control law mainly accounts for an inversion of the diffusion dynamics.

When the diffusive effect can be neglected with respect to the convective one, one may be interested in using a standard prediction-based control law to compensate only for the dominant delay dynamics and, e.g., decrease the controller complexity. However, a well-known fact about predictionbased techniques [START_REF] Palmor | Stability properties of Smith dead-time compensator controllers[END_REF] is that they suffer from being sensitive to various disturbances such as delay mismatch or, to a lesser extent, plant parameters uncertainties [START_REF] Mondié | Delay robustness of closed loop finite assignment for input delay systems[END_REF][17] [START_REF] Zhong | Robust Control of Time-Delay Systems[END_REF]. Thus, there is no guarantee a priori that a prediction-based feedback is robust to diffusion in the actuation path.

For this reason, in this paper, we consider a convection/diffusion PDE-ODE cascade coupled with a predictionbased controller. To analyze the resulting closed-loop robustness, we follow the distributed parameters systems methodology introduced in [START_REF] Krstic | Boundary Control of PDEs: a Course on Backstepping Designs[END_REF] which has proved to be effective to conclude on the robustness of prediction-based controller D. Bresch-Pietri (corresponding author) is with CNRS at the Department of Automatic Control, Gipsa-lab, 11 rue des Mathématiques, 38000 Grenoble, France Email: delphine.bresch-pietri@gipsa-lab.fr M. Krstic is with the Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla CA 92093, USA to several classes of uncertainties [START_REF] Bresch-Pietri | Adaptive trajectory tracking despite unknown input delay and plant parameters[END_REF], [START_REF] Bresch-Pietri | Adaptive control scheme for uncertain time-delay systems[END_REF], [START_REF] Bresch-Pietri | Delay-adaptive control for nonlinear systems[END_REF]. In details, focusing on the sole convection effect in the actuation path, we propose to apply to the distributed input a backstepping transformation which has been developed for input-delay systems. This allows us to map the original system into a target system which is suitable for Lyapunov stability analysis. We investigate two distinct frameworks, full-state feedback and boundary feedback, and show that, in both cases, exponential stabilization is obtained in the sense of the L p norm, provided that the diffusion parameter is small enough. This is the main contribution.

The paper is organized as follows. In Section II, we start by presenting the problem under consideration. Then, in Section III and IV, we tackle the cases of full-state and boundary control respectively and present our robustness results along with their proofs. This robustness property is illustrated with a numerical example in Section V. We conclude by sketching directions of future work.

II. PROBLEM STATEMENT AND CONTROL DESIGN

Consider the following ODE-PDE cascade

         Ẋ(t) =AX(t) + Bu(0,t) u t (x,t) =u x (x,t) + εu xx (x,t) u(1,t) =U(t) εu x (0,t) =0 (1) in which X ∈ R n , the control U is scalar, u(•, 0) ∈ C (R, H 1 ([0, 1]
)) and ε ≥ 0 is a small parameter. We consider in the following that the ODE state X is measured and that the pair (A, B) is controllable. In the sequel, we will consider K such that A + BK is Hurwitz.

When ε = 0, system (1) is an input-delay system, with a unitary delay. It is well-known that the following control law compensates for the delay

U(t) =K e A X(t) + t t-1 e A(t-s) BU(s)ds =K e A X(t) + 1 0 e A(1-x) Bu(x,t)dx (2)
namely that, after one unit of time, the control law (2) yields the delay-free dynamics Ẋ(t) = (A + BK)X(t) which is exponentially stable. It is also worth noticing that, when ε = 0, one can solve explicitly the transport PDE and obtain u(x,t) = U(t + x -1) which explains the second expression provided in [START_REF] Bekiaris-Liberis | Nonlinear Control Under Nonconstant Delays[END_REF].

In this paper, we aim at studying the stability of the closedloop system (1)-( 2) in presence of diffusion, that is, for ε ≥ 0. As, when ε = 0, the two expressions provided in [START_REF] Bekiaris-Liberis | Nonlinear Control Under Nonconstant Delays[END_REF] are not equivalent, we consider them separately in the following.

III. FULL-STATE FEEDBACK

In this section, we consider that u(•,t) is measured. This corresponds to a dynamical system where convection/diffusion arises from a physical transport process, as opposed to computational or communication-based delay. Consequently, we select the control law as

U(t) =K e A X(t) + 1 0 e A(1-x) Bu(x,t)dx (3)
We have the following robustness result. Theorem 1: Consider the closed-loop system consisting of the plant (1) and the control law [START_REF] Bekiaris-Liberis | Predictor-feedback stabilization of multi-input nonlinear systems[END_REF]. Define the following functional

Γ p (t) = |X(t)| 2p + 1 0 u(x,t) 2p dx (4)
for a given p ∈ N * . There exists ε p > 0 such that, if ε < ε p , then there exist ρ > 0 (independent of p) and R p > 0 such that

Γ p (t) ≤ R p Γ p (0)e -2pρt , t ≥ 0 (5)
This result states that the desired exponential convergence is achieved, provided that diffusion is small enough. It is worth noticing that this result holds for any p ∈ N * , but that the corresponding diffusion coefficient bound ε p decreases with p. Thus, one cannot obtain exponential stability in the infinity norm (see Section III-B for more details on this point).

We now provide the proof of this result.

A. Proof of Theorem 1 1) Backstepping transformation: We consider the following backstepping transformation of the distributed state

w(x,t) =u(x,t) - x 0 Ke A(x-y) Bu(y,t)dy -Ke Ax X(t) (6)
Lemma 1: The backstepping transformation (6), along with the control law (3), transform the plant (1) into the following target system

         Ẋ =(A + BK)X(t) + Bw(0,t) w t (x,t) =w x (x,t) + εw xx (x,t) + εg(x,t) w(1,t) =0 εw x (0,t) = -ε(KBw(0,t) + K(A + BK)X(t)) (7) in which g(x,t) =Ke Ax A((A + BK)X(t) + Bw(0,t)) (8)
Proof: Taking a time-derivative of ( 6) and using an integration by parts, one obtains

w t =u t -KBu(x,t) -K x 0 e A(x-y) (ABu + εBu xx )dy -Ke Ax AX (9)
Similarly, one obtains for the first and second-order spatial derivatives

w x =u x -KBu(x,t) -K x 0
Ae A(x-y) Bu(y,t)dy -Ke Ax AX (10)

w xx =u xx -KBu x (x,t) -KABu(x,t) -K x 0 A 2 e A(x-y) Bu(y,t)dy -Ke Ax A 2 X (11)
Thus, it follows that

w t =w x + εw xx + ε -K x 0 e A(x-y) Bu xx dy + KBu x (x,t) (12) + KABu(x,t) + K x 0 A 2 e A(x-y) Bu(y,t)dy + Ke Ax A 2 X
which gives the expected result with an integration by parts. The boundary conditions in [START_REF] Bresch-Pietri | Delay-adaptive control for nonlinear systems[END_REF] follow straightforwardly respectively from the definition of the control law (3) and the backstepping transformation [START_REF] Bresch-Pietri | Adaptive trajectory tracking despite unknown input delay and plant parameters[END_REF], and from equation ( 10) evaluated at x = 0.

2) Stability study in L p -norm: Define the following Lyapunov functional candidate

V p (t) =(X T PX) p + b 1 0 e 2px w(x,t) 2p dx (13)
in which b is a positive constant, p ∈ N * and P is the solution to the Lyapunov equation associated to A + BK and a given positive definite matrix Q. For the sake of conciseness, in the sequel, we denote w(t) 2p 2p = 1 0 e 2px w(x,t) 2p dx. Taking a time-derivative of ( 13), one gets

Vp (t) = -p(X T PX) p-1 X T QX + 2p(X T PX) p-1 X T PBw(0,t) -bw(0,t) 2p -2pb w(t) 2p 2p + 2pbε 1 0 e 2px w 2p-1 w xx dx + 2pbε 1 0 e 2px g(x,t)w 2p-1 (x,t)dx (14)
Performing an integration by parts, one obtains

Vp (t) = -p(X T PX) p-1 X T QX + 2p(X T PX) p-1 X T PBw(0,t) -bw(0,t) 2p -2pb w(t) 2p 2p + εb 2pw(0,t) 2p -2pw(0,t) 2p-1 w x (0,t) + (2p) 2 w(t) 2p 2p - 1 0 2p(2p -1)e 2px w(x,t) 2p-2 w x (x,t) 2 dx + 2pbε 1 0 e 2px g(x,t)w(x,t) 2p-1 dx (15)
Further, with an integration by parts and applying Young inequality, one can observe that 2p Consequently, one gets

Vp (t) ≤ -pλ (Q)(X T PX) p-1 |X| 2 + 2p(X T PX) p-1 |X||PBw(0,t)| -bw(0,t) 2p -2pb w(t) 2p 2p + εb 2pw(0,t) 2p -2pw(0,t) 2p-1 w x (0,t) -2(2p -1)w(0,t) 2p -2p(2p -1) w(t) 2p 2p + (2p) 2 w(t) 2p 2p + 2pbε 1 0 e 2px g(x,t)w(x,t) 2p-1 dx (18)
From the definition of the function g in ( 8) and applying Young inequality, one obtains the existence of a positive constant M > 0 (which does not depend on p) such that

2p 1 0 e 2px g(x,t)w(x,t) 2p-1 dx ≤M 2p w(t) 2p 2p + w(0,t) 2p + |X(t)| 2p (19)
Similarly, one gets

2p(X T PX) p-1 |X||PBw(0,t)| ≤(2p -1) λ (Q) 4
λ (P) λ (P) (X T PX) p-1 |X(t)| 2p 2p-1 + 4 λ (Q) λ (P) λ (P) 2p-1 |PB| 2p w(0,t) 2p ≤ pλ (Q) 2 (X T PX) p-1 |X(t)| 2 + 4 λ (Q) λ (P) λ (P) 2p-1 |PB| 2p w(0,t) 2p (20)
in which λ (•) and λ (•) denote respectively the minimum and maximum eigenvalues of a given matrix, and

-2pw(0,t) 2p-1 w x (0,t) ≤2pKBw(0,t) 2p + (2p -1)w(0,t) 2p + |K(A + BK)| 2p |X(t)| 2p (21)
Plugging these inequalities into [START_REF] Palmor | Stability properties of Smith dead-time compensator controllers[END_REF], one gets

Vp (t) ≤ - pλ (Q) 4 (X T PX) p-1 |X| 2 - pλ (Q)λ (P) p-1 4 -bε(M + |K(A + BK)| 2p ) |X| 2p -2pb(1 -ε(M + 1)) w(t) 2p 2p -b - 4 λ (Q) λ (P) λ (P) 2p-1 |PB| 2p -bε(1 + M + 2p|KB|) w(0,t) 2p (22)
Assume that b >2 4λ (P)

λ (Q)λ (P) 2p-1 |PB| 2p (23)
and that ε < ε p with

ε p = min pλ (Q)λ (P) p-1 4b(M + |K(A + BK)| 2p ) , 1 2(M + 1 + 2p|KB|) (24)
then, one obtains

Vp (t) ≤ -2p min λ (Q) 8λ (P) , 1 2 V p (t) (25)
which gives

V p (t) ≤ e -2p min λ (Q) 8λ (P) , 1 2 t V p (0) (26)
3) Equivalence between V p and Γ p : To conclude the proof, it now remains to show that V p and Γ p are equivalent, i.e., that there exist c 1 , c 2 > 0 such that c 1 V p ≤ Γ p ≤ c 2 V p . This fact actually follows straightforwardly by applying Young inequality to (6) and its inverse u(x,t) =w(x,t) +

x 0 Ke (A+BK)(x-y) Bw(y,t)dy + Ke (A+BK)x X(t) (27)
However, note that these constants c 1 and c 2 depend on p. Then, it follows from (26) that ( 5) holds with R p = c 2 c 1 and ρ = min λ (Q) 8λ (P) , 1 2 .

B. Comments on the norm selection

It is worth noting that, even if our robustness result holds for all L p norm, it does not extend to the infinity norm. Indeed, the bound ε p defined in (24) depends on the value of p and tends to zero as p tends to infinity. Thus, according to our analysis at least, prediction-based control is not robust to diffusion in infinity norm. Interestingly, the same conclusion can be obtained by studying the system via Input-to-State Stability and attempting to apply the Small Gain Theorem.

IV. BOUNDARY FEEDBACK

In this section, as the distributed input is seldom measured in application, we consider the more realistic case in which the actuator state is unknown. Consequently, we choose the control law as

U(t) =K e A X(t) + t t-1 e A(t-s) BU(s)ds (28)
In this case, a similar result to the one provided in the previous section holds. Theorem 2: Consider the closed-loop system consisting of the plant (1) and the control law (28). Define the following functional

ϒ p (t) =|X(t)| 2p + t t-1 U(s) 2p ds ( 29
)
+ t t-1 U(s) 2p ds + t t-1 Ü(s) 2p ds + 1 0 u(x,t) 2p dx
for a given p ∈ N. For any initial condition

U 0 : s ∈ [-1, 0] → U(s) ∈ C 2 ([-1, 0], R)
, there exists ε p > 0 such that, if ε < ε p , then there exist ρ > 0 (independent of p) and R p > 0 such that

ϒ p (t) ≤ R p ϒ p (0)e -2pρt , t ≥ 0 (30)
The definition (29) of the functional ϒ is due to the specific form of Lyapunov functional considered in the following proof.

Proof of Theorem 2:

To account for the fact that the control law explicitly neglects diffusion in the actuator path, we define the distributed estimate û(x,t) = U(t + x -1), x ∈ [0, 1]. Thus, the controller (28) can be reformulated as

U(t) =K e A X(t) + 1 0 e A(1-x) B û(x,t)dx (31)
and, introducing the estimation error ũ = uû, one obtains

                   Ẋ(t) =AX(t) + B[ ũ(0,t) + û(0,t)] ũt (x,t) = ũx (x,t) + ε( ũxx (x,t) + ûxx (x,t)) ũ(1,t) =0 ε ũx (0,t) = -ε ûx (0,t) ût (x,t) = ûx (x,t) û(1,t) =U(t) (32)
We now consider the following backstepping transformation ŵ(x,t) = û(x,t) -x 0 Ke A(x-y) B û(y,t)dy -Ke Ax X(t) (33)

Following the same steps as in the proof of Lemma 1, one can prove that this transformation, along with the control law (31), map the previous system into

                   Ẋ(t) =(A + BK)X(t) + B[ ũ(0,t) + ŵ(0,t)] ũt (x,t) = ũx (x,t) + ε( ũxx (x,t) + h(x,t)) ũ(1,t) =0 ε ũx (0,t) = -ε[ ŵx (0,t) + KB ŵ(0,t) + K(A + BK)X(t)] ŵt (x,t) = ŵx (x,t) ŵ(1,t) =0 (34)
in which, considering the inverse backstepping transformation (similar to the one defined in the previous section in (27))

h(x,t) = ŵxx (x,t) + KB ŵx (x,t) + K(A + BK)B ŵ(x,t) + K x 0 (A + BK) 2 e (A+BK)(x-y) B ŵ(y,t)dy + Ke (A+BK)x (A + BK) 2 X(t) (35)
In the following analysis, we also need the dynamics of the first and second spatial-derivatives of the backstepping variable which is simply

ŵxt (x,t) = ŵxx (x,t) ŵx (1,t) =0 ŵxxt (x,t) = ŵxxx (x,t) ŵxx (1,t) =0 (36)
We are now ready to start the Lyapunov analysis. We consider the following functional candidate

V p (t) =(X T PX) p + b 1 0 e 2px ũ(x,t) 2p dx + b 1 0 e 2px ( ŵ(x,t) 2p + ŵx (x,t) 2p + ŵxx (x,t) 2p )dx (37)
in which b > 0, p ∈ N * and P is the solution to the Lyapunov equation associated to A +BK and a given positive definite matrix Q. Taking a time-derivative and performing integrations by parts, one obtains

Vp (t) = -p(X T PX) p-1 X T QX -b ũ(0,t) 2p -2pb ũ(t) 2p 2p + 2p(X T PX) p-1 X T PB[ ũ(0,t) + ŵ(0,t)] -b( ŵ(0,t) 2p + ŵx (0,t) 2p + ŵxx (0,t) 2p ) -2pb( ŵ(t) 2p 2p + ŵx (t) 2p 2p + ŵxx (t) 2p 2p ) + 2pbε 1 0 e 2px ũ(x,t) 2p-1 ( ũxx + h(x,t))dx (38)
in which, using integrations by parts,

1 0 e 2px ũ(x,t) 2p-1 ũxx dx = -ũ(0,t) 2p-1 ũx (0,t) + ũ(0,t) 2p + 2p ũ(t) 2p 2p -(2p -1) 1 0 e 2px ũ2p-2 ũ2 x dx (39)
Following similar steps as in ( 16)-( 17), one gets 2p 1 0 e 2px ũ(x,t) 2p-1 ũxx dx ≤ -2p ũ(0,t) 2p-1 ũx (0,t) + ũ(0,t) 2p + 2p ũ(t) Plugging together (38), ( 40) and (41), one obtains, using Young inequality,

Vp (t) ≤ -p λ (Q) 4 (X T PX) p-1 |X| 2 -b - 8λ (P) λ (Q)λ (P) 2p-1 |PB| 2p -bε|KB| ŵ(0,t) 2p -p λ (Q)λ (P) p-1 4 -εb( M + |K(A + BK)| 2p ) |X(t)| 2p -b - 8λ (P) λ (Q)λ (P) 2p-1 |PB| 2p -bε(3 + |KB|) ũ(0,t) 2p -b(1 -ε) ŵx (0,t) 2 -2pb(1 -ε(1 + M)) ũ(t) 2p 2p -2pb(1 -ε)( ŵ(t) 2p 2p + ŵx (t) 2p 2p + ŵxx (t) 2p 2p ) (42)
Consequently, choosing b > 2 8λ (P)

λ (Q)λ (P) 2p-1 |PB| 2p (43)
and if ε < ε p with

ε p = min 1 2(3 + |KB|) , pλ (Q)λ (P) p-1 4b( M + |K(A + BK)| 2p ) (44)
one obtains that

Vp (t) ≤ -2p min λ (Q) 8λ (P) , 1 2 V p (t) (45)
It now remains to show that ϒ p and V p are equivalent. This can be done as previously by considering the inverse backstepping transformation and applying Young inequality.

It is also necessary to notice that ûx (x,t) = U(t + x -1) and ûxx (x,t) = Ü(t + x -1) which explains the presence of integrals of U and Ü in (29). This concludes the proof.

The same comments as in the previous section regarding the infinity norm apply here.

V. NUMERICAL SIMULATIONS

To illustrate the merits of the proposed controller, we consider a numerical example inspired from the fluid flow system considered, e.g., in [START_REF] Bresch-Pietri | Prediction-based control of moisture in a convective flow[END_REF] [START_REF] Hernandez Perez | A time-delay approach for modeling and control of mist in a poiseuille flow[END_REF]. It consists in an horizontal tube equipped with inlet and outlet fans and a mist injector. The mist injection is fixed and the control objective is to regulate the moisture at the output of the tube, that is, the ODE state X is the moisture while the control U is the fan power. A schematic view of the set-up is given in Fig. 1.

As studied in [START_REF] Bresch-Pietri | Prediction-based control of moisture in a convective flow[END_REF] [10], a stable first-order delay equation encapsulates the main features of the dynamics under consideration. A more detailed modeling enables to conclude that diffusion also occurs inside the tube, but is negligible compared to convection. The moisture measured at the output thus satisfies (1) with A = -1/τ 0 and B = K P /τ 0 and A + BK is Hurwitz for any constant K < 0. We chose K = -10 in simulation. One can easily observe1 that the results presented in this paper straightforwardly extend to the case of a non-unitary delay D. Thus, for simulation, we take2 D = 6s, τ 0 = 9s, ε = 0.04s -1 . This choice leads to a Peclet number of εD = 0.24. As only boundary measurements are available in this example, we implement the controller (28).

We compare the obtained closed-loop performances with the full-state diffusion/convection compensation solution proposed in [START_REF] Krstic | Compensating actuator and sensor dynamics governed by diffusion PDEs[END_REF], [START_REF] Susto | Control of PDE-ODE cascades with neumann interconnections[END_REF] as a reference. Corresponding results are depicted in Fig. 2. First, one can observe that convergence is indeed obtained, even if the Peclet number is quite high in this example. Also, one can notice that, when using the proposed prediction-based controller, transient performances are slightly degraded as an overshoot is generated. However, overall, this controller reasonably performs compared to a strategy tailored to handle diffusion and using the distributed state knowledge.

VI. CONCLUSION

In this paper, we proved that prediction-based controller for linear systems subject to a single input-delay are robust to diffusion in the actuation path. Future works should investigate the extension of this property to systems subject to multiple input-delays [START_REF] Bekiaris-Liberis | Predictor-feedback stabilization of multi-input nonlinear systems[END_REF], [START_REF] Tsubakino | Predictor-feedback for multi-input lti systems with distinct delays[END_REF].

Fig. 1 .

 1 Fig. 1. Schematic view of the flowtube system under consideration as an illustrative example.

  Evolution of the moisture, the fan power and the parameter estimate.

Fig. 2 .

 2 Fig. 2. Numerical results for the closed-loop dynamics of the moisture modeled by (1) with A = -1/τ 0 and B(θ ) = K P /τ 0 , for two different controllers. The initial condition corresponds to a moisture (relative humidity) equilibrium of 50% that one wants to regulate to a new operating condition corresponding to 60% of moisture. The prediction-based controller (28) has been implemented with a feedback gain K = -10. The nominal diffusion/convection controller from[START_REF] Krstic | Compensating actuator and sensor dynamics governed by diffusion PDEs[END_REF],[START_REF] Susto | Control of PDE-ODE cascades with neumann interconnections[END_REF] has been implemented with a similar feedback gain K = -10.

Performing, from an unnormalized process Dv t = v x + av xx , a simple change of variable u(x,t) = v(x, Dt).

Note that the convection speed is actually time-varying, but we neglect those variations here for the sake of simplicity


Christophe Chalons A
email: christophe.chalons@uvsq.fr.

AND Paola Goatin B
email: paola.goatin@inria.fr.

Luis M Villada
email: lvillada@ubiobio.cl.

HIGH ORDER NUMERICAL SCHEMES FOR ONE-DIMENSION NON-LOCAL CONSERVATION LAWS

Keywords: 2000 Mathematics Subject Classification. 35L65, 90B20, 76T20, 65M08, 65M60 Scalar conservation laws, non-local flux, sedimentation models, traffic flow model, Galerkin discontinuous schemes, finite volume schemes

This paper focuses on the numerical approximation of the solutions of non-local conservation laws in one space dimension. These equations are motivated by two distinct applications, namely a traffic flow model in which the mean velocity depends on a weighted mean of the downstream traffic density, and a sedimentation model where either the solid phase velocity or the solid-fluid relative velocity depends on the concentration in a neighborhood. In both models, the velocity is a function of a convolution product between the unknown and a kernel function with compact support. It turns out that the solutions of such equations may exhibit oscillations that are very difficult to approximate using classical first-order numerical schemes. We propose to design Discontinuous Galerkin (DG) schemes and Finite Volume WENO (FV-WENO) schemes to obtain high-order approximations. As we will see, the DG schemes give the best numerical results but their CFL condition is very restrictive. On the contrary, FV-WENO schemes can be used with larger time steps. We will see that the evaluation of the convolution terms necessitates the use of quadratic polynomials reconstructions in each cell in order to obtain the high-order accuracy with the FV-WENO approach. Simulations using DG and FV-WENO schemes are presented for both applications.

Introduction

This paper is concerned with the design of numerical schemes for the one-dimensional Cauchy problem for non-local scalar conservation law of the form ρ t + (f (ρ)V (ρ * ω η )) x = 0, x ∈ R, t > 0, ρ(x, 0) = ρ 0 (x)

x ∈ R,

where the unknown density ρ depends on the space variable x and the time variable t, ρ → V (ρ) is a given velocity function, and ρ → g(ρ) = f (ρ)V (ρ) is the usual flux function for the corresponding local scalar conservation law. Here, (1.1) is non-local in the sense that the velocity function V is evaluated on a "neighborhood" of x ∈ R defined by the convolution of the density ρ and a given kernel function ω η with compact support. In this paper, we are especially interested in two specific forms of (1.1) which naturally arise in traffic flow modelling [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF] and sedimentation problems [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF]. They are given as follows.

A non-local traffic flow model. In this context, we follow [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF] and consider (1.1) as an extension of the classical Lighthill-Whitham-Richards traffic flow model, in which the mean velocity is assumed to be a non-increasing function of the downstream traffic density and where the flux function is given by

f (ρ) = ρ, V (ρ) = 1 -ρ, (ρ * ω η )(x, t) = x+η x
ω η (y -x)ρ(y, t)dy.

(1.2)

Four different nonnegative kernel functions will be considered in the numerical section, namely ω η (x) = 1/η, ω η (x) = 2(η -x)/η 2 , ω η (x) = 3(η 2 -x 2 )/(2η 3 ) and ω η (x) = 2x/η 2 with support on [0, η] for a given value of the real number η > 0. Notice that the well-posedness of this model together with the design of a first and a second order FV approximation have been considered in [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF].

A non-local sedimentation model. Following [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF], we consider under idealized assumptions that equation (1.1) represents a one-dimensional model for the sedimentation of small equal-sized spherical solid particles dispersed in a viscous fluid, where the local solids column fraction ρ = ρ(x, t) is a function of depth x and time t. In this context, the flux function is given by

f (ρ) = ρ(1 -ρ) α , V (ρ) = (1 -ρ) n , (ρ * ω η )(x, t) = 2η -2η
ω η (y)ρ(x + y, t)dy,

where n ≥ 1 and the parameter α satisfies α = 0 or α ≥ 1. The function V is the so-called hindered settling factor and the convolution term ρ * ω η is defined by a symmetric, nonnegative, and piecewise smooth kernel function ω η with support on [-2η, 2η] for a parameter η > 0 and R ω η (x)dx = 1. More precisely, the authors define in [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF] a truncated parabola K by

K(x) = 3 8 1 - x 2 4
for |x| < 2, K(x) = 0 otherwise, and set ω η (x) := η -1 K(η -1 x).

(1.4)

Conservation laws with non-local terms arise in a variety of physical and engineering applications: besides the above cited ones, we mention models for granular flows [START_REF] Amadori | An integro-differential conservation law arising in a model of granular flow[END_REF], production networks [START_REF] Gugat | Analysis of a system of nonlocal conservation laws for multi-commodity flow on networks[END_REF], conveyor belts [START_REF] Göttlich | Modeling, simulation and validation of material flow on conveyor belts[END_REF], weakly coupled oscillators [START_REF] Amadori | On the global well-posedness of {BV} weak solutions to the kuramotosakaguchi equation[END_REF], laser cutting processes [START_REF] Colombo | Nonlocal systems of balance laws in several space dimensions with applications to laser technology[END_REF], biological applications like structured populations dynamics [START_REF] Perthame | Transport equations in biology[END_REF], crowd dynamics [START_REF] Carrillo | An improved version of the Hughes model for pedestrian flow[END_REF][START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF] or more general problems like gradient constrained equations [START_REF] Amorim | On a nonlocal hyperbolic conservation law arising from a gradient constraint problem[END_REF]. While several analytical results on non-local conservation laws can be found in the recent literature (we refer to [START_REF] Amorim | On the numerical integration of scalar nonlocal conservation laws[END_REF] for scalar equations in one space dimension, [START_REF] Ambrosio | Hamiltonian ODEs in the Wasserstein space of probability measures[END_REF][START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Piccoli | Transport equation with nonlocal velocity in Wasserstein spaces: convergence of numerical schemes[END_REF] for scalar equations in several space dimensions and [START_REF] Aggarwal | Nonlocal systems of conservation laws in several space dimensions[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF][START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF] for multi-dimensional systems of conservation laws), few specific numerical methods have been developed up to now. Finite volume numerical methods have been studied in [START_REF] Aggarwal | Nonlocal systems of conservation laws in several space dimensions[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF][START_REF] Kurganov | Non-oscillatory central schemes for a traffic flow model with Arrehenius look-ahead dynamics[END_REF][START_REF] Piccoli | Time-evolving measures and macroscopic modeling of pedestrian flow[END_REF]. At this regard, it is important to notice that the lack of Riemann solvers for non-local equations limits strongly the choice of the scheme. At the best of our knowledge, two main approaches have been proposed in the literature to treat nonlocal problems: first and second order central schemes like Lax-Friedrichs or Nassyau-Tadmor [START_REF] Aggarwal | Nonlocal systems of conservation laws in several space dimensions[END_REF][START_REF] Amorim | On the numerical integration of scalar nonlocal conservation laws[END_REF][START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF][START_REF] Kurganov | Non-oscillatory central schemes for a traffic flow model with Arrehenius look-ahead dynamics[END_REF] and Discontinuous Galerkin (DG) methods [START_REF] Göttlich | Discontinuous Galerkin Method for Material Flow Problems[END_REF]. In particular, the comparative study presented in [START_REF] Göttlich | Discontinuous Galerkin Method for Material Flow Problems[END_REF] on a specific model for material flow in two space dimensions, involving density gradient convolutions, encourages the use of DG schemes for their versatility and lower computational cost, but further investigations are needed in this direction. Besides that, the computational cost induced by the presence of non-local terms, requiring the computation of quadrature formulas at each time step, motivate the development of high order algorithms.

The aim of the present article is to conduct a comparison study on high order schemes for a class of non-local scalar equations in one space dimension, focusing on equations of type (1.1). In Section 2 we review DG and FV-WENO schemes for classical conservation laws. These schemes will be extended to the non-local case in Sections 3 and 4. Section 5 is devoted to numerical tests.

A review of Discontinuous Galerkin and Finite Volume WENO schemes for local conservation laws

The aim of this section is to introduce some notations and to briefly review the DG and FV-WENO numerical schemes to solve the classical local nonlinear conservation law

ρ t + g(ρ) x = 0, x ∈ R, t > 0, ρ(x, 0) = ρ 0 (x), x ∈ R. (2.1)
We first consider {I j } j∈Z a partition of R. The points x j will represent the centers of the cells

I j = [x j-1 2 , x j+ 1 2
], and the cell sizes will be denoted by

∆x j = x j+ 1 2 -x j-1 2
. The largest cell size is h = sup j ∆x j . Note that, in practice, we will consider a constant space step so that we will have h = ∆x.

2.1.

The Discontinuous Galerkin approach. In this approach, we look for approximate solutions in the function space of discontinuous polynomials

V h := V k h = {v : v| I j ∈ P k (I j )}
, where P k (I j ) denotes the space of polynomials of degree at most k in the element I j . The approximate solutions are sought under the form

ρ h (x, t) = k l=0 c (l) j (t)v (l) j (x), v (l) j (x) = v (l) (ζ j (x)),
where c (l) j are the degrees of freedom in the element I j . The subset {v (l) j } l=0,...,k constitutes a basis of P k (I j ) and in this work we will take Legendre polynomials as a local orthogonal basis of P k (I j ), namely

v (0) (ζ j ) = 1, v (1) (ζ j ) = ζ j , v (2) (ζ j ) = 1 2 3ζ 2 j -1 , . . . , ζ j := ζ j (x) =
x -x j ∆x/2 , see for instance [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF][START_REF] Qiu | Runge-Kutta discontinuous Galerkin method using WENO limiters[END_REF]. Multiplying (2.1) by v h ∈ V h and integrating over I j gives

I j ρ t v h dx - I j g(ρ)v h,x dx + g(ρ(•, t))v h (•) x j+ 1 2 x j-1 2 = 0, (2.2)
and the semi-discrete DG formulation thus consists in looking for ρ h ∈ V h , such that for all v h ∈ V h and all j,

I j ρ h t v h dx - I j g(ρ h )v h,x dx + ĝj+ 1 2 v - h (x j+ 1 2 ) -ĝj-1 2 v + h (x j-1 2 ) = 0, (2.3)
where

ĝj+ 1 2 = ĝ(ρ h,- j+ 1 2 , ρ h,+ j+ 1 2
) is a consistent, monotone and Lipschitz continuous numerical flux function. In particular, we will choose to use the Lax-Friedrichs flux

ĝ(a, b) := g(a) + g(b) 2 + α a -b 2 , α = max u |g (u)|.
Let us now observe that if v h is the l-th Legendre polynomial, we have

v + h (x j-1 2 ) = v (l) (ζ j (x j-1 2
)) = (-1) l , and

v + h (x j+ 1 2 ) = v (l) (ζ j (x j+ 1 2
)) = 1, ∀j , l = 0, 1, . . . , k. Therefore, replacing ρ h (x, t) by

ρ h j (x, t) and v h (x) by v (l) (ζ j (x)) in (2.
3), the degrees of freedom c (l)

j (t) + 1 a l - I j g(ρ h j ) d dx v (l) (ζ j (x))dx + ĝ(ρ h,- j+ 1 2 , ρ h,+ j+ 1 2 ) -(-1) l ĝ(ρ h,- j-1 2 , ρ h,+ j-1 2 ) = 0, (2.4)
with

a l = I j (v (l) (ζ j (x))) 2 dx = ∆x 2l + 1 , l = 0, 1, . . . , k.
On the other hand, the initial condition can be obtained using the L 2 -projection of ρ 0 (x), namely c (l)

j (0) = 1 a l I j ρ 0 (x)v (l) (ζ j (x))dx = 2l + 1 2 1 -1 ρ 0 ∆x 2 y + x j v (l) (y)dy, l = 0, . . . , k.
The integral terms in (2.4) can be computed exactly or using a high-order quadrature technique after a suitable change of variable, namely

I j g(ρ h j ) d dx v (l) (ζ j (x))dx = 1 -1 g ρ h j ∆x 2 y + x j , t (v (l) ) (y)dy.
In this work, we will consider a Gauss-Legendre quadrature with N G = 5 nodes for integrals on [-1, 1]

1 -1 g(y)dy = N G e=1 w e g(y e ),
where y e are the Gauss-Legendre quadrature points such that the quadrature formula is exact for polynomials of degree until 2N G -1 = 9 [START_REF] Abramowitz | Handbook of mathematical functions[END_REF].

The semi-discrete scheme (2.4) can be written under the usual form

d dt C(t) = L(C(t)),
where L is the spatial discretization operator defined by (2.4). In this work, we will consider a time-discretisation based on the following total variation diminishing (TVD) third-order Runge-Kutta method [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF],

C (1) = C n + ∆tL(C n ), C (2) = 3 4 C n + 1 4 C (1) + 1 4 ∆tL(C (1) ), (2.5
) 2) ).

C n+1 = 1 3 C n + 2 3 C (2) + 2 3 ∆tL(C (
Other TVD or strong stability preserving SSP time discretization can be also used [START_REF] Gottlieb | High order strong stability preserving time discretizations[END_REF]. The CFL condition is ∆t ∆x max

ρ |g (ρ)| ≤ C CFL = 1 2k + 1 ,
where k is the degree of the polynomial, see [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF]. The scheme (2.4) and (2.5) will be denoted RKDG.

2.2. Generalized slope limiter. It is well-known that RKDG schemes like the one proposed above may oscillate when sharp discontinuities are present in the solution. In order to control these instabilities, a common strategy is to use a limiting procedure. We will consider the so-called generalized slope limiters proposed in [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF]. With this in mind and

ρ h j (x) = k l=0 c (l) j v (l) (ζ j (x)) ∈ P k (I j ), we first set u - j+ 1 2 := ρj + mm(ρ h j (x j+ 1 2 ) -ρj , ∆ + ρj , ∆ -ρj ) and u + j-1 2 := ρj -mm(ρ j -ρ h j (x j-1 2 ), ∆ + ρj , ∆ -ρj ),
where ρj is the average of ρ h on I j , ∆ + ρj = ρj+1 -ρj , ∆ -ρj = ρj -ρj-1 , and where mm is given by the minmod function limiter

mm(a 1 , a 2 , a 3 ) = s • min j |a j | if s = sign(a 1 ) = sign(a 2 ) = sign(a 3 ) 0 otherwise,
or by the TVB modified minmod function

mm(a 1 , a 2 , a 2 ) = a 1 if |a 1 | ≤ M b h 2 , mm(a 1 , a 2 , a 3 ) otherwise, (2.6)
where M b > 0 is a constant. According to [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF][START_REF] Qiu | Runge-Kutta discontinuous Galerkin method using WENO limiters[END_REF], this constant is proportional to the secondorder derivative of the initial condition at local extrema. Note that if M b is chosen too small, the scheme is very diffusive, while if M b is too large, oscillations may appear.

The values u -

j+ 1 2 and u + j-1 2
allow to compare the interfacial values of ρ h j (x) with respect to its local cell averages. Then, the generalized slope limiter technique consists in replacing ρ h j on each cell I j with ΛΠ h defined by

ΛΠ h (ρ h j ) =      ρ h j if u + j-1 2 = ρ h j (x j-1 2
) and u -

j+ 1 2 = ρ h j (x j+ 1 2 ), ρj + (x -x j ) ∆x/2 mm(c (1)
j , ∆ + ρj , ∆ -ρj ) otherwise.
Of course, this generalized slope limiter procedure has to be performed after each inner step of the Runge-Kutta scheme (2.5).

2.3.

The Finite Volume WENO approach. In this section, we solve the nonlinear conservation law (2.1) by using a high-order finite volume WENO scheme [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. Let us denote by ρ(x j , t) the cell average of the exact solution ρ(•, t) in the cell I j :

ρ(x j , t) := 1 ∆x I j ρ(x, t)dx.
The unknowns are here the set of all {ρ j (t)} j∈Z which represent approximations of the exact cell averages ρ(x j , t). Integrating (2.1) over I j we obtain

d dt ρ(x j , t) = - 1 ∆x g(ρ(x j+ 1 2 , t)) -g(ρ(x j-1 2 , t)) , ∀j ∈ Z.
This equation is approximated by the semi-discrete conservative scheme

d dt ρj (t) = - 1 ∆x ĝj+ 1 2 -ĝj-1 2 , ∀j ∈ Z, (2.7)
where the numerical flux ĝj+ . In order to obtain a (2k -1)th-order WENO approximation, we first compute k reconstructed values ρ(r)

j+ 1 2 = k-1 i=0 c i r ρi-r+j , r = 0, . . . , k -1,
that correspond to k possible stencils S r (j) = {x j-r , . . . , x j-r+k-1 } for r = 0, . . . , k -1. The coefficients c i r are chosen in such a way that each of the k reconstructed values is kth-order accurate [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. Then, the (2k -1)th-order WENO reconstruction is a convex combination of all these k reconstructed values and defined by

ρ l j+ 1 2 = k-1 r=0 w r ρ(r) j+ 1 2 ,
where the positive nonlinear weights w r satisfy k-1 r=0 w r = 1 and are defined by

w r = α r k-1 s=0 α s , α r = d r ( + β r ) 2 .
Here d r are the linear weights which yield the (2k -1)th-order accuracy, β r are called the "smoothness indicators" of the stencil S r (j), which measure the smoothness of the function ρ in the stencil, and is a small parameter used to avoid dividing by zero (typically = 10 -6 ). The exact form of the smoothness indicators and other details about WENO reconstructions can be found in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. The reconstruction of ρ r j-1 2 is obtained in a mirror symmetric fashion with respect to x j-1

2

. The semi-discrete scheme (2.7) is then integrated in time using the (TVD) third-order Runge-Kutta scheme (2.5), under the CFL condition

∆t ∆x max ρ |g (ρ)| ≤ C CFL < 1.
.

Construction of DG schemes for non-local problems

We now focus on the non-local equation (1.1), for which we set R(x, t) := (ρ * ω η )(x, t). Since (l)

ρ h (x, t)| I j = k l=0 c (l) j (t)v (l) (ζ j (x)) ∈ P k (I j ),
j (t) + 1 a l - I j f (ρ h )V (R) d dx v (l) (ζ j (x))dx + fj+ 1 2 -(-1) l fj-1 2 = 0, (3.1)
where fj+ 1 2 is a consistent approximation of f (ρ)V (R) at interface x j+1/2 . Here, we consider again the Lax-Friedrichs numerical flux defined by . Since R is defined by a convolution, we naturally set R h,-

fj+ 1 2 = 1 2 f (ρ h,- j+ 1 2 )V (R h,- j+ 1 2 ) + f (ρ h,+ j+ 1 2 )V (R h,+ j+ 1 2 ) + α(ρ h,- j+ 1 2 -ρ h,+ j+ 1 2 ) , (3.2)
with α := max ρ |∂ ρ (f (ρ)V (ρ))| and where R h,-
j+ 1 2 = R h,+ j+ 1 2 = R(x j+ 1 2 , t) =: R j+ 1 2
, so that (3.2) can be written as

fj+ 1 2 := f (ρ h,- j+ 1 2 , ρ h,+ j+ 1 2 , R j+ 1 2 ) = 1 2 (f (ρ h,- j+ 1 2 ) + f (ρ h,+ j+ 1 2 ))V (R j+ 1 2 ) + α(ρ h,- j+ 1 2 -ρ h,+ j+ 1 2 ) . (3.3)
Next, we propose to approximate the integral term in (3.1) using the following high-order Gauss-Legendre quadrature technique,

I j f (ρ h )V (R) d dx v (l) (ζ j (x))dx = 1 -1 f (ρ h ∆x 2 y + x j , t )V R ∆x 2 y + x j , t (v (l) ) (y)dy = N G e=1 w e f (ρ h (x e , t))V (R (x e , t)) (v (l) ) (y e ), (3.4)
where we have set xe = ∆x 2 y e +x j ∈ I j , y e being the Gauss-Legendre quadrature points ensuring that the quadrature formula is exact for polynomials of order less or equal to 2N G -1.

It is important to notice that the DG formulation for the non-local conservation law (3.1) requires the computation of the extra integral terms R j+ 1 2 in (3.3) and R (x e , t) in (3.4) for each quadrature point, which increases the computational cost of the strategy. For ρ h (x, t)| I j ∈ P k (I j ), we can compute these terms as follows for both the LWR and sedimentation non-local models considered in this paper.

Non-local LWR model. For the non-local LWR model, we impose the condition N ∆x = η for some N ∈ N, so that we have

R j+ 1 2 = x j+ 1 2 +η x j+ 1 2 ω η (y -x j+ 1 2 )ρ h (y, t)dy = N i=1 I j+i ω η (y -x j+ 1 2 )ρ h j+i (y, t)dy = N i=1 k l=0 c (l)
j+i (t)

I j+i ω η (y -x j+ 1 2 )v (l) (ζ j+i (y))dy = N i=1 k l=0 c (l) j+i (t) ∆x 2 1 -1 ω η ∆x 2 y + (i - 1 2
)∆x v (l) (y)dy

Γ i,l = N i=1 k l=0 c (l) j+i Γ i,l ,
while for each quadrature point xe we have

R (x e , t) = xe+η xe ω η (y -xe )ρ h (y, t)dy = x j+ 1 2 xe ω η (y -xe )ρ h j (y, t)dy Γa + N -1 i=1 I j+i ω η (y -xe )ρ h j+i (y, t)dy Γ i b + xe+η x j+N -1 2 ω η (y -xe )ρ h j+N (y, t)dy Γc .
The three integrals Γ a , Γ i b and Γ c can be computed with the same change of variable as before, namely

Γ a = ∆x 2 1 ye ω η ∆x 2 (y -y e ) k l=0 c (l) j (t)v (l) (y)dy = k l=0 c (l) j (t) ∆x 2 1 ye ω η ∆x 2 (y -y e ) v (l) (y)dy Γ (e) 0,l = k l=0 c (l) j Γ (e) 0,l , Γ i b = ∆x 2 1 -1 ω η ∆x 2 (y -y e ) + i∆x k l=0 c (l) j+i (t)v (l) (y)dy = k l=0 c (l) j+i ∆x 2 1 -1 ω η ∆x 2 (y -y e ) + i∆x v (l) (y)dy Γ (e) i,l = k l=0 c (l) j+i Γ (e)
i,l ,

Γ c = ∆x 2 ye -1 ω η ∆x 2 (y -y e ) + N ∆x k l=0 c (l) j+N v (l) (y)dy = k l=0 c (l) j+N (t) ∆x 2 ye -1 ω η ∆x 2 (y -y e ) + N ∆x v (l) (y)dy Γ (e) N,l = k l=0 c (l) j+N Γ (e)
N,l .

Finally we can compute R (x e , t) as

R (x e , t) = N i=0 k l=0 c (l) j+i (t)Γ (e)
i,l , in order to evaluate (3.4).

Non-local sedimentation model. Considering now the non-local sedimentation model, we impose N ∆x = 2η for some N ∈ N, so that we have Remark. In order to compute integral terms in (3.1) as accurately as possible, the integrals R j+ 1 2 and R (x e , t) above, and in particular, the coefficients Γ i,l , must be calculated exactly or using a suitable quadrature formula accurate to at least O(∆x l+p ) where p is the degree of the convolution term ω η . It is important to notice that the coefficients can be precomputed and stored in order to save computational time. Finally the semi-discrete scheme (3.1) can be discretized in time using the (TVD) third-order Runge-Kutta scheme (2.5), under the CFL condition ∆t ∆x max

R j+ 1 2 = x j+ 1 2 +2η x j+ 1 2 -2η ω η (y -x j+ 1 2 )ρ h (y, t)dy = N i=-N +1 k l=0 c (l) j+i (t)Γ i,
ρ |∂ ρ (f (ρ)V (ρ))| ≤ C CFL = 1 2k + 1 ,
where k is the degree of the polynomial.

Construction of FV schemes for non-local conservation laws

Let us now extend the FV-WENO strategy of Section 2.3 to the non-local case. We first integrate (1.1) over I j to obtain

d dt ρ(x j , t) = - 1 ∆x f (ρ(x j+ 1 2 , t))V (R(x j+ 1 2 , t)) -f (ρ(x j-1 2 , t))V (R(x j-1 2 , t)) , ∀j ∈ Z,
so that the semi-discrete discretization can be written as

d dt ρj (t) = - 1 ∆x fj+ 1 2 -fj-1 2 , ∀j ∈ Z, (4.1)
where the use of the Lax-Friedrichs numerical flux gives

fj+ 1 2 = f (ρ l j+ 1 2 , ρ r j+ 1 2 , R j+ 1 2 ) = 1 2 (f (ρ l j+ 1 2 ) + f (ρ r j+ 1 2 ))V (R j+ 1 2 ) + α(ρ l j+ 1 2 -ρ r j+ 1 2
) .

Recall that ρ l At this stage, it is crucial to notice that in the present FV framework, the approximate solution is piecewise constant on each cell I j , so that a naive evaluation of the convolution terms may lead to a loss of high-order accuracy. Let us illustrate this. Considering for instance the non-local LWR model and using that ρ(x, t) is piecewise constant on each cell naturally leads to

R j+ 1 2 = x j+ 1 2 +η x j+ 1 2 ω η (y -x j+ 1 2 )ρ(y, t)dy = N i=1 I j+i ω η (y -x j+ 1 2 )ρ(y, t)dy = ∆x N i=1 ρj+i I j+i ω η (y -x j+ 1 2
)dy, which does not account for the high-order WENO reconstruction. In order to overcome this difficulty, we propose to approximate the value of ρ(x, t) using quadratic polynomials in each cell. This strategy is detailed for each model in the following two subsections. for the non-local LWR model. In order to compute the integral

R j+ 1 2 = x+η x j+ 1 2 ω η (y -x j+ 1 2 )ρ(y, t)dy,
we propose to consider a reconstruction of ρ(x, t) on I j by taking advantage of the high-order WENO reconstructions ρ r j-1 2 and ρ l j+ 1 2 at the boundaries of I j , as well as the approximation of the cell average ρn j . More precisely, we propose to define a polynomial P j (x) of degree 2 on I j by

P j (x j-1 2 ) = ρ r j-1 2 , P j (x j+ 1 2 ) = ρ l j+ 1 2 , 1 ∆x I j P j (x)dx = ρn j ,
which is very easy to handle. In particular, we have

P j (x) := a j,0 + a j,1 x -x j ∆x/2 + a j,2 3 x -x j ∆x/2 2 -1 /2, x ∈ I j , with a j,0 = ρn j , a j,1 = 1 2 ρ l j+ 1 2 -ρ r j-1 2 , a j,2 = 1 2 ρ l j+ 1 2 + ρ r j-1 2 -ρn j .
Observe that we have used the same polynomials as in the DG formulation, i.e., P j (x) = 2 l=0 a j,l v (l) (ζ j (y)). With this, R j+ 1 2 can be computed as

R j+ 1 2 = N i=1 I j+i ω η (y -x j+ 1 2 )P j+i (y)dy = N i=1 I j+i ω η (y -x j+ 1 2 ) 2 l=0 a j+i,l v (l) (ζ j+i (y))dy = N i=1 2 l=0 a j+i,l I j+i ω η (y -x j+ 1 2 )v (l) (ζ j+i (y))dy = N i=1 2 l=0 a j+i,l ∆x 2 1 -1 ω η ∆x 2 y + (i - 1 2 )∆x v (l) (y)dy = N i=1 2 l=0 a j+i,l Γ i,l ,
where the coefficients Γ l,i are computed exactly or using a high-order quadrature approximation.

R j+ 1 2 = 2η -2η ω η (y)ρ(x j+ 1 2 + y, t)dy = x j+ 1 2 +2η x j+ 1 2 -2η ω η (y -x j+ 1 2
)ρ(y, t)dy.

Considering again the assumption N ∆x = 2η, we get

R j+ 1 2 = N i=-N +1 I j+i ω η (y -x j+ 1 2 )P j+i (y)dy = N i=-N +1 2 l=0 a j+i,l I j+i ω η (y -x j+ 1 2 )v (l) (ζ j+i (y))dy = N i=-N +1 2 l=0 a j+i,l ∆x 2 1 -1 ω η ∆x 2 y + (i - 1 2 )∆x v (l) (y)dy = N i=-N +1 2 l=0 a j+i,l Γ i,l .
To conclude this section, let us remark that the coefficients Γ i,l are computed for l = 0, . . . , k in the DG formulation, where k is the degree of the polynomials in P k (I j ), while in the FV formulation, the coefficients Γ i,l are computed only for l = 0, . . . , 2 due to the quadratic reconstruction of the unknown in each cell. The semi-discrete scheme (4.1) is then integrated in time using the (TVD) third-order Runge-Kutta scheme (2.5), under the CFL condition

∆t ∆x max ρ |∂ ρ (f (ρ)V (ρ))| ≤ C CFL < 1.
.

Numerical experiments

In this section, we propose several test cases in order to illustrate the behaviour of the RKDG and FV-WENO high-order schemes proposed in the previous Sections 3 and 4 for the numerical approximation of the solutions of the non-local traffic flow and sedimentation models on a bounded interval I = [0, L] with boundary conditions.

We consider periodic boundary conditions for the traffic flow model, i.e. ρ(0, t) = ρ(L, t) for t ≥ 0 and zero-flux boundary conditions for sedimentation model, in order to simulate a batch sedimentation process. More precisely, we assume that ρ(x, t) = 0 for x ≤ 0 and ρ(x, t) = 1 for x ≥ 1. Given a uniform partition of [0, L] {I j } M j=1 with ∆x = L/M , in order to compute the numerical fluxes fj+1/2 for j = 0, . . . , M + 1, we define ρ n j in the ghost cells as follow: for the traffic model,

ρ n -1 := ρ n M -1 , ρ n 0 := ρ n M , ρ n M +i := ρ n i , for i = 1, . . . , N ;
and for the sedimentation model ρ n i := 0, for i = -N, . . . , 0, and ρ n M +i := 1, for i = 1, . . . , N.

A key element of this section will be the computation of the Experimental Order of Accuracy (EOA) of the proposed strategies, which is expected to coincide with the theoretical order of accuracy given by the high-order reconstructions involved in the corresponding numerical schemes. Let us begin with a description of the practical computation of the EOA. Regarding the RKDG schemes, if ρ ∆x (x, T ) and ρ ∆x/2 (x, T ) ∈ V k h are the solutions computed with M and 2M mesh cells respectively, the L 1 -error is computed by

e(∆x) = ρ ∆x (x, T ) -ρ ∆x/2 (x, T ) L 1 = M j=1 I 2j-1 |ρ ∆x j (x, T ) -ρ ∆x/2 2j-1 (x, T )|dx + I 2j |ρ ∆x j (x, T ) -ρ ∆x/2 2j (x, T )|dx,
where the integrals are computed with a high-order quadrature formula. As far as the FV schemes are concerned, the L 1 -error is computed as

e(∆x) = ρ ∆x (x, T ) -ρ ∆x/2 (x, T ) L 1 = ∆x M j=1 |ρ ∆x j (x j , T ) -ρj (x j )|dx,
where ρj (x) is a third-degree polynomial reconstruction of ρ ∆x/2 j (x, T ) at point x j , i.e.,

ρj (x j ) = 9 16 ρ ∆x/2 2j + ρ ∆x/2 2j-1 - 1 16 ρ ∆x/2 2j+1 + ρ ∆x/2 2j-2 .
In both cases, the EOA is naturally defined by γ(∆x) = log 2 (e(∆x)/e(∆x/2)) .

In the following numerical tests, the CFL number is taken as C CFL = 0.2, 0.1, 0.05 for RKDG1, RKDG2 and RKDG3 schemes respectively, and C CFL = 0.6 for FV-WENO3 and FV-WENO5 schemes. For RKDG3 and FV-WENO5 cases ∆t is further reduced in the accuracy tests. 2) at time T = 0.1 with η = 0.1 and w η (x) = 3(η 2 -x 2 )/(2η 3 ). We set ∆x = 1/800 and compare the numerical solutions obtained with the FV-WENO3, FV-WENO5, and for RKDG1, RKDG2 and RKDG3 we use a generalized slope limiter (2.6) with M b = 35. The results displayed in Fig. 1 are compared with respect to the reference solution which was obtained with FV-WENO5 scheme and ∆x = 1/3200. In Fig1 (b) we observe that RKDG schemes are much more accurate than FV-WENO3 and even FV-WENO5.

5.2. Test 1b: non-local LWR model. Now, we consider an initial condition with a small perturbation ρ 0 (x) = 0.35 -(x -0.5) exp (-2000(x -0.5) 2 ) and an increasing kernel function ω η = 2x/η 2 , with η = 0.05 and periodic boundary conditions. According to [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF], the non-local LWR model is not stable with increasing kernels, in the sense that oscillations develop in short time. Fig. 2 displays the numerical solution with different RKDG and FV-WENO schemes with ∆x = 1/400 at time T = 0.3. The profile provided by the first order scheme proposed in [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF] is also included. The reference solution is computed with a FV-WENO5 scheme with ∆x = 1/3200. We observe that the numerical solutions obtained with the high-order schemes provide better approximations of the oscillatory shape of the solution than the first-order scheme, and that RKDG schemes give better approximations than the FV-WENO schemes. x ≤ 0, 0.5 0 < x < 1, 1

x ≥ 1, with zero-flux boundary conditions in the interval [0, 1], and compute the solution at time T = 1, with parameters α = 1, n = 3 and a = 0.025. We set ∆x = 1/400 and compute the solutions with different RKDG and FV-WENO schemes, including the first-order Lax-Friedrichs scheme used in [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF]. The results displayed in Fig. 3 are compared to a reference solution computed with FV-WENO5 and ∆x = 1/3200. Compared to the reference solution, we observe that RKDG1 is more accurate than FV-WENO3, and FV-WENO5 more accurate than RKDG3 and RKDG2 (Fig. 3). In particular, we observe that the numerical solutions obtained with the high-order schemes provide better approximations of the oscillatory shape of the solution than the first order scheme. These oscillations can possibly be explained as layering phenomenon in sedimentation [START_REF] Siano | Layered sedimentation in suspensions of monodisperse spherical colloidal particles[END_REF], which denotes a traveling staircases pattern, looking as several distinct bands of different concentrations. We observe in Fig. 4, where the evolution of ρ h (•, t) is displayed for t ∈ [0, 5], that this layering phenomenon is observed with high order schemes, in this case RKDG1 and FV-WENO5, instead of Lax-Friedrichs scheme. ). We compare solutions computed with FV Lax-Friedrichs, FV-WENO and RKDG schemes, using ∆x = 1/400. Reference solution is computed with FV-WENO5 with ∆x = 1/3200. 5.4. Test 3: Experimental Order of Accuracy for local conservation laws. In this subsection we compute the EOA for local conservation laws. More precisely, we first consider the advection equation ρ t +ρ x = 0, and the initial datum ρ 0 (x) = 0.5+0.4 sin(πx) for x ∈ [-1, 1] with periodic boundary conditions. We compare the numerical approximations obtained with ∆x = 1/M and M = 20, 40, 80, 160, 320, 640 at time T = 1. RKDG schemes are used without limiting function and L 1 -errors are collected in Table 1. We observe that the correct EOA is obtained for the different schemes, moreover, we observe that the L 1 -errors obtained with FV-WENO3 and RKDG1 are comparable, and also the L 1 -errors for FV-WENO5 and RKDG3. We also consider the nonlinear local LWR model ρ t + (ρ(1 -ρ)) x = 0, with initial datum ρ 0 (x) = 0.5 + 0.4 sin(πx) for x ∈ [-1, 1] and periodic boundary conditions. A shock wave appears at time T = 0.3. We compare the numerical approximations obtained with ∆x = 1/M and M = 20, 40, 80, 160, 320, 640 at time T = 0.15. RKDG schemes are used without limiters and L 1 -errors are given in Table 2. We again observe that the correct EOA is obtained for the different schemes, and the L 1 -errors obtained with FV-WENO3 and RKDG1 are comparable, and also the L 1 -errors for FV-WENO5 and RKDG3. We observe that the same kind of layering phenomenon reported in [START_REF] Siano | Layered sedimentation in suspensions of monodisperse spherical colloidal particles[END_REF] seems to be observed with high order schemes RKDG1 and FV-WENO5 instead of LxF scheme.

For RKDG schemes, we obtain the correct EOA. For the FV-WENO schemes, the EOA is also correct thanks to the in-cell quadratic reconstructions used to compute the non-local terms.

For the non-local sedimentation model, we compute the solution of (1.1)-(1.3)-(1.4) with initial data ρ 0 (x) = 0.8 sin(πx) 10 for x ∈ [0, 1] with η = 0.025, α = 1, V (ρ) = (1 -ρ) 3 and ∆x = 1/M with M = 100, 200, 400, 800, 1600, 3200 at final time T = 0.04. The results are given in Table 3.

k = 1 k = 2 k = 3 FV-WENO3 FV-WENO5 1/∆x L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) 20
(x) = 0.5 + 0.4 sin(πx) on x ∈ [-1, 1]. k = 1 k = 2 k = 3 WENO3 WENO5 1/∆x L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x) 20

Conclusion

In this paper we developed high order numerical approximations of the solutions of nonlocal conservation laws in one space dimension motivated by application to traffic flow and sedimentation models. We propose to design Discontinuous Galerkin (DG) schemes which can be applied in a natural way and Finite Volume WENO (FV-WENO) schemes where we have used quadratic polynomial reconstruction in each cell to evaluate the convolution term in order to obtain the high-order accuracy. The numerical solutions obtained with high-order schemes provide better approximations of the oscillatory shape of the solutions than the first order schemes. We also remark that DG schemes are more accurate but expensive, while FV-WENO are less accurate but also less expensive. The present work thus establishes the preliminary basis for a deeper study and optimum programming of the methods, so that efficiently plots could then be considered. w η (x) = 1/η w η (x) = 2(η -x)/η 2 w η (x) = 3(η 2 -x 2 )/(2η 3 ) 1/∆x L 1 -error γ(∆x) L 1 -error γ(∆x) L 10 . Experimental order of accuracy at T = 0.04 with f (ρ) = ρ(1 -ρ), V (ρ) = (1 -ρ) 3 , ω η (x) := η -1 K(η -1 x) with η = 0.025.

1 2:= ĝ(ρ l j+ 1 2 , ρ r j+ 1 2 ) 1 2 and ρ r j+ 1 2 2 ,

 122112 is the Lax-Friedrichs flux and ρ l j+ are some left and right high-order WENO reconstructions of ρ(x j+ 1 t) obtained from the cell averages {ρ j (t)} j∈Z . Let us focus on the definition of ρ l j+ 1 2

  is a weak solution of the non-local problem (1.1), the coefficients c

  (l) j (t) can be calculated by solving the following differential equation, d dt c

j+ 1 2 and R h,+ j+ 1 2

 21 are the left and right approximations of R(x, t) at the interface x j+ 1 2

  l , and for each quadrature point xe , R (x e , t) = xe+2η xe-2η ω η (y -xe )ρ h (y, t)dy = y e ) -N ∆x v (l) (y)dy, Γ y e ) + i∆x v (l) (y)dy, Γ y e ) + N ∆x v (l) (y)dy.

j+ 1 2 and ρ r j+ 1 2 2 .

 212 are the left and right WENO high-order reconstructions at point x j+1

4. 1 .

 1 Computation of R j+ 1 2

4. 2 . 1 2

 21 Computation of R j+ for non-local sedimentation model. As far as the non-local sedimentation model (1.1)-(1.3) is concerned, we have

5. 1 .

 1 Test 1a: non-local LWR model. We consider the Riemann problem ρ(x, 0) = 0.95 x ∈ [-0.5, 0.4], 0.05 otherwise, with absorbing boundary conditions and compute the numerical solution of (1.1)-(1.

Figure 1 .

 1 Figure 1. Test 1a. (a) Solution of Riemann problem with w η (x) = 3(η 2x 2 )/(2η 3 ) at T = 0.1. We compare solutions computed with FV-WENO and RKDG schemes using ∆x = 1/800. (b) Zoomed view of the left discontinuity in Fig (a). Reference solution is computed with FV-WENO5 with ∆x = 1/3200.

5. 3 .

 3 Test 2: non-local sedimentation model. We now solve (1.1)-(1.3)-(1.4) with the piecewise constant initial datum ρ(x,

Figure 2 .

 2 Figure 2. Test1b. Numerical solution non-local LWR model with ω η = 2x/η 2 , η = 0.05 and ρ 0 (x) = 0.35 -(x -0.5) exp (-2000(x -0.5) 2). We compare solutions computed with FV Lax-Friedrichs, FV-WENO and RKDG schemes, using ∆x = 1/400. Reference solution is computed with FV-WENO5 with ∆x = 1/3200.

5. 5 .

 5 Test 4: Experimental Order of Accuracy for the non-local problems. We now consider the non-local LWR and sedimentation models. Considering the non-local LWR model, we compute the solution of (1.1)-(1.2) with initial data ρ 0 (x) = 0.5 + 0.4 sin(πx) for x ∈ [-1, 1], periodic boundary conditions, with η = 0.1 and ∆x = 1/M and M = 80, 160, 320, 640, 1280 at final time T = 0.15. The results are given for different kernel functions w η (x) in Table3.

Figure 3 .

 3 Figure 3. Test2. Riemann problem for non-local sedimentation model. Comparing numerical solution computed with FV Lax-Friedrichs, FV-WENO schemes and RKDG schemes using limiter function with M b = 50 using ∆x = 1/400. Reference solution is computed with FV-WENO5 with ∆x = 1/3200.

Figure 4 .

 4 Figure 4. Test 2. Nonlocal sedimentation model. History solution for t ∈ [0, 5]computed with LxF (left) RKDG1 (center) and WENO 5 (right) with ∆x = 1/400. We observe that the same kind of layering phenomenon reported in[START_REF] Siano | Layered sedimentation in suspensions of monodisperse spherical colloidal particles[END_REF] seems to be observed with high order schemes RKDG1 and FV-WENO5 instead of LxF scheme.

Table 1 .

 1 Test 3: Advection equation. Experimental order of accuracy for FV-WENO schemes and RKDG schemes without Limiter at T = 1 with initial condition ρ 0

		1.80E-03	-	9.70E-05	-	1.00E-06	-	3.72E-02	-	4.63E-04	-
	40	4.30E-04 2.06 1.28E-05 2.92 6.04E-08 4.05 1.17E-02 1.66 1.41E-05 5.03
	80	1.06E-04 2.02 1.63E-06 2.96 3.92E-09 3.94 2.85E-03 2.04 4.33E-07 5.02
	160 2.64E-05 2.00 2.06E-07 2.98 2.58E-10 3.92 4.24E-04 2.74 1.35E-08 5.00
	320 6.61E-06 2.00 2.59E-08 2.99 1.51E-11 3.32 3.06E-05 3.79 1.15E-11 5.20

Table 2 .

 2 Test

		2.03E-03	-	1.88E-04	-	5.03E-06	-	6.92E-03	-	2.80E-04	-
	40	4.98E-04 2.02 3.06E-05 2.62 2.81E-07 4.16 2.37E-03 1.77 1.41E-05 4.31
	80	1.23E-04 2.01 4.71E-06 2.69 1.72E-08 4.03 5.95E-04 1.99 6.14E-07 4.52
	160 3.08E-05 2.00 7.11E-07 2.72 1.05E-09 4.02 8.20E-05 2.86 2.62E-08 4.55
	320 7.71E-06 2.00 1.06E-07 2.74 6.62E-11 3.99 5.47E-06 3.90 8.61E-10 4.92

3: Classical LWR equation. Experimental order of accuracy for FV-WENO schemes and RKDG schemes without Limiter at T = 0.15 with initial condition ρ 0 (x) = 0.5 + 0.4 sin(πx) on x ∈ [-1, 1].

Table 3 .

 3 Test 4: Non-local LWR model. Experimental order of accuracy. Initial condition ρ 0 (x) = 0.5 + 0.4 sin(πx), η = 0.1, numerical solution at T = 0.15 for FV-WENO and RKDG schemes without generalized slope limiter. L 1 -err γ(∆x) L 1 -err γ(∆x) L 1 -err γ(∆x)L 1 -err γ(∆x) L 1 -err γ(∆x) 98E-05 2.0 4.64E-07 2.6 3.95E-08 3.5 6.02E-05 2.1 1.01E-07 4.7 400 4.95E-06 2.0 7.75E-08 2.6 3.37E-09 3.6 7.87E-06 2.9 3.22E-09 4.9 800 1.23E-06 2.0 1.23E-08 2.6 2.47E-10 3.8 5.50E-07 3.8 1.11E-10 4.8 1600 3.09E-07 2.0 2.02E-09 2.6 1.22E-11 4.2 4.64E-08 3.7 4.02E-12 4.7

	1 -error	γ(∆x)

Table 4 .

 4 Test 4: Non-local sedimentation problem. Initial condition ρ 0 (x) = 0.8 sin(πx)

Acknowledgements

This research was conducted while LMV was visiting Inria Sophia Antipolis Méditerranée and the Laboratoire de Mathématiques de Versailles, with the support of Fondecyt-Chile project 11140708 and "Poste Rouge" 2016 of CNRS -France.

The authors thank Régis Duvigneau for useful discussions and advices.