Conservation laws with discontinuous flux function on networks: a splitting algorithm

1.   Introduction

Traffic flow models based on scalar conservation laws with continuous flux functions are widely used in the literature. For a general presentation of the models, we refer to the books [11,12,23] and the references therein. Extensions to road traffic networks have been also established. We mention in particular the contributions [6,15], where the authors introduce the coupled network problem and show the existence of solutions. Within this article, we are concerned with the special case of scalar conservation laws with discontinuous flux in the unknown that are motivated in the traffic flow theory by the observation of a gap between the free flow and the congested flow regime [4,5,8]. This phenomenon generates an interesting dynamical behavior called zero waves, i.e., waves with infinite (negative) speed but zero wave strength, and has been investigated in recent years either from a theoretical or numerical point of view, see for instance [2,19,20,22,24] or more generally [1,3,7,13].

To the best of our knowledge, the study of scalar conservation laws with discontinuous flux functions on networks is still missing in the traffic flow literature. However, in the context of supply chains with discontinuous flux such considerations have been already done [10,14]. We remark that supply chain models differ essentially from traffic flow models due to simpler dynamics and different coupling conditions.

In this work, we aim to derive a traffic network model, where the dynamics on each road are governed by a scalar conservation law with discontinuous flux function in the unknown. For simplicity, we restrict to piecewise linear flux functions. Special emphasis is put on the coupling at junction points to ensure a unique admissible weak solution. In particular, we focus on dispersing junctions where the number of incoming roads does not exceed the number of outgoing roads and merging with two incoming and one outgoing road. The latter type of junction can be extended to the case of multiple incoming roads and a single outgoing one. In order to construct a suitable numerical scheme that is not based on regularization techniques we adapt the splitting algorithm originally introduced in [22]. Therein, the discontinuous flux is decomposed into a Lipschitz continuous flux and a Heaviside flux such that a two-point monotone flux scheme, e.g., Godunov, can be employed in an appropriate manner. This algorithm has been studied in [22] for the case of a single road only. However, in the network case, multiple roads with possibly disjunctive flux functions need to be considered at a junction point to ensure mass conservation. Hence, the key challenge is to determine the correct flux through the junction in an appropriate manner. Therefore, a detailed case distinction in accordance with the theoretical investigations is provided for the different types of junctions. The numerical results validate the proposed algorithm for some relevant network problems.

The paper is organized as follows: in Section 2 we discuss the basic model and Riemann problems which permit to derive an exact solution. We extend the modeling framework to networks in Section 3 and focus on the coupling conditions. In Section 4, we introduce how the splitting algorithm [22] can be extended to also deal with the different types of junctions. Finally, we present a suitable discretization and numerical simulations in Section 5.

2.   LWR model with discontinuous flux on a single edge

In this section, we briefly recall the case of the Lighthill-Whitham-Richards (LWR) model [18,21] on a single edge with a flux function having a single decreasing jump at $u^* \in [0, u_{{\rm{max}}}]$ that has been intensively studied in e.g., [19,22,24]. The shape of this flux function is inspired by empirical research which suggest that the fundamental diagram is of a reverse $\lambda$ shape and consist of a free flow area and a drop to a congested area, see [4,5,8].

Following [22], we consider the scalar conservation law Eq (2.1),

More precisely the flux function is defined as follows Eq (2.2),

We denote $f_1(u^*) = f(u^*-)$ and $f_2(u^*) = f(u^*+)$, where $f(u^*-) > f(u^*+)$. In order to assign a single numerical value to $f(u^*)$ we arbitrarily define $f(u^*) = f(u^*-)$. One example is shown in Figure 1. The jump has magnitude

As usual we require for the flux function $f : [0, u_{\rm{max}}] \rightarrow \mathbb{R}^+$ with $f(0) = f(u_{{\rm{max}}}) = 0$ and in addition, we assume $f_1$ to be linear increasing and $f_2$ linear decreasing.

As in [22], the multivalued version of $f$ is defined by Eq (2.4),

Finally, we have to discuss the imposed boundary conditions at $x = b$ in Eq (2.1). Due to the flux discontinuity we have a non-standard boundary condition $\mathcal{F}(t) \in \tilde f (s(t))$. If the boundary condition at the right is $s(t) \neq u^*$, then the additional boundary condition is redundant, i.e., $\mathcal{F}(t) = f(s(t))$. Otherwise, if $s(t) = u^*$, we have to choose between two possible flux values $\tilde f(u^*) \in \{f(u^*+), f(u^*-)\}$. Then, the function $\mathcal{F}$ assigns a flux value to $u(b,t) = s(t)$ according to [22]

The state of traffic ahead of $x = b$ is additional information that is determined independently from the other data of the problem.

Remark 2.1. [22,Remark 1.3] We note that for the boundary condition at the left end the state of traffic ahead of $x = a$ is known from the start values. Hence, no additional boundary conditions are necessary.

The following assumptions are important for the proof of existence and uniqueness of solutions.

Assumption 2.2. [22,Assumption 1.1] The initial data satisfies $u_0(x) \in [0, u_{\mathit{{\rm{max}}}}]$ for $x \in (a,b)$ and $u_0 \in \mathit{{\rm{BV}}}([a,b]).$ The boundary data satisfies $r(t),s(t) \in [0, u_{\mathit{{\rm{max}}}}]$ for $t \in [0,T]$, and $r,s \in \mathit{{\rm{BV}}}([0,T])$.

A weak solution is intended in the following sense:

Definition 2.3 [22,Definition 1.1] A function $u \in L^{\infty} (\Pi_T)$ is called a weak solution to the initial-boundary value problem (2.1) if there exists a function $v \in L^{\infty}(\Pi_T)$ satisfying

such that for each $\psi \in C_0^1((a,b) \times[0,T)),$

As usual, weak solutions do not lead to a unique solution and additional criteria are necessary to rule out physically incorrect solutions. In particular, the discontinuity of the flux prohibits from directly using the classical approaches. Note that in [22] an adapted version of Oleinik's entropy condition [9] is used to single out the correct solution, while in [24] the convex hull construction [17] is used to construct solutions to Riemann problems.

Here, we will concentrate on the convex hull construction. For completeness we will shortly recall the solutions to Riemann problems considered in [24] as they are essential in order to construct a Riemann solver at a junction.

2.1.   Riemann solutions

We consider a Riemann problem with initial data $(u_L,u_R)$ and its exact solution. For simplicity we will focus on a piecewise linear flux function as in [24] and in Figure 1. The exact solution is then derived with a regularization of the flux and the convex hull construction. As regularizations are not unique we choose a simpler regularization than the one in [24]. In fact, we use the regularization introduced in [19] for a piecewise quadratic flux function and apply it to the piecewise linear flux function. In the following we simply recall the results from [19] in case of a linear flux function for the sake of completeness and in order to understand better the upcoming analysis. Note that the the results also coincide with the ones of [24].

We consider the following flux function by Eq (2, 6),

with the regularized flux function given by Eq (2.7),

We define $u_\epsilon : = u^* + \epsilon$. An example of the regularization is shown in Figure 2. Even though the regularization differs slightly from the one used in [24], the different cases to consider and the approach to solve the Riemann problems are completely analogous. We refer to [24,Section 3] for further details. We define $h: [0, u_{\rm{max}}] \rightarrow \mathbb{R}$ to be a function connecting the points $(u^*, f(u^*+))$ and the initial value lying on Flux 2. In addition, $\gamma$ denotes the intersection of $h$ and Flux 1, see also Figure 3. The point $\gamma$ distinguishes the cases 3 and 4 in the following discussion.

Case 1: Either $u_L < u^*, u_R < u^*$ or $u_L > u^*, u_R > u^*$

This case corresponds to the classical case of solving Riemann problems, where the solution consists of a single rarefaction wave or shock, see [17].

Case 2: $u_R < u^* < u_L$

By using the smallest convex hull approach the solution consists of a contact line following $f(u)$, connecting $u_R$ and $u^*$ and a shock, connecting $u^*$ and $u_L$. The speed of the contact discontinuity is given by $d_1>0$ and the speed of the shock can be calculated with the Rankine-Hugoniot condition, i.e.,

where we recall that $f(u^*) = f(u^*-)$ holds. The exact solution is then given by Eq (2.8),

Case 3: $\gamma < u_L < u^* < u_R$

Here, the solution is given by a shock connecting $u_L$ and $u^*$ and contact discontinuity reconnecting it to $u_R$. The speed of the shock equals:

Note that due to $u_L<u_R$ the convex hull construction chooses the flux $f(u^*+)$ at $u^*$. The contact discontinuity connecting $u^*$ and $u_R$ moves at the speed $e_1<0$. Hence,

Case 4: $u_L < u^* < u_R$ and $\gamma \geq u_L$

In this case, we get only one shock connecting $u_L$ and $u_R$ as a solution due to the condition $\gamma \geq u_L$. The speed of the shock can be calculated by the Rankin-Hugoniot condition. The exact solution is given by Eq (2.10),

Remark 2.4. As aforementioned in [19] Riemann solutions for piecewise quadratic discontinuous flux functions are derived. They also cover the case of a piecewise linear flux function if the quadratic terms are zero. For a general quadratic discontinuous flux, the solutions are more involved since no contact discontinuities occur.

Up to now, we have not addressed the case, where one of the boundary conditions equals the critical density $u^*$. In this situation so called zero waves can occur, see [24]. Those waves travel with speed $\mathcal{O}(\frac{1}{\epsilon})$ from right to left and have strength $\mathcal{O}(\epsilon)$, depending on the regularization parameter $\epsilon$. Although they have no strength, information can be exchanged between Riemann problems lying next to each other influencing the solution. For more details on the origin of zero waves in the context of traffic flow we refer to [24]. In particular, they arise when double Riemann problems are considered. For a detailed discussion of these double Riemann problems we refer to [24,Section 4].

3.   LWR model with discontinuous flux on networks

Next, we focus on networks where we allow for discontinuous flux functions. The key idea is to consider the regularized flux function $f^\epsilon$ and take the limit $\epsilon \rightarrow 0$ to obtain results for the original problem at intersections. In doing so, the standard approaches can be directly applied to the regularized flux function $f^\epsilon$.

We start with a short introduction to the network setting. For more details on traffic flow network models we refer the reader to [11] and the references therein.

Let $G = (\mathcal{V}, \mathcal{E})$ be a directed finite graph consisting of a nonempty set $\mathcal{V}$ of vertices and a nonempty set of edges $\mathcal{E}$ representing the roads. We denote by $\delta_v^-$ and $\delta_v^+$ the set of all incoming and outgoing edges for every vertex $v \in \mathcal{V}$. Every edge $e \in \mathcal{E}$ is interpreted as an interval $I_e = (a_e,b_e) \subset \mathbb{R}$ representing the spatial extension of the road $e$. For all edges $e$, which do not discharge into a vertex, i.e. $e \notin \cup_{v \in \mathcal{V}} \delta_v^-$, we set an infinite road length by $b_e = +\infty$. Analogously, for all edges $e$, which do not originate from a vertex, i.e. $e \notin \cup_{v \in \mathcal{V}} \delta_v^+$, we set an infinite road length by $a_e = -\infty$.

We call the couple $(\mathcal{I}, \mathcal{V})$ with $\mathcal{I} = \{I_e: e\in \mathcal{E}\}$ a road network.

In order to derive the network solution, we restrict to the description of a single junction $v\in\mathcal{V}$. Hence, from now on we consider a fixed junction with $n$ incoming and $m$ outgoing roads. Then, we need additional information how vehicles are distributed from one road to another. Hence, we need to define a distribution matrix $A \in [0,1]^{n \times m}$ as in [11]:

To conserve the mass we assume $\sum_{j = 1}^m \beta_{i,j} = 1$ for every $i \in \{1, \dots, n\}$. Furthermore, some technical assumptions on the distribution matrix are necessary which can be found in ref. [11] (Eq (4.2.4)). These assumptions are necessary to obtain unique solutions at the junction for the case $n\leq m$. Following ref. [11] (definition 4.2.4) and definition 2.3, weak solutions at a junction are intended in the following sense

Definition 3.1 (Weak solution at a junction). Let $v \in \mathcal{V}$ be a vertex with $n$ incoming and $m$ outgoing roads, represented by the intervals $I_1^{{\rm{in}}}, \dots I_n^{{\rm{in}}}$ and $I_1^{{\rm{out}}}, \dots I_n^{{\rm{out}}}$. A weak solution at the junction $v$ is a collection $u$ and $w$ of functions $u_i^{{\rm{in}}}: I_i^{{\rm{in}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}$, $w_i^{{\rm{in}}}: I_i^{{\rm{in}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}$ with $w_i^{{\rm{in}}}(x,t)\in \tilde f(u_i^{{\rm{in}}}(x,t))$ a.e., where $i = 1, \dots,n$, $u_j^{{\rm{out}}}$ : $I_j^{{\rm{out}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}$ and $w_j^{{\rm{out}}}: I_j^{{\rm{out}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}$ with $w_j^{{\rm{out}}}(x,t)\in \tilde f(u_j^{{\rm{out}}}(x,t))$ a.e., where $j = 1, \dots,m$, respectively, such that

for every collection of test functions $\phi_i^{{\rm{in}}} \in C_0^1((a_i^{{\rm{in}}},b_i^{{\rm{in}}}] \times [0,T]), i = 1, \dots,n$ and $\phi_j^{{\rm{out}}} \in C_0^1((a_j^{{\rm{out}}},b_j^{{\rm{out}}}] \times [0,T]), j = 1, \dots,m$, satisfying

for $i \in {1,\dots n} \rm{ and } j \in {1, \dots,m}.$

Additionally, in order to get unique solutions, we will consider the following concept of admissible solutions, which adapts ref. [11] (rule (A) and (B), p. 81) to the discontinuous setting:

Definition 3.2 (Admissbile Weak Solution). We call $u$ and $w$ (as defined in definition 3.1) an admissible weak solution at a junction $v \in \mathcal{V}$ with the corresponding distribution matrix $A$ if

1. $u_i^{{\rm{in}}}(\cdot,t) \in {\rm{BV}}(I_i^{{\rm{in}}}), u_j^{{\rm{out}}}(\cdot,t) \in {\rm{BV}}(I_j^{{\rm{out}}}) \;{\rm{ for\; all }}\; i,j,$

2. $u$ and $w$ constitute a weak solution at the junction $v$,

3. ${w_j^{{\rm{out}}}}(a_j^{{\rm{out}}} +,t)) = \sum_{i = 1}^n \beta_{i,j} {w_i^{{\rm{in}}}}(b_i^{{\rm{in}}} -,t)) \;{\rm{ for\; all }}\; j = 1, \dots, m \rm{ and } t \geq 0,$

4. $\sum_{i = 1}^n {w_i^{{\rm{in}}}}(b_i^{{\rm{in}}} -,t))$ is maximal subject to 2. and 3.

In particular, the maximization of the inflow with respect to the distributions parameters and the technical assumption [11] guarantee the uniqueness of solutions for a continuous flux.

If $n > m$, that means more incoming than outgoing roads, we need so-called right of way parameters which prescribe in which order the vehicles are going to drive. We consider the special case $n = 2$ and $m = 1$. Let $C \in { \mathbb{R}}$ denote the amount of vehicles that can move on the road. The right of way parameter $q \in (0,1)$ describes the following situation: If not all vehicles can move to the outgoing road, then $qC$ is the amount of vehicles that enter from the first and $(1-q)C$ the amount from the second road. If there are more roads, we need to define several right of way parameters.

For solving the maximization problems imposed by the definition 3.2 so-called supply and demand functions can be used, see [16]. The demand describes the maximal flux the incoming road wants to send. In contrast, the supply describes the maximal flux the outgoing road is able to absorb. The definition of the supply and demand functions of the regularized function is straightforward. As $f(u^*) \in \{ f(u^*+), f(u^*-)\}$ is not unique, we have to consider the limit for $\epsilon \rightarrow 0$ to determine the correct values of the supply and demand for the discontinuous flux function. Nevertheless, it becomes apparent that the maximal flux in both cases is given by $f(u^*-)$. This leads to the following definition:

Definition 3.3. For a network with flux function $f$ having a discontinuity at $u^*$, the demand is given by

On the contrary, the supply reads as

and

Remark 3.4. If we consider the regularized flux function $f^\epsilon$, the definition of supply and demand function is completely analogous by replacing $f$ with $f^\epsilon$.

In order to show existence and uniqueness in the discontinuous case we need to define an additional function. For a regularized flux function we notice that for every flux value, we get two different density values, see left picture in Figure 4. As different density values lead to different solutions, we need to be able to distinguish them and choose the correct solution. In the continuous or regularized case a mapping usually called $\tau$ is used to determine the density values for which the flux values coincide, see also left picture in Figure 4 and ref. [11] (definition 3.2.6) for the precise definition of $\tau$. This situation becomes more complicated for the discontinuous flux functions. For $\epsilon \rightarrow 0$ and $f(u) > f(u^*+)$, the corresponding density value given by $\tau$ can converge to the discontinuity area. This situation is pictured in Figure 4 on the right. Therefore, we need to adapt the definition of the mapping $\tau$ to also cover this case and determine the corresponding density value. This mapping is defined in the following and is denoted by $\eta$. For readability, we also define as an abbreviation $f^{-1}_+: = f^{-1}(f(u^*+))$.

Definition 3.5. Let the function $\eta : [0, u_{\rm{max}}] \rightarrow [0, u_{\rm{max}}]$ satisfy:

1. $f(\eta(u)) = f(u)$ for every $u \in \bigl([0, { f^{-1}_+}] \cup [u^*,u_{\rm{max}}]\bigr)$ and $\eta(u) \neq u$ for every $u \in \bigl([0, { f^{-1}_+}] \cup (u^*,u_{\rm{max}}]\bigr)$.

2. For $u \in \bigl({ f^{-1}_+}, u^*\bigr)$ it holds, $\eta(u) = u^*$ and $f(\eta(u)) = f(u)$.

Note that if $u \in \bigl({ f^{-1}_+}, u^*\bigr)$, the flux value is given by $f(u)\in(f(u^*+),f(u^*-))$.

Remark 3.6. We note that the mapping $\tau$ from ref. [11] (definition 3.2.6) is in principle given by the first statement in definition 3.5 with the separated interval $[0, { f^{-1}_+}] \cup [u^*,u_{\rm{max}}]$ being replaced by $[0,u_{\rm{max}}]$. Hence, the only difference between $\tau$ and $\eta$ lies in the treatment of $u \in \bigl({ f^{-1}_+}, u^*\bigr)$.

3.1.   Riemann solutions at a junction

Now, we present a Riemann solver for two types of junctions. First, we consider a junction with $n$ incoming and $m$ outgoing roads with $n\leq m$, second, a junction with $n = 2$ incoming and $m = 1$ outgoing roads. For the continuous case the existence of solutions satisfying in particular the definition 3.2 has been shown in [11]. Since the continuous flux is not multivalued, the definition in 3.2 coincides with ref. [11] (Definition 4.2.4). As we have seen in Section 2.1, the solution to a Riemann problem may produce more than one wave. In order to have a well-defined network problem it is important that the waves produced at the junction have negative speed on the incoming and positive speed on the outgoing roads. This is ensured by the following theorem:

Theorem 3.7. Let $w \in \mathcal{V}$ be a junction with $n$ incoming and $m$ outgoing roads with $n\leq m$. For every constant initial datum $u_{i,0}^{\mathit{{\rm{in}}}} , u_{j,0}^{\mathit{{\rm{out}}}} \in [0, u_{\mathit{{\rm{max}}}}],\ i = 1,\ldots,n,\ j = 1,\ldots,m$, there exists a unique admissible weak solution $u$ in the sense of definition 3.2. Moreover, for the incoming road the solution is given by the wave $(u_{i,0}^{\mathit{{\rm{in}}}} , u_{i}^{\mathit{{\rm{in}}}}),\ i = 1,\ldots,n$ and for every outgoing road the solution is the wave $(u_{j,0}^{\mathit{{\rm{out}}}},u_{j}^{\mathit{{\rm{out}}}}), j = 1,\ldots,m$, where

and

Proof. Using the definition of the supply and demand functions in definition 3.3 and the results from [12,section 5.2.3] we can follow the proof of [12,theorem 5.1.2] and uniquely determine the inflows which maximize the flux through junctions subject to the distribution parameters. It remains to show that by the choice of the density values the correct waves are induced.

We start with considering the outgoing roads. If $u_{j,0}^{{\rm{out}}}\in [0,u^*)$, the choice $u_j^{{\rm{out}}}\in[0,u^*]$ produces the classical case of a shock or contact discontinuity, both with the positive speed $d_1$, since the flux is a linear increasing function. If $u_{j,0}^{{\rm{out}}} = u^*$, we need to know if the traffic ahead is congested or not in order to determine the correct flux value of $u^*$ and the corresponding density. In particular, we are now concerned with a double Riemann problem by $u_j^{{\rm{out}}}$, $u_{j,0}^{{\rm{out}}} = u^*$ and the information about the traffic ahead. By the analysis done in [24,Section 4] the choices of $u_j^{{\rm{out}}}$ give rise to waves with positive velocities, even though the velocities differ depending on the traffic situation. Note that the free flowing case corresponds to case 1 in [24,Section 4] and the congested case to case 2. In particular, in the latter case $f_j^{{\rm{out}}}\leq f(u^*+) = f(u^*)$ holds such that the solution is a shock with velocity always greater or equal to zero, compare also case 4 of Section 2.1.

For $u_{j,0}^{{\rm{out}}}\in (u^*,u_{{\rm{max}}}]$ we need to choose the density $u_{j}^{{\rm{out}}} \in \{ u_{j,0}^{{\rm{out}}} \} \cup [0, \eta(u_{j,0}^{{\rm{out}}}))$ which corresponds to the incoming flux. Here, we are in case 4 of Section 2.1, as $u_{j,0}^{{\rm{out}}}<\gamma$, such that the solution is a shock with positive velocity.

On the contrary, considering the incoming roads and $u_{i,0}^{\rm{in}} \in [0, { f^{-1}_+}]$, the choice $u_i^{{\rm{in}}}\in \{ u_{i,0}^{{\rm{in}}} \} \cup (\eta(u_{i,0}^{{\rm{in}}}),u_{{\rm{max}}}]$ also corresponds to case 4, but this time the shock has a negative velocity.

Now, let $u_{i,0}^{\rm{in}} \in (u^*, u_{\rm{max}}]$ and $f_i^{{\rm{in}}}$ denote the uniquely determined inflow. If $f_i^{{\rm{in}}}<f(u^*+)$, the choice $u_i^{{\rm{in}}}\in (u^*, u_{\rm{max}}]$ corresponds to the solution of a Riemann problem with a linear decreasing function and hence either a shock or contact discontinuity with negative velocity. On the contrary if $f_i^{{\rm{in}}}\geq f(u^*+)$, we choose $u_i^{{\rm{in}}} = u^*$. Here, the convex hull reconstruction induces a shock as a solution with velocity

Further, if $u_{i,0}^{\rm{in}} = u^*$, we would normally also need to consider a double Riemann problem for which we would need the information of boundary data of the incoming roads. Nevertheless, the solutions at the junction are the same. As can be seen in [24,section 4] cases 2 and 4, if $f_i^{{\rm{in}}}<f(u^*+)$ the choice $u_i^{{\rm{in}}}\in (u^*, u_{\rm{max}}]$ induces a contact discontinuity with negative velocity as the solution at the junction. If $f_i^{{\rm{in}}}\geq f(u^*+)$ we choose $u_i^{{\rm{in}}} = u^*$ which gives a constant solution.

Now, let us turn to the remaining case of $u_{i,0}^{\rm{in}} \in ({ f^{-1}_+}, u^*)$. Again the solution depends on the inflow. If $f_i^{{\rm{in}}}<f(u^*+)$, the choice of $u_{i}^{{\rm{in}}}\in\{u_{i,0}^{{\rm{in}}}\} \cup (u^*,u_{{\rm{max}}}]$ induces negative waves. Nevertheless, depending on $u_{i,0}^{{\rm{in}}}> \gamma$ or $u_{i,0}^{{\rm{in}}}\leq \gamma$ we are either in case 3 or case 4 of Section 2.1. Even though case 3 induces two waves, all have a negative velocity. If $f_i^{{\rm{in}}}\geq f(u^*+)$, choosing $u_i^{{\rm{in}}} = u^*$ gives a shock with velocity

as $f_i^{{\rm{in}}}$ is bounded from above by $f(u_{i,0}^{{\rm{in}}})$ due to the definition of the demand.

Hence, the choices of the densities induce the correct waves.

Now, we consider the case of more incoming than outgoing roads. Exemplary, we study the 2–to–1 situation, even though the results can be easily extended to the $n$–to–1 case. We assume that a right of way parameter $q$ is given and the unique density of the discontinuous flux can be determined according to the following theorem:

Theorem 3.8. Let $w \in \mathcal{V}$ be a junction with $n = 2$ incoming and $m = 1$ outgoing roads with a right way parameter $q \in (0, 1)$. For every constant initial datum $u_{i,0}^{\mathit{{\rm{in}}}} , u_{1,0}^{\mathit{{\rm{out}}}} \in [0, u_{\mathit{{\rm{max}}}}]$ with $i \in \{1,2\}$, there exists a unique admissible weak solution $u$ in the sense of definition 3.2. Moreover, for every incoming road $i$ the solution is given by the wave $(u_{i,0}^{\mathit{{\rm{in}}}} , u_{i}^{\mathit{{\rm{in}}}})$ and for the outgoing edge the solution is the wave $(u_{1,0}^{\mathit{{\rm{out}}}},u_{1}^{\mathit{{\rm{out}}}})$, where

and

Proof. Following [11,Section 3.2.2], the flux values at the junction can be calculated with the following steps:

1. Calculate the maximal possible flux $f_{\rm{max}} = \min\{ D(u_{1,0}^{\rm{in}}) + D(u_{2,0}^{\rm{in}}), S(u_{1,0}^{\rm{out}})\}$.

2. Consider the right of way parameter and the flux maximization and calculate the intersection $P = (qf_{\rm{max}}, (1-q)f_{\rm{max}})$.

3. If $P$ is in the feasible area $\Omega = \{ z\in [0, D(u_{1,0}^{in})] \times [0, D(u_{2,0}^{\rm{in}})]: z_1 + z_2 \leq f_{\rm{max}}\}$, set $(z_1,z_2) = P$. If not, determine the point in $\Omega \cap \{(z_1,z_2): z_1+z_2 = f_{{\rm{max}}}\}$ which is the closest to the line $z_2 = (1-q)/q z_1$.

4. The flux values are given by $f^{\rm{in}}_1 = z_1,\ f^{\rm{in}}_2 = z_2,\ f^{\rm{out}}_1 = z_1+z_2$.

Completely analogous to theorem 3.7 we can show that the choice of the densities admits the correct wave speeds.

The Riemann solutions proposed in theorem 3.7 and theorem 3.8 are the key ingredients for the splitting algorithm on networks in the next section.

4.   A splitting algorithm on networks

Different problems might occur when designing a numerical scheme for a conservation law with discontinuous flux. However, the main difficulties are induced by the zero waves. Since these waves have infinite speed, the regular CFL condition is scaled by the regularization parameter $\epsilon$. This leads to very inefficient step sizes and high computational effort. Possible ways to avoid using the regular CFL condition are implicit methods [20] or the splitting algorithm introduced in [22], see also [2]. In this section we recall the main ideas of the splitting algorithm and also present its extension to networks.

We consider a flux function $f$ with discontinuity at $u^*$, which satisfies assumption 2.2. The jump has magnitude $\alpha = f(u^*-)-f(u^*+)$. The main idea of the splitting algorithm is to split the discontinuous function $f$ into two parts, namely $p$ and $g$, such that $f(u) = p(u) + g(u)$. Here, $g$ is a piecewise constant function defined by

where $H$ denotes the Heaviside function. This function extracts the jump from the function $f$. In addition, we define $p(u) = f(u) - g(u).$ Hence, $p$ is continuous and can be used for numerical algorithms. An example can be seen in Figure 5.

4.1.   A single road

This case has been already treated in [22] and will be the basis for the splitting algorithm on networks. When solving the scalar conservation law (2.1) on a single road, the boundary value in the case $s(t) = u^*$ is not unique. In the continuous case, the flux value is determined by $\mathcal{F}(t)$ in Eq (2.5). When splitting the function $f$, the discontinuity is shifted to the function $g$. So analogously, we define the multivalued version of $g$ in order to be able to assign a flux value in the critical case, i.e.,

Furthermore, we define $\mathcal{G}(t) \in \tilde g(s(t)).$ It holds,

Additionally to the assumptions 2.2, we assume:

Assumption 4.1. [22,Assumption 1.1] The initial data satisfies $u_0(x) \in [0, u_{\mathit{{\rm{max}}}}]$ for $x \in (a,b),\ u_0 \in \mathit{{\rm{BV}}}([a,b])$ and $g(u_0) \in \mathit{{\rm{BV}}}([a,b]).$ We also assume that $\mathcal{G}(t) \in [-\alpha,0]$ for $t \in [0,T]$, and $\mathcal{G} \in \mathit{{\rm{BV}}}([0,T])$.

Remark 4.2. We emphasize that the original splitting algorithm for a single road [22] is not limited to piecewise linear discontinuous flux functions. Another prominent example might be concave piecewise quadratic flux functions with discontinuity again at $u^*$, cf. [19,22,24].

Then, we are able to handle the flux function $p$ and the jump $g$ separately. As usual, we discretize the domain $\Pi_T$ by a spatial mesh and a temporal mesh. Let $\Delta x$ denote the spatial and $\Delta t$ the temporal step size. We define the grid constant via

For an integer ${K} \in \mathbb{N}$ the grid size is given by $\Delta x = \frac{b-a}{{K}+1}$ and the grid points by $x_{{k}} = a+{k} \Delta x$. Further, we define the ${K}$ disjunct intervals $I_{{k}} = [x_{{k}} - \frac{\Delta x}{2}, x_{{k}} + \frac{\Delta x}{2})$ and $\mathcal{{K}} = \{ 1, \cdots {K}\}$ and $\mathcal{{K}}^+ = \{ 0, \cdots {K}+1\}$. The time steps are defined by $t^n = n \Delta t$, for $n = 0,1, \dots,N$. The value $N \in \mathbb{N}$ is chosen such that the time horizon fulfills $T \in [t^N, t^N+\Delta t).$

As the algorithm splits the function $f$ into two parts, $f(u) = p(u) + g(u)$, every step of the algorithm consists of two half steps. Here, the first half step corresponds to approximately solving $u_t + g(u)_x = 0$ denoted by $U_{{k}}^{n+\frac{1}{2}} \approx u(x_{{k}},t^n).$ Consequently, the second step approximates $u_t + p(u)_x = 0$ and is denoted by $U_{{k}}^n \approx u(x_{{k}},t^n).$

We denote the backward spatial difference by $\Delta_- U_{{k}}^n = U_{{k}}^n - U_{{k}-1}^n$ and the forward spatial difference by $\Delta_+ U_{{k}}^n = U_{{k}+1}^n - U_{{k}}^n$. The initial values are denoted by $U_{{k}}^0 = u_0(x_{{k}})$, the boundary values by

The function $g$ covers the issue of $s(t^n) = u^*$. We define

That means, we can describe the boundary value $\mathcal{F}(t)$ via $\mathcal{G}(t)$. For the splitting algorithm in its original fashion [22] we additionally need the following function.

Definition 4.3 [22,Eq (3.7)] Let $G(z) : = z - \lambda g(z)$ with grid constant $\lambda$ and $\tilde G$ its multivalued version. The function $\tilde G$ is strictly increasing and has a unique inverse, $z \rightarrow \tilde G^{-1}(z)$, which is a single valued function. It is Lipschitz continuous and nondecreasing. It holds, $0 \leq \tilde G^{-1}(z) \leq z$. The functions are given by

The splitting algorithm [22] can be then expressed as

Note that the first half step, which includes the first two equations, is implicit. Nevertheless, instead of solving a nonlinear system of equations, the equation can be solved backwards in space starting with ${k} = {K}$. The last equation in (4.4) uses the Godunov scheme for the second half step which flux function is denoted by $p^g$.

We note that for the implicit equation a CFL condition is not needed, but it is required for the third step. As $p$ is continuous, we can use the standard CFL condition of the Godunov scheme. Note that similar to [22] also other two-point monotone schemes, e.g., the Lax-Friedrichs method, can be used.

As shown in [22,Theorem 5.1] the splitting algorithm (4.4) converges to a weak solution of Eq (2.1). However, obtaining a similar statement about weak entropy solutions is still an open problem.

4.2.   Networks

The key idea to numerically solve such discontinuous conservation laws on networks is to use the splitting algorithm only on the roads and determine the correct in- and outflows at the boundaries by the help of the Riemann solver established in, e.g., theorem 3.7. As the splitting algorithm works with flux values, there is no need to compute the exact densities at the junction. Instead we need to know how the solution at the junction influences the flux values. The algorithm for a single junction is depicted in algorithm 1. The general description of the algorithm allows for either junctions with given distribution or right of way parameters. For simplicity, we assume that each road is represented by the same interval $(a,b)$, such that we have $I_1^{{\rm{in}}} = \dots = I_n^{{\rm{in}}} = I_1^{\rm{out}} = \dots = I_m^{\rm{out}} = (a,b)$.

Remark 4.4. Note that this simplification enables the use of the same grid points on each road which spares further sub- or superindices. However, the algorithm can be easily adapted to different road lengths.

We assume in the following that the space and time grid is the same as in the previous subsection. The approximate solutions are denoted by $U_{i,k}^{n,{\rm{in}}}\approx u_{i,k}^{n,{\rm{in}}}(x_k,t^n)$ for $i = 1,\dots,n$ and $U_{j,k}^{n,{\rm{out}}}\approx u_{j,k}^{n,{\rm{out}}}(x_{{k}},t^n)$ for $j = 1,\dots,m$. The half steps are denoted accordingly. We assume that the initial values $u^{\rm{in}}_{i,0}(x)$ and $u^{\rm{out}}_{j,0}(x)$ and the boundary values $u^{\rm{in}}_{i,0}(a,t)$ and $u^{\rm{out}}_{j,0}(b,t)$ are given. Then, we can directly set the starting values as $U_{i,k}^{0,{\rm{in}}} = u^{\rm{in}}_{i,0}(x_k)$ for $i = 1,\dots,n$ and $U_{j,k}^{0,{\rm{out}}} = u^{\rm{out}}_{j,0}(x_k)$ for $j = 1,\dots,m$. The left boundary values for the incoming roads are given by $U_{i,0}^{n,{\rm{in}}} = u^{\rm{in}}_{i,0}(a,t^n)$ and the right boundary values for the outgoing roads by $U_{j,0}^{n,{\rm{out}}} = u^{\rm{out}}_{j,0}(b,t^n)$. The missing boundary data, i.e. the right one for the incoming road and the left one for the outgoing road, are determined by the flux values at the junction.

The overall strategy of the splitting algorithm on a network consists of three important steps:

1. Solve the optimization problem induced by definition 3.2 (in particular item 4) at the junction to calculate the flux values $f_i^{\rm{in}}, f_j^{\rm{out}}$ with the discontinuous flux $f$.

Here, it is crucial to use the discontinuous flux function $f$ when calculating demand and supply at the junction, which are necessary for the optimization problem, since otherwise not only the flux values can be different but also the decision whether supply or demand is active. Therefore, in line 6, we need to compute first the fluxes at the junction with the original discontinuous flux function $f$.

These flux values bring us now to:

$2.$ Consider the corresponding density values at the junction applying either theorem 3.7 or theorem 3.8 depending on the type of junction. If the density value is greater or equal than $u^*$, i.e., the road is congested, the flux value needs to be adjusted. Hence, we determine adjusted flux values $f_\rm{adj.}^\rm{in/out}$.

Using the calculated (unadjusted) flux values from step one, we can determine the densities at the junction with the help of the appropriate Riemann solver (theorem 3.7 and Eqs (3.4)–(3.5) or theorem 3.8 and Eqs (3.6)–(3.7)) at the junction (line 7). Then, these density values can be used to calculate the corresponding flux value of $p$ to get the boundary values and decide about an adjustment. This means that if the density given by the Riemann solver suggests that the traffic ahead is free flowing, no adjustment needs to be made - as the flux value already lies on $p$. If the traffic ahead is congested we take the flux value on the curve $p$ for the corresponding density value and not its original flux value on the curve $f$. Both the first step, i.e solving the optimization problem at the junction, and the second step, i.e., how to adjust the flux, strongly depend on the junction type. Hence, we postpone a detailed discussion to the following subsections, in which we will consider the 1–to–1, 1–to–2 and 2–to–1 case in more detail.

In addition, the first steps give us all the ingredients for the final step:

$3.$ Away from the junction, the splitting algorithm can be used. At the junction itself, the adjusted fluxes are used to determine the missing boundary data.

The adjusted flux values from the previous step are important for the second half step (line 17 or 28) of the splitting algorithm which uses a Godunov type scheme based on $p$. Here, the missing fluxes in and out of the junction the Godunov scheme needs are simply given by these adjusted flux values. The flux values can be used directly as boundary conditions, see lines 15 and 26, and the exact density values are not needed for the scheme.

Further boundary data is needed in the first half step of the algorithm, lines 12–13 and 23–24. Here, we start with $k = K$ to avoid solving a nonlinear equation system. For example, the value $g^{n,{\rm{out}}}_{j,K+1}$ is determined by the density value at the boundary. Nevertheless, for every outgoing road the boundary values are known and the value $g^{n,{\rm{out}}}_{j,K+1}$ can be determined as before in Eq (4.3). However, for the incoming road, i.e., for $g^{n,{\rm{in}}}_{i,K+1}$, the density value at the junction is necessary. From the adjustment made in step two we know whether the traffic ahead is congested or not. In addition, instead of working with the exact density values at the junction directly, we can also use the adjusted flows computed in the second step (or line 8). Here, we can simply compute the boundary value $g^{n,{\rm{in}}}_{i,K+1}$ by the difference between the fluxes calculated with $f$ (before adjustment) and the ones with $p$ (after adjustment), i.e.,

Note that the definition of $g^{n,{\rm{in}}}_{i,K+1}$ in Eq (4.5) is slightly different in comparison to Eq (4.3) or for the outgoing roads. Nevertheless, if the density value on the incoming road is smaller than $u^*$, $g^{n,{\rm{in}}}_{i,K+1}$ equals zero (no adjustment needs to be made) and if it is greater $u^*$, it equals $-\alpha$. Only in the case the density value equals $u^*$ the situation is more involved, as the traffic on the incoming road can be congested and $f(u^*-)>f^{{\rm{in}}}>f(u^*+)$ can occur such that we need $g^{n,{\rm{in}}}_{i,K+1}\in (-\alpha,0)$.

Furthermore, we can decrease the computational costs of the algorithm: In the second step (or in line 7 of algorithm 1) the density at the junction is computed. This can be very expensive and hence we aim to avoid this. In the third step we have seen that for the missing boundary data only the adjusted flux values are necessary and not the densities at the junction themselves. Therefore, studying first each junction type in detail allows to determine the corresponding flux values based on the density values and the supply and demand functions and the intermediate expensive step for the computation of the exact densities can be skipped. Hence, we combine the first and second step of the strategy in one single step. In the following, we will study a 1–to–1 junction in detail and present the tailored algorithm. As the strategy is completely analogous for the 1–to–2 junction and 2–to–1 junction, we will only present the algorithms and discuss important properties. The algorithms can then be used to replace the lines 6 to 8 in algorithm 1.

Remark 4.5. The extension to $n$–to–1 or 1–to–$m$ junctions follows immediately from the upcoming discussion. For more complex junctions, it is possible to adapt the proposed strategy by studying the different solutions of the optimization problem at the junction.

Remark 4.6. Further note that the presented strategy and also algorithm 1 can be used for arbitrary junctions and nonlinear discontinuous flux functions once an appropriate Riemann solver similar to theorem 3.7 and 3.8 is established.

4.2.1.   1–to–1 junction

First, we consider a junction with one incoming and one outgoing road in detail. Let $u_{1,0}^{\rm{in}}, u_{1,0}^{\rm{out}}$ be constant initial values. The flux values are given by $f^{{\rm{in}}}_1 = f^{{\rm{out}}}_1 = \min\{D(u_{1,0}^{\rm{in}}),S(u_{1,0}^{\rm{out}})\}$. We distinguish the following cases:

Case A: demand and supply are equal

There are two different situations depending on the density value of the incoming road where demand and supply can be equal.

1. $\bf{u_{1,0}^{\rm{in}}} \bf{\geq u^*}$: Here, the traffic on the outgoing road needs to be free-flowing or equal to $u^*$ and the traffic ahead of the outgoing road is free-flowing. Supply and demand equal the maximal possible flux $f(u^*-)$. For these flux values the solution of the Riemann solver on the incoming and outgoing road is given by $u^*$. As corresponding flux values already lie on $p$, we do not need to adjust them.

2. $\bf{u_{1,0}^{\rm{in}}} \bf{< u^*}$: Demand and supply can also be equal, if the number of vehicles the demand wants to send is equal to the number of vehicles the supply can take. Therefore, $u_{1,0}^{\rm{in}} < \gamma$ and $u_{1,0}^{{\rm{out}}}>u^*$ need to hold. The solution proposed by the Riemann solver at the junction is the initial condition for each road. Hence, the traffic on the outgoing road stays congested and we need to adjust the flux value on the outgoing road, i.e., $f_{1,\rm{adj}}^{\rm{out}} = f_1^{\rm{out}} + \alpha$. The flux on the incoming road stays free-flowing.

Case B: supply is restrictive

If the supply is restrictive, i.e. $f^{{\rm{in}}}_1 = f^{{\rm{out}}}_1 = S(u_{1,0}^{\rm{out}})<D(u_{1,0}^{\rm{in}})$, the outgoing road cannot take the whole amount of cars, the demand wants to send. From this follows that the traffic on the outgoing road is congested, i.e. $u_{1,0}^{\rm{out}} > u^*$ or $u_{1,0}^{\rm{out}} = u^*$ and the traffic ahead being congested. This has an effect on the incoming road. Again, we face two situations:

1. $\bf{u_{1,0}}^{\rm{in}} \bf{< u^*}$: Initially, the traffic on the incoming road is free-flowing and the traffic on the outgoing road is congested. As the supply is restrictive, the traffic on the incoming road congests as well, the solution of the Riemann solver is on both roads given by $u_{1,0}^{\rm{out}}$. Therefore, we have to adjust the flux value, $f_{1,\rm{adj}}^{\rm{in}} = f_1^{\rm{in}}+\alpha.$ In addition, we need to set ${g_{1,K+1}^{n,{\rm{in}}}} = -\alpha$. The outgoing flux needs to be also adjusted, i.e., $f_{1,\rm{adj}}^{\rm{out}} = f_1^{\rm{out}}+\alpha.$

2. $\bf{u_{1,0}}^{\rm{in}} \bf{> u^*}$: If the traffic on the incoming road is congested as well, the solution of the Riemann solver on both roads is again given by $u_{1,0}^{\rm{out}}$. This leads to the same adjustments as in the previous case.

Case C: demand is restrictive

If the demand is restrictive, the outgoing road is able to take the whole amount of vehicles the demand wants to send. Here, we only have one possible situation for the initial condition on the incoming road:

1. $\bf{u_{1,0}}^{\rm{in}} \bf{\leq u^*}$: The solution of the Riemann problem on each road is now given by $u_{1,0}^{\rm{in}}$, such that no adjustment is necessary.

Note that we have the same flux function on each road. Hence, in the case ${u_{1,0}}^{\rm{in}} {> u^*}$ either the supply is restrictive or demand and supply are equal.

The whole procedure is summarized in algorithm 2.

Remark 4.7. Theoretically, the Riemann solver in theorem 3.7 coincides for a 1–to–1 junction with the Riemann solver on a single road. Hence, the procedure described in algorithm 2 leads to the same solution. In contrast to that on the numerical level, the splitting algorithm for a 1–to–1 junction does not exactly coincide with the splitting algorithm used for a single road [22]. The reason for the computational difference is that the splitting algorithm for the 1–to–1 junction considers the exact solution of the Riemann problem at the junction point and hence uses exact values for $g_{1,K+1}^{n,{\rm{in}}}$ while for a single road this value is only approximated using $g_{1,K+1}^{n,{\rm{out}}}$.

4.2.2.   1–to–2 junction

The difference to the 1–to–1 junction is now that we have to consider two supplies. So the case distinctions to determine the flux values are a bit more complex. However, the procedure itself does not change. The details can be seen in algorithm 3.

We remark that if the inflow on the first road equals the demand and the restriction given by at least one of the supplies, the latter road needs to be congested. The Riemann solver states congestion such that the flux needs to be adjusted. Then, the incoming road either stays free flowing or the solution is given by $u^*$. Here, no adjustment is necessary. Note that the conditions $S_i\neq f(u^*-)$ are necessary for the specific case of $\beta_i = 1$.

If at least one of the supply restrictions is active, we need to adjust the corresponding flux values as in the 1–to–1 case and the incoming road. Nevertheless, here an interesting case (which is not possible in the 1–to–1 situation) can occur. The solution of the Riemann problem at the incoming road can be given by $u^*$ and the inflow is in $f^{{\rm{in}}}_1\in(f(u^*+),f(u^*-))$, see also the proof of theorem 3.7. We need to adjust the flux to $p(u^*) = f(u^*-)$. Hence, in this case if we use Eq (4.5) the adjustment is between $0$ and $-\alpha$. Finally, as in the 1–to–1 case there is no need for an adjustment if the demand is restrictive.

4.2.3.   2–to–1 junction

Recall that for two incoming and one outgoing roads, the maximal possible flux on the outgoing road is given by $f_{\rm{max}} = \min\{ D(u_{1,0}^{\rm{in}}) + D(u_{2,0}^{\rm{in}}), S(u_{1,0}^{\rm{out}})\}$. The corresponding flux values on the incoming roads are always smaller or equal to the actual demand. The algorithm for the corresponding adjusted flux values is listed in algorithm 4.

As before, no adjustment is needed when the demand on both roads is smaller than the supply. If the supply is active we might need to adjust the outgoing road and in most cases at least one incoming road. As in the 1–to–2 case, $f_i^{{\rm{in}}}$ can be greater than $f(u^*+)$ such that the Riemann solver of theorem 3.8 gives $u^*$ as a solution with the corresponding flux $p(u^*) = f(u^*-)$ leading to an adjustment smaller than $\alpha$.

5.   Numerical simulation

In this section, we present some numerical examples to compare the splitting algorithm on networks with the Riemann solver. Further, we compute the solution using a regularized flux and the Godunov scheme. We consider the following flux function:

The corresponding regularization is given by Eq (2.7).

We consider in particular the 1–to–2 and 2–to–1 situations. For our comparison, we choose constant initial data on each road. The junction is always located at $x = 0$, such that the incoming roads have negative $x$-coordinates and the outgoing ones positive. All roads are assumed to have the same length of 2. Furthermore, we compute the sum of the $L^1$ errors on each road at $t = 1$ to compare the numerical approaches. As the splitting algorithm requires $\lambda\leq 1$ as CFL condition, but the regularization approach $\lambda\leq \epsilon$, we compare different CFL constants and regularization parameters with each other.

5.1.   1–to–2 junction

In both scenarios the supply of the first outgoing road is restrictive. The parameter settings are as follows:

$(1)$ $u_{1,0}^{\rm{in}} = 0.4$, $u_{1,0}^{\rm{out}} = 0.9$, $u_{2,0}^{\rm{out}} = 0.7$, $\beta = (0.75, 0.25)$.

The exact solution is given by

Apparently, the solution induces two waves on the incoming road.

$(2)$ $u_{1,0}^{\rm{in}} = 0.4$, $u_{1,0}^{\rm{out}} = 0.7$, $u_{2,0}^{\rm{out}} = 0.2$, $\beta = (0.5, 0.5)$.

The exact solution is given by

This example generates $f_1^{{\rm{in}}} = 0.3$ such that we are in the specific case that the flux has to be adjusted, but only to $0.5$.

In Table 1, the $L^1$-errors are compared for different $\Delta x$ and the numerical convergence rate (CR) is determined by a least square fitting method. For the splitting algorithm we are allowed to choose a rather large CFL constant of 0.75 while for a direct comparison to the regularized approach $\lambda = 0.1$. is chosen. Furthermore, the regularization approach is also computed with a smaller regularization parameter. For the regularization approach we choose $\lambda = \epsilon$.

We can see that the error terms obtained by the splitting algorithm are the lowest and so the computational costs. Obviously, the error terms increase with a lower CFL due to numerical diffusion. For a direct comparison with the regularized approach, the CFL condition should be the same. Meaning that the second column in Table 1 for the splitting algorithm should be compared with the first one of the regularized approach. In this case, we can see that the splitting algorithm also performs better in both examples. By choosing a smaller regularization parameter the performance of the regularized approach increases, but also the computational costs. To obtain similar error terms as for the splitting algorithm the regularization parameter needs to be further reduced at very high computational costs.

For $\lambda = 0.1$ and $\Delta x = 0.01$ the approximate solutions and the exact solution are displayed in figure 6 for the first example. On the incoming road the splitting algorithm approximates the solution much better, while for the regularized approach an even smaller regularization parameter would be needed for a correct approximation. On the outgoing roads both schemes capture the right waves. Note that we zoomed in on the second outgoing road to make the small differences visible.

5.2.   2–to–1 junction

Here, we consider two scenarios, where in the first scenario the demand is restrictive while in the second one the supply. The parameter settings are as follows:

$(1)$ $u_{1,0}^{\rm{in}} = 0.2$, $u_{2,0}^{\rm{in}} = 0.25$, $u_{1,0}^{\rm{out}} = 0.3$, $q = 0.75$.

The exact solution is given by

$(2)$ $u_{1,0}^{\rm{in}} = 0.6$, $u_{2,0}^{\rm{in}} = 0.7$, $u_{1,0}^{\rm{out}} = 0.4$, $q = 0.8$.

The exact solution is given by

In particular, the flux value for the first incoming road needs to be adjusted from $0.4$ to $0.5$ in the splitting algorithm, while on the second road it is adjusted as usual from $0.1$ to $0.35 = 0.1+\alpha$. Note that due to the high velocity on the first incoming road we evaluate the error at $t = 0.5$.

Considering the $L^1$ errors in Table 2 similar observations as in the 1–to–2 case can be made. The splitting algorithm performs best and has the lowest computational costs. In the first example the solution of the splitting algorithm coincides with the regularized approach as the flows are equal at the junction due to the active demand, see second column of the splitting algorithm and first one of the regularized approach in Table 2. Nevertheless, in the second example the splitting algorithm performs much better for the same CFL constant. In addition, the regularization parameter needs to be even smaller than 0.01 to reach the precision of the faster splitting algorithm.

For $\lambda = 0.1$ and $\Delta x = 0.01$ the approximate solutions and the exact solution of the second example are displayed in Figure 7. In particular, on the first incoming road the splitting algorithm approximates the solution much better. For the considered resolution the difference between the exact solution and its approximation by the splitting algorithm is not visible on the first road. On the other roads the approximate solutions are very similar. Here, we zoomed in to display the differences to the exact solution.

6.   Conclusions

We have presented a Riemann solver at a junction for conservation laws with discontinuous flux. We have adapted the splitting algorithm of [22] to networks and demonstrated its validity in comparison with the exact solution. We have also pointed out that the splitting algorithm on networks is faster and more accurate than the approach with a regularized flux. Future research includes the investigation of other network models, where the flux is discontinuous.

Acknowledgments

J.F. was supported by the German Research Foundation (DFG) under grant HE 5386/18-1 and S.G. under grant GO 1920/10-1.

Conflict of interest

The authors declare there is no conflict of interest.