Analysis of the implicit upwind finite volume scheme with rough coefficients

TL;DR

This paper analyzes the implicit upwind finite volume scheme for the linear continuity equation with low regularity coefficients, establishing a convergence rate of at least 1/2 in a weak sense.

Contribution

It provides the first convergence rate analysis of the scheme in the DiPerna-Lions setting with rough coefficients and unstructured meshes.

Findings

Convergence rate of at least 1/2 in logarithmic Kantorovich-Rubinstein distance.

Error bounds applicable to low regularity velocity fields.

Validation of the scheme's effectiveness in rough coefficient scenarios.

Abstract

We study the implicit upwind finite volume scheme for numerically approximating the linear continuity equation in the low regularity DiPerna-Lions setting. That is, we are concerned with advecting velocity fields that are spatially Sobolev regular and data that are merely integrable. We prove that on unstructured regular meshes the rate of convergence of approximate solutions generated by the upwind scheme towards the unique distributional solution of the continuous model is at least 1/2. The numerical error is estimated in terms of logarithmic Kantorovich-Rubinstein distances and provides thus a bound on the rate of weak convergence.

Equations

$$\{\begin{array}{rlrl}{\partial }_t\rho +\nabla \cdot \left(u\rho \right)& =f& & \text{in }(0,T)\times \Omega ,\\ \rho (0,\cdot )& ={\rho }^0& & \text{in }\Omega .\end{array}$$

$$u\cdot \nu =0\text{on }\partial \Omega ,$$

$${\int }_{\Omega }f(t,x)dx=0\text{for all }t\in (0,T).$$

$${\int }_{\Omega }\rho (t,x)dx={\int }_{\Omega }{\rho }^0(x)dx\text{for all }t\in (0,T).$$

$$u\in L^1((0,T);W^{1,p}(\Omega ))\text{and}(\nabla \cdot u{)}^{-}\in L^1((0,T);L^{\infty }(\Omega ))$$

$$\Vert \rho {\Vert }_{L^{\infty }(L^q)}\le {\Lambda }^{1-\frac{1}{q}}\left(\Vert {\rho }^0{\Vert }_{L^q}+\Vert f{\Vert }_{L^1(L^q)}\right),$$

$$\rho \in L^{\infty }((0,T);L^q(\Omega ))\text{where}\text{nicefrac}1p+\text{nicefrac}1q=1.$$

$$h:=\underset{K\in \mathcal{T}}{\max }\mathrm{diam}K.$$

$$\Vert \psi {\Vert }_{L^1(\partial K)}$$

$$\Vert \psi -{\psi }_K{\Vert }_{L^1(K)}$$

$$\frac{|\partial K|}{|K|}\le \frac{C}{h},$$

$${\rho }_K^0:=\text{vbox}-{\int }_K{\rho }^{0}dx.$$

$$f_K^n:=\text{vbox}-{\int }_{t^n}^{t^{n+1}}\text{vbox}-{\int }_{K}fdxdt.$$

$$u_{KL}^n:=\text{vbox}-{\int }_{t^n}^{t^{n+1}}\text{vbox}-{\int }_{K\mid L}u\cdot {\nu }_{KL}d{\mathcal{H}}^{d-1}dt.$$

$$\frac{\rho_{K}^{n+1}-\rho_{K}^{n}}{{\delta}}+\sum_{L\sim K}\frac{\lvert K{\scriptstyle\mid}L\rvert}{\lvert K\rvert}\;\big{(}u_{KL}^{n+}\,\rho_{K}^{n+1}-u_{KL}^{n-}\,\rho_{L}^{n+1}\big{)}=f_{K}^{n}$$

$$\rho_{{\delta},h}(t,x):=\rho^{n}_{K}\qquad\text{for a.e.\ }(t,x)\in\bigl{[}t^{n},t^{n+1}\bigr{)}\times K.$$

$$\frac{{\rho }_K^{n+1}-{\rho }_K^n}{\delta }+\sum _{L\sim K}\frac{|K\mid L|}{|K|}{u}_{KL}^n\frac{{\rho }_K^{n+1}+{\rho }_L^{n+1}}{2}+\sum _{L\sim K}\frac{|K\mid L|}{|K|}\left|u_{KL}^n\right|\frac{{\rho }_K^{n+1}-{\rho }_L^{n+1}}{2}=f_K^n,$$

$$\frac{q-1}{q}{\int }_I\Vert (\nabla \cdot u{)}^{-}{\Vert }_{L^{\infty }}dt\le \frac{\kappa -1}{\kappa }\text{for all intervals I of length }{\delta }_{\max }(\kappa ).$$

$${\rho }^0\in L^q(\Omega ),\text{and}f\in L^1((0,T);L^q(\Omega ))\cap L^{\infty }((0,T);W^{-1,1}(\Omega )),$$

$$u\in L^1((0,T);W^{1,p}(\Omega ))\text{with}(\nabla \cdot u{)}^{-}\in L^1((0,T);L^{\infty }(\Omega )),$$

$$u\in L^{\infty }((0,T)\times \Omega ).$$

$$\inf_{\pi\in\Pi\left(\rho(t,\cdot),\rho_{{\delta},h}(t,\cdot)\right)}\iint\log\bigg{(}\frac{|x-y|}{r}+1\bigg{)}d\pi(x,y)\lesssim 1+\frac{\sqrt{h\,T\,\|u\|_{L^{\infty}}}+\sqrt{{\delta}\,T\phantom{\|}\!}\,\|u\|_{L^{\infty}}}{r}\!$$

$${\rho }_{\delta ,h}⟶\rho \text{weakly with rate at most h+δ},$$

$$\underset{h\to 0}{\lim }{h}^{-\epsilon }\Vert \rho -{\rho }_h{\Vert }_{L^1((0,T)\times \Omega )}\gtrsim 1.$$

$$\Vert \nabla {\rho }_{\delta ,h}{\Vert }_{L^1((0,T)\times \Omega )}\lesssim \sqrt{\frac{T}{hU}}\text{and}\Vert {\partial }_t{\rho }_{\delta ,h}{\Vert }_{L^1((0,T)\times \Omega )}\lesssim \sqrt{\frac{T}{\delta }}.$$

$${\int }_{\Omega }{\rho }_{1}dx={\int }_{\Omega }{\rho }_{2}dx,$$

$$\pi [A\times \Omega ]={\int }_A{\rho }_{1}dx\text{and}\pi [\Omega \times A]={\int }_A{\rho }_{2}dx\text{for all measurable }A\subset \Omega .$$

$$\int_{\Omega}\big{(}\zeta(x)+\theta(y)\big{)}\,d\pi(x,y)=\int_{\Omega}\zeta\,\rho_{1}\,dx+\int_{\Omega}\theta\,\rho_{2}\,dx\qquad\text{for all }\zeta,\theta\in C(\overline{\Omega}),$$

$$D_r({\rho }_1,{\rho }_2):=\underset{\pi \in \Pi ({\rho }_1,{\rho }_2)}{\inf }\iint \log \left(\frac{|x-y|}{r}+1\right)d\pi (x,y).$$

$$D_r({\rho }_1,{\rho }_2)=\underset{\zeta }{\sup }\left\{{\int }_{\Omega }\zeta ({\rho }_1-{\rho }_2)dx:|\zeta (x)-\zeta (y)|\le \log \left(\frac{|x-y|}{r}+1\right)\right\},$$

$$D_r({\rho }_1,{\rho }_2)=\iint \log \left(\frac{|x-y|}{r}+1\right)d{\pi }_{\mathrm{opt}}(x,y)={\int }_{\Omega }{\xi }_{\mathrm{opt}}({\rho }_1-{\rho }_2)dx.$$

Full text

1. Introduction

The continuity equation is one of the most fundamental and at the same time most elementary partial differential equation with applications in a wide range of problems from physics, engineering, biology or social science. It describes the conservative transport of a quantity (e.g., a mass or number density, temperature, concentration, or tracer) by a given velocity field. In many important examples, the velocity field itself is related to the actual configuration of this quantity by means of a momentum-type equation, thermodynamic law or other basic principles. In this paper, we are interested in the purely linear model in which such a feedback of the actual configuration on the velocity field is neglected. For the sake of a larger applicability, we will allow for external sources and sinks.

We are thus concerned with the linear inhomogeneous continuity equation in a bounded Lipschitz domain $\Omega$ in $\mathbf{R}^{d}$ and time interval $(0,T)$, that is,

[TABLE]

Here, $\rho$ is the evolving (scalar) quantity with initial configuration $\rho^{0}$, $u$ the velocity field and $f$ is the source-sink distribution. We will often refer to $\rho$ as a density even if it may take negative values.

We will assume in this paper that there is no loss of mass by transport across the boundary, that is, we suppose that the velocity field is tangential at $\partial\Omega$,

[TABLE]

if $\nu$ denotes the outer unit normal vector along the boundary. Under this hypothesis we can apply a standard transformation which ensures that the total sources and sinks are balanced. We will accordingly demand that

[TABLE]

The previous two assumptions together imply that the continuity equation is indeed conservative in the sense that

[TABLE]

It is obvious that this identity is formally true and it indeed holds under the assumptions of this paper which will be specified in the following.

In many relevant applications, for instance, in turbulent transport of mass or heat, neither the velocity field nor the transported density are expected to be regular functions. The minimal mathematical requirement for guaranteeing well-posedness of the Cauchy problem (1) is that $u$ has Sobolev or bounded variation ($BV$) regularity in the spatial variable. This is the setting studied in the ground-breaking papers of DiPerna and Lions [9] and Ambrosio [2]. We will focus on the case of Sobolev vector fields, which is the setting originally studied by DiPerna and Lions, and will thus assume that

[TABLE]

for some $p\in(1,\infty]$. Here the superscript minus sign indicates the negative part of the divergence. While the first hypothesis implies the so-called renormalization property, which, in a certain sense, is the verification of the chain rule for solutions of (1), the second hypothesis yields the validity of the a priori estimate

[TABLE]

provided that the right-hand side is finite, where $\Lambda:=\exp\left(\|(\nabla\cdot u)^{-}\|_{L^{1}(L^{\infty})}\right)$ is the compressibility constant associated with $u$. Here and in the following, we use the shorter notation $L^{r}(X)$ instead of $L^{r}((0,T);X)$ for any Banach space $X$ and $r\in[1,\infty]$. In principle, the results in [9, 2] allow for the consideration of densities that are merely summable. In this case, the relevant notion of solutions is that of renormalized solutions. In the present paper, we are interested in the smaller class of distributional solutions, which are well-defined if

[TABLE]

This integrability assumption is consistent with the a priori estimate (5).

The aim of this paper is to provide an error estimate for the implicit upwind finite volume approximation on unstructured meshes of distributional solutions to the continuity equation (1) in the DiPerna–Lions setting [9]. The upwind scheme is the most classical stable monotone and mass preserving numerical method for hyperbolic conservation laws, whose underlying idea is the numerical approximation of the upstream transport of cell averages. Convergence rates for upwind methods were intensively studied in the past forty years, starting with the pioneering work by Kuznetsov in 1976 [13]. In the linear case with Lipschitz regular vector fields and $BV$ or $H^{1}$ initial data, the optimal results on unstructured meshes are due to Johnson and Pitkäranta [12]111In fact, Johnson and Pitkäranta consider the more general discrete Galerkin approximation., Vila and Villedieu [23], Merlet and Vovelle [15, 14], Delarue and Lagoutière [7], and Aguillon and Boyer [1]. All these works provide $\nicefrac{{1}}{{2}}$ convergence rates in terms of the mesh size $h$. Concerning the optimality of these results, we refer to the papers by Peterson [17] and Tang and Teng [22]. The case of lower regularity data was treated only very recently. Under a one-sided Lipschitz condition on the vector field, Delarue, Lagoutière and Vauchelet estimated the approximation error by $\mathcal{O}(h^{\nicefrac{{1}}{{2}}})$ for arbitrary measure-valued solutions [8]. Using similar (probabilistic) techniques, the authors of the present paper derived the same rate for distributional solutions in the class (6) under the assumption (4), cf. [18]. Unfortunately, both works [8] and [18] are restricted to Cartesian meshes.

In the present work, convergence rates for upwind methods on unstructured meshes in the low regularity framework are obtained for the first time. As in the classical setting, we prove a $\mathcal{O}(h^{\nicefrac{{1}}{{2}}})$ error bound in the case of rough densities (6) and Sobolev vector fields (4). Less importantly, we treat the implicit scheme instead of the explicit scheme studied earlier in [18], but we’re convinced that the explicit one could be handled in a similar fashion. The step from Cartesian to arbitrary meshes appears to be substantial: While on Cartesian meshes the upwind scheme allows for an interpretation as a finite difference method, the scheme is genuinely of finite volume type on unstructured meshes. In our previous work [18] (and that of Delarue, Lagoutière and Vauchelet [8]) the Cartesian mesh geometry allowed for a comparison of the Lagrangian flow associated with the continuous problem with a stochastic flow associated with a probabilistic (Markov chain) interpretation of the upwind scheme. On unstructured meshes this comparison seems to fail. Moreover, an adaptation of the construction used for the transport equation with Lipschitz regular vector field by Delarue and Lagoutière [7] to the continuity equation in the low regularity regime is not apparent. To overcome this difficulty, in the present work we choose to work mostly with the Eulerian specification of the problem rather than the Lagrangian one considered in [8, 18]. Instead of comparing the flows, we directly estimate the distance of continuous and approximate solution. For this purpose, we derive a number of new optimal stability estimates for continuity equations that built up on and further extend techniques recently established in [19, 20]. The Lagrangian formulation of transport enters our analysis only through the superposition principle.

The focus on density functions of low regularity comes along with a change in topology: The results established in [8, 18] and the present paper quantify the rate of weak convergence (of measures). This is contrasted by the “classical” setting with $BV$ or $H^{1}$ densities considered in [13, 12, 23, 15, 14, 7], where optimal rates in strong Lebesgue norms can be proved. In [18] we show that this change of topology is not at all a pathology of the applied method. In fact, quite elementary examples indicate that, firstly, no rates exist that are uniform in the initial data and, secondly, rate $\nicefrac{{1}}{{2}}$ weak convergence is (almost) sharp in the sense that for any small $\varepsilon$ there exists configurations that generate rates of at most $\nicefrac{{1}}{{2}}+\varepsilon$. In this respect, the qualitative convergence results in strong norms for upwind schemes with rough coefficients obtained by Walkington [25] and Boyer [4] are optimal, too.

The paper is organized as follows: Section 2 contains the precise definition of the implicit upwind finite volume scheme. In Section 3 we present and discuss our main results. Section 4 provides an overview on relevant facts about Kantorovich–Rubinstein distances. In Section 5 we derive stability and weak $BV$ estimates for the upwind scheme. The final Section 6 is devoted to the error analysis. We conclude this paper with a short appendix on generalized means.

2. The upwind scheme: I. Definition

The upwind scheme is a classical finite volume method for approximating hyperbolic conservation laws. It approximates, roughly speaking, the evolution of volume averages by means of the flux over their boundaries, where it only uses the upstream values to update neighboring cell averages in each time step. See also the monograph [11] for references and further details.

In this section, we describe the upwind scheme for the continuity equation with arbitrary source-sink distribution (1). Since there is no mass flux across the domain boundary, cf. (2), there is no need to restrict the geometry of $\Omega$ in what follows.

We consider as tessellation $\mathcal{T}$ of $\Omega$ a finite disjoint polyhedral covering of $\Omega$. That is a set of finite many closed connected sets $K$ with disjoint interiors and such that $\overline{\Omega}=\bigcup_{K\in\mathcal{T}}K$. Moreover, for each $K\in\mathcal{T}$ exists a polyhedral set $K^{\prime}$ such that $K=K^{\prime}\cap\overline{\Omega}$. These sets $K$ are called control volumes or simply cells. The interior boundary of each cell, that is $\partial K\setminus\partial\Omega$, is the union of finitely many flat closed and connected $d-1$ dimensional faces. We write $K\sim L$ whenever $K$ and $L$ are two neighboring cells and we denote by $K{\scriptstyle\mid}L$ the joint edge $K\cap L$. The unit normal vector on the edge $K{\scriptstyle\mid}L$ pointing from $K$ to $L$ will be denoted by $\nu_{KL}$, so that $\nu_{KL}=-\nu_{LK}$. The definition of the scheme involves the edge size $|K{\scriptstyle\mid}L|$ and the volume $|K|$, where, by a common abuse of notation, we have used $|\cdot|$ for both, the $d-1$ dimensional Hausdorff measure and the $d$ dimensional Lebesgue measure. Finally, the mesh size of the tessellation $\mathcal{T}$ is defined as the maximal cell diameter,

[TABLE]

For our convergence analysis, we have to assume a certain regularity of the mesh which ensures that standard geometric constants can be chosen independently of the tessellation, and in particular, independently of the mesh size $h$. To be more specific, we require that the constants $C$ in the trace and Poincaré estimates

[TABLE]

are uniform in $K\in\mathcal{T}$ and $h>0$. Hereby, \psi_{K}:=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.7027pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.80164pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.67079pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.17912pt}}\!\int_{K}\psi\,dx is the average of $\psi$ on $K$. For a proof of the above trace and Poincaré estimates with implicit constants, we refer to [10, Chapter 4.3 and 4.5]. It is worth to note that for the particular choice $\psi\equiv 1$, the trace estimate implies the isoperimetric property

[TABLE]

which in turn guarantees that the volume of each cell $K\in\mathcal{T}$ is of order $h^{d}$, and its surface area is of order $h^{d-1}$, uniformly in $h$.

As we are concerned with an implicit scheme, the time step size ${\delta}$ can be chosen independently of the mesh size. For convenience, we assume that ${\delta}$ is fixed in each step, and thus, we can write $t^{n}=n\,{\delta}$ for the $n$-th time step. The final time step is the largest integer $N$ satisfying $N\,{\delta}\leq T$. We will sometimes write $\llbracket 0,N\rrbracket:=\{0,1,\dots,N\}$.

The upwind scheme approximates solutions by cell averages. On each cell $K\in\mathcal{T}$, the initial datum $\rho^{0}$ is thus approximated by its average

[TABLE]

Similarly, in each time interval $[t^{n},t^{n+1})$ and cell $K$, the source term is replaced by

[TABLE]

The scheme takes into account the net fluxes across the cell faces. The average normal velocity from one cell $K$ to a neighboring cell $L\sim K$ is defined by

[TABLE]

Here, $\mathcal{H}^{d-1}$ denotes the $d-1$ dimensional Hausdorff measure. Notice that these quantities are well-defined thanks to the trace theorem for Sobolev functions (see, for instance, [10, Chapter 4.3]) and assumption (4). By the sign convention of cell normals it holds $u_{KL}^{n}=-u_{LK}^{n}$. We need to distinguish between the fluxes inwards and outwards each reference cell. We thus write $u_{KL}^{n\pm}=\left(u_{KL}^{n}\right)^{\pm}$ with $(q)^{+}:=\max\{0,q\}$ and $(q)^{-}:=\max\{0,-q\}$ denoting the positive and the negative part of a quantity $q\in\mathbf{R}$, respectively.

The implicit upwind finite volume scheme for the linear continuity equation (1) now reads

[TABLE]

for every $n\in\llbracket 0,N\rrbracket$ and $K\in\mathcal{T}$. Each $\rho_{K}^{n}$ can be thus thought of as the approximate volume average of the exact solution at time $t^{n}$. The approximate solution $\rho_{{\delta},h}$ is accordingly given by

[TABLE]

In the case $n=0$, we also write $\rho_{h}^{0}=\rho_{{\delta},h}^{0}=\rho_{{\delta},h}(0,\cdot)$. The approximate source term $f_{{\delta},h}$ is defined analogously.

For later purposes it is beneficial to remark that the upwind scheme (12) can equivalently be formulated as

[TABLE]

which follows from the identities u_{KL}^{n+}=\frac{1}{2}\big{(}\lvert u_{KL}^{n}\rvert+u_{KL}^{n}\big{)} and u_{KL}^{n-}=\frac{1}{2}\big{(}\lvert u_{KL}^{n}\rvert-u_{KL}^{n}\big{)}. Well-posedness (cf. Lemma 1) and quantified stability (cf. Lemma 2) analogous to (5) follow under the additional condition that ${\delta}\leq{\delta}_{\max}$, where for some $\kappa>1$ the maximal time step size ${\delta}_{\max}={\delta}_{\max}(\kappa)$ is such that

[TABLE]

A similar condition on the time step size was introduced earlier by Boyer, see [4, Eq. (3.1)]. We introduce the constant $\kappa$ in order to quantify how close the approximate solutions $\rho_{\delta,h}$ get to satisfying the a priori estimate (5). In fact, we are able to prove a substitute for (5) satisfied by $\rho_{\delta,h}$ in Lemma 2 below, in which under condition (15) on the maximal time step size the exponent $1-\nicefrac{{1}}{{q}}$ on the compressibility constant in (5) is replaced by $\kappa(1-\nicefrac{{1}}{{q}})$. Notice that in the case of divergence-free vector fields, we can set $\kappa=1$.

3. Main results

Our main result is an estimate on the numerical error generated by the implicit upwind finite volume scheme (12) for the continuity equation (1). Before stating the result, we recall or specify the underlying hypotheses. We suppose that the initial configuration $\rho^{0}$ and the source-sink distribution $f$ are integrable functions such that

[TABLE]

for some $q\in(1,\infty)$. Here $W^{-1,1}(\Omega)$ is the space that is dual to the homogeneous Lipschitz space $\dot{W}^{1,\infty}(\Omega)$. For the advecting velocity field $u$ we suppose that

[TABLE]

where $p\in(1,\infty)$ is such that $\nicefrac{{1}}{{p}}+\nicefrac{{1}}{{q}}=1$. Slightly modifying the arguments of DiPerna and Lions [9], one can show that under these assumptions, the Cauchy problem for the continuity equation (1) is well-posed in the class of functions $\rho$ with $\rho\in L^{\infty}((0,T);L^{q}(\Omega))$. Notice that for the purpose of well-posedness one could drop the assumption that $f\in L^{\infty}((0,T);W^{-1,1}(\Omega))$. This assumption, however, is crucial for the purpose of optimal stability estimates that will enter our analysis. Likewise, in [9] the vector field $u$ can possibly be unbounded, but for our numerical analysis it is important for $u$ to be in addition uniformly bounded in time and space

[TABLE]

Let us now give our precise result.

Theorem 1**.**

Suppose that $\delta_{\max}$, $\rho^{0}$, $f$, and $u$ are given such that (15), (16), (17) and (18) hold. Let $\rho\in L^{\infty}((0,T);L^{q}(\Omega))$ denote the exact solution to the continuity equation (1) and for ${\delta}\in(0,{\delta}_{\max}\wedge 1)$ and $h\in(0,1)$ let $\rho_{{\delta},h}$ denote the approximate solution associated via (13) to the upwind scheme (12) with coefficients (9), (10) and (11). Suppose that the mesh is non-degenerate in the sense of (7). Then for any $r>0$ it holds

[TABLE]

uniformly in $t\in(0,T)$.

Before discussing the result, let us briefly comment on the notation. Here and in the following we write $a\lesssim b$ whenever there is a constant $C$ independent of $h$ and ${\delta}$ such that $a\leq Cb$. Notice that we kept the supremum norm of the velocity field and the total time in the right-hand side of (19) for the purpose of dimensional consistency: Since $r$ has the dimension of a length, the term on the right-hand side is dimension-free. Other terms dependent on the velocity fields, the data or the solution have been absorbed into the implicit constant inside “$\lesssim$”.

In our main estimate (19), $\Pi\left(\rho(t,\cdot),\rho_{{\delta},h}(t,\cdot)\right)$ is the set of all joint measures with marginals $(\rho-\rho_{{\delta},h})^{+}$ and $(\rho-\rho_{{\delta},h})^{-}$, which is non-empty because $\rho(t,\cdot)$ and $\rho_{{\delta},h}(t,\cdot)$ have same total mass, see Equation (34) on page 34 below. The quantity on the left-hand side of (19) is a Kantorovich–Rubinstein distance, that originates from the theory of optimal mass transportation. Roughly speaking, in the original context this distance describes the minimal total cost that is necessary for transferring the configuration $\rho$ into the configuration $\rho_{{\delta},h}$ if $\log(z/r+1)$ is the cost for shipping a unit volume over the distance $z$. Notice that the cost function (and thus the Kantorovich–Rubinstein distance) is singular in the limit $r\to 0$. We will give a precise definition of the marginal conditions and some important features of Kantorovich–Rubinstein distances in Section 4 below.

The result in Theorem 1 can be interpreted as follows: Since Kantorovich–Rubinstein distances metrize weak convergence (of measures) [24, Theorem 7.12], by choosing $r=\sqrt{h}+\sqrt{{\delta}}$ one deduces from (19) that

[TABLE]

as $h\to 0$ and ${\delta}\to 0$. The result shows the classical $\nicefrac{{1}}{{2}}$ convergence rate for upwind schemes on unstructured meshes found earlier in [13, 12, 23, 15, 14, 7], just that the strong norms considered in the classical setting are traded for weak convergence measures in the DiPerna–Lions setting. In case of the explicit scheme with Cartesian meshes under a CFL condition of the form ${\delta}\,\|u\|_{L^{\infty}}\leq h$, we obtained an analogous result in our previous work [18]: We proved that the rate of weak convergence is at least ${\nicefrac{{1}}{{2}}}$. A simple example moreover shows that this bound is optimal! Indeed, considering the one-dimensional setting with constant velocity, we prove that, on the one hand, no convergence rates can be obtained in strong norms: For every small $\varepsilon>0$ there exist an initial configuration such that

[TABLE]

For these data we show on the other hand that the rate of weak convergence measured in the $1$-Wasserstein distance (i.e., the Kantorovich–Rubinstein distance with linear cost function) is at most $\nicefrac{{1}}{{2}}+\varepsilon$, which almost matches the $\mathcal{O}(h^{\nicefrac{{1}}{{2}}})$ error bound—though for a different weak metric. We believe that similar calculations also apply for the implicit scheme yielding the optimality of Theorem 1. We plan to address this issue in the future. Strong convergence without rates of the implicit upwind scheme (and the more general discrete Galerkin approximation) in the setting of this paper was proved earlier by Walkington [25] and Boyer [4].

The order of the upwind scheme is formally $1$. The loss in the convergence rate from $1$ to $\nicefrac{{1}}{{2}}$ is caused by numerical diffusion, which is analytically manifested in “weak $BV$” estimates (cf. Proposition 1). These estimates have in the (heuristic) case $|u_{KL}^{n}|\sim\|u\|_{L^{\infty}}\sim U$ the form

[TABLE]

Such estimates are standard tools in the convergence analysis of numerical schemes for hyperbolic equations under classical regularity assumptions, see also [11, Chapters 5–7].

4. Transport distance with logarithmic cost function

In this section we provide the rigorous definition of the Kantorovich–Rubinstein distance appearing in our error estimate (19) and we collect those of its properties which will be relevant in the subsequent analysis. For a general introduction into the theory of optimal transportation, we refer to Villani’s monograph [24].

Given two nonnegative distributions $\rho_{1}$ and $\rho_{2}$ on $\Omega$ with same total mass, i.e.,

[TABLE]

the set $\Pi(\rho_{1},\rho_{2})$ consists of all those joint measures $\pi$ on the product space $\Omega\times\Omega$ which have the marginals $\rho_{1}$ and $\rho_{2}$, that is,

[TABLE]

This is equivalent to the requirement that

[TABLE]

where $C(\overline{\Omega})$ is the set of all functions that are continuous up to the boundary of $\Omega$. For any positive number $r$, we then define

[TABLE]

Functionals of this type were originally introduced by Kantorovich to compute the minimal cost for transferring goods from producers to consumers. In this context, $d_{r}(x,y):=\log({|x-y|}/{r}+1)$ plays the role of a cost function. The measures $\pi$ are usually referred to as transport plans. It is not difficult to see that the infimum in (22) is in fact attained, see, for instance, Theorem 1.3 in [24]. The corresponding minimizer is unique because $d$ is strictly concave [24, Theorem 2.45] and will in the sequel be denoted by $\pi_{{\mathrm{opt}}}$ and called optimal transport plan.

Instead of working with (22) directly, we will mostly consider the dual formulation

[TABLE]

cf. [24, Theorem 1.14], which admits a maximizer $\zeta_{{\mathrm{opt}}}$, the so-called Kantorovich potential, cf. [24, Exercise 2.35]. In particular, it holds that

[TABLE]

The duality formula is known as the Kantorovich–Rubinstein theorem and has a number of important implications. First, $D_{r}(\rho_{1},\rho_{2})$ is a transshipment cost which only depends on the difference $\rho_{1}-\rho_{2}$. We may therefore extend the definition of $D_{r}(\rho_{1},\rho_{2})$ to densities of the same mass (20) that are not necessarily nonnegative, so that

[TABLE]

Because $d_{r}(x,y)$ defines a distance on $\Omega$, $D_{r}(\rho_{1},\rho_{2})$ becomes a distance on the space of functions with same (finite) mass [24, Theorem 7.3]. We will accordingly refer to $D_{r}(\rho_{1},\rho_{2})$ as a Kantorovich–Rubinstein distance. An immediate consequence of this observation is the validity of the triangle inequality

[TABLE]

There is a second triangle-type inequality that we will make use of later and which immediately follows from the dual formulation (23), namely

[TABLE]

By the sublinearity of the logarithm it is $d_{r}(x,y)\leq\frac{\lvert x-y\rvert}{r}$ and therefore

[TABLE]

where the latter is the $1$-Wasserstein distance in its dual formulation.

One of the fundamental properties of Kantorovich–Rubinstein distances is of topological nature and plays a central role in the interpretation of our main result in Theorem 1: Kantorovich–Rubinstein distances metrize weak convergence (of measures) in the sense that

[TABLE]

as $k\to\infty$, for any sequence $(\rho_{k})_{k\in\mathbf{N}}$ of densities of same mass as $\rho$, cf. [24, Theorem 7.12].

We will finally mention a crucial relation between optimal transport plan $\pi_{{\mathrm{opt}}}$ and Kantorovich potential $\zeta_{{\mathrm{opt}}}$. The minimizer $\pi_{{\mathrm{opt}}}$ is concentrated on the set

[TABLE]

cf. [24, Exercise 2.37], which in turn yields a formula for the derivative of the Lipschitz function $\zeta_{{\mathrm{opt}}}$. Indeed, it holds that

[TABLE]

for $\pi_{{\mathrm{opt}}}$-a.e. $(x,y)\in\Omega\times\Omega$. The global estimate

[TABLE]

is a consequence of the characteristic Lipschitz condition in (23).

5. The upwind scheme: II. Properties and estimates

Let us start by quoting a result on existence, uniqueness, conservativity and monotonicity of the upwind scheme, for which we refer to [11].

Lemma 1**.**

Under the assumption of Theorem 1 has the implicit upwind finite volume scheme (12) a unique solution. This solution is mass preserving in the sense that

[TABLE]

Moreover, if $\rho_{h}^{0}$ and $f_{{\delta},h}$ are both nonnegative so is $\rho_{{\delta},h}$.

The next lemma provides numerical stability for the implicit upwind scheme. The numerical stability estimate (30) below is the discrete counterpart of the a priori estimate (5). In addition to that we achieve control of spatial and temporal discrete gradients—a manifestation of the numerical diffusion introduced by the scheme—see (31) below. We will see later in Proposition 1 that the latter convert into weak $BV$ estimates.

Lemma 2** (Stability and energy estimate).**

Suppose that $\rho_{h}^{0}$ and $f_{{\delta},h}$ are nonnegative. Let $\rho_{{\delta},h}$ be the solution to the implicit upwind scheme (12). Then for any $q\in(1,\infty)$, any $\kappa>1$ and any ${\delta}\leq{\delta}_{\max}(\kappa)$ as defined in (15), it holds

[TABLE]

where \Lambda_{{\delta},h}:=\exp\big{(}\lVert\left(\nabla\cdot u\right)_{{\delta},h}^{-}\rVert_{L^{1}(L^{\infty})}\big{)} with $(\nabla\cdot u)_{{\delta},h}$ defined analogously to $f_{{\delta},h}$ in (10) and (13).

In addition, for any $\bar{q}\in(1,\min\{q,2\}]$, the following spatial and temporal discrete gradient bounds hold

[TABLE]

where $C_{\bar{q}}$ is a numerical factor with $C_{\bar{q}}\to\infty$ as $\bar{q}\to 1$.

As it will become clear in the proof, the stability estimate (30) is also valid in the limiting case $q=1$. However, it is not clear to us how to extend (31) to that case. This, in fact, is the reason why we have to restrict ourselves to the setting with $q>1$.

Proof.

From the monotonicity of the scheme stated in Lemma 1 and the nonnegativity assumption on the data $\rho^{0}_{h}$ and $f_{{\delta},h}$, we deduce that the solution $\rho_{{\delta},h}$ is nonnegative, too. For the proof, it will be convenient to use the formulation (14) of the upwind scheme. Multiplication of (14) by $\lvert K\rvert$ yields

[TABLE]

Let us denote the terms in the above identity as $\operatorname{I}^{n}_{K}+\operatorname{II}^{n}_{K}+\operatorname{III}^{n}_{K}=\operatorname{IV}^{n}_{K}$. Our derivation of (30) mimics the one of (5) in the continuous setting. First, we test the equation with $(\rho_{K}^{n+1})^{q-1}$ and sum over $K$. For the first term $\operatorname{I}^{n}_{K}$, we use the Hölder inequality to obtain

[TABLE]

Next, by recalling that $u_{KL}^{n}=-u_{LK}^{n}$ we symmetrize the second term $\operatorname{II}^{n}_{K}$

[TABLE]

Let us introduce the $q$-mean

[TABLE]

and note that $\theta_{2}(a,b)$ is the arithmetic mean (see Appendix A for some of its properties). By the definition of the $q$-mean, we have the identity

[TABLE]

In particular, $\operatorname{II}^{n}$ can be decomposed into the sum $\operatorname{II}^{n}_{1}+\operatorname{II}^{n}_{2}$, where

[TABLE]

To estimate the term $\operatorname{II}^{n}_{1}$, we do another symmetrization and use the divergence theorem in every cell $K$, so that

[TABLE]

where we have set $\lambda^{n}:={\delta}\,\|(\nabla\cdot u_{{\delta},h}(t^{n}))^{-}\|_{L^{\infty}}$ for abbreviation. Let us now estimate the remainder term $\operatorname{II}^{n}_{2}$. In Appendix A, we derive the following elementary estimate between $\theta_{q}$ and $\theta_{2}$:

[TABLE]

Applying this estimate inside of $\operatorname{II}^{n}_{2}$, we arrive at

[TABLE]

Likewise, summation in $K$ and symmetrization leads to a similar bound on $\operatorname{III}^{n}_{K}$, namely

[TABLE]

Finally, the term obtained from $\operatorname{IV}^{n}_{K}$ after testing by $(\rho_{K}^{n+1})^{q-1}$ and summation in $K$ is estimated by the Hölder inequality as

[TABLE]

A combination of all the estimates so far gives

[TABLE]

where $c_{q}$ is obtained as

[TABLE]

Dropping for the moment the second term on the left hand side of (32) and dividing by $\lVert\rho_{{\delta},h}^{n+1}\rVert_{L^{q}}^{q-1}$ gives the bound

[TABLE]

From the choice of ${\delta}_{\max}(\kappa)$ in (15) it follows $\,\tfrac{q-1}{q}\lambda^{n}\leq\tfrac{\kappa-1}{\kappa}$ and therefore

[TABLE]

An iteration of (33) thus generates the bound (30).

In order to establish the energy estimates, we reconsider term $\operatorname{I}^{n}$. By the convexity of the map $a\mapsto a^{\bar{q}}$, we notice that for $\bar{q}\in(1,2]$ it holds

[TABLE]

Applying this estimate with $a=\rho_{K}^{n+1}$ and $b=\rho_{K}^{n}$ yields

[TABLE]

If we use this bound on $\operatorname{I}^{n}$, we obtain instead of (32) that

[TABLE]

It is furthermore convenient to rewrite the third term on the left-hand side by applying the following elementary inequality for any $\bar{q}\in(1,2]$ and any $a,b>0$

[TABLE]

with $a=\rho_{K}^{n+1}$ and $b=\rho_{L}^{n+1}$. Therewith, summation over $n$, then yields

[TABLE]

Applying the just proven stability estimate (30) yields (31). ∎

Before stating and proving the afore mentioned weak $BV$ estimates for the upwind scheme, we note two obvious relations between exact and approximated data. From the definitions in Section 2, it immediately follows that the initial total masses are identical in the continuous and the discrete models:

[TABLE]

Similarly, as for the continuous problem, the approximate net source is vanishing at any time $t\in(0,T)$, because with $n\in\llbracket 0,N-1\rrbracket$ such that $t\in[t^{n},t^{n+1})$ it holds

[TABLE]

We will see in Lemma 1 below that these two facts together entail for any $t\in(0,T)$ that

[TABLE]

Let us now prove some basic estimates between discretized and continuous versions of various quantities.

Lemma 3**.**

The following estimates hold: $\|f_{{\delta},h}\|_{L^{1}(L^{q})}\leq\|f\|_{L^{1}(L^{q})}$, $\|\rho_{h}^{0}\|_{L^{q}}\leq\|\rho^{0}\|_{L^{q}}$ and $\lVert(\nabla\cdot u)_{{\delta},h}^{-}\rVert_{L^{1}(L^{\infty})}\leq\lVert\left(\nabla\cdot u\right)^{-}\rVert_{L^{1}(L^{\infty})}.$

Proof.

The estimate for the initial data is a straight-forward consequence of Jensen’s inequality. Indeed,

[TABLE]

By essentially the same reasoning, for any source term of the form $f(t,x)=f^{1}(t)f^{2}(x)$, it holds that

[TABLE]

if by subscript ${\delta}$ and $h$ we denote the discretization in time and space, respectively. It remains to conclude with an approximation argument: Thanks to the density of smooth functions in $L^{1}(L^{q})$, it is enough to prove the statement for continuous functions. Moreover, the Stone–Weierstraß theorem enables us to furthermore approximate continuous functions uniformly by functions of the form $f^{1}(t)f^{2}(x)$—for which the estimate is shown above. The proof for $\nabla\cdot u$ follows from a combination of the Jensen inequality applied to the convex function $a\mapsto(a)^{-}$ with the argument above. ∎

As a consequence, under the assumptions of Theorem 1, the expressions in the right-hand sides of (30) and (31) are both $\mathcal{O}(1)$, for instance,

[TABLE]

Let us now establish the weak $BV$ estimates in space and time, which will occur in later estimates and are a manifestation of the numerical diffusion.

Proposition 1** (Weak $BV$ estimates).**

Suppose the assumptions of Theorem 1 hold and let $\rho_{{\delta},h}$ be the solution to the implicit upwind scheme (12). Then, it holds

[TABLE]

Notice that these estimates do not have any counterparts in the continuous model (1). In particular, no compactness estimates can be inferred.

Proof.

The solution to the upwind scheme can be split into $\rho_{{\delta},h}=(\rho_{{\delta},h})_{+}-(\rho_{{\delta},h})_{-}$ where $(\rho_{{\delta},h})_{\pm}$ is the nonnegative discrete solution with data $(\rho_{h}^{0})^{\pm}$ and $(f_{{\delta},h})^{\pm}$ (first discretized then decomposed). Hence, once we have established the estimates (37) and (36) for nonnegative data, the estimate follows for general data just by the triangle inequality and the above observation. Therefore, let $\rho_{{\delta},h}$ be the nonnegative solution of the upwind scheme to the nonnegative data $\rho_{h}^{0}$ and $f_{{\delta},h}$. The term (37) is estimated by applying the Cauchy-Schwarz inequality, whereby we smuggle in an additional weight $\left((\rho_{K}^{n+1}+\rho_{L}^{n+1})/2\right)^{\bar{q}-2}$ for some $\bar{q}\in(1,\min\{q,2\}]$:

[TABLE]

After summing over $n$ and another Hölder inequality in time, the term $\operatorname{I}^{n}$ can be directly estimated by the energy estimate (31) from Lemma 2 for $\bar{q}\in(1,\min\{q,2\}]$. For the term $\operatorname{II}^{n}$, we observe that $\left((a+b)/2\right)^{2-\bar{q}}\leq a^{2-\bar{q}}+b^{2-\bar{q}}$ for any $a,b>0$ and use the regularity assumption on the mesh (8) to estimate

[TABLE]

Now, we choose $\bar{q}=q$ if $q\in(1,2)$ and obtain after summation over $n$:

[TABLE]

Since $q\geq 1$, we have that $2-q<q$ and can estimate $\|\rho_{{\delta},h}\|_{L^{\infty}(L^{2-q})}$ by $\|\rho_{{\delta},h}\|_{L^{\infty}(L^{q})}$ up to a factor depending on $\lvert\Omega\rvert$. In the case $q\geq 2$, we choose $\bar{q}=2$ in (31) and deduce the analogous result. Notice that the solution is uniformly bounded in any $L^{\infty}(L^{\bar{q}})$ with $\bar{q}\leq q$ thanks to the stability estimate (30) and the bounds from Lemma 3. A combination of the previous estimates yields (37).

The estimate (36) follows along the same lines by first applying the Cauchy-Schwarz inequality with the same weight and then the a priori estimates (30) and (31). ∎

6. Error estimates and proof of Theorem 1

In this section, we present the proof of our main result Theorem 1. We will see that there are two classes of discretization errors contributing to estimate (19). The first class consists of errors introduced by the discretization of time and space, and thus by the corresponding finite volume approximations of data and coefficients. These errors are $\mathcal{O}(h+{\delta})$. The second class is caused by the discretization of the scheme, which is usually referred to as the truncation error. These latter errors are all $\mathcal{O}(h^{\nicefrac{{1}}{{2}}}+{\delta}^{\nicefrac{{1}}{{2}}})$, and are related to the phenomenon of numerical diffusion, see, for instance, Section 2.4 in [18].

6.1. Discretization of data

We begin the error analysis by addressing the various errors caused by the discretization of time and space. The first one concerns the discretization of the time steps for the continuous problem.

Lemma 4** (Discretization of time).**

For $n\in\llbracket 0,N-1\rrbracket$ and any $t\in[t^{n},t^{n+1})$ it holds

[TABLE]

The next error is caused by the spatial discretization of the data. Here $f_{h}$ is defined analogously to $\rho_{h}^{0}$, thus still continuous in time. The statement of this lemma is an immediate consequence of the stability estimates for continuity equations established in [19]. The argument for the initial data is already given in [18, Lemma 8].

Lemma 5** (Spatial discretization of data).**

Let $\rho^{h}$ be the solution of the continuity equation (1) with initial datum $\rho^{0}_{h}$ and source-sink distribution $f_{h}$. Then it holds for any $t\in[0,T]$ that

[TABLE]

Thirdly, we consider the error caused by discretizing data and coefficients in the time variable. Here the subscript ${\delta}$ refers to discretization by averaging over \bigl{[}t^{n},t^{n+1}\bigr{)}.

Lemma 6** (Temporal discretization of data).**

Let $\rho^{{\delta}}$ be the solution of the continuity equation (1) with driving vector field $u_{{\delta}}$ and source-sink distribution $f_{{\delta}}$. Then it holds for any $\ell\in\llbracket 0,N\rrbracket$ that

[TABLE]

The argument for this last error estimate substantially differs from those for Lemmas 4 and 5 (and those from [20, 18]). In fact, to control the errors caused by the discretization of time we have to subtly change from Eulerian to Lagrangian coordinates.

We postpone the proofs of Lemmas 4 to 6 until Subsection 6.3 below. In view of the triangle inequality (25), these first results imply that for any $t\in(0,T)$ it holds

[TABLE]

where $\ell\in\llbracket 0,N\rrbracket$ is such that $t\in\left[t^{\ell},t^{\ell+1}\right)$ and $\rho^{\delta,h}$ is the unique distributional solution to the continuity equation (1) with initial datum $\rho_{h}^{0}$, source-sink distribution $f_{\delta,h}$ and velocity $u_{\delta}$. The remaining error $D_{r}(\rho^{\delta,h}(t^{\ell}),\rho_{\delta,h}(t^{\ell}))$ governs the convergence rate. Its treatment will be illustrated in what follows.

6.2. The error caused by the scheme

By the virtue of eliminating the discretization errors in the previous subsection, we tacitly assume from here on that the continuous problem is solved with data $\rho_{h}^{0}$, $f_{{\delta},h}$ and $u_{{\delta}}$ and that $t=t^{\ell}$ for some $\ell\in\llbracket 0,N\rrbracket$, i.e., we assume that $\rho(t)=\rho^{\delta,h}(t^{\ell})$. We will see that the truncation error caused by the scheme is $\mathcal{O}(\sqrt{h}+\sqrt{{\delta}})$. Our goal in the following is to portray the main steps in the estimate of the discrete rate of change of the distance $D_{r}(\rho(t^{n}),\rho_{\delta,h}(t^{n}))$. In fact, it turns out that instead of analyzing $\rho_{\delta,h}$, it is more convenient to study the piecewise linear temporal approximate solution defined by

[TABLE]

Since $\widehat{\rho}_{{\delta},h}(t^{n})=\rho_{{\delta},h}(t^{n})$ for any $n\in\llbracket 0,N\rrbracket$, there is no additional error term to consider by replacing the piecewise constant (in time) approximation by the piecewise linear approximation. The advantage of considering $\widehat{\rho}_{{\delta},h}$ instead of $\rho_{{\delta},h}$ is that the former is (weakly) differentiable. Indeed, in a first step, we formally compute the rate of change of the Kantorovich–Rubinstein distance between $\rho$ and $\widehat{\rho}_{{\delta},h}$:

[TABLE]

Here, $\zeta_{{\mathrm{opt}}}=\zeta_{{\mathrm{opt}}}(t)$ denotes the Kantorovich potential corresponding to $D_{r}(\rho,\widehat{\rho}_{{\delta},h})$ at time $t$. By construction, it holds $\partial_{t}\widehat{\rho}_{{\delta},h}={\delta}^{-1}(\rho_{h}^{n+1}-\rho_{h}^{n})$ where we have set $\rho_{h}^{n}=\rho_{h}(t^{n})$ for any $n$. The time derivative of $\rho$ does in general not exist. However, using the continuity equation and formally integrating by parts, the above formula may be rewritten as

[TABLE]

In fact, arguing as in [19, Lemma 1], this identity can be established rigorously. We now use formulation (14) of the upwind scheme and notice that the source term drops out because $f=f_{{\delta},h}$ by the virtue of Lemma 5. After integration over \bigl{[}t^{n},t^{n+1}\bigr{)}, we thus find

[TABLE]

where

[TABLE]

with (\zeta_{{\mathrm{opt}}})_{K}^{n}:=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.7027pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.80164pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.67079pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.17912pt}}\!\int_{t^{n}}^{t^{n+1}}\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.7027pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.80164pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.67079pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.17912pt}}\!\int_{K}\zeta_{{\mathrm{opt}}}\,dx\,dt.

Up to a shift in the time variable, the first error term $\operatorname{I}^{n}$ can be controlled by the techniques developed in [19] to establish stability estimates for continuity equations. The time shift can then be compensated with the help of the temporal weak $BV$ estimate (36) from Proposition 1.

Lemma 7** (Estimate of $\operatorname{I}^{n}$).**

[TABLE]

The second error term $\operatorname{II}^{n}$ has to be rewritten using an integration by parts. An important ingredient in the following bound is the spatial weak $BV$ estimate (37) from Proposition 1.

Lemma 8** (Estimate of $\operatorname{II}^{n}$).**

[TABLE]

The last error term $\operatorname{III}^{n}$ is of purely diffusive origin and is controlled again thanks to (37) from Proposition 1.

Lemma 9** (Estimate of $\operatorname{III}^{n}$).**

[TABLE]

The proof of Theorem 1 now follows by a combination of the Lemmas 4–9. Indeed, summing over $n\in\llbracket 0,\ell-1\rrbracket$ in (39), using the fact that $\rho(0)=\rho^{{\delta},h}(0)=\rho_{{\delta},h}(0)$ by the hypothesis of this subsection, and using Lemmas 7–9, we find that

[TABLE]

if $h$ is sufficiently small. Plugging this estimate into (38) then yields the statement in Theorem 1.

6.3. The proofs of the Lemmas 4–9

To simplify the notation in the following, we occasionally write $\psi_{t}$ for $\psi(t,\cdot)$ for a given function $\psi\in L^{1}((0,T);X)$.

We start with the short proof of Lemma 4.

Proof of Lemma 4.

We choose $n\in\llbracket 0,N\rrbracket$ such that $t\in[t^{n},t^{n+1})$ and write $\zeta_{t}=\zeta_{{\mathrm{opt}}}(t,\cdot)$ if $\zeta_{{\mathrm{opt}}}$ is the Kantorovich potential corresponding to $D_{r}(\rho_{t},\rho_{t^{n}})$, so that $D_{r}(\rho_{t},\rho_{t^{n}})=\int_{\Omega}\zeta_{t}(\rho_{t}-\rho_{t^{n}})\,dx$. Using the distributional formulation of (1), we then find

[TABLE]

In view of (29), the first term on the right-hand side is easily controlled by

[TABLE]

For the second term we have by using (27)

[TABLE]

A combination of the previous estimates yields the statement of the lemma. ∎

The proofs of Lemmas 5–9 have the flavor of stability estimates for continuity equations as those that were recently established in [19, 20]. Their main ingredients are new estimates on the rate of change of the Kantorovich–Rubinstein distance in which Crippa–De Lellis-type inequalities play a central role. Such inequalities, first derived in [6] and adapted to the context of Kantorovich–Rubinstein distances in [5] and [16], provide a way of bounding integrals of difference quotients by $L^{p}$ norms of gradients. At the heart of their proofs is a standard tool from harmonic analysis: the maximal function. The maximal function of a given function $v$ on $\Omega$ is defined by

[TABLE]

One of the fundamental estimates in the theory of singular integrals is the $L^{p}$ bound on the maximal function,

[TABLE]

which is valid for any $p\in(1,\infty]$. In the sequel, we will apply this estimate to the velocity gradient $\nabla u$. The fact that (40) fails in the case $p=1$ is the reason why or theory does not extend to the setting where $u\in L^{1}(W^{1,1})$. Moreover, we need a pointwise estimate which is of Morrey-type, namely

[TABLE]

Hereby, $\bar{v}$ denotes a Sobolev extension of $v$ on $\mathbf{R}^{d}$. The fundamental inequality for maximal functions (40) can be found in any standard reference on harmonic analysis, see, for instance, [21, p. 5, Theorem 1]. The Morrey-type estimate (41) is rather elementary. Its proof is, for instance, contained in [10, p. 143, Theorem 3]. The Sobolev extension theorem is given in [10, p. 135, Theorem 1].

Proof of Lemma 5.

For notational convenience in the proof, we omit the time dependence of the occurring functions, if the context allows it. We first notice that the rate of change of the Kantorovich distance between the two solutions $\rho$ and $\rho^{h}$ which are both advected by the same velocity field takes the form

[TABLE]

The proof for the homogeneous continuity equation is contained in [19] and is easily adapted to the case that includes sources and sinks. The second term on the right-hand side is estimated by $D_{r}(f,f_{h})$ due to the dual variational formula (23). For the first term we claim that

[TABLE]

We postpone the proof of (42) until later and continue with the estimate for $D_{r}\left(\rho,\rho^{h}\right)$. Integration in time and the a priori estimate (5) imply that

[TABLE]

for some $C<\infty$. The first two terms on the right-hand side are controlled in a similar way. It is thus enough to focus on one of them, say $D(\rho^{0},\rho_{h}^{0})$. The estimate of this term was already given in [18, Lemma 8]. For the convenience of the reader, we provide the proof here again; this time, however, we present an argument based on (23). If $\zeta_{{\mathrm{opt}}}$ denotes the associated Kantorovich potential, it holds

[TABLE]

where in the second equality we have used the symmetry of the averaging operator $(\cdot)_{h}$ defined by averaging over each cell $K\in\mathcal{T}$. For every $x\in K$ we use the Lipschitz property of $\zeta_{{\mathrm{opt}}}$, cf. (23), in the sense that

[TABLE]

It immediately follows that $D(\rho^{0},\rho_{h}^{0})\leq{h}/{r}\|\rho^{0}\|_{L^{1}}$.

It remains to prove (42). We use the marginal condition (21) and the formula (28) for the gradient of the Kantorovich potential to rewrite

[TABLE]

Here $\pi_{{\mathrm{opt}}}$ is the optimal transport plan that corresponds to the marginals $(\rho-\rho^{h})^{+}$ and $(\rho-\rho_{h})^{-}$, cf. (24). At this point, we need an estimate of Crippa–De Lellis-type to control the expression on the right-hand side by $\|\nabla u\|_{L^{p}}$. Dropping $r$ and making use of the Morrey-type inequality (41), we find that

[TABLE]

where $\bar{u}$ denotes a Sobolev extension of $u$ to all of $\mathbf{R}^{d}$. In view of the marginal condition (21), the term on the right-hand side becomes

[TABLE]

An application of Hölder’s inequality and of the fundamental estimate (42) thus yield

[TABLE]

The statement in (42) now follows from the continuity of the extension operator, the a priori estimate (5) for the continuity equation and the Young-type estimates in Lemma 3. ∎

The proof of Lemma 6 relies in addition on the superposition principle, which represents the solution of the continuity equation (1) in terms of a Lagrangian flow, and some properties of this flow. We will see that it is enough to consider the classical setting with spatially smooth velocity fields. For further references and the convenience of the reader, we recall its definition and basic properties. For $0\leq s\leq t\leq T$ a mapping $\Phi_{t,s}:\Omega\to\Omega$ is called Lagrangian flow if for every $x\in\Omega$ the mapping $t\mapsto\Phi_{t,s}(x)$ is an (integral) solution to the ordinary differential equation

[TABLE]

The no-flux boundary condition (2) ensures that the flow $\Phi_{t,s}$ maps $\Omega$ to $\Omega$. By the compressibility assumption (4) (which can be retained under approximation) it follows that the Jacobian given by

[TABLE]

is bounded below by $\Lambda^{-1}$. In particular, it holds that

[TABLE]

The constant $\Lambda$ is accordingly referred to as the compressibility constant of $\Phi$. It is moreover well-known that the family $\{\Phi_{t,s}\}_{0\leq s\leq t\leq T}$ is an inhomogeneous semigroup, that is

[TABLE]

The theory for Lagrangian flows has been extended to the low regularity framework by DiPerna and Lions [9]. See also the recent contribution of Ambrosio, Colombo and Figalli [3] for the theory maximal flows on bounded domains.

The superposition principle for a solution $\rho$ to the inhomogeneous continuity equation (1) reads

[TABLE]

where we have set $\phi_{t}:=\Phi_{t,0}$ and where the push forward of a map $\phi:\Omega\to\Omega$ is defined for any Borel function $\xi$ and Borel measure $\mu$ by

[TABLE]

We now proceed with the proof of Lemma 6.

Proof of Lemma 6.

In view of the stability results obtained in [9] or [19] and by a standard approximation argument, it is enough to consider vector fields that are smooth in the spatial variable.

Let $\Phi^{{\delta}}$ be the Lagrangian flow for the vector field $u_{{\delta}}$. Based on the superposition principle (45), we have the following representations

[TABLE]

Let us introduce some convenient abbreviations for the homogeneous and inhomogeneous components of the solutions $\rho$ and $\rho^{\delta}$, respectively

[TABLE]

where for $s\in[t^{n},t^{n+1})$, we set $I^{\delta}_{s}:=[t^{n},t^{n+1})$. These definitions imply $\rho_{t}=\rho_{t}^{\hom}+F_{t}$ and likewise $\rho_{t}^{\delta}=\rho_{t}^{{\delta},\hom}+F_{t}^{{\delta}}$. Now, by using the triangle inequality (26) for the Kantorovich–Rubinstein distance, we estimate

[TABLE]

By the boundedness assumption (18) on the vector field and the definition of the flows (43), it follows for any $s\in[0,T]$ and $t_{0},t_{1}\in[s,T]$ the basic estimate

[TABLE]

We first turn to the estimate of D_{r}\big{(}\rho_{t^{\ell}}^{\hom},\rho_{t^{\ell}}^{{\delta},\hom}\big{)}. Denoting by $\zeta_{{\mathrm{opt}}}$ the corresponding optimal Kantorovich–Rubinstein potential, we have the estimate

[TABLE]

where in the last inequality we have used the Lipschitz property of $\zeta_{{\mathrm{opt}}}$ in (23). It is thus enough to control the quantity on the right-hand side. Using the triangle inequality in $\mathbf{R}^{d}$, the concavity of the logarithm, and the definition of the flows, we first estimate

[TABLE]

for any $n\in\llbracket 0,N\rrbracket$. The numerator in the integrand is controlled as follows: In view of the definition of $u_{{\delta}}$ and the Morrey-type estimate (41), we have that

[TABLE]

Using the estimate (48), we get $|\phi_{\sigma}-\phi_{\tau}^{\delta}|\leq|\phi_{t^{n}}-\phi_{t^{n}}^{\delta}|+2{\delta}\|u\|_{L^{\infty}}$. Combining the previous estimates thus yields

[TABLE]

It remains to notice that by recalling (46) and applying Hölder’s inequality it follows

[TABLE]

The terms on the right-hand side are estimated with the help of the fundamental inequality for maximal functions (40) and the compressibility property (44). The second term on the right-hand side of (49) is treated similarly by using in addition Lemma 3. Hence, summing over $n\in\llbracket 0,\ell-1\rrbracket$ yields

[TABLE]

Let us now focus on the term $D_{r}(F_{t^{\ell}},F_{t^{\ell}}^{\delta})$. If $\zeta_{{\mathrm{opt}}}$ is the Kantorovich potential corresponding to this transport distance, then we obtain analogously to the above computation

[TABLE]

The estimation of the expression on the right proceeds very analogously to the estimations above. We will thus only sketch it. In a first step, we obtain, similarly as in (49) that

[TABLE]

for any $n\in\llbracket 0,N\rrbracket$ with $s\leq t^{n}$. By using the same arguments as above and summing over all $n\in\llbracket k,\ell-1\rrbracket$, where $k$ is such that $s,\sigma\in\left[t^{k-1},k\right)$, we find that

[TABLE]

We observe that thanks to (48), we have the estimate

[TABLE]

Plugging this bound into (50) and integrating over $\sigma$ and $s$ thus yields

[TABLE]

which is what we aimed to show. ∎

Proof of Lemma 7.

We split the term $\operatorname{I}^{n}$ into the sum $\operatorname{I}_{1}^{n}+\operatorname{I}^{n}_{2}$ with

[TABLE]

We claim that

[TABLE]

The proof of the first estimate in (51) is a direct consequence of estimate (42) that was established in the proof of Lemma 5. (Here, it is important to notice that $\zeta_{{\mathrm{opt}}}$ is a Kantorovich potential associated with $D_{r}(\rho,\widehat{\rho}_{{\delta},h})$.) In fact, applying (42) together with the triangle inequality for $\|\cdot\|_{L^{q}}$ yields

[TABLE]

Summing over $n$ and invoking the a priori estimates (5) and (35) gives the result.

For the second statement in (51), we notice that because of

[TABLE]

we have

[TABLE]

Using the Lipschitz bound for $\zeta_{{\mathrm{opt}}}$ (29) and summing over $n$ gives

[TABLE]

The conclusion is an immediate consequence of the temporal weak $BV$ estimate (36) from Proposition 1. ∎

Proof of Lemma 8.

For the proof, we write $\zeta:=\zeta_{{\mathrm{opt}}}$ and $u^{n}:=u(t^{n})$. We first notice that the anti-symmetry of $u_{KL}^{n}$ implies that

[TABLE]

Moreover, the divergence theorem yields

[TABLE]

Regarding the other term in $\operatorname{II}^{n}$, we have by integration by parts on each cell $K$ that

[TABLE]

It thus follows that $\operatorname{II}^{n}$ can be rewritten as

[TABLE]

where $\zeta^{n}=\zeta^{n}(x)$ denotes the time average of $\zeta$ over $\left(t^{n},t^{n+1}\right)$. We can furthermore decompose this expression

[TABLE]

and we set $\operatorname{II}^{n}=:\operatorname{II}^{n}_{1}+\operatorname{II}^{n}_{2}+\operatorname{II}^{n}_{3}$, accordingly. To estimate $\operatorname{II}_{1}^{n}$, we notice that we can smuggle in the constant function $\zeta_{K}^{n}$ in each boundary integral, leading to

[TABLE]

Thanks to the Lipschitz property of $\zeta_{{\mathrm{opt}}}$, cf. (29), it holds that $|\zeta^{n}-\zeta_{K}^{n}|\leq h/r$ uniformly in $K$. We now combine the trace and Poincaré estimate (7), to the effect that

[TABLE]

We are thus left with

[TABLE]

Summing over $n$ and using Hölder’s inequality and the a priori bound (35) yields that

[TABLE]

To estimate $\operatorname{II}^{n}_{2}$, we first notice that for symmetry reasons, we have

[TABLE]

Using again the Lipschitz property for $\zeta_{{\mathrm{opt}}}$, cf. (29), the latter is bounded as follows:

[TABLE]

We apply the spatial weak $BV$ estimate (37) from Proposition 1 showing

[TABLE]

It remains to investigate $\operatorname{II}^{n}_{3}$. Once again, we use the Lipschitz bound of $\zeta_{{\mathrm{opt}}}$ in (29), Hölder’s inequality and bound $\nabla\cdot u$ by $\nabla u$ up to a dimension dependent constant to estimate

[TABLE]

Summing over $n$ and using the a priori estimate (30) yields

[TABLE]

This concludes the proof of Lemma 8. ∎

Proof of Lemma 9.

From the Lipschitz bound (29) we deduce for any two neighboring cells $K$ and $L$ that

[TABLE]

Therewith, we estimate after doing a symmetrization of the sum:

[TABLE]

The proof concludes by summing over $n\in\llbracket 0,N-1\rrbracket$ and applying the spatial weak $BV$ estimate (37) from Proposition 1. ∎

Appendix A The $q$-mean

We briefly describe some helpful estimates for the $q$-mean defined for $q>1$ by

[TABLE]

1. (i)

The function $\theta_{q}$ has the following integral representation

[TABLE] 2. (ii)

The function $\theta_{q}$ is $1$-homogeneous: for any $c>0$ it holds $\theta_{q}(c\,a,c\,b)=c\,\theta_{q}(a,b)$. 3. (iii)

For any positive numbers $a\neq b$ is $q\mapsto\theta_{q}(a,b)$ strictly increasing. 4. (iv)

The function $(a,b)\mapsto\theta_{q}(a,b)$ is concave for $q\in(1,2)$ and convex for $q\in(2,\infty)$. 5. (v)

For any $a,b>0$ it holds

[TABLE]

Proof.

For the identity (i) let $t_{q}(a,b;s):=\left((1-s)a^{q-1}+sb^{q-1}\right)^{\frac{1}{q-1}}$ denote the integrand on the right hand side. Then, by a straightforward calculation it follows $\partial_{s}t_{q}(a,b;s)=\frac{b^{q-1}-a^{q-1}}{q-1}t(s)^{2-q}$ and hence $T_{q}(a,b;s):=\frac{q-1}{q}\frac{t(s)^{q}}{b^{q-1}-a^{q-1}}$ is its primitive from which (i) follows.

The $1$-homogeneity as stated in (ii) follows immediately from the definition.

For proving (iii), we note that $t_{q}(a,b;s)$ is the $\ell^{q-1}$ norm on the probability space \big{(}\{0,1\},(1-s)\delta_{0}+s\delta_{1}\big{)} for a function $f:\{0,1\}\to\mathbf{R}$ taking values $f(0)=a$ and $f(1)=b$. Hence, the statement is a consequence of the ordering of the $\ell^{p}$ spaces, which extends to any value $p\in\mathbf{R}$.

The property (iv) follows by calculating the Hessian of $(a,b)\mapsto t_{s}(a,b):$

[TABLE]

Now, one immediately recovers that $(a,b)\mapsto t_{q}(a,b;s)$ is negative semidefinite for $q\in[0,2)$ and positive semidefinite for $q>2$.

For (v), we can assume by symmetry and $1$-homogeneity that $a\in(0,1)$ and $b=1$. From (iv), we have that the mapping $(0,1)\ni a\mapsto\theta_{q}(a,1)$ is concave for $q\in(1,2)$ and convex for $q>2$. Let us first assume $q\in(1,2)$, then by convexity of $a\mapsto\theta_{2}(a,1)-\theta_{q}(a,1)$, we estimate using the secant inequality between the points [math] and $1$:

[TABLE]

For $q>2$, we apply the same argument to the convex function $a\mapsto\theta_{q}(a,1)-\theta_{2}(a,1)$. ∎