Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity

We consider an extension of the traffic flow model proposed by Lighthill, Whitham and Richards, 
in which the mean velocity depends on a weighted mean of the downstream traffic density. 
We prove well-posedness and a regularity result for entropy weak solutions of the corresponding Cauchy problem, 
and use a finite volume central scheme to compute approximate solutions. 
We perform numerical tests to illustrate the theoretical results and to investigate the limit 
as the convolution kernel tends to a Dirac delta function.

1. Introduction. Macroscopic traffic flow models provide nowadays a validated and powerful approach to simulate and manage traffic on road networks (we refer the interested reader to [25] for an overview of modelling approaches and practical applications and to [11] for a detailed review of the mathematical theory of macroscopic models on networks). These models are based on hyperbolic equations derived from fluid dynamics, and describe the spatio-temporal evolution of macroscopic quantities like vehicle density and mean velocity. One of the seminal macroscopic models was introduced in the mid 1950s by Lighthill and Whitham [21] and Richards [23], who proposed to complete the one-dimensional mass conservation equation ∂ t ρ(t, x) + ∂ x f (t, x) = 0 with a closure relation between speed and density, leading to the fundamental diagram f (t, x) = f (ρ(t, x)) = ρ(t, x)v(ρ(t, x)), which can be derived from empirical speed-density or flow-density data. The classical LWR model therefore reads ∂ t ρ(t, x) + ∂ x (ρ(t, x)v(ρ(t, x))) = 0, x ∈ R, t > 0, where ρ(t, x) ≥ 0 is the traffic density and the non-increasing function v : R + → R + is the mean velocity. Thanks to its simplicity and its established analytical and numerical properties, this model is widely used in modern traffic engineering, see e.g. [6,14] for some recent applications. Nevertheless, like every model describing complex systems, it suffers form some limitations that limit its applicability. In particular, like all classical macroscopic models, it allows for speed discontinuities, which translate in infinite acceleration values. This restrains the use of these equations in combination with consumption and atmospheric and noise pollution models, which require an accurate estimation of the acceleration terms. Aiming to overcome this difficulty, a prototype of non-local traffic flow models has been recently introduced in [5], where the speed is assumed to depend on a weighted mean of the downstream traffic conditions. Due to the action of a convolution product, the speed becomes a Lipschitz function of the space and time variables, guaranteeing bounded acceleration.
In this paper, we provide a generalization and some extensions of the results presented in [5]: instead of focussing on the linear decreasing function v(ρ) = 1 − ρ, we consider general non-increasing velocities, we provide a regularity result and study a higher order finite volume central scheme for a better approximation of the solutions. More precisely, we consider the following mass conservation equation for traffic flow with non-local mean velocity: In (1), we consider a speed function v : where ρ max denotes the maximal car density. We denote the downstream convolution product as This term is intended to reproduce the fact that drivers adapt their velocity to the downstream traffic, assigning a greater importance to closer vehicles. Setting We will consider solutions ρ = ρ(t, x) satisfying the following definition (see [4,5,8,17]): for all ϕ ∈ C 1 c (R 2 ; R + ) and κ ∈ R. Conservation laws with non-local in space terms arose recently in several application domains, such as granular flows [2], sedimenatation [4], pedestrian flows [7], conveyor belts [13] and aggregation phenomena [16]. Also, other traffic flow models with non-local terms have been proposed and studied in [15,19,24]. Some general analytical results on non-local conservation laws, proving existence and eventually uniqueness of solutions, can be found in [3] for scalar equations in one space dimension, in [8] for equations in several space dimension and in [1,9,10] for multi-dimensional systems of conservation laws.
Unlike all the above mentioned results, in the case of equation (1) we deal with a very specific (non-increasing, discontinuous) convolution kernel. Despite its lack of regularity at the origin, the monotonicity of the kernel allows to obtain refined L ∞ and BV estimates on Lax-Friedrichs approximate solutions, which as a by-product allow to prove the existence of weak entropy solutions by a standard compactness argument. In particular, we prove that a maximum principle holds, which guarantees the positivity of solutions and a uniform upper bound, given by the maximum of the initial density. We stress that these properties, which are of utmost importance in traffic flow applications, depend heavily on the non-increasing monotonicity of the kernel, as numerical counterexamples show, see [5,Section 5].
The paper is organized as follows: Section 2 is devoted to the analytical results, summarized in Theorem 2.1 and some intermediate relevant statements. Section 3 describes a generalisation of the Nessyahu-Tadmor finite volume central scheme and Section 4 collects the results of the numerical tests. The final Appendix contains some interesting details of the proof of the BV-estimates.
2. Well-posedness and regularity. The main results are reported in the following theorem: admits a unique weak entropy solution in the sense of Definition 1.1. Moreover it holds min Uniqueness of solutions is proved in [5]. The proof of existence is based on the convergence of a sequence of approximate solutions as in [5], see also [3,4,7,8,9] for related results. Here we briefly describe the finite volume scheme used to construct approximate solutions and we state its key proprieties: maximum principle, bounded total variation and discrete entropy inequalities, which allow to prove the convergence to a weak entropy solution.
Let us consider a space step ∆x such that η = N ∆x, for some N ∈ N, and a time step ∆t subject to a CFL condition which will be specified later. For j ∈ Z and n ∈ N, let x j+1/2 = j∆x be the cells interfaces, x j = (j − 1/2)∆x the cells centers and t n = n∆t the time mesh. To construct a finite volume approximate solution we take a piece-wise constant approximation of the initial datum ρ 0 , we denote w k η := w η (k∆x) for k = 0, . . . , N − 1 and Note that the discretization chosen for w η implies so it slightly underestimates traffic mean velocity and may introduce unphysical negative velocities. Other discretizations are possible, see for example [4, eq. (3. 2)], but our choice is simpler to implement and does not sensibly affect the results. We consider the following modified Lax-Friedrichs flux adapted to (1): with λ = ∆t/∆x. Easy calculations, detailed in [12], show that the numerical scheme is not monotone w. r. to the variables ρ j+k , k = 2, . . . , N − 2. Therefore, classical convergence arguments do not apply in this context. To prove the convergence of the numerical approximations, we will rely on the following results, see [5,12] for details of the proofs.

Proposition 1. (Maximum principle and L
Then the finite volume approximation ρ n j , j ∈ Z and n ∈ N, constructed using scheme (7), satisfies the bounds ρ m ≤ ρ n j ≤ ρ M for all j ∈ Z and n ∈ N, under the CFL condition This is an important property for traffic flow applications. In particular, we expect solutions to (1) not to exceed the maximal density bound ρ max . It would be interesting to see if different convolution choices are able to preserve this upper bound property while allowing density to increase with respect to the initial values.
, and let ρ ∆x be given by (7). If α ≥ v max + A w η (0)∆x (and α ≥ 2A w η (0)∆x) and the CFL condition ∆t ≤ ∆x/(α + 2A w η (0)∆x)) holds, then for every T > 0 the following discrete space BV estimate holds In particular, solutions to (1) are not total variation decreasing, unlike the solutions of the corresponding classical (non non-local) conservation law. Moreover, it is worth noticing that the bound provided by (9) depends on the L ∞ -norm of the derivatives of the velocity function v. Also, unlike the case of linear velocity studied in [5], the scheme is not monotonicity preserving for general functions v. We refer the reader to the Appendix for a sketch of the proof.
Following [3, Proposition 2.8] and [5, Section 3.3], we derive a discrete entropy inequality for the approximate solution generated by (7). Let us denote Proposition 3. (Discrete entropy inequalities) Let ρ n j , j ∈ Z, n ∈ N, be given by (7). Then, if α ≥ 1 and the CFL condition (8) holds, for all j ∈ Z, n ∈ N we have The proof of convergence of ρ ∆x to an entropy weak solution in the sense of Definition 1.1 follows closely [5,Section 4] and is detailed in [12,Section 2.5].
Concerning the regularity of the solutions of (1), (2), we recover the same result as in [4,Lemma 5.5].
Proposition 4. (Lipschitz regularity) For any T > 0, the numerical solution generated by the scheme (7) converges to a Lipschitz continuous function ρ, provided ρ 0 is also Lipschitz continuous.
The proof is detailed in [12, Section 2.6].

3.
A finite volume central scheme. In this section, we describe a central scheme for (1) inspired by [18], which is an extension of the second-order Nessyahu-Tadmor scheme [22]. At each time step, the solution is first approximated by a piece-wise linear function, then evolved to the next time step according to the integral form of conservation law.
Denoting by the flux function in (1) can be written as With the notations introduced at the beginning of Section 2, let us assume that at time t = t n we have the solution approximated by its cell averages We construct the piecewise linear interpolant approximation given bỹ Formally, we have a second-order approximation provided the slopes s n j are at least first-order approximations of the derivatives ρ x (t n , x j ). We use a generalized minmod reconstruction (as in [18]) with: where the minmod function is defined by and the parameter θ influences the numerical viscosity of the scheme: larger values of θ correspond to less dissipative but more oscillatory approximations, see [20]. We now evolve the approximated solution (11) to the next time level t = t n+1 by integrating equation (10) over the space-time volumes [x j , x j−1 ] × [t n , t n+1 ], which gives: ρ n+1 The first integral on the right-hand side of (13) can be computed exactly: The flux integral in (13) should be computed using the approximate solution of the initial value problem (10), (11) in the time interval (t n , t n+1 ) with initial data obtained at t = t n . Given CFL condition we compute the flux integrals in (13) by the mid-point quadrature: where we approximate the intermediate time level values of ρ and R with the corresponding Taylor expansions, namely: By (11) we have:ρ n (x j ) =ρ n j and by (10) we evaluate the time derivative ρ t as: The space derivative F x in (17) is computed using the minmod limiter: We compute the terms in (16) by the composite trapezoidal and mid-point rules, noting that the first and the last intervals are halved. Recalling that N := η/∆x we get ∂yF (ρ(t n , y), R(t n , y)) w(y − xj) dy = x j +η x j F (ρ(t n , y), R(t n , y)) w (y − xj) dy − F (ρ(t n , y), R(t n , y)) w(y − xj) F (ρ(t n , x j+k ), R(t n , x j+k ))w (k∆x).
Since the derived scheme uses alternating, staggered grids, we have to distinguish between the odd and even time steps. Formulas (14), (11)-(12), (15)-(16) describe the odd steps. The even steps are obtained by shifting the indexes in the aforementioned equations by 1 2 and the computational domain should be extended by ∆x is constant near the boundaries, so we use these constant values on an extended domain.
4. Numerical tests. In this section, we perform numerical simulations with different choices for the speed law v and the convolution kernel w η . More precisely, we consider the following velocity functions, see [11]: with v max = 1 and ρ max = 1, and the kernels w η ∈ C 1 ([0, η]; R + ): We recall that in [5] the authors considered the Greenshield model with n = 1 (i.e. linear decreasing velocity). For the tests presented below, the space domain is given by the interval [−1, 1], the space discretization mesh ∆x = 0.002 and η = 0.1, where not specified otherwise. Absorbing conditions are imposed at the boundaries. At the right boundary, we add N = η/∆x ghost cells and define ρ n j = ρ n 2 ∆x for every j = 2 ∆x +1, . . . , 2 ∆x +N , thus extending the solution constantly equal to the last value inside the domain. At the left boundary, we just need to add one ghost cell, as in classical problems.

4.1.
Monotonicity of the solution. In Figure 1 we compare the solutions of the Riemann problem with initial datum with constant kernel w η = 1/η. Numerical simulations confirm that the scheme generally is not monotonicity preserving, as shown in the case of the Greenberg's and the Underwood's velocities, which are non-linear, see the Appendix for more details on this aspect. As a consequence, the corresponding total variations ( Figure  2) are not constant.

4.2.
Limit η 0. We consider equation (1) equipped with Greenshield's velocity (n = 1) in the situation studied in [18, Section 3, Example 1]: a red traffic light is located at x = −0.1 and it turns green at the initial time t = 0, so the initial datum is: First of all, we compare the approximate solutions obtained with the Lax-Friedrichs and the Nessyahu-Tadmor schemes at time T = 0.5 for the mesh sizes ∆x = 0.01 and ∆x = 0.001, see Figure 3. In all cases, we remark the higher precision obtained by the Nessyahu-Tadmor scheme.  We now investigate the convergence to the solution of the classical conservation law ∂ t ρ + ∂ x (ρ(1 − ρ)) = 0 as η 0. To this end, we compare the solutions of (1) obtained using scheme (14) with different values of η = 0.1, 0.01, 0.001 with the local solution at time T = 0.5 for ∆x = 0.0001 (Figure 4). Contrarily to the results obtained in the case of the Arrhenius look-ahead model, the solutions to the non-local LWR model seem to converge to the solution of the classical LWR equation as the support of the convolution kernel shrinks to a point, as shown in Figure 4. In particular, we remark that the solutions of the non-local model (1) display sharp discontinuities in the density profile (one of which being a decreasing jump that could look unphysical, but satisfies the entropy condition (3)). On the contrary, the corresponding velocities are Lipschitz continuous, with a Lipschitz constant that depends on the L ∞ -norm of w η . Indeed, it is easy to show that there holds the estimate   We also compute numerical convergence orders and L 1 -errors of non-local solutions ρ η computed with scheme (14), θ = 1, for η 0, with respect to the local solutionρ computed with classical Lax-Friedrichs method as reference solution, as e(η) = ρ η (T, ·) −ρ(T, ·) L 1 , and γ(η) = log 10 e(η) e(η/10) , for η = 0.1, 0.01, 0.001 (see Table 1).
Appendix: Sketch of the proof of Proposition 2.
We observe here that the coefficient in (20d) is not positive in general, hence the scheme is not monotonicity preserving unless v = 0 as was the case in [5].
Details of the above calculations can be found in [12].