On the implementation of a finite volumes scheme with monotone transmission conditions for scalar conservation laws on a star-shaped network

TL;DR

This paper validates a finite volume numerical scheme with monotone transmission conditions for scalar conservation laws on star-shaped networks by comparing with explicit solutions and Riemann problem solutions, supporting its accuracy and applicability.

Contribution

It provides the first validation of a finite volume scheme with monotone transmission conditions for networked scalar conservation laws through explicit and Riemann problem comparisons.

Findings

Scheme accurately reproduces explicit solutions

Effective in solving Riemann problems on networks

Supports further validation of the finite volume method

Abstract

In this paper we validate the implementation of the numerical scheme proposed in [3]. The validation is made by comparison with an explicit solution here obtained, and the solutions of Riemann problems for several networks. We then perform some simulations in order to qualitatively validate the model under consideration. Such results represent also a first step for the validation of the finite volumes scheme introduced in [9].

Tables

Table 1. Table 1 . Relative 𝐋 1 superscript 𝐋 1 \mathbf{L}^{1} -error at time t = 2.4 𝑡 2.4 t=2.4 computed in Section  5 .

Number		Rate of		Rate of		Rate of
of cells	$E_{{𝐋}^1}^{s,\Gamma }$	convergence	$E_{{𝐋}^1}^{s,𝖨}$	convergence	$E_{{𝐋}^1}^{s,𝖩}$	convergence
per arc
$60$	$6.5374\times {10}^{-2}$	-	$2.1928\times {10}^{-1}$	-	$1.4155\times {10}^{-2}$	-
$120$	$3.4281\times {10}^{-2}$	$0.9313$	$1.1380\times {10}^{-1}$	$0.9464$	$7.8554\times {10}^{-3}$	$0.8496$
$600$	$7.6754\times {10}^{-3}$	$0.9302$	$2.4933\times {10}^{-2}$	$0.9441$	$1.9468\times {10}^{-3}$	$0.8625$
$1200$	$4.8890\times {10}^{-3}$	$0.8800$	$1.6393\times {10}^{-2}$	$0.8831$	$1.0579\times {10}^{-3}$	$0.8660$
$6000$	$1.9875\times {10}^{-3}$	$0.7721$	$7.1933\times {10}^{-3}$	$0.7575$	$2.5294\times {10}^{-4}$	$0.8738$
$12000$	$1.6804\times {10}^{-3}$	$0.7034$	$6.3143\times {10}^{-3}$	$0.6815$	$1.3587\times {10}^{-4}$	$0.8774$

Equations

$${\rho }_{h,t}+f_h({\rho }_h{)}_x=0,\text{for }t>0,x\in {\Omega }_h,h\in \mathsf{H},$$

$${\rho }_h(0,x)={\rho }_{h,0}(x),x\in {\Omega }_h,$$

$$\sum _{i\in \mathsf{I}}{f}_i\left({\rho }_i(t,0^{-})\right)=\sum _{j\in \mathsf{J}}{f}_j\left({\rho }_j(t,0^{+})\right).$$

$$\int_{0}^{\infty}\int_{\Omega_{h}}\Bigl{(}|\rho_{h}-k|\partial_{t}\phi+\mathrm{sign}\left(\rho_{h}-k\right)\left(f_{h}(\rho_{h})-f_{h}(k)\right)\partial_{x}\phi\Bigr{)}\,{\,\rm{d}}x{\,\rm{d}}t+\int_{\Omega_{h}}\left|\rho_{h,0}(x)-k\right|\phi(0,x)\,{\,\rm{d}}x\geq 0;$$

$$\sum _{i\in \mathsf{I}}{f}_i({\rho }_i(t,0-))=\sum _{j\in \mathsf{J}}{f}_j({\rho }_j(t,0+)).$$

$$\{\begin{array}{ll}{\partial }_t\rho +{\partial }_{x}f(\rho )=0,& t>0,x\in \text{mathbbm}R,\\ \rho (0,x)=\{\begin{array}{ll}a& \text{if }x<0,\\ b& \text{if }x>0,\end{array}& x\in \text{mathbbm}R,\end{array}$$

$$G(a,b)=\{\begin{array}{ll}{\min }_{s\in [a,b]}\{f(s)\}& \text{if }a\le b,\\ {\max }_{s\in [b,a]}\{f(s)\}& \text{if }a\ge b.\end{array}$$

$$\{\begin{array}{l}{\partial }_t\rho +{\partial }_{x}f(\rho )=0,\text{for }(t,x)\in \text{mathbbm}{R}_{+}\times \text{mathbbm}{R}_{-},\\ \rho (t,0)={\rho }_b(t),\\ \rho (0,x)={\rho }_0(x),\end{array}$$

$$\{\begin{array}{ll}{\partial }_t{\rho }_h+{\partial }_x{f}_h({\rho }_h)=0,& \text{on }]0,T[\times {\Omega }_h,\\ {\rho }_h(t,0)=v_h(t),& \text{on }]0,T[,\\ {\rho }_h(0,x)={\rho }_{h,0}(x),& \text{on }{\Omega }_h,\end{array}$$

$$f_h\left({\rho }_h(t,0-)\right)=G_h\left({\rho }_h(t,0-),v_h(t)\right),\text{if }h\in \mathsf{I};$$

$$f_h\left({\rho }_h(t,0+)\right)=G_h\left(v_h(t),{\rho }_h(t,0+)\right),\text{if }h\in \mathsf{J}.$$

$$v_h(t)=p(t),\text{for a.e. }t\in \text{mathbbm}{R}_{+},\text{for all }h\in \mathsf{H}.$$

$$\sum _{i\in \mathsf{I}}{G}_i\left({\rho }_i(t,0^{-}),p(t)\right)=\sum _{j\in \mathsf{J}}{G}_j\left(p(t),{\rho }_j(t,0^{+})\right),\text{for a.e. t>0}.$$

$$\sum _{i=1}^m{\int }_{-(M-Lt)}^0|{\rho }_i(t,x)-{\rho }_i^{\prime }(t,x)|dx+\sum _{j=m+1}^{m+n}{\int }_0^{M-Lt}|{\rho }_j(t,x)-{\rho }_j^{\prime }(t,x)|dx$$

$$\le \sum _{i=1}^m{\int }_{-M}^0|{\rho }_{i,0}(x)-{\rho }_{i,0}^{\prime }(x)|dx+\sum _{j=m+1}^{m+n}{\int }_0^M|{\rho }_{j,0}(x)-{\rho }_{j,0}^{\prime }(x)|dx.$$

$$f(\rho )\equiv f_h(\rho )=\rho (1-\rho )$$

$${\rho }_{1,0}(x)={\chi }_{[-1/2,0]}(x),{\rho }_{2,0}(x)=3/4{\chi }_{[-1/4,0]}(x),{\rho }_{3,0}(x)=0.$$

$$G_1(1,p_1)+G_2(3/4,p_1)=G_3(p_1,0),$$

$${\mathfrak{R}}_{O,1}=\frac{1}{2}\left(1-\frac{x}{t}\right),\text{for }-t\le x\le -\frac{\sqrt{2}}{2}t.$$

$${\mathfrak{S}}_{O,2}:\dot{x}(t)=\sigma \left(\frac{3}{4},p_1\right),x(0)=0.$$

$${\mathfrak{R}}_{O,3}=\frac{1}{2}\left(1-\frac{x}{t}\right),\text{for }0\le x\le t.$$

$${\mathfrak{S}}_{C,2}:\dot{x}(t)=\sigma \left(0,p_1\right),x(t_C)=x_C,$$

$${\mathfrak{S}}_{D,1}:\dot{x}(t)=\sigma \left(0,{\mathfrak{R}}_{O,1}\left(t,x(t)\right)\right),x(1/2)=-1/2.$$

$${\mathfrak{S}}_{E,1}:\dot{x}(t)=\sigma \left(0,p_1\right),x(t_E)=x_E.$$

$$G_1\left(p_1,p_2\right)+G_2\left(0,p_2\right)=G_3\left(p_2,\frac{1}{2}\right),\text{for }t>t_F.$$

$${\mathfrak{R}}_{F,1}(t,x)=\frac{1}{2}\left(1-\frac{x}{t-t_F}\right),\text{for}-\frac{\sqrt{2}}{2}(t-t_F)<x\le 0.$$

$${\mathfrak{S}}_{H,3}:\dot{x}(t)=\sigma \left(0,{\mathfrak{R}}_{O,3}\left(t,x(t)\right)\right),x(t_H)=0.$$

$$\bigcup _{\ell \le -1}(x_{\ell }^i,x_{\ell +1}^i),i\in \mathsf{I},\bigcup _{\ell \ge 0}(x_{\ell }^j,x_{\ell +1}^j),j\in \mathsf{J},$$

$$\Delta t\underset{h}{\max }\{L_h\}\le \frac{\Delta x}{2},$$

$${\rho }_{\ell +\frac{1}{2}}^{h,0}=\frac{1}{\Delta x}{\int }_{x_{\ell }^h}^{x_{\ell +1}^h}{\rho }_h^0(x)dx,$$

$$p^s\in [0,{\rho }_{\max }]\text{s.t. }\sum _{i\in \mathsf{I}}{G}_i({\rho }_{-\frac{1}{2}}^{i,s},p^s)=\sum _{j\in \mathsf{J}}{G}_j(p^s,{\rho }_{\frac{1}{2}}^{j,s}),$$

Full text

Such results represent also a first step for the validation of the finite volumes scheme introduced in [9].

The author acknowledges the support of the Région Bourgogne Franche-Comté, projet 2017-2020 “Analyse mathématique et simulation numérique d’EDP issus de problèmes de contrôle et du trafic routier" and of the Université de Bourgogne Franche-Comté, projet Chrysalide 2017 “Contrôle, analyse numérique et applications d’équations hyperboliques sur un réseau". The author is member of the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM)

1. Introduction

We investigate from the numerical point of view the model developed in [3] by Andreianov, Coclite and Donadello, called here ACD, which describes the evolution of traffic at a junction consisting of $m$ incoming and $n$ outgoing arcs. Incoming arcs are parametrized by $x\in\mathbbm{R}_{-}$ and numbered by the index $i\in\mathsf{I}=\{1,\ldots,m\}$, while outgoing arcs are parametrized by $x\in\mathbbm{R}_{+}$ and numbered by the index $j\in\mathsf{J}=\{m+1,\ldots,m+n\}$ in such a way that the junction is always located at $x=0$. We denote the generic arc by $\Omega_{h}$, $h\in\mathsf{H}=\{1,\ldots,m+n\}$, and the network by $\Gamma=\Pi_{h\in\mathsf{H}}\Omega_{h}$.

We describe the evolution of traffic on each arc by the Lighthill-Whitham-Richards (LWR) model [13, 15], namely by a scalar conservation law of the form

[TABLE]

where $\rho_{h}$ is the density and $f_{h}$ is the flux on the $h$-th arc. We assume that the arcs have a common maximal density $\rho_{\max}>0$ and the fluxes are all bell-shaped (unimodal), Lipschitz and non-degenerate nonlinear i.e.

1. (F)

for all $h\in\mathsf{H}$, $f_{h}\in\mathbf{Lip}\left([0,\rho_{\max}];\mathbbm{R}_{+}\right)$ with $\|f_{h}^{\prime}\|_{\infty}\leq L_{h}$, $f_{h}(0)=0=f_{h}(\rho_{\max})$, and there exists $\rho_{h,c}\in\left(0,\rho_{\max}\right)$ such that $f_{h}^{\prime}(\rho)\leavevmode\nobreak\ (\rho_{h,c}-\rho)>0$ for a.e. $\rho\in[0,\rho_{\max}]$, 2. (NLD)

for all $h\in\mathsf{H}$, $f^{\prime}_{h}$ is not constant on any non-trivial subinterval of $[0,\rho_{\max}]$.

We augment (1) with the initial conditions

[TABLE]

where $\rho_{h,0}\in\mathbf{L}^{\infty}(\Omega_{h};[0,\rho_{\max}])$, $h\in\mathsf{H}$. We also impose the conservation of the total density at the junction, i.e. for a.e. $t\in\mathbbm{R}_{+}$

[TABLE]

Notice that the previous equation makes sense as the assumption (F) ensures the existence of strong traces [14, 16].

In [3] the authors prove the well-posedness of solutions obtained as vanishing viscosity limits for the Cauchy problem (1)-(2). Their result relies upon the explicit characterization of the class of admissible weak solution at the junction in terms of vanishing viscosity germ, see [1, 4, 5]. It represents a generalized study of the model investigated in [8], in which the authors establish the existence of weak solutions as limit of vanishing viscosity approximations. Such results are relevant in the perspective of a theoretical analysis of PDEs on networks, in particular, they allow the extension of the analogy between vanishing viscosity and numerical scheme both in the network case. Furthermore, their analysis is applicable to general junction solvers enjoying enjoying the order-preservation property, see for instance [10].

The aim of this paper is to validate the implementation of the numerical scheme proposed in [3] by comparison with an explicit solution here computed and with the solutions of Riemann problems both for merge, divide and 2-2 cases. Moreover, we show the consistency of the scheme with respect to the case of a network with no discontinuity at the junction, namely taking the same flux on each arc, and finally a convergence analysis is also performed. These results are used in the validation of the finite volumes scheme with point constraints at the junction introduced in [9].

The paper is organized as follows. In Section 2 we briefly recall the main theoretical results for ACD. In Section 3 we compute an explicit solution for the problem in the case of a merge consisting of two incoming and one outgoing arcs. In Section 4 we present the numerical scheme. Finally, Section 5 is devoted to the validation of the scheme.

2. Well-posedness of ACD in the frame of admissible solutions

The well-posedness for the general Cauchy problem (1)- (2)- (3) is established in [3] in the frame of admissible solution.

We recall some definitions.

Definition 2.1**.**

A function $\rho_{h}\in\mathbf{L}^{\infty}(\mathbbm{R}_{+}\times\Omega_{h};[0,\rho_{\max}])$, $h\in\mathsf{H}$, is a weak solution of (5) if

• •

for every $k\in[0,\rho_{\max}]$ and every nonnegative test function $\phi\in\mathcal{C}^{\infty}(\mathbbm{R}\times\Omega_{h};\mathbbm{R})$ with compact support

[TABLE]

• •

for a.e. $t>0$, it holds

[TABLE]

We remind the formulation of the Bardos-LeRoux-Nédélec boundary condition for conservation laws in terms of the Godunov numerical flux (see [7, 11]), which will be useful for the definition of admissible solution at the junction.

Definition 2.2**.**

The Godunov flux related to a flux $f$ satisfying (F) is the function which associates to any couple $(a,b)\in[0,\rho_{\max}]^{2}$ the value $f(\rho(t,0^{-}))=f(\rho(t,0^{+}))$ (the equality holds due to the Rankine-Hugoniot condition), where $\rho$ is the Kruzhkov [12] entropy solution to the Riemann problem

[TABLE]

see Figure 1. The analytical expression of the Godunov flux is given by

[TABLE]

We denote by $G_{h}$ the Godunov flux associated with the flux $f_{h}$, $h\in\mathsf{H}$.

Consider the initial boundary value problem

[TABLE]

and assume $\rho$ is a Kruzkov entropy solution in the interior of the half plane $\mathbbm{R}_{+}\times\mathbbm{R}_{-}$. Then $\rho$ satisfies the boundary condition in the sense of Bardos-LeRoux-Nédélec (see [7]) if and only if $f(\rho(t,0-))=G(\rho(t,0-),\rho_{b}(t))$.

Fix an initial condition $\vec{\rho}_{0}=\left(\rho_{1,0},\cdots,\rho_{m+n,0}\right)\in\mathbf{L}^{\infty}(\mathbbm{R}_{+}\times\Gamma;[0,\rho_{\max}]^{m+n})$. We look for a function $\vec{\rho}=\left(\rho_{1},\cdots,\rho_{m+n}\right)$ such that for every $h\in\mathsf{H}$, $\rho_{h}\in\mathbf{L}^{\infty}(\mathbbm{R}_{+}\times\Omega_{h},[0,\rho_{\max}])$ is a weak entropy solution of

[TABLE]

where $\vec{v}:\mathbbm{R}_{+}\to[0,\rho_{\max}]^{m+n}$ is to be fixed in the sequel in order to ensure the conservation at the junction.

The condition (4) is equivalent to ask for the traces $\rho_{h}(t,0\pm)$ to satisfy the boundary condition in the sense of Bardos-LeRoux-Nédélec

[TABLE]

In order to describe the solution of (1) we postulate (see [1, 6])

[TABLE]

The criterion for the choice of $p$ is equivalent to the condition (4), indeed, due to (6) and (7), we can express (4) in the following way

[TABLE]

We can now give the definition of admissible solution.

Definition 2.3**.**

Given an initial condition $\vec{\rho}_{0}\in\mathbf{L}^{\infty}(\Gamma;[0,\rho_{\max}]^{m+n})$, we say that $\vec{\rho}=\left(\rho_{1},\cdots,\rho_{m+n}\right)$ in $\mathbf{L}^{\infty}(\mathbbm{R}_{+}\times\Gamma;[0,\rho_{\max}]^{m+n})$ is an admissible solution for the Cauchy problem at the junction (1) associated with $\vec{\rho}_{0}$, if there exists a function $p\in\mathbf{L}^{\infty}(\mathbbm{R}_{+}:[0,\rho_{\max}])$ such that, for any $h\in\mathsf{H}$, $\rho_{h}$ is a weak solution of (5) with $v_{h}$, $h\in\mathsf{H}$ chosen to fulfill (8), and such that $\vec{\rho}$, $p$ fulfill (4).

The authors provide a characterization of vanishing viscosity limits for the problem (1) in terms of $m+n$ Dirichlet problems on $\Omega_{h}$, $h\in\mathsf{H}$ coupled by a transmission condition at the junction. To this aim, they introduce the vanishing viscosity germ (see [4, 5, 1]), which can be identified by the set of all possible stationary admissible solutions to (1) that are constant on each road of $\Gamma$.

The main result of [3] is summarized in the following theorem.

Theorem 2.4** (Theorem 3.1 in [3]).**

For any given initial condition $\vec{\rho}_{0}=\left(\rho_{0,1},\ldots,\rho_{0,m+n}\right)$ in $\mathbf{L}^{\infty}(\Gamma;\mathbbm{R}^{m+n})$ the problem (1) admits a unique admissible solution $\vec{\rho}$ in $\mathbf{L}^{\infty}(\mathbbm{R}_{+}\times\Gamma;[0,\rho_{\max}]^{m+n})$.

Moreover, if $\vec{\rho}$ and $\vec{\rho}^{\prime}$ are two admissible solutions corresponding respectively to the initial condition $\vec{\rho}_{0}$ and $\vec{\rho}_{0}^{\prime}$, then for all $M>0$ and $t<M/L$, where $L=\max\left\{\|f_{h}^{\prime}\|_{\mathbf{L}^{\infty}([0,\rho_{\max}];\mathbbm{R})}\,|\,h=1,\ldots,m+n\right\}$,

[TABLE]

In particular, the map that associates to $\vec{\rho}_{0}$ the unique corresponding admissible profile $\vec{\rho}(t)$, is non-expansive with respect to the $\mathbf{L}^{1}$ distance for all $t>0$.

3. An explicit admissible solution at a merge

In this section we compute an explicit solution for the problem consisting of two incoming and one outgoing arcs. We consider

[TABLE]

as the flux for each arc. As initial condition, we choose

[TABLE]

The exact solution is obtained by an explicit analysis of the wave-front interactions, with computer assisted computation of the front slopes and interaction times.

At time $t=0$, let $p_{1}\approx 0.85$ be the solution of

[TABLE]

then, on $\Omega_{1}$ a rarefaction $\mathfrak{R}_{O,1}$ starts from $O(0,0)$ and its values are given by

[TABLE]

On $\Omega_{2}$ starts the backward shock $\mathfrak{S}_{O,2}$ given by

[TABLE]

On $\Omega_{3}$ a rarefaction starts from $O(0,0)$ and its values are given by

[TABLE]

On $\Omega_{2}$, let $C(x_{C},t_{C})$ be the point where the shock $\mathfrak{S}_{B,2}:\,x(t)=-\frac{1}{4}+\frac{t}{4}$ originated from $B(-1/4,0)$ interacts with the shock $\mathfrak{S}_{O,2}$. As a result, from $C$ starts a shock given by

[TABLE]

and reaches the junction $x=0$ at time $t=t_{F}=3/2$ that corresponds to the time at which the second incoming arc becomes empty. On $\Omega_{1}$, in $D(-1/2,1/2)$, the stationary shock $\mathfrak{S}_{A,1}$ originated from $A(-1/2,0)$ interacts with the rarefaction $\mathfrak{R}_{O,1}$. As a result, from $D$ starts a shock $\mathfrak{S}_{D,1}$ given by

[TABLE]

Let $E(x_{E},t_{E})$ be the intersection between $\mathfrak{S}_{D,1}$ and $x(t)=-(\sqrt{2}/2)\,t$. From this point starts a forward shock

[TABLE]

Let $p_{2}=1/2$ be the solution of

[TABLE]

Therefore, a rarefaction appears on $\Omega_{1}$:

[TABLE]

Let $G$ be the point where $\mathfrak{S}_{E,1}$ and $\mathfrak{R}_{F,1}$ interact. From this point starts a forward shock $\mathfrak{S}_{G,1}$, with left state $\rho=0$, which reaches the junction at time $t_{H}\approx 2.75$, then $\Omega_{1}$ is empty. Finally, on $\Omega_{3}$ at time $t_{H}$ starts a shock which interacts with the rarefaction $\mathfrak{R}_{O,3}$ generating the additional shock

[TABLE]

4. Finite volumes numerical scheme

We fix a constant space step $\Delta x$. For $\ell\in\mathbbm{Z}$ and $h\in\mathsf{H}$, we set $x_{\ell}^{h}=\ell\Delta x$. We define the cell centers $x_{\ell+\frac{1}{2}}^{h}=(\ell+\frac{1}{2})\Delta x$ for $\ell\in\mathbbm{Z}$ and consider the uniform spatial mesh on each $\Omega_{h}$

[TABLE]

so that the position of the junction $x=0$ corresponds to $x_{0}^{h}$ for each edge. Then we fix a constant time step $\Delta t$ satisfying the CFL condition

[TABLE]

and for $s\in\mathbbm{N}$ we define the time discretization $t^{s}=s\Delta t$. At each time $t^{s}$, $\rho_{\ell+\frac{1}{2}}^{h,s}$ represents an approximation of the main value of the solution on the interval $[x_{\ell}^{h},x_{\ell+1}^{h})$, $\ell\in\mathbbm{Z}$, along the $h$-th arc. We initialize the scheme by discretizing the initial conditions

[TABLE]

for all $h\in\mathsf{H}$ and for $\ell\leq-1$ if $h\in\mathsf{I}$, $\ell\geq 0$ if $h\in\mathsf{J}$.

For each $s\in\mathbbm{N}$, at all cell interfaces $x_{\ell}^{h}$ with $\ell\neq 0$ we consider the standard Godunov flux $G_{h}$ corresponding to the flux $f_{h}$. At the junction $x_{0}^{h}$ we take on each arc $\Omega_{h}$ the Godunov flux corresponding to the admissible solution of the Riemann problem at the junction, defined as in [3], which we compute by means of a two-step procedure:

1. (i)

find

[TABLE] 2. (ii)

compute

[TABLE]

where

[TABLE]

The choice of the Godunov’s flux is motivated by the fact that all admissible stationary solutions are exact solutions for such scheme. However, one can use any other numerical flux that is monotone, consistent and Lipschitz continuous.

A convergence result for the scheme (13)-(14)-(15) can be found in [3].

5. Validation of the numerical scheme

We propose here to validate the numerical scheme (13)-(14)-(15) making a comparison with the explicit solution computed in the Section 3; by comparison with the solution of Riemann problems and by showing the consistence of the scheme with respect to a network consisting of one single arc.

We consider the explicit solution to (1) constructed in Section 3. The setup for the simulation is as follows. We consider $[-3/5,0]$ as domain of computation for the incoming arcs and $[0,3/5]$ for the outgoing one, and $\Delta x=0.5\times 10^{-4}$, $\Delta t=0.25\times 10^{-4}$ as space and time step, respectively. A qualitative comparison between the numeric solution $x\mapsto\rho_{\Delta}(t,x)$ and the explicit solution $x\mapsto\rho(t,x)$ at different fixed times $t$ is shown in Figure 2.

Additionally, we perform a convergence analysis for this test. We introduce the relative $\mathbf{L}^{1}$-error respectively for the whole network, for the incoming and for the outgoing arcs at a given time $t^{s}$ as follows

[TABLE]

Table 1 depicts the relative $\mathbf{L}^{1}$-error with respect to the space step at the fixed time $t=2.4$. The time step is fixed to $\Delta t=0.25\times 10^{-4}$. Since we are dealing with a first order scheme approximating discontinuous solutions, the sub-linear convergence rate found results expected.

5.1. Riemann problem for a $2$-$1$ merge.

We consider a network consisting of three edges and one junction, with two incoming and one outgoing arcs. We consider $[-\frac{1}{2},0]$ as domain of computation for the incoming arcs and $[0,\frac{1}{2}]$ for the outgoing one, and we take a normalized flux $f(\rho)=\rho\left(1-\rho\right)$ for each arc.

In Figures 3 we present a qualitative comparison between the numerically computed solution and the explicitly one at time $t=\frac{1}{2}$, corresponding to the Riemann problems having $\vec{\rho}_{0,a}=\left(1/4,1/3,4/5\right)$ and $\vec{\rho}_{0,b}=\left(1/4,2/3,1/5\right)$ as initial conditions, respectively.

We observe good agreement between the exact and the numeric solution. The parameters for the simulations are $\Delta x=10^{-4}$ and $\Delta t=0.5\,\Delta x=0.5\times 10^{-4}$.

5.2. Riemann problem for a $1$-$2$ divide.

We consider a network consisting of three edges and one junction, with one incoming and two outgoing arcs. We consider $[-\frac{1}{2},0]$ as domain of computation for the incoming arc and $[0,\frac{1}{2}]$ for the outgoing ones, and we take a normalized flux $f(\rho)=\rho\left(1-\rho\right)$ for each arc. We consider as initial conditions for the Riemann problems $\vec{\rho}_{0,a}=\left(1/4,2/3,4/5\right)$ and $\vec{\rho}_{0,b}=\left(3/4,1/3,4/5\right)$, respectively.

Figure 4 shows a qualitative comparison between the numerically computed solution and the explicitly one at time $t=\frac{1}{2}$, corresponding to the above initial conditions. Also in this case, we can observe good agreements between the exact solution and its numerical approximation. The parameters for the computed solution are $\Delta x=10^{-4}$ and $\Delta t=0.5\,\Delta x=0.5\times 10^{-4}$.

5.3. Riemann problem for $2$-$2$ network.

We consider here a junction with two incoming and two outgoing arcs. We consider $[-1/6,0]$ as domain of computation for the incoming arcs and $[0,1/6]$ for the outgoing ones, and we take a normalized flux $f(\rho)=\rho\left(1-\rho\right)$ for each arc. As initial condition for the Riemann problems on the network we choose $\vec{\rho}_{0,a}=\left(1/4,1/5,2/3,5/6\right)$ and $\vec{\rho}_{0,b}=\left(3/4,1/5,1/3,1/6\right)$.

Also in this case, in Figures 5 and 6 we find a good agreement between the exact and the numerical solutions. The parameters for the simulation are $\Delta x=10^{-4}$ and $\Delta t=0.5\,\Delta x=0.5\times 10^{-4}$.

5.4. A network with no discontinuity at junction.

In this section we consider a network consisting of an arc without any discontinuity at the junction, namely, we assume that the flux on the incoming arc coincides with the flux on the outgoing arc, and therefore, this setting is equivalent to the case of a single arc. This simulation exploits the idea consisting in solving two scalar conservation laws on half-space coupled by an ad hoc transmission condition at the interface [1]. We remark that, in the case of discontinuous flux, the transmission condition can be interpreted in terms of a flux constraint at the interface [2].

We apply the scheme to two different domains, one including the junction located at $x=0$, and one not. More in details, we choose $[-1/2,1/2]$ and $[0,1]$ as domain of computation, and $f(\rho)=\rho(1-\rho)$ as flux along the arcs. We consider the initial densities $\rho_{1}^{0}=0.75\chi_{[-1/4,0]}$ and $\rho_{2}^{0}=0.75\chi_{[0,1/4]}$ for the first simulation, and $\rho_{2}^{0}=0.75\chi_{[0,1/4]}$ for the second one. The parameters of computation are $\Delta x=10^{-4}$ and $\Delta t=0.5\,\Delta x=0.5\times 10^{-4}$.

Figure 7 displays the comparison between the profiles of solutions in the two simulations at times $t=2/5$ and $t=3/4$. We can observe the same qualitative behavior shifted of $|x|=1/2$.