Well-posedness theory for nonlinear scalar conservation laws on networks

Ulrik Skre Fjordholm, Markus Musch and Nils Henrik Risebro

Abstract

We consider nonlinear scalar conservation laws posed on a network. We establish $L^{1}$ stability, and thus uniqueness, for weak solutions satisfying the entropy condition. We apply standard finite volume methods and show stability and convergence to the unique entropy solution, thus establishing existence of a solution in the process. Both our existence and stability/uniqueness theory is centred around families of stationary states for the equation. In one important case – for monotone fluxes with an upwind difference scheme – we show that the set of (discrete) stationary solutions is indeed sufficiently large to suit our general theory. We demonstrate the method’s properties through several numerical experiments.

1 Introduction

Partial differential equations (PDEs) on networks have a large number of applications, including fluid flow in pipelines, traffic flow on a network of roads, blood flow in vessels, data networks, irrigation channels and supply chains. A treatment of this wide range of applications can be found in the review articles [5, 11] and the references therein. In this paper we will focus on scalar, one-dimensional conservation laws

u_{t}+f(u)_{x}=0

(1.1)

on a network. Here, $u=u(x,t)$ is the unknown conserved variable and $f$ is a scalar flux function defined either on $\mathbb{R}$ or some subinterval. We aim to make sense of the conservation law on a directed graph and obtain existence, uniqueness, stability and approximability results.

Consider a network represented by a connected and directed graph. We tag the edges of this graph with an index $k$ and impose on each edge a scalar conservation law

\begin{split}u_{t}^{k}+f^{k}(u^{k})_{x}=0&,\qquad x\in D^{k},\ t>0\\ u^{k}(x,0)=\bar{u}^{k}(x)&,\qquad x\in D^{k}\end{split}

(1.2)

for some spatial domain $D^{k}\subset\mathbb{R}$ . (Here and in the remainder, a superscript $k$ will refer to an edge or a vertex.) We may think of edges as pipes or roads and the vertices as intersections, with the convention that the direction of travel is in the positive $x$ -direction, as shown in Figure 1.

Figure 1: A star shaped network with two ingoing and three outgoing edges.

Example 1.1.

A common one-dimensional model for traffic flow on a single road is described by (1.1) with $f(u)=s_{\max}u(u_{\max}-u)$ . Here $u=u(x,t)$ denotes the number of cars per length unit of road (density of vehicles), $u_{\max}$ denotes the maximal density of the road, and $s_{\max}$ the maximum speed. Commonly, the length unit is scaled so that $u_{\max}=1$ . This model is often referred to as the Lighthill–Whitham model [18].

A traffic model on a network might incorporate different speed limits $s_{\max}^{k}$ on each road section $k$ , or different numbers of lanes. In order to model the latter, each unknown $u^{k}$ can be scaled according to the capacity (proportional to the number of lanes) $\alpha^{k}$ on road $k$ , and the guiding principle to obtain the correct scaling should be the conservation of the total mass (i.e., the total number of cars) $\sum_{k}\int_{D^{k}}u^{k}(x,t)\,dx$ . Thus, if $\alpha^{k}$ denotes the capacity and $v^{k}(x,t)$ denotes the density of vehicles per lane, then the correct scaling is $u^{k}=\alpha^{k}v^{k}$ , and $u^{k}$ should satisfy (1.2) with

f^{k}(u^{k})=\alpha^{k}f\big{(}u^{k}/\alpha^{k}\big{)}

and e.g. $f(v)=v(1-v)$ .

In this paper we will be interested in uniqueness and stability for nonlinear scalar conservation laws on a network, as well as in constructing a numerical approximation and proving convergence of the numerical scheme. As opposed to many existing results, where the flux function on each edge is the same [8, 13], we want to allow for a different flux function $f^{k}$ on each edge $D^{k}$ of the network. Assuming that each flux $f^{k}$ is continuous and independent of the space variable, the problem can be seen as a PDE with a discontinuous flux, with the points of discontinuity sitting on the vertices of the graph. In fact, if our network would be the trivial network with only one ingoing and one outgoing edge then this would be exactly the case of a conservation law on the real line with a flux function with one discontinuity located at the vertex. Because of the connection to the theory for conservation laws with discontinuous fluxes (see e.g. [2]), we will borrow several ideas from this theory.

It is well-established that nonlinear hyperbolic conservation laws develop shocks in finite time. Therefore, solutions are always understood in the weak sense. Unfortunately, weak solutions to nonlinear hyperbolic conservation laws turn out to be non-unique, and additional conditions, usually referred to as entropy conditions, are imposed to select a unique solution. If the flux function is continuous then the theory of entropy solutions is covered by Kruzkhov’s theory [16]. For conservation laws with discontinuous fluxes the choice of entropy conditions is not obvious, and different physical models might lead to different entropy conditions. Although suitable entropy conditions can yield uniqueness, different entropy conditions are known to yield different solutions; see [2] and references therein. This study of different entropy conditions for conservation laws with discontinuous fluxes culminated in the paper by Andreianov, Karlsen and Risebro [2]. The authors relate the question of admissibility of a solution to the properties of a set of constant solutions, a so-called germ. Inspired by the entropy theory of Andreianov, Karlsen and Risebro, we investigate so-called stationary and discrete stationary solutions for our graph problem and thus derive an entropy theory for conservation laws on networks. Although our theory is influenced by the theory in [2], we have strived to make this paper as self-contained as possible.

In [13] the authors show uniqueness and existence to the Riemann problem as well as existence of a weak solution of the Cauchy problem on a network of roads in the case that the flux functions on each edge are identical. In [7, 6, 1, 8] the authors show well-posedness results for vanishing viscosity solutions. In [12] the authors investigate two entropy type conditions. However, in none of the existing literature can one find a general entropy condition which leads to uniqueness and stability of solutions. In the present work we aim to address this deficiency in the existing theory of conservation laws on networks.

The second important question to address is existence of a solution. Our approach will be to construct an approximation of the exact entropy solution by constructing a finite volume scheme. We will prove convergence to an entropy solution, thereby also proving existence of a solution.

In this article we focus on monotone fluxes – that is, each flux $f^{k}$ is either increasing or decreasing. For non-monotone fluxes only some of the techniques in the present paper are applicable; a traffic flow model, in which fluxes are assumed to be strictly concave, will be treated in a forthcoming work [10].

This article is structured as follows: In Section 2 we define our mathematical framework. We show uniqueness of entropy solutions to our problem in Section 3. In Section 4 we define a finite difference scheme appropriate for our problem, and in Section 5 we prove that our numerical scheme converges towards the unique entropy solution. In Section 6 we show that a class of monotone flux functions fits in our general scheme. Numerical experiments for the monotone case are presented in Section 7.

While the theory outlined in Sections 2 through 4 holds for conservation laws with general flux functions, the convergence theory in Sections 5 and 7 focuses on monotone flux functions and upwind numerical fluxes.

2 Entropy solutions on networks

Consider a network (or directed graph) of vertices and edges; for simplicity we will assume that the network contains a single vertex, along with $N_{{\mathrm{in}}}\in\mathbb{N}$ edges entering and $N_{{\mathrm{out}}}\in\mathbb{N}$ edges exiting the vertex (see Figure 1). (The generalization to general networks will follow analogously, due to the finite speed of propagation of the equations considered.) We think of the $N_{\mathrm{in}}$ edges as being to the left of the vertex and the $N_{\mathrm{out}}$ edges to its right. The ingoing edges will be labelled $k\in{{\mathscr{I}}_{\mathrm{in}}}\coloneqq\{-{N_{\mathrm{in}}},\dots,-1\}$ and the outgoing edges $k\in{{\mathscr{I}}_{\mathrm{out}}}\coloneqq\{1,\dots,{N_{\mathrm{out}}}\}$ . We denote $N\coloneqq{N_{\mathrm{in}}}+{N_{\mathrm{out}}}$ and we let ${\mathscr{I}}\coloneqq{{\mathscr{I}}_{\mathrm{in}}}\cup{{\mathscr{I}}_{\mathrm{out}}}$ denote the set of all edge indices. Placing the vertex at the origin $x=0$ , the incoming edges have coordinates $x\in\mathbb{R}_{-}=(-\infty,0)$ , while the outgoing edges have coordinates $x\in\mathbb{R}_{+}=(0,\infty)$ ; we will denote the $k$ -th edge by

D^{k}=\begin{cases}\mathbb{R}_{-}&\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},\\ \mathbb{R}_{+}&\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}.\end{cases}

On each edge $D^{k}$ we now impose the scalar conservation law (1.1), resulting in the $N$ distinct equations

\begin{split}u^{k}_{t}+f^{k}(u^{k})_{x}=0&\qquad\text{for }x\in D^{k},\ k\in{\mathscr{I}}.\end{split}

(2.1)

The collection of functions ${\mathbf{u}}=(u^{k})_{k\in{\mathscr{I}}}$ can be thought of as a function ${\mathbf{u}}\colon\Omega\to\mathbb{R}$ , where

\Omega\coloneqq\bigcup_{k\in{\mathscr{I}}}D^{k}\times\{k\}.

On the Borel $\sigma$ -algebra ${\mathscr{B}}(\Omega)\simeq\big{\{}\prod_{k\in{\mathscr{I}}}A^{k}\times\{k\}\ :\ A^{k}\in{\mathscr{B}}(D^{k})\big{\}}$ on $\Omega$ we define the measure $\lambda={\mathcal{L}}\times\#$ , where ${\mathcal{L}}$ is the one-dimensional Lebesgue measure and $\#$ is the counting measure; thus, the integral of ${\mathbf{u}}=(u^{k})_{k\in{\mathscr{I}}}$ is

\int_{\Omega}{\mathbf{u}}\,d\lambda=\sum_{k\in{\mathscr{I}}}\int_{D^{k}}u^{k}(x)\,dx.

(2.2)

The set of $L^{\infty}$ -bounded, real-valued functions on $\Omega$ will be denoted by $L^{\infty}(\Omega;\lambda)$ . We define the total variation of a function ${\mathbf{u}}\in L^{\infty}(\Omega;\lambda)$ as the sum of the variations of its components:

\operatorname{TV}({\mathbf{u}})\coloneqq\int_{\Omega}\Bigl{|}\frac{d{\mathbf{u}}}{dx}\Bigr{|}\,d\lambda=\sum_{k\in{\mathscr{I}}}\int_{D^{k}}\Bigl{|}\frac{du^{k}}{dx}(x)\Bigr{|}\,dx.

(2.3)

2.1 Weak solutions

Definition 2.1 (Weak Solution).

We say that a function ${\mathbf{u}}\in L^{\infty}\big{(}\mathbb{R}_{+};L^{\infty}(\Omega;\lambda)\big{)}$ is a weak solution of (2.1) with initial data $\bar{\mathbf{u}}\in L^{\infty}(\Omega;\lambda)$ if

\begin{split}\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}u^{k}\varphi^{k}_{t}+f^{k}(u^{k})\varphi^{k}_{x}\,dx\,dt+\sum_{k\in{\mathscr{I}}}\int_{D^{k}}\bar{u}^{k}(x)\varphi^{k}(x,0)\,dx=0\end{split}

(2.4)

for all $\varphi^{k}\in C_{c}^{\infty}\big{(}\overline{D^{k}}\times[0,\infty)\big{)}$ satisfying $\varphi^{k}(0,t)\equiv\varphi^{1}(0,t)$ for all $k\in{\mathscr{I}}$ .

Weak solutions automatically satisfy a Rankine–Hugoniot condition at the intersection:

Proposition 2.2 (Rankine–Hugoniot condition).

Let $(u^{k})_{k\in{\mathscr{I}}}$ be a weak solution of (2.1) such that $f^{k}\circ u^{k}(\cdot,t)$ has a strong trace at $x=0$ for every $k\in{\mathscr{I}}$ and a.e. $t>0$ . Then

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}(u^{k})(0,t)=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}(u^{k})(0,t)\qquad\text{for a.e.~{}}t>0.

(2.5)

Proof.

Define

\theta_{\varepsilon}(x)=\begin{cases}\frac{1}{\varepsilon}(\varepsilon+x)&\text{if }x\in[-\varepsilon,0]\\ \frac{1}{\varepsilon}(\varepsilon-x)&\text{if }x\in[0,\varepsilon]\\ 0&\text{if }|x|>\varepsilon.\end{cases}

(2.6)

We define $\Phi(x,t)\coloneqq\theta_{\varepsilon}(x)\psi(t)$ where $\psi\in C_{c}^{\infty}([0,\infty))$ . The partial derivatives of $\Phi$ are

\Phi_{x}(x,t)=\begin{cases}\frac{1}{\varepsilon}\psi(t)&\text{if }x\in[-\varepsilon,0]\\ -\frac{1}{\varepsilon}\psi(t)&\text{if }x\in[0,\varepsilon]\\ 0&\text{if }|x|>\varepsilon\end{cases}\qquad\text{and}\qquad\Phi_{t}(x,t)=\theta_{\varepsilon}(x)\psi^{\prime}(t).

By a density argument, $\Phi$ qualifies as an admissible test function. Thus, we can insert $\Phi$ into the weak formulation (2.4) to get

	$\displaystyle 0={}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}u^{k}\Phi_{t}^{k}+f^{k}(u^{k})\Phi_{x}^{k}\,dx\,dt+\sum_{k\in{\mathscr{I}}}\int_{D^{k}}\bar{u}^{k}(x)\Phi^{k}(x,0)\,dx$
	$\displaystyle={}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}u^{k}\theta_{\varepsilon}(x)\psi^{\prime}(t)\,dt\,dx$
		$\displaystyle+\frac{1}{\varepsilon}\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}\cap(-\varepsilon,\varepsilon)}\operatorname{sgn}(k)f^{k}(u^{k})\psi(t)\,dx\,dt$
		$\displaystyle-\frac{1}{\varepsilon}\sum_{k\in{\mathscr{I}}}\int_{D^{k}\cap(-\varepsilon,\varepsilon)}\bar{u}^{k}(x,0)\psi(t)\,\,dt$
	$\displaystyle\to{}$	$\displaystyle-\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\operatorname{sgn}(k)f^{k}(u^{k})\psi(t)\,dt$

as $\varepsilon\to 0$ , which is equivalent to (2.5). ∎

Definition 2.3 (Stationary Solution).

A stationary solution of (2.1) is a weak solution of (2.1) which is constant in time and is a strong solution on each edge $D^{k}$ . We see from (2.4) and (2.5) that the stationary solutions are precisely those satisfying $u^{k}(x,t)\equiv c^{k}\in\mathbb{R}$ for $x\in D^{k},\ t\geqslant 0$ and $k\in{\mathscr{I}}$ , and where $c^{k}$ satisfy the Rankine–Hugoniot condition

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}(c^{k})=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}(c^{k}).

(2.7)

Thus, we can identify each stationary solution with a vector ${\mathbf{c}}=(c^{k})_{k\in{\mathscr{I}}}\in\mathbb{R}^{N}$ .

Remark 2.4.

Note that if we only required stationary solutions to be weak solutions on each edge $D^{k}$ then they could exhibit arbitrarily many jump discontinuities. More precisely, if $f$ is not injective then a “stationary weak solution” could jump arbitrarily often between values $u^{k}=c^{k,\pm}$ , where $f(c^{k,-})=f(c^{k,+})$ .

2.2 Entropy conditions

Next, we formulate conditions that will single out a unique weak solution.

Definition 2.5 (Kruzkov entropy pairs).

The Kruzkov entropy pairs are the pairs of functions $\eta_{c}(u)=|u-c|$ , $q_{c}^{k}(u)=\operatorname{sgn}(u-c)\big{(}f^{k}(u)-f^{k}(c)\big{)}$ for $c\in\mathbb{R}$ .

The Kruzkov entropy pairs lead to a consistency condition on sets of stationary solutions:

Definition 2.6.

A subset $\mathscr{G}\subset\mathbb{R}^{N}$ consisting of stationary solutions of (2.1) is mutually consistent if

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}q_{c^{k}}^{k}\big{(}\tilde{c}^{k}\big{)}\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}q_{c^{k}}^{k}\big{(}\tilde{c}^{k}\big{)}

(2.8)

for every pair ${\mathbf{c}},\tilde{{\mathbf{c}}}\in\mathscr{G}$ , where $q_{c}^{k}$ are the Kruzkov entropy flux functions.

The set of stationary solutions $\mathscr{G}$ will determine what class of initial data we can consider:

Definition 2.7.

Let $\mathscr{G}\subset\mathbb{R}^{N}$ . We let $L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ be the set of $L^{\infty}$ -bounded data with values lying in the orthogonal convex hull^*^**This nomenclature is somewhat non-standard. The authors found several similar, but non-equivalent, definitions of the orthogonal convex hull in the literature. of $\mathscr{G}$ ,

\begin{split}L^{\infty}_{\mathrm{oco}}(\mathscr{G})=\Bigl{\{}&{\mathbf{u}}\in L^{\infty}(\Omega;\lambda)\ :\ \exists\ {\mathbf{c}},{\mathbf{d}}\in\mathscr{G}\text{ s.t. }[c^{k},d^{k}]\subset\pi^{k}(\mathscr{G})\ \forall\ k\in{\mathscr{I}}\\ &\text{and }c^{k}\leqslant u(x,k)\leqslant d^{k}\ \forall\ (x,k)\in\Omega\text{ and }c^{k}\leqslant d^{l}\ \forall\ k,l\in{\mathscr{I}}\Bigr{\}}\end{split}

(2.9)

where $\pi^{k}$ is the projection $\pi^{k}({\mathbf{c}})=c^{k}$ . (Note that we are ignoring the vertex index $k=0$ .)

The condition $[c^{k},d^{k}]\subset\pi^{k}(\mathscr{G})\ \forall\ k\in{\mathscr{I}}$ asserts that for every $\alpha\in[c^{k},d^{k}]$ , there is some $\tilde{\mathbf{c}}\in\mathscr{G}$ such that $\tilde{c}^{k}=\alpha$ . This property is needed in the doubling-of-variables argument in the proof of stability and uniqueness.

Example 2.8.

If $N_{\mathrm{in}}=N_{\mathrm{out}}=1$ and $f^{k}(u)=f(u)=u^{2}$ then the stationary solutions are all ${\mathbf{c}}\in\mathbb{R}^{2}$ of the form ${\mathbf{c}}=(c,c)$ or ${\mathbf{c}}=(c,-c)$ for $c\in\mathbb{R}$ . Both $\mathscr{G}_{1}=\{(c,c):c\in\mathbb{R}\}$ and $\mathscr{G}_{2}=\{(c,-c):c\geqslant 0\}$ (as well as any subset of these) are mutually consistent, as is $\mathscr{G}_{1}\cup\mathscr{G}_{2}$ . Note that no point of the form $(-c,c)$ for $c>0$ can be added to any of these sets and remain mutually consistent. Similarly, no set containing both $(-c,c)$ and $(-d,d)$ for distinct $c,d>0$ can be mutually consistent. The set $\mathscr{G}_{1}$ stands out as the smallest closed set which is such that both components span all of $\mathbb{R}$ , i.e. $\pi^{k}\mathscr{G}_{1}=\mathbb{R}$ for $k=-1,1$ . It is readily checked that $L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{1})=L^{\infty}(\Omega;\lambda)$ , and that any subset of $\mathscr{G}_{1}$ will yield a strictly smaller set of initial data. It is similarly straightforward to check that $\mathscr{G}_{2}$ generates a very restrictive set of initial data:

L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{2})=\big{\{}u\in L^{\infty}(\Omega;\lambda)\ :\ u^{-1}\equiv c,\ u^{1}\equiv-c\text{ for some }c\geqslant 0\big{\}}.

Thus, $\mathscr{G}_{1}$ is the smallest mutually consistent set of stationary solutions that allow for initial data in all of $L^{\infty}(\Omega;\lambda)$ .

Similar considerations hold when $f^{k}(u^{k})=\alpha^{k}f(u^{k}/\alpha^{k})$ for $\alpha_{-1}\neq\alpha_{1}$ (see Example 1.1).

Definition 2.9 (Entropy Solution).

Let $\mathscr{G}\subset\mathbb{R}^{N}$ be a mutually consistent set of stationary solutions of (2.1) and let $\bar{{\mathbf{u}}}\in L^{\infty}(\Omega;\lambda)$ . We say that a function ${\mathbf{u}}\in L^{\infty}(\mathbb{R}_{+},L^{\infty}(\Omega;\lambda))$ is an entropy solution of (2.1) with respect to $\mathscr{G}$ with initial data $\bar{\mathbf{u}}$ if

\begin{split}\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}\eta_{c^{k}}(u^{k})\varphi^{k}_{t}+q_{c^{k}}^{k}(u^{k})\varphi^{k}_{x}\,dx\,dt&\\ +\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\eta_{c^{k}}(\bar{u}^{k}(x))\varphi^{k}(x,0)\,dx&\geqslant 0\end{split}

(2.10)

for every ${\mathbf{c}}\in\mathscr{G}$ and every $0\leqslant\varphi\in C_{c}^{\infty}\big{(}\Omega\times[0,\infty)\big{)}$ satisfying $\varphi^{k}(0,t)\equiv\varphi^{1}(0,t)$ for all $k\in{\mathscr{I}}$ .

Audusse and Perthame [3] considered an entropy condition similar to (2.10), but in the context of spatially dependent, discontinuous flux functions.

We show first that entropy solutions are invariant in the set $L^{\infty}_{\mathrm{oco}}$ from Definition 2.7.

Lemma 2.10.

Let $\mathscr{G}\subset\mathbb{R}^{N}$ be a mutually consistent set of stationary solutions of (2.1) and let ${\mathbf{u}}\in L^{\infty}(\mathbb{R}_{+},L^{\infty}(\Omega;\lambda))$ be an entropy solution w.r.t. $\mathscr{G}$ with initial data $\bar{\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ . Then ${\mathbf{u}}(t)\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ for a.e. $t>0$ .

Proof.

Select ${\mathbf{c}},{\mathbf{d}}\in\mathscr{G}$ such that ${\mathbf{c}}\leqslant\bar{{\mathbf{u}}}\leqslant{\mathbf{d}}$ (cf. (2.9)). Add inequality (2.10) and equation (2.4) for both $c^{k}$ and $u^{k}$ to get

	$\displaystyle\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}\big{(}c^{k}-u^{k}\big{)}^{+}\varphi^{k}_{t}+H\big{(}c^{k}-u^{k}\big{)}\big{(}f^{k}(c^{k})-f^{k}(u^{k})\big{)}\varphi^{k}_{x}\,dx\,dt$
	$\displaystyle+\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\underbrace{\big{(}c^{k}-\bar{u}^{k}(x)\big{)}^{+}}_{=\,0}\varphi^{k}(x,0)\,dx$	$\displaystyle\geqslant 0$

(where $\cdot^{+}=\max(\cdot,0)$ and $H=\operatorname{sgn}^{+}$ is the Heaviside function). Replacing $\varphi$ by a sequence of approximations of the identity function on the set $\Omega\times[0,T]$ yields

\sum_{k\in{\mathscr{I}}}\int_{D^{k}}\big{(}c^{k}-u^{k}(x,T)\big{)}^{+}\,dx\leqslant 0

for a.e. $T>0$ , whence ${\mathbf{u}}(T)\geqslant{\mathbf{c}}$ . It follows similarly that ${\mathbf{u}}(T)\leqslant{\mathbf{d}}$ , and hence, ${\mathbf{u}}(T)\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ . ∎

The above lemma enables us to show that entropy solutions have strong traces.

Theorem 2.11.

Let $\mathscr{G}\subset\mathbb{R}^{N}$ be a mutually consistent set of stationary solutions of (2.1) and let ${\mathbf{u}}\in L^{\infty}(\mathbb{R}_{+},L^{\infty}(\Omega;\lambda))$ be an entropy solution w.r.t. $\mathscr{G}$ with initial data $\bar{\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ . Then the functions $q^{k}\circ u^{k}$ and $f^{k}\circ u^{k}$ admit strong traces on $\{x=0\}$ , for any $k\in{\mathscr{I}}$ .

Proof.

It follows from Lemma 2.10 that $u^{k}$ is a Kruzkhov entropy solution on $D^{k}$ , for any $k\in{\mathscr{I}}$ . We can therefore apply [19, Theorem 1.4] to obtain the desired conclusion. ∎

Proposition 2.12.

Let $\mathscr{G}\subset\mathbb{R}^{N}$ be a set of stationary solutions of (2.1). Let ${\mathbf{u}}$ be an entropy solution of (2.1) w.r.t. $\mathscr{G}$ such that $q_{c^{k}}^{k}\circ u^{k}(\cdot,t)$ has a strong trace at $x=0$ for every $k\in{\mathscr{I}}$ and a.e. $t>0$ . Then

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}q_{c^{k}}^{k}(u^{k})(0,t)\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}q_{c^{k}}^{k}(u^{k})(0,t)\qquad\text{for a.e.~{}}t>0

(2.11)

for every ${\mathbf{c}}\in\mathscr{G}$ .

Proof.

We take a positive test function $0\leqslant\psi\in C_{c}^{\infty}((0,\infty))$ . As in the proof of Proposition 2.2 we define $\Phi(x,t)\coloneqq\theta_{\varepsilon}(x,t)\psi(t)$ where $\theta_{\varepsilon}$ is given by (2.6). Now we insert $\Phi$ as test function into the entropy inequality (2.10) to get

	$\displaystyle 0\leqslant{}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}\eta_{c^{k}}(u^{k})\theta_{\varepsilon}(x)\psi^{\prime}(t)\,dx\,dt$
		$\displaystyle-\frac{1}{\varepsilon}\sum_{k\in{\mathscr{I}}}\biggl{(}\begin{aligned} &\int_{0}^{\infty}\int_{D^{k}\cap(-\varepsilon,\varepsilon)}\operatorname{sgn}(k)q_{c^{k}}^{k}(u^{k})\psi(t)\,dx\,dt\\ &+\int_{D^{k}}\eta_{c^{k}}(\bar{u}^{k}(x))\theta_{\varepsilon}(x)\psi(t)\,dx\biggr{)}\end{aligned}$
	$\displaystyle\to{}$	$\displaystyle-\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\operatorname{sgn}(k)q_{c^{k}}^{k}(u^{k}(0,t))\psi(t)\,dt$

as $\varepsilon\to 0$ , which shows the desired inequality. ∎

Definition 2.13 (Maximality).

A mutually consistent set of stationary solutions $\mathscr{G}$ is called maximal if there is no mutually consistent set of stationary solutions $\widetilde{\mathscr{G}}$ having $\mathscr{G}$ as a strict subset.

Lemma 2.14.

Let $\mathscr{G}$ be a maximal, mutually consistent set of stationary solutions and let ${\mathbf{d}}\in\mathbb{R}^{N_{\mathrm{in}}+N_{\mathrm{out}}}$ satisfy the Rankine–Hugoniot condition (2.7) as well as

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}q_{c^{k}}^{k}\big{(}d^{k}\big{)}\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}q_{c^{k}}^{k}\big{(}d^{k}\big{)}\qquad\forall\ {\mathbf{c}}\in\mathscr{G}.

Then ${\mathbf{d}}\in\mathscr{G}$ .

Proof.

By assumption, $\mathscr{G}\cup\{{\mathbf{d}}\}$ is a mutually consistent set of stationary solutions. But $\mathscr{G}$ is maximal, whence $\mathscr{G}\cup\{{\mathbf{d}}\}=\mathscr{G}$ . ∎

Remark 2.15.

By Proposition 2.12, we can take the vector ${\mathbf{d}}$ to be the trace of an entropy solution at the vertex. This observation will be important in the proof of the stability result in Section 3,

3 Stability and Uniqueness

Theorem 3.1 (Entropy Solutions are $L^{1}$ stable).

Let $\mathscr{G}\subset\mathbb{R}^{N}$ be a mutually consistent, maximal set of stationary solutions. Let ${\mathbf{u}},\,{\mathbf{v}}$ , be entropy solutions of (2.1) w.r.t. $\mathscr{G}$ with initial data $\bar{\mathbf{u}},\,\bar{\mathbf{v}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})\cap L^{1}(\Omega;\lambda)$ , and assume each function $u^{k},v^{k}$ has a strong trace at $x=0$ . Then

\sum_{k\in{\mathscr{I}}}\big{\|}u^{k}(t)-v^{k}(t)\big{\|}_{L^{1}(D^{k})}\leqslant\sum_{k\in{\mathscr{I}}}\big{\|}\bar{u}^{k}-\bar{v}^{k}\big{\|}_{L^{1}(D^{k})}

for every $t>0$ . In particular, there exists at most one entropy solution for given initial data.

Proof.

From Lemma 2.10 it follows that ${\mathbf{u}}(t),{\mathbf{v}}(t)\in L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ for a.e. $t>0$ . Let $k\in{{\mathscr{I}}_{\mathrm{in}}}$ ; the case $k\in{{\mathscr{I}}_{\mathrm{out}}}$ will follow analogously. The first step is a standard doubling of variables argument on each edge $k\in{{\mathscr{I}}_{\mathrm{in}}}$ by selecting $\varphi^{k}\in C_{c}^{\infty}(D^{k}\times[0,\infty))$ and $\varphi_{l}\equiv 0$ for $l\neq k$ . The doubling of variables technique on a single edge gives:

\begin{split}\int_{D^{k}}\int_{0}^{\infty}\big{|}u^{k}(x,t)-v^{k}(x,t)\big{|}\varphi_{t}+q_{v}^{k}(u)\varphi_{x}\,dt\,dx&\\ +\int_{D^{k}}\big{|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{|}\varphi(x,0)\,dx&\geqslant 0.\end{split}

(3.1)

Next, for general $\varphi^{k}\in C_{c}^{\infty}\big{(}\overline{D^{k}}\times[0,\infty)\big{)}$ , we cut off the functions near $x=0$ and couple the terms (3.1) on each edge together by utilizing (2.11). For $h>0$ we define

\mu_{h}(x)\coloneqq\begin{cases}0&x\in(-\infty,-2h)\\ \frac{1}{h}(x+2h)&x\in[-2h,-h)\\ 1&x\in[-h,0]\end{cases}

and

\Psi_{h}(x)\coloneqq 1-\mu_{h}(x).

The derivative of $\Psi_{h}$ reads

\Psi_{h}^{\prime}(x)=\begin{cases}0&x\in(-\infty,-2h)\\ -\frac{1}{h}&x\in[-2h,-h)\\ 0&x\in[-h,0].\end{cases}

Define $\varphi^{k}(x,t)\coloneqq\xi^{k}(x,t)\Psi_{h}(x)$ for a function $\xi^{k}\in C_{c}^{\infty}\big{(}\overline{D^{k}}\times[0,\infty)\big{)}$ . We insert $\varphi$ into equation (3.1) to get

	$\displaystyle\int_{D^{k}}\int_{0}^{\infty}\big{\|}u^{k}(x,t)-v^{k}(x,t)\big{\|}\xi_{t}^{k}\Psi_{h}+q^{k}_{v^{k}}(u^{k})\xi_{x}^{k}\Psi_{h}\,dt\,dx$
	$\displaystyle+\int_{D^{k}}\int_{0}^{\infty}q^{k}_{v^{k}}(u^{k})\xi^{k}\Psi_{h}^{\prime}\,dt\,dx+\int_{D^{k}}\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}\xi^{k}\Psi_{h}\,dx$	$\displaystyle\geqslant 0.$

Sending $h\downarrow 0$ we get

	$\displaystyle\int_{D^{k}}\int_{0}^{\infty}\big{\|}u^{k}(x,t)-v^{k}(x,t)\big{\|}\xi_{t}^{k}+q^{k}_{v^{k}}(u^{k})\xi_{x}^{k}\,dt\,dx+\int_{D^{k}}\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}\xi^{k}\,dx$
	$\displaystyle+\lim_{h\downarrow 0}\int_{-2h}^{-h}\int_{0}^{\infty}q^{k}_{v^{k}}(u^{k})\xi^{k}\Psi_{h}^{\prime}\,dt\,dx$	$\displaystyle\geqslant 0.$

Since the traces of $q^{k}(u^{k})$ and $q^{k}(v^{k})$ exist, we get

\displaystyle-\lim_{h\downarrow 0}\frac{1}{h}\int_{0}^{T}\int_{-2h}^{-h}q^{k}_{v^{k}}(u^{k})\xi^{k}\,dx\,dt=-\int_{0}^{T}q_{v^{k-}}^{k}(u^{k})\xi^{k}(0,t)\,dt.

We therefore obtain

\begin{split}\int_{D^{k}}\int_{0}^{\infty}\big{|}u^{k}(x,t)-v^{k}(x,t)\big{|}\xi_{t}+q^{k}_{v^{k}}(u^{k})\xi_{x}\,dt\,dx&\\ +\int_{D^{k}}\big{|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{|}\,dx-\int_{0}^{T}q_{v^{k}}^{k}(u^{k})\xi(0,t)\,dt&\geqslant 0.\end{split}

(3.2)

By an analogous argument we get

	$\displaystyle\int_{D^{k}}\int_{0}^{\infty}\big{\|}u^{k}(x,t)-v^{k}(x,t)\big{\|}\xi_{t}+q^{k}_{v^{k}}(u^{k})\xi_{x}\,dt\,dx$
	$\displaystyle+\int_{D^{k}}\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}\,dx+\int_{0}^{T}q_{v^{k}}^{k}(u^{k})\xi(0,t)\,dt$	$\displaystyle\geqslant 0$

for $k\in{{\mathscr{I}}_{\mathrm{out}}}$ . Fix $s>0$ , let $r,\kappa>0$ , and let $\alpha_{r}:\mathbb{R}_{-}\to\mathbb{R}$ and $\beta_{\kappa}\colon\mathbb{R}_{+}\to\mathbb{R}$ be given by

\alpha_{r}(x)=\begin{cases}0&x\in(-\infty,-r-1]\\ x+r+1&x\in(-r-1,-r)\\ 1&x\in[-r,0)\end{cases}

\beta_{\kappa}(t)=\begin{cases}1&t\in[0,s]\\ \frac{1}{\kappa}(\kappa+s-t)&t\in(s,s+\kappa)\\ 0&t\in[s+\kappa,\infty).\end{cases}

Via a standard regularization argument one can check that $\varphi(x,t)=\alpha_{r}(x)\beta_{\kappa}(t)$ is an admissible test function. We compute the partial derivatives of $\varphi$ :

\varphi_{t}(x,t)=\begin{cases}0&t\in[0,s]\\ -\frac{1}{\kappa}\alpha_{r}(x)&t\in(s,s+\kappa]\\ 0&t\in(s+\kappa,\infty)\end{cases}

and

\varphi_{x}(x,t)=\begin{cases}0&x\in(-\infty,-r-1)\\ \beta_{\kappa}(t)&x\in(-r-1,-r)\\ 0&x\in(-r,0).\end{cases}

We insert this into (3.2) to get

	$\displaystyle-\frac{1}{\kappa}\int_{s}^{s+\kappa}\int_{-r-1}^{0}\big{\|}u^{k}(x,t)-v^{k}(x,t)\big{\|}\alpha_{r}(x)\,dx\,dt$
	$\displaystyle+\int_{0}^{s+\kappa}\int_{-r-1}^{-r}q_{v^{k}}^{k}(u^{k})\beta_{\kappa}(t)\,dx\,dt+\int_{-r-1}^{0}\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}\alpha_{r}(x)\,dx$
	$\displaystyle-\int_{0}^{s+\kappa}q_{v^{k}}^{k}(u^{k}(0,t))\beta_{\kappa}(t)\,dt$	$\displaystyle\geqslant 0.$

Letting $\kappa\to 0$ and $r\to\infty$ , we get

\big{\|}u^{k}(x,t)-v^{k}(x,t)\big{\|}_{L^{1}(D^{k})}\leqslant\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}_{L^{1}(D^{k})}-\int_{0}^{s}q_{v^{k}}^{k}(u^{k}(0,t))\,dt.

An analogous inequality holds for $k\in{{\mathscr{I}}_{\mathrm{out}}}$ . We sum over all edges to get

	$\displaystyle\sum_{k\in{\mathscr{I}}}\big{\\|}u^{k}(x,t)-v^{k}(x,t)\big{\\|}_{L^{1}(D^{k})}$
	$\displaystyle\leqslant\sum_{k\in{\mathscr{I}}}\int_{D^{k}}\big{\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\|}\,dx+\underbrace{\int_{0}^{s}\sum_{k\in{\mathscr{I}}}\operatorname{sgn}(k)q^{k}_{v^{k}}(u^{k})}_{\begin{subarray}{c}\text{$\leqslant 0$ by \eqref{eq:RHentr} and}\\ \text{ maximality of $\mathscr{G}$}\end{subarray}}$
	$\displaystyle\leqslant\sum_{k\in{\mathscr{I}}}\big{\\|}\bar{u}^{k}(x)-\bar{v}^{k}(x)\big{\\|}_{L^{1}(D^{k})}.$

∎

4 Numerical approximation

In this section we construct a finite volume numerical approximation for (2.1) and prove stability and convergence properties of the method. The numerical method is rather standard for hyperbolic conservation laws, but an important feature of the method is that the vertex is discretized as a separate control volume. Although this control volume vanishes as the mesh parameter ${\Delta x}$ is passed to zero, its presence will ensure that entropy is correctly dissipated at the vertex, even in the limit ${\Delta x}\to 0$ .

4.1 A finite volume method on networks

Let ${\Delta t},{\Delta x}>0$ be given discretization parameters. We define the index sets

\displaystyle D_{\mathrm{disc}}^{+}\coloneqq\mathbb{N},\qquad D_{\mathrm{disc}}^{-}\coloneqq-\mathbb{N},\qquad D_{\mathrm{disc}}^{k}\coloneqq D_{\mathrm{disc}}^{\operatorname{sgn}(k)},\qquad D_{\mathrm{disc}}^{0}\coloneqq\{0\}.

For $n\in\mathbb{N}_{0}$ we discretize^†^††In numerical experiments, the timestep ${\Delta t}$ is chosen dynamically at each step $n$ in order to comply with the CFL condition derived in Section 4. We use a uniform timestep for the sake of simplicity only. $t^{n}=n{\Delta t}$ , and for $k\in{\mathscr{I}}$ and $i\in\mathbb{Z}$ we let $x_{i+{\nicefrac{{1}}{{2}}}}=({i+{\nicefrac{{1}}{{2}}}}){\Delta x}$ , and partition the physical domain into cells

{\mathscr{C}}_{i}^{k}=D^{k}\cap\big{(}x_{i-{\nicefrac{{1}}{{2}}}},\,x_{i+{\nicefrac{{1}}{{2}}}}\big{)}.

We define the mesh size at the vertex by ${\Delta x}_{0}\coloneqq\sum_{k\in{\mathscr{I}}}|{\mathscr{C}}_{0}^{k}|=N{\Delta x}/2$ , where $|A|$ denotes the Lebesgue measure of $A\subset\mathbb{R}$ . We make the finite volume approximation

	$\displaystyle u_{i}^{k,n}$	$\displaystyle\approx\frac{1}{{\Delta x}}\int_{{\mathscr{C}}_{i}^{k}}u^{k}(x,t^{n})\,dx\qquad\text{for }i\in D^{k}_{\mathrm{disc}},$
	$\displaystyle u_{0}^{n}$	$\displaystyle\approx\frac{1}{{\Delta x}_{0}}\sum_{k\in{\mathscr{I}}}\int_{{\mathscr{C}}_{0}^{k}}u^{k}(x,t^{n})\,dx.$

Fix some $i\in D_{\mathrm{disc}}^{k}$ , let $\varphi^{k}(x,t)=\frac{1}{{\Delta t}{\Delta x}}\mathbbm{1}_{{\mathscr{C}}_{i}^{k}}(x)\mathbbm{1}_{[t^{n},\,t^{n+1})}(t)$ and $\varphi^{l}\equiv 0$ for $l\neq k$ , and (after an approximation procedure) insert these into (2.4). We then obtain the numerical method

	$\frac{u_{i}^{k,n+1}-u_{i}^{k,n}}{{\Delta t}}+\frac{F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,n}}{{\Delta x}}=0$	(4.1a)
where $F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}=F^{k}\big{(}u^{k,n}_{i},u^{k,n}_{i+1}\big{)}$ is an approximation of the mean flux through $x_{i+{\nicefrac{{1}}{{2}}}}$ over the time interval $[t^{n},t^{n+1})$ ,
	$F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}\approx\frac{1}{{\Delta t}}\int_{t^{n}}^{t^{n+1}}f^{k}\big{(}u^{k}(x_{i+{\nicefrac{{1}}{{2}}}},t)\big{)}\,dt.$
For the special cell $i=0$ we let $\varphi^{k}(x,t)=\frac{1}{{\Delta t}{\Delta x}_{0}}\mathbbm{1}_{{\mathscr{C}}_{0}^{k}}(x)\mathbbm{1}_{[t^{n},\,t^{n+1})}(t)$ for $k\in{\mathscr{I}}$ to obtain
	$\frac{u_{0}^{n+1}-u_{0}^{n}}{{\Delta t}}+\frac{1}{{\Delta x}_{0}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{\nicefrac{{1}}{{2}}}^{k,n}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k,n}\biggr{)}=0.$	(4.1b)

We will use the notational convention that $u_{0}^{k,n}\equiv u_{0}^{n}$ for all $k\in{\mathscr{I}}$ . (We postpone the definition of the initial data $u_{i}^{k,0}$ until Section 4.3.)

Given a numerically computed solution $({\mathbf{u}}_{i}^{n})_{i,n}$ , we define the piecewise constant function

{\mathbf{u}}_{\Delta t}(x,k,t)=u_{i}^{k,n}\qquad\text{for }x\in{\mathscr{C}}_{i}^{k},\ t\in[t^{n},t^{n+1}).

(4.2)

We remark that the integral of ${\mathbf{u}}_{\Delta t}$ w.r.t. the measure $\lambda$ (cf. (2.2)) can be written

\int_{\Omega}{\mathbf{u}}_{\Delta t}(\cdot,t)\,d\lambda=\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}u_{i}^{k,n}\,{\Delta x}+u_{0}^{n}\,{\Delta x}_{0}

(4.3)

for any $t\in[t^{n},t^{n+1})$ , and the total variation of ${\mathbf{u}}_{\Delta t}$ (cf. (2.3)) can be written

\begin{split}\operatorname{TV}({\mathbf{u}}_{\Delta t}(\cdot,t))=&\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{|}u_{i+1}^{k,n}-u_{i}^{k,n}\big{|}+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{|}u_{i}^{k,n}-u_{i-1}^{k,n}\big{|}\\ =&\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{|}u_{i}^{k,n}-u_{i-1}^{k,n}\big{|}+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{|}u_{i+1}^{k,n}-u_{i}^{k,n}\big{|}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\big{|}u_{0}^{n}-u_{-1}^{k,n}\big{|}+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\big{|}u_{0}^{n}-u_{1}^{k,n}\big{|}.\end{split}

(4.4)

Note also that a numerical method of the form (4.1) is conservative in the sense that the total mass $\sum_{k\in{\mathscr{I}}}\int_{\Omega}{\mathbf{u}}_{\Delta t}\,d\lambda$ is independent of $n$ :

	$\displaystyle\int_{\Omega}{\mathbf{u}}_{\Delta t}(\cdot,t^{n+1})\,d\lambda={}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}u_{i}^{k,n+1}{\Delta x}+u_{0}^{n+1}{\Delta x}_{0}$
	$\displaystyle={}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}u_{i}^{k,n}{\Delta x}-{\Delta t}\big{(}F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\big{)}+u_{0}^{n}{\Delta x}_{0}$
		$\displaystyle{{}-{\Delta t}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{{\nicefrac{{1}}{{2}}}}^{k}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k}\biggr{)}$
	$\displaystyle={}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}u_{i}^{k,n}{\Delta x}+u_{0}^{n}{\Delta x}_{0}$
	$\displaystyle={}$	$\displaystyle\int_{\Omega}{\mathbf{u}}_{\Delta t}(\cdot,t^{n})\,d\lambda.$

As a shorthand for the scheme (4.1) we define the functions

	$G^{k}\big{(}u_{i-1}^{k},u_{i}^{k},u_{i+1}^{k}\bigr{)}\coloneqq u_{i}^{k}-\frac{{\Delta t}}{{\Delta x}}\big{(}F^{k}\bigl{(}u_{i}^{k},u_{i+1}^{k}\bigr{)}-F^{k}\bigl{(}u_{i-1}^{k},u_{i}^{k}\bigr{)}\big{)}$	(4.5a)
for $k\in{\mathscr{I}}$ and
	$\begin{split}&G^{0}\bigl{(}u_{-1}^{-{N_{\mathrm{in}}}},\dots,u_{-1}^{-1},u_{0},u_{1}^{1},\dots,u_{1}^{{N_{\mathrm{out}}}}\bigr{)}\\ &\quad\coloneqq u_{0}-\frac{{\Delta t}}{{\Delta x}_{0}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F^{k}\bigl{(}u_{0},u_{1}^{k}\bigr{)}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F^{k}\bigl{(}u_{-1}^{k},u_{0}\bigr{)}\biggr{)},\end{split}$	(4.5b)

enabling us to write (4.1) in the update form

\begin{split}u_{i}^{k,n+1}&=G^{k}\big{(}u_{i-1}^{k,n},u_{i}^{k,n},u_{i+1}^{k,n}\bigr{)}\qquad\qquad\text{for }i\in D_{\mathrm{disc}}^{k},\ k\in{\mathscr{I}}\\ u_{0}^{n+1}&=G^{0}\bigl{(}u_{-1}^{-{N_{\mathrm{in}}},n},\dots,u_{-1}^{-1,n},u_{0}^{n},u_{1}^{1,n},\dots,u_{1}^{{N_{\mathrm{out}}},n}\bigr{)}.\end{split}

(4.6)

As a shorthand for (4.6), we will sometimes use the notation

u_{i}^{k,n+1}=G^{k}\big{(}{\mathbf{u}}_{i-1}^{n},{\mathbf{u}}_{i}^{n},{\mathbf{u}}_{i+1}^{n}\big{)}\qquad\text{for }i\in D_{\mathrm{disc}}^{k},\ k\in{\mathscr{I}}_{0},

(4.6’)

where ${\mathbf{u}}_{i}^{n}$ is the vector containing all numerical values at index $i$ at time $n$ .

Definition 4.1 (Monotone scheme).

The difference scheme (4.6’) is monotone if

{\mathbf{u}}^{n}\leqslant{\mathbf{v}}^{n}\quad\Rightarrow\quad{\mathbf{u}}^{n+1}\leqslant{\mathbf{v}}^{n+1},

where ${\mathbf{u}}^{n}\leqslant{\mathbf{v}}^{n}$ means that every component $u_{i}^{k,n}$ of ${\mathbf{u}}^{n}$ is no greater than the corresponding component of ${\mathbf{v}}^{n}$ .

We state a straightforward CFL-type condition which ensures monotonicity of the numerical scheme.

Proposition 4.2.

Consider a consistent finite volume method (4.1), where $F^{k}$ is nondecreasing in the first variable and nonincreasing in the second. Then the method is monotone under the CFL condition

{\Delta t}\max_{k,u,v}\Bigl{|}\frac{\partial F^{k}}{\partial u}(u,v)\Bigr{|}\leqslant{\Delta x}/2,\qquad{\Delta t}\max_{k,u,v}\Bigl{|}\frac{\partial F^{k}}{\partial v}(u,v)\Bigr{|}\leqslant{\Delta x}/2.

(4.7)

Proof.

We can calculate the derivatives to the update functions to get

	$\displaystyle\frac{\partial G^{k}}{\partial u_{i-1}^{k}}$	$\displaystyle=\frac{{\Delta t}}{{\Delta x}}\frac{\partial F_{i-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i-1}^{k}},\qquad\frac{\partial G^{k}}{\partial u_{i+1}^{k}}=-\frac{{\Delta t}}{{\Delta x}}\frac{\partial F_{i+{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i+1}^{k}}$
	$\displaystyle\frac{\partial G^{k}}{\partial u_{i}^{k}}$	$\displaystyle=1-\frac{{\Delta t}}{{\Delta x}}\biggl{(}\frac{\partial F_{i+{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i}^{k}}-\frac{\partial F_{i-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i}^{k}}\biggr{)},$

for each $k\in{\mathscr{I}}$ , and

$\displaystyle\frac{\partial G^{0}}{\partial u_{-1}^{k}}$	$\displaystyle=\frac{{\Delta t}}{{\Delta x}_{0}}\frac{\partial F_{-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{-1}^{k}}$	$\displaystyle\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},$
$\displaystyle\frac{\partial G^{0}}{\partial u_{1}^{k}}$	$\displaystyle=-\frac{{\Delta t}}{{\Delta x}_{0}}\frac{\partial F_{{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{1}^{k}}$	$\displaystyle\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}},$
$\displaystyle\frac{\partial G^{0}}{\partial u_{0}}$	$\displaystyle=1-\frac{{\Delta t}}{{\Delta x}_{0}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\frac{\partial F_{{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{0}^{n}}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\frac{\partial F_{-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{0}^{n}}\biggr{)}$

on the vertex. We would like these derivatives to be non-negative. The monotonicity of $F^{k}$ guarantees that the first, second, fourth and fifth expressions are non-negative. Applying monotonicity of $F^{k}$ to the third and sixth terms, we get

\frac{\partial G^{k}}{\partial u_{i}^{k}}=\frac{{\Delta t}}{{\Delta x}}\biggl{(}\frac{{\Delta x}}{{\Delta t}}-\Bigl{|}\frac{\partial F_{i+{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i}^{k}}\Bigr{|}-\Bigl{|}\frac{\partial F_{i-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{i}^{k}}\Bigr{|}\biggr{)}\geqslant 0

(by (4.7)) and

	$\displaystyle\frac{\partial G^{0}}{\partial u_{0}}$	$\displaystyle=\frac{{\Delta t}}{{\Delta x}_{0}}\biggl{(}\frac{{\Delta x}_{0}}{{\Delta t}}-\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\Bigl{\|}\frac{\partial F_{{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{0}^{n}}\Bigr{\|}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\Bigl{\|}\frac{\partial F_{-{\nicefrac{{1}}{{2}}}}^{k}}{\partial u_{0}^{n}}\Bigr{\|}\biggr{)}$
(using ${\Delta x}_{0}=N{\Delta x}/2$ )
		$\displaystyle\geqslant\frac{{\Delta t}}{{\Delta x}_{0}}\biggl{(}N\frac{{\Delta x}}{2{\Delta t}}-{N_{\mathrm{out}}}\max_{k,u,v}\Bigl{\|}\frac{\partial F^{k}}{\partial u}(u,v)\Bigr{\|}-{N_{\mathrm{in}}}\max_{k,u,v}\Bigl{\|}\frac{\partial F^{k}}{\partial v}(u,v)\Bigr{\|}\biggr{)}$
		$\displaystyle\geqslant 0$

by (4.7). ∎

4.2 Discrete stationary solutions

In the same way that stationary solutions are essential for the well-posedness of entropy solutions (cf. Section 3), they are essential to the stability and convergence of numerical methods on networks. Asserting that a numerical solution is both constant in time and on each edge yields the following definition.

Definition 4.3 (Discrete Stationary Solution).

Consider a consistent, conservative numerical method (4.1). A discrete stationary solution for (4.1) is a vector

{\mathbf{c}}^{\mathrm{disc}}\coloneqq(c^{-{N_{\mathrm{in}}}},\dots,c^{{N_{\mathrm{out}}}})\in\mathbb{R}^{N+1}

satisfying the Rankine–Hugoniot condition

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}(c^{k})=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}(c^{k})

(4.8)

as well as the conditions


$\displaystyle F^{k}(c^{k},c^{0})$	$\displaystyle=f^{k}(c^{k})\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},$	(4.9a)
$\displaystyle F^{k}(c^{0},c^{k})$	$\displaystyle=f^{k}(c^{k})\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}.$	(4.9b)

In the remainder, sets of discrete stationary solutions will be denoted with a superscript, $\mathscr{G}^{0}$ , to signal that they also include a value at the vertex $i=0$ .

Notation 4.4.

We will sometimes index a discrete stationary solution as

	${\mathbf{c}}_{i}=\begin{cases}\big{(}c^{-{N_{\mathrm{in}}}},\dots,c^{-1}\big{)}&i<0\\ c^{0}&i=0\\ \big{(}c^{1},\dots,c^{N_{\mathrm{out}}}\big{)}&i>0\end{cases}$	(4.10a)
for $i\in\mathbb{Z}$ and, by extension,
	$c_{i}^{k}=\begin{cases}c^{k}&i\neq 0\\ c^{0}&i=0.\end{cases}$	(4.10b)

Using the notation (4.6’), it is readily checked that discrete stationary solutions are precisely those that are constant on each edge and satisfy

{\mathbf{c}}_{i}=G^{k}({\mathbf{c}}_{i-1},{\mathbf{c}}_{i},{\mathbf{c}}_{i+1})\qquad\forall\ i\in D_{\mathrm{disc}}^{k},\ k\in{\mathscr{I}}_{0}.

Remark 4.5.

The conditions (4.9) say that the numerical fluxes at the vertex reduce to the upwind flux on the in edges and the downwind flux on the out edges. This can be interpreted as information only flowing into the vertex, not out of it. This is consistent with the interpretation of the vertex as a stationary shock.

Remark 4.6.

Discrete stationary solutions ${\mathbf{c}}=\big{(}c^{-{N_{\mathrm{in}}}},\dots,c^{{N_{\mathrm{out}}}}\big{)}$ fulfil a discrete version of the Rankine–Hugoniot type condition (2.7),

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F^{k}\big{(}c^{k},c^{0}\big{)}=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F^{k}\big{(}c^{0},c^{k}\big{)}.

Lemma 4.7.

Consider a consistent, conservative numerical scheme (4.1). Let ${\mathbf{c}}=(c^{k})_{k\in{\mathscr{I}}}$ be a stationary solution for (1.1) and let $c^{0}\in\mathbb{R}$ . Then the vector ${\mathbf{c}}_{\mathrm{disc}}=(c^{-{N_{\mathrm{in}}}},\dots,c^{-1},c^{0},c^{1},\dots,c^{{N_{\mathrm{out}}}})\in\mathbb{R}^{N+1}$ is a discrete stationary solution if and only if

c^{0}\in\bigcap_{k\in{{\mathscr{I}}_{\mathrm{in}}}}(H^{k})^{-1}(\{f^{k}(c^{k})\})\mathrel{\bigcap}\bigcap_{k\in{{\mathscr{I}}_{\mathrm{out}}}}(J^{k})^{-1}(\{f^{k}(c^{k})\})

where

\displaystyle H^{k}(c)

\displaystyle\coloneqq F^{k}(c^{k},c)\quad\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},\qquad J^{k}(c)\coloneqq F^{k}(c,c^{k})\quad\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}.

Proof.

We can rewrite conditions (4.9a) and (4.9b) as

		$\displaystyle\eqref{eq:consistencyOutEdge}\quad\Leftrightarrow\quad H^{k}(c^{0})=f^{k}(c^{k})\quad\Leftrightarrow\quad c^{0}\in(H^{k})^{-1}(\{f^{k}(c^{k})\})$
for $k\in{{\mathscr{I}}_{\mathrm{in}}}$ , and
		$\displaystyle\eqref{eq:consistencyInEdge}\quad\Leftrightarrow\quad J^{k}(c^{0})=f^{k}(c^{k})\quad\Leftrightarrow\quad c^{0}\in(J^{k})^{-1}(\{f^{k}(c^{k})\})$

for $k\in{{\mathscr{I}}_{\mathrm{out}}}$ . Hence, if (4.9a), (4.9b) are satisfied then $c^{0}$ must lie in all of the sets on the right hand side, and hence in their intersection. Conversely, if $c^{0}$ lies in the intersection, then (4.9a), (4.9b) are satisfied. ∎

4.3 A stability framework for general consistent, conservative, monotone methods

We set out to prove an $L^{\infty}$ bound, $L^{1}$ contractiveness and Lipschitz continuity in time for solutions computed with a general consistent, conservative, monotone finite volume method on a network. Our starting point will be a class of discrete stationary solutions $\mathscr{G}_{\mathrm{disc}}^{0}\subset\mathbb{R}^{N+1}$ for a conservative finite volume method (4.1). We take initial data $\bar{u}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}}^{0})$ (cf. (2.9)), we let ${\mathbf{c}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ be as specified in (2.9), and consider the finite volume method (4.1) initialized by

u_{i}^{k,0}=\frac{1}{{\Delta x}}\int_{{\mathscr{C}}_{i}^{k}}\bar{u}^{k}(x)\,dx,\qquad u_{0}^{0}=c^{0}.

(4.11)

(The value $c^{0}$ is chosen for convenience, and any value in $[c^{0},d^{0}]$ will have the desired effect.)

Lemma 4.8.

Consider monotone numerical flux functions $F^{k}$ ( $k\in{\mathscr{I}}$ ). Let ${\mathbf{c}},{\mathbf{d}}$ be discrete stationary solutions satisfying $c^{k}\leqslant d^{l}$ for all $k,l\in{\mathscr{I}}$ (cf. Definition 2.7). Then $c^{0},d^{0}$ can be modified such that ${\mathbf{c}},{\mathbf{d}}$ remain discrete stationary solutions and such that $c^{0}\leqslant d^{0}$ .

Proof.

Define

I^{k}(c^{k})\coloneqq\begin{cases}\big{(}F^{k}(c^{k},\cdot)\big{)}^{-1}\big{(}\{f^{k}(c^{k})\}\big{)}&\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}}\\ \big{(}F^{k}(\cdot,c^{k})\big{)}^{-1}\big{(}\{f^{k}(c^{k})\}\big{)}&\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}.\end{cases}

Since all $F^{k}$ are monotone, each $I^{k}(c^{k})$ is a connected interval which contains $c^{k}$ , and moreover, Lemma 4.7 says that $c^{0}\in\bigcap_{k\in{\mathscr{I}}}I^{k}(c^{k})$ . This implies that $[\![c^{0},c^{k}]\!]\subset I^{k}(c^{k})$ , where $[\![a,b]\!]=[\min(a,b),\,\max(a,b)]$ . Since $\bigcap_{k\in{\mathscr{I}}}[\![c^{0},c^{k}]\!]$ is nonempty, the number

\tilde{c}^{0}\coloneqq\min\biggl{(}\bigcap_{k\in{\mathscr{I}}}[\![c^{0},c^{k}]\!]\biggr{)}

exists and satisfies $\tilde{c}^{0}\leqslant\max_{k\in{\mathscr{I}}}c^{k}$ . Appealing again to Lemma 4.7, ${\mathbf{c}}$ remains a discrete stationary solution if $c^{0}$ is replaced by $\tilde{c}^{0}$ . In a similar way we replace $d^{0}$ by

\tilde{d}^{0}\coloneqq\max\biggl{(}\bigcap_{k\in{\mathscr{I}}}[\![d^{0},d^{k}]\!]\biggr{)},

which satisfies $\tilde{d}^{0}\geqslant\min_{k\in{\mathscr{I}}}d^{k}$ . By our hypothesis, it follows that $\tilde{c}^{0}\leqslant\tilde{d}^{0}$ . ∎

Proposition 4.9.

Consider a consistent, conservative, monotone finite volume method (4.1), (4.11) with a set of discrete stationary states $\mathscr{G}_{\mathrm{disc}}^{0}$ . For any initial data $\bar{\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}}^{0})$ , the numerical solution is uniformly $L^{\infty}$ bounded.

Proof.

Pick discrete stationary states ${\mathbf{c}},{\mathbf{d}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ as in (2.9). It is clear that the initial data defined in (4.11) satisfy $c^{k}\leqslant u_{i}^{k,0}\leqslant d^{k}$ for all $i\in D_{\mathrm{disc}}^{k}$ and $k\in{\mathscr{I}}_{0}$ . If the same holds at some time step $n\in\mathbb{N}_{0}$ then (using the notation (4.6’), (4.10))

u_{i}^{k,n+1}=G^{k}\bigl{(}{\mathbf{u}}_{i-1}^{n},{\mathbf{u}}_{i}^{n},{\mathbf{u}}_{i+1}^{n}\bigr{)}\geqslant G^{k}\bigl{(}{\mathbf{c}}_{i-1},{\mathbf{c}}_{i},{\mathbf{c}}_{i+1}\bigr{)}=c_{i}^{k}

for all $i\in D_{\mathrm{disc}}^{k}$ and $k\in{\mathscr{I}}_{0}$ , and similarly, $u_{i}^{k,n+1}\leqslant d_{i}^{k}$ . ∎

Definition 4.10 ( $L^{1}$ contractive method).

A numerical method (4.6’) is $L^{1}$ contractive if

\displaystyle\bigl{\|}{\mathbf{u}}_{\Delta t}(\cdot,t)-{\mathbf{v}}_{\Delta t}(\cdot,t)\bigr{\|}_{L^{1}(\Omega;\lambda)}\leqslant\|\bar{\mathbf{u}}-\bar{\mathbf{v}}\|_{L^{1}(\Omega;\lambda)}

for all $t\geqslant 0$ , where ${\mathbf{u}}_{\Delta t}$ and ${\mathbf{v}}_{\Delta t}$ are the projection of the numerical solution (cf. (4.2)) computed with initial data $\bar{\mathbf{u}},\bar{\mathbf{v}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}}^{0})\cap L^{1}(\Omega;\lambda)$ , respectively. (See (4.3) for the integral of ${\mathbf{u}}_{\Delta t}$ , ${\mathbf{v}}_{\Delta t}$ w.r.t. $\lambda$ .)

We state the well known Crandall–Tartar lemma which we will use in the following proof. Here and below, we use the notation $a\vee b=\max(a,b)$ .

Theorem 4.11 (Crandall–Tartar: [9, Proposition 1]).

Let $(\Omega,\lambda)$ be a measure space. Let $C\subset L^{1}(\Omega;\lambda)$ have the property that $f,g\in C$ implies $f\vee g\in C$ . Let $V\colon C\to L^{1}(\Omega;\lambda)$ satisfy $\int_{\Omega}V(f)\,d\lambda=\int_{\Omega}f\,d\lambda$ for $f\in C$ . Then the following three properties of $V$ are equivalent:

(a)

$f,g\in C$ and $f\leqslant g$ a.e. implies $V(f)\leqslant V(g)$ a.e.,
(b)

$\int_{\Omega}(V(f)-V(g))^{+}\leqslant\int_{\Omega}(f-g)^{+}$ for $f,g\in C$ ,
(c)

$\int_{\Omega}\big{|}V(f)-V(g)\big{|}\leqslant\int_{\Omega}|f-g|$ for $f,g\in C$ .

We can now prove $L^{1}$ -contractivity of our method.

Theorem 4.12.

Every conservative, consistent monotone method (4.1), (4.11) is $L^{1}$ -contractive.

Proof.

Let $C={C_{\Delta x}}$ be the set of piecewise constant functions,

{C_{\Delta x}}=\Bigl{\{}{\mathbf{u}}\in L^{1}\cap L^{\infty}(\Omega;\lambda)\ :\ {\mathbf{u}}(x)=\sum_{k\in{\mathscr{I}}_{0}}\sum_{i\in D_{\mathrm{disc}}^{k}}u_{i}^{k}\mathbbm{1}_{{\mathscr{C}}_{i}^{k}}\text{ for }u_{i}^{k}\in\mathbb{R}\Bigr{\}}.

We define the operator $V:{C_{\Delta x}}\to{C_{\Delta x}}$ mapping a numerical solution to the next time step,

	$\displaystyle V({\mathbf{u}})\coloneqq{}$	$\displaystyle\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\mathbbm{1}_{{\mathscr{C}}_{i}^{k}}\Bigl{(}u_{i}^{k}-\frac{{\Delta t}}{{\Delta x}}\big{(}F^{k}\big{(}u_{i}^{k},u_{i+1}^{k}\big{)}-F^{k}\big{(}u_{i-1}^{k},u_{i}^{k}\big{)}\big{)}\Bigr{)}$
		$\displaystyle+\sum_{k\in{\mathscr{I}}}\mathbbm{1}_{{\mathscr{C}}_{0}^{k}}\biggl{(}u_{0}^{0}-\frac{{\Delta t}}{{\Delta x}_{0}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F^{k}(u_{0},u_{1}^{k})-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F^{k}(u_{-1}^{k},u_{0})\biggr{)}\biggr{)}.$

By the definition (2.2) of the measure $\lambda$ (cf. also (4.3)), we have $\int_{\Omega}V({\mathbf{u}})\,d\lambda=\int_{\Omega}{\mathbf{u}}\,d\lambda$ for all ${\mathbf{u}}\in{C_{\Delta x}}$ . We apply the Crandall–Tartar lemma to conclude $L^{1}$ -contractivity of the numerical solution operator $V$ . ∎

From $L^{1}$ -contractivity we get continuity in time as a corollary:

Corollary 4.13.

Consider a consistent, conservative and monotone method (4.1). Let ${\mathbf{u}}_{\Delta t}$ be an approximate solution computed with this method and assume that all numerical fluxes $F^{k}$ are Lipschitz continuous in both arguments. Then computed solutions are uniformly $L^{1}$ Lipschitz continuous in time, i.e.,

	$\displaystyle\big{\\|}{\mathbf{u}}_{\Delta t}(t^{n+1})-{\mathbf{u}}_{\Delta t}(t^{n})\big{\\|}_{L^{1}(\Omega;\lambda)}$	$\displaystyle\leqslant\big{\\|}{\mathbf{u}}_{\Delta t}(t^{1})-{\mathbf{u}}_{\Delta t}(t^{0})\big{\\|}_{L^{1}(\Omega;\lambda)}$
		$\displaystyle\leqslant{\Delta t}\big{(}C\operatorname{TV}({\mathbf{u}}^{0})+\bar{M}\big{)},$

where the constants $C$ and $\bar{M}$ depend on the flux functions and on the initial data.

Proof.

We compute

	$\displaystyle\bigl{\\|}$	$\displaystyle{\mathbf{u}}_{\Delta t}(t^{n+1})-{\mathbf{u}}_{\Delta t}(t^{n})\bigr{\\|}_{L^{1}(\Omega;\lambda)}$
		$\displaystyle=\big{\\|}V({\mathbf{u}}_{\Delta t}(t^{n}))-V({\mathbf{u}}_{\Delta t}(t^{n-1}))\big{\\|}_{L^{1}(\Omega;\lambda)}$
(using Theorem 4.11(c))
		$\displaystyle\leqslant\big{\\|}{\mathbf{u}}_{\Delta t}(t^{n})-{\mathbf{u}}_{\Delta t}(t^{n-1})\big{\\|}_{L^{1}(\Omega;\lambda)}\leqslant\dots\leqslant\big{\\|}{\mathbf{u}}_{\Delta t}(t^{1})-{\mathbf{u}}_{\Delta t}(t^{0})\big{\\|}_{L^{1}(\Omega;\lambda)}$
		$\displaystyle={\Delta x}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\|u_{i}^{k,1}-u_{i}^{0}\|+{\Delta x}_{0}\|u_{0}^{k,1}-u_{0}^{0}\|$
		$\displaystyle={\Delta t}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}F_{i+{\nicefrac{{1}}{{2}}}}^{k,0}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,0}\big{\|}+{\Delta t}\biggl{\|}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{\nicefrac{{1}}{{2}}}^{k,0}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k,0}\biggr{\|}$
		$\displaystyle={\Delta x}\lambda\sum_{\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,0}-u_{i-1}^{k,0}\big{\|}$
		$\displaystyle\quad+{\Delta t}\biggl{\|}\!\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{\nicefrac{{1}}{{2}}}^{k,0}-F^{k,0}\big{(}u_{0}^{0},u_{0}^{0}\big{)}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k,0}-F^{k,0}(u_{0}^{0},u_{0}^{0})\\ &+\overbrace{\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}\big{(}u_{0}^{0}\big{)}}^{\eqqcolon f_{\mathrm{out}}(u_{0}^{0})}-\overbrace{\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}\big{(}u_{0}^{0}\big{)}}^{\eqqcolon f_{\mathrm{in}}(u_{0}^{0})}\biggr{\|}\end{aligned}$
		$\displaystyle\leqslant{\Delta t}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}L^{k}\bigl{(}\big{\|}u_{i}^{k,0}-u_{i-1}^{k,0}\big{\|}+\big{\|}u_{i+1}^{k,0}-u_{i}^{k,0}\big{\|}\bigr{)}$
		$\displaystyle\quad+{\Delta t}\underbrace{\Biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}L^{k}\big{\|}u_{1}^{k,0}-u_{0}^{0}\big{\|}+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}L^{k}\big{\|}u_{0}^{0}-u_{-1}^{k,0}\big{\|}+\overbrace{\big{\|}f_{\mathrm{out}}\big{(}u_{0}^{0}\big{)}-f_{\mathrm{in}}\big{(}u_{0}^{0}\big{)}\big{\|}}^{\leqslant M}\Biggr{)}}_{\eqqcolon\bar{M}}$
		$\displaystyle\leqslant{\Delta t}(C\operatorname{TV}({\mathbf{u}}^{0})+\bar{M}),$

where we collect all constants into the global constant $C$ . We can bound $\big{|}f_{\mathrm{out}}\big{(}u_{0}^{0}\big{)}-f_{\mathrm{in}}\big{(}u_{0}^{0}\big{)}\big{|}\leqslant\bar{M}$ with a constant $\bar{M}\in\mathbb{R}$ since $f_{\mathrm{in}},f_{\mathrm{out}}$ are continuous and $u_{0}^{0}\in L^{\infty}$ . ∎

5 Convergence of finite volume schemes

We are now in place to prove convergence in the case where the flux functions $f^{k}$ are strictly monotone. We do this by using the upwind method where the numerical flux functions are defined by

F^{k}(u,v)=\begin{cases}f^{k}(u)&\text{if $f^{k}$ is increasing,}\\ f^{k}(v)&\text{if $f^{k}$ is decreasing.}\end{cases}

We shall show that the set of discrete approximations is compact in $L^{\infty}\big{(}[0,\infty);L^{1}(\Omega;\lambda)\big{)}$ , and that any limit limit is an entropy solution. In particular, this convergence result establishes existence of an entropy solution. We show convergence to a weak solution by proving a Lax–Wendroff type theorem:

Theorem 5.1 (Lax–Wendroff theorem).

Fix $T>0$ . Assume that each flux function $f^{k}$ is locally Lipschitz continuous and strictly monotone. Let $\mathscr{G}_{\mathrm{disc}}^{0}$ be a class of discrete stationary solutions for the upwind method and let ${\mathbf{u}}_{\Delta t}$ be computed from the upwind method with initial data $\bar{\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}}^{0})\cap L^{1}(\Omega;\lambda)$ . Consider a subsequence $\big{(}{\mathbf{u}}_{{\Delta t}_{\ell}}\big{)}_{\ell\in\mathbb{N}}$ such that ${\Delta t}_{\ell}\to 0$ and ${\mathbf{u}}_{{\Delta t}_{\ell}}\to{\mathbf{u}}$ in $L^{\infty}([0,T];L^{1}(\Omega;\lambda))$ as $k\to\infty$ . Then the limit ${\mathbf{u}}$ is the unique entropy solution to (2.1), that is, ${\mathbf{u}}$ satisfies (2.10).

Proof.

We write ${\Delta x}$ and ${\Delta t}$ rather than ${\Delta x}_{\ell}$ and ${\Delta t}_{\ell}$ , and we shall show that $u$ satisfies the entropy condition (2.10) for every ${\mathbf{c}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ . Choosing stationary solutions ${\mathbf{c}},{\mathbf{d}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ such that ${\mathbf{c}}\leqslant{\mathbf{u}}_{\Delta t}\leqslant{\mathbf{d}}$ (cf. Proposition 4.9) in particular shows that $u$ is a weak solution.

Let ${\mathbf{c}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ and consider the Crandall–Majda numerical entropy fluxes

\displaystyle Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}=F^{k}\bigl{(}u_{i}^{k,n}\vee c^{k}_{i},u_{i+1}^{k,n}\vee c^{k}_{i}\bigr{)}-F^{k}\bigl{(}u_{i}^{k,n}\wedge c^{k},u_{i+1}^{k,n}\wedge c^{k}\bigr{)}

for $i=0,1,\dots$ when $k\in{{\mathscr{I}}_{\mathrm{out}}}$ , and for $i=\dots,-2,-1$ when $k\in{{\mathscr{I}}_{\mathrm{in}}}$ , and

\displaystyle Q_{-{\nicefrac{{1}}{{2}}}}^{n}=\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}Q_{-{\nicefrac{{1}}{{2}}}}^{k,n},\qquad Q_{{\nicefrac{{1}}{{2}}}}^{n}=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}Q_{{\nicefrac{{1}}{{2}}}}^{k,n}

(cf. Notation 4.4 for the definition of $c_{i}^{k}$ ). Recalling the definition (4.5) of the update functions $G^{k}$ , we see that

	$\displaystyle G^{k}\bigl{(}u_{i-1}^{k,n}\vee c^{k},u_{i}^{k,n}\vee c^{k},u_{i+1}^{k,n}\vee c^{k}\bigr{)}-G^{k}\bigl{(}u_{i-1}^{k,n}\wedge c^{k},u_{i}^{k,n}\wedge c^{k},u_{i+1}^{k,n}\wedge c^{k}\bigr{)}$
	$\displaystyle\qquad=\bigl{\|}u_{i}^{k,n}-c^{k}\bigr{\|}-\frac{{\Delta t}}{{\Delta x}}\bigl{(}Q_{i+{\nicefrac{{1}}{{2}}}}^{k}-Q_{i-{\nicefrac{{1}}{{2}}}}^{k}\bigr{)},$

for $k\in{\mathscr{I}}$ and $i\in D_{\mathrm{disc}}^{k}$ . Hence,

\begin{split}\big{|}u_{i}^{k,n+1}-c^{k}\big{|}&=u_{i}^{k,n+1}\vee c^{k}-u_{i}^{k,n+1}\wedge c^{k}\\ &=G^{k}\bigl{(}u_{i-1}^{k,n},u_{i}^{k,n},u_{i+1}^{k,n}\bigr{)}\vee c^{k}-G^{k}\bigl{(}u_{i-1}^{k,n},u_{i}^{k,n},u_{i+1}^{k,n}\bigr{)}\wedge c^{k}\\ &\leqslant G^{k}\bigl{(}u_{i-1}^{k,n}\vee c^{k},u_{i}^{k,n}\vee c^{k},u_{i+1}^{k,n}\vee c^{k}\bigr{)}\\ &\quad-G^{k}\bigl{(}u_{i-1}^{k,n}\wedge c^{k},u_{i}^{k,n}\wedge c^{k},u_{i+1}^{k,n}\wedge c^{k}\bigr{)}\\ &=\bigl{|}u_{i}^{k,n}-c^{k}\bigr{|}-\frac{{\Delta t}}{{\Delta x}}\bigl{(}Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-Q_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\bigr{)}.\end{split}

(5.1)

Similarly, we find that

\bigl{|}u_{0}^{n+1}-c^{0}\bigr{|}-\bigl{|}u_{0}^{n}-c^{0}\bigr{|}+\frac{{\Delta t}}{{\Delta x}_{0}}\bigl{(}Q_{\nicefrac{{1}}{{2}}}^{n}-Q_{-{\nicefrac{{1}}{{2}}}}^{n}\bigr{)}\leqslant 0

(5.2)

We choose $T=M{\Delta t}$ for a natural number $M$ , multiply the above $N+1$ inequalities with a test function $\varphi$ and sum up to get

	$\displaystyle 0$	$\displaystyle\geqslant\sum_{n=0}^{M}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\varphi_{i}^{k,n}\biggl{(}\bigl{(}\big{\|}u_{i}^{k,n+1}-c^{k}\big{\|}-\big{\|}u_{i}^{k,n}-c^{k}\big{\|}\bigr{)}+\frac{{\Delta t}}{{\Delta x}}\bigl{(}Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-Q_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\bigr{)}\biggr{)}$
		$\displaystyle\quad+\sum_{n=0}^{M}\varphi_{0}^{n}\biggl{(}\frac{N}{2}\bigl{(}\big{\|}u_{0}^{n+1}-c^{0}\big{\|}-\big{\|}u_{0}^{n}-c^{0}\big{\|}\bigr{)}+\frac{{\Delta t}}{{\Delta x}}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}Q_{{\nicefrac{{1}}{{2}}}}^{k,n}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}Q_{-{\nicefrac{{1}}{{2}}}}^{k,n}\biggr{)}\biggr{)},$

where $\varphi_{i}^{k,n}=\varphi^{k}(x_{i},t^{n})$ . After summation by parts we get

	$\displaystyle 0$	$\displaystyle\geqslant-\sum_{n=1}^{M}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,n}-c^{k}\big{\|}\bigl{(}\varphi_{i}^{k,n}-\varphi_{i}^{k,n-1}\bigr{)}-\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,0}-c^{k}\big{\|}\varphi_{i}^{k,0}$
		$\displaystyle\quad-\frac{{\Delta t}}{{\Delta x}}\sum_{n=0}^{M}\biggl{(}\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\bigl{(}\varphi_{i}^{k,n}-\varphi_{i-1}^{k,n}\bigr{)}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}\bigl{(}\varphi_{i+1}^{k,n}-\varphi_{i}^{k,n}\bigr{)}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\varphi_{-1}^{k,n}Q_{-{\nicefrac{{1}}{{2}}}}^{k,n}-\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\varphi_{1}^{k,n}Q_{{\nicefrac{{1}}{{2}}}}^{k,n}\biggr{)}\end{aligned}$
		$\displaystyle\quad-\frac{N}{2}\big{\|}u_{0}^{n}-c^{0}\big{\|}\varphi_{0}^{0}-\frac{N}{2}\sum_{n=1}^{M}\big{\|}u_{0}^{n}-c^{0}\big{\|}\bigl{(}\varphi_{0}^{n}-\varphi_{0}^{n-1}\bigr{)}$
		$\displaystyle\quad+\frac{{\Delta t}}{{\Delta x}}\sum_{n=0}^{M}\biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}Q_{{\nicefrac{{1}}{{2}}}}^{n}\varphi_{0}^{n}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}Q_{-{\nicefrac{{1}}{{2}}}}^{n}\varphi_{0}^{n}\biggr{)}.$

After shifting the $i$ index on the second line we get

	$\displaystyle 0$	$\displaystyle\geqslant-\underbrace{{\Delta t}{\Delta x}\sum_{n=1}^{M}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,n}-c^{k}\big{\|}\biggl{(}\frac{\varphi_{i}^{k,n}-\varphi_{i}^{k,n-1}}{{\Delta t}}\biggr{)}}_{=\,A_{1}}$
		$\displaystyle\quad-\underbrace{{\Delta x}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,0}-c^{k}\big{\|}\varphi_{i}^{k,0}}_{=\,A_{2}}$
		$\displaystyle\quad-\underbrace{{\Delta t}{\Delta x}\sum_{n=0}^{M}\biggl{(}\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}\biggl{(}\frac{\varphi_{i+1}^{k,n}-\varphi_{i}^{k,n}}{{\Delta x}}\biggr{)}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\biggl{(}\frac{\varphi_{i}^{k,n}-\varphi_{i-1}^{k,n}}{{\Delta x}}\biggr{)}\biggr{)}\end{aligned}}_{=\,A_{3}}$
		$\displaystyle\quad-\underbrace{{\Delta x}\frac{N}{2}\big{\|}u_{0}^{0}-c^{0}\big{\|}\varphi_{0}^{0}-{\Delta t}{\Delta x}\frac{N}{2}\sum_{n=1}^{M}\big{\|}u_{0}^{n}-c^{0}\big{\|}\biggl{(}\frac{\varphi_{0}^{n}-\varphi_{0}^{n-1}}{{\Delta t}}\biggr{)}}_{=\,A_{4}}.$

Taking limits we get for ${\Delta t},{\Delta x}\to 0$

\displaystyle A_{1}\to\sum_{k\in{\mathscr{I}}}\int_{0}^{\infty}\int_{D^{k}}\big{|}u^{k}-c^{k}\big{|}\varphi_{t}^{k}\,dx\,dt,

and for ${\Delta x}\to 0$

\displaystyle A_{2}\to\sum_{k\in{\mathscr{I}}}\big{|}u_{0}^{k}-c^{k}\big{|}(x)\varphi^{k}(x,0)\,dx,\qquad A_{4}\to 0.

Thus, we are left with $A_{3}$ . Since the scheme is the upwind method, we can write

	$\displaystyle A_{3}$	$\displaystyle={\Delta t}{\Delta x}\sum_{n=0}^{M}\biggl{(}\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}q_{c^{k}}^{k}\bigl{(}u_{i}^{k,n}\bigr{)}\biggl{(}\frac{\varphi_{i+1}^{k,n}-\varphi_{i}^{k,n}}{{\Delta x}}\biggr{)}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}q_{c^{k}}^{k}\bigl{(}u_{i}^{k,n}\bigr{)}\biggl{(}\frac{\varphi_{i}^{k,n}-\varphi_{i-1}^{k,n}}{{\Delta x}}\biggr{)}\biggr{)}\end{aligned}$
		$\displaystyle\to\sum_{k\in{\mathscr{I}}}\int_{0}^{T}\int_{D^{k}}q_{c^{k}}^{k}\bigl{(}u^{k}(x,t)\bigr{)}\varphi_{x}(x,t)\,dx\,dt$

as ${\Delta t},{\Delta x}\to 0$ , due to the a.e. pointwise convergence of ${\mathbf{u}}_{{\Delta t}}$ to ${\mathbf{u}}$ . ∎

To show compactness we want to apply Helly’s theorem:

Theorem 5.2 (Helly’s theorem, [14, p. 437]).

Let $A\subset C([a,b])$ be a set of functions on an interval $[a,b]\subset\mathbb{R}$ for which there exists some $M>0$ such that

\operatorname{TV}(v)+\|v\|_{\infty}\leqslant M\qquad\forall\ v\in A.

Then $A$ is relatively compact with respect to uniform convergence.

Now we have everything in place to proof a compactness theorem.

Theorem 5.3 (Compactness and Convergence to a Weak Solution).

Fix $T>0$ . Assume that each flux function $f^{k}$ is locally Lipschitz continuous and strictly monotone. Let $\mathscr{G}_{\mathrm{disc}}^{0}$ be a set of discrete stationary states for the upwind method. Let ${\mathbf{u}}_{\Delta t}$ be computed from the upwind method with initial data $\bar{\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}}^{0})\cap L^{1}(\Omega;\lambda)$ , and assume that $\operatorname{TV}(\bar{\mathbf{u}})<\infty$ . Then the numerical solution $\{{\mathbf{u}}_{\Delta t}\}_{{\Delta t}>0}$ converges in $C([0,T],L_{\mathrm{loc}}^{1}(\Omega;\lambda))$ to a weak solution ${\mathbf{u}}$ .

Proof.

We first show convergence of the sequence of functions ${\mathbf{g}}_{\Delta t}\colon\Omega\times[0,T]\to\mathbb{R}$ ,

{\mathbf{g}}_{\Delta t}(x,k,t)\coloneqq f^{k}(u^{k}_{\Delta t}(x,t)).

The sequence ${\mathbf{g}}_{\Delta t}$ is uniformly $L^{\infty}$ bounded, by Proposition 4.9, and it is Lipschitz continuous in time:

	$\displaystyle\int_{\Omega}\big{\|}{\mathbf{g}}_{\Delta t}(t^{n+1})-{\mathbf{g}}_{\Delta t}(t^{n})\big{\|}\,d\lambda$	$\displaystyle\leqslant C_{f}\int_{\Omega}\big{\|}{\mathbf{u}}_{\Delta t}(t^{n+1})-{\mathbf{u}}_{\Delta t}(t^{n})\big{\|}\,d\lambda$
		$\displaystyle\leqslant C_{f}(C\operatorname{TV}(\bar{\mathbf{u}})+\bar{M}){\Delta t},$

by Corollary 4.13. We can bound the total variation of ${\mathbf{g}}_{\Delta t}$ by

	$\displaystyle\operatorname{TV}({\mathbf{g}}_{\Delta t}(\cdot,t))$	$\displaystyle=\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}f^{k}(u_{i}^{k,n})-f^{k}(u_{i-1}^{k,n})\big{\|}$
		$\displaystyle\quad+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\big{\|}f^{k}(u_{0}^{n})-f^{k}(u_{-1}^{k,n})\big{\|}$
		$\displaystyle\leqslant\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}f^{k}(u_{i}^{k,n})-f^{k}(u_{i-1}^{k,n})\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle=\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle=\frac{{\Delta x}}{{\Delta t}}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,n+1}-u_{i}^{k,n}\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle\leqslant{\Delta x}(C\operatorname{TV}(\bar{\mathbf{u}})+\bar{M})+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}.$

Applying Ascoli’s compactness theorem together with Helly’s theorem, we get the existence of a subsequence ${\Delta t}_{\ell}\to 0$ such that ${\mathbf{g}}_{{\Delta t}_{\ell}}\to{\mathbf{g}}$ in $C([0,T],L^{1}_{\mathrm{loc}}(\Omega;\lambda))$ for some function ${\mathbf{g}}$ . The strict monotonicity of $f^{k}$ implies that

{\mathbf{u}}_{\Delta t}(x,k,t)=(f^{k})^{-1}\big{(}{\mathbf{g}}_{\Delta t}(x,k,t)\big{)},

and hence, also ${\mathbf{u}}_{{\Delta t}_{\ell}}$ converges in $C([0,T],L^{1}_{\mathrm{loc}}(\Omega;\lambda))$ to some function ${\mathbf{u}}$ . Theorem 5.1 implies that ${\mathbf{u}}$ is the entropy solution; since this solution is unique (Theorem 3.1), the entire sequence $\{{\mathbf{u}}_{\Delta t}\}_{{\Delta t}>0}$ must converge to ${\mathbf{u}}$ . ∎

6 Discrete Stationary Solutions for Monotone Flux Functions

So far we have shown that if a sufficiently large class of stationary and discrete stationary solutions exists, then our equations on the network are well posed and the finite volume numerical approximations converge to the entropy solution. In this section we show that such classes exist in the case where either all fluxes $f^{k}$ are strictly increasing or all are strictly decreasing. We henceforth assume that all fluxes are increasing; one can attain analogous results for decreasing fluxes following the same arguments. In the following we want to investigate the sets of discrete stationary solutions implied by the upwind method.

We define

f_{\mathrm{in}}(u)\coloneqq\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}(u),\qquad f_{\mathrm{out}}(u)\coloneqq\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}(u)\qquad\text{for }u\in\mathbb{R}.

It is clear that $f_{\mathrm{in}},f_{\mathrm{out}}$ are monotone by the monotonicity of their summand components. In particular, the two functions are invertible.

For the upwind method the conditions (4.9a) and (4.9b) become


$\displaystyle f^{k}(c^{k})$	$\displaystyle=f^{k}(c^{k})\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},$	(6.1a)
$\displaystyle f^{k}(c^{0})$	$\displaystyle=f^{k}(c^{k})\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}.$	(6.1b)

This is equivalent to


$\displaystyle c^{k}$	$\displaystyle=c^{k}\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{in}}},$	(6.2a)
$\displaystyle c^{0}$	$\displaystyle=c^{k}\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}},$	(6.2b)

due to the invertibility of the flux functions $f^{k}$ . It is obvious as well, that for two different discrete stationary solutions ${\mathbf{c}},{\mathbf{d}}$ satisfying $c^{k}\leqslant d^{k}$ for $k\in{\mathscr{I}}$ , we also have $c^{0}\leqslant d^{0}$ . Henceforth, we denote

\mathscr{G}_{\mathrm{disc}}^{0}\coloneqq\bigl{\{}\text{all discrete stationary states for the upwind method}\bigr{\}}

and we let

\mathscr{G}_{\mathrm{disc}}\coloneqq\bigl{\{}\big{(}c^{-N_{\mathrm{in}}},\dots,c^{-1},c^{1},\dots,c^{N_{\mathrm{out}}}\big{)}\mid{\mathbf{c}}\in\mathscr{G}_{\mathrm{disc}}^{0})\bigr{\}}.

Although it might be too difficult to find a full characterization of the set $L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ of admissible initial data, we will be able to characterize large subsets of $L^{\infty}_{\mathrm{oco}}(\mathscr{G})$ . Let

I_{\mathrm{in}}\coloneqq f_{\mathrm{in}}^{-1}(R_{\mathrm{in}}\cap R_{\mathrm{out}}),\qquad I_{\mathrm{out}}\coloneqq f_{\mathrm{out}}^{-1}(R_{\mathrm{in}}\cap R_{\mathrm{out}})

where

R_{\mathrm{in}}\coloneqq f_{\mathrm{in}}(\mathbb{R}),\qquad R_{\mathrm{out}}\coloneqq f_{\mathrm{out}}(\mathbb{R}).

By the continuity of $f_{\mathrm{in}},f_{\mathrm{out}}$ , the sets $I_{\mathrm{in}},I_{\mathrm{out}}$ are closed intervals.

Theorem 6.1.

We have $\mathscr{L}\subset L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}})$ , where

\mathscr{L}\coloneqq\Bigl{\{}{\mathbf{u}}\in L^{\infty}(\Omega;\lambda)\mid u^{k}(x)\in I_{\mathrm{in}}\ \forall\ k\in{{\mathscr{I}}_{\mathrm{in}}},\ u^{k}(x)\in I_{\mathrm{out}}\ \forall\ k\in{{\mathscr{I}}_{\mathrm{out}}}\Bigr{\}}.

In particular, if $f_{\mathrm{in}},f_{\mathrm{out}}$ have the same range $R_{\mathrm{in}},R_{\mathrm{out}}$ , then $L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}})=L^{\infty}(\Omega;\lambda)$ .

Proof.

Let ${\mathbf{u}}\in\mathscr{L}$ . Since ${\mathbf{u}}\in L^{\infty}(\Omega;\lambda)$ , and $I_{\mathrm{in}},I_{\mathrm{out}}$ are closed, we also have

\underline{c}_{\mathrm{in}}\coloneqq\inf_{\begin{subarray}{c}x\in D^{k}\\ k\in{{\mathscr{I}}_{\mathrm{in}}}\end{subarray}}u^{k}(x)\in I_{\mathrm{in}},\qquad\overline{c}_{\mathrm{in}}\coloneqq\sup_{\begin{subarray}{c}x\in D^{k}\\ k\in{{\mathscr{I}}_{\mathrm{in}}}\end{subarray}}u^{k}(x)\in I_{\mathrm{in}}

and likewise for $\underline{c}_{\mathrm{out}},\overline{c}_{\mathrm{out}}$ . By continuity of $f_{\mathrm{in}},f_{\mathrm{out}}$ , there are $\underline{d}_{\mathrm{in}}\in I_{\mathrm{in}}$ and $\underline{d}_{\mathrm{out}}\in I_{\mathrm{out}}$ satisfying $\underline{d}_{\mathrm{in}}\leqslant\underline{c}_{\mathrm{in}}$ and $\underline{d}_{\mathrm{out}}\leqslant\underline{c}_{\mathrm{out}}$ so that $f_{\mathrm{in}}(\underline{d}_{\mathrm{in}})=f_{\mathrm{out}}(\underline{d}_{\mathrm{out}})$ , that is, the vector $\underline{{\mathbf{d}}}\coloneqq\bigl{(}\underline{d}_{\mathrm{in}},\dots,\underline{d}_{\mathrm{in}},\underline{d}_{\mathrm{out}},\dots,\underline{d}_{\mathrm{out}}\bigr{)}$ is a stationary solution. This stationary solution clearly satisfies (6.2), whence ${\mathbf{d}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ .

In a similar way one finds a stationary solution $\overline{{\mathbf{d}}}\coloneqq\bigl{(}\overline{d}_{\mathrm{in}},\dots,\overline{d}_{\mathrm{in}},\overline{d}_{\mathrm{out}},\dots,\overline{d}_{\mathrm{out}}\bigr{)}\in\mathscr{G}_{\mathrm{disc}}^{0}$ which bounds ${\mathbf{u}}$ from above. Since now

\underline{d}_{\mathrm{in}}\leqslant u^{k}(x)\leqslant\overline{d}_{\mathrm{in}}\;\forall\ k\in{{\mathscr{I}}_{\mathrm{in}}}\qquad\text{and}\qquad\underline{d}_{\mathrm{out}}\leqslant u^{k}(x)\leqslant\overline{d}_{\mathrm{out}}\;\forall\ k\in{{\mathscr{I}}_{\mathrm{out}}}

we conclude that ${\mathbf{u}}\in L^{\infty}_{\mathrm{oco}}(\mathscr{G}_{\mathrm{disc}})$ . ∎

Proposition 6.2.

Consider a conservation law on a network with strictly increasing fluxes $f^{k}$ . Let $\mathscr{G}_{\mathrm{disc}}^{0}$ denote the set of all discrete stationary solutions for the upwind method. Then the set

\mathscr{G}_{\mathrm{disc}}\coloneqq\Big{\{}\big{(}c^{-{N_{\mathrm{in}}}},\dots,c^{-1},c^{1},\dots,c^{{N_{\mathrm{out}}}}\big{)}\mid{\mathbf{c}}\in\mathscr{G}^{0}_{\mathrm{mon}}\Big{\}}

is a mutually consistent, maximal set of stationary solutions.

Proof.

Every ${\mathbf{c}}\in\mathscr{G}_{\mathrm{disc}}$ is a stationary solution due to (4.8).

To prove mutual consistency of $\mathscr{G}_{\mathrm{disc}}$ we plug a discrete stationary solution ${\mathbf{d}}\in\mathscr{G}_{\mathrm{disc}}^{0}$ into (5.1) to get for $n\in\mathbb{N}$ ,

Q_{-{\nicefrac{{1}}{{2}}}}^{k,n}\geqslant Q_{-{\nicefrac{{3}}{{2}}}}^{k,n}\;\text{for }k\in{\mathscr{I}}_{\mathrm{in}}\quad\text{and}\quad Q_{{\nicefrac{{3}}{{2}}}}^{k,n}\geqslant Q_{\nicefrac{{1}}{{2}}}^{k,n}\;\text{for }k\in{\mathscr{I}}_{\mathrm{out}}.

Since we are using the upwind scheme, this reduces to

q_{c^{k}}^{k}(d^{k})\geqslant q_{c^{0}}^{k}(d^{0})\;\text{for }k\in{\mathscr{I}}_{\mathrm{in}}\quad\text{and}\quad q_{c^{0}}^{k}(d^{0})\geqslant q^{k}_{c^{k}}(d^{k})\text{for }k\in{\mathscr{I}}_{\mathrm{out}}.

In the same manner, plugging $d^{0}$ into (5.2) gives us

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}Q_{-{\nicefrac{{1}}{{2}}}}^{k,n}\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}Q_{\nicefrac{{1}}{{2}}}^{k,n}.

Combining these two observations, we get

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}q_{c^{k}}^{k}(d^{k})\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}q_{c^{0}}^{k}(d^{0})\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}q_{c^{0}}^{k}(d^{0})\geqslant\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}q_{c^{k}}^{k}(d^{k}).

As ${\mathbf{c}},{\mathbf{d}}$ were arbitrary, it follows that $\mathscr{G}_{\mathrm{disc}}$ is mutually consistent.

If for some vector ${\mathbf{d}}$ , $\mathscr{G}_{\mathrm{disc}}\cup\{{\mathbf{d}}\}$ is mutually consistent, we know that

\sum_{k\in{\mathscr{I}}_{\mathrm{in}}}q^{k}_{c^{k}}\big{(}d^{k}\big{)}\geqslant\sum_{k\in{\mathscr{I}}_{\mathrm{out}}}q^{k}_{c^{k}}\big{(}d^{k}\big{)}.

We choose $c^{k}=d^{k}\ \forall\ k\in{\mathscr{I}}_{\mathrm{in}}$ and $c^{0}=(f_{\mathrm{out}})^{-1}(\sum_{k\in{\mathscr{I}}_{\mathrm{in}}}f^{k}(c^{k}))$ . Since all $f^{k}$ are monotonically increasing, the entropy flux reduces to $q^{k}_{c^{k}}(d^{k})=\big{|}f^{k}(c^{k})-f^{k}(d^{k})\big{|}$ and thus,

0\geqslant\sum_{k\in{\mathscr{I}}_{\mathrm{out}}}\big{|}f^{k}(c^{0})-f^{k}(d^{k})\big{|},

which implies $d^{k}=c^{0}$ for $k\in{\mathscr{I}}_{\mathrm{out}}$ , and thus, ${\mathbf{d}}\in\mathscr{G}_{\mathrm{disc}}$ . In other words, $\mathscr{G}_{\mathrm{disc}}$ is maximal. ∎

Remark 6.3.

Due to Proposition 6.2 we now know that the doubling of variables argument in Section 3 can be applied in this case and thus, the stability result holds.

Example 6.4.

Consider the traffic flow Example 1.1 with $f^{k}(v)=\alpha^{k}f(v/\alpha^{k})$ for $\alpha^{k}>0$ and $f$ some fixed flux, say, the Burgers flux $f(v)=v^{2}/2$ . Restrict the solutions to the monotone region $u^{k}\geqslant 0$ . Consider the Godunov method with the numerical flux $F(a,b)=f(a)$ . The Rankine–Hugoniot condition reads

\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\frac{(c^{k})^{2}}{\alpha^{k}}=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\frac{(c^{k})^{2}}{\alpha^{k}}.

For the discrete stationary solutions, the additional conditions (4.9) read

\displaystyle f^{k}(c^{0})=f^{k}(c^{k})\qquad\text{for }k\in{{\mathscr{I}}_{\mathrm{out}}}

(the condition for $k\in{{\mathscr{I}}_{\mathrm{in}}}$ automatically holds). Since $f^{k}$ is strictly monotone on $\mathbb{R}_{+}$ , we get $c^{k}=c^{0}$ for all $k\in{{\mathscr{I}}_{\mathrm{out}}}$ . Thus, the set of all discrete stationary solutions is

\mathscr{G}_{\mathrm{disc}}^{0}=\Big{\{}{\mathbf{c}}\in\mathbb{R}^{N+1}\ :\ c^{k}=c^{0}\ \forall\ k\in{{\mathscr{I}}_{\mathrm{out}}},\ {\textstyle\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\frac{(c^{k})^{2}}{\alpha^{k}}=\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\frac{(c^{0})^{2}}{\alpha^{k}}}\Big{\}}.

7 Numerical Examples

We show numerical experiments for some example cases including results for linear and nonlinear as well as convex and concave fluxes. In all experiments we use a CFL number of ${\nicefrac{{1}}{{2}}}$ – that is, ${\Delta t}$ is chosen so that there is equality in (4.7). In all experiments we compute the experimental order of convergence (EOC) as $p\approx\dfrac{\log(e^{j+1}/e^{j})}{\log({\Delta x}_{j+1}/{\Delta x}_{j})}$ on a series of successive grids with $2^{j}$ cells, where $e^{j}$ denotes the $L^{1}$ error on grid level $j$ . The error is computed as the $L^{1}$ difference to a high-resolution reference solution. All errors and EOC are displayed in Table 1.

Figure 2: A network with a periodic edge.

Example 7.1 (Burgers’ equation with roundabout).

In this example we include a roundabout – an edge whose endpoints meet at the same vertex, as shown in Figure 2. We also include an ingoing edge and two outgoing edges, amounting to a total of two ingoing and three outgoing edges. As initial data we choose constants on the roundabout and the outgoing edges and two different constants on the independent ingoing edge. After a while the shock in the initial data on the independent ingoing edge will reach the edge and create new Riemann problems. We choose the initial data

\displaystyle\bar{u}^{-1}(x)=\begin{cases}2&\text{if }0\leqslant x<0.5,\\ \sqrt{2}&\text{if }0.5\leqslant x,\end{cases}\qquad\bar{u}^{-2}=\bar{u}^{1}=\bar{u}^{2}=\bar{u}^{3}\equiv 1.

We take all edges to have length $1$ and choose zero Neumann boundary data on the outer boundaries. On the vertex we set $u_{0}^{0}=1$ . On the ingoing edge with index $-1$ we have a travelling shock wave

u^{-1}(x,t)=\begin{cases}2&\text{if }0\leqslant x<\frac{1}{2-\sqrt{2}}t\\ \sqrt{2}&\text{if }\frac{1}{2-\sqrt{2}}t\leqslant x\end{cases}

which will hit the vertex at $t^{*}=1-\frac{1}{\sqrt{2}}$ . To compute the solution after $t^{*}$ we compute the new vertex value $c^{0}=\sqrt{\nicefrac{{5}}{{3}}}$ and therefore get the Riemann problem

\displaystyle u^{k}(x,t^{*})=\begin{cases}\sqrt{\nicefrac{{5}}{{3}}}&\text{if }x=0,\\ 1&\text{if }0<x\end{cases}

for $k=1,2,3$ , which results in a travelling shock wave with speed $s=\frac{1}{\sqrt{\nicefrac{{5}}{{3}}}-1}$ . At time $t^{**}\coloneqq\sqrt{\nicefrac{{5}}{{3}}}-\frac{1}{\sqrt{2}}$ the travelling shock wave which originated on the roundabout edge hits the vertex once again, resulting in a new set of Riemann problems on the outgoing edge. This process will continue in a periodic fashion.

We compute up to time $T=0.5$ . A plot of the exact and approximate solution to this example at two different times is shown in Figure 3. The accuracy and order of convergence of the numerical approximation are shown in Table 1.

Refer to caption — Figure 3: Initial state and state at $t=0.7$ of a Burgers-type equation with travelling shock wave which hits the vertex at time $t=1-\frac{1}{\sqrt{2}}$ . Here, the graph includes a periodic edge.

Example 7.2.

We construct an example where we take the flux function from the traffic flow example in [13], $f(u)=4u(1-u)$ , but allow for different fluxes on different edges, $f^{k}(u)=\alpha^{k}f\big{(}\frac{u}{\alpha^{k}}\big{)}$ for $\alpha^{k}>0$ , and compute on a star shaped graph with two ingoing edges and three outgoing edges like in Figure 1. The initial data is chosen so that all fluxes are strictly increasing over the range of $\bar{u}$ ; thus, the fluxes $f^{k}$ are in effect monotonously increasing functions. We choose constant solutions on the two ingoing roads and constant initial data on the outgoing roads which are chosen such that on one road a shock will develop, on one road the solution will stay constant over time and on one road a rarefaction wave will develop.

Solving the conditions (4.8), (4.9) for $c^{0}$ yields

c^{0}=\frac{3\pm\sqrt{9-4\bigl{(}\frac{1}{\alpha_{1}}+\frac{1}{\alpha_{2}}+\frac{1}{\alpha_{3}}\bigr{)}\big{(}u^{-1}-\frac{1}{\alpha_{-1}}(u^{-1})^{2}+u^{-2}-\frac{1}{\alpha_{-2}}(u^{-2})^{2}\big{)}}}{2\big{(}\frac{1}{\alpha_{1}}+\frac{1}{\alpha_{2}}+\frac{1}{\alpha_{3}}\big{)}}.

For the incoming edges to have a monotonically increasing flux we impose $u_{-i}\leqslant\frac{1}{2}\alpha_{-i},$ for $i=1,2$ and for outgoing edges $c^{0}\leqslant\frac{1}{2}\min\{\alpha_{1},\alpha_{2},\alpha_{3}\}$ . We choose $\alpha_{-1}=\alpha_{-2}=1$ , $\alpha_{1}=\alpha_{2}=4$ and $\alpha_{3}=2$ with initial data

\displaystyle\bar{u}^{-1}=\bar{u}^{-2}\equiv 0.5,\qquad\bar{u}^{1}\equiv 0,\qquad\bar{u}^{2}\equiv\frac{1}{2}\big{(}3-\sqrt{7}\big{)},\qquad\bar{u}^{3}\equiv 1.

This gives us $c^{0}=\frac{1}{2}(3-\sqrt{7})$ . On the outer boundary we choose zero Neumann boundary conditions. For $u^{3}$ we will get a shock

u^{3}(x,t)=\begin{cases}\frac{1}{2}(3-\sqrt{7})&\text{if }x<st,\\ 0&\text{if }x\geqslant st,\end{cases}

with speed $s=\frac{(3-\sqrt{7})(2-\frac{3-\sqrt{7}}{4})-3}{\frac{3-\sqrt{7}}{2}-1},$ and a rarefaction wave for $u^{1}$ of the form

u^{1}(x,t)=\begin{cases}\frac{3-\sqrt{7}}{2}&\text{if }x<2(\sqrt{7}-1)t,\\ 1-\frac{x}{4t}&\text{if }2(\sqrt{7}-1)t\leqslant x<4t,\\ 0&\text{if }4t\leqslant x.\end{cases}

On edge $2$ we get the constant solution $u^{2}\equiv\frac{1}{2}(3-\sqrt{7})$ .

We compute up to time $T=0.2$ . A plot of the exact and approximate solution at two different timepoints is shown in Figure 4. Accuracy and order of convergence of the numerical approximation are shown in Table 1.

In addition to the examples described above we show errors and experimental order of convergence (EOC) for several additional examples in Table 1.

Example 7.3 (EOC: Linear advection).

We consider a linear advection equation with two ingoing edges and three outgoing edges as in Figure 1 with initial data

\displaystyle\bar{u}^{-1}(x,t)=\begin{cases}2&0\leqslant x<0.8,\\ 1&0.8\leqslant x,\end{cases}\qquad\bar{u}^{-2}\equiv 1,\qquad\bar{u}^{1}=\bar{u}^{2}=\bar{u}^{3}\equiv\frac{2}{3},

and Dirichlet boundary conditions adapted to the edge values. We initialize the vertex node by $u_{0}^{0}=\frac{2}{3}$ . We compute up to time $T=0.5$ .

Example 7.4 (EOC: Burgers’ equation with elementary waves).

We choose $\bar{u}^{-1}=\bar{u}^{-2}\equiv 1$ as initial data on the ingoing roads and $\bar{u}^{1}\equiv 0$ , $\bar{u}^{2}\equiv\sqrt{\nicefrac{{2}}{{3}}}$ and $\bar{u}^{3}\equiv 2$ on the outgoing edges of a star shaped graph as in Figure 1. The conditions on the numerical flux imply then $3(c^{0})^{2}=2\Leftrightarrow c^{0}=\sqrt{\nicefrac{{3}}{{2}}}$ . Thus, we get the following Riemann problems on the outgoing roads:

\displaystyle\bar{u}^{1}(x)=\begin{cases}\sqrt{\nicefrac{{2}}{{3}}}&x=0,\\ 0&x>0,\end{cases}\qquad\bar{u}^{2}=\begin{cases}\sqrt{\nicefrac{{2}}{{3}}}&x=0,\\ \sqrt{\nicefrac{{2}}{{3}}}&x>0,\end{cases}\qquad\bar{u}^{3}=\begin{cases}\sqrt{\nicefrac{{2}}{{3}}}&x=0,\\ 2&x>0,\end{cases}

with zero Neumann boundary conditions at the outer edges. The solution to these problems are a shock, a constant solution and a rarefaction wave, respectively. We compute up to time $T=0.3$ .

Example 7.5 (EOC: Burgers’ equation with travelling shock).

We consider a Burgers-type equation with two ingoing edges and three outgoing edges as in Figure 1 with initial data

\displaystyle\bar{u}^{-1}(x)=\begin{cases}2&0\leqslant x<0.8,\\ 1&0.8\leqslant x,\end{cases}\qquad\bar{u}^{-2}\equiv 1,\qquad\bar{u}^{1}=\bar{u}^{2}=\bar{u}^{3}\equiv\sqrt{\frac{2}{3}},

with Dirichlet boundary conditions of the same value as the associated edge. On the vertex node the initial condition is chosen as $u_{0}^{0}=\sqrt{\frac{2}{3}}$ . We compute up to $T=0.5$ .

	Example 7.3		Example 7.5		Example 7.4		Example 7.1		Example 7.2
Grid level	$L^{1}$ error	EOC	$L^{1}$ error	EOC	$L^{1}$ error	EOC	$L^{1}$ error	EOC	$L^{1}$ error	EOC
3	0.10877	-	0.11630	-	0.14459	-	0.07087	-	0.09904	-
4	0.05496	0.98	0.07136	0.70	0.08016	0.85	0.0546	0.38	0.04913	1.01
5	0.03649	0.59	0.04372	0.71	0.04651	0.79	0.03117	0.81	0.02844	0.79
6	0.02629	0.47	0.02255	0.96	0.02711	0.78	0.01903	0.71	0.01627	0.81
7	0.01830	0.52	0.01360	0.73	0.01495	0.86	0.01115	0.77	0.00919	0.82
8	0.01255	0.54	0.00653	1.06	0.00925	0.69	0.00644	0.79	0.00527	0.80
9	0.00883	0.51	0.00325	1.01	0.00480	0.95	0.00330	0.96	0.00268	0.98
10	0.00625	0.50	0.00160	1.02	0.00295	0.70	0.00173	0.93	0.00150	0.84
11	0.00442	0.50	0.00086	0.90	0.00152	0.96	0.00085	1.03	0.00084	0.84
12	0.00312	0.50	0.00040	1.10	0.00081	0.91	0.00042	1.02	0.00047	0.84

Table 1:

L^{1}

errors and estimated orders of convergence (EOC) for a selection of examples.

7.1 Comments on the experiments

Convergence order estimates for finite volume methods for nonlinear scalar conservation laws are due to Kuznetsov [17] for the continuous flux case and due to Badwaik, Ruf [4] for the case of monotone fluxes with points of discontinuity. In both of those cases the analytically proven convergence rate is at least $\sqrt{{\Delta x}}$ . Our numerical experiments indicate the same lower bound on the convergence rate for our numerical methods on graphs. Considering the fact that $f_{\mathrm{in}}$ and $f_{\mathrm{out}}$ from Section 6 are monotone it might be possible to generalize the result of Badwaik and Ruf to networks.

8 Summary and Outlook

In conclusion we have defined a framework for the analysis and numerical approximation of conservation laws on networks. We extended the concepts well known from the conventional case such as weak solution, entropy solution and monotone methods to make sense on a directed graph. We defined a reasonable entropy condition under which we have shown stability and uniqueness of an analytic solution. Existence is shown by convergence of a conservative, consistent, monotone difference scheme. In an upcoming work [10] we want to address convergence of a numerical method where the fluxes $f^{k}$ are not monotone but concave, as is usually found in traffic flow models. This includes deriving a sufficiently large set of stationary and discrete stationary solutions for this case. Further, we want to extend our model to include boundary conditions and derive a convergence order estimate for numerical approximations. As for future work, a generalization to systems of conservation laws would be highly desirable. One could also try to construct numerical schemes for equations incorporating diffusive fluxes like it was done in [15] on the line. Generalized models would span more complex scenarios such as blood circulation [5] in a network of veins or a river delta by the means of Euler equations and shallow water equations, respectively.

References

[1] B. Andreianov, G. M. Coclite, and C. Donadello. Well-posedness for vanishing viscosity solutions of scalar conservation laws on a network. Discrete and Continuous Dynamical Systems, 37(11):5913–5942, 2017.
[2] B. Andreianov, K. H. Karlsen, and N. H. Risebro. A theory of $L^{1}$ -dissipative solvers for scalar conservation laws with discontinuous flux. Archive for Rational Mechanics and Analysis, 201:27–86, 2011.
[3] E. Audusse and B. Perthame. Uniqueness for scalar conservation laws with discontinuous flux via adapted entropies. Proceedings of the Royal Society of Edingburgh, 135:253–266, 2005.
[4] J. Badwaik and A. M. Ruf. Convergence rates of monotone schemes for conservation laws with discontinuous flux. SIAM Journal on Numerical Analysis, 58(1):607–629, 2020.
[5] A. Bressan, S. Canic, M. Garavello, M. Herty, and B. Piccoli. Flows on networks: recent results and perspectivees. EMS Surveys in Mathematical Sciences, 1:47–111, 2014.
[6] G. M. Coclite and L. di Ruvo. Vanishing viscosity for traffic on networks with degenerate diffusivity. Mediterranean Journal of Mathematics, 16, 2019.
[7] G. M. Coclite and C. Donadello. Vanishing viscosity on a star-shaped graph under general transmission conditions at the node. Networks and heterogeneous media, 15:197–213, 2020.
[8] G. M. Coclite and M. Garavello. Vanishing viscosity for traffic on networks. SIAM Journal on Mathematical Analysis, 42(4):1761–1783, 2010.
[9] M. G. Crandall and L. Tartar. Some relations between nonexpansive and order preserving mappings. Proceedings of the American Mathematical Society, 78(3):385–390, 1980.
[10] U. S. Fjordholm, M. Musch, and N. H. Risebro. Well-posedness of traffic flow models on networks. in preparation, 2021.
[11] M. Garavello. A review of conservation laws on networks. Networks and Heterogeneous Media, 5:565–581, 2010.
[12] M. Garavello and B. Piccoli. Entropy type conditions for riemann solvers at nodes. Advances in differential Equations, 16(1/2):113–144, 2011.
[13] H. Holden and N. H. Risebro. A mathematical model of traffic flow on a network of unidirectional roads. SIAM Journal on Mathematical Analysis, 26:999–1017, 1995.
[14] H. Holden and N. H. Risebro. Front Tracking for Hyperbolic Conservation Laws, volume 52 of Applied Mathematical Sciences. Springer, 2015.
[15] K. H. Karlsen, N. H. Risebro, and J. D. Towers. Upwind difference approximations for degenerate parabolic convection-diffusion equations with a discontinuous coefficient. IMA Journal of Numerical Analysis, 22:623–664, 2002.
[16] S. N. Kruzkov. First order quasilinear equaitons in several independent variables. Math. USSR Sb., 10(2):217–243, 1970.
[17] N. Kuznetsov. Accuracy of some approximate methods for computing the weak solutions of a first-order quasi-linear equation. USSR Computational Mathematics and Mathematical Physics, 16(6):105–119, 1976.
[18] M. J. Lighthill and G. B. Whitham. On kinematic waves II. A theory of traffic flow on long crowded roads. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 229:317–345, 1955.
[19] E. Y. Panov. Existence of strong traces for quasi-solutions of multidimensional conservation laws. Journal of Hyperbolic Differential Equations, 4(4):729–770, 2007.

	$\displaystyle\bigl{\\|}$	$\displaystyle{\mathbf{u}}_{\Delta t}(t^{n+1})-{\mathbf{u}}_{\Delta t}(t^{n})\bigr{\\|}_{L^{1}(\Omega;\lambda)}$
		$\displaystyle=\big{\\|}V({\mathbf{u}}_{\Delta t}(t^{n}))-V({\mathbf{u}}_{\Delta t}(t^{n-1}))\big{\\|}_{L^{1}(\Omega;\lambda)}$
(using Theorem 4.11(c))
		$\displaystyle\leqslant\big{\\|}{\mathbf{u}}_{\Delta t}(t^{n})-{\mathbf{u}}_{\Delta t}(t^{n-1})\big{\\|}_{L^{1}(\Omega;\lambda)}\leqslant\dots\leqslant\big{\\|}{\mathbf{u}}_{\Delta t}(t^{1})-{\mathbf{u}}_{\Delta t}(t^{0})\big{\\|}_{L^{1}(\Omega;\lambda)}$
		$\displaystyle={\Delta x}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\|u_{i}^{k,1}-u_{i}^{0}\|+{\Delta x}_{0}\|u_{0}^{k,1}-u_{0}^{0}\|$
		$\displaystyle={\Delta t}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}F_{i+{\nicefrac{{1}}{{2}}}}^{k,0}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,0}\big{\|}+{\Delta t}\biggl{\|}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{\nicefrac{{1}}{{2}}}^{k,0}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k,0}\biggr{\|}$
		$\displaystyle={\Delta x}\lambda\sum_{\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,0}-u_{i-1}^{k,0}\big{\|}$
		$\displaystyle\quad+{\Delta t}\biggl{\|}\!\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}F_{\nicefrac{{1}}{{2}}}^{k,0}-F^{k,0}\big{(}u_{0}^{0},u_{0}^{0}\big{)}-\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}F_{-{\nicefrac{{1}}{{2}}}}^{k,0}-F^{k,0}(u_{0}^{0},u_{0}^{0})\\ &+\overbrace{\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}f^{k}\big{(}u_{0}^{0}\big{)}}^{\eqqcolon f_{\mathrm{out}}(u_{0}^{0})}-\overbrace{\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}f^{k}\big{(}u_{0}^{0}\big{)}}^{\eqqcolon f_{\mathrm{in}}(u_{0}^{0})}\biggr{\|}\end{aligned}$
		$\displaystyle\leqslant{\Delta t}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}L^{k}\bigl{(}\big{\|}u_{i}^{k,0}-u_{i-1}^{k,0}\big{\|}+\big{\|}u_{i+1}^{k,0}-u_{i}^{k,0}\big{\|}\bigr{)}$
		$\displaystyle\quad+{\Delta t}\underbrace{\Biggl{(}\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}L^{k}\big{\|}u_{1}^{k,0}-u_{0}^{0}\big{\|}+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}L^{k}\big{\|}u_{0}^{0}-u_{-1}^{k,0}\big{\|}+\overbrace{\big{\|}f_{\mathrm{out}}\big{(}u_{0}^{0}\big{)}-f_{\mathrm{in}}\big{(}u_{0}^{0}\big{)}\big{\|}}^{\leqslant M}\Biggr{)}}_{\eqqcolon\bar{M}}$
		$\displaystyle\leqslant{\Delta t}(C\operatorname{TV}({\mathbf{u}}^{0})+\bar{M}),$

	$\displaystyle 0$	$\displaystyle\geqslant-\underbrace{{\Delta t}{\Delta x}\sum_{n=1}^{M}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,n}-c^{k}\big{\|}\biggl{(}\frac{\varphi_{i}^{k,n}-\varphi_{i}^{k,n-1}}{{\Delta t}}\biggr{)}}_{=\,A_{1}}$
		$\displaystyle\quad-\underbrace{{\Delta x}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,0}-c^{k}\big{\|}\varphi_{i}^{k,0}}_{=\,A_{2}}$
		$\displaystyle\quad-\underbrace{{\Delta t}{\Delta x}\sum_{n=0}^{M}\biggl{(}\begin{aligned} &\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i+{\nicefrac{{1}}{{2}}}}^{k,n}\biggl{(}\frac{\varphi_{i+1}^{k,n}-\varphi_{i}^{k,n}}{{\Delta x}}\biggr{)}\\ &+\sum_{k\in{{\mathscr{I}}_{\mathrm{out}}}}\sum_{i\in D_{\mathrm{disc}}^{k}}Q_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\biggl{(}\frac{\varphi_{i}^{k,n}-\varphi_{i-1}^{k,n}}{{\Delta x}}\biggr{)}\biggr{)}\end{aligned}}_{=\,A_{3}}$
		$\displaystyle\quad-\underbrace{{\Delta x}\frac{N}{2}\big{\|}u_{0}^{0}-c^{0}\big{\|}\varphi_{0}^{0}-{\Delta t}{\Delta x}\frac{N}{2}\sum_{n=1}^{M}\big{\|}u_{0}^{n}-c^{0}\big{\|}\biggl{(}\frac{\varphi_{0}^{n}-\varphi_{0}^{n-1}}{{\Delta t}}\biggr{)}}_{=\,A_{4}}.$

	$\displaystyle\operatorname{TV}({\mathbf{g}}_{\Delta t}(\cdot,t))$	$\displaystyle=\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}f^{k}(u_{i}^{k,n})-f^{k}(u_{i-1}^{k,n})\big{\|}$
		$\displaystyle\quad+\sum_{k\in{{\mathscr{I}}_{\mathrm{in}}}}\big{\|}f^{k}(u_{0}^{n})-f^{k}(u_{-1}^{k,n})\big{\|}$
		$\displaystyle\leqslant\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}f^{k}(u_{i}^{k,n})-f^{k}(u_{i-1}^{k,n})\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle=\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}F_{i+{\nicefrac{{1}}{{2}}}}^{k,n}-F_{i-{\nicefrac{{1}}{{2}}}}^{k,n}\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle=\frac{{\Delta x}}{{\Delta t}}\sum_{k\in{\mathscr{I}}}\sum_{i\in D_{\mathrm{disc}}^{k}}\big{\|}u_{i}^{k,n+1}-u_{i}^{k,n}\big{\|}+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}$
		$\displaystyle\leqslant{\Delta x}(C\operatorname{TV}(\bar{\mathbf{u}})+\bar{M})+N\\|{\mathbf{g}}_{\Delta t}\\|_{\infty}.$

Well-posedness theory for nonlinear scalar conservation laws on networks

Abstract

1 Introduction

Example 1.1.

2 Entropy solutions on networks

2.1 Weak solutions

Definition 2.1 (Weak Solution).

Proposition 2.2 (Rankine–Hugoniot condition).

Proof.

Definition 2.3 (Stationary Solution).

Remark 2.4.

2.2 Entropy conditions

Definition 2.5 (Kruzkov entropy pairs).

Definition 2.6.

Definition 2.7.

Example 2.8.

Definition 2.9 (Entropy Solution).

Lemma 2.10.

Proof.

Theorem 2.11.

Proof.

Proposition 2.12.

Proof.

Definition 2.13 (Maximality).

Lemma 2.14.

Proof.

Remark 2.15.

3 Stability and Uniqueness

Theorem 3.1 (Entropy Solutions are L1L^{1} stable).

Proof.

4 Numerical approximation

4.1 A finite volume method on networks

Definition 4.1 (Monotone scheme).

Proposition 4.2.

Proof.

4.2 Discrete stationary solutions

Definition 4.3 (Discrete Stationary Solution).

Notation 4.4.

Remark 4.5.

Remark 4.6.

Lemma 4.7.

Proof.

4.3 A stability framework for general consistent, conservative, monotone methods

Lemma 4.8.

Proof.

Proposition 4.9.

Proof.

Definition 4.10 (L1L^{1} contractive method).

Theorem 4.11 (Crandall–Tartar: [9, Proposition 1]).

Theorem 4.12.

Proof.

Corollary 4.13.

Proof.

5 Convergence of finite volume schemes

Theorem 5.1 (Lax–Wendroff theorem).

Proof.

Theorem 5.2 (Helly’s theorem, [14, p. 437]).

Theorem 5.3 (Compactness and Convergence to a Weak Solution).

Proof.

6 Discrete Stationary Solutions for Monotone Flux Functions

Theorem 6.1.

Proof.

Proposition 6.2.

Proof.

Remark 6.3.

Example 6.4.

7 Numerical Examples

Example 7.1 (Burgers’ equation with roundabout).

Example 7.2.

Example 7.3 (EOC: Linear advection).

Example 7.4 (EOC: Burgers’ equation with elementary waves).

Example 7.5 (EOC: Burgers’ equation with travelling shock).

7.1 Comments on the experiments

8 Summary and Outlook

References

Theorem 3.1 (Entropy Solutions are $L^{1}$ stable).

Definition 4.10 ( $L^{1}$ contractive method).