A Godunov type scheme for a class of LWR traffic flow models with non-local flux

1.   Introduction

Over recent years, non-local conservation laws gained growing interest for a wide field of applications such as supply chains [3], sedimentation [4], conveyor belts [12], crowd motion [8] or traffic flow [5,10]. In the latter case, the well-known Lighthill-Whitham-Richards (LWR) model [15,16] has been extended by considering non-local velocity terms depending on the downstream traffic so that drivers adapt their velocity to the mean traffic in front, see [5,10].

The well-posedness of special non-local flux problems has been investigated in for example [1,2,7,9]. However, only a few numerical schemes have been applied so far to solve these type of equations. The most common approach are first order LxF type schemes [2,4,5,7,10], while recently second- and higher-order schemes have been introduced [6,11]. We remark that also these higher-order methods rely on LxF type numerical flux functions, which imply the same drawbacks known from local conservation laws. Certainly, the LxF type scheme offers a powerful tool to numerically analyze non-local flux problems but typically leads to approximate solutions with strong diffusive behavior. As we are interested in a more accurate approach, we present a Godunov type scheme for a class of scalar conservation laws with non-local flux. In addition, by deriving several properties of the scheme, we prove the well-posedness for these special non-local conservation laws, which in contrast to other models [5,7,10] focus on a non-local mean velocity of the downstream traffic. Furthermore, the Godunov type scheme approach allows for physically reasonable solutions meaning that a maximum principle is satisfied and negative velocities as well as negative fluxes are avoided.

This work is organized as follows: In Section 2, we present the considered class of non-local conservation laws for traffic flow. Afterwards, we introduce the Godunov type scheme and derive important properties of the scheme such as $L^\infty$ and bounded variation (BV) estimates in Section 3. Those are also used to show the well-posedness of the proposed traffic model. In Section 4, we present numerical examples, which demonstrate the better accuracy of the Godunov type scheme in comparison to the widely used LxF type scheme and also provide a comparison to the model in [7].

2.   Modeling

We briefly present an already existing traffic flow model with non-local flux originally introduced in [5,7,10]. Based on the modeling ideas therein, we propose an adapted model and show its well-posedness. The key difference appears in the flux function, where instead of a mean downstream density a mean downstream velocity is considered.

2.1.   An existing model with mean downstream density

The model considered in [7] is given by a scalar conservation law of the form

where

For initial conditions

the existence and uniqueness of weak entropy solutions is stated in [7,Theorem 1] if the following hypotheses are satisfied:

Note that a whole family of kernel functions $w_\eta$ is considered to also analyze the limit behaviour of the model ($\eta\rightarrow0$ and $\eta\rightarrow\infty$).

2.2.   Model considering a mean downstream velocity

The non-local model (1) to (3) can be applied in the context of traffic flow and chooses the velocity based on a mean downstream traffic density. In contrast to this approach, it is also reasonable to assume that drivers adapt their speed based on a mean downstream velocity, anticipating the future space in front of them, see Figure 1.

Therefore, we consider a slightly different model compared to (1) to (3), namely

where

and we have given initial conditions

For simplicity, let us also define

In (4) and (5), we assume the same hypotheses as for (1) and (2) and one additional restriction:

As we will see in Section 3, the reformulation of the original model (1) to (3) keeps the main properties and allows for a straightforward application of a Godunov type scheme.

Remark 1. We note that in the case of a linear velocity function $v(\rho)$, e.g. $v(\rho) = 1-\rho$, the model given by (4) and (5) coincides with (1) and (2).

The weak entropy solutions of problem (4) to (6) are intended in the following sense:

Definition 2.1. [13,Definition 1] A function $\rho \in (L^1\cap L^\infty \cap {\rm{BV}})(\mathbb{R^+}\times\mathbb{R};\mathbb{R})$ is a weak entropy solution if

for all $\phi \in C_c^1(\mathbb{R}^2;\mathbb{R}^+)$ and $\kappa \in I = [0, \rho_{\max}]$.

Our main result concerning the new model is given by the following theorem, which states the well-posedness of problem (4) to (6).

Theorem 2.2. Let $\rho_0 \in {\rm{BV}}(\mathbb{R};[0, \rho_{\max}])$ and hypotheses (H2) hold. Then, the Cauchy problem

admits a unique weak entropy solution in the sense of Definition 2.1 and

The proof consists of two parts: existence and uniqueness of entropy solutions. While the uniqueness proof follows from the Lipschitz continuous dependence of weak entropy solutions on the initial data, the existence proof is based on a construction of a converging sequence of approximate solutions defined by a Godunov type scheme.

2.3.   Uniqueness of entropy solutions

One part of the proof to Theorem 2.2 is to show uniqueness of entropy solutions for the model (4) to (6). Therefore, we prove the Lipschitz continuous dependence of weak entropy solutions with respect to the initial data. Here, we follow [5,7,10] and use Kruzkov's doubling of variables technique [13]. Note that in the following $\Vert \cdot \Vert$ denotes $\Vert \cdot \Vert_{L^\infty}$.

Theorem 2.3. Under hypotheses (H2), let $\rho$ and $\sigma$ be two entropy solutions of (4) to (6) with initial data $\rho_0$ and $\sigma_0$, respectively. Then, for any $T>0$, there holds

with $K$ given by (12).

Proof. The functions $\rho$ and $\sigma$ are weak entropy solutions of

respectively, and $V$, $U$ are bounded measurable functions and Lipschitz continuous w.r.t. $x$ since $\rho, \sigma\in (L^1\cap L^\infty \cap BV)(\mathbb{R}^+\times \mathbb{R};\mathbb{R})$. Using the classical doubling of variables technique, we get the following inequality:

where $\rho_x$ must be understood in the sense of measures. Applying the mean value theorem and using the properties of the kernel function, we deduce

Using the Leibniz integral rule and again the mean value theorem, we can also obtain for a.e. $x\in\mathbb{R}$

If we plug (10) and (11) into (9), we obtain

with

By Gronwall's lemma we get (8) and for $\sigma_0 = \rho_0$ the uniqueness of entropy solutions.

The existence of weak entropy solutions is now proved in Section 3. We therefore introduce a Godunov type scheme used to construct approximate solutions.

3.   A Godunov type scheme

The main new contribution of this work is to develop a suitable Godunov type numerical scheme for the non-local model (4) to (6). We derive $L^\infty$ and BV bounds for the approximate solutions and further provide, due to a discrete entropy inequality, all ingredients to finally prove the well-posedness result Theorem 2.2.

3.1.   Numerical scheme

We take a space step $h$ such that for the size of the (downstream) kernel $\eta = Nh$ with $N\in \mathbb{N}$ holds. The time step is denoted by $\tau$ and we set $\lambda = \tau/h$. We define the space grid by $x_{j+\frac{1}{2}} = (j+\frac{1}{2})h$ being the cell interfaces and $x_j = jh$ the cell centers for $j\in \mathbb{Z}, $ see Figure 2. The finite volume approximate solution $\bar{\rho}(t, x)$ is denoted by $\rho_j^n$ for $(t, x)\in [t^n, t^{n+1}[\times ]x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]$.

Within the proposed scheme, we intend to mimic the numerical flux function of the Godunov scheme for local conservation laws, i.e., minimizing or maximizing the flux within $[\rho_{j}^n, \rho_{j+1}^n]$ or $[\rho_{j+1}^n, \rho_{j}^n]$ depending on whether $\rho_j^n \le \rho_{j+1}^n$ or $\rho_{j+1}^n \le \rho_{j}^n$, respectively. As the convolution term $V(t, x)$ is a non-local velocity function (see (7)), there is no straightforward way of adapting this result to the non-local case. Therefore, we first examine the flux at the interface $x_{j+\frac{1}{2}}$, where the convolution term for $t\in [t^n, t^{n+1}[$ is then approximated by

with

which is motivated by

Remark 2. For (14) we follow [4,Equation (3.2)] to satisfy

if $W_0 = 1$ and $\gamma_k$ is computed exactly or by an appropriate quadrature formula (such that all $\gamma_k$ are non-negative and $\sum_{k = 0}^{N-1} \gamma_k = W_0 = 1$ holds).

An example for the computation of $V_{j+\frac{1}{2}}^n$ can be seen in Figure 2. The corresponding values of $\rho$, which are used in the computation of $V_{j+\frac{1}{2}}^n$ for the case $N = 2$, are gray-shaded.

Based on the approximate convolution term $V_{j+\frac12}^n$, we adapt the numerical flux function of the Godunov scheme as follows:

Summarizing, the entire Godunov type scheme is initialized by the initial data $\rho_0$ as

and can be computed using (13) and (16) by the finite volume scheme

Remark 3. Analogously to the scheme (17) and (18) with (13), a Godunov type scheme for the model (1) to (3) considered in [7] can be derived. For this, similar properties can be shown analogously to the following sections. One major advantage of our Godunov type scheme unlike the LxF type scheme used in [5,7,10] is that the numerical fluxes $F_{j+\frac{1}{2}}^n = V_{j+\frac{1}{2}}^n g(\rho_j^n)$ are always non-negative.

3.2.   Maximum principle

The approximate solutions constructed by the Godunov type scheme (18) satisfy a strict maximum principle:

Theorem 3.1. Let hypotheses (H2) hold. For a given initial datum $\rho^0_j, \ j\in\mathbb{Z}$ with $\rho_M = \sup_{j\in\mathbb{Z}} \rho^0_j$ and $\rho_m = \inf_{j\in\mathbb{Z}} \rho^0_j$, the approximate solutions constructed by the scheme (18) satisfy the bounds

if the following Courant-Friedrichs-Levy (CFL) condition holds:

Proof. We prove the claim per induction. For $n = 0$ the claim is obvious, so we suppose

holds for a fixed $n\in \mathbb{N}$.

Before considering $\rho_j^{n+1}$ we show some general inequalities:

where we used the monotonicity of $w_\eta$ and $v$, and the mean value theorem.

Analogously, we obtain from (20)

By multiplying inequality (21) by $g(\rho_M)$, subtracting $V_{j+\frac{1}{2}}^n g(\rho_j^n)$ and applying the mean value theorem, we obtain

Therefore, under the CFL condition (19), we have

Analogously, we obtain

to show

which gives us the claim.

The maximum principle ensures that the numerical solution to (4) to (6) is bounded from above by $\rho_M = \sup_{j\in\mathbb{Z}} \rho^0_j \in [0, \rho_{\max}]$ and hence does not exceed the maximal density $\rho_{\max}$. In addition, the scheme is positivity preserving as the solution stays non-negative, since $\rho_m = \inf_{j\in\mathbb{Z}} \rho^0_j\geq 0$.

3.3.   BV estimates

Next, we derive a BV estimate for the approximate solutions constructed by the Godunov type scheme (18). Similar to the LxF type scheme analyzed in [5,7,10], BV estimates cannot be derived using the standard general approaches. In particular, the Godunov type scheme also does not fit into the classical assumptions of total variation diminishing (TVD) schemes, as the total variation may slightly increase (as it is the same for the analytical solution). Nevertheless, the numerical scheme has a bounded total variation. Further, to finally also prove the existence of solutions to the model (4) to (6), we also need to provide a bound on the (discrete) variation in space and time. We begin with the BV estimate in space:

Theorem 3.2. Let hypotheses (H2) hold, $\rho_0\in BV(\mathbb{R};[0, \rho_{\max}])$ and let $\bar \rho$ be given by (18). If the CFL condition (19) holds, then for every $T>0$ the following discrete space BV estimate is satisfied:

with $C(w_\eta, v, g, \rho_{\max}) = w_\eta(0)( \Vert v'\Vert \Vert g \Vert \rho_{\max} + \Vert v\Vert \Vert g'\Vert)$.

Proof. Let us define

Then, we obtain

where $\xi_{j+\frac{1}{2}}$ is between $\rho_{j}^n$ and $\rho_{j+1}^n$. With (20) we derive

where $\zeta_{j+\frac{1}{2}}$ is again between $\rho_{j}^n$ and $\rho_{j+1}^n$. Thus, we have

Due to the CFL condition (19) and hypotheses (H2), all terms before the differences are positive and we get

Rearranging the indices we obtain

Using inequality (21), for which

holds, and with $\gamma_0\leq hw_\eta(0)$, we obtain

Therefore, we recover the following estimate for the total variation

We are now able to also provide an estimate for the discrete total variation in space and time:

Theorem 3.3. Let hypotheses (H2) hold, $\rho_0\in BV(\mathbb{R};[0, \rho_{\max}])$ and let $\bar \rho$ be given by (18). If the CFL condition (19) holds, then for every $T>0$ the following discrete space and time total variation estimate is satisfied:

with $C(w_\eta, v, g, \rho_{\max}) = w_\eta(0)( \Vert v'\Vert \Vert g \Vert \rho_{\max} + \Vert v\Vert \Vert g'\Vert)$.

Proof. We fix $T\in \mathbb{R}^+$. If $T\leq \tau$, then $TV(\bar{\rho}; \mathbb{R}\times [0, T])\leq T\cdot TV(\rho_0)$. For $T>\tau$ let $M\in \mathbb{N}\setminus \{0\}$ such that $M\tau < T\leq (M+1) \tau$. Then

If we consider the scheme (18), we obtain

Taking absolute values yields

Summing over $j$ and rearranging the indices gives us

so that we have

Therefore, we recover

as desired.

3.4.   Discrete entropy inequality

As another desirable property and final ingredient regarding the proof of Theorem 2.2, we next show that the approximate solutions obtained by the Godunov type scheme (18) fulfill a discrete entropy inequality. Therefore, we follow [2,5,7,10] and define

with $ a\wedge b = \max(a, b)$ and $ a\vee b = \min(a, b)$.

Theorem 3.4. Let $\rho_j^n, \ j\in\mathbb{Z}, \ n\in\mathbb{N}$ be given by (18), and let the CFL condition (19) and hypotheses (H2) hold. Then we have

for all $j\in\mathbb{Z}, \ n\in\mathbb{N}$ and $\kappa\in I = [0, \rho_{\max}]$.

Proof. The proof closely follows [2,5,7]. We set

which is a monotone non-decreasing function with respect to each variable under the CFL condition (19) since we have

Moreover, we have the identity

Then, by monotonicity, the definition of the scheme (18) and by using (for the last inequality) the non-negativity of $(a, b)\mapsto ({\rm {sgn}}(a+b)-{\rm {sgn}}(a))(a+b)$, we get (24):

3.5.   Proof of Theorem 2.2

Since we have already shown uniqueness of weak entropy solutions to the model (4) to (6), it remains to finalize the existence proof. Similar to [5,Section 4] and [7,Proof of Theorem 1], the convergence of the approximate solutions constructed by the Godunov type scheme (18) towards the unique weak entropy solution can be proven by applying Helly's theorem. The latter can be applied due to Theorems 3.1 and 3.3 and states that there exists a sub-sequence of the constructed $\bar \rho$ that converges to some $\rho \in (L^1 \cap L^\infty \cap BV)(\mathbb{R}^+ \times \mathbb{R};I)$. Following a Lax-Wendroff type argument, see [14,Theorem 12.1], one can show that the limit function $\rho$ is a weak entropy solution of (4) to (6) in the sense of Definition 2.1. Together with the uniqueness result in Theorem 2.3, this concludes the proof of Theorem 2.2.

4.   Numerical examples

In this section, we present some numerical examples demonstrating the advantages of the Godunov type scheme in comparison to the widely used LxF type scheme. The latter will be briefly introduced in the following section. In addition, we also comment on the differences between the model considered in this work (4) to (6) and the earlier one (1) to (3).

4.1.   A LxF type scheme

The common scheme used so far for the problem (1) to (3) is a LxF type scheme, where the downstream velocity of the convolution term is computed by

with $w_\eta^k = w_\eta(kh)$. This discretization of $w_\eta$ implies (for $W_0 = 1$)

Thus, the approximation (25) of the convolution term slightly overestimates the traffic density and therewith underestimates the velocity. Further, unphysical densities beyond $\rho_{\max}$ or negative velocities are possible. This can be both avoided by discretizing the kernel function as proposed in (14). In the following subsections we will use this discretization to avoid that the accuracy studies are biased by different quadrature rule errors.

The numerical flux function of the LxF scheme is given by

with $\alpha\geq 0$ being the viscosity coefficient. This leads to the scheme

For the corresponding CFL condition and restrictions on $\alpha$, we refer to [7,Proposition 2].

Remark 4. Note that the LxF type scheme can be adapted to the model (4) and vice versa the Godunov type scheme to model (1), where in comparison the LxF type scheme adds more diffusion to the numerical solution (see also Section 4.2).

4.2.   Accuracy of the Godunov type scheme

In the following sections, we consider the solution of model (4) to (6) and model (1) to (3) in the case of traffic flow. So we have $g(\rho) = \rho$, a velocity function $v(\rho)$ specified below, and a road of length $L = 1$. We are interested in the solution at a certain final time $T$ for the initial conditions

For simplicity, we use periodic boundary conditions in all our examples. To compute the $L^1$ error of an approximate solution $\rho_h$ with step size $h$ compared to a reference solution $\rho_{\tilde{h}}$ with step size $\tilde h$ at time $T$, we apply

Now, in order to compare the accuracy of the Godunov scheme and the LxF scheme, we first use the linear velocity function

as for this choice both models coincide (see Remark 1) and the different discretization schemes both have been well analyzed. We consider the final time $T = 0.1$ and the kernel function $w_\eta(x) = 3(\eta^2-x^2)/(2\eta^3)$ with $\eta = 0.1$. We compare the $L^1$ distances of the LxF and the Godunov type scheme to a reference solution computed with the LxF type scheme and $\tilde h = 0.02\cdot 2^{-9}$. The spatial step size is given by $h = 0.02\cdot 2^{-n}$ with $n\in\{0, \dots, 6\}$. The time step parameter $\tau$ is given by the minimum of the CFL condition (19) and the CFL condition of the LxF type scheme.

The results for this first test case are given in Table 1. Obviously, the $L^1$ errors of the Godunov type scheme are significantly smaller than the ones of the LxF scheme.

The better accuracy of the Godunov type scheme can also be seen directly in Figure 3. We notice that in the presence of the two discontinuities, the Godunov type scheme in particular shows a better resolution of the solution structure than the LxF scheme at the left-hand side, while the resolution for the rarefaction wave close to the jump on the right-hand side is quite similar for the short time period $T = 0.1$.

If we consider a longer time period, i.e. $T = 1$, the good performance of the Godunov type scheme can be observed in Figure 4. Here, the reference solution is computed with the LxF type scheme with a spatial step size of $\tilde h = 0.02\cdot 2^{-7}$.

In addition, we consider the accuracy for a non-linear velocity function,

We choose the constant kernel $w_\eta(x) = \frac{1}{\eta}$ with $\eta = 0.1$ and consider again the solution of the initial conditions (27) but at time $T = 0.05$ on a road of length $L = 1$. The spatial step size is given as above by $h = 0.02\cdot 2^{-n}$ with $n\in\{0, \dots, 6\}$. The reference solution is computed by the LxF scheme adapted to the model (4) and with $\tilde h = 0.02\cdot 2^{-9}$.

The results for the non-linear test case can be seen in Table 2. Similar to the linear test case, the $L^1$ errors of the Godunov type scheme are significantly smaller than the ones of the LxF scheme. In addition, the better resolution of the Godunov type scheme can be seen in Figure 5. Again, similar to the linear case, the Godunov type scheme in particular shows a better resolution of the solution structure than the LxF scheme at the left-hand side, while the resolution for the rarefaction wave close to the jump on the right-hand side is quite similar.

4.3.   Comparison of the models

Next, we aim to discuss the differences between the models (4) to (6) and (1) to (3). To see the different dynamics within the two models, we have to choose a non-linear velocity function and we choose the same non-linear velocity $v(\rho) = 1-\rho^5$ as before with the same parameter $\eta = 0.1$ for the constant kernel function and final time $T = 0.05$.

For a fair comparison of the evolution of densities, we apply the Godunov type scheme to both models. Figure 6 shows the results for a spatial step size $h = 0.01$ and a time step size $\tau$ determined by the CFL condition. It can be observed that the approximate solution of the earlier model (1) to (3) contains rather large oscillations while resolving the traffic jam. In contrast to that, the model (4) to (6) resolves the traffic jam in a more monotone way.

4.4.   Limit $\eta\rightarrow 0$

Finally, we take a look at the limit case $\eta\rightarrow0$ and investigate numerically whether the approximate solutions constructed by the proposed Godunov type scheme converge towards the solution of the (local) LWR traffic model. Note that this property is far from obvious since the constants in Theorems 2.3 and 3.3 blow up (see also [7]). For a numerical investigation, we consider the same (non-linear) scenario as above with a fixed space step size $h = 0.5\cdot 10^{-4}$ and vary $\eta \in \{10^{-1}, $ $10^{-2}, $ $10^{-3}, $ $10^{-4}\}$. As final time we take $T = 0.05$.

To evaluate the convergence, we compute the $L^1$ distances between the approximate solutions obtained for the proposed Godunov type scheme applied to (4) to (6) and the result of a classical Godunov scheme for the corresponding local LWR problem. Obviously, the corresponding $L^1$ distances shown in Table 3 demonstrate the convergence towards the solution of the local traffic model. The results are further illustrated in Figure 7.

5.   Conclusion

In this work, we have presented a Godunov type scheme for a class of non-local conservation laws. For this novel scheme we provide $L^\infty$ and BV bounds as well as a discrete entropy inequality. Based on these results, we also proved the well-posedness, i.e., existence and uniqueness of weak entropy solutions. The proposed Godunov type scheme can be adapted to other classes of non-local conservation laws and is very promising as it adds less numerical diffusion to the solution compared to the LxF type scheme. Certainly, this advantage can also be exploited by using the underlying numerical flux function within higher-order methods.

In future work we aim at constructing several higher order methods based on the presented Godunov type scheme. In addition, the considered model with mean downstream velocity may be advantageous in the context of networks. Here we aim at investigating appropriate coupling conditions and suitable discretization schemes.