A differential model for growing sandpiles on networks

TL;DR

This paper models the equilibrium of granular material on networks using a differential system of Monge-Kantorovich type, establishing existence, uniqueness, and numerical methods for solutions.

Contribution

It introduces a differential model for sandpiles on networks and proves well-posedness using viscosity solutions, with numerical approximation and experiments.

Findings

Existence and uniqueness of solutions are proven.

Numerical schemes for the model are developed.

Numerical experiments validate the approach.

Abstract

We consider a system of differential equations of Monge-Kantorovich type which describes the equilibrium configurations of granular material poured by a constant source on a network. Relying on the definition of viscosity solution for Hamilton-Jacobi equations on networks, recently introduced by P.-L. Lions and P. E. Souganidis, we prove existence and uniqueness of the solution of the system and we discuss its numerical approximation. Some numerical experiments are carried out.

Equations

$$\begin{array}[]{l}-D\left(v(x)\,\,Du(x)\right)=f(x),\\[6.0pt] \ (|Du(x)|-1)v(x)=0,\\[6.0pt] \ |Du(x)|\leq 1,\end{array}$$

$$|Du|\le 1$$

$$\begin{array}[]{ll}-D\left(v(x)\eta(x)\,\,Du(x)\right)=f(x),\\[4.0pt] (\eta(x)|Du(x)|-1)v(x)=0\,,\\[4.0pt] \eta(x)|Du(x)|\leq 1\end{array}$$

$${\pi }_j(0),{\pi }_j({\ell }_j)\in \mathcal{V},e_j={\pi }_j((0,{\ell }_j))and{\overline{e}}_j={\pi }_j([0,{\ell }_j]).$$

$$\text{Dist}(x,y):=\underset{\mathcal{P}(x,y)}{\inf }\{|{\pi }_{j_1}^{-1}(x)-{\pi }_{j_1}^{-1}(x_1)|+\sum _{i=2}^n{\ell }_{j_i}+|{\pi }_{j_{n+1}}^{-1}(y)-{\pi }_{j_{n+1}}^{-1}(x_n)|\}.$$

$$u_j(t):=u({\pi }_j(t))\text{and}{u}_j(t):=u_{{\bar{e}}_j}({\pi }_j(t)),t\in [0,{\ell }_j],j\in \mathcal{J}.$$

$${\int }_{\mathcal{N}}u(x)dx=\sum _{j\in \mathcal{J}}{\int }_{e_j}u(x)dx=\sum _{j\in \mathcal{J}}{\int }_0^{{\ell }_j}{u}_j(t)dt,$$

$$D_{j}u(x)=u_j^{\prime }({\pi }_j^{-1}(x)),\text{ if }x\in e_j,$$

$$D_{j}u(x_{i})=\left\{\begin{array}[]{ll}\lim\limits_{h\to 0^{+}}(u_{j}(h)-u_{j}(0))/h\,,&\hbox{if }x_{i}=\pi_{j}(0)\\[4.0pt] \lim\limits_{h\to 0^{+}}(u_{j}(\ell_{j}-h)-u_{j}(\ell_{j}))/h\,,&\hbox{if }x_{i}=\pi_{j}(\ell_{j})\,.\end{array}\right.$$

$$({\eta }_j{)}_{j\in \mathcal{J}}\in \prod _{j=1}^{M}C([0,{\ell }_j])\text{ and }\min \{{\eta }_j(x);x\in {\bar{e}}_j,j\in \mathcal{J}\}>0.$$

$${\cal D}(x,y):=\inf_{{\cal P}(x,y)}\left\{\big{|}\int_{\pi_{j_{1}}^{-1}(x)}^{t_{1}}\frac{1}{\eta_{j_{1}}(s)}ds\big{|}\ +\sum_{i=2}^{n}\big{|}\int_{t_{i-1}}^{t_{i}}\frac{1}{\eta_{j_{i}}(s)}ds\big{|}+\big{|}\int_{t_{n}}^{\pi_{j_{n+1}}^{-1}(y)}\frac{1}{\eta_{j_{n+1}}(s)}ds\big{|}\right\}$$

$$d(x):=\underset{y\in \partial \mathcal{N}}{\min }\mathcal{D}(x,y),x\in \mathcal{N}.$$

$$|Du(x)|-\frac{1}{\eta (x)}=0,x\in \stackrel{\circ }{N},$$

$$|D_j\varphi (x)|-\frac{1}{{\eta }_j(x)}\le 0;$$

$$|D_j\varphi (x)|-\frac{1}{{\eta }_j(x)}\ge 0;$$

$$\underset{j\in {\mathrm{Inc}}_i}{\max }\left\{|D_j\varphi (x_i)|-\frac{1}{{\eta }_j(x_i)}\right\}\ge 0;$$

$$S_j(u):=\{t\in (0,{\ell }_j):u_j\text{ is not differentiable at }t\},j\in \mathcal{J},$$

$${\sigma }_{ij}(u):=\mathrm{sgn}[D_{j}u(x_i)],i\in \mathcal{I},j\in {\mathrm{Inc}}_i,$$

$${\mathrm{Inc}}_i^{\pm }(u):=\{j\in {\mathrm{Inc}}_i:{\sigma }_{ij}(u)=\pm 1\}\text{and}{N}_i^{\pm }(u)=\#{\mathrm{Inc}}_i^{\pm }(u),i\in \mathcal{I}.$$

$${\sigma }_{ij}(d)+{\sigma }_{ik}(d)=0.$$

$$S(d):=\{{\pi }_j(S_j(d));j\in \mathcal{J}\}\cup \{x_i\in \mathcal{V}:d\text{ has a local maximum at }{x}_i\}.$$

$$\#\mathcal{I}-\#\mathcal{J}\le \#S(d)\le \#\mathcal{J}-\#{\mathcal{I}}_T.$$

$${\Sigma }_j(d):=S_j(d)\cup T_j(d),j\in \mathcal{J},$$

$$T_j(d):=\{{\pi }_j^{-1}(x_i);x_i\in {\overline{e}}_j\text{ s.t. }j\in {\mathrm{Inc}}_i^{-}(d)\},j\in \mathcal{J}.$$

$$P_j(t):=t+{\tau }_j(t){\eta }_j(t)d_j^{\prime }(t),t\in [0,{\ell }_j],$$

$${\tau }_j(t):=\min \{s\ge 0:t+s{\eta }_j(t)d_j^{\prime }(t)\in {\Sigma }_j(d)\}$$

$$d(y)=\mathcal{D}(x,y)+d(x).$$

$$f\in L^{\infty }(\mathcal{N}),f\ge 0,\text{meas}(\mathrm{supp}(f))>0,$$

$$-$$

$$\eta (x)|Du(x)|\le 1,$$

$$\eta (x)|Du(x)|=1,$$

Full text

** AMS subject classification:**

35R02, 35C15, 47J20, 49L25, 35M33.

Keywords:

networks, granular matter, Monge-Kantorovich system, viscosity solutions.

1 Introduction

In this paper we shall analyze the differential system

[TABLE]

on a network ${\cal N}$, i.e. a collection of vertices joined by non self-intersecting edges, and we shall provide a characterization of the solution $(u,v)$.

System (1.1), when considered on a bounded domain $\Omega$ of $\mathbb{R}^{n}$, arises in several different frameworks. For example, it characterizes optimal plans in the mass transfer problem, and it is also related to the behavior of the solution of the $p$-Laplace equation as $p\to\infty$. Another illustrating example for system (1.1) is a mathematical model for granular matter. The deposition of homogeneous granular matter such as sand, when being poured onto objects from sources above, is well understood ([4, 9, 12]). In this framework system (1.1) is related to equilibrium configurations which occur on flat and bounded domains without sides (such as tables). In these configurations the granular matter can only form heaps with local steepness not exceeding the angle of repose $\alpha$ characteristic of the matter. Consequently, if we denote by $u$ the height of the standing layer at $x$ and normalize the angle $\alpha$ in such a way that $\tan(\alpha)=1$, then $u$ must satisfy the equation

[TABLE]

inside the domain and vanish on the boundary. A selection criterium among all the possible admissible configurations is given by the maximal volume solution of (1.2). This extremal solution coincides with the distance function from the boundary and it is characterized as the unique viscosity solution of (1.2) satisfying $u=0$ at the boundary. Eventually matter poured by the source on the standing layer $u$ would roll at a speed proportional to the slope $Du$. Denoted with $v$ the height of the rolling layer, since the matter is conserved inside the domain, it follows that $v$ satisfies a conservation law given by first equation in (1.1). The analysis of (1.1) performed in the one-dimensional case in [9] has been extended to bi-dimensional domains in [5] (see also [6, 7, 8]) where existence, uniqueness and representation formula for the solution of system (1.1) completed with the appropriate boundary conditions are proved.

In this paper, we are interested in the case where the matter is poured on a network. More precisely, we suppose that the edges of the network are bounded on both sides by sufficiently high walls and that the sand can run out of the network only at the boundary vertices, while at the other vertices is interchanged between the incident edges. Pouring sand in a network, several sand heaps will start to grow, each two of them separated by at least one boundary vertex. At the equilibrium any additional sand portion which violates the angle of repose is forced to leave the network at the boundary vertices. The profile of the standing layer $u$ is then described by a continuous function on the network, which vanishes at the boundary points, maximizes the volume functional, and satisfies the eikonal equation almost everywhere in the edges. As in the case of a bounded domain $\Omega$, the distance function from the boundary satisfies all these properties and it can be characterized as the viscosity solution of the corresponding eikonal equation with vanishing boundary condition. Assuming that inside the edges the matter is conserved, the height of the rolling layer $v$ satisfies a conservation law as in (1.1). Hence we have the following system on the network

[TABLE]

where $f$ is the source term and $\eta$ represents a non constant angle of repose. It is important to observe that the notion of viscosity solution for the Hamilton-Jacobi equations involves also the internal vertices (see Definition 3.1). Hence the eikonal equation in (1.3) has to be completed only with a condition at the boundary vertices, i.e. $u=0$ on $\partial{\cal N}$. Instead, we need to add to the conservation law in (1.3) transmission conditions at the internal vertices, expressing the conservation of the total matter interchanged among the edges incident a same vertex (see (4.7)).

A preliminary and fundamental step for the analysis of (1.3) requires the study of the regularity of the distance function from the network boundary and the characterization of its singular set. Then, for the system (1.3) completed with the mentioned appropriate boundary and transmissions conditions, we prove existence through a representation formula, and uniqueness (to be intended for $v$ on the whole network and for $u$ on the set $\{v\neq 0\}$, exactly as in the case of a bounded domain $\Omega$).

The paper is organized as follows. In Section 2 we introduce some basic notations and definitions. Section 3 is devoted to the study of the eikonal equation on the network. In Section 4 we prove existence of a solution of (1.3) on the network through its representation formula, while in Section 5 we give the uniqueness result. A finite-difference approximation scheme for the problem and some numerical examples illustrating the theory are described in Sections 6 and 7 respectively.

2 Preliminary definitions and notations

This section is devoted to the definitions and notations that we shall use in the sequel. These definitions are nowadays classical but not necessarily standard in the literature and we recall them for the readers convenience. A network ${\cal N}$ is a finite collection of $N$ distinct points $x_{i}$ in $\mathbb{R}^{n}$, also denoted vertices or nodes, and $M$ non self-intersecting distinct curves $e_{j}$ in $\mathbb{R}^{n}$, also denoted edges or arcs, whose endpoints are vertices of ${\cal N}$. Setting ${\cal V}:=\{x_{i}\}_{i\in\mathcal{I}}$, $\mathcal{I}=\{1,\dots,N\}$, and ${\cal E}:=\{e_{j}\}_{j\in\mathcal{J}}$, $\mathcal{J}=\{1,\dots,M\}$, we have ${\cal N}=({\cal V},{\cal E})$. We also assume that ${\cal N}$ is connected and that there are no loops.

To each $e_{j}\in{\cal E}$ is associated a diffeomorphism $\pi_{j}:[0,\ell_{j}]\to\mathbb{R}^{n}$, $\ell_{j}\in(0,\infty)$, such that

[TABLE]

For $x_{i}\in{\cal V}$, ${\rm Inc}_{i}:=\{j\in\mathcal{J}:x_{i}\in\overline{e}_{j}\}$ is the set of all the indices of the edges having an endpoint at $x_{i}$, and when $j,k\in{\rm Inc}_{i}$, we shall say that $e_{j}$ and $e_{k}$ are incidents. We define the boundary of the network as a subset $\partial{\cal N}:=\{x_{i}\}_{i\in\mathcal{I}_{B}}$ of ${\cal V}$, the boundary vertices, for a given nonempty set $\mathcal{I}_{B}\subset\mathcal{I}$, and we set ${\overset{\circ}{\cal N}}:={\cal N}\setminus\partial{\cal N}$. We always assume $x_{i}\in\partial{\cal N}$ whenever $\#{\rm Inc}_{i}=1$. The complement set $\{x_{i}\}_{i\in\mathcal{I}_{T}}$, $\mathcal{I}_{T}:=\mathcal{I}\setminus\mathcal{I}_{B}$ , shall be called the set of transition vertices.

The $\pi_{j}$’s allows us to endow the network with the following natural metric (see [3, 10]). Given $x,y\in{\cal N}$, let ${\cal P}(x,y)$ denote a path connecting them along ${\cal N}$, i.e. a finite sequence of closed edges and sub-edges $(\bar{e}_{j_{1}}\cap\bar{e}(x),\bar{e}_{j_{2}},\dots,\bar{e}_{j_{n}},\bar{e}_{j_{n+1}}\cap\bar{e}(y))$ such that $x\in\bar{e}_{j_{1}}$, $y\in\bar{e}_{j_{n+1}}$, $\bar{e}_{j_{i}}\cap\bar{e}_{j_{i+1}}=\{x_{i}\}\subset{\cal V}$, $i=1,\dots,n$, and $e(x)$ and $e(y)$ are the sub-edges with endpoints $x$, $x_{1}$ and $y$, $x_{n}$ respectively. Then,

[TABLE]

The distance (2.1) makes ${\cal N}$ a compact metric network.

To each function $u$ defined on ${\cal N}$ and each function $u=(u_{\bar{e}_{j}})_{j\in\mathcal{J}}$ defined on $\prod_{j=1}^{M}\bar{e}_{j}$, with $u_{\bar{e}_{j}}$ defined on $\bar{e}_{j}$, we associate the projection $(u_{j})_{j\in\mathcal{J}}$ defined on the parameters’ space as

[TABLE]

(2.2) allows us to make no difference between $u:{\cal N}\to\mathbb{R}$ and $u=(u_{\bar{e}_{j}})_{j\in\mathcal{J}}$ in the sequel. For the sake of simplicity, we shall also improperly write $u_{j}(x)$ whenever $x\in\bar{e}_{j}$, instead of $u_{j}(\pi_{j}^{-1}(x))$.

Next, the integral of $u$ on ${\cal N}$ is naturally defined as

[TABLE]

while the Lebesgue spaces $L^{p}({\cal N})$, $p\in[1,\infty]$, are defined as the direct product spaces $\prod_{j=1}^{M}L^{p}(0,\ell_{j})$ endowed with the norm $\|u\|_{p}:=\sum_{j\in\mathcal{J}}\|u_{j}\|_{L^{p}(0,\ell_{j})}$. The space $C({\cal N})$ of continuous functions on ${\cal N}$ is the space of $u$ such that $(u_{j})_{j\in\mathcal{J}}\in\prod_{j=1}^{M}C([0,\ell_{j}])$ and $u_{j}(x_{i})=u_{k}(x_{i})$ for all $j,k\in{\rm Inc}_{i}$ and all $i\in\mathcal{I}$.

Concerning derivatives, with $Du(x)$, $x\in{\cal N}$, we denote $(D_{j}u(x))_{j\in\mathcal{J}}$, where

[TABLE]

while, if $x=x_{i}\in{\cal V}$ and $j\in{\rm Inc}_{i}$, $D_{j}u(x_{i})$ is the internal oriented derivative of $u$ at $x_{i}$ along the arc $e_{j}$, i.e.

[TABLE]

Then, the space $C^{1}({\cal N})$ consists of all the functions $u\in C({\cal N})$ such that $(u_{j})_{j\in\mathcal{J}}\in\prod_{j=1}^{M}C^{1}([0,\ell_{j}])$ and it is endowed with the norm $\|u\|_{C^{1}}=\max_{k=0,1}\|D^{k}u\|_{\infty}$. Observe that no continuity condition at the vertices is prescribed for the derivatives of a function $u\in C^{1}({\cal N})$.

Finally, as for the Lebesgue spaces $L^{p}({\cal N})$, the Sobolev spaces $W^{1,p}({\cal N})$, $p\in[1,\infty]$, is the product space $\prod_{j=1}^{M}W^{1,p}(0,\ell_{j})$ endowed with the norm $\|u\|_{1,p}:=\sum_{j\in\mathcal{J}}\|u_{j}\|_{W^{1,p}(0,\ell_{j})}$. It is worth noticing that, with the above definition, $u\in W^{1,p}({\cal N})$ implies $(u_{j})_{j\in\mathcal{J}}\in\prod_{j=1}^{M}C([0,\ell_{j}])$ but not $u\in C({\cal N})$. This is not the standard definition of the Sobolev space, but it is the convenient one for the problem we are concerned here (see [3]).

To conclude this preliminary section, we observe that the diffeomorphisms $\pi_{j}$ induce necessarily an orientation on the edges $e_{j}$. However, all the results that will follow are independent on that orientation as well as on the $\pi_{j}$ themselves. Indeed, changing the family $(\pi_{j})_{j\in\mathcal{J}}$ leads simply to turn $({\cal N},\text{Dist})$ into an isomorphic compact metric network.

3 A weighted distance on ${\cal N}$ and its singular set

In view of the problem motivating the present analysis, it is natural to assume that the network is not homogeneous. Therefore, we introduce a measure of the capacity of each edge of the network to transport and allocate matter (i.e. a non constant angle of repose) through a function $\eta$ satisfying

[TABLE]

Then, with the notations of the previous section, we define a new metric on ${\cal N}$ taking into account the heterogeneity of the edges, as

[TABLE]

where, for a given path ${\cal P}(x,y)$, we have set $\pi_{j_{i}}^{-1}(x_{i-1})=t_{i-1}$ and $\pi_{j_{i}}^{-1}(x_{i})=t_{i}$, for simplicity. Since the network is finite, the number of paths connecting $x$ to $y$ and composed of distinct edges is finite too. Therefore, the infimum in (3.2) is finite and attained. With (3.2), we also define the usual distance function from the boundary $\partial{\cal N}$

[TABLE]

We shall call all paths realizing ${\cal D}(x,y)$ and $d(x)$, geodesic paths. Obviously, ${\cal D}(x,y)$ is equivalent to the metric (2.1) induced by the $(\pi_{j})_{j\in\mathcal{J}}$. However, the geodesic paths given by (3.2) are not necessary the same ones given by (2.1) (see the numerical tests in Section 7).

The remaining of this section is devoted to prove that (3.3) is a viscosity solution of the eikonal equation

[TABLE]

according to Definition 3.1 below. There are several frameworks in which a viscosity solution theory for Hamilton-Jacobi equations on networks has been developed [1, 10, 13]. Here we shall consider the theory recently introduced in [11], since it allows us to deal with non continuous hamiltonians, as the one in (3.4). This theory has been developed for a flat junction-type network, but extends easily to the network under consideration, yelding that a viscosity solution $u$ of (3.4) is a viscosity solution in each edge $e_{j}$ and a constrained supersolution at $x_{i}$, for all $i\in{\cal I}_{T}$. We shall give also the definitions and the results necessary to obtain the regularity properties of $d$ we shall need in the sequel. All the proofs are postponed to the Appendix A since they are quite classical.

Definition 3.1**.**

Given $u\in C({\cal N})$,

• (i)

$u$* is a (viscosity) subsolution of (3.4) if for any $x\in e_{j}$, $j\in\mathcal{J}$, and for any test function $\phi\in C^{1}({\cal N})$ for which $(u-\phi)$ attains a local maximum at $x$, we have*

[TABLE]

• (ii)

$u$* is a (viscosity) supersolution of (3.4) if the following holds:*

$\bullet$* for $x\in e_{j}$, $j\in\mathcal{J}$ and for any test function $\phi\in C^{1}({\cal N})$ such that $(u-\phi)$ attains a local minimum at $x$, we have*

[TABLE]

$\bullet$* for $x_{i}\in{\cal V}$, $i\in\mathcal{I}_{T}$, and for any test function $\phi\in C^{1}({\cal N})$ such that $(u-\phi)$ attains a local minimum at $x_{i}$, we have*

[TABLE]

• (iii)

$u$* is a (viscosity) solution of (3.4) if it is both a viscosity subsolution and a viscosity supersolution of (3.4).*

Definition 3.2**.**

Given a function $u\in W^{1,\infty}({\cal N})$, we set

[TABLE]

and, whenever $D_{j}u(x_{i})$ exists and is not zero,

[TABLE]

[TABLE]

The set $S_{j}(u)$ is the set of singular points of $u$ inside the edge $e_{j}$, while $\sigma_{ij}(u)$ is the slope of $u$ at the vertex $x_{i}$ along the arc $e_{j}$. We observe that $\sigma_{ij}(u)=1$ (respectively $\sigma_{ij}(u)=-1$) if and only if the graph of $u_{j}$ leaves $x_{i}$ “uphill” (respectively “downhill”) the vertex $x_{i}$. Moreover, $\sigma_{ij}(u)$ does not depend on the orientation of $e_{j}$ induced by $\pi_{j}$.

The next proposition states that if $u$ is a viscosity solution of the eikonal equation (3.4), then each edge contains no or exactly one singular point.

Proposition 3.3**.**

Let $u$ be a viscosity solution of (3.4). Then, $u\in(W^{1,\infty}\cap C)({\cal N})$, satisfies the eikonal equation a.e. over ${\cal N}$ and

• (i)

$u$* does not attain a local minimum on ${\overset{\circ}{\cal N}}$;*

• (ii)

$u$* attains a local maximum at $x\in e_{j}$ if and only if $\pi_{j}^{-1}(x)\in S_{j}(u)$;*

• (iii)

for all $i\in\mathcal{I}$ and $j\in\mathcal{J}$, $\sigma_{ij}(u)$ is well defined and $\#S_{j}(u)\in\{0,1\}$; moreover, if $j\in{\rm Inc}_{i}\cap{\rm Inc}_{k}$, $i,k\in\mathcal{I}$, it holds:

$\bullet$* $\#S_{j}(u)=0$ if and only if $\sigma_{ij}(u)+\sigma_{kj}(u)=0$,*

$\bullet$* $\#S_{j}(u)=1$ if and only if $\sigma_{ij}(u)=\sigma_{kj}(u)=1$ ;*

• (iv)

$\sum_{i\in\mathcal{I}}N_{i}^{\pm}(u)\mp\sum_{j\in\mathcal{J}}\#S_{j}(u)=\#\mathcal{J}$* .*

Proposition 3.4**.**

For any fixed $x_{0}\in{\cal N}$, the function ${\cal D}(x_{0},\cdot)$ is a viscosity solution of (3.4) in ${\overset{\circ}{\cal N}}\setminus\{x_{0}\}$ and $d$ is the unique viscosity solution of (3.4) with $d=0$ on $\partial{\cal N}$.

We are now in a position to provide a complete description of the singular set of the distance function $d$. It is worth to recall that in the case of a smooth domain $\Omega\subset\mathbb{R}^{n}$, the singular set of the euclidian distance from $\partial\Omega$ is the set of points where this function is not differentiable. Its closure coincides with the set of points having multiple geodesics connecting them to $\partial\Omega$.

In the case of a network, the structure of the singular set is determined as well by the structure of the network (see Proposition 3.3 (iv) and (3.13) below). However, the two characterizations do not apply and do not coincide, in general, as they are: there could be points of the network connected to the boundary by more than one geodesic path and where the distance from boundary is differentiable. Indeed, as proved above, if $d$ is not differentiable at $x\in{\cal E}$, then $x$ is a local maximum point for $d$ and there are at least two geodesic paths connecting $x$ to $\partial{\cal N}$. On the other hand, if $x_{i}\in{\cal V}$, with $\#{\rm Inc}_{i}>2$, is a transition vertex connected to $\partial{\cal N}$ by two distinct geodesic paths, and $x_{i}$ is not a maximum point of $d$, then there exists at least one edge incident to $x_{i}$ such that (all or some of) the points belonging to it can be also connected to $\partial{\cal N}$ by the same geodesic paths. At the same time the differentiability of $d$ on that points is not a priori excluded. It is also easy to observe that if $x_{i}\in{\cal V}$, $i\in\mathcal{I}_{T}$, is not a maximum point of $d$ on ${\cal N}$, then there exist $j,k\in\mathcal{J}$ such that $j\in{\rm Inc}_{i}^{+}(d)$, $k\in{\rm Inc}_{i}^{-}(d)$ and so

[TABLE]

Note that by Proposition 3.3, $\sigma_{ij}(d)$ and hence ${\rm Inc}_{i}^{\pm}(d)$ are defined for all $i\in\mathcal{I}$ and $j\in\mathcal{J}$. Furthermore, conditions (3.11) can be interpreted as a weak differentiability (classical if $\eta_{j}(x_{i})=\eta_{k}(x_{i})$) of $d$ at $x_{i}$ along the couple $(\bar{e}_{j},\bar{e}_{k})$.

In light of this observations, the natural definition of the singular set of $d$ to the case of a network, conciliating the two characterizations above, is the following

[TABLE]

Using claim (iv) of Proposition 3.3 and the facts that $N_{i}^{-}(d)=0$ while $N_{i}^{+}(d)\geq 1$ if $i\in\mathcal{I}_{B}$, $N_{i}^{-}(d)\geq 2$ while $N_{i}^{+}(d)=0$ if $x_{i}$ is a transition vertex where $d$ has a local maximum, and $N_{i}^{\pm}(d)\geq 1$ on the remaining vertices, it is easily seen that

[TABLE]

Finally, we shall define a normal distance to $S(d)$, selecting on each $\bar{e}_{j}$ the point nearest to $S(d)$. So, let introduce the projection set in $\bar{e}_{j}$

[TABLE]

where

[TABLE]

By Proposition 3.3 again, given $j\in\mathcal{J}$, one set between $S_{j}(d)$ and $T_{j}(d)$ is a singleton and the other one is empty. Therefore, $\Sigma_{j}(d)$ is also a singleton and, for all $j\in\mathcal{J}$, we can define the projection of $t\in[0,\ell_{j}]$ onto the projection set $\Sigma_{j}(d)$ as

[TABLE]

where

[TABLE]

measure the distance of $\pi_{j}(t)$ to $\pi_{j}(P_{j}(t))$ along the direction of $D_{j}d(\pi_{j}(t))$ in the metric (2.1). We have that $(\tau_{j})_{j\in\mathcal{J}}\in\Pi_{j=1}^{M}C([0,\ell_{j}])$ and $(\tau_{j})_{j\in\mathcal{J}}$ is zero on $\Pi_{j=1}^{M}\Sigma_{j}(d)$. However, it is not possible to define from $(\tau_{j})_{j\in\mathcal{J}}$ a continuous function on ${\cal N}$. This is one of the major differences with the normal distance to the cut locus defined in [5, 6]. On the other hand, it is easy to see that for any $j\in\mathcal{J}$ and $x\in\bar{e}_{j}$, setting $t=\pi_{j}^{-1}(x)$ and $y=\pi_{j}(t+r\,\eta_{j}(t)\,d^{\prime}_{j}(t))\in\bar{e}_{j}$, with $r\in[0,\tau_{j}(t)]$, it holds

[TABLE]

Furthermore, one can readily iterate the previous projection procedure to prove that for all $x\in{\cal N}$ there exists at least one $y\in S(d)$ such that (3.16) holds true. As a consequence, any geodesic path realizing $d(x)$ does not contain points of the singular set $S(d)$ except possibly $x$, $d$ is piecewise $C^{1}$ along the geodesic paths, and the only points of non differentiability are transition vertices of the path, where $d$ satisfies (3.11).

4 The sandpiles problem on ${\cal N}$ : existence of a solution

Let $\eta$ be given as in (3.1), and let the (constant in time) matter source be represented by a function $f$ satisfying

[TABLE]

where ${\rm supp}(f)$ stands for the usual essential support. In addition, to every subset $A_{i}$ of ${\rm Inc}_{i}$, $i\in\mathcal{I}_{T}$, we associate a fixed collection of positive coefficients $(C_{ij})_{j\in A_{i}}$ such that $\sum_{j\in A_{i}}C_{ij}=1$.

We are now in a position to consider system (1.3) on the network ${\cal N}$, i.e.

[TABLE]

endowed with the Dirichlet homogeneous boundary condition at the boundary vertices

[TABLE]

to complete (4.3) and (4.4), and the transmission conditions at each $x_{i}$, $i\in\mathcal{I}_{T}$,

[TABLE]

to complete the conservation law (4.2).

The $(C_{ij})_{j\in{\rm Inc}_{i}^{-}(u)}$ in (4.7) are the positive coefficients associated to the subset ${\rm Inc}_{i}^{-}(u)$ of ${\rm Inc}_{i}$, and condition (4.7) amounts to impose that the mass $\sum_{k\in{\rm Inc}_{i}^{+}(u)}v_{k}(x_{i})$ of the rolling layer $v$ entering in a given transition vertex $x_{i}$ is released in each of the $N_{i}^{-}(u)$ outgoing edges according to the distribution coefficients $C_{ij}$. It is worth noticing here that it is possible to consider $C_{ij}=1/{N_{i}^{-}(u)}$ for all $j\in{\rm Inc}_{i}^{-}(u)$, which corresponds to assume that all the mass entering in a vertex is uniformly distributed in the outgoing arcs. Indeed, a remarquable consequence of the uniqueness result we shall prove is that the $N_{i}^{-}(u)$ are invariant for all $i\in\mathcal{I}_{T}$ such that $v(x_{i})\neq 0$ and depend uniquely on the structure of the network, since here ${\rm Inc}_{i}^{-}(u)={\rm Inc}_{i}^{-}(d)$ (see Corollary 5.4). Furthermore, since $\eta$ is not a priori continuous on the vertices, it has to be expected that the $v$ component of the solution is also not continuous on the vertices. This is the reason why (4.4) has to be solved on the set of $x\in{\overset{\circ}{\cal N}}$ where $v(x)\neq 0$ instead of $v(x)>0$, thus including the transition vertices $x_{i}$ where $v_{j}(x_{i})$ is zero for some of the $j\in{\rm Inc}_{i}$. A natural definition of weak solution of (4.2)–(4.7) is then the following.

Definition 4.1**.**

We say that $(u,v)$ is a solution of (4.2)–(4.7) if

• (i)

$(v_{j})_{j\in\mathcal{J}}\in\prod_{j=1}^{M}C([0,\ell_{j}])$* and $v\geq 0$ in ${\cal N}$;*

• (ii)

$u\in(W^{1,\infty}\cap C)({\cal N})$, $u\geq 0$ in ${\cal N}$, $\eta(x)\,|Du(x)|\leq 1$ a.e. in ${\overset{\circ}{\cal N}}$ and $u$ is a viscosity solution of

[TABLE]

• (iii)

For every function $\psi\in(W^{1,\infty}\cap C)({\cal N})$ such that $\psi=0$ on $\partial{\cal N}$, it holds

[TABLE]

• (iv)

$u$* and $(u,v)$ satisfy the boundary and the transition conditions (4.6) and (4.7), respectively .*

Before proving the existence result for (4.2)–(4.7), it is useful to analyse the consistency and well-posedness of the transition condition (4.7). First of all, it is worth noticing that if $(u,v)$ is a classical solution of (4.2), then $(u,v)$ satisfies (4.8) and the conservation of the flux at each transition vertices $x_{i}$, i.e.

[TABLE]

and vice-versa. However, (4.9) is not sufficient to make the problem well posed since the values $v_{j}(x_{i})$ for each $j\in{\rm Inc}_{i}$ are not univocally determined by (4.9), and more specific conditions has to be considered. Next, if $(u,v)$ is a solution in the sense of the definition above, and if $j\in{\rm Inc}_{i}$, $i\in\mathcal{I}_{T}$, is such that $\sigma_{ij}(u)$ is not defined, it holds necessarily $v_{j}(x_{i})=0$. Indeed, assuming by contradiction that $v_{j}(x_{i})>0$, the continuity of $v_{j}$ on $[0,\ell_{j}]$ implies the existence of a sub-interval of $[0,\ell_{j}]$, with one endpoint in $\pi_{j}^{-1}(x_{i})$, along which $v_{j}>0$. Hence, $u$ is a viscosity solution of the eikonal equation on that sub-interval and Proposition 3.3 assures the existence of $\sigma_{ij}(u)$. On the other hand, if $\sigma_{ij}(u)$ is well defined, it determines if $j$ is either in ${\rm Inc}_{i}^{+}(u)$ or in ${\rm Inc}_{i}^{-}(u)$. If eventually ${\rm Inc}_{i}^{+}(u)=\emptyset$, we use the classical convention that the sum in the r.h.s. of (4.7) is zero. In any case, the transition condition (4.7) is meaningful and implies

[TABLE]

i.e., by Proposition 3.3 again, the conservation of the flux (4.9) at each transition vertices.

We now prove the existence result giving an explicit representation formula for a solution of the problem. This formula generalizes the one in [9] (see also [5, 6]) and at the same time takes into account the transmission of the matter through the transition vertex.

Theorem 4.2**.**

A solution of (4.2)–(4.7) in the sense of Definition 4.1, is given by the pair $(d,v^{f})$, with $d$ the distance function defined in (3.3) and the projections $v_{j}^{f}$, $j\in\mathcal{J}$, of $v^{f}$ given by

[TABLE]

where $\chi_{T_{j}(d)}$ denotes the characteristic function of $T_{j}(d)$ and $x_{i}=\pi_{j}(P_{j}(t))$. Moreover, $v^{f}$ is zero on the singular set $S(d)$ defined in (3.12), $v^{f}\in W^{1,\infty}({\cal N})$ and $(d,v^{f})$ satisfies (4.2) pointwise on ${\cal E}\setminus\{\pi_{j}(S_{j}(d));j\in\mathcal{J}\}$.

Let us observe that the first term in the r.h.s. of (4.11) is non negative and takes into account the matter poured by the source vertically onto each edge (see [5, 9]). Concerning the second term, for a fixed $j\in\mathcal{J}$, if $T_{j}(d)$ is empty, or equivalently if $S_{j}(d)$ is a singleton, $\chi_{T_{j}(d)}\equiv 0$ and the term makes sense giving no contribution. If $T_{j}(d)$ is not empty, then $T_{j}(d)=\{P_{j}(t)\}\subset\{0,\ell_{j}\}$ for all $t\in[0,\ell_{j}]$, $S_{j}(d)$ is empty, $d$ is strictly monotone on $e_{j}$ and $\pi_{j}(P_{j}(t))$ is the endpoint $x_{i}$ of $e_{j}$ such that $\sigma_{ij}(d)=-1$. Again, the second term makes sense and it gives a positive contribution to $v^{f}_{j}(t)$ iff $x_{i}$ is not a maximum point for $d$ over ${\cal N}$ and therefore ${\rm Inc}_{i}^{+}(d)$ is not empty (i.e. if there are edges $e_{k}$ “ingoing downhill” into $x_{i}$). Resuming, the second term in the r.h.s. of (4.11) adds to the rolling layer due to the source $f$, the rolling layer coming from the ingoing edges.

Proof.

Thanks to Propositions 3.3 and 3.4, we need only to prove that $(d,v^{f})$ satisfies conditions (i) and (iii) in Definition 4.1, the transition condition (4.7), and to check that $(v^{f}_{j})_{j\in\mathcal{J}}$ is zero on $S(d)$.

The positiveness of $v^{f}$ in ${\cal N}$ follows by the definition itself. Next, if $t\in S_{j}(d)$, as already observed $T_{j}(d)$ is empty and $\tau_{j}(t)=0$, so that both terms in (4.11) are zero. Hence, $(v^{f}_{j})_{j\in\mathcal{J}}$ is zero on the maximum points of $d$ belonging to the edges, i.e. on $\prod_{j=1}^{M}S_{j}(d)$. If $x_{i}\in{\cal V}$ is a maximum point of $d$ on ${\cal N}$, then for all $j\in{\rm Inc}_{i}$ it holds that $S_{j}(d)$ is empty, $T_{j}(d)=\{\pi_{j}^{-1}(x_{i})\}$ and $\tau_{j}(t)=0$ for $t=\pi_{j}^{-1}(x_{i})$. Hence, the first term in (4.11) is zero. Since the set ${\rm Inc}_{i}^{+}(d)$ is necessarily empty too, the second term in (4.11) is also zero and $v^{f}$ results to be zero on the maximum that $d$ attains on ${\cal V}$.

Concerning the transition conditions (4.7), recall that $\sigma_{ij}(d)$ is always defined. Moreover, if $x_{i}$ is a transition node and $j\in{\rm Inc}_{i}^{-}(d)$, then $t=\pi_{j}^{-1}(x_{i})\in T_{j}(d)$, $\tau_{j}(t)=0$, $P_{j}(t)=t$ and (4.11) reduces to (4.7). The conservation of the flux (4.9) follows too.

It remains to obtain the regularity of $(v^{f}_{j})_{j\in\mathcal{J}}$ and to prove (4.8). Let $j\in\mathcal{J}$ be fixed. Assume that $S_{j}(d)=\emptyset$ and that $d$ is increasing along $e_{j}$. Then, $T_{j}(d)=\{\ell_{j}\}$, $P_{j}(t)=\ell_{j}$, $\eta_{j}(t)\,d^{\prime}_{j}(t)=1$ for all $t\in[0,\ell_{j})$ and $\tau_{j}(t)=\ell_{j}-t$ for all $t\in[0,\ell_{j}]$. Hence, with $i\in\mathcal{I}_{T}$ s.t. $j\in{\rm Inc}_{i}^{-}(d)$, using (4.7), (4.11) becomes

[TABLE]

In particular $v_{j}^{f}$ is continuous on $[0,\ell_{j}]$ and belongs to $W^{1,\infty}(0,\ell_{j})$ with $(v^{f}_{j})^{\prime}(t)=-f_{j}(t)$ a.e. $t\in(0,\ell_{j})$. Therefore, for a test function $\psi$ as in Definiton 4.1, and for $k\in\mathcal{I}$ s.t. $j\in{\rm Inc}_{k}^{+}(d)$, we have

[TABLE]

If $S_{j}(d)=\emptyset$ and $d$ is decreasing along $e_{j}$, it is easily seen that (4.13) still holds true by reflection.

Now, assume that $T_{j}(d)=\emptyset$ and $S_{j}(d)=\{\bar{t}\,\}$. Then, $\chi_{T_{j}(d)}\equiv 0$, $P_{j}(t)=\bar{t}$, $\eta_{j}(t)\,d_{j}^{\prime}(t)=1$ for $t\in(0,\bar{t}\,)$ and $\eta_{j}(t)\,d_{j}^{\prime}(t)=-1$ for $t\in(\bar{t},\ell_{j})$. Hence, for $t\in[0,\bar{t})$

[TABLE]

while, for $t\in(\bar{t},\ell_{j}]$,

[TABLE]

Therefore, $(v^{f}_{j})^{\prime}(t)=-f_{j}(t)$ a.e. $t\in(0,\bar{t})$, $(v^{f}_{j})^{\prime}(t)=f_{j}(t)$ a.e. $t\in(\bar{t},\ell_{j})$, $v^{f}_{j}\in W^{1,\infty}(0,\ell_{j})$, and we found again that $v^{f}_{j}(\bar{t})=0$. Denoting $x_{i}$ and $x_{k}$ the endpoints of $e_{j}$ and computing as before, we get

[TABLE]

It is worth noticing that in this case $j\in{\rm Inc}_{i}^{+}(d)\cap{\rm Inc}_{k}^{+}(d)$.

Finally, summing up (4.13) and (4.16) with respect to $j\in\mathcal{J}$, taking into account that $\psi\in C({\cal N})$ is zero on $\partial{\cal N}$ and that $(d,v^{f})$ satisfies (4.10), we obtain

[TABLE]

and (4.8) is proved. ∎

5 Uniqueness on ${\rm supp}(v^{f})$

This section is devoted to the proof of a uniqueness result for (4.2)–(4.7) over ${\rm supp}(v^{f})$, where for ${\rm supp}(v^{f})$ we intend $\Pi_{j=1}^{M}\pi_{j}({\rm supp}(v^{f}_{j}))$. In order to illustrate the complexity of the uniqueness problem, in the fashion of [6, 9], we introduce the function

[TABLE]

and the space $X:=\{u\in(W^{1,\infty}\cap C)({\cal N})\,:\,\eta(x)|Du(x)|\leq 1\text{ a.e. }x\in{\cal N}\}$.

It is easily seen that $X=\text{Lip}^{1}({\cal N}):=\{u\in C({\cal N})\,:\,|u(x_{1})-u(x_{2})|\leq{\cal D}(x_{1},x_{2})\,,\forall\ x_{1},x_{2}\in{\cal N}\}$ and that the distance function $d$ is the maximal nonnegative function in $X_{0}:=\{u\in X\,:\,u=0\text{ on }\partial{\cal N}\}$. Concerning $u^{f}$ we have the following.

Lemma 5.1**.**

*The function $u^{f}$ belongs to $X$, satisfies $0\leq u^{f}\leq d$ in ${\cal N}$ and it is the smallest nonnegative function among the nonnegative functions $u\in X$ such that $u=d$ on ${\rm supp}(f)$. Moreover, *

• (i)

$u^{f}=d$* in ${\rm supp}(v^{f})$ ;*

• (ii)

$u^{f}=d$* in ${\cal N}$ if and only if $S(d)\subset{\rm supp}(f)$ .*

Proof.

The function $u^{f}$ is a nonnegative and continuous function over ${\cal N}$ by definition. Furthermore

[TABLE]

implying $u\in X$. Indeed, let assume $u^{f}(x_{1})>0$ (otherwise the claim (5.1) is obvious) and let $y\in{\rm supp}(f)$ be a point realizing the maximum for $u^{f}(x_{1})$. Then, it holds

[TABLE]

If $d(y)-{\cal D}(x_{2},y)>0$, (5.1) follows by the triangular inequality for ${\cal D}$. Otherwise, $d(y)\leq{\cal D}(x_{2},y)$ and

[TABLE]

Next, if $x\in{\cal N}$ is such that $u^{f}(x)>0$ and $y\in{\rm supp}(f)$ realizes the maximum for $u^{f}(x)$, then $u^{f}(x)=d(y)-{\cal D}(x,y)\leq d(x)$. Therefore $u^{f}\leq d$ over ${\cal N}$. On the other hand, if $x\in{\rm supp}(f)$, then $u^{f}(x)\geq d(x)$ by definition again, and therefore $u^{f}=d$ on ${\rm supp}(f)$.

Consider now a nonnegative function $u\in X$ satisfying $u=d$ on ${\rm supp}(f)$. Take any $x\in{\cal N}$ such that $u^{f}(x)>0$ and let $y\in{\rm supp}(f)$ realize the maximum for $u^{f}(x)$. Then,

[TABLE]

and the first claim is totally proved.

To prove (i), we recall that $u^{f},d\in C({\cal N})$. Therefore, it is sufficient to obtain (i) on $\text{int}({\rm supp}(v^{f}))$. Let $x_{0}\in\text{int}({\rm supp}(v^{f}))$. We claim that there exists $x_{1}\in{\rm supp}(f)$ such that

[TABLE]

and $u^{f}(x_{0})=d(x_{0})$ follows by the definition and the properties of $u^{f}$.

Since $v^{f}(x_{0})\neq 0$, there exists at least one $j\in\mathcal{J}$ such that $x_{0}\in\overline{e}_{j}$ and $v^{f}_{j}(\pi_{j}^{-1}(x_{0}))>0$. For this $j$ fixed, two different cases can occur that give rise to an iteration procedure leading to (5.2).

First case: $S_{j}(d)=\{\bar{t}\,\}$ and $T_{j}(d)=\emptyset$. Then, the restriction of $d$ to ${\overline{e}}_{j}$ has a global maximum at $\pi_{j}(\bar{t})$ and $v^{f}_{j}(\bar{t})=0$ (see Theorem 4.2). Set $t_{0}=\pi_{j}^{-1}(x_{0})$. By formulae (4.14) and (4.15)

[TABLE]

Since $v^{f}_{j}(t_{0})>0$, it follows that $\bar{t}\neq t_{0}$ and there exists a set $A\subset(0,|\bar{t}-t_{0}|)$ with positive measure where the source $f_{j}$ is a.e. positive. Let $r_{1}\in A$ and set $x_{1}=\pi_{j}(t_{0}+r_{1}\,\eta_{j}(t_{0})\,d_{j}^{\prime}(t_{0}))$. Then, $x_{1}\in{\rm supp}(f)\cap e_{j}$, $x_{1}\neq x_{0}$, $x_{1}\neq\pi_{j}(\bar{t})$ and (5.2) follows in this case by (3.16).

Second case: $S_{j}(d)=\emptyset$. Then, $T_{j}(d)=\{\pi_{j}^{-1}(x_{i})\}$, with $x_{i}$ the endpoints of $e_{j}$ such that $i\in\mathcal{I}_{T}$ and $j\in{\rm Inc}_{i}^{-}(d)$, (see (3.15)). By formula (4.11)

[TABLE]

Again, since $v^{f}_{j}(t_{0})>0$, (at least) one of the two terms in the r.h.s. of (5.3) has to be positive. If the integral term is positive, we can proceed similarly to the first case to get (5.2). Otherwise there exists $k\in{\rm Inc}_{i}^{+}(d)$ such that $v^{f}_{k}(\pi_{k}^{-1}(x_{i}))>0$. If $S_{k}(d)$ is not empty, we can argue as in the first case along $\overline{e}_{k}$ with $x_{i}$ (endpoint of $e_{k}$) instead of $x_{0}$, to conclude that there exists $x_{1}\in\text{int}({\rm supp}(f))\cap e_{k}$ such that $d(x_{1})=d(x_{i})+{\cal D}(x_{i},x_{1})$. Furthermore, $d(x_{i})=d(x_{0})+{\cal D}(x_{0},x_{i})$ because either $x_{0}=x_{i}$ (this is the case for instance when $\tau_{j}(t_{0})=0$) or $x_{0}\neq x_{i}$ and one geodesic path from $x_{i}$ to $\partial{\cal N}$ has to pass through $x_{0}$ (otherwise $S_{j}(d)$ should not be empty). Since the distance $d$ is increasing along $\bar{e}_{k}$ from $x_{i}$ to $x_{1}$, a geodesic path from $x_{1}$ to $\partial{\cal N}$ has also to pass through $x_{i}$ and $x_{0}$, so that : ${\cal D}(x_{0},x_{1})={\cal D}(x_{0},x_{i})+{\cal D}(x_{i},x_{1})$. Hence,

[TABLE]

and the claim (5.2) follows once again.

Iteration procedure. If there is no $k\in{\rm Inc}_{i}^{+}(d)$ such that $v^{f}_{k}(\pi_{k}^{-1}(x_{i}))>0$ and $S_{k}(d)\neq\emptyset$, we apply the arguments of the second case to $x_{i}$ and $e_{k}$, instead of $x_{0}$ and $e_{j}$ respectively. Note that $d$ is increasing along $e_{k}$ from $x_{i}$ to the other endpoint of $e_{k}$. Therefore, $\tau_{k}(\pi_{k}^{-1}(x_{i}))>0$ and the first integral term in (5.3) is positive iff ${\rm supp}(f)\cap e_{k}\neq\emptyset$. If the latter holds true, we can proceed similarly to the first case to get (5.2). Otherwise, we iterate the procedure. Since ${\cal N}$ is finite and ${\rm supp}(f)\neq\emptyset$, after a finite number of steps we arrive necessarily to the source and obtain the claim (5.2).

It remains to prove (ii). Let assume that $S(d)\subset{\rm supp}(f)$. By the previous results, $u^{f}=d$ on $S(d)$. Let, $x\in{\overset{\circ}{\cal N}}\setminus S(d)$. Using the properties of $u^{f}$ and (3.16), there exists $y\in S(d)$ s.t.

[TABLE]

i.e. $u^{f}(x)=d(x)$. On the other hand, let assume that $u^{f}=d$ over ${\cal N}$ and that $S(d)\not\subset{\rm supp}(f)$. Take $x$ in $S(d)$ such that $x\notin{\rm supp}(f)$. Since $S(d)\cap\partial{\cal N}=\emptyset$, $u^{f}(x)=d(x)>0$, so that there exists $y\in{\rm supp}(f)$ realizing the maximum for $u^{f}(x)$. Let $z\in\partial{\cal N}$ such that $d(x)={\cal D}(x,z)$. Then, it holds

[TABLE]

giving $d(y)={\cal D}(y,z)={\cal D}(x,y)+d(x)$. Since $x\neq y$, the latter identity implies that a geodesic path from $y$ to $\partial{\cal N}$ pass through $x\in S(d)$ and this cannot be true (see the properties of $d$ in Section 3). Therefore, (ii) is totally proved. ∎

From Lemma 5.1 it follows that all the nonnegative functions $u\in X$ such that $u=d$ on ${\rm supp}(f)$, also satisfy $u=u^{f}=d$ on ${\rm supp}(v^{f})$. Therefore, if in addition $u=0$ on $\partial{\cal N}$, these functions are all good candidates to be the first component of the solution of (4.2)–(4.7), with $u^{f}$ the minimal one, since together with $v^{f}$ they satisfy (4.8). However, the transmission condition (4.7) is satisfied by each of $(u,v^{f})$ on the transition vertices $x_{i}$ where $v^{f}(x_{i})>0$, but nothing can be infered for the remaining transition vertices. Next lemma shows that if $(u,v)$ satisfies all the requirements of Definition 4.1 except the transmission condition, then the $u$-component is identified on ${\rm supp}(v^{f})$ as equal to $d$, but again nothing can be deduced about the $v$-component. The transmission condition (4.7) has to be henceforth a key tool for the uniqueness result, as we shall see in Theorem 5.3.

Lemma 5.2**.**

If $(u,v)$ satisfies the points (i), (ii) and (iii) of Definition 4.1 and the boundary condition (4.6), then the couple $(d,v)$ satisfies also (4.8) and $u=d=u^{f}$ on ${\rm supp}(v^{f})$.

Proof.

The proof will follow by the two identities below

[TABLE]

and by Lemma 5.1. Thanks to the properties of $u$ and $d$ we are allowed to use the test function $\psi=u-d$ in (4.8) for $(u,v)$ to obtain

[TABLE]

where the negative sign is due to the fact that $f\geq 0$ and $u\leq d$ in ${\cal N}$. Moreover, $|Dd(x)|=\frac{1}{\eta(x)}$ a.e. $x\in{\cal N}$ and $|Du(x)|=\frac{1}{\eta(x)}$ a.e. on $\{x\in{\overset{\circ}{\cal N}}\,:\,v(x)\neq 0\}$, so that

[TABLE]

The latter together with (5.5) give us

[TABLE]

and (5.4) follows. ∎

Theorem 5.3**.**

*If $(u,v)$ is a solution of (4.2)–(4.7) in the sense of Definition 4.1, then $u=d=u^{f}$ on ${\rm supp}(v^{f})$ and $v=v^{f}$ on $\Pi_{i=1}^{M}\bar{e}_{j}$. Moreover, if $S(d)\subset{\rm supp}(f)$, then $(u,v)=(d,v^{f})$ on ${\cal N}\times\Pi_{i=1}^{M}\bar{e}_{j}$. *

Proof.

We shall prove that $v=v^{f}$ on $\Pi_{i=1}^{M}\bar{e}_{j}$ by several steps, using the fact that $(d,v)$ satisfies (4.8). First, let us show that $v$ is zero on each maximum point of $d$ over ${\cal N}$, as it is the case for $v^{f}$ (see Theorem 4.2). For $j\in\mathcal{J}$ such that $S_{j}(d)=\{\bar{t}\,\}$ is not empty, let $n\in\mathbb{N}$ be sufficiently large so that $[\bar{t}-\frac{1}{n},\bar{t}+\frac{1}{n}]\subset(0,\ell_{j})$. Let $\psi_{n}=(\psi_{j,n})_{j\in\mathcal{J}}$ be a sequence of test functions uniformly bounded in $L^{\infty}({\cal N})$, that are zero on each edges except on $e_{j}$ and satisfying

[TABLE]

Then, taking into account the monotonicity property of $d$, (4.8) for $(d,v)$ gives us

[TABLE]

Passing to the limit as $n\to\infty$ we obtain $v_{j}(\bar{t})=0$. Hence,

[TABLE]

Let now $x_{i}\in{\cal V}$, $i\in\mathcal{I}_{T}$, be a maximum point for $d$ in ${\cal N}$. The argument is similar to the previous one, except that we have to take into account all the edges incident to $x_{i}$. Assume, without loss of generality, that $\pi_{j}^{-1}(x_{i})=\ell_{j}$ for all $j\in{\rm Inc}_{i}$. Hence, $\eta_{j}(t)d^{\prime}_{j}(t)=1$ for $t\in(0,\ell_{j})$ and all $j\in{\rm Inc}_{i}$. Let $n\in\mathbb{N}$ be sufficiently large so that $\ell_{j}-\frac{1}{n}>0$ for all $j\in{\rm Inc}_{i}$ and consider a sequence of test functions $\psi_{n}$ uniformly bounded in $L^{\infty}({\cal N})$ and satisfying

[TABLE]

Then, (4.8) for $(d,v)$ becomes

[TABLE]

Passing to the limit as $n\to\infty$ we obtain $\sum_{j\in{\rm Inc}_{i}}v_{j}(\ell_{j})=0$ and, since $v\geq 0$ in ${\cal N}$, $v_{j}(\ell_{j})=0$ for all $j\in{\rm Inc}_{i}$, i.e. $v(x_{i})=0$.

In order to prove that $v=v^{f}$ everywhere else in ${\cal N}$, we introduce the following partition of ${\cal E}$

[TABLE]

where

[TABLE]

${\cal E}_{0}$ contains all the edges $e_{j}$ such that $\bar{e}_{j}\cap S(d)\neq\emptyset$. In particular, if the edge $e_{j}$ has two boundary vertices as endpoints, then $e_{j}\in{\cal E}_{0}^{\prime}$. Furthermore, the two sets giving ${\cal E}_{0}$ are disjoints and they can not be both empty. Therefore the partition is well defined and finite.

Let assume that there exists $e_{j}\in{\cal E}_{0}^{\prime}$. Let $S_{j}(d)=\{\bar{t}\}$ and take $t\in(0,\bar{t}\,)$. Choose $n\in\mathbb{N}$ sufficiently large so that $[t-\frac{1}{n},\bar{t}+\frac{1}{n}]\subset(0,\ell_{j})$ and a sequence $\psi_{n}$ of test functions uniformly bounded in $L^{\infty}({\cal N})$, that are zero on each edges of the network except on $e_{j}$, with ${\rm supp}(\psi_{j,n})=[t-\frac{1}{n},\bar{t}+\frac{1}{n}]$ and $\psi_{j,n}=1$ on $[t,\bar{t}]$. Proceeding as above and taking into account (5.6), (4.8) for $(d,v)$ gives us

[TABLE]

Passing to the limit as $n\to\infty$, we get

[TABLE]

i.e. $v_{j}$ is given by (4.14) for $t\in(0,\bar{t}]$ as $v_{j}^{f}$. Formula (4.15) for $t\in[\bar{t},\ell_{j})$ can be obtained similarly. Hence $v=v^{f}$ holds on each $e_{j}\in{\cal E}_{0}^{\prime}$ and by continuity on $\bar{e}_{j}$.

Let assume now that there exists $e_{j}\in{\cal E}_{0}^{\prime\prime}$ and let $x_{i}\in{\cal V}$, $i\in\mathcal{I}_{T}$, be the endpoint of $e_{j}$ where $d$ has a maximum. Then, $e_{k}\in{\cal E}_{0}^{\prime\prime}$ for all $k\in{\rm Inc}_{i}$. Assume, without loss of generality, that $\pi_{k}^{-1}(x_{i})=\ell_{k}$ for all $k\in{\rm Inc}_{i}$. Let $t\in(0,\ell_{j}\,)$ be fixed and choose a sequence $\psi_{n}$ of test functions uniformly bounded in $L^{\infty}({\cal N})$, such that ${\rm supp}(\psi_{n})\subset\cup_{k\in{\rm Inc}_{i}}\bar{e}_{k}$ and for $n$ sufficiently large satisfies on $\bar{e}_{j}$

[TABLE]

while on $\bar{e}_{k}$, $k\in{\rm Inc}_{i}$, $k\neq j$,

[TABLE]

Observing that $\eta_{k}(r)d_{k}^{\prime}(r)=1$ on $(0,\ell_{k})$ for all $k\in{\rm Inc}_{i}$, by (4.8) for $(d,v)$ we obtain

[TABLE]

Passing again to the limit as $n\to\infty$ and recalling that $v(x_{i})=v^{f}(x_{i})=0$, we get

[TABLE]

i.e. $v_{j}$ is given by formula (4.12) in $[0,\ell_{j}]$ as $v_{j}^{f}$ and $v=v^{f}$ follows on $\bar{e}_{j}$, for each $e_{j}\in{\cal E}_{0}^{\prime\prime}$.

We are now ready to iterate the procedure. Assume that ${\cal E}_{1}$ is not empty, otherwise the proof is complete, and fix $e_{j}\in{\cal E}_{1}$. By the definition of the partition of ${\cal E}$ and the previous steps, there exists $i\in\mathcal{I}_{T}$ such that $j\in{\rm Inc}_{i}$ and

[TABLE]

Moreover, any such index $k$ in (5.9) belongs to ${\rm Inc}_{i}^{+}(d)$. Indeed, either $e_{k}\in{\cal E}_{0}^{\prime}$ so that $S_{k}(d)\neq\emptyset$, or $e_{k}\in{\cal E}_{0}^{\prime\prime}$ and $x_{i}$ is not a maximum point of $d$ on ${\cal N}$, otherwise $e_{j}$ should also belong to ${\cal E}_{0}^{\prime\prime}$. On the other hand, $j\in{\rm Inc}_{i}^{-}(d)$ since ${\cal E}\setminus{\cal E}_{0}$ does not contain any singular point of $d$ and $i\in{\cal I}_{T}$. Obviously, the same holds true for all $e_{j}\in{\cal E}_{1}$ such that $j\in{\rm Inc}_{i}$.

Assume again that $\pi_{k}^{-1}(x_{i})=\ell_{k}$ for all $k\in{\rm Inc}_{i}$. Let $t\in(0,\ell_{j})$ be fixed and choose as before a sequence $\psi_{n}$ of test functions uniformly bounded in $L^{\infty}({\cal N})$, such that ${\rm supp}(\psi_{n})\subset\cup_{k\in{\rm Inc}_{i}}\bar{e}_{k}$ and satisfying (5.7) and (5.8). Observing that, for $n$ large enough, $\eta_{k}(r)d_{k}^{\prime}(r)=-1$ on the support of $\psi_{k,n}$ for all $k\in{\rm Inc}_{i}^{+}(d)$, while $\eta_{k}(r)d_{k}^{\prime}(r)=1$ for all $k\in{\rm Inc}_{i}^{-}(d)$, and proceeding as before, (4.8) for $(d,v)$ gives

[TABLE]

Passing again to the limit as $n\to\infty$, we get

[TABLE]

Recalling that $j\in{\rm Inc}_{i}^{-}(d)$, the continuity of $v_{j}$ allows us to pass to the limit $t\to\ell_{j}$ in (5.10) to have

[TABLE]

i.e. the conservation of the flux (4.10) for $(d,v)$ in $x_{i}$. Plugging (5.11) into (5.10), the latter becomes

[TABLE]

It remains to prove that $v_{j}(\ell_{j})=v_{j}^{f}(\ell_{j})$, which implies that $v_{j}=v_{j}^{f}$ on $\bar{e}_{j}$ by formula (4.12). Recall that $x_{i}=\pi_{k}(\ell_{k})$ for all $k\in{\rm Inc}_{i}$. We distinguish the following two cases.

$(i)$ If $v^{f}(x_{i})=0$, then $v^{f}_{k}(\ell_{k})=0$ for all $k\in{\rm Inc}_{i}$ and in particular for all $k\in{\rm Inc}_{i}^{+}(d)$. By (5.9), the conservation of the flux (5.11) and the positivity of $v$, it follows that $v(x_{i})=0$ too.

$(ii)$ If $v^{f}(x_{i})\neq 0$, taking into account the transmission condition for $(d,v^{f})$, there exists (at least one) $k\in{\rm Inc}_{i}^{+}(d)$ such that $v_{k}^{f}(x_{i})>0$, while $v_{j}^{f}(x_{i})>0$ for all $j\in{\rm Inc}_{i}^{-}(d)$. For those indices, $\sigma_{ij}(u)$ and $\sigma_{ik}(u)$ are well defined with $\sigma_{ij}(u)=\sigma_{ij}(d)$ and $\sigma_{ik}(u)=\sigma_{ik}(d)$, since $u=d$ on ${\rm supp}(v^{f})$. Hence, $\{k\in{\rm Inc}_{i}^{+}(d)\,:\,v_{k}^{f}(x_{i})>0\}\subseteq{\rm Inc}_{i}^{+}(u)$ and ${\rm Inc}_{i}^{-}(d)\subseteq{\rm Inc}_{i}^{-}(u)$. Consequently, by (5.9) and the transmission condition (4.7) for $(u,v)$, it follows that $v_{j}(x_{i})>0$ for all $j\in{\rm Inc}_{i}^{-}(u)$. Moreover, if there exists $k\in{\rm Inc}_{i}^{+}(d)$ such that $v_{k}^{f}(x_{i})=0$, then $v_{k}(x_{i})=0$ and either $\sigma_{ik}(u)$ is not defined or $\sigma_{ik}(u)=1$, i.e. $k\in{\rm Inc}_{i}^{+}(u)$ (see (4.7) again). Resuming, $\{k\in{\rm Inc}_{i}^{+}(d)\,:\,v_{k}^{f}(x_{i})>0\}\subseteq{\rm Inc}_{i}^{+}(u)\subseteq{\rm Inc}_{i}^{+}(d)$, while ${\rm Inc}_{i}^{-}(d)={\rm Inc}_{i}^{-}(u)$ and the transmission condition for $(u,v)$ again gives us $v_{j}(\ell_{j})=v_{j}^{f}(\ell_{j})$.

Iterating similar arguments on the remaining ${\cal E}_{m}$, after a finite number of step we get the claim. ∎

Corollary 5.4**.**

If $(u,v)$ is a solution of (4.2)–(4.7) in the sense of Definition 4.1, then for all $i\in\mathcal{I}_{T}$ such that $v(x_{i})\neq 0$, the sets ${\rm Inc}_{i}^{\pm}(u)$ are not empty and satisfy

[TABLE]

Remark 5.5

It is possible to consider in the model an additional source term located at some of the transition vertices of the networks, $g:\{x_{i}\}_{i\in\mathcal{I}_{T}}\to[0,\infty)$. Since the additional sand poured by the source $g$ only influences the total mass rolling in the vertices, the sandpiles differential model is given by the same Monge-Kantorovich system discussed above but with the transition condition (4.7) replaced by

[TABLE]

with again $\sum_{j\in{\rm Inc}_{i}^{-}(u)}C_{ij}=\sum_{j\in{\rm Inc}_{i}^{-}(u)}K_{ij}=1$. The existence and uniqueness results still holds true with formula (4.11) replaced by

[TABLE]

6 An approximation scheme for the sandpiles problem

In this section we consider an approximation scheme to compute the solution $(d,v^{f})$ in (3.3) and (4.11).

Given the positive integers $M_{j}$, $j\in\mathcal{J}$, we define on each parameter’s interval $[0,\ell_{j}]$ the locally uniform partition $t_{m}^{j}:=m\,h_{j}$, $m=0,\dots,M_{j}+1$, with space step $h_{j}:=\ell_{j}/(M_{j}+1)$. The corresponding spatial grid on $e_{j}$ is then $\mathcal{G}_{j}:=\{x_{m}^{j}=\pi_{j}(t_{m}^{j})\,;\,m=1,\dots,M_{j}\}$, on $\bar{e}_{j}$ is $\bar{\mathcal{G}}_{j}:=\{x_{m}^{j}=\pi_{j}(t_{m}^{j})\,;\,m=0,\dots,M_{j}+1\}$, while $\mathcal{G}^{h}:=\{\mathcal{G}_{j}\,;\,j\in\mathcal{J}\}\cup{\cal V}$, with $h:=\max_{j\in\mathcal{J}}\{h_{j}\}$, is the grid on ${\cal N}$. We also set $\partial\mathcal{G}^{h}:=\partial{\cal N}$.

Next, for $x_{1},x_{2}\in\mathcal{G}^{h}$, we say that $x_{1}$ and $x_{2}$ are adjacent and we write $x_{1}\sim x_{2}$ if there exists $j\in\mathcal{J}$ such that $\text{Dist}(x_{1},x_{2})=h_{j}$. We call a discrete path ${\cal P}^{h}(x,y)$ connecting $x\in\mathcal{G}^{h}$ to $y\in\mathcal{G}^{h}$ any finite set $\{x_{0}=x,x_{1},\dots,x_{n}=y\}$ with $n\geq 1$, $x_{m}\in\mathcal{G}^{h}$ and $x_{m}\sim x_{m+1}$, $m=0,\dots,n-1$.

In order to compute an approximation $d^{h}$ of the distance $d$ in (3.3), we consider the following finite difference scheme for the eikonal equation (3.4) with homogeneous Dirichlet boundary condition

[TABLE]

It is easily seen that problem (6.1) admits a unique solution $\delta^{h}:\mathcal{G}^{h}\to\mathbb{R}$ given by

[TABLE]

where the minimum in (6.2) is taken over all discrete path ${\cal P}^{h}(x,y)$ and all $y\in\partial\mathcal{G}^{h}$. Indeed, a function $u^{h}:\mathcal{G}^{h}\to\mathbb{R}$ that is zero on $\partial\mathcal{G}^{h}$, satisfies the discrete eikonal equation in (6.1) iff

[TABLE]

and there exists at least one $y(x)\sim x$ for which the equality in (6.3) holds true. The discrete function in (6.2) satisfies the previous requirements and it is actually the unique solution. To prove the uniqueness claim, denote hereafter $z(x)$ the grid point which realizes the equality in (6.3) for $\delta^{h}$. Then, if $x\in\mathcal{G}^{h}\setminus\partial\mathcal{G}^{h}$ is such that $y(x)\sim x$ belongs to $\partial\mathcal{G}^{h}\cap\bar{e}_{j}$, $j\in\mathcal{J}$, it follows that $z(x)\in\partial\mathcal{G}^{h}\cap\bar{e}_{j}$ as well and $u^{h}(x)=\delta^{h}(x)$ since

[TABLE]

Now, let assume that $\sigma:=\max\{|u^{h}(x)-\delta^{h}(x)|\,;\,x\in\mathcal{G}^{h}\}>0$ and let $A^{\pm}:=\{x\in\mathcal{G}^{h}:u^{h}(x)-\delta^{h}(x)=\pm\sigma\}$. If $x\in A^{+}$, by (6.4), $z(x)\notin\partial\mathcal{G}^{h}$, but $z(x)\in A^{+}$ since

[TABLE]

Consequently, $u^{h}(x)-u^{h}(z(x))=\delta^{h}(x)-\delta^{h}(z(x))>0$, i.e. $u^{h}(x)>u^{h}(z(x))$. The latter implies that $u^{h}$ can not attain the minimum over the finite set $A^{+}$. Hence $A^{+}=\emptyset$. Changing the role between $u^{h}$ and $\delta^{h}$ it can be proved that $A^{-}=\emptyset$, so that $u^{h}\equiv\delta^{h}$ follows by contradiction.

It is also worth noticing that for any $x\in\mathcal{G}^{h}$ there may exist (at most) one grid point $y$ adjacent to $x$ such that $\delta^{h}(x)=\delta^{h}(y)$. Whenever the latter holds, we extend the grid $\mathcal{G}^{h}$ adding the new grid points $(x+y)/2$ and defining $\delta^{h}((x+y)/2):=\delta^{h}(x)+h/2$. The approximation $d^{h}$ of $d$ we shall consider here is then the continuous linear interpolation over ${\cal N}$ of the values attained by $\delta^{h}$ on the enlarged grid. For the sake of simplicity, we shall use the same notations as before for the enlarged grid and its grids points but we replace the uniform space step $h_{j}$ on $e_{j}$ with the non-uniform space steps

[TABLE]

Thanks to the previous procedure, the $d^{h}$ function shares the same nice properties of the distance (3.3) discussed in Section 3. Indeed, let define the forward difference quotients of $d^{h}$ on each edge $e_{j}$

[TABLE]

the corresponding signs $\varsigma_{m}^{j}:={\rm sgn}[\Delta_{m}^{j}(d^{h})]$ and the set of critical points of $d^{h}$ inside the edge $e_{j}$

[TABLE]

Then, with the same definition (3.9), for the slope of $d^{h}$ at the vertices $x_{i}$, and (3.10), it is straightforward to prove that $d^{h}$, $\sigma_{ij}(d^{h})$ and $S_{j}(d^{h})$ satisfy properties (i) to (iv) in Proposition 3.3. In particular, the critical points of $d^{h}$ are maximum points and they are attained uniquely in grid points.

Next, in order to obtain an approximation of formula (4.11), using the same definition (3.14)-(3.15) for the projection set in $\bar{e}_{j}$ of $d^{h}$, we set, for $m=0,\dots,M_{j}+1$,

[TABLE]

[TABLE]

and we define the projection of $t_{m}^{j}$ onto $\Sigma_{j}(d^{h})$ as

[TABLE]

With the latter we can finally approximate the function $v^{f}$ over $\Pi_{j=1}^{M}\bar{\mathcal{G}}_{j}$ by means of

[TABLE]

if $\varsigma_{m}^{j}=1$, and

[TABLE]

if $\varsigma_{m}^{j}=-1$. The last step is to define $v^{f,h}$ as the continuous linear interpolation of the values $v^{f,h}_{j}(t_{m}^{j})$ over $\Pi_{i=1}^{M}\bar{e}_{j}$.

Theorem 6.1**.**

$(d^{h},v^{f,h})$* converges as $h\to 0$ to the solution $(d,v^{f})$ of the sand piles problem, uniformly in ${\cal N}\times\Pi_{j=1}^{M}\overline{e}_{j}$.*

The proof of the above convergence result is quite standard and we leave it to the reader. Indeed, the convergence of $d^{h}$ toward $d$ can be easily proved by a combination of classical arguments in viscosity solution theory and the Comparison Principle in [11]. Once the uniform convergence of $d^{h}$ is obtained, the convergence of $v^{f,h}$ toward $v^{f}$ easily follows by the comparison of the explicit formulas (4.11) and (6.5)-(6.6) observing that the singular set of $d^{h}$ converges to the singular set of $d$.

7 SPNET and numerical tests

In this section we first briefly introduce SPNET (Sand Piles on NETworks), an easy-to-use program written in C we developed for the numerical approximation of the sand pile problem. The interested reader can download the software at http://www.dmmm.uniroma1.it/$\sim$fabio.camilli/spnet.html. Next, we shall consider three numerical tests showing the features of the proposed method and providing empirical convergence analysis.

SPNET takes in input .net files, which are simple text files containing lists of vertices and edges, formatted as follows:

#SPNET

#v1 x1 y1 type1

...

#vN xN yN typeN

#e1 start1 end1 n1 f1(t) eta${}^{-1}_{1}$(t)

...

#eM startM endM nM fM(t) eta${}^{-1}_{M}$(t)

where

• (i)

#SPNET is just a header to recognize .net files;

• (ii)

#vi xi yi typei defines a vertex with typei equal to b (for boundary) or t (for transition) with coordinates (xi,yi);

• (iii)

#ej startj endj nj fj(t) eta${}^{-1}_{j}$(t) defines an edge connecting the vertices with indices startj and endj. The edge is parametrized from the vertex startj to the vertex endj using nj discretization nodes. Finally, fj(t) and eta${}_{j}^{-1}$(t) are respectively the sand source (4.1) and the inverse of the spatial inhomogeneity (3.1) on the current edge, defined as symbolic analytic functions of the parameter $t$ ranging in $[0,1]$.

The program checks for syntax errors in the input .net file, then quickly computes the solution pair $(d^{h},v^{f,h})$ and plots the results in a gnuplot window (http://www.gnuplot.info). The current view can be saved into a .pdf or .svg file, and some information about the solution can be printed on the screen.

Let us spend some words on the actual implementation. After allocating proper data structures for the network, the program first computes the approximation $d^{h}$ of the distance function $d$. This is readily done using the scheme (6.1) in a fast marching fashion for all the edges, suitably modified to handle the transition vertices. More precisely, one has to correctly propagate information to all the incoming edges when a transition vertex is encountered (for further details we refer the reader to [14]). The second step consists in computing all relevant data derived from $d^{h}$, including its slopes, singular set and projections. Finally the edges are dynamically processed, starting from those containing points of the singular set, i.e. using the partition $({\cal E}_{0},{\cal E}_{1},\dots)$ of ${\cal E}$ (see Theorem 5.3). The approximation $v^{f,h}$ of the rolling component $v^{f}$ is then computed recursively from ${\cal E}_{0}$ up to the boundary $\partial{\cal N}$, using the values obtained in the previous iterations, according to formula (6.5)-(6.6). In all the following examples the amount of rolling sand at the transition vertices is assigned along the incoming edges by prescribing the coefficients $C_{ij}$ equal to $1/N_{i}^{-}(d)$ (uniform repartition of the incoming sand), but this is not a restriction for the code that indeed can be used in the general setting.

**Test 1. ** We start considering a simple example of planar network composed of four vertices $\{x_{i}\}_{i=0,1,2,3}$, with $\mathcal{I}_{T}=\{0\}$ and $\mathcal{I}_{B}=\{1,2,3\}$, and three edges $\{e_{j}\}_{j=1,2,3}$ connecting each one $x_{0}$ to $x_{j}$, $j\neq 0$, with parameter’s sets $[0,\ell_{j}]=[0,1/2]$ for $j=1,2,$ and $[0,\ell_{3}]=[0,1]$, (see Figure 1 (a)). Moreover, we assume $\eta\equiv 1$ on the whole network, so that the metric (3.2) coincide with (2.1), and the sand source $f$ is given edgewise by (see Figure 1 (b))

[TABLE]

In this situation, the solution of the sandpile problem can be computed explicitly. In particular, the distance $d$ from the boundary is given by

[TABLE]

where the value $t=\frac{1}{4}\in[0,\ell_{3}]$ gives the unique singular point of the distance $d$ over the network. It follows that ${\rm Inc}_{0}^{+}=\{3\}$, ${\rm Inc}_{0}^{-}=\{1,2\}$, and we get $v^{f}$ by (4.11) as

[TABLE]

where $\frac{7}{32}$ is the contribution of the rolling layer on $e_{3}$ to $x_{0}$, which is uniformly split in $e_{1}$, $e_{2}$. Moreover, since $f\neq 0$ all over the network, by Theorem 5.3 the couple $(d,v^{f})$ is the unique solution of the sandpile problem.

In Figure 2 we show the numerical solution $(d^{h},v^{f,h})$ computed by SPNET using a uniform discretization step $h=10^{-2}$ for all the edges.

Note that the component $v^{f,h}$ is multivalued at $x_{0}$ as $v^{f}$.

We remark that since in this test the sand source $f$ is a linear function on each edge, the trapezoidal quadrature rule in (6.5)-(6.6) computes the exact values at the grids points. Since $d$ is also a piecewise linear function, the only source of error is given by the wrong localization of the singular point and we never see an error if we compare the exact and the approximate solution only on the grid points. Hence, for a more fair comparison, we evaluate the error on the whole network introducing in this way an additional interpolation error for $v^{f,h}$. To sample both the $L^{\infty}$ and $L^{1}$ errors, we use a very fine grid of $10^{4}$ nodes per edge.

In Figure 3 we show the errors (in logarithmic scale) against a uniform discretization step $h$ ranging from $10^{-1}$ to $10^{-3}$.

We readily observe for the distance $d^{h}$ an $L^{\infty}$ error of order $1$ and an $L^{1}$ error of order $2$ (see respectively Figures 3 (a) and 3 (b)). Moreover, we confirm that both errors vanish (up to machine precision and round-off errors) for all the grids containing the singular point.

The behavior of the errors for the rolling component $v^{f,h}$ is instead slightly different. Indeed, we see in Figures 3 (c) and 3 (d) that both are at worst of order $1$, but become of order $2$ (they are not almost zero as for $d$!) if the corresponding grids contain the singular point. As explained, in this situation the errors in $v^{f,h}$ are zero only on the network grid, elsewhere we pay for the interpolation.

**Test 2. ** We now extend the previous network adding two additional edges $e_{4}$ and $e_{5}$, connecting respectively $x_{2}$ to $x_{1}$ and $x_{1}$ to $x_{3}$. Moreover, we set $\mathcal{I}_{T}=\{0,1,2\}$ and $\mathcal{I}_{B}=\{3\}$ (see Figure 4 (a)). We assume the sand source $f$ and the spatial dishomogeneity $1/\eta$ given edgewise in normalized coordinates $t\in[0,1]$ respectively by (see Figure 4 (b))

[TABLE]

In Figure 5 we show the numerical solution $(d^{h},v^{f,h})$ computed by SPNET using a uniform discretization step $h=10^{-2}$ for all the edges.

Note that $v^{f,h}$ is multivalued at $x_{0}$ and $x_{3}$ and continuous at the other vertices. Moreover, $v^{f,h}$ is zero on $e_{1}$, $e_{5}$ and part of $e_{4}$, and we loose uniqueness of the solution pair $(d,v^{f})$.

**Test 3. ** We finally consider a more complex network composed of 6 vertices and 9 edges, which is a two-level pre-fractal for the Sierpiński triangle. The extremal vertices $x_{0}$, $x_{1}$ and $x_{2}$ are the boundary vertices, whereas the internal ones are the transition vertices (see Figure 6 (a)).

Again, we assume $\eta\equiv 1$ on the whole network and we take a sand source $f$ similar to the one in Test 2, supported in a sub-interval of the edges $e_{7}$ and $e_{8}$ (see also Figure 6 (b))

[TABLE]

In Figure 7 we show the numerical solution $(d^{h},v^{f,h})$ computed by SPNET using a uniform discretization step $h=10^{-2}$ for all the edges. We observe a more rich structure. The distance $d^{h}$ has three singular points and $v^{f,h}$ is continuous at the vertices $x_{0}$ and $x_{4}$ and multivalued at the remaining vertices. Moreover, $v^{f,h}$ is zero on $e_{9}$ where uniqueness fails.

Appendix A Appendix A

Proof of Proposition 3.3  We shall prove only claims (i) and (ii) since the remaining (iii) and (iv) are straightforward consequences. The statements (i) and (ii) follow on ${\cal E}$ observing that if $w$ is a viscosity solution of

[TABLE]

with $g\in C([a,b])$, $g>0$, then one of the following cases is true:

• (i)

$w$ is a piecewise $C^{1}$ function over $[a,b]$ with exactly one singular point $\bar{t}\in(a,b)$, $w^{\prime}>0$ on $[a,\bar{t})$ and $w^{\prime}<0$ on $(\bar{t},b]$;

• (ii)

$w\in C^{1}([a,b])$ and either $w^{\prime}>0$ or $w^{\prime}<0$ on $[a,b]$.

In fact, $w$ cannot have minimum point inside $(a,b)$. Otherwise, if $w$ has a local minimum at $t_{0}\in(a,b)$, the constant function $\phi(t)\equiv w(t_{0})$ is a test function such that $(w-\phi)$ has a local minimum at $t_{0}$. Since $\phi^{\prime}(t_{0})\equiv 0$ and $g$ is positive, we get a contradiction with (3.6), the definition of supersolution at $t_{0}$. On the other hand, if $w$ has two maximum points inside $(a,b)$, say $t_{1}$ and $t_{2}$, then, because of the previous property, $w$ is constant in the interval $(t_{1},t_{2})$ and we get again a contradiction by taking a test function constant over $(t_{1},t_{2})$.

Finally, the fact that $u$ does not atteint a local minimum in a transition vertex too follows applying similar arguments. $\Box$

Proof of Proposition 3.4  Let $x_{0}\in{\cal N}$ be fixed and $u(\cdot):={\cal D}(x_{0},\cdot)$. The continuity of $u$ is a straightforward consequence of the definition of ${\cal D}$ itself. Next, we will show only that $u$ satisfies the supersolution condition (3.7) at $x_{i}$, $i\in{\cal I}_{T}$, since the proofs of (3.5) and (3.6) can be obtained with the same type of reasoning. Moreover, (3.5) and (3.6) follow with equality.

Let $\phi\in C^{1}({\cal N})$ be a test function such that $(u-\phi)$ has a local minimum at the transition vertex $x_{i}$. If in turn $x_{0}$ is a transition vertex, we assume that $x_{i}\neq x_{0}$. Let ${\cal P}(x_{0},x_{i})=(\bar{e}_{j_{1}}\cap\bar{e}(x_{0}),\dots,\bar{e}_{j_{n+1}})$ be a geodesic path realizing ${\cal D}(x_{0},x_{i})$. Then, $(u_{j_{n+1}}-\phi_{j_{n+1}})$ has a local minimum at $\pi_{j_{n+1}}^{-1}(x_{i})$ and

[TABLE]

for all $x\in e_{j_{n+1}}$ in a neighborhood of $x_{i}$. Since for those $x$, ${\cal P}(x_{0},x)=(\bar{e}_{j_{1}}\cap\bar{e}(x_{0}),\dots,\bar{e}_{j_{n+1}}\cap\bar{e}(x))$ is a geodesic path realizing ${\cal D}(x_{0},x)$, it holds that ${\cal D}(x_{0},x_{i})-{\cal D}(x_{0},x)={\cal D}(x,x_{i})$ and

[TABLE]

Recalling (2.3), the latter gives us easily

[TABLE]

and (3.7) follows.

The proof that $d(x)=\min_{y\in\partial{\cal N}}{\cal D}(x,y)$ is a viscosity solution of (3.4) is quite standard. Indeed, $d$ is a supersolution since it is the minimum of the finite number of supersolutions ${\cal D}(\cdot,y)$, $y\in\partial{\cal N}$. Moreover, it is a subsolution since it is the maximum of a family of Lipschitz continuous subsolution and therefore we can apply the Perron method (see [2, Thm. 2.12]). Finally, to prove that $d$ is a supersolution at the transition vertices $x_{i}$, let $\phi$ be a test function such that $(d-\phi)$ has a local minimum at $x_{i}$ and $\overline{y}\in\partial{\cal N}$ such that $d(x_{i})={\cal D}(x_{i},\overline{y})$. Hence, for $x$ in a neighborhood of $x_{i}$ we have

[TABLE]

and therefore ${\cal D}(\cdot,\overline{y})-\phi(\cdot)$ has a minimum at $x_{i}$. Following the arguments above, we get (3.7). Since the uniqueness of the solution of (3.4) with homogeneous Dirichlet boundary condition is a consequence of the comparison theorem in [11], the proof is complete. $\Box$

Acknowledgment This research was initiated while the third author was visiting the Department of Basic and Applied Sciences for Engineering at “Sapienza” University of Rome. The third author acknowledges the support of the french “ANR blanche” project Kibord : ANR-13-BS01-0004 and of “Progetto Gnampa 2015: Network e controllabilità”.