A variational finite volume scheme for Wasserstein gradient flows

TL;DR

This paper introduces a novel variational finite volume scheme for approximating Wasserstein gradient flows, ensuring energy decay, non-negativity, and convergence, applicable to various energies including the linear Fokker-Planck equation.

Contribution

It develops a new discretization method that preserves the variational structure of Wasserstein gradient flows at the discrete level, with proven convergence and robustness.

Findings

Scheme guarantees non-negativity and energy decay.

Proven convergence for linear Fokker-Planck equation.

Numerical results show first-order accuracy and robustness.

Abstract

We propose a variational finite volume scheme to approximate the solutions to Wasserstein gradient flows. The time discretization is based on an implicit linearization of the Wasserstein distance expressed thanks to Benamou-Brenier formula, whereas space discretization relies on upstream mobility two-point flux approximation finite volumes. Our scheme is based on a first discretize then optimize approach in order to preserve the variational structure of the continuous model at the discrete level. Our scheme can be applied to a wide range of energies, guarantees non-negativity of the discrete solutions as well as decay of the energy. We show that our scheme admits a unique solution whatever the convex energy involved in the continuous problem, and we prove its convergence in the case of the linear Fokker-Planck equation with positive initial density. Numerical illustrations show that it is first order accurate in both time and space, and robust with respect to both the energy and the initial profile.

Figures

a](#S4.F2.sf1) it is noticeable that scheme ([36](#S2.E36)) scheme dissipates the energy faster than the FV one, being indeed a bit more diffusive. This is an expected behavior since the scheme is built to maximize the decrease of the energy and this is actually one of the main strength of the approach. In [Figure 2b](#S4.F2.sf2), one can see that the two dissipations tends to the real one when finer mesh and smaller time step are used, for both schemes, despite the fact that ([36](#S2.E36)) still dissipates faster. In the end, in [Figure 2c](#S4.F2.sf3) it is remarkable that for a very small

Tables

Table 1. Table 1. Time-space convergence for the two schemes. Integration on the time step [ 0 , 0.25 ] . 0 0.25 [0,0.25].

		FV				LJKO
h	dt	${ϵ}_{L^{\infty }}$	$r$	${ϵ}_{L^1}$	$r$	${ϵ}_{L^{\infty }}$	$r$	${ϵ}_{L^1}$	$r$
0.2986	0.0500	0.1634	/	0.0350	/	0.1463	/	0.0334	/
0.1493	0.0250	0.0856	0.932	0.0176	0.997	0.0651	1.169	0.0145	1.120
0.0747	0.0125	0.0434	0.979	0.0087	1.015	0.0449	0.535	0.0066	1.134
0.0373	0.0063	0.0218	0.996	0.0043	1.009	0.0297	0.598	0.0033	1.007
0.0187	0.0031	0.0109	0.999	0.0022	1.004	0.0174	0.770	0.0017	0.943
0.0093	0.0016	0.0054	1.000	0.0011	1.001	0.0095	0.870	0.0009	0.947

Table 2. Table 2. Time-space convergence for scheme ( 36 ). Integration on the time step [ 0.5 , 0.25 ] . 0.5 0.25 [0.5,0.25].

		LJKO
h	dt	${ϵ}_{L^{\infty }}$	$r$	${ϵ}_{L^1}$	$r$
0.2986	0.0500	0.1186	/	0.0216	/
0.1493	0.0250	0.0618	0.9411	0.0109	0.9857
0.0747	0.0125	0.0307	1.0110	0.0053	1.0311
0.0373	0.0063	0.0152	1.0116	0.0026	1.0213
0.0187	0.0031	0.0076	1.0078	0.0013	1.0119
0.0093	0.0016	0.0038	1.0042	0.0006	1.0062

Equations

$$\{\begin{array}{ll}{\partial }_t\varrho -\nabla \cdot (\varrho \nabla \frac{\delta \mathcal{E}}{\delta \rho }[\varrho ])=0& \text{in }{Q}_T=\Omega \times (0,T),\\ \varrho \nabla \frac{\delta \mathcal{E}}{\delta \rho }[\varrho ]\cdot \mathbf{n}=0& \text{on }{\Sigma }_T=\partial \Omega \times (0,T),\\ \varrho (\cdot ,0)={\rho }^0& \text{in}\Omega .\end{array}$$

$$W_2^2(\rho ,\mu )=\underset{\gamma \in \Gamma (\rho ,\mu )}{\min }{\iint }_{\Omega \times \Omega }|\mathbf{y}-\mathbf{x}{|}^{2}d\gamma (\mathbf{x},\mathbf{y}),$$

$$\Gamma(\rho,\mu)=\Big{\{}\gamma\in\mathbb{P}(\Omega\times\Omega):\gamma^{1}=\rho,\gamma^{2}=\mu\Big{\}},$$

$${\partial }_t\varrho =\Delta \varrho +\nabla \cdot (\varrho \nabla V)\text{in}{Q}_T,$$

$$\mathcal{E}(\rho )={\int }_{\Omega }[\rho \log \frac{\rho }{e^{-V}}-\rho +e^{-V}]\mathrm{d}\mathbf{x}.$$

$$\{\begin{array}{l}{\rho }_{\tau }^0={\rho }^0,\\ {\rho }_{\tau }^n\in {\text{argmin}}_{\rho }\frac{1}{2\tau }{W}_2^2(\rho ,{\rho }_{\tau }^{n-1})+\mathcal{E}(\rho ).\end{array}$$

$$\underset{\rho ,\mathbf{v}}{\inf }\frac{1}{2}{\int }_{t^{n-1}}^{t^n}{\int }_{\Omega }\rho |\mathbf{v}{|}^2\mathrm{d}\mathbf{x}\mathrm{d}t+\mathcal{E}(\rho (t^n)),$$

$$\{\begin{array}{ll}{\partial }_t\rho +\nabla \cdot (\rho \mathbf{v})=0& \text{in}\Omega \times (t^{n-1},t^n),\\ \rho \mathbf{v}\cdot \mathbf{n}=0& \text{on }\partial \Omega \times (t^{n-1},t^n),\\ \rho (t^{n-1})={\rho }_{\tau }^{n-1}& \text{in}\Omega .\end{array}$$

$$\begin{gathered} \underset{\rho ,\mathbf{m}}{\inf }\underset{\varphi }{\sup }{\int }_{t^{n-1}}^{t^n}{\int }_{\Omega }\frac{|\mathbf{m}{|}^2}{2\rho }\mathrm{d}\mathbf{x}\mathrm{d}t+{\int }_{t^{n-1}}^{t^n}{\int }_{\Omega }(\rho {\partial }_t\varphi +\mathbf{m}\cdot \nabla \varphi )\mathrm{d}\mathbf{x}\mathrm{d}t \\ +{\int }_{\Omega }[\varphi (t^{n-1}){\rho }_{\tau }^{n-1}-\varphi (t^n)\rho (t^n)]\mathrm{d}\mathbf{x}+\mathcal{E}(\rho (t^n)). \end{gathered}$$

$$\underset{\varphi }{\sup }\underset{\rho }{\inf }{\int }_{t^{n-1}}^{t^n}{\int }_{\Omega }({\partial }_t\varphi -\frac{1}{2}|\nabla \varphi {|}^2)\rho \mathrm{d}\mathbf{x}\mathrm{d}t+{\int }_{\Omega }[\varphi (t^{n-1}){\rho }_{\tau }^{n-1}-\varphi (t^n)\rho (t^n)]\mathrm{d}\mathbf{x}+\mathcal{E}(\rho (t^n)).$$

$$\underset{\varphi (t^{n-1})}{\sup }{\int }_{\Omega }\varphi (t^{n-1}){\rho }_{\tau }^{n-1}\mathrm{d}\mathbf{x}+\underset{\rho (t^n)}{\inf }\left[\mathcal{E}(\rho (t^n))-{\int }_{\Omega }\varphi (t^n)\rho (t^n)\mathrm{d}\mathbf{x}\right],$$

$$\{\begin{array}{ll}-{\partial }_t\varphi +\frac{1}{2}|\nabla \varphi {|}^2\le 0& \text{in}\Omega \times (t^{n-1},t^n),\\ \varphi (t^n)\le \frac{\delta \mathcal{E}}{\delta \rho }[\rho (t^n)]& \text{in}\Omega .\end{array}$$

$$G_{\varphi (t^n)}(\rho )=\mathcal{E}(\rho )-{\int }_{\Omega }\varphi (t^n)\rho \mathrm{d}\mathbf{x}$$

$$={\int }_{\Omega }\varphi (t^n)(\rho -\bar{\rho })\mathrm{d}\mathbf{x}+\mathcal{E}(\bar{\rho })-{\int }_{\Omega }\varphi (t^n)\bar{\rho }\mathrm{d}\mathbf{x}=G_{\varphi (t^n)}(\bar{\rho }).$$

$$\{\begin{array}{l}{\partial }_t\varphi -\frac{1}{2}|\nabla \varphi {|}^2=0,\\ {\partial }_t\rho -\nabla \cdot (\rho \nabla \varphi )=0,\end{array}\text{in}\Omega \times (t^{n-1},t^n),\text{with}\{\begin{array}{l}\rho (t^{n-1})={\rho }_{\tau }^{n-1},\\ \varphi (t^n)=\frac{\delta \mathcal{E}}{\delta \rho }[\rho (t^n)],\end{array}\text{in}\Omega .$$

$${\Vert \rho -\mu \Vert }_{{\dot{H}}_{\rho }^{-1}}=W_2(\rho ,\mu )+o(W_2(\rho ,\mu )),\forall \rho ,\mu \in \mathbb{P}(\Omega ).$$

$${\Vert h\Vert }_{{\dot{H}}_{\rho }^{-1}}=\left\{\underset{\phi }{\sup }{\int }_{\Omega }h\phi \mathrm{d}\mathbf{x}|{\Vert \phi \Vert }_{{\dot{H}}_{\rho }^1}\le 1\right\},\text{with}{\Vert \phi \Vert }_{{\dot{H}}_{\rho }^1}={\left({\int }_{\Omega }\rho |\nabla \phi {|}^2\mathrm{d}\mathbf{x}\right)}^{1/2},$$

$$\{\begin{array}{ll}\rho -\mu -\nabla \cdot (\rho \nabla \psi )=0& \text{in}\Omega ,\\ \nabla \psi \cdot \mathbf{n}=0& \text{on}\partial \Omega .\end{array}$$

$${\int }_{\Omega }(\rho -\mu )\phi \mathrm{d}\mathbf{x}=-{\int }_{\Omega }\nabla \cdot (\rho \nabla \psi )\phi \mathrm{d}\mathbf{x}={\int }_{\Omega }\rho \nabla \psi \cdot \nabla \phi \mathrm{d}\mathbf{x}\le \Vert \psi {\Vert }_{{\dot{H}}_{\rho }^1}\Vert \phi {\Vert }_{{\dot{H}}_{\rho }^1},$$

$${\rho }_{\tau }^n\in \underset{\rho \in \mathbb{P}(\Omega )}{\mathrm{argmin}}\frac{1}{2\tau }{\Vert \rho -{\rho }_{\tau }^{n-1}\Vert }_{{\dot{H}}_{\rho }^{-1}(\Omega )}^2+\mathcal{E}(\rho ),n\ge 1.$$

$$\underset{\rho ,\mathbf{m}}{\inf }{\int }_{\Omega }\frac{|\mathbf{m}{|}^2}{2\tau \rho }\mathrm{d}\mathbf{x}+\mathcal{E}(\rho ),\text{subject to: }\{\begin{array}{rlrl}& \rho -{\rho }_{\tau }^{n-1}+\nabla \cdot \mathbf{m}=0& & \text{in }\Omega ,\\ & \mathbf{m}\cdot \mathbf{n}=0& & \text{on }\partial \Omega .\end{array}$$

$$\underset{\rho ,\mathbf{m}}{\inf }\underset{\varphi }{\sup }{\int }_{\Omega }\frac{|\mathbf{m}{|}^2}{2\tau \rho }\mathrm{d}\mathbf{x}-{\int }_{\Omega }(\rho -{\rho }_{\tau }^{n-1})\varphi \mathrm{d}\mathbf{x}+{\int }_{\Omega }\mathbf{m}\cdot \nabla \varphi \mathrm{d}\mathbf{x}+\mathcal{E}(\rho ),$$

$$\underset{\varphi }{\sup }{\int }_{\Omega }{\rho }_{\tau }^{n-1}\varphi \mathrm{d}\mathbf{x}+\underset{\rho }{\inf }{\int }_{\Omega }(-\varphi -\frac{\tau }{2}|\nabla \varphi {|}^2)\rho \mathrm{d}\mathbf{x}+\mathcal{E}(\rho ).$$

$${\varphi }_{\tau }^n+\frac{\tau }{2}|\nabla {\varphi }_{\tau }^n{|}^2\le \frac{\delta \mathcal{E}}{\delta \rho }[{\rho }_{\tau }^n],$$

$$\{\begin{array}{l}{\varphi }_{\tau }^n+\frac{\tau }{2}|\nabla {\varphi }_{\tau }^n{|}^2=\frac{\delta \mathcal{E}}{\delta \rho }[{\rho }_{\tau }^n],\\ \frac{{\rho }_{\tau }^n-{\rho }_{\tau }^{n-1}}{\tau }-\nabla \cdot ({\rho }_{\tau }^n\nabla {\varphi }_{\tau }^n)=0,\end{array}$$

$$⟨\mathbf{h},\mathbf{\varphi }{⟩}_{\mathcal{T}}=\sum _{K\in \mathcal{T}}{h}_K{\varphi }_K{m}_K,\forall \mathbf{h}={(h_K)}_{K\in \mathcal{T}},\mathbf{\varphi }={({\varphi }_K)}_{K\in \mathcal{T}},$$

$${\Psi }_{\mathcal{T}}^{\ast }(\mathbf{\rho };\mathbf{\varphi })=\frac{1}{2}\sum _{\begin{array}{c}\sigma \in \Sigma \\ \sigma =K|L\end{array}}{a}_{\sigma }{\rho }_{\sigma }{\left({\varphi }_K-{\varphi }_L\right)}^2\ge 0,$$

$${\rho }_{\sigma }=\{\begin{array}{ll}{\rho }_K& \text{if}{\varphi }_K>{\varphi }_L,\\ {\rho }_L& \text{if}{\varphi }_K<{\varphi }_L,\end{array}\forall \sigma =K|L\in \Sigma .$$

$$F_{K\sigma }+F_{L\sigma }=0,\forall \sigma =K|L\in \Sigma .$$

$${\mathbb{F}}_{\mathcal{T}}=\left\{{\boldsymbol{F}}=\left(F_{K\sigma},F_{L\sigma}\right)_{\sigma=K|L\in\Sigma}\in\mathbb{R}^{2\Sigma}\;\middle|\;\text{\eqref{eq:cons_sig} holds}\right\}.$$

$${\mathbb{R}}_0^{\mathcal{T}}=\left\{\mathbf{h}={\left(h_K\right)}_{K\in \mathcal{T}}\in {\mathbb{R}}^{\mathcal{T}}|⟨\mathbf{h},1{⟩}_{\mathcal{T}}=0\right\}$$

Full text

1. A strategy to approximate Wasserstein gradient flows

1.1. Generalities about Wasserstein gradient flows

Given a convex and bounded open subset $\Omega$ of $\mathbb{R}^{d}$, a strictly convex and proper energy functional $\mathcal{E}:L^{1}(\Omega;\mathbb{R}_{+})\to[0,+\infty]$, and given an initial density $\rho^{0}\in L^{1}(\Omega;\mathbb{R}_{+})$ with finite energy, i.e. such that $\mathcal{E}(\rho^{0})<+\infty$, we want to solve problems of the form:

[TABLE]

Equation 1 expresses the continuity equation for a time evolving density $\varrho$, starting from the initial condition $\rho^{0}$, convected by the velocity field $-\nabla\frac{\delta\mathcal{E}}{\delta\rho}[\varrho]$. The mixed boundary condition the system is subjected to represents a no flux condition across the boundary of the domain for the mass: the total mass is therefore preserved.

It is now well understood since the pioneering works of Otto [33, 50, 51] that equations of the form of (1) can be interpreted as the gradient flow in the Wasserstein space w.r.t. the energy $\mathcal{E}$ [2]. A gradient flow is an evolution stemming from an initial condition and evolving at each time following the steepest decreasing direction of a prescribed functional. Consider the space $\mathbb{P}(\Omega)$ of nonnegative measures defined on the bounded and convex domain $\Omega$ with prescribed total mass that are absolutely continuous w.r.t. the Lebesgue measure (hence $\mathbb{P}(\Omega)\subset L^{1}(\Omega;\mathbb{R}_{+})$). The Wasserstein distance $W_{2}$ between two densities $\rho,\mu\in\mathbb{P}(\Omega)$ is the cost to transport one into the other in an optimal way with respect to the cost given by the squared euclidean distance, namely the optimization problem

[TABLE]

with the set $\Gamma(\rho,\mu)$ of admissible transport plans given by

[TABLE]

where $\gamma^{1},\gamma^{2}$ denote the first and second marginal measure, respectively.

A typical example of problem entering the framework of (1) is the linear Fokker-Planck equation

[TABLE]

complemented with no-flux boundary conditions and an initial condition. In (3), $V\in W^{1,\infty}(\Omega)$ denotes a Lipschitz continuous exterior potential. In this case, the energy functional is

[TABLE]

The potential $V$ is defined up to an additive constant, which can be adjusted so that the densities $e^{-V}$ and $\rho^{0}$ have the same mass. Beside this simple example studied for instance in [33, 10], many problems have been proven to exhibit the same variational structure. Porous media flows [51, 37, 15], magnetic fluids [50], supraconductivity [4, 3], crowd motions [45], aggregation processes in biology [22, 9], semiconductor devices modelling [35], or multiphase mixtures [18, 32] are just few examples of problems that can be represented as gradient flows in the Wasserstein space. Designing efficient numerical schemes for approximating their solutions is therefore a major issue and our leading motivation.

1.2. JKO semi-discretization

An intriguing question is how to solve numerically a gradient flow. Problem (1) can of course be directly discretized and solved using one of the many tools available nowadays for the numerical approximation of partial differential equations. The development of energy diminishing numerical methods based on classical ODE solvers for the march in time has been the purpose of many contributions in the recent past, see for instance [8, 16, 17, 13, 54, 49, 19]. Nevertheless, the aforementioned methods disregard the fact that the trajectory aims at optimizing the energy decay, in opposition to methods based on minimizing movement scheme (often called JKO scheme after [33]). This scheme can be thought as a generalization to the space ${\mathbb{P}}(\Omega)$ (the mass being defined by the initial data $\rho^{0}$) equipped with the metric $W_{2}$ of the backward Euler scheme and writes:

[TABLE]

The parameter $\tau$ is the time discretization step. Scheme (5) generates a sequence of measures $\left(\rho_{\tau}^{n}\right)_{n\geq 1}$. Using this sequence it is possible to construct a time dependent measure by gluing them together in a piecewise constant (in time) fashion: $\rho_{\tau}(t)=\rho_{\tau}^{n}$, for $t\in(t^{n-1}=(n-1)\tau,t^{n}=n\tau]$. Under suitable assumptions on the functional $\mathcal{E}$, it is possible to prove the uniform convergence in time of this measure to weak solutions $\varrho$ of (1) (see for instance [2] or [53]).

Lagrangian numerical methods appear to be very natural (especially in dimension 1) to approximate the Wasserstein distance and thus the solution to (5). This was already noticed in [36], and motivated numerous contributions, see for instance [43, 12, 44, 34, 23, 20, 38]. In our approach, we rather consider an Eulerian method based on Finite Volumes for the space discretization. The link between monotone Finite Volumes and optimal transportation was simultaneously highlighted by Mielke [46] and Maas [41, 29, 25, 42, 30]. But these works only focuses on the space discretization, whereas we are interested in the fully discrete setting. Moreover, the approximation based on upstream mobility we propose in Section 2.3 does not enter their framework. Last but not least, let us mention the so-called ALG2-JKO scheme [7, 14] where the optimization problem (5) is discretized and then solved thanks to an augmented Lagrangian iterative method. Our approach is close to the one of [7], with the goal to obtain a faster numerical solver.

Thanks to formal calculations, let us highlight the connection of the minimization problem involved at each step of (5) with a system coupling a forward in time conservation law with a backward in time Hamilton-Jacobi (HJ) equation. The problem can be rewritten thanks to Benamou-Brenier dynamic formulation of optimal transport [6] as

[TABLE]

where the density and velocity curves satisfy weakly

[TABLE]

The next value $\rho_{\tau}^{n}$ is chosen equal to $\rho(t^{n})$ for the optimal $\rho$ in (6)–(7). Using the momentum ${\boldsymbol{m}}=\rho{\boldsymbol{v}}$ instead of ${\boldsymbol{v}}$ as a variable, and incorporating the constraint (7) in (6) yields the saddle-point problem

[TABLE]

We will refer to (8) as the primal problem. The dual problem is obtained by exchanging inf and sup in (8). Strong duality can be proven and the problem hence does not change. Optimizing first w.r.t. ${\boldsymbol{m}}$ leads to ${\boldsymbol{m}}=-\rho{\nabla}\phi$, so that the dual problem writes

[TABLE]

Because of the first term in (9), the infimum is equal to $-\infty$ unless $-\partial_{t}\phi+\frac{1}{2}|{\nabla}\phi|^{2}\leq 0$ a.e. in $\Omega\times(t^{n-1},t^{n})$, with equality $\rho$-almost everywhere since $\rho\geq 0$. Moreover, optimizing w.r.t. $\rho(t^{n})$ provides that $\phi(t^{n})\leq\frac{\delta\mathcal{E}}{\delta\rho}[\rho(t^{n})]$ with equality $\rho(t^{n})$-almost everywhere. Hence the dual problem can be rewritten as

[TABLE]

subject to the constraints

[TABLE]

On the one hand the monotonicity of the backward HJ equation $-\partial_{t}\phi+\frac{1}{2}|{\nabla}\phi|^{2}=f$ with respect to its right-hand side $f\leq 0$ implies that given $\phi(t^{n})$, the solution (which exists) of $-\partial_{t}\phi+\frac{1}{2}|{\nabla}\phi|^{2}=0$ gives a bigger value at $\phi(t^{n-1})$ and thus a better competitor for (11). On the other given $\phi$, among all the $\rho$ satisfying the constraint, the competitor $\bar{\rho}=\frac{\delta\mathcal{E}}{\delta\rho}^{-1}(\phi(t^{n}))$ is optimal. Indeed let $G_{\phi(t^{n})}(\rho)=\mathcal{E}(\rho)-\int_{\Omega}\phi(t^{n})\rho$, the convexity of $\mathcal{E}$ at $\bar{\rho}$ implies

[TABLE]

At the end of the day, the primal-dual optimality conditions of problem (5) finally amounts to the mean field game

[TABLE]

The optimal $\rho_{\tau}^{n}$ of (5) is then equal to $\rho(t^{n})$. The no-flux boundary condition reduces to ${\nabla}\phi\cdot{\boldsymbol{n}}=0$ on $\partial\Omega\times(t^{n-1},t^{n})$.

The approximation of the system (12) is a natural strategy to approximate the solution to (1). This approach was for instance at the basis of the works [7, 21]. These methods require a sub-time stepping to solve system (12) on each interval $(t^{n-1},t^{n})$, yielding a possibly important computational cost. The avoidance of this sub-time stepping is the main motivation of the time discretization we propose now.

1.3. Implicit linearization of the Wasserstein distance and LJKO scheme

Let us introduce in the semi-discrete in time setting the time discretization to be used in the fully discrete setting later on. The following ansatz is at the basis of our approach: when $\tau$ is small, $\rho_{\tau}^{n}$ is close to $\rho_{\tau}^{n-1}$. Then owing to [55, Section 7.6] (see also [52]), the Wasserstein distance between two densities $\rho$ and $\mu$ of ${\mathbb{P}}(\Omega)$ is close to some weighted $H^{-1}$ distance, namely

[TABLE]

In the above formula, we denoted by

[TABLE]

so that $\|\rho-\mu\|_{\dot{H}^{-1}_{\rho}}=\|\psi\|_{\dot{H}^{1}_{\rho}}$ with $\psi$ solution to

[TABLE]

Indeed, in view of (14)–(15), there holds

[TABLE]

with equality if $\varphi=\psi/\|\psi\|_{\dot{H}^{1}_{\rho}}$. Equation (15) can be thought as a linearization of the Monge-Ampère equation.

In view of (13), a natural idea is to replace the Wasserstein distance by the weighted $\dot{H}_{\rho}^{-1}$ norm in (5), leading to what we call the implicitly linearized JKO (or LJKO) scheme:

[TABLE]

The choice of an implicit weight $\rho$ in (16) appears to be particularly important when $\{\rho^{n-1}_{\tau}=0\}$ has a non-empty interior set, which can not be properly invaded by the $\rho_{\tau}^{n}$ if one chooses the explicit (but computationally cheaper) weight $\rho_{\tau}^{n-1}$ as in [48]. Our time discretization is close to the one that was proposed very recently in [40] where the introduction on inner time stepping was also avoided. In [40], the authors introduce a regularisation term based on Fisher information, which mainly amounts to stabilize the scheme thanks to some additional non-degenerate diffusion. In our approach, we manage to avoid this additional stabilization term by taking advantage of the monotonicity of the involved operators.

At each step $n\geq 1$, (16) can be formulated as a constrained optimization problem. To highlight its convexity, we perform the change of variables $(\rho,\psi)\mapsto(\rho,{\boldsymbol{m}}=-\rho{\nabla}\psi)$, in analogy with (6), and rewrite step $n$ as:

[TABLE]

Incorporating the constraint in the above formulation yields the following inf-sup problem:

[TABLE]

the supremum w.r.t. $\phi$ being $+\infty$ unless the constraint is satisfied. Problem (18) is strictly convex in $(\rho,{\boldsymbol{m}})$ and concave (since linear) in $\phi$. Exploiting Fenchel-Rockafellar duality theory it is possible to show that strong duality holds, so that (18) is equivalent to its dual problem where the inf and the sup have been swapped. Optimizing w.r.t. to ${\boldsymbol{m}}$ yields the optimality condition ${\boldsymbol{m}}=-\tau\rho{\nabla}\phi$, hence the problem reduces to

[TABLE]

The problem is now strictly convex in $\rho$ and concave in $\phi$. Optimizing w.r.t. $\rho$ leads to the optimality condition

[TABLE]

with equality on $\{\rho_{\tau}^{n}>0\}$. In the above formula, $\phi_{\tau}^{n}$ denote the optimal $\phi$ realizing the sup in (19). Similarly to was was done in the previous section for the JKO scheme, it is possible to show again that saturating inequality (20) on $\{\rho_{\tau}^{n}=0\}$ is optimal since the mapping $f\mapsto\phi$ solution to $\phi+\frac{\tau}{2}|{\nabla}\phi|^{2}=f$ is monotone. Finally, the optimality conditions for the LJKO problem (16) write

[TABLE]

set on $\Omega$, complemented with homogeneous Neumann boundary condition ${\nabla}\phi_{\tau}^{n}\cdot{\boldsymbol{n}}=0$ on $\partial\Omega$. We can interpret (21) as the one step resolvent of the mean-field game (12). Both the forward in time continuity equation and the backward in time HJ equation are discretized thanks to one step of backward Euler scheme.

1.4. Goal and organisation of the paper

As already noted, most of the numerical methods based on backward Euler scheme disregard the optimal character of the trajectory $t\mapsto\varrho(t)$ of the exact solution to (1). Rather than discretizing directly the PDE (1), which can be thought as the Euler-Lagrange equation for the steepest descent of the energy, we propose to first discretize w.r.t. space the functional appearing in the optimization problem (16), and then to optimize. The corresponding Euler-Lagrange equations will then encode the optimality of the trajectory. The choice of the LJKO scheme (16) rather than the classical JKO scheme (5) is motivated by the fact that solving (21) is computationally affordable. Indeed, it merely demands to approximate two functions $\rho_{\tau}^{n},\phi_{\tau}^{n}$ rather than time depending trajectories in function space as for the JKO scheme (12). This allows in particular to avoid inner time stepping as in [7, 21], making our approach much more tractable to solve complex problems.

Two-Point Flux Approximation (TPFA) Finite Volumes are a natural solution for the space discretization. The are naturally locally conservative thus well-suited to approximate the conservation laws. Moreover, they naturally transpose to the discrete setting the monotonicity properties of the continuous operators. Monotonicity was crucial in the derivation of the optimality conditions (21), as it will also be the case in the fully discrete framework later on. This led us to use upstream mobilities in the definition of the discrete counterpart of the squared $\dot{H}^{1}_{\rho}$ norm. The system (21) thus admits a discrete counterpart (36). The derivation of the fully discrete Finite Volume scheme based on the LJKO time discretization is performed in Section 2, where we also establish the well-posedness of the scheme, as well as the preservation at the discrete level of fundamental properties of the continuous model, namely the non-negativity of the densities and the decay of the energy along time. In Section 3, we show that our scheme converges in the case of the Fokker-Planck equation (3) under the assumption that the initial density is bounded from below by a positive constant. Even though we do not treat problem (1) in its full generality, this results shows for the consistency of the scheme. Finally, Section 4 is devoted to numerical results, where our scheme is tested on several problems including systems.

2. A variational Finite Volume scheme

The goal of this section is to define the fully discrete scheme to solve (1), and to exhibit some important properties of the scheme. But at first, let us give some assumptions and notations on the mesh.

2.1. Discretization of $\Omega$

The domain $\Omega\subset\mathbb{R}^{d}$ is assumed to be polygonal if $d=2$ or polyhedral if $d=3$. The specifications on the mesh are classical for TPFA Finite Volumes [27]. More precisely, an admissible mesh of $\Omega$ is a triplet $\left(\mathcal{T},\overline{\Sigma},{({\boldsymbol{x}}_{K})}_{K\in\mathcal{T}}\right)$ such that the following conditions are fulfilled.

1. (i)

Each control volume (or cell) $K\in\mathcal{T}$ is non-empty, open, polyhedral and convex. We assume that $K\cap L=\emptyset$ if $K,L\in\mathcal{T}$ with $K\neq L$, while $\bigcup_{K\in\mathcal{T}}\overline{K}=\overline{\Omega}$. The Lebesgue measure of $K\in\mathcal{T}$ is denoted by $m_{K}>0$. 2. (ii)

Each face $\sigma\in\overline{\Sigma}$ is closed and is contained in a hyperplane of $\mathbb{R}^{d}$, with positive $(d-1)$-dimensional Hausdorff (or Lebesgue) measure denoted by $m_{\sigma}=\mathcal{H}^{d-1}(\sigma)>0$. We assume that $\mathcal{H}^{d-1}(\sigma\cap\sigma^{\prime})=0$ for $\sigma,\sigma^{\prime}\in\overline{\Sigma}$ unless $\sigma^{\prime}=\sigma$. For all $K\in\mathcal{T}$, we assume that there exists a subset $\overline{\Sigma}_{K}$ of $\overline{\Sigma}$ such that $\partial K=\bigcup_{\sigma\in\overline{\Sigma}_{K}}\sigma$. Moreover, we suppose that $\bigcup_{K\in\mathcal{T}}\overline{\Sigma}_{K}=\overline{\Sigma}$. Given two distinct control volumes $K,L\in\mathcal{T}$, the intersection $\overline{K}\cap\overline{L}$ either reduces to a single face $\sigma\in\overline{\Sigma}$ denoted by $K|L$, or its $(d-1)$-dimensional Hausdorff measure is [math]. 3. (iii)

The cell-centers $({\boldsymbol{x}}_{K})_{K\in\mathcal{T}}\subset\Omega$ are pairwise distinct and are such that, if $K,L\in\mathcal{T}$ share a face $K|L$, then the vector ${\boldsymbol{x}}_{L}-{\boldsymbol{x}}_{K}$ is orthogonal to $K|L$ and has the same orientation as the normal ${\boldsymbol{n}}_{KL}$ to $K|L$ outward w.r.t. $K$.

Cartesian grids, Delaunay triangulations or Voronoï tessellations are typical examples of admissible meshes in the above sense. Since no boundary fluxes appear in our problem, the boundary faces $\Sigma_{\rm ext}=\{\sigma\subset\partial\Omega\}$ are not involved in our computations. Nonzeros fluxes may only occur across internal faces $\sigma\in\Sigma=\overline{\Sigma}\setminus\Sigma_{\rm ext}$. We denote by $\Sigma_{K}=\overline{\Sigma}_{K}\cap\Sigma$ the internal faces belonging to $\partial K$, and by $\mathcal{N}_{K}$ the neighboring cells of $K$, i.e., $\mathcal{N}_{K}=\{L\in\mathcal{T}\;|\;K|L\in\Sigma_{K}\}$. For each internal face $\sigma=K|L\in\Sigma$, we refer to the diamond cell $\Delta_{\sigma}$ as the polyhedron whose edges join ${\boldsymbol{x}}_{K}$ and ${\boldsymbol{x}}_{L}$ to the vertices of $\sigma$. The diamond cell $\Delta_{\sigma}$ is convex if ${\boldsymbol{x}}_{K}\in K$ and ${\boldsymbol{x}}_{L}\in L$. Denoting by $d_{\sigma}=|{\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{L}|$, the measure $m_{\Delta_{\sigma}}$ of $\Delta_{\sigma}$ is then equal to $m_{\sigma}d_{\sigma}/d$, where $d$ stands for the space dimension. The transmissivity of the face $\sigma\in\Sigma$ is defined by $a_{\sigma}=m_{\sigma}/d_{\sigma}$.

The space $\mathbb{R}^{\mathcal{T}}$ is equipped with the scalar product

[TABLE]

which mimics the usual scalar product on $L^{2}(\Omega)$.

2.2. Upstream weighted dissipation potentials

Since the LJKO time discretization presented in Section 1.3 relies on weighted $\dot{H}^{1}_{\rho}$ and $H^{-1}_{\rho}$ norms, we introduce the discrete counterparts to be used in the sequel. As it will be explained in what follows, the upwinding yields to problems to introduce discrete counterparts to the norms. To bypass this difficulty, we adopt a formalism based on dissipation potentials inspired from the one of generalized gradient flows introduced by Mielke in [46]. This framework was used for instance to study the convergence of the semi-discrete in space squareroot Finite Volume approximation of the Fokker-Planck equation, see [31].

Let ${\boldsymbol{\rho}}=\left(\rho_{K}\right)_{K\in\mathcal{T}}\in\mathbb{R}_{+}^{\mathcal{T}}$, and let ${\boldsymbol{\phi}}=\left(\phi_{K}\right)_{K\in\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$, then we define the upstream weighted discrete counterpart of $\frac{1}{2}\|\phi\|_{\dot{H}^{1}_{\rho}}^{2}$ by

[TABLE]

where $\rho_{\sigma}$ denotes the upwind value of ${\boldsymbol{\rho}}$ on $\sigma\in\Sigma$:

[TABLE]

Because of the upwind choice of the mobility (23), the functional (22) is not symmetric, i.e., $\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}};{\boldsymbol{\phi}})\neq\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}};-{\boldsymbol{\phi}})$ in general, which prohibits to define a semi-norm from $\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}};\cdot)$. But one easily checks that ${\boldsymbol{\phi}}\mapsto\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}},{\boldsymbol{\phi}})$ is convex, continuous thus lower semi-continuous (l.s.c.) and proper.

Let us now turn to the definition of the discrete counterpart of $\frac{1}{2}\|\cdot\|^{2}_{\dot{H}^{-1}_{\rho}}$. To this end, we introduce the space ${\mathbb{F}}_{\mathcal{T}}\subset\mathbb{R}^{2\Sigma}$ of conservative fluxes. An element ${\boldsymbol{F}}$ of ${\mathbb{F}}_{\mathcal{T}}$ is made of two outward fluxes $F_{K\sigma},F_{L\sigma}$ for each $\sigma=K|L\in\Sigma$, and one flux $F_{K\sigma}$ per boundary face $\sigma\in\Sigma_{K}$. We impose the conservativity across each internal faces

[TABLE]

In what follows, we denote by $F_{\sigma}=|F_{K\sigma}|=|F_{L,\sigma}|$. There are no flux across the boundary faces. The space ${\mathbb{F}}_{\mathcal{T}}$ is then defined as

[TABLE]

Now, we define the subspace

[TABLE]

and

[TABLE]

where the minimization over ${\boldsymbol{F}}$ is restricted to the linear subspace of ${\mathbb{F}}_{\mathcal{T}}$ such that

[TABLE]

In (25), $\rho_{\sigma}$ denotes the upwind value w.r.t. ${\boldsymbol{F}}$, i.e.,

[TABLE]

In the case where some $\rho_{\sigma}$ vanish, we adopt the following convention in (25) and in what follows:

[TABLE]

Remark that this condition is similar to the one implicitly used in (8) and (17). Summing (26) over $K\in\mathcal{T}$ and using the conservativity across the edges (24), one notices that there is no ${\boldsymbol{F}}\in{\mathbb{F}}_{\mathcal{T}}$ satisfying (26) unless ${\boldsymbol{h}}\in\mathbb{R}^{\mathcal{T}}_{0}$. But when ${\boldsymbol{h}}\in\mathbb{R}^{\mathcal{T}}_{0}$, the minimization set in (25) is never empty. Note that $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})$ may take infinite values when ${\boldsymbol{\rho}}$ vanishes on some cells, for instance $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})=+\infty$ if $h_{K}>0$ and $\rho_{K}=0$ for some $K\in\mathcal{T}$.

Formula (25) deserves some comments. This sum is built to approximate $\int_{\Omega}\frac{|{\boldsymbol{m}}|^{2}}{2\rho}{\rm d}{\boldsymbol{x}}$. The flux $F_{\sigma}$ approximates $|{\boldsymbol{m}}\cdot{\boldsymbol{n}}_{\sigma}|$, and thus encodes the information on ${\boldsymbol{m}}$ only in the one direction (normal to the face $\sigma$) over $d$. But on the other hand, the volume $d_{\sigma}m_{\sigma}$ is equal to $dm_{\Delta_{\sigma}}$ which allows to hope that the sum is a consistent approximation of the integral. This remark has a strong link with the notion of inflated gradients introduced in [24, 26]. The convergence proof carried out in Section 3 somehow shows the non-obvious consistency of this formula.

At the continuous level, the norms $\|\cdot\|_{\dot{H}^{1}_{\rho}}$ and $\|\cdot\|_{H^{-1}_{\rho}}$ are in duality. This property is transposed to the discrete level in the following sense.

Lemma 2.1**.**

Given ${\boldsymbol{\rho}}\geq{\bf 0}$, the functionals ${\boldsymbol{h}}\mapsto\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})$ and ${\boldsymbol{\phi}}\mapsto\Psi^{*}_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\phi}})$ are one another Legendre transforms in the sense that

[TABLE]

In particular, both are proper convex l.s.c. functionals. Moreover, if $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})$ is finite, then there exists a discrete Kantorovitch potential ${\boldsymbol{\phi}}$ solving

[TABLE]

such that

[TABLE]

Proof.

Let ${\boldsymbol{\rho}}\geq{\bf 0}$ be fixed. Incorporating the constraint (26) in (25), and using the definition of $\rho_{\sigma}$ and the twice conservativity constraint (24), we obtain the saddle point primal problem

[TABLE]

The functional in the right-hand side is convex and coercive w.r.t. ${\boldsymbol{F}}$ and linear w.r.t. ${\boldsymbol{\phi}}$, so that strong duality holds. We can exchange the sup and the inf in the above formula to obtain the dual problem, and we minimize first w.r.t. ${\boldsymbol{F}}$, leading to

[TABLE]

Substituting $F_{K\sigma}$ by $\rho_{\sigma}\frac{\phi_{K}-\phi_{L}}{d_{\sigma}}$ in the dual problem leads to (28), while the constraint (26) turns to (29). The fact that $\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}},\cdot)$ is also the Legendre transform of $\Psi_{\mathcal{T}}({\boldsymbol{\rho}},\cdot)$ follows from the fact that it is convex l.s.c., hence equal to its relaxation.

When $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})$ is finite, then the supremum in (28) is achieved, ensuring the existence of the corresponding discrete Kantorovitch potentials ${\boldsymbol{\phi}}$. Finally, multiplying (29) by the optimal $\phi_{K}$ and by summing over $K\in\mathcal{T}$ yields $\langle{\boldsymbol{h}},{\boldsymbol{\phi}}\rangle_{\mathcal{T}}=2\Psi^{*}_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\phi}})$. Substituting this relation in (28) shows the relation $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{h}})=\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}};{\boldsymbol{\phi}})$. ∎

Our next lemma can be seen as an adaptation to our setting of a well known properties of optimal transportation, namely $\rho\mapsto\frac{1}{2}W_{2}^{2}(\rho,\mu)$ is convex, which is key in the study of Wasserstein gradient flows.

Lemma 2.2**.**

Let ${\boldsymbol{\mu}}\in\mathbb{R}^{\mathcal{T}}_{+}$, the function ${\boldsymbol{\rho}}\mapsto\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\mu}}-{\boldsymbol{\rho}})$ is proper and convex on $({\boldsymbol{\mu}}+\mathbb{R}_{0}^{\mathcal{T}})\cap\mathbb{R}_{+}^{\mathcal{T}}.$

Proof.

The function ${\boldsymbol{\rho}}\mapsto\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\mu}}-{\boldsymbol{\rho}})$ is proper since it is equal to [math] at ${\boldsymbol{\rho}}={\boldsymbol{\mu}}$. Then it follows from (28) that

[TABLE]

Since ${\boldsymbol{\rho}}\mapsto\Psi_{\mathcal{T}}^{*}({\boldsymbol{\rho}};{\boldsymbol{\phi}})$ is linear, $\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\mu}}-{\boldsymbol{\rho}})$ is defined as the supremum of linear functions, whence it is convex. ∎

2.3. A variational upstream mobility Finite Volume scheme

The finite volume discretization replaces the functions $\rho_{\tau}^{n},\phi_{\tau}^{n}$ at time step $n\geq 1$ defined on $\Omega$ with the vectors ${\boldsymbol{\rho}}^{n}\in\mathbb{R}_{+}^{\mathcal{T}}$ and ${\boldsymbol{\phi}}^{n}\in\mathbb{R}^{\mathcal{T}}$. In each cell $K$, the restriction of each of these functions is approximated by a single real number $\rho_{K}^{n},\phi_{K}^{n}$, which can be thought as its mean value located in the cell center ${\boldsymbol{x}}_{K}$. Given ${\boldsymbol{\rho}}^{0}\in\mathbb{R}_{+}^{\mathcal{T}}$, the space ${\mathbb{P}}_{\mathcal{T}}$ which is the discrete counterpart of ${\mathbb{P}}(\Omega)$ is then defined by

[TABLE]

It is compact. The energy $\mathcal{E}$ is discretized into a strictly convex functional $\mathcal{E}_{\mathcal{T}}\in C^{1}(\mathbb{R}_{+}^{\mathcal{T}};\mathbb{R}_{+})$ that we do not specify yet. We refer to Sections 3 and 4 for explicit examples.

We have introduced all the necessary material to introduce our numerical scheme, which combines upstream weighted Finite Volumes for the space discretization and the LJKO time discretization:

[TABLE]

A further characterization of the scheme is needed for its practical implementation, but the condensed expression (32) already provides crucial informations gathered in the following theorem. Note in particular that our scheme automatically preserves mass and the positivity since the solutions $\left({\boldsymbol{\rho}}^{n}\right)_{n\geq 1}$ belong to ${\mathbb{P}}_{\mathcal{T}}$.

Theorem 2.3**.**

For all $n\geq 1$, there exists a unique solution ${\boldsymbol{\rho}}^{n}\in{\mathbb{P}}_{\mathcal{T}}$ to (32). Moreover, energy is dissipated along the time steps. More precisely,

[TABLE]

Proof.

The functional ${\boldsymbol{\rho}}\mapsto\frac{1}{\tau}\Psi_{\mathcal{T}}({\boldsymbol{\rho}};{\boldsymbol{\rho}}^{n-1}-{\boldsymbol{\rho}})+\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}})$ l.s.c. and strictly convex on the compact set ${\mathbb{P}}_{\mathcal{T}}$ in view of Lemma 2.2 and of the assumptions on $\mathcal{E}_{\mathcal{T}}$. Moreover, it is proper since ${\boldsymbol{\rho}}^{n-1}$ belongs to its domain. Therefore, it admits a unique minimum on ${\mathbb{P}}_{\mathcal{T}}$. The energy / energy dissipation estimate (33) is obtained by choosing ${\boldsymbol{\rho}}={\boldsymbol{\rho}}^{n-1}$ as a competitor in (32). ∎

In view of (31), and after rescaling the dual variable ${\boldsymbol{\phi}}\leftarrow\frac{{\boldsymbol{\phi}}}{\tau}$, solving (32) amounts to solve the saddle point problem

[TABLE]

which is equivalent to its dual problem

[TABLE]

Our strategy for the practical computation of the solution to (32) is to solve the system corresponding to the optimality conditions of (35). So far, we did not take advantage of the upwind choice of the mobility (23) (we only used the linearity of $({\boldsymbol{\rho}},{\boldsymbol{\phi}})\mapsto\left(\rho_{\sigma}\right)_{\sigma\in\Sigma}$ in the proofs of Lemmas 2.1 and 2.2, which also holds true for a centered choice of the mobilities). The upwinding will be key in the proof of the following theorem, which, roughly speaking, states that there is no need of a Lagrange multiplier for the constraint ${\boldsymbol{\rho}}\geq{\bf 0}$.

Theorem 2.4**.**

The unique solution $({\boldsymbol{\rho}}^{n},{\boldsymbol{\phi}}^{n})$ to system

[TABLE]

where $\rho_{\sigma}^{n}$ denotes the upwind value, i.e.,

[TABLE]

is a saddle point of (35).

System (36) is the discrete counterpart of (21), whose derivation relied on the monotonicity of the inverse of the operator $\phi\mapsto\phi+\frac{\tau}{2}|{\nabla}\phi|^{2}$. Before proving Theorem 2.4, let us show that the space discretization preserves this property at the discrete level. To this end, we introduce the functional ${\boldsymbol{\mathcal{G}}}=\left(\mathcal{G}_{K}\right)_{K}\in C^{1}(\mathbb{R}^{\mathcal{T}};\mathbb{R}^{\mathcal{T}})$ defined by

[TABLE]

Lemma 2.5**.**

Given ${\boldsymbol{f}}\in\mathbb{R}^{\mathcal{T}}$, there exists a unique solution to ${\boldsymbol{\mathcal{G}}}({\boldsymbol{\phi}})={\boldsymbol{f}}$, and it satisfies

[TABLE]

Moreover, let ${\boldsymbol{\phi}},\widetilde{{\boldsymbol{\phi}}}$ be the solutions corresponding to ${\boldsymbol{f}}$ and $\widetilde{{\boldsymbol{f}}}$ respectively, then

[TABLE]

Proof.

Given ${\boldsymbol{f}}\geq\widetilde{{\boldsymbol{f}}}$ and ${\boldsymbol{\phi}},\widetilde{{\boldsymbol{\phi}}}$ corresponding solutions, let $K^{*}$ be the cell such that

[TABLE]

Then, for all the neighboring cells $L$ of $K^{*}$, it holds $\phi_{K^{*}}-\tilde{\phi}_{K^{*}}\leq\phi_{L}-\tilde{\phi}_{L}$ and therefore $\phi_{K^{*}}-\phi_{L}\leq\tilde{\phi}_{K^{*}}-\tilde{\phi}_{L}$ which implies

[TABLE]

Recall ${\boldsymbol{f}}\geq\widetilde{{\boldsymbol{f}}}$ so $\mathcal{G}_{K^{*}}({\boldsymbol{\phi}})\geq\mathcal{G}_{K^{*}}(\tilde{{\boldsymbol{\phi}}})$ together with (39) it yields $\phi_{K^{*}}\geq\tilde{\phi}_{K^{*}}$. Finally as in $K^{*}$ the difference $\phi_{K}-\tilde{\phi}_{K}$ is minimal, we obtain $\phi_{K}\geq\tilde{\phi}_{K}$ for all $K\in\mathcal{T}$. The uniqueness of the solution ${\boldsymbol{\phi}}$ of ${\boldsymbol{\mathcal{G}}}({\boldsymbol{\phi}})={\boldsymbol{f}}$ follows directly. The maximum principle (37) is also a straightforward consequence of (38) as one can compare ${\boldsymbol{\phi}}$ to $(\min{\boldsymbol{f}}){\bf 1}$ and $(\max{\boldsymbol{f}}){\bf 1}$ which are fixed points of ${\boldsymbol{\mathcal{G}}}$. Finally, existence follows from Leray-Schauder fixed-point theorem [39] as the bounds (37) are uniform whatever $\tau\geq 0$. ∎

With Lemma 2.5 at hand, we can now prove Theorem 2.4.

Proof of Theorem 2.4.

Uniqueness of the solution ${\boldsymbol{\rho}}^{n}$ to (32) was already proved in Theorem 2.3. Owing to (33), $\Psi_{\mathcal{T}}({\boldsymbol{\rho}}^{n};{\boldsymbol{\rho}}^{n-1}-{\boldsymbol{\rho}}^{n})$ is finite. So Lemma 2.1 ensures the existence of a discrete Kantorovitch potential ${\boldsymbol{\phi}}^{n}$ satisfying (after a suitable rescaling by $\tau^{-1}$)

[TABLE]

The above condition is the optimality condition w.r.t. ${\boldsymbol{\phi}}$ in (35). The optimality condition w.r.t. ${\boldsymbol{\rho}}$ writes

[TABLE]

where the Lagrange multiplier ${\boldsymbol{\pi}}^{n}=\left(\pi_{K}^{n}\right)_{K}\leq{\bf 0}$ for the constraint ${\boldsymbol{\rho}}\geq{\bf 0}$ satisfies $\langle{\boldsymbol{\pi}}^{n},{\boldsymbol{\rho}}^{n}\rangle_{\mathcal{T}}=0$. To prove Theorem 2.4, we thus has to prove that one can set ${\boldsymbol{\pi}}^{n}={\bf 0}$ and that the solution to system (36) is unique.

Let $({\boldsymbol{\rho}}^{n},\widetilde{{\boldsymbol{\phi}}}^{n})$ be a solution to (40)–(41) with some possibly non-zero Lagrange multiplier ${\boldsymbol{\pi}}^{n}\leq{\bf 0}$. Denote by

[TABLE]

Then we deduce from (40) and from the upstream choice of the mobility that

[TABLE]

Hence the following alternative holds for all the neighbours $L\in\mathcal{N}_{K}$ of $K\in\mathcal{Z}^{n}$:

[TABLE]

Let ${\boldsymbol{\phi}}^{n}$ be such that ${\boldsymbol{\phi}}^{n}\geq\widetilde{{\boldsymbol{\phi}}}^{n}$ with $\phi_{K}^{n}=\widetilde{\phi}_{K}^{n}\;\text{on}\;\mathcal{P}^{n}$, then it follows from the upwinding that (40) holds for all $K\in\mathcal{T}$, while (41) holds for $K\in\mathcal{P}^{n}$.

Reproducing the proof of Lemma 2.5, there exists a unique solution $\left(\overline{\phi}_{K}^{n}\right)_{K\in\mathcal{Z}^{n}}$ to the system

[TABLE]

and, since $\frac{\partial\mathcal{E}_{\mathcal{T}}}{\partial\rho_{K}}({\boldsymbol{\rho}}^{n})\geq\frac{\partial\mathcal{E}_{\mathcal{T}}}{\partial\rho_{K}}({\boldsymbol{\rho}}^{n})+\pi_{K}^{n}$, the monotonicity of ${\boldsymbol{\mathcal{G}}}$ yields $\overline{\phi}_{K}^{n}\geq\widetilde{\phi}_{K}^{n}$ for all $K\in\mathcal{Z}^{n}$. Then define ${\boldsymbol{\phi}}^{n}\in\mathbb{R}^{\mathcal{T}}$ by

[TABLE]

it follows from (42) that $({\boldsymbol{\rho}}^{n},{\boldsymbol{\phi}}^{n})$ fulfills (40). Moreover, it follows from the upwinding in (41) and from (42) that Equation (41) for $K\in\mathcal{P}^{n}$ does not depend on $\left({\boldsymbol{\phi}}_{K}^{n}\right)_{K\in\mathcal{Z}^{n}}$, so that $({\boldsymbol{\rho}}^{n},{\boldsymbol{\phi}}^{n})$ satisfies also (41), but with ${\boldsymbol{\pi}}^{n}={\bf 0}$. Finally, owing to Lemma 2.5, ${\boldsymbol{\phi}}^{n}$ is unique, concluding the proof of Theorem 2.4. ∎

2.4. Comparison with the classical backward Euler discretization

The scheme (32) is based on a “first discretize then optimize” approach. We have built a discrete counterpart of $\frac{1}{2}W_{2}^{2}$ and a discrete energy $\mathcal{E}_{\mathcal{T}}$, then the discrete dynamics is chosen in an optimal way by (32). In opposition, the continuous equation (1) can be thought as the Euler-Lagrange optimality condition for the steepest descent of the energy. A classical approach to approximate the optimal dynamics is to discretize directly (1), leading to what we call a “first optimize then discretize” approach. It is classical for the semi-discretization in time of (1) to use a backward Euler scheme. If one combines this technic with upstream weighted Finite Volumes, we obtain the following fully discrete scheme:

[TABLE]

This scheme has no clear variational structure in the sense that, to our knowledge, $\check{\boldsymbol{\rho}}^{n}$ is no longer the solution to an optimization problem. However, it shares some common features with our scheme (32): it is mass and positivity preserving as well as energy diminishing.

Proposition 2.6**.**

Given ${\boldsymbol{\rho}}^{n-1}\in{\mathbb{P}}_{\mathcal{T}}$, there exists at least one solution $(\check{\boldsymbol{\rho}}^{n},\check{\boldsymbol{\phi}}^{n})\in{\mathbb{P}}_{\mathcal{T}}\times\mathbb{R}^{\mathcal{T}}$ to system (43), which satisfies

[TABLE]

Proof.

Summing (43) over $K\in\mathcal{T}$ provides directly the conservation of mass, i.e., $\langle\check{\boldsymbol{\rho}}^{n},{\bf 1}\rangle_{\mathcal{T}}=\langle{\boldsymbol{\rho}}^{n-1},{\bf 1}\rangle_{\mathcal{T}}$. Assume for contradiction that $\mathcal{K}^{n}=\left\{K\in\mathcal{T}\;\middle|\;\check{\rho}_{K}^{n}<0\right\}\neq\emptyset,$ then choose $K^{\star}\in\mathcal{K}^{n}$ such that $\check{\phi}_{K^{\star}}^{n}\geq\check{\phi}_{K}^{n}$ for all $K\in\mathcal{K}^{n}$. Then it follows from the upwind choice of the mobility in (43) that

[TABLE]

so that $\check{\rho}_{K^{\star}}^{n}\geq\rho_{K^{\star}}^{n-1}\geq 0$, showing a contradiction. Therefore, $\mathcal{K}^{n}=\emptyset$ and $\check{\boldsymbol{\rho}}^{n}\geq{\bf 0}$. These two a priori estimates (mass and positivity preservation) are uniform w.r.t. $\tau\geq 0$, thus they are sufficient to prove the existence of a solution $(\check{\boldsymbol{\rho}}^{n},\check{\boldsymbol{\phi}}^{n})$ to (43) thanks to a topological degree argument [39].

Let us now turn to the derivation of the energy / energy dissipation inequality (44). Multiplying (43) by $\check{\phi}_{K}^{n}$ and summing over $K\in\mathcal{T}$ provides

[TABLE]

The definition of $\check{\boldsymbol{\phi}}^{n}$ and the convexity of $\mathcal{E}_{\mathcal{T}}$ yield $\langle\check{\boldsymbol{\rho}}^{n}-{\boldsymbol{\rho}}^{n-1},\check{\boldsymbol{\phi}}^{n}\rangle_{\mathcal{T}}\geq\mathcal{E}_{\mathcal{T}}(\check{\boldsymbol{\rho}}^{n})-\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{n-1}).$ Thus to prove (44), it remains to check that

[TABLE]

In view of (29), $\tau\check{\boldsymbol{\phi}}^{n}$ is a discrete Kantorovitch potential sending ${\boldsymbol{\rho}}^{n-1}$ on $\check{\boldsymbol{\rho}}^{n}$ for the mobility corresponding to $\check{\boldsymbol{\rho}}^{n}$. Therefore (45) holds as a consequence of (30). ∎

Next proposition provides a finer energy / energy dissipation estimate than (33), which can be thought as discrete counterpart to the energy / energy dissipation inequality (EDI) which is a characterization of generalized gradient flows [2, 46].

Proposition 2.7**.**

Given ${\boldsymbol{\rho}}^{n-1}\in{\mathbb{P}}_{\mathcal{T}}$, let ${\boldsymbol{\rho}}^{n}$ be the unique solution to (32) and let $\check{\boldsymbol{\rho}}^{n}$ be a solution to (43), then

[TABLE]

where $\check{\boldsymbol{\phi}}^{n}$ is defined by $m_{K}\check{\phi}_{K}^{n}=\frac{\partial\mathcal{E}_{\mathcal{T}}}{\partial\rho_{K}}(\check{\boldsymbol{\rho}}^{n})$ for all $K\in\mathcal{T}$.

Proof.

Since $\check{\boldsymbol{\rho}}^{n}$ belongs to ${\mathbb{P}}_{\mathcal{T}}$, it is an admissible competitor for (32), thus

[TABLE]

Combining this with (44) and bearing in mind that $\frac{1}{\tau}\Psi_{\mathcal{T}}({\boldsymbol{\rho}}^{n};{\boldsymbol{\rho}}^{n-1}-{\boldsymbol{\rho}}^{n})=\tau\Psi_{\mathcal{T}}^{*}\left({\boldsymbol{\rho}}^{n};{\boldsymbol{\phi}}^{n}\right)$ thanks to (30), we obtain the desired inequality (46). ∎

3. Convergence in the Fokker-Planck case

In this section, we investigate the limit of the scheme when the time step $\tau$ and the size of the mesh $h_{\mathcal{T}}$ tend to [math] in the specific case of the Fokker-Planck equation (3). The size of the mesh is defined by $h_{\mathcal{T}}=\max_{K\in\mathcal{T}}h_{K}$ with $h_{K}={\rm diam}(K)$. To this end, we consider a sequence $\left(\mathcal{T}_{m},\overline{\Sigma}_{m},\left({\boldsymbol{x}}_{K}\right)_{K\in\mathcal{T}_{m}}\right)_{m\geq 1}$ of admissible discretizations of $\Omega$ in the sense of Section 2.1 and a sequence $\left(\tau_{m}\right)_{m\geq 1}$ of time steps such that $\lim_{m\to\infty}\tau_{m}=\lim_{m\to\infty}h_{\mathcal{T}_{m}}=0.$ We also make the further assumptions on the mesh sequence: there exists $\zeta>0$ such that, for all $m\geq 1$,

[TABLE]

Let $T>0$ be an arbitrary finite time horizon, then we assume for the sake of simplicity that $\tau_{m}=T/N_{m}$ for some integer $N_{m}$ tending to $+\infty$ with $m$. For the ease of reading, we remove the subscript $m\geq 1$ when it appears to be unnecessary for understanding.

Given $V\in C^{2}(\overline{\Omega})$, we define the discrete counterpart of the energy (4) by

[TABLE]

where $V_{K}=V({\boldsymbol{x}}_{K})$ for all $K\in\mathcal{T}$. In view of the above formula, there holds

[TABLE]

Given an initial condition $\varrho^{0}\in{\mathbb{P}}(\Omega)$ with positive mass, i.e. $\int_{\Omega}\varrho^{0}{\rm d}{\boldsymbol{x}}>0$, and such that $\mathcal{E}(\varrho^{0})<\infty$, it is discretized into ${\boldsymbol{\rho}}^{0}=\left(\rho_{K}^{0}\right)_{K\in\mathcal{T}}$ defined by

[TABLE]

Note that the energy $\mathcal{E}_{\mathcal{T}}$ is not in $C^{1}(\mathbb{R}^{\mathcal{T}}_{+})$ since its gradient blows up on $\partial\mathbb{R}_{+}^{\mathcal{T}}$. However, the functional $\mathcal{E}_{\mathcal{T}}$ is continuous and strictly convex on $\mathbb{R}^{\mathcal{T}}_{+}$, hence the scheme (32) still admits a unique solution ${\boldsymbol{\rho}}^{n}$ for all $n\geq 1$ thanks to Theorem 2.3, since its proof does not use the differentiability of the energy. Thanks to the conservativity of the scheme and definition (49) of ${\boldsymbol{\rho}}^{0}$, one has

[TABLE]

Let us show that ${\boldsymbol{\rho}}^{n}>{\bf 0}$ for all $n\geq 1$. To this end, we proceed as in [53, Lemma 8.6].

Lemma 3.1**.**

Assume that $\varrho^{0}$ has positive mass, then the iterated solutions $\left({\boldsymbol{\rho}}^{n}\right)_{n\geq 1}$ to scheme (32) satisfy ${\boldsymbol{\rho}}^{n}>{\bf 0}$ for all $n\geq 1$. Moreover, there exists a unique sequence $\left({\boldsymbol{\phi}}^{n}\right)_{\geq 1}$ of discrete Kantorovitch potentials such that the following optimality conditions are satisfied for all $K\in\mathcal{T}$ and all $n\geq 1$:

[TABLE]

Proof.

Define $\overline{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\varrho^{0}{\rm d}{\boldsymbol{x}}$ and $\overline{{\boldsymbol{\rho}}}=\overline{\rho}{\bf 1}\in{\mathbb{P}}_{\mathcal{T}}$, and by ${\boldsymbol{\rho}}^{n}_{\epsilon}=\left(\rho_{K,\epsilon}^{n}\right)_{K\in\mathcal{T}}=\epsilon\overline{{\boldsymbol{\rho}}}+(1-\epsilon){\boldsymbol{\rho}}^{n}\in{\mathbb{P}}_{\mathcal{T}}$ for some arbitrary $\epsilon\in(0,1)$. Since ${\boldsymbol{\rho}}^{n}$ is optimal in (32), there holds

[TABLE]

The convexity of ${\boldsymbol{\rho}}\mapsto\Psi_{\mathcal{T}}({\boldsymbol{\rho}},{\boldsymbol{\rho}}^{n-1}-{\boldsymbol{\rho}})$ implies that

[TABLE]

while the boundedness of $V$ provides

[TABLE]

Therefore, the right-hand side in (52) can be overestimated by

[TABLE]

for some $C$ depending on ${\boldsymbol{\rho}}^{n},{\boldsymbol{\rho}}^{n-1}$ and $V$ but not on $\epsilon$. Setting $\mathcal{Z}^{n}=\{K\in\mathcal{T}\;|\;\rho_{K}^{n}=0\}$ and $\mathcal{P}^{n}=\{K\in\mathcal{T}\;|\;\rho_{K}^{n}>0\}={(\mathcal{Z}^{n})}^{c}$, we have

[TABLE]

and, thanks to the convexity of $\rho\mapsto\rho\log\rho$ and to the monotonicity of $\rho\mapsto\log\rho$,

[TABLE]

Then dividing by $\epsilon$ and letting $\epsilon$ tend to [math], we obtain that

[TABLE]

which is only possible if $\mathcal{Z}^{n}=\emptyset$, i.e., ${\boldsymbol{\rho}}^{n}>{\bf 0}$. This implies that $\mathcal{E}_{\mathcal{T}}$ is differentiable at ${\boldsymbol{\rho}}^{n}$, hence the optimality conditions (36) hold, which rewrites as (50)–(51) thanks to (48). By the way, the uniqueness of the discrete Kantorovitch potential ${\boldsymbol{\phi}}^{n}$ for all $n\geq 1$ is provided by Theorem 2.4. ∎

Lemma 3.1 allows to define two functions $\rho_{\mathcal{T},\tau}$ and $\phi_{\mathcal{T},\tau}$ by setting

[TABLE]

It follows from the conservativity of the scheme and definition (49) of ${\boldsymbol{\rho}}^{0}$ that

[TABLE]

so that $\rho_{\mathcal{T},\tau}(\cdot,t)$ belongs to ${\mathbb{P}}(\Omega)$ for all $t\in(0,T)$.

The goal of this section is to prove the following theorem.

Theorem 3.2**.**

Assume that $\varrho^{0}\geq\rho_{\star}$ for some $\rho_{\star}\in(0,+\infty)$ and that $\mathcal{E}(\varrho^{0})<+\infty$, and let $\left(\mathcal{T}_{m},\overline{\Sigma}_{m},\left({\boldsymbol{x}}_{K}\right)_{K\in\mathcal{T}_{m}}\right)_{m\geq 1}$ be a sequence of admissible discretizations of $\Omega$ such that $h_{\mathcal{T}_{m}}$ and $\tau_{m}$ tend to [math] while conditions (47) hold. Then up to a subsequence, $\left(\rho_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$ tends in $L^{1}(Q_{T})$ towards a weak solution $\varrho\in L^{\infty}((0,T);L^{1}(\Omega))\cap L^{2}((0,T);W^{1,1}(\Omega))$ of (3) corresponding to the initial data $\varrho^{0}$.

The proof is based on compactness arguments. At first in Section 3.1, we derive some a priori estimates on the discrete solution. These estimates will be used to obtain some compactness on $\rho_{\mathcal{T}_{m},\tau_{m}}$ and $\phi_{\mathcal{T}_{m},\tau_{m}}$ in Section 3.2. Finally, we identify the limit value as a weak solution in Section 3.3.

3.1. Some a priori estimates

First, let us show that if the continuous initial energy $\mathcal{E}(\varrho^{0})$ is bounded, then so does its discrete counterpart $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{0})$.

Lemma 3.3**.**

Given $\varrho^{0}\in{\mathbb{P}}(\Omega)$ such that $\mathcal{E}(\varrho^{0})<+\infty$, and let ${\boldsymbol{\rho}}^{0}$ be defined by (49), then there exists $C_{1}$ depending only on $\Omega$, $V$ and $\varrho^{0}$ (but not on $\mathcal{T}$) such that $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{n})\leq C_{\ref{cte:L2H1_phi_1}}$ for all $n\geq 0$.

Proof.

It follows from (33) that $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{n})\leq\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{0})$ for all $n\geq 1$. Rewriting $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{0})$ as

[TABLE]

with

[TABLE]

we deduce from the definition (49) of ${\boldsymbol{\rho}}^{0}$ and Jensen’s inequality that

[TABLE]

Since $V$ is continuous, there exists $\widetilde{{\boldsymbol{x}}}_{K}\in K$ such that $\int_{K}e^{-V}{\rm d}{\boldsymbol{x}}=m_{K}e^{-V(\widetilde{{\boldsymbol{x}}}_{K})}$. Therefore,

[TABLE]

Similarly, it follows from the mean value theorem that there exists $\check{\boldsymbol{x}}_{K}\in K$ such that $m_{K}V(\check{\boldsymbol{x}}_{K})\rho_{K}^{0}=\int_{K}\varrho^{0}V{\rm d}{\boldsymbol{x}}.$ Hence,

[TABLE]

Combining (54)–(56) in (53) shows that $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{0})\leq\mathcal{E}(\varrho^{0})+C$ for some $C$ depending only on $V$, $\Omega$ and $\varrho^{0}$. ∎

Our next lemma shows that if $\varrho^{0}$ is bounded away from [math], then so does $\rho_{\mathcal{T},\tau}$.

Lemma 3.4**.**

Using the convention $\log(0)=-\infty$, one has

[TABLE]

In particular, if $\varrho^{0}\geq\rho_{\star}$ for some $\rho_{\star}\in(0,+\infty)$, then there exists $\alpha>0$ depending only on $V$ and $\rho_{\star}$ (but not on $\mathcal{T},\tau$ and $n$) such that ${\boldsymbol{\rho}}^{n}\geq\alpha{\bf 1}$ for all $n\geq 1$.

Proof.

It follows directly from (50) that $\log(\rho_{K}^{n})+V_{K}\geq\phi_{K}^{n}$ for all $K\in\mathcal{T}$. Let $K_{\star}\in\mathcal{T}$ be such that $\phi_{K_{\star}}^{n}\leq\phi_{K}^{n}$ for all $K\in\mathcal{T}$, then the conservation equation (51) ensures that $\rho_{K_{\star}}^{n}\geq\rho_{K_{\star}}^{n-1}$. On the other hand, since

[TABLE]

the discrete HJ equation (50) provides that

[TABLE]

Assume now that $\varrho^{0}\geq\rho_{\star}$, then for all $K\in\mathcal{T}$ and all $n\geq 0$,

[TABLE]

Therefore, we obtain the desired inequality with $\alpha=\rho_{\star}e^{-\|V^{+}\|_{\infty}-\|V^{-}\|_{\infty}}$. ∎

Our third lemma deals with some estimates on the discrete gradient of the discrete Kantorovitch potentials $\left({\boldsymbol{\phi}}^{n}\right)_{n}$.

Lemma 3.5**.**

Let $({\boldsymbol{\rho}}^{n},{\boldsymbol{\phi}}^{n})$ be the iterated solution to (36), then

[TABLE]

Moreover, if $\varrho^{0}\geq\rho_{\star}\in(0,+\infty)$, then there exists $C_{2}$ depending on $\Omega$, $V$ and $\varrho^{0})$ such that

[TABLE]

Proof.

Since $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}})\geq 0$ for all ${\boldsymbol{\rho}}\in{\mathbb{P}}_{\mathcal{T}}$, summing (33) over $n\in\{1,\dots,N\}$ yields

[TABLE]

Thanks to (30), the left-hand side rewrites

[TABLE]

so that it only remains to use Lemma 3.3 to recover (57).

Finally, if $\varrho^{0}$ is bounded from below by some $\rho_{\star}>0$, then Lemma 3.4 shows that $\rho_{K}^{n}\geq\alpha$ for some $\alpha$ depending only on $\rho_{\star}$ and $V$. Therefore, since $\rho_{\sigma}^{n}$ is either equal to $\rho_{K}^{n}$ or to $\rho_{L}^{n}$ for $\sigma=K|L\in\Sigma$, then (58) holds with $C_{\ref{cte:L2H1_phi_2}}=\frac{C_{\ref{cte:L2H1_phi_1}}}{\alpha}$. ∎

The discrete solution $\rho_{\mathcal{T},\tau}$ is piecewise constant on the cells. To study the convergence of the scheme, we also need a second reconstruction $\rho_{\Sigma,\tau}$ of the density corresponding to the edge mobilities. It is defined by

[TABLE]

Lemma 3.6**.**

There exists $C_{3}$ depending only on $\zeta$ and $\varrho^{0}$ such that

[TABLE]

Moreover, there exists $C_{4}$ depending only on $\zeta,V$ and $\varrho^{0}$ such that

[TABLE]

Proof.

Since $t\mapsto\rho_{\Sigma,\tau}(\cdot,t)$ is piecewise constant, it suffices to check that the above properties at each $t^{n}$, $1\leq n\leq N$. In view of the definition of $\rho_{\Sigma,\tau}$, one has

[TABLE]

The first term can easily be overestimated by $\int_{\Omega}\rho_{\mathcal{T},\tau}({\boldsymbol{x}},t^{n}){\rm d}{\boldsymbol{x}}=\int_{\Omega}\varrho^{0}{\rm d}{\boldsymbol{x}}$. Since $\rho_{\sigma}^{n}\leq\rho_{K}^{n}+\rho_{L}^{n}$, the second term in the above expression can be overestimated by

[TABLE]

Using the regularity property of the mesh (47c), we obtain that

[TABLE]

so that (60) holds with $C_{\ref{cte:rho_Sig}}=(1+\zeta)\int_{\Omega}\varrho^{0}{\rm d}{\boldsymbol{x}}$.

Reproducing the above calculations, one gets that

[TABLE]

Since $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{n})\leq\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{0})\leq C_{\ref{cte:L2H1_phi_1}}$ and since $V$ is uniformly bounded, we obtain that (61) holds with $C_{\ref{cte:H_Sig}}=(1+\zeta)\left(C_{\ref{cte:L2H1_phi_1}}+\|(1-V)^{+}\|_{\infty}\right)$. ∎

The last lemma of this section can be thought as a discrete $\left(L^{\infty}((0,T);W^{1,\infty}(\Omega))\right)^{\prime}$ estimate on $\partial_{t}\rho_{\mathcal{T},\tau}$. This estimate will be used to apply a discrete nonlinear Aubin-Simon lemma [5] in the next section.

Lemma 3.7**.**

Let $\varphi\in C^{\infty}_{c}(Q_{T})$, then define $\varphi_{K}^{n}=\frac{1}{m_{K}}\int_{K}\varphi({\boldsymbol{x}},t^{n}){\rm d}{\boldsymbol{x}}$ for all $K\in\mathcal{T}$. There exists $C_{5}$ depending only on $\zeta,T,\varrho^{0},d$, such that

[TABLE]

Proof.

Multiplying (51) by $\varphi_{K}^{n}$ and summing over $K\in\mathcal{T}$ and $n\in\{1,\dots,N\}$ yields

[TABLE]

Applying Cauchy-Schwarz inequality on the right-hand side then provides

[TABLE]

The first term in the right-hand side is bounded thanks to Lemma 3.5. On the other hand, the regularity of $\varphi$ ensures that there exists $\widetilde{{\boldsymbol{x}}}_{K}\in K$ such that $\varphi({\boldsymbol{x}}_{K},t^{n})=\varphi_{K}^{n}$ for all $K\in\mathcal{T}$. Thanks to the regularity assumptions (47a)–(47b) on the mesh, there holds

[TABLE]

Hence, the second term of the right-hand side in (62) can be overestimated by

[TABLE]

the last inequality being a consequence of Lemma 3.6. Combining all this material in (62) shows the desired estimate with $C_{\ref{cte:drho_dt}}=(1+2\zeta(1+\zeta))\sqrt{C_{\ref{cte:L2H1_phi_1}}C_{\ref{cte:rho_Sig}}Td}$. ∎

3.2. Compactness of the approximate solution

The goal of this section is to show enough compactness in order to be able to pass to the limit $m\to\infty$. For the sake of readability, we remove the subscript $m$ unless necessary.

Owing to Lemma 3.3, one has $\mathcal{E}_{\mathcal{T}}({\boldsymbol{\rho}}^{n})\leq C_{\ref{cte:L2H1_phi_1}}$ for all $n\in\{1,\dots,N\}$. Proceeding as in the proof of Lemma 3.6, this allows to show that

[TABLE]

for some $C_{6}$ depending only on $\varrho^{0}$, $\zeta$ and $V$. Combining de La Vallée Poussin’s theorem with Dunford-Pettis’ one [56, Ch. XI, Theorem 3.6], there exists $\varrho\in L^{\infty}((0,T);L^{1}(\Omega))$ such that, up to a subsequence,

[TABLE]

Since $\rho\mapsto\rho\log\rho$ is convex, $f\mapsto\iint_{Q_{T}}f\log f{\rm d}{\boldsymbol{x}}{\rm d}t$ is l.s.c. for the weak convergence in $L^{1}(Q_{T})$ (see for instance [11, Corollary 3.9]), so that (63) yields

[TABLE]

Moreover, since $\rho_{\mathcal{T},\tau}\geq\alpha$ thanks to Lemma 3.4, then $\varrho\geq\alpha$ too.

Our goal is to show that $\varrho$ is the unique weak solution to the Fokker-Planck equation (3) corresponding to the initial data $\varrho^{0}$. Even though the continuous problem is linear, (64) is not enough to pass to the limit in our nonlinear scheme. Refined compactness have to be derived in this section so that one can identify $\varrho$ as the solution to (3) in the next section. To show enhanced compactness (and most of all the consistency of the scheme in the next section), we have to assume that the initial data is bounded away from [math].

Proposition 3.8**.**

Assume that $\varrho^{0}\geq\rho_{\star}\in(0,+\infty)$, then, up to a subsequence,

[TABLE]

Proof.

Our proof of (66)–(67) relies on ideas introduced in [47] that were adapted to the discrete setting in [5]. Define the two convex and increasing conjugated functions defined on $\mathbb{R}_{+}$:

[TABLE]

then the following inequality holds for any measurable functions $f,g:Q_{T}\to\mathbb{R}$:

[TABLE]

Now, notice that since $\rho_{\mathcal{T},\tau}$ is bounded from below thanks to Lemma 3.4 and bounded in $L^{1}(Q_{T})$, then $\log\rho_{\mathcal{T},\tau}$ is bounded in $L^{p}(Q_{T})$ for all $p\in[1,\infty)$ and $\Upsilon(|\log(\rho_{\mathcal{T},\tau})|)$ is bounded in $L^{1}(Q_{T})$. As a consequence, there exists $\ell\in L^{\infty}((0,T);L^{p}(\Omega))$ such that

[TABLE]

Since $f\mapsto\iint_{Q_{T}}\Upsilon(|f|)$ is convex thus l.s.c. for the weak convergence, we infer that $\Upsilon(|\ell|)$ belongs to $L^{1}(Q_{T})$. Moreover, in view of (65), $\Upsilon^{*}(\varrho)$ belongs also to $L^{1}(Q_{T})$. Therefore, thanks to (69), the function $\varrho\ell$ is in $L^{1}(Q_{T})$.

Define the quantities

[TABLE]

and by $r_{\mathcal{T},\tau}\in L^{1}(Q_{T})$ the function defined

[TABLE]

Thanks to Lemma 3.5, ${\|r_{\mathcal{T},\tau}\|}_{L^{1}(Q_{T})}\leq\frac{1}{2}C_{\ref{cte:L2H1_phi_2}}\tau$. As a consequence, $r_{\mathcal{T}_{m},\tau_{m}}$ tends to [math] in $L^{1}(Q_{T})$ as $m$ tends to $+\infty$.

Let ${\boldsymbol{\xi}}\in\mathbb{R}^{d}$ be arbitrary, we denote by $\Omega_{\boldsymbol{\xi}}=\{{\boldsymbol{x}}\in\Omega\;|\;{\boldsymbol{x}}+{\boldsymbol{\xi}}\in\Omega\}$. Then using (50) and the triangle inequality, we obtain that for all $m\geq 1$, there holds

[TABLE]

where, denoting by $V_{\mathcal{T}}({\boldsymbol{x}})=V_{K}$ if ${\boldsymbol{x}}\in K$, we have set

[TABLE]

Since $\left(r_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$ and $\left(V_{\mathcal{T}_{m}}\right)_{m\geq 1}$ are compact in $L^{1}(Q_{T})$ and $L^{1}(\Omega)$ respectively, it follows from the Riesz-Frechet-Kolmogorov theorem (see for instance [11, Exercise 4.34]) that there exists $\omega\in C(\mathbb{R}_{+};\mathbb{R}_{+})$ with $\omega(0)=0$ such that

[TABLE]

On the other hand, the function $\phi_{\mathcal{T},\tau}$ belongs to $L^{1}((0,T);BV(\Omega))$ and the integral in time of its total variation in space can be estimated as follows:

[TABLE]

with $C_{7}=\sqrt{d|\Omega|TC_{\ref{cte:L2H1_phi_2}}}$. This implies in particular that $A_{2,m}({\boldsymbol{\xi}})\leq C_{\ref{cte:BV}}|{\boldsymbol{\xi}}|$ for all $m\geq 1$. Combining this estimate with (71) in (50) yields

[TABLE]

The combination of (72) with Lemma 3.7 is exactly what one needs to reproduce the proof of [5, Proposition 3.8], which shows that the product of the weakly convergent sequences ${(\rho_{\mathcal{T}_{m},\tau_{m}})}_{m}$ and ${(\log\rho_{\mathcal{T}_{m},\tau_{m}})}_{m}$ converges towards the product of their weak limits:

[TABLE]

Let us now identify $\ell$ as $\log(\varrho)$ thanks to Minty’s trick. Let $\kappa>0$ and $\varphi\in C^{\infty}_{c}(Q_{T};\mathbb{R}_{+})$ be arbitrary, then thanks to (73),

[TABLE]

As a consequence, $\left(\varrho-\kappa\right)\left(\ell-\log\kappa\right)\geq 0$ a.e. in $Q_{T}$ for all $\kappa>0$, which holds if and only if $\ell=\log\varrho$. To finalize the proof of (66)–(67), define

[TABLE]

Then (73) implies that

[TABLE]

As a consequence, $c_{m}$ tends to [math] almost everywhere in $Q_{T}$, which implies the $\rho_{\mathcal{T}_{m},\tau_{m}}$ tends almost everywhere towards $\varrho$ (up to a subsequence). Then (66)–(67) follow from Vitali’s convergence theorem (see for instance [56, Chap. XI, Theorem 3.9]).

Finally, one has $\phi_{\mathcal{T},\tau}=\log\rho_{\mathcal{T},\tau}+V_{\mathcal{T}}-r_{\mathcal{T},\tau}$. In view of the above discussion, the right-hand side converges strongly in $L^{1}(Q_{T})$ up to a subsequence towards $\log\varrho+V$, then so does the left-hand side. This provides (68) and concludes the proof of Proposition 3.8. ∎

Next lemma shows that $\rho_{\Sigma,\tau}$ shares the same limit $\varrho$ as $\rho_{\mathcal{T},\tau}$.

Lemma 3.9**.**

Assume that $\varrho^{0}\geq\rho_{\star}\in(0,+\infty)$, then

[TABLE]

Proof.

Thanks to Lemma 3.6, it follows from the de La Vallée-Poussin and Dunford Pettis theorems that $\left(\rho_{\Sigma_{m},\tau_{m}}\right)_{m\geq 1}$ is relatively compact for the weak topology of $L^{1}(Q_{T})$. Combining this with (64), we infer that, up to a subsequence, $\left(\rho_{\Sigma_{m},\tau_{m}}-\rho_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$ converges towards some $w$ weakly in $L^{1}(Q_{T})$. Thanks to Vitali’s convergence theorem, it suffices to show that from any subsequence of $\left(\rho_{\Sigma_{m},\tau_{m}}-\rho_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$, one can extract a subsequence that tends to [math] a.e. in $Q_{T}$ (so that the whole sequence converges towards $w=0$), or equivalently

[TABLE]

since both $\left(\rho_{\Sigma_{m},\tau_{m}}\right)_{m\geq 1}$ and $\left(\rho_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$ are bounded away from [math] thanks to Lemma 3.4. Bearing in mind the definition (59) of $\rho_{\Sigma_{m},\tau_{m}}$, and one has

[TABLE]

Using (50) and the triangle inequality, one gets that

[TABLE]

with

[TABLE]

and

[TABLE]

Using again that $dm_{\Delta_{\sigma}}=d_{\sigma}m_{\sigma}\leq\zeta h_{\mathcal{T}}m_{\sigma}$ thanks to (47a), one has

[TABLE]

Since $|r_{K}^{n}-r_{L}^{n}|\leq r_{K}^{n}+r_{L}^{n}$, the regularity assumption (47c) on the mesh implies that

[TABLE]

Since $V$ is Lipschitz continuous, $|V_{K}-V_{L}|\leq\|{\nabla}V\|_{\infty}d_{\sigma}\leq\zeta\|{\nabla}V\|_{\infty}h_{\mathcal{T}}$ for all $\sigma=K|L\in\Sigma$ thanks to (47a). Therefore,

[TABLE]

so that (74) holds, concluding the proof of Lemma 3.9. ∎

3.3. Convergence towards the unique weak solution

Our next lemma is a important step towards the identification of the limit $\varrho$ as a weak solution to the continuous Fokker-Planck equation (3). Define the vector field ${\boldsymbol{F}}_{\Sigma,\tau}:Q_{T}\to\mathbb{R}^{d}$ by

[TABLE]

Lemma 3.10**.**

Assume that $\varrho^{0}\geq\rho_{\star}\in(0,+\infty)$, then, up to a subsequence, the vector field ${\boldsymbol{F}}_{\Sigma_{m},\tau_{m}}$ converges weakly in $L^{1}(Q_{T})^{d}$ towards $-{\nabla}\varrho-\varrho{\nabla}V$ as $m$ tends to $+\infty$. Moreover, $\sqrt{\varrho}$ belongs to $L^{2}((0,T);H^{1}(\Omega))$, while $\varrho$ belongs to $L^{2}((0,T);W^{1,1}(\Omega))$.

Proof.

Let us introduce the inflated discrete gradient ${\boldsymbol{G}}_{\Sigma,\tau}$ of $\phi_{\mathcal{T},\tau}$ defined by

[TABLE]

so that ${\boldsymbol{F}}_{\Sigma,\tau}=-\rho_{\Sigma,\tau}{\boldsymbol{G}}_{\Sigma,\tau}$. Thanks to Lemma 3.5,

[TABLE]

thus we know that, up to a subsequence, ${\boldsymbol{G}}_{\Sigma,\tau}$ converges weakly towards some ${\boldsymbol{G}}$ in $L^{2}(Q_{T})^{d}$ as $m$ tends to $+\infty$. Since $\phi_{\mathcal{T},\tau}$ tends to $\log\varrho+V$, cf. (68), then the weak consistency of the inflated gradient [24, 26] implies that ${\boldsymbol{G}}={\nabla}(\log\varrho+V)$.

Define now ${\boldsymbol{H}}_{\Sigma,\tau}=\sqrt{\rho_{\Sigma,\tau}}{\boldsymbol{G}}_{\Sigma,\tau}$, then using again Lemma 3.5,

[TABLE]

so that there exists ${\boldsymbol{H}}\in L^{2}(Q_{T})^{d}$ such that, up to a subsequence, ${\boldsymbol{H}}_{\Sigma,\tau}$ tends to ${\boldsymbol{H}}$ weakly in $L^{2}(Q_{T})^{d}$. But since $\sqrt{\rho_{\Sigma,\tau}}$ converges strongly towards $\sqrt{\varrho}$ in $L^{2}(Q_{T})$, cf. Lemma 3.6, and since ${\boldsymbol{G}}_{\Sigma,\tau}$ tends weakly towards ${\nabla}(\log\varrho+V)$ in $L^{2}(Q_{T})^{d}$, we deduce that ${\boldsymbol{H}}_{\Sigma,\tau}$ tends weakly in $L^{1}(Q_{T})^{d}$ towards $\sqrt{\varrho}{\nabla}(\log\varrho+V)=2{\nabla}\sqrt{\varrho}+\sqrt{\varrho}{\nabla}V={\boldsymbol{H}}$. In particular, $\sqrt{\varrho}$ belongs to $L^{2}((0,T);H^{1}(\Omega))$. Now, we can pass in the limit $m\to+\infty$ in ${\boldsymbol{F}}_{\Sigma,\tau}=-\sqrt{\rho_{\Sigma,\tau}}{\boldsymbol{H}}_{\Sigma,\tau}$, leading to the desired result. ∎

In order to conclude the proof of Theorem 3.2, it remains to check that any limit value $\varrho$ of the scheme is a solution to the Fokker-Planck equation (3) in the distributional sense.

Proposition 3.11**.**

Let $\varrho$ be a limit value of $\left(\rho_{\mathcal{T}_{m},\tau_{m}}\right)_{m\geq 1}$ as described in Section 3.2, then for all $\varphi\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))$, one has

[TABLE]

Proof.

Given $\varphi\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))$, we denote by $\varphi_{K}^{n}=\varphi({\boldsymbol{x}}_{K},t^{n})$. Then multipying (51) by $-\varphi_{K}^{n-1}$ and summing over $K\in\mathcal{T}$ and $n\in\{1,\dots,N\}$ leads to

[TABLE]

where we have set

[TABLE]

and

[TABLE]

Since $\rho_{\mathcal{T},\tau}$ converges in $L^{1}(Q_{T})$ towards $\varrho$, cf. Proposition 3.8, and since $\varphi$ is smooth,

[TABLE]

It follows from the definition (49) of $\rho_{K}^{0}$ that the piecewise constant function $\rho_{\mathcal{T}}^{0}$, defined by $\rho_{\mathcal{T}}^{0}({\boldsymbol{x}})=\rho_{K}^{0}$ if ${\boldsymbol{x}}\in\mathcal{T}$, converges in $L^{1}(\Omega)$ towards $\varrho^{0}$. Therefore, since $\varphi$ is smooth,

[TABLE]

Let us define

[TABLE]

Then it follows from Lemma 3.10 that

[TABLE]

To conclude the proof of Proposition 3.11, it only remains to check that

[TABLE]

Since $\varphi$ is smooth and since $d_{\sigma}{\boldsymbol{n}}_{KL}={\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{L}$ thanks to the orthogonality condition satisfied by the mesh,

[TABLE]

for some $C_{\varphi}$ depending only on $\varphi$. Therefore,

[TABLE]

Applying Cauchy-Schwarz inequality, one gets that

[TABLE]

thanks to Lemma 3.6. ∎

4. Numerical results

To check the correctness and reliability of our formulation we performed some numerical tests. Before that, we are going to present some details on the solution of the nonlinear system involved in the scheme.

4.1. Newton method

Due to the explicit formulation of the optimality condition of the saddle point problem (35), it appears extremely convient to use a Newton method for their solution. Given ${\boldsymbol{u}}^{n-1}=({\boldsymbol{\phi}}^{n-1},{\boldsymbol{\rho}}^{n-1})\in\mathbb{R}^{2\mathcal{T}}$ solution of the scheme at the time step $n-1$, the Newton method aims at constructing a sequence of approximations of ${\boldsymbol{u}}^{n}$ as ${\boldsymbol{u}}^{n,k+1}={\boldsymbol{u}}^{n,k}+\boldsymbol{d}^{k}$, $\boldsymbol{d}^{k}=(\boldsymbol{d}_{{\boldsymbol{\phi}}}^{k},\boldsymbol{d}_{{\boldsymbol{\rho}}}^{k})$ being the Newton direction, solution to the block-structured system of equations

[TABLE]

In the above linear system, ${\boldsymbol{f}}_{{\boldsymbol{\phi}}}^{k}$ and ${\boldsymbol{f}}_{{\boldsymbol{\rho}}}^{k}$ are the discrete HJ and continuity equations evaluated in ${\boldsymbol{u}}^{n,k}$, and ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\phi}}}^{k}$, ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\rho}}}^{k}$, ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\phi}}}^{k}$ and ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k}$ are the four blocks of the Hessian matrix ${\boldsymbol{J}}^{k}$ of the discrete functional in (35) evaluated in ${\boldsymbol{u}}^{n,k}$. The sequence converges to the unique solution ${\boldsymbol{u}}^{n}$ as soon as the initial guess is sufficiently close to it, which is ensured for a sufficiently small time step by taking ${\boldsymbol{u}}_{0}^{n}={\boldsymbol{u}}^{n-1}$. The algorithm stops when the $\ell^{\infty}$ norm of the discrete equations is smaller than a prescribed tolerance or if the maximum number of iterations is reached. It is possible to implement an adaptative time stepping: if the Newton method converges in few iterations the time step $\tau$ increases; if it reaches the maximum number of iterations the time step is decreased and the method restarted. Issues could arise if the iterate ${\boldsymbol{u}}^{n,k}$ reaches negative values, especially if the energy is not defined for negative densities. To avoid this problem two possible strategies may be implemented: the iterate may be projected on the set of positive measure by taking ${\boldsymbol{u}}^{n,k}=({\boldsymbol{u}}^{n,k})^{+}$; the method may be restarted with a smaller time step.

In case of a local energy functional, as it is the case for the Fokker-Planck and many more examples, the block ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k}$ is diagonal and therefore straightforward to invert. System (76) can be rewritten in term of the Schur complement and solved for $\boldsymbol{d}_{{\boldsymbol{\phi}}}^{k}$ as

[TABLE]

while $\boldsymbol{d}_{{\boldsymbol{\rho}}}^{k}=({\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k})^{-1}\,({\boldsymbol{f}}_{{\boldsymbol{\rho}}}^{k}-{\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\phi}}}^{k}\,\boldsymbol{d}_{{\boldsymbol{\phi}}}^{k})$.

Proposition 4.1**.**

The Schur complement ${\boldsymbol{S}}^{k}={\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\phi}}}^{k}-{\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\rho}}}^{k}\,({\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k})^{-1}\,{\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\phi}}}^{k}$ is symmetric and negative definite.

Proof.

${\boldsymbol{S}}^{k}$ is symmetric since ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\phi}}}^{k}$ and ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k}$ are, while ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\rho}}}^{k}=({\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\phi}}}^{k})^{T}$. The matrix ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k}$ is positive definite since the problem is strictly convex, whereas ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\phi}}}^{k}$ is negative definite if $\rho_{K}^{n,k}>0,\forall K\in\mathcal{T}$, since the problem is strictly concave, but it is semi-negative definite if the density vanishes somewhere. Therefore, it is sufficient to show that the matrix ${\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\rho}}}^{k}=({\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\phi}}}^{k})^{T}={\boldsymbol{M}}+{\boldsymbol{A}}^{k}$ is invertible. ${\boldsymbol{M}}$ is a diagonal matrix such that $({\boldsymbol{M}})_{K,K}=m_{K}$, whereas

[TABLE]

and, for $L\neq K$,

[TABLE]

Therefore the columns of ${\boldsymbol{A}}^{k}$ sum up to [math], so that $({\boldsymbol{J}}_{{\boldsymbol{\phi}},{\boldsymbol{\rho}}}^{k})$ is a column M-matrix [28] and thus invertible. ∎

In case the matrix ${\boldsymbol{J}}_{{\boldsymbol{\rho}},{\boldsymbol{\rho}}}^{k}$ is simple to invert it is then possible to decrease the computational complexity of the solution of system (76). Moreover, it is possible to exploit for the solution of system (77) symmetric solvers, which are computationally more efficient than non-symmetric ones.

4.2. Fokker-Planck equation

We first tackle the gradient flow of the Fokker-Planck energy, namely eq. 3. In section 3 we showed the $L^{1}$ convergence of the scheme. Consider the specific potential $\rho V(\mathbf{x})=-\rho gx$: for this case it is possible to design an analytical solution and test the convergence of the scheme. Consider the domain $\Omega=[0,1]^{2}$, the time interval $[0,0.25]$ and the following analytical solution of the Fokker-Planck equation (built from a one-dimensional one):

[TABLE]

where $\alpha=\pi^{2}+\frac{g^{2}}{4}$. On the domain $\Omega=[0,1]^{2}$, the function $\rho(x,y,t)$ is positive and satisfies the mixed boundary conditions $(\nabla\rho+\rho\nabla V)\cdot{\boldsymbol{n}}|_{\partial\Omega}=0$. We want to exploit the knowledge of this exact solution to compute the error we commit in the spatial and time integration. Consider a sequence of meshes $\left(\mathcal{T}_{m},\overline{\Sigma}_{m},{({\boldsymbol{x}}_{K})}_{K\in\mathcal{T}_{m}}\right)$ with decreasing mesh size $h_{\mathcal{T}_{m}}$ and a sequence of decreasing time steps $\tau_{m}$ such that $\frac{h_{\mathcal{T}_{m+1}}}{h_{\mathcal{T}_{m}}}=\frac{\tau_{m+1}}{\tau_{m}}$. In particular, we used a sequence of Delaunay triangular meshes such that the mesh size halves at each step, obtained subdividing at each step each triangle into four using the edges midpoints. Three subsequent partitioning of the domain are shown in figure 1.

Let us introduce the following mesh-dependent errors:

[TABLE]

where $\rho({\boldsymbol{x}}_{K},n\tau_{m})$ is the value in the cell center of the triangle $K$ of the analytical solution at time $n\tau_{m}$, $n$ running from 0 to the total number of time steps $N_{m}$. The upstream Finite Volume scheme with implicit Euler discretization of the temporal derivative is known to exhibit order one of convergence applied to this problem, both in time and space. This means that the $L^{\infty}((0,T);L^{1}(\Omega))$ and $L^{1}((0,T);L^{1}(\Omega))$ errors halve whenever $h_{\mathcal{T}}$ and $\tau$ halve. We want to inspect whether scheme (36) recovers the same behavior.

For the sequence of meshes and time steps, for $m$ going from one to the total number of meshes, we computed the solution to the linear Fokker-Planck equations and the errors, using both the Finite Volume scheme and scheme (36). The results are shown in Table 1. For each mesh size and time step $m$, it is represented the error together with the rate with respect to the previous one. Scheme (36) exhibits the same order of convergence of the FV scheme. It is noticeable that the rate of convergence of scheme (36) senses a big drop and then recovers order one, especially in the $L^{\infty}((0,T);L^{1}(\Omega))$ error. This is due to the fact that the initial condition $\rho({\boldsymbol{x}}_{K},0)$ is too close to zero, and in particular equal to zero on the set ${1}\times[0,1]$, and scheme (36) tends to be repulsed away from zero due to the singularity of the gradient of the first variation of the energy. In Table 2 we repeated the convergence test for the time interval $[0.05,0.25]$: the convergence profile sensibly improves.

To further investigate and compare the behavior of scheme (36) with the FV, we computed also the energy decay along the trajectory. We call dissipation the difference $\mathcal{E}(\rho)-\mathcal{E}(\rho^{\infty})$, where $\rho^{\infty}$ is the final equilibrium condition, the long time behavior. Since we are discretizing a gradient flow, its dissipation is a useful criteria to assess the goodness of the scheme. The long time value of the energy is equal to:

[TABLE]

It is possible to define the equilibrium solution also on the discrete dynamics on the grid. Namely, the equilibrium solution for the dynamic defined on the grid is

[TABLE]

as it can be easily checked to be the unique minimizer of the discrete energy $\mathcal{E}_{h}=\sum_{K\in\mathcal{T}}E(\rho_{K})m_{K}$ subject to the constraint of the conservation of the mass,

[TABLE]

with $\lambda$ lagrange multiplier associated with the constraint. $M$ is the constant that makes $(\rho_{K}^{\infty})_{K\in\mathcal{T}}$ have the same total mass:

[TABLE]

It is immediate to observe that this is indeed the equilibrium solution in the FV scheme, since with such density the potential in (43) is constant:

[TABLE]

For the scheme (36) instead, as it appears clear from 2.1, whenever $\rho_{K}^{n}=\rho_{K}^{n-1},\forall K\in\mathcal{T}$, as it is the case for an equilibrium solution, the potential is constant. From the potential equation one gets again

[TABLE]

In Figure 2 it is represented the semilog plot of the dissipation of the system, computed for the scheme (36) and the FV one, $\mathcal{E}^{h}(\rho_{K})-\mathcal{E}^{h}(\rho_{K}^{\infty})$, and the real solution, $\mathcal{E}(\rho)-\mathcal{E}(\rho^{\infty})$. In Figure 2a it is noticeable that scheme (36) scheme dissipates the energy faster than the FV one, being indeed a bit more diffusive. This is an expected behavior since the scheme is built to maximize the decrease of the energy and this is actually one of the main strength of the approach. In Figure 2b, one can see that the two dissipations tends to the real one when finer mesh and smaller time step are used, for both schemes, despite the fact that (36) still dissipates faster. In the end, in Figure 2c it is remarkable that for a very small time step the dissipation of the two schemes tends to coincide, as it is expected. For the time parameter going to zero the two schemes coincide.

4.3. Porous medium equation

Consider again the domain $\Omega=[0,1]^{2}$. The porous medium equation,

[TABLE]

has been proven in [51] to be a gradient flow in Wasserstein space with respect to the energy

[TABLE]

for a given $m$ strictly greater than one. Our aim is to show that scheme (36) works regardless of the positivity assumption on the initial measure $\rho_{0}$. For this reason, we use a confining potential $V({\boldsymbol{x}})=\frac{1}{2}||{\boldsymbol{x}}-0.5||^{2}_{2}$. The equilibrium solution of the gradient flow should then be the Barenblatt profile \rho^{\infty}({\boldsymbol{x}})=\big{(}\frac{m-1}{2m}\max(1-||{\boldsymbol{x}}||^{2})\big{)}^{\frac{1}{m-1}}.

In Figure 3 the evolution of an initial density with compact support is shown for the case $m=4$. As expected, the solution converges towards the Barenblatt profile.

4.4. Salinity intrusion problem

We want to show now that scheme (36) can be used for the solutions of systems of equations of the type of eq. 1. We consider the problem of salinity intrusion in an unconfined aquifer. Under the assumption that the two fluids, the fresh and the salty water, are immiscible and the domains occupied by each fluid are separated by a sharp interface, the problem can be modeled via the system of equations

[TABLE]

completed with the no-flux boundary conditions

[TABLE]

and initial conditions $f(t=0)=f_{0},g(t=0)=g_{0}$, with $f_{0},g_{0}\in L^{\infty}(\Omega),f_{0},g_{0}\geq 0$. The quantities $f$, $g$, and $b$ represent respectively the thickness of the fresh layer of water, the thickness of the salty water layer and the height of the bedrock. Therefore the quantity $b+g$ represents the height of the sharp interface separating the two fluids. The parameter $\nu=\frac{\rho_{f}}{\rho_{s}}$ is the ratio between the constant mass density of the fresh and salt water. Equation 78 has been proven in [37] to be a Wasserstein gradient flow with respect to the energy

[TABLE]

Figure 4 represents the bedrock $b$ and the evolution of the surfaces of salt water, $b+g$, and of fresh water, $b+g+f$ (see [1] for a full description of the test case). Also this case is not covered from the theoretical analysis we performed on the convergence of the scheme but still scheme (36) works. As already said, from numerical evidences the scheme works under much more general and mild hypotheses.

Acknowledgements. CC acknowledges the support of the Labex CEMPI (ANR-11-LABX-0007-01). GT acknowledges that this project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754362. We also thanks Guillaume Carlier and Quentin Mérigot for fruitful discussions.