Lagrangian discretization of crowd motion and linear diffusion

TL;DR

This paper introduces a Lagrangian discretization method for modeling crowd motion and linear diffusion, using optimal transport to enforce density constraints and demonstrating convergence and numerical implementation in 2D.

Contribution

It presents a novel Lagrangian discretization scheme for crowd motion and diffusion equations, with convergence proofs and practical 2D numerical examples.

Findings

Convergence of discrete measures to continuous PDE solutions in 1D.

Feasibility of 2D numerical implementation.

Effective enforcement of density constraints via optimal transport.

Abstract

We study a model of crowd motion following a gradient vector field, with possibly additional interaction terms such as attraction/repulsion, and we present a numerical scheme for its solution through a Lagrangian discretization. The density constraint of the resulting particles is enforced by means of a partial optimal transport problem at each time step. We prove the convergence of the discrete measures to a solution of the continuous PDE describing the crowd motion in dimension one. In a second part, we show how a similar approach can be used to construct a Lagrangian discretization of a linear advection-diffusion equation, interpreted as a gradient flow in Wasserstein space. We provide also a numerical implementation in 2D to demonstrate the feasibility of the computations.

Figures

Tables

Table 1. Table 1: Error err h subscript err ℎ \mathrm{err}_{h} between the exact and numeric solution to the crowd motion model as a function of the space-discretization h ℎ h .

$h$	$\frac{1}{20}$	$\frac{1}{30}$	$\frac{}{40}$	$\frac{1}{50}$	$\frac{1}{100}$	$\frac{1}{200}$
${\mathrm{err}}_h$	$5.24\cdot {10}^{-2}$	$3.06\cdot {10}^{-2}$	$2.15\cdot {10}^{-2}$	$1.70\cdot {10}^{-2}$	$4.96\cdot {10}^{-3}$	$2.80\cdot {10}^{-3}$

Equations

$${\partial }_t\rho -\Delta \rho -\mathrm{div}(\rho \nabla V)=0$$

$${\partial }_t\rho -\Delta {\rho }^m-\mathrm{div}(\rho \nabla V)=0,m\in (1,+\infty )$$

$$\{\begin{array}{l}{\partial }_t\rho +\mathrm{div}(\rho v)=0\\ v=-\nabla p-\nabla V\\ 0\le \rho \le 1,\text{leavevmode}p\ge 0,\text{leavevmode}p(1-\rho )=0.\end{array}$$

$$F(\rho )=\int \rho \log \rho$$

$$F(\rho )=\{\begin{array}{ll}0& if\rho \le 1\\ +\infty & if\; not\end{array}$$

$$F_{\epsilon }(\mu ):=\underset{\rho }{\inf }F(\rho )+\frac{1}{2\epsilon }{W}_2^2(\mu ,\rho ).$$

$$F_{\epsilon }(x_1,\dots ,x_N)=F_{\epsilon }\left(\frac{1}{N}\sum _{1\le i\le N}{\delta }_{x_i}\right)$$

$$\{\begin{array}{l}\frac{1}{N}{\dot{x}}_i(t)=-{\nabla }_{x_i}{F}_{{\epsilon }_N}(x_1(t),\dots ,x_N(t))-\frac{1}{N}\nabla V(x_i(t)),\\ X^N(0)=X_0^N,\end{array}$$

$$X_0^N\in \underset{X\in {\mathbb{R}}^{Nd}}{\mathrm{argmin}}{W}_2^2({\rho }_0,{\mu }_X),$$

$$K=\{\rho \in {\mathcal{P}}_{\text{ac}}(\Omega ):\rho \le 1\},$$

$$\{\begin{array}{ll}{\partial }_t\rho +\mathrm{div}(\rho v)=0,& \rho \in K,\\ v={\Pi }_{T_{\rho }K}(-\nabla V),& \end{array}$$

$$T_{\rho }K=\left\{v\in L^2(\rho ;{\mathbb{R}}^d):\mathrm{div}(v)\ge 0\text{ on }\{\rho =1\},v\cdot n\le 0\text{ on }\partial \Omega \cap \{\rho =1\}\right\}.$$

$$\{\begin{array}{ll}{\partial }_t\rho +\mathrm{div}(\rho v)=0,& \rho \in K,\\ v=-\nabla V-w,& w\in N_{\rho }K,\end{array}$$

$$N_{\rho }K=\left\{\nabla p:p\in H^1(\Omega ),p\ge 0,p(1-\rho )=0\right\}.$$

$$\{\begin{array}{ll}{\partial }_t\rho +\mathrm{div}(\rho v)=0,& \rho \in K,\\ v=-\nabla V-\nabla p,& \\ p\ge 0,p(1-\rho )=0,& \end{array}$$

$$\{\begin{array}{l}{\partial }_t\rho \in -\partial E(\rho ),\\ \rho (0)={\rho }_0,\end{array}$$

$$E(\rho )=\int V\mathrm{d}\rho +F(\rho ),whereF(\rho )=i_K(\rho )=\{\begin{array}{ll}0& if\rho \in K\\ +\infty & if\; not.\end{array}$$

$${\mathcal{P}}_N({\mathbb{R}}^d)=\left\{\frac{1}{N}\sum _{i=1}^N{\delta }_{x_i}:x_i\in {\mathbb{R}}^d\right\}.$$

$$E_{\epsilon }(\rho )={\int }_{\Omega }V\mathrm{d}\rho +\frac{1}{2\epsilon }{\mathsf{d}}_K^2(\rho )={\int }_{\Omega }V\mathrm{d}\rho +\frac{1}{2\epsilon }\underset{\sigma \in K}{\min }{W}_2^2(\rho ,\sigma ).$$

$$E_{\epsilon }(X)=E_{\epsilon }({\mu }_X)=\frac{1}{N}\sum _{i=1}^{N}V(x_i)+\frac{1}{2\epsilon }{\mathsf{d}}_K^2({\mu }_X).$$

$${\mathcal{D}}_N=\{(x_1,\dots ,x_N)\in {\mathbb{R}}^{Nd}\mid x_i=x_{j}for\; somei\ne j\}.$$

$$W_2^2({\mu }_X,{\mu }_{X+V})=\frac{1}{N}\Vert V{\Vert }^2.$$

$$\begin{cases}\frac{1}{N}\dot{X}^{N}(t)=-\nabla E^{N}\bigl{(}X^{N}(t)\bigr{)},\\ X^{N}(0)=X^{N}_{0},\end{cases}$$

$$\underset{\sigma \in K}{\min }{W}_2^2(\mu ,\sigma ).$$

$$\frac{\partial E^{N}}{\partial x_{i}}(X)=\frac{1}{N}\nabla V(x_{i})+\frac{1}{N\varepsilon_{N}}\bigl{(}x_{i}-\beta_{i}(X)\bigr{)},$$

$$\begin{cases}\displaystyle\dot{x}_{i}^{N}(t)=-\nabla V\bigl{(}x_{i}^{N}(t)\bigr{)}+\frac{1}{\varepsilon_{N}}\bigl{[}\beta_{i}\bigl{(}X^{N}(t)\bigr{)}-x_{i}^{N}(t)\bigr{]},\\ X^{N}(0)=X^{N}_{0}.\end{cases}$$

$$X^{\prime }(t)\in -{\partial }^{+}{E}_{\epsilon }(X(t)),$$

$$|X_1(t)-X_2(t)|\ge e^{-Ct}|X_1(0)-X_2(0)|,$$

$$\frac{1}{{\epsilon }_N}{W}_2^2({\rho }_0,{\mu }_N(0))\le C,\underset{N\to \infty }{\lim }{\epsilon }_N=0.$$

$$\frac{1}{{\epsilon }_N^2}{\int }_0^T{W}_2^2\left({\sigma }_N,\frac{1}{N}\sum _{i=1}^N{\delta }_{{\beta }_i(X^N)}\right)\mathrm{d}t\le C,$$

$$M_{N}(t)=\sum_{i=1}^{N}\dot{x}_{i}^{N}(t)\frac{1}{N}\delta_{x_{i}^{N}(t)}=-\frac{1}{N}\sum_{i=1}^{N}\left[\nabla V\bigl{(}x_{i}^{N}(t)\bigr{)}+\frac{1}{\varepsilon_{N}}\bigl{(}x_{i}^{N}(t)-\beta_{i}^{N}(t)\bigr{)}\right]\delta_{x_{i}^{N}(t)}.$$

Full text

1 Introduction

In this paper we present an approximation scheme to solve evolution PDEs which have a gradient-flow structure in the space of probability measures $\mathcal{P}(\Omega)$ endowed with the Wasserstein distance $W_{2}$. Here, $\Omega\subset\mathbb{R}^{d}$ is a given compact domain where the evolution takes place, and the PDE will naturally be complemented by no-flux boundary conditions. The approximation that we present is Lagrangian in the sense that the evolving measure $\rho_{t}$ will be approximated by an empirical measure of the form $\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}(t)}$ and we will look for the evolution of the points $x_{i}$. We will use this approximation to provide an efficient numerical method, based on the most recent developments in semi-discrete optimal transport [6, 13, 12, 3]. Here, “semi-discrete” refers to the fact that the discretization of the diffusion effects in the evolution equation involves computation of the optimal transport plans between an empirical measure and diffuse measures.

Starting from the pioneering work of Otto and Jordan-Kinderlherer-Otto [18, 9] it is well-known that some linear and non-linear diffusion equations can be expressed in terms of a gradient flow in the space $W_{2}(\Omega)$. More precisely, the Fokker-Planck equation

[TABLE]

is the gradient flow of the energy $E(\rho):=\int\rho\log\rho+\int Vd\rho$, and the porous-medium equation

[TABLE]

is the gradient flow of the energy $E(\rho):=\int\frac{1}{m-1}\rho^{m}+\int Vd\rho$, Recently, also the limit case $m=\infty$ has been considered, in the framework of crowd motion [16]. In this case, the functional is $E(\rho):=\int Vd\rho$ if $\rho$ satisfies the constraint $\rho\leq 1$, and $E(\rho)=+\infty$ otherwise; the corresponding PDE is

[TABLE]

One can see the appearance of a pressure $p$ accounting for the constraint $\rho\leq 1$. We will come back later to the precise meaning and formulations of this last equation.

Approximation by empirical measures

Since any probability measure can be approximated by empirical measures, it is tempting to perform an approximation scheme just by considering the gradient flow of one of the above energy functionals $E$ on the set $\mathcal{P}_{N}(\Omega)$ of uniform measures on $N$ atoms, and then let $N\to\infty$. Unfortunately, the domain of the above functionals is reduced to absolutely continuous measures, and its intersection with $\mathcal{P}_{N}(\Omega)$ is empty. The main idea and novelty of this paper is to write $E(\rho)=F(\rho)+\int V\mathop{}\!\mathrm{d}\rho$ where $F$ is the entropy or the congestion constraint,

[TABLE]

and to replace $F$ by its Moreau-Yosida regularization

[TABLE]

The energies $F_{\varepsilon}$ are finite and well-defined for arbitrary probability measures $\mu$, and converge to $F$ as $\varepsilon\to 0$. More importantly, we will see that it is possible to compute very efficiently $F_{\varepsilon}(\mu)$ when $\mu\in\mathcal{P}_{N}(\Omega)$. The evolution of the discrete measures is then dealt with by keeping track of the positions of the particles $X=(x_{1},\dotsc,x_{N})\in\mathbb{R}^{Nd}$ in the support of the associated measure $\mu_{X}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}}$. Thanks to the correspondence between $X$ and $\mu_{X}$, we can think of $F_{\varepsilon}$ as an energy on the space of particle positions too, given by

[TABLE]

The discrete gradient flow then takes the form of a system of ODEs

[TABLE]

for a suitable choice of $\varepsilon_{N}\to 0$. The particles are only coupled by the forces $-\nabla_{x_{i}}F_{\varepsilon_{N}}(x_{1},\ldots,x_{N})$ due to diffusion or to the congestion constraint. Finally, we note that the initial condition can be selected by optimal quantization of the initial density $\rho_{0}\in\mathscr{P}(\Omega)$,

[TABLE]

granting, in many situations, an initial error of approximately $W_{2}^{2}(\rho_{0},\mu_{X_{0}^{N}})\lesssim(1/N)^{1/d}$.

Convergence

In the paper we will present the approximation scheme and prove, in these two cases, the convergence of the curves of empirical measures to the solution of the corresponding PDE, under the assumption that a certain bound on the approximate solutions themselves is satisfied. This assumption is unnatural, and it would be desirable to remove it, or replace it with an assumption on the approximation of the initial data. Yet, this seems to be a non-trivial problem, which is closely related to the general question of the convergence of gradient flows once the functionals $\Gamma$-converge. We refer to [19, 22] as classical papers on this question. In these papers, a semi-continuity property on the slopes of the functionals is required, which is in the same spirit of the bounds we need. These required bounds are stronger in the crowd motion case, as the equation is non-linear and stronger compactness is needed, while sligthly weaker in the Fokker-Planck case, which is indeed a linear equation. We show that such bounds can be obtained in dimension 1 (see §4). However, we insist that these bounds can be verified numerically in general, which makes the scheme we propose interesting for the approximation in arbitrary dimension.

Comparison to existing Lagrangian schemes

The optimal transport interpretation of advection-diffusion equations by Otto and Jordan-Kinderlherer-Otto [18, 9] has already led to many Lagrangian schemes:

• •

In dimension $d=1$ it is quite easy to construct such Lagrangian schemes. This is due to the fact that the Wasserstein space $(\mathscr{P}(\mathbb{R}),W_{2})$ can be isometrically embedded into $L^{2}([0,1])$ through the inverse cumulative distribution function. Discretizing probability densities in a lagrangian way then amounts to discretizing inverse cdfs, which can be done using finite elements of order one [2] or two [15], leading respectively to piecewise constant or piecewise linear densities.

• •

In dimension $d\geq 2$, one cannot isometrically embed the Wasserstein space in a $L^{p}$ space, making the choice of discretization less canonical. There exists discretization based on piecewise constant densities over Voronoi cells [5], or using Gaussian mixtures (“blobs”) [4]. Evans, Savin and Gangbo [8] have proposed a general way to rewrite some Wasserstein gradient flows as gradient flows in the space of diffeomorphisms (the idea is to write $\rho_{t}=s_{t\#}\rho_{0}$ where $s_{t}$ is a time-dependent diffeomorphism), which is used to construct a Lagrangian discretization in [10].

Our scheme could be considered as a variant of the scheme of [5] involving Voronoi cells but it is, to the best of our knowledge, the first one to use Laguerre cells (see §5.1), which are objects intrinsically related to the Wasserstein structure of the equation. The use of Laguerre cells, easy to handle via modern computational geometry tools, is now state-of-the-art when the transport cost is quadratic. However, we stress that we only use quadratic transport costs in the computation of the Moreau-Yosida regularization of the diffusion or congestion term, which means that we could a priori think of attacking with the same ideas other PDEs, which have a gradient structure for other distances. We refer for instance to [1] for PDEs induced by a cost functions of the form $c(x,y)=|x-y|^{q}$ or [17] for the so-called relativistic heat equation.

2 Lagrangian discretization of crowd motion

Formulation of the continuous problem

We fix a compact domain $\Omega\subset\mathbb{R}^{d}$ and a potential $V\in C^{1}(\mathbb{R}^{d})$ bounded from below, e.g. $V\geq 0$. The crowd is described by a probability measure $\rho$ in $\Omega$. Each agent tries to follow the gradient vector field $-\nabla V$ while ensuring that the probability density satisfies the density constraint $\rho\leq 1$. Therefore, we introduce the constraint set

[TABLE]

and we assume that $\lvert\Omega\rvert\geq 1$ so that $K$ is non-empty. There are a few possible ways to express this idea of constrained motion, which are at least formally equivalent. The first version is straightforward, but a bit problematic to formulate rigorously. The mass evolution is expressed by a continuity equation where the driving vector field is the projection of $-\nabla V$ onto the tangent cone of the set $K$:

[TABLE]

where

[TABLE]

Note that the choice of the boundary conditions on $v\cdot n$ is not relevant at this stage, since anyway the equation $\partial_{t}\rho+\operatorname{div}(\rho v)=0$ is already intended with no-flux boundary conditions $\rho v\cdot n=0$, i.e. $v\cdot n=0$ on $\{\rho>0\}$. Choosing to impose $v\cdot n\leq 0$ on the part of the boundary where the density $\rho$ is saturated allows for an easier expression of the normal cone below. Indeed, projecting to the tangent cone amounts to subtracting a vector from the normal cone, dual to the tangent cone, thereby leading to

[TABLE]

where

[TABLE]

This can be seen by observing that $T_{\rho}K$ contains all divergence-free vector fields and hence, by Helmoltz decomposition, any vector $w$ in the dual cone $N_{\rho}K$ is such that $\rho w$ is the gradient of a scalar function. Once we identify that this function should vanish on $\{\rho<1\}$, the constraints on the pressure $p$ are obtained by dualizing those defining $T_{\rho}K$. More explicitly,

[TABLE]

where $p$ is to be thought of as a pressure field enforcing the density constraint. Finally, solutions to (2.2) also formally coincide with the gradient flow

[TABLE]

with respect to the Wasserstein distance of the energy $E:\mathscr{P}(\Omega)\to\mathbb{R}\cup\{\infty\}$ given by

[TABLE]

The equation (2.5) is a differential inclusion, meaning that a priori the element of the subdifferential which is selected is not known. This translates into the fact that the pressure only belongs to a cone and is not an explicit function of the density, as it happens in other equations, for instance of porous medium type. However, in many differential inclusion problems there is indeed uniqueness, corresponding to the idea that only choice of the pressure will preserve the constraint: at a very formal level, we can say that even if we only imposed $\nabla\cdot v\geq 0$ on the saturated region — and this was done in order not to increase the density in the future after saturation — we should also have $\nabla\cdot v=0$, at least for $t>0$, otherwise the density would violate the constraint in the past. This means that $p$ should be uniquely determined as the solution of $-\Delta p=\Delta V$ on the saturated region $\{\rho=1\}$, with Dirichlet boundary condition on the part of the boundary of this region contained in the interior of $\Omega$ and suitable non-homogeneous Neumann condition $\nabla p\cdot n=-\nabla V\cdot n$ on the part contined in $\partial\Omega$.

Discretization

As explained in the introduction, our strategy for the numerical solution of the crowd motion is to employ a Lagrangian discretization of (2.5), meaning that we consider the time evolution of a probability measure which remains uniform over a set containing $N$ points:

[TABLE]

Since the intersection between the constraint set $K=\{\rho\in\mathscr{P}_{\text{ac}}(\Omega):\rho\leq 1\}$ and $\mathscr{P}_{N}(\mathbb{R}^{d})$ is empty, we are forced to replace the constraint $\rho\in K$ with a penalization, therefore considering the regularized energy given by

[TABLE]

Note that $\frac{1}{2\varepsilon}\mathsf{d}_{K}^{2}$ is the Moreau-Yosida regularization of the convex indicator function $F=i_{K}$ of the set $K$.

The evolution of the discrete measures is dealt with by keeping track of the positions of the particles $X=(x_{1},\dotsc,x_{N})\in\mathbb{R}^{Nd}$, to which corresponds the associated measure $\mu_{X}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}}$. Thanks to the correspondence between $X$ and $\mu_{X}$, we can think of $E_{\varepsilon}$ as an energy on the space of particle positions too, given by

[TABLE]

Assume for a moment that the point set $X\in\mathbb{R}^{Nd}$ does not belong to the set

[TABLE]

Then it is easy to see that for small perturbations $V=(v_{1},\ldots,v_{N})\in\mathbb{R}^{Nd}$, the optimal transport map between $\mu_{X}$ and $\mu_{X+V}$ simply maps $x_{i}$ to $x_{i}+v_{i}$, thus showing that the map $X\mapsto\mu_{X}$ is locally isometric, i.e. if $\lVert V\rVert$ is small enough,

[TABLE]

This suggests to replace the continuous Wasserstein gradient flow (2.5) by the discrete gradient flow (2.7) with respect to the Euclidean metric:

[TABLE]

where $E^{N}=E_{\varepsilon_{N}}$ for a suitable choice of $\varepsilon_{N}\to 0$. The initial condition can be selected by optimal quantization as in (1.3).

Given $\mu\in\mathscr{P}_{N}(\mathbb{R}^{d})$, let $\sigma\in K$ be the projection of $\mu$ onto $K$, i.e. a minimizer of

[TABLE]

Its existence follows by compactness, while its uniqueness and continuity with respect to $\mu$ are guaranteed by Proposition 5.2 in [7]. Let also $T:\Omega\to\mathbb{R}^{d}$ be the (unique) optimal transport map from $\sigma$ to $\mu$. The cell $L_{i}=T^{-1}(x_{i})$ represents the part of the mass of $\sigma$ which is attached to the particle $x_{i}$. Denoting by \beta_{i}(X)=\mathop{}\mkern-3.0mu\mathchoice{\hbox to0.0pt{\displaystyle\vbox{\hbox to6.14195pt{\hss$\textstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\textstyle\vbox{\hbox to6.14195pt{\hss$\scriptstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\scriptstyle\vbox{\hbox to4.78084pt{\hss$\scriptscriptstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\scriptscriptstyle\vbox{\hbox to4.78084pt{\hss$\scriptscriptstyle{-}$\hss}}\hss}}\mkern-3.0mu\int_{L_{i}}x\mathop{}\!\mathrm{d}\sigma(x)=\int_{L_{i}}x\mathop{}\!\mathrm{d}\sigma(x)/\int_{L_{i}}\mathop{}\!\mathrm{d}\sigma(x) the barycenter of the cell $L_{i}$, subsection 5.1 gives

[TABLE]

and therefore (2.7) becomes more explicitly

[TABLE]

Moreover, subsection 5.1 shows that $\frac{1}{2\varepsilon}\mathsf{d}^{2}_{K}$ is $\frac{1}{2N\varepsilon}$-concave, which proves that the vector field $-\nabla E_{\varepsilon}(X)$ is well-defined a.e. and provides several useful properties of the flow of this vector field. In particular, following [21] (slightly adapting the proof of Proposition 2.3), one can prove the existence, for every initial datum, of a solution of

[TABLE]

where $\partial^{+}$ is the superdifferential of semiconcave functions. Solutions of this ODE satisfy a reverse Gronwall inequality which provides

[TABLE]

opposite to what happens for the gradient flows of semiconvex functions: this is not useful for uniqueness purposes, but states that solutions cannot concentrate too much, and in particular we obtain that for a.e. initial datum the flow avoids for a.e. time the non-differentiability set. In particular, we have existence of solutions of $X^{\prime}(t)=-\nabla E_{\varepsilon}(X(t)),$ in the almost everywhere sense and for almost every initial datum. Therefore it is always possible to find $\mu_{N}(0)$ and $X_{N}$ that satisfy the hypothesis of the following theorem. However, we insist that the goal of the present paper is to show approximation results, and the existence proof for the “discrete” problem with a finite number of particles is not the core of our analysis, which explains why we do not provide more details about the existence for a.e. intial datum.

Theorem 2.1** (Convergence of the discrete scheme).**

For every $N\in\mathbb{N}$, let $\varepsilon_{N}\in(0,\infty)$ and $\mu_{N}(0)\in\mathscr{P}_{N}(\mathbb{R}^{d})$ be such that

[TABLE]

Let $X_{N}\in C^{1}([0,T],\mathbb{R}^{Nd})$ be a solution of (2.7) or, equivalently, (2.8) and let $\mu_{N}:[0,T]\to\mathscr{P}_{N}(\mathbb{R}^{d})$ be the corresponding curve of measures. Assume that

[TABLE]

for some constant independent of $N$ and that $\rho_{0}\in K$. Then, as $N\to\infty$, and up to subsequences, $\mu_{N}\to\rho$ in C^{0}\bigl{(}[0,T];W_{2}(\mathbb{R}^{d})\bigr{)}, where $\rho$ is a weak solution to (2.4).

Proof.

For simplicity of notation, let us write $\beta_{i}^{N}(t)$ in place of \beta_{i}\bigl{(}X^{N}(t))\bigr{)}. Define the time dependent vector valued measures

[TABLE]

Define also the space-time measures $\mu_{N}\in\mathscr{M}_{+}([0,T]\times\Omega)$ and $M_{N}\in\mathscr{M}([0,T]\times\Omega;\mathbb{R}^{d})$ given by

[TABLE]

By construction they satisfy the equation

[TABLE]

in a weak sense because, for any $\zeta\in C^{1}_{c}(\mathbb{R}^{d})$, one has

[TABLE]

This also means that $(\mu_{N},M_{N})$ solve the continuity equation

[TABLE]

in distributional sense.

The first step is showing that the measures $\mu_{N}$ admit a limit in some sense. This is a consequence of the energy estimate

[TABLE]

because then the Benamou-Brenier formula for the $W_{2}$ distance

[TABLE]

shows that the functions [0,T]\to\bigl{(}\mathscr{P}(\Omega),W_{2}\bigr{)}:t\mapsto\mu_{N}(t) are equi-continuous, since they are all $1/2$-Hölder with the same constant. Ascoli-Arzelà then ensures that $\mu_{N}\to\rho$ in $C([0,T],W_{2}(\Omega))$, up to a subsequence. In particular, $\mu_{N}\rightharpoonup\rho$ in $\mathscr{M}_{+}([0,T]\times\Omega)$.

Next, we show that also the family of measures $M_{N}$ admits a limit. Indeed,

[TABLE]

and by compactness in the space of measures, they admit a weak limit $M_{N}\rightharpoonup M$ in $\mathscr{M}([0,T]\times\Omega;\mathbb{R}^{d})$.

In particular, the weak convergence of $\mu_{N}$ and $M_{N}$ is sufficient to pass to the limit (2.11) and infer that

[TABLE]

To show that $M\ll\rho$, we will use a few properties of the Benamou-Brenier functional $\mathcal{B}_{2}:\mathcal{M}([0,T]\times\Omega)\times\mathcal{M}([0,T]\times\Omega,\mathbb{R}^{d})\to\mathbb{R}\cup\{+\infty\}$. These properties are found in Proposition 5.18 of [20]: (i) If $M\ll\mu$, $\mathcal{B}_{2}(\mu,M)=\int_{[0,T]\times\Omega}\left\lvert\frac{\mathop{}\!\mathrm{d}M}{\mathop{}\!\mathrm{d}\mu}\right\rvert^{2}\mathop{}\!\mathrm{d}\mu$. (ii) The functional $\mathcal{B}_{2}$ is lower semi-continuous wrt narrow convergence. (iii) $\mathcal{B}_{2}(\mu,M)<+\infty$ only if $\mu\geq 0$ and $M\ll\mu$. In our case, using the fact that $M_{N}\ll\mu_{N}$ and previous computations, we get that

[TABLE]

is uniformly bounded. Then, by lower semi-continuity, $\mathcal{B}_{2}(\rho,M)$ is finite. This implies, by the third property of $\mathcal{B}_{2}$ that $M\ll\rho$.

Let now $\sigma_{N}$ be the projection of $\mu_{N}$ on $K$ and let $T_{N}:\Omega\to\mathbb{R}^{d}$ be the optimal transport map $(T_{N})_{\#}\sigma_{N}=\mu_{N}$. Notice that

[TABLE]

therefore $\sigma_{N}(t)$ converges to the same limit $\rho(t)$ as $\mu_{N}(t)$. We used that $\rho_{0}\in K$ so that $\mathop{}\!\mathrm{d}_{K}(\rho_{0})=0$. In particular this means that $\rho(t)\in K$ for all $t$.

For $\xi\in C^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$ and omitting time dependence for brevity, we have

[TABLE]

By Brenier’s Theorem and the particular structure of the optimal partial transport problem, we have that

[TABLE]

where $\varphi_{N}:\Omega\to\mathbb{R}$ is a semi-concave function satisfying $\varphi_{N}\leq 0$ and $(1-\sigma_{N})\varphi_{N}=0$; see [16, Lemma 3.1] or Proposition 5.1 of this paper. If we introduce the pressure field

[TABLE]

we have that

[TABLE]

We must show that $p_{N}\rightharpoonup p$ in L^{2}\bigl{(}[0,T];H^{1}(\Omega)\bigr{)} to some admissible pressure field. This follows from the equi-boundedness

[TABLE]

together with a Poincaré inequality based on the fact that $|\{p_{N}=0\}|\geq|\Omega|-1$. This allows to transform $L^{2}$ bounds on the gradients into full $H^{1}$ bounds. We now prove the aforementioned equi-boundedness. For every $t\in[0,T]$, which we omit for brevity of notation, we have that

[TABLE]

where in the last step we use that $\int_{L_{i}}(y-\beta_{i})\mathop{}\!\mathrm{d}\sigma_{N}(y)=0$.

The first term can be treated with the bound given by Assumption (2.9). For the second term, notice that

[TABLE]

by (2.12). In conclusion,

[TABLE]

The next step is to show that $p(1-\rho)=0$. The difficulty is that both $\sigma_{N}$ and $p_{N}$ are converging weakly, which is not sufficient in order to pass to the limit the nonlinear relation $p_{N}(1-\sigma_{N})=0$.

For $0\leq t_{0}<t_{1}\leq T$, let us introduce the average pressure

[TABLE]

Define also the measures $\lambda_{N}\in\mathscr{M}_{+}([0,T])$ given by

[TABLE]

Their total masses

[TABLE]

are uniformly bounded, as previously shown; therefore, up to a subsequence, they converge weakly to a measure $\lambda\in\mathscr{M}_{+}([0,T])$. Note that the second equality in the previous formula uses that $p_{N}(1-\sigma_{N})=0$.

For every $N$, we split in two pieces the following identity:

[TABLE]

In order to deal with the first integral, we observe that we have strong convergence $p_{N}^{t_{0},t_{1}}\xrightarrow{L^{2}(\Omega)}p^{t_{0},t_{1}}$ since $p_{N}^{t_{0},t_{1}}$ is bounded in $H^{1}$ (as it is obtained as an average, for fixed $t_{0},t_{1}$ of a time-dependent function on which we have $L^{2}_{t}H^{1}_{x}$ bounds), and that we have $\sigma_{N}(t_{0})\rightharpoonup\rho(t_{0})$ as $N\to\infty$ (this convergence is a weak convergence of measures, but it is also weak in $L^{2}$ because of the $L^{\infty}$ bounds on the densities). Hence, the first integral converges to

[TABLE]

At any Lebesgue point $t_{0}$ of the map $[0,T]\to L^{2}(\Omega):t\mapsto p(t)$ we have $p^{t_{0},t_{1}}\xrightarrow[t_{1}\to t_{0}]{}p(t_{0})$, hence

[TABLE]

Employing Section 2, the second integral can be estimated as

[TABLE]

where $\omega$ is a continuity modulus for the curve $t\mapsto\sigma_{N}(t)$ in the Wasserstein space $W_{2}$. Note that the continuity of the curve $t\mapsto\sigma_{N}(t)$ comes from the continuity of $t\mapsto\mu_{N}(t)$ (a consequence of (2.13)) and the continuity of the projection operator $\mu\mapsto\arg\min_{\sigma}W_{2}^{2}(\mu,\sigma)$; see, for instance, [7]).

When $N\to\infty$, for almost every $t_{0}$ and $t_{1}$ we have

[TABLE]

which tends to a finite constant when $t_{1}\to t_{0}$ for a.e. $t_{0}$. This latter fact comes from differentiation of measures: this quantity represents the density of the measure $\lambda$ w.r.t. the Lebesgue measure on the real line, and the limit exists and is equal to the density of the absolutely continuous part of $\lambda$ a.e. If we also consider the factor $\omega(t_{1}-t_{0})$, we see that the second integral goes to [math] for almost every $t_{0}$ when taking the limits $N\to\infty$ and $t_{1}\to t_{0}$, in this order.

Summing up, we have shown that

[TABLE]

for almost every $t_{0}\in[0,T]$, which proves $p(1-\rho)=0$, by the positivity of $p$.

We can finally show that $M=(-\nabla V-\nabla p)\rho$. Fix $\xi\in C^{1}([0,T]\times\Omega;\mathbb{R}^{d})$. By (2.14), we know that

[TABLE]

The first term passes to the limit because

[TABLE]

For the second term,

[TABLE]

therefore

[TABLE]

The following lemma is borrowed from [16, Lemma 3.5] (but was first presented in other papers, such as [14]).

Lemma 2.2**.**

Let $\mu_{0},\mu_{1}\in\mathscr{P}(\Omega)$ be probability measures with densities bounded by $1$. Then for all $f\in H^{1}(\Omega)$ we have that

[TABLE]

Remark 2.3*.*

The convergence result can be generalized to handle PDEs involving other terms such as self-interaction involving a $\mathcal{C}^{1,1}$ kernel $W$:

[TABLE]

which can be regarded as the gradient flow, in $(\mathscr{P}(\Omega),W_{2})$ of the energy

[TABLE]

provided the Kernel $W$ is even. In this case, the discretized system becomes

[TABLE]

The only difference with respect to Theorem 2.1 will the presence of a modified velocity field $\nabla V_{N}$ instead of $\nabla V$, with $V_{N}=V+W*\mu_{N}$. In the proof of Theorem 2.1 we used the weak convergence of $\mu_{N}$ to handle the term $\int\nabla V\cdot\xi d\mu_{N}$; we will now also need the uniform convergence $\nabla V_{N}\to\nabla(V+W*\rho)$ (which is a consequence of the regularity assumption on $W$) to handle the same term. Also note that in the definition of the flow one can omit the term $\nabla W(x_{i}-x_{j})$ for $i=j$, as it is usually done, since anyway $\nabla W(0)=0$.

Theorem 2.4** (Convergence of the discrete scheme in 1D).**

Let $\Omega\subset\mathbb{R}$ be an interval. For every $N\in\mathbb{N}$, let $\varepsilon_{N}=1/N$ and $\mu_{N}(0)\in\mathscr{P}_{N}(\mathbb{R})$ be such that

[TABLE]

Let $X_{N}\in C^{1}([0,T],\mathbb{R}^{N})$ be a solution of (2.7) or, equivalently, (2.8) and let $\mu_{N}:[0,T]\to\mathscr{P}_{N}(\mathbb{R})$ be the corresponding curve of measures. Then, as $N\to\infty$, and up to subsequences, $\mu_{N}\to\rho$ in C^{0}\bigl{(}[0,T];W_{2}(\mathbb{R})\bigr{)}, where $\rho$ is a weak solution to (2.4).

Proof.

This is an immediate consequence of Theorem 2.1 and Section 4, which allows to verify the assumption

[TABLE]

3 Lagrangian discretization of linear diffusion

The previously presented Lagrangian scheme can be adapted to solve also the advection-diffusion equation

[TABLE]

on a bounded domain $\Omega$, with no-flux boundary conditions. This equation arises as the gradient flow with respect to $W_{2}$ of the energy

[TABLE]

We adopt the same Lagrangian discretization as before. For the atomic measures $\mu_{N}$, the entropy (the second term in the energy) is identically $+\infty$, therefore we need to substitute it with a similar functional, in the same manner that we replaced the hard constraint $\rho\in K$ with a penalization. To this end, we consider its Moreau-Yosida regularization

[TABLE]

and the new energy becomes $E_{\varepsilon}(\rho)=\int_{\Omega}V\mathop{}\!\mathrm{d}\rho+H_{\varepsilon}(\rho).$ Letting $F_{\varepsilon}(x_{1},\ldots,x_{N})=H_{\varepsilon}(\mu_{X})$, the discrete measure $\mu_{N}(t)$ represented by the particles $X_{N}(t)$ can then evolve according to the system of ODE as before, namely

[TABLE]

Theorem 3.1** (Convergence of the discrete scheme).**

For every $N\in\mathbb{N}$, let $\varepsilon_{N}\in(0,\infty)$ and $\mu_{N}(0)\in\mathscr{P}_{N}(\mathbb{R}^{d})$ such that

[TABLE]

Let $X_{N}\in C^{1}([0,T],\mathbb{R}^{Nd})$ be a solution of (3.3) and let $\mu_{N}:[0,T]\to\mathscr{P}_{N}(\mathbb{R}^{d})$ be the corresponding curve of measures. Assume that

[TABLE]

where $\sigma_{N}$ is the minimizer in the definition of $H_{\varepsilon}(\mu_{N})$ (see (3.2)). Then, as $N\to\infty$, and up to subsequences, $\mu_{N}\to\rho$ in C^{0}\bigl{(}[0,T];W_{2}(\mathbb{R}^{d})\bigr{)}, where $\rho$ is a weak solution to (3.1).

Proof.

Define as before the vector measures

[TABLE]

Together with $\mu_{N}$, they solve the continuity equation

[TABLE]

Moreover,

[TABLE]

where we used that the domain is bounded to see that the entropy is bounded from below by some constant, implying E^{N}\bigl{(}X^{N}(0)\bigr{)}\geq C. This shows that the functions [0,T]\to\bigl{(}\mathscr{P}(\Omega),W_{2}\bigr{)}:t\mapsto\mu_{N}(t) are equi-continuous, because they are $1/2$-Hölder with the same constant. Ascoli-Arzelà then ensures that $\mu_{N}\rightharpoonup\rho$, up to a subsequence.

The rest of the proof is similar to the previous one with the following modifications.

The measure $\sigma_{N}$ minimizing

[TABLE]

satisfies

[TABLE]

where $\varphi_{N}$ is the optimal potential from $\sigma_{N}$ to $\mu_{N}$ and $c_{N}$ is a normalization constant. This optimality condition can be recovered from the first variation of the objective functional. This implies that $\mathscr{L}^{d}$-almost everywhere on $\Omega$ the density $\sigma_{N}$ is strictly positive and

[TABLE]

Passing to the limit $M_{N}$ in order to get $M=-\nabla V\rho-\nabla\rho$ is now easier because

[TABLE]

as before. Setting $p_{N}=\log\sigma_{N}$ (which is the term which plays a similar role to that of the pressure in the previous section), we have for any $\xi\in\mathcal{C}^{0}_{c}(\Omega)$,

[TABLE]

The first step in the above equation is justified because

[TABLE]

The last term tends to [math] by writing, as in (2.15),

[TABLE]

The first term tends to [math] by assumption, and the second term is $O(\varepsilon_{N})$ because of (2.16). ∎

Theorem 3.2** (Convergence of the discrete scheme in 1D).**

Let $\Omega\subset\mathbb{R}$ be a bounded interval. For every $N\in\mathbb{N}$, take a number $\varepsilon_{N}>0$ and $\mu_{N}(0)\in\mathscr{P}_{N}(\mathbb{R})$ such that

[TABLE]

Let $X_{N}\in C^{1}([0,T],\mathbb{R}^{N})$ be a solution of (3.3) and let $\mu_{N}:[0,T]\to\mathscr{P}_{N}(\mathbb{R})$ be the corresponding curve of measures. Then, as $N\to\infty$, and up to subsequences, $\mu_{N}\to\rho$ in C^{0}\bigl{(}[0,T];W_{2}(\mathbb{R})\bigr{)}, where $\rho$ is a weak solution to (3.1).

Proof.

This is an immediate consequence of Theorem 3.1 and Section 4, which provide

[TABLE]

4 Bounds in 1D

In this section we prove that, for both the crowd motion and the linear diffusion discretizations, in one dimension there are bounds on the quantities which are relevant for the application of Theorems 2.1 and 3.1. The results come from a static analysis, in the sense that the evolution equations do not play any role in the estimates.

We begin with the easier case of the crowd motion.

Proposition 4.1**.**

Let $\Omega\subset\mathbb{R}$ be an interval. Let $x_{1},\dots,x_{N}\in\mathbb{R}$, let $\mu_{N}$ be the corresponding atomic measure, $\sigma_{N}$ its $W_{2}$-projection on $\{\rho\in\mathscr{P}(\Omega):\rho\leq 1\}$, $T$ the optimal transport map between $\sigma_{N}$ and $\mu_{N}$. Define the Laguerre cells and their barycenters as

[TABLE]

[TABLE]

Then, choosing $\varepsilon=1/N$,

[TABLE]

Proof.

Each Laguerre cell $L_{i}$ is an interval of length $1/N$ and its barycenter is the midpoint. Moreover, $\sigma_{N}$ has constant density $1$ on every cell $L_{i}$, therefore

[TABLE]

which gives the claim. ∎

We now pass to the case which is relevant for linear diffusion.

Proposition 4.2**.**

Let $\Omega\subset\mathbb{R}$ be a bounded interval. Let $x_{1},\dots,x_{N}\in\mathbb{R}$, let $\mu_{N}$ be the corresponding atomic measure, and define, for $\varepsilon>0$:

[TABLE]

Let $\beta_{1},\dots,\beta_{N}$ be the barycenters of the Laguerre cells $L_{1},\dots,L_{N}$ of $\sigma_{N}$. Then we have

[TABLE]

where $C(\Omega)$ only depends on the length $|\Omega|$ of $\Omega$.

The proof of this proposition relies on the next lemma, which is a particular case of the main theorem of [11].

Lemma 4.3**.**

Let $\sigma=e^{-(y-x)^{2}/(2\varepsilon)}\mathop{}\!\mathrm{d}y$ be a Gaussian measure, let $L$ be an interval such that $\sigma(L)=\frac{1}{N}$. Then, the variance of N\sigma\mathop{\raisebox{-0.5468pt}{\reflectbox{\rotatebox[origin={br}]{-90.0}{\lnot}}}}L is upper bounded by the variance of $\sigma$:

[TABLE]

Proof of Section 4.

Let $\ell_{i}=\lvert L_{i}\rvert$ denote the length of the $i$-th Laguerre cell. We fix a parameter $\bar{\ell}\in(0,1)$ to be specified later and divide the cells in two groups:

• •

short cells: $S=\{i:\ell_{i}<\bar{\ell}\}$;

• •

long cells: $L=\{i:\ell_{i}\geq\bar{\ell}\}$.

Notice that $\lvert S\rvert\leq N$ and $\lvert L\rvert\leq\lvert\Omega\rvert/\bar{\ell}$. By subsection 5.1, we know that the restriction of $\sigma_{N}$ to the Laguerre cell $L_{i}$ is proportional to a Gaussian of the form $\exp(-\frac{1}{2\varepsilon}(x-x_{i})^{2})$. Therefore, using Section 4 to estimate the contribution of the long cells, we get

[TABLE]

where in the last step we chose $\bar{\ell}=\sqrt{\varepsilon}$. ∎

5 Numerical scheme

5.1 Computation of the Moreau-Yosida regularization

Let $F:\mathscr{P}_{\text{ac}}(\Omega)\to\mathbb{R}\cup\{+\infty\}$, which we assume to be lower-semicontinuous with respect to the Wasserstein metric $W_{2}$. We consider the Moreau-Yosida regularization

[TABLE]

and we assume that for every $X\in\mathbb{R}^{Nd}$, the minimization problem defining $F_{\varepsilon}(X)$ admits a unique solution. This assumption is satisfied in the two relevant cases for this paper, since the projection onto measures with bounded densities is always unique, see [7], and the minimizer in the entropy case is unique because of strict convexity. We let

[TABLE]

Our first proposition gives an explicit formulation for the gradient of $F_{\varepsilon}$ given a solution $\sigma$ of the minimization problem defininig $F_{\varepsilon}$. We recall that a function $F$ on $\mathbb{R}^{k}$ is $\lambda$-semi-concave if and only if $F-\lambda\lVert\cdot\rVert^{2}$ is concave.

Proposition 5.1**.**

$F_{\varepsilon}$* is $\frac{1}{2N\varepsilon}$-semi-concave on $\mathbb{R}^{Nd}$ and continuously differentiable on $\mathbb{R}^{Nd}\setminus\mathcal{D}_{N}$. Given $X=(x_{1},\ldots,x_{N})\in\mathbb{R}^{Nd}\setminus\mathcal{D}_{N}$, we let $\sigma$ the unique minimizer in (5.1), $T$ the unique optimal transport map between $\sigma$ and $\mu_{X}$, and $L_{i}=T^{-1}(x_{i})$. Then,*

[TABLE]

Proof.

First, let us underline that we work under the assumption that the optimal $\sigma$ in the minimization problem defining $F_{\varepsilon}(X)$ is unique for every $X$. This uniqueness implies continuity of the map $X\mapsto\sigma$. Let $X\in\mathbb{R}^{Nd}$, $\sigma$ the unique minimizer in (5.1) and $T$ the unique optimal transport map between $\sigma$ and $\mu_{X}$. Given $Y\in\mathbb{R}^{Nd}$,

[TABLE]

By construction, one can decompose $\sigma=\sum_{1\leq i\leq N}\sigma_{i}$ where $\sigma_{i}\geq 0$, $T(\sigma_{i})=x_{i}$, and such that $\sigma_{i}(\Omega)=1/N$. Considering the transport which maps $\sigma_{i}$-a.e. point to $y_{i}$, one gets

[TABLE]

where we have set $\beta_{i}(X)=N\int x\mathop{}\!\mathrm{d}\sigma_{i}$ and $\beta(X)=(\beta_{1}(X),\ldots,\beta_{N}(X))$. This inequality can be rewritten as

[TABLE]

where $G_{\varepsilon}(X)=F_{\varepsilon}(X)-\frac{1}{2N\varepsilon}\lVert X\rVert^{2}$. This shows that $-\frac{1}{N\varepsilon}\beta(X)$ belongs to the superdifferential of $G_{\varepsilon}$ at $X$, so that the superdifferential of $G_{\varepsilon}$ is never empty, also showing that $G_{\varepsilon}$ is concave. If $X\not\in\mathcal{D}_{N}$, then \sigma_{i}=\sigma\mathop{\raisebox{-0.5468pt}{\reflectbox{\rotatebox[origin={br}]{-90.0}{\lnot}}}}T^{-1}(x_{i}) and the point $\beta_{i}(X)$ is uniquely defined (we use here the hypothesis on the uniqueness of the minimal $\sigma$ in (5.1)). Using the stability of optimal transport maps, we get that $X\in\mathbb{R}^{Nd}\setminus D\mapsto\beta(X)$ is continuous on $\mathbb{R}^{Nd}\setminus D_{N}$, which shows that $G_{\varepsilon}\in\mathcal{C}^{1}(\mathbb{R}^{Nd}\setminus D_{N})$ and that

[TABLE]

The properties of $F_{\varepsilon}$ can be deduced from those of $G_{\varepsilon}$. ∎

The next two propositions explain how to compute the optimal $\sigma$ in the definition of the Moreau-Yosida regularization in the crowd motion and linear diffusion. Using Kantorovich duality, this problem can be reformulated as the computation of a Kantorovich potential satisfying a finite-dimensional non-linear system, (5.2) or (5.3).

Given $x_{1},\ldots,x_{N}\in\mathbb{R}^{d}$ and $\psi\in\mathbb{R}^{N}$ we define the Laguerre cell of the point $x_{i}$ with respect to $\Omega$ as

[TABLE]

Note that the next proposition is a special case of the characterization of Wasserstein projections on $K$, proven in [16, Lemma 3.1]. It is illustrated in Figure 1.

Proposition 5.2**.**

Consider $F:\mathscr{P}(\Omega)\to\mathbb{R}$ defined by $F(\mu)=0$ if $\mu\in K$ and $F(\mu)=+\infty$ otherwise, where $K$ is defined in (2.1). Then for all $X\in\mathbb{R}^{Nd}\setminus\mathcal{D}_{N}$, there exists $\psi\in\mathbb{R}^{N}_{-}$ such that

[TABLE]

Given such a $\psi$, define $\varphi=\min(\min_{i}\lVert\cdot-x_{i}\rVert^{2}-\psi_{i},0)$ and $\sigma=\mathbf{1}_{\{\varphi<0\}\cap\Omega}$.

• (a)

$\sigma\in\mathscr{P}_{\text{ac}}(\Omega)$* is the Wasserstein projection of $\mu_{X}$ on $K$,*

• (b)

$\varphi\leq 0,$* $\varphi(1-\sigma)=0$*

• (c)

$(\varphi,\psi)$* is an admissible pair of Kantorovich potential in the transport from $\sigma$ to *$\mu_{X}$

Proof.

By Kantorovich duality, one can write for any $\sigma\in\mathscr{P}_{\text{ac}}(\Omega)$,

[TABLE]

thus giving

[TABLE]

Switching the minimum and the maximum, we get the following dual problem

[TABLE]

where, we set

[TABLE]

and

[TABLE]

With similar arguments as in [12, Theorem 1.1], one can prove that $D$ is concave, $\mathcal{C}^{1}$, and that its partial derivatives are

[TABLE]

It is easy to see that the maximum is attained in the dual problem, thus proving the existence of $\psi\leq 0$ satisfying (5.2). Define $\sigma$ and $\varphi$ as in the statement. Then, the property $\varphi(1-\sigma)$ is obvious. In addition,

[TABLE]

so that the pair $(\varphi,\psi)$ is admissible in the dual Kantorovich problem. It is also optimal in the optimal transport problem between $\sigma$ and $\mu_{X}$ by construction, since $\varphi$ coincides with the $c$-transform of $\psi$ on the support of $\sigma$. This shows that

[TABLE]

Since the converse inequality holds by weak duality, we get strong duality, and in particular $\sigma$ is the solution to the primal problem. ∎

The following proposition, dealing with the linear diffusion case, is obtained in a very similar manner (one can for instance use Proposition 8.6 in [20] to get the optimality condition for the dual problem). We also refer the reader to Theorems 3.1–3.2 in [3], where similar results are shown for more general functionals.

Proposition 5.3**.**

Let $H:\mathscr{P}(\Omega)\to\mathbb{R}$ be Boltzmann’s functional, and $X\in\mathbb{R}^{Nd}\setminus D_{N}$. Then, there exists $\psi\in\mathbb{R}^{N}$ such that

[TABLE]

Given such a $\psi$, define $\varphi=\min_{i}\lVert\cdot-x_{i}\rVert^{2}-\psi_{i}$ and $\sigma=e^{-\frac{\varphi}{2\varepsilon}}\mathbf{1}_{\Omega}$. Then,

• (a)

$\sigma\in\mathscr{P}_{\text{ac}}(\Omega)$* is the unique minimizer of $\min_{\mathscr{P}_{\text{ac}}(\Omega)}\frac{1}{2\varepsilon}W_{2}^{2}(\cdot,\mu_{X})+H(\cdot)$.*

• (b)

$\frac{1}{2\varepsilon}\varphi+\log(\sigma)=0$**

• (c)

$(\varphi,\psi)$* is a pair of optimal Kantorovich potentials in the transport from $\sigma$ to $\mu_{X}$.*

Remark 5.4*.*

In practice, equations (5.2) and (5.3) are solved using the same damped Newton algorithm as in [12]. The cells $L_{i}(\psi)$ are computed using computational geometry techniques ensuring a near-linear computational time in $2$D. The integrals are computed exactly in the crowd motion case, and using quadratures ensuring a negligible numerical error in the linear diffusion case.

Remark 5.5*.*

The computation of the Moreau-Yosida regularization of the congestion constraint and the entropy is implemented in the open-source library sd-ot, which is available at https://github.com/sd-ot. The numerical schemes for crowd motion and linear diffusion are implemented as examples in the Python package.

5.2 Numerical experiments (crowd motion)

In this paragraph, we consider $\Omega\subseteq\mathbb{R}^{2}$ a compact domain, $V:\Omega\to\mathbb{R}$ a potential, and we define as usual the congestion term $F:\mathscr{P}(\Omega)\to\mathbb{R}$ by $F(\mu)=0$ if $\mu$ has density $\leq 1$ and $+\infty$ if not. We consider the discretization of the crowd motion model explained above: an initial point set $X^{0}=(x_{1}^{0},\ldots,x_{N}^{0})$ is evolved through the ODE system

[TABLE]

which we discretize using a simple explicit Euler scheme:

[TABLE]

Propositions 5.1–5.1 can be used to compute the gradient of the regularized congestion term $F_{\varepsilon}$. Figure 1 illustrates this computation by showing a point set $X=(x_{1},\ldots,x_{N})$, the projected measure $\sigma\in\mathscr{P}_{\text{ac}}(\Omega)$, $\sigma\leq 1$ and the gradient $(\nabla_{x_{i}}F_{\varepsilon}(X))_{1\leq i\leq N})$.

Radial case

As a first test case, we consider a simple problem with radial symmetry, introduced in [16, Section 5], and whose solution is explicit. The domain is the set $\Omega=\{x\in\mathbb{R}^{2}\mid x_{2}\geq\lvert x_{1}\rvert,\lVert x\rVert\leq R\}$, and the potential is given by $V(x)=\lVert x\rVert$. In our experiment, we assume that $R=2$ and $\alpha=\frac{1}{\pi}$, so that $\rho^{0}=\alpha\mathbf{1}_{\Omega}$ is a probability measure. As shown in [16], the evolution of the crowd is then given by

[TABLE]

where $b$ is a solution of

[TABLE]

Given $h>0$, we denote $N_{h}=\mathrm{Card}(\Omega\cap h\mathbb{Z}^{2})$ and we let $x_{1}^{0},\ldots,x_{N_{h}}^{0}$ be the an arbitrary numbering of the points in the intersection $\Omega\cap h\mathbb{Z}^{2}$. In all experiments, we set $\tau=\frac{h}{2},$ $\varepsilon=h$ and $T=1$. Figure 2 displays the evolution of the Laguerre cells at six time steps. To get error estimates, we measure the Wasserstein distance between:

• •

$\bar{\rho}_{t}=(x\mapsto\lVert x\rVert)_{\#}\rho_{t}\in\mathscr{P}(\mathbb{R})$, which is the distribution of distances from the origin, computed from the exact solution $\rho_{t}$;

• •

$\bar{\mu}^{k}=\frac{1}{N_{h}}\sum_{1\leq i\leq N_{h}}\delta_{\lVert x_{i}^{k}\rVert}\in\mathscr{P}(\mathbb{R})$ the distribution of distances from the origin, computed on the discrete solution $\mu^{k}=\frac{1}{N_{h}}\sum_{1\leq i\leq N_{h}}\delta_{x_{i}}$.

The relation between $h$ and $\mathrm{err}_{h}=\max_{0\leq k\leq\frac{T}{\tau}}W_{2}(\bar{\rho}_{k\tau},\bar{\mu}^{k}),$ as reported in Table 1, suggests a near-linear convergence rate.

Bimodal case

In this case, is obtained as the union of three squares $\Omega=\Omega_{\ell}\cup\Omega_{r}\cup\Omega_{c}$: two “rooms” $\Omega_{\ell}$ and $\Omega_{r}$ joined by a corridor $\Omega_{c}$, where

[TABLE]

The crowd is initially located in the left room $\Omega_{\ell}$ and the potential $V$ is constructed as the distance function to the two corners $\{(\frac{7}{3}\alpha,\alpha),(\frac{7}{3}\alpha,0)\}$ of the right room $\Omega_{r}$. More precisely, $V$ is obtained as the solution to the following Eikonal equation, which is computed using a fast marching method:

[TABLE]

Given $h>0$, we denote $N_{h}=\mathrm{Card}(\Omega_{\ell}\cap h\mathbb{Z}^{2})$ and we let $x_{1}^{0},\ldots,x_{N_{h}}^{0}$ be the list of points in $\Omega_{\ell}\cap h\mathbb{Z}^{2}$, so that the crowd is initially located on the left square $\Omega_{\ell}$. Here, we set the final time to $T=3$, and as before, we have $\varepsilon=h$, $\tau=\frac{h}{2}$.

Visualization

In Figures 3 and 4, we visualize the distribution of the crowd at the timesteps $t_{i}=\frac{i/5}{T}$ for $0\leq i\leq 5$, for two space discretizations $h=\frac{\alpha}{30}$ and $h=\frac{\alpha}{80}$. In Figure 5.4, we highlight some features of the Lagrangian trajectory. First, given a particle $x_{i}^{0}$, we can compute the minimum time required for the particle to enter the right room:

[TABLE]

The exit time $\tau_{i}$ is displayed as a function of the particle coordinate $x_{i}^{0}$ at time $t=0$ in the first row of Figure 5.4. This figure shows (as one could expect) that the exit time is not proportional to the distance to the “door” $\Omega_{\ell}\cap\Omega_{c}$, as people in front tend to escape faster than those on side of the door. Finally, the next two rows of show the trajectory of the particles. The trajectories seem to be regular in time, but also seem to depend continuously on the initial condition (except near the non-differentiability locus of the potential $V$).

Remark 5.6* (On the assumption (2.9)).*

In the 2D cases treated here, we cannot guarantee that our numerical solutions converge to a solution of the crowd motion equation since convergence requires the estimate (2.9). However, it is quite easy to see that this estimate holds if one is able to prove that the diameter of the Laguerre cells is bounded uniformly by $C(1/N)^{1/d}$ – this is what is established in 1D in Section 4. Figure 3 (and in fact all our simulations) suggest that such an estimate is satisfied in practice.

Remark 5.7* (Lagrangian interpretation).*

The Eulerian crowd-motion equation (2.4) can be turned into a Lagrangian equation by introducing the map $s_{t}\in\mathrm{L}^{2}(\rho_{0},\Omega)$, which describes the displacement of the crowd from its position at time $t=0$ (more precisely, $s_{t}(x)$ is the position at time $t$ which was at $x$ at time [math]). Formally, $s$ should satisfy the following system

[TABLE]

and we expect that the numerical solution, shown in Figure 5.4, provides a piecewise-constant (in space) approximation to $s$. We note however that the system (5.5) has not been studied; showing existence of solutions to this system would require to better understand the regularity of the pressure $p$ appearing in (2.4).

5.3 Numerical experiments (diffusion)

In this paragraph, we consider $\Omega\subseteq\mathbb{R}^{2}$ a compact domain, and we let $F$ be Boltzmann’s functional. We consider the discretization of the heat equationl explained above: an initial point set $X^{0}=(x_{1}^{0},\ldots,x_{N}^{0})$ is evolved through the ODE system (3.3) (with $V=0$), which we discretize again using a simple explicit Euler scheme:

[TABLE]

In the numerical example presented in Figures 6,7 and 8, the initial density is uniform over a disk $D$, and approximated by the uniform measure over $h\mathbb{Z}^{2}\cap D$. Despite the lack of convergence result in 2D, one can observe the consistency between the simulations with $h=\frac{1}{30}$ and $h=\frac{1}{80}$. Some of the cells near the boundary of the disk $D$ are very elongated; however this does not a priori prevent Assumption (3.4) to hold, since there are few elongated cells and the assumed bound is on a mean quantity.

Remark 5.8* (Lagrangian interpretation).*

The trajectories we construct (displayed in Figure 8) should not be interpreted as realizations of solutions of the stochastic ODE associated with the heat equation. As in remark 5.2, we expect (but do not prove) that our numerical solutions actually approximate the solution to a Lagrangian equation which can be derived from the heat equation, namely

[TABLE]

In contrast with Remark 5.2, the existence of solutions to (5.6) has been established in an article of Evans, Gangbo and Savin [8], assuming that the initial density $\rho_{0}$ is bounded from above and below. Their result can also be extended to some non-linear diffusion equations, under assumptions on the nonlinearity. This Lagrangian point of view has already been used to construct numerical schemes for nonlinear diffusion equations, see [10].

Acknowledgements

This work has been supported by Agence nationale de la recherche (ANR-16-CE40-0014 - MAGA - Monge-Ampère et Géométrie Algorithmique). The authors would like to thank the anonymous referees for several remarks which lead to an improvement of the manuscript.