A Lagrangian scheme for the solution of nonlinear diffusion equations using moving simplex meshes

TL;DR

This paper introduces a Lagrangian numerical scheme using moving simplex meshes for solving nonlinear degenerate diffusion equations, effectively capturing solution support growth and maintaining energy estimates.

Contribution

It presents a novel Lagrangian scheme based on gradient flow structure and finite linear maps, with proven energy consistency and applicability to complex nonlinear diffusion equations.

Findings

Scheme effectively tracks support growth in porous medium equations

Energy estimates are preserved in the discrete scheme

Numerical experiments demonstrate scheme's robustness and accuracy

Abstract

A Lagrangian numerical scheme for solving nonlinear degenerate Fokker-Planck equations in space dimensions $d\ge2$ is presented. It applies to a large class of nonlinear diffusion equations, whose dynamics are driven by internal energies and given external potentials, e.g. the porous medium equation and the fast diffusion equation. The key ingredient in our approach is the gradient flow structure of the dynamics. For discretization of the Lagrangian map, we use a finite subspace of linear maps in space and a variational form of the implicit Euler method in time. Thanks to that time discretisation, the fully discrete solution inherits energy estimates from the original gradient flow, and these lead to weak compactness of the trajectories in the continuous limit. Consistency is analyzed in the planar situation, $d=2$. A variety of numerical experiments for the porous medium equation indicates that the scheme is well-adapted to track the growth of the solution's support.

Equations

$${\partial }_t\rho$$

$$\rho (\cdot ,0)$$

$$P(r)=r{h}^{\prime }(r)-h(r)\text{for all r≥0},$$

$$\displaystyle\partial_{t}\rho+\nabla\cdot\big{(}\rho\,\mathbf{v}[\rho]\big{)}=0,$$

$$\displaystyle\mathbf{v}[\rho]=-\nabla\big{(}h^{\prime}(\rho)+V\big{)}.$$

$${\partial }_t{G}_t=\mathbf{v}[{\rho }_t]\circ G_t,G(0,\cdot )=G^0,$$

$${\rho }_t=(G_t{)}_{\#}\overline{\rho }$$

$${\partial }_t{G}_t=-\nabla \left[h^{\prime }\left(\frac{\overline{\rho }}{\det \mathrm{D}{G}_t}\right)\right]\circ G_t-\nabla V\circ G_t.$$

$$\displaystyle\mathcal{E}(\rho)=\int_{\mathbb{R}^{d}}\big{[}h(\rho(x))+V(x)\rho(x)\big{]}\,\mathrm{d}x,$$

$$\mathbf{E}(G|\overline{\rho }):=\mathcal{E}(G_{\#}\overline{\rho })={\int }_K\left[\tilde{h}\left(\frac{\det \mathrm{D}G}{\overline{\rho }}\right)+V\circ G\right]\overline{\rho }\mathrm{d}\omega ,\tilde{h}(s):=sh(s^{-1}),$$

$${\int }_{X}f(x)\mathrm{d}{\mu }_n(x)\to {\int }_{X}f(x)\mathrm{d}\mu (x)$$

$${\int }_{X}g\circ T(x)\mathrm{d}\mu (x)={\int }_{Y}g(y)\mathrm{d}\nu (y).$$

$${\mathcal{P}}_2^{\text{ac}}({\mathbb{R}}^d)=\left\{\rho \in L^1({\mathbb{R}}^d);\rho \ge 0,{\int }_{{\mathbb{R}}^d}\rho (x)\mathrm{d}x=1,{\int }_{{\mathbb{R}}^d}\Vert x{\Vert }^2\rho (x)\mathrm{d}x<\infty \right\}.$$

$$G_{\#}\overline{\rho }=\frac{\overline{\rho }}{\det \mathrm{D}G}\circ G^{-1}.$$

$${\mathrm{W}}_2({\rho }_0,{\rho }_1):=\inf {\left\{{\int }_{{\mathbb{R}}^d}\Vert T(x)-x{\Vert }^2{\rho }_0(x)\mathrm{d}x;T:{\mathbb{R}}^d\to {\mathbb{R}}^d\text{measurable},T_{\#}{\rho }_0={\rho }_1\right\}}^{\frac{1}{2}}.$$

$${\mathrm{W}}_2({\rho }_0,{\rho }_1)=\inf {\left\{{\int }_K\Vert G(\omega )-G_0(\omega ){\Vert }^2\overline{\rho }(\omega )\mathrm{d}\omega ;G:K\to {\mathbb{R}}^d\text{measurable},G_{\#}\overline{\rho }={\rho }_1\right\}}^{\frac{1}{2}},$$

$${\rho }_{\tau }^n:=\underset{\rho \in {\mathcal{P}}_2^{\text{ac}}({\mathbb{R}}^d)}{\mathrm{argmin}}{\mathcal{E}}_{\tau }(\rho ;{\rho }_{\tau }^{n-1}),{\mathcal{E}}_{\tau }(\rho ,\hat{\rho }):=\frac{1}{2\tau }{\mathrm{W}}_2(\rho ,\hat{\rho }{)}^2+\mathcal{E}(\rho ).$$

$$\displaystyle G_{\tau}^{n}:=\operatorname*{argmin}_{G\in L^{2}(K\to{\mathbb{R}^{d}};\overline{\rho})}\mathbf{E}_{\tau}\big{(}G;G_{\tau}^{n-1}\big{)},\quad\mathbf{E}_{\tau}(G;\hat{G})=\frac{1}{2\tau}\int_{K}\|G-\hat{G}\|^{2}\,\,\mathrm{d}\overline{\rho}+\mathbf{E}(G|\overline{\rho}).$$

$${\overline{\rho }}_{\mathcal{T}}^m:=\frac{{\mu }_{\mathcal{T}}^m}{|{\Delta }_m|}\text{for the simplex masses}{\mu }_{\mathcal{T}}^m:={\int }_{{\Delta }_m}\overline{\rho }(\omega )\mathrm{d}\omega .$$

$$G(\omega )=A_m\omega +b_m\text{for all }\omega \in {\Delta }_m,$$

$${△}^d:=\left\{\xi =({\xi }_1,\dots ,{\xi }_d)\in {\mathbb{R}}_{\ge 0}^d;\sum _{j=1}^d{\xi }_j\le 1\right\},$$

$$r_m(\xi )$$

$$q_m(\xi )$$

$$\displaystyle\det A_{m}=\frac{\det\mathrm{D}q_{m}}{\det\mathrm{D}r_{m}}=\frac{\det Q_{\mathscr{T}}^{m}[G]}{2|\Delta_{m}|}\quad\text{where}\quad Q_{\mathscr{T}}^{m}[G]:=\big{(}G_{m,1}-G_{m,0}\big{|}\cdots\big{|}G_{m,d}-G_{m,0}\big{)}.$$

$$\displaystyle\mathbf{E}(G|\overline{\rho}_{\mathscr{T}})=\sum_{\Delta_{m}\in\mathscr{T}}\mu_{\mathscr{T}}^{m}\big{[}\mathbb{H}_{\mathscr{T}}^{m}(G)+\mathbb{V}_{\mathscr{T}}^{m}(G)\big{]}$$

$${\mathbb{H}}_{\mathcal{T}}^m(G):=\tilde{h}\left(\frac{\det A_m}{{\overline{\rho }}_{\mathcal{T}}^m}\right)=\tilde{h}\left(\frac{\det Q_{\mathcal{T}}^m[G]}{2{\mu }_{\mathcal{T}}^m}\right)$$

$$\displaystyle\mathbb{V}_{\mathscr{T}}^{m}(G)=\fint_{\Delta_{m}}V(A_{m}\omega+b_{m})\,\mathrm{d}\omega=\fint_{\mathord{\bigtriangleup}}V\big{(}r_{m}(\omega)\big{)}\,\mathrm{d}\omega.$$

$$\Vert G-G^{\ast }{\Vert }_{L^2(K;{\overline{\rho }}_{\mathcal{T}})}^2={\int }_K\Vert G-G^{\ast }{\Vert }^2{\overline{\rho }}_{\mathcal{T}}\mathrm{d}\omega =\sum _{{\Delta }_m\in \mathcal{T}}{\mu }_{\mathcal{T}}^m{\mathbb{L}}_{\mathcal{T}}^m(G,G^{\ast }).$$

$${\mathbb{L}}_{\mathcal{T}}^m(G,G^{\ast })$$

$$={\text{fint}}_{△}\Vert r_m(\omega )-r_m^{\ast }(\omega ){\Vert }^2\mathrm{d}\omega$$

$$=\frac{2}{(d+1)(d+2)}\sum _{0\le i\le j\le d}(G_{m,i}-G_{m,i}^{\ast })\cdot (G_{m,j}-G_{m,j}^{\ast }).$$

$$\displaystyle G_{\boxplus}^{n}:=\operatorname*{argmin}_{G\in\mathcal{A}_{\mathscr{T}}}\mathbf{E}_{\boxplus}\big{(}G;G_{\boxplus}^{n-1}\big{)}\quad\text{where}\quad\mathbf{E}_{\boxplus}(G;G^{*})=\frac{1}{2\tau}\|G-G^{*}\|_{L^{2}(K;\overline{\rho}_{\mathscr{T}})}^{2}+\mathbf{E}(G|\overline{\rho}_{\mathscr{T}}).$$

$$S_c:=\left\{G\in {\mathcal{A}}_{\mathcal{T}};{\mathbf{E}}_{⊞}(G;G^{\ast })\le c\right\}\text{with}c:=\mathbf{E}(G^{\ast }|{\overline{\rho }}_{\mathcal{T}})$$

Full text

DM was supported by the DFG Collaborative Research Center TRR 109, “Discretization in Geometry and Dynamics”. BD and DSMcC were supported by the Leverhulme Trust research project grant “Novel discretisations for higher-order nonlinear PDE” (RPG-2015-69).

1. Introduction

We study a Lagrangian discretization of the following type of initial value problem for a nonlinear Fokker–Planck equation:

[TABLE]

This problem is posed for the time-dependent probability density function $\rho\colon\mathbb{R}_{\geq 0}\times\mathbb{R}^{d}\to\mathbb{R}_{\geq 0}$, with initial condition $\rho^{0}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$. We assume that the pressure $P\colon\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ can be written in the form

[TABLE]

for some non-negative and convex $h\in C^{1}(\mathbb{R}_{\geq 0})\cap C^{\infty}(\mathbb{R}_{>0})$ and that $V\in C^{2}(\mathbb{R}^{d})$ is a non-negative potential without loss of generality.

Problem (1.1) encompasses a large class of diffusion equations, such as the heat equation ($P(r)=r,V=0$), the porous medium equation ($P(r)=r^{m}/(m-1),m>1,V=0$) and the fast diffusion equation ($P(r)=r^{m}/(m-1),m<1,V=0$), and extends to related problems with given external potentials $V$. To motivate our discretization, we first briefly recall the Lagrangian form of the dynamics: rewriting (1.1) as a transport equation, we obtain

[TABLE]

More general systems can be written in this form. For instance, interaction potentials leading to aggregation equations, relativistic heat equations, $p$-Laplacian equations and Keller-Segel type models can also be included; see Carrillo and Moll [12] and the references therein for a good account of models enjoying this structural form. Here, we reduce to models with nonlinear degenerate diffusion, i.e. $h(0)=h^{\prime}(0)=0$, and confining potentials in order to explore our new discretization.

The system (1.3) naturally induces a Lagrangian representation of the dynamics, which can be summarized as follows. Below, we use the notation $G_{\#}\overline{\rho}$ for the push-forward of $\overline{\rho}$ under a map $G\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$; the definition is recalled in (2.1).

Lemma 1.1**.**

Assume that $\rho\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}_{\geq 0}$ is a smooth positive solution of (1.1). Let $\overline{\rho}\colon\mathbb{R}^{d}\to\mathbb{R}_{\geq 0}$ be a given reference density, and let $G^{0}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ be a given map such that $G^{0}_{\#}\overline{\rho}=\rho^{0}$. Further, let $G\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the flow map associated to (1.3b), satisfying

[TABLE]

where $\rho_{t}:=\rho(t,\cdot)$ and $G_{t}:=G(t,\cdot)\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$. Then

[TABLE]

at any $t\in[0,T]$.

In short, the solution $G$ to (1.4) is a Lagrangian map for the solution $\rho$ to (1.1). This fact is an immediate consequence of (1.3a); for convenience of the reader, we recall the proof in Appendix A. Notice that (1.5) can be substituted for $\rho$ in the expression (1.3b) for the velocity, which makes (1.4) an autonomous evolution equation for $G$:

[TABLE]

A more explicit form of (1.6) is derived in (5.2).

There is a striking structural relation between (1.1) and (1.6): it is well-known (see Otto [30] or Ambrosio, Gigli and Savaré [1]) that (1.1) is a gradient flow for the relative Renyi entropy functional

[TABLE]

with respect to the $L^{2}$-Wasserstein metric on the space $\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ of probability densities on ${\mathbb{R}^{d}}$. It appears to be less well known (see Evans et al. [18], Carrillo and Moll [12], or Carrillo and Lisini [11]) that also (1.6) is a gradient flow, namely for the functional

[TABLE]

on the Hilbert space $L^{2}(K\to{\mathbb{R}^{d}};\overline{\rho})$ of maps from $K$ to ${\mathbb{R}^{d}}$, where $\overline{\rho}$ is a reference measure supported on $K\subset\mathbb{R}^{d}$. We shall discuss these gradient flow structures in more detail in Section 2 below.

The particular spatio-temporal discretization of the initial value problem (1.1) that we study in this paper is based on these facts. Instead of numerically integrating (1.1a) to obtain the density $\rho$ directly, we approximate the Lagrangian map $G$ — using a finite subspace of linear maps in space — which then allows us to recover $\rho$ a posteriori via (1.5). Our approximation for $G$ is constructed by solving a succession of minimization problems that are naturally derived from the gradient flow structure behind (1.4). Via variational methods, this provides compactness estimates on the discrete solutions.

The use of a finite subspace of linear maps for the Lagrangian maps has a geometric interpretation: the induced densities are piecewise constant on triangles whose vertices move in time. In the sequel, we further assume that $\lim_{s\to\infty}sh^{\prime\prime}(s)=+\infty$ in order to prevent the collapse of the images of our Lagrangian maps $G_{t}$, as proven below. However, we do not yet know how to prevent (globally) the possible intersection of the images of the approximated Lagrangian maps.

This approach is alternative to the one developed by Carrillo et al. [14, 12], where $G$ is obtained by directly solving the PDE (1.6) numerically with finite differences or Galerkin approximation via finite element methods. In fact, the main difference between these two strategies can be summarized as follows: while Carrillo et al. [14, 12] follows the strategy minimize first then discretize, our present approach is to discretize first then minimize. In other words, in Carrillo et al. [14, 12] one minimizes first, obtaining as Euler-Lagrange equations the implicit Euler discretization of (1.6) in [14] approximated by the explicit Euler method in [12], and then discretized in space. In the present approach, we discretize first approximating the space of Lagrangian maps by a suitable finite subspace of linear maps, and then we minimize obtaining a nonlinear system of equations to find the approximated Lagrangian map within that set of linear maps.

Let us mention that other numerical methods have been developed to conserve particular properties of solutions of the gradient flow (1.1). Finite volume methods preserving the decay of energy at the semi-discrete level, along with other important properties like non-negativity and mass conservation, were proposed in the papers [4, 9]. Particle methods based on suitable regularisations of the flux of the continuity equation (1.1) have been proposed in the papers [16, 23, 24, 32]. A particle method based on the steepest descent of a regularized internal part of the energy $\mathcal{E}$ in (1.7) by substituting particles by non-overlapping blobs was proposed and analysed in Carrillo et al. [10, 13]. Moreover, the numerical approximation of the JKO variational scheme has already been tackled by different methods using pseudo-inverse distributions in one dimension (see [5, 8, 20, 35]) or solving for the optimal map in a JKO step (see [3, 22]). Finally, note that gradient-flow-based Lagrangian methods in one dimension for higher-order, drift diffusion and Fokker-Planck equations have recently been proposed in the papers [17, 27, 28, 29].

There are two main arguments in favour of our taking this indirect approach of solving (1.6) instead of solving (1.1). The first is our interest in structure preserving discretizations: the scheme that we present builds on the non-obvious “secondary” gradient flow representation of (1.1) in terms of Lagrangian maps. The benefits are monotonicity of the transformed entropy functional $\mathbf{E}$ and an $L^{2}$-control on the metric velocity for our fully discrete solutions, that eventually lead to weak compactness of the trajectories in the continuous limit. We remark that our long-term goal is to design a numerical scheme that makes full use of the much richer “primary” variational structure of (1.1) in the Wasserstein distance, that is reviewed in Section 2 below. However, despite significant effort in the recent past — see, e.g., the references [3, 6, 13, 14, 17, 19, 22, 25, 31, 35] — it has not been possible so far to preserve features like metric contractivity of the flow under the discretization, except in the rather special situation of one space dimension (see Matthes and Osberger [25]). This is mainly due to the non-existence of finite-dimensional submanifolds of $\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ that are complete with respect to generalized geodesics.

The second motivation is that Lagrangian schemes are a natural choice for numerical front tracking, see, e.g., Budd [7] for first results on the numerical approximation of self-similar solutions to the porous medium equation. We recall that due to the assumed degeneracy $P^{\prime}(0)=0$ of the diffusion in (1.1), solutions that are compactly supported initially remain compactly supported at all times. A numerically accurate calculation of the moving edge of support is challenging, since the solution can have a very complex behavior near that edge, like the waiting time phenomenon (see Vazquez [33]). Our simulation results for $\partial_{t}\rho=\Delta(\rho^{3})$ — that possess an analytically known, compactly supported, self-similar Barenblatt solution — indicated that our discretization is indeed able to track the edge of support quite accurately.

This work is organized as follows. In Section 2 we present an overview of previous results in gradient flows pertaining our work. Section 3 is devoted to the introduction of the linear set of Lagragian maps and the derivation of the numerical scheme. Section 4 shows the compactness of the approximated sequences of discretizations and we give conditions leading to the eventual convergence of the scheme towards (1.1). Section 5 deals with the consistency of the scheme in two dimensions while Section 6 gives several numerical tests showing the performance of this scheme.

2. Gradient flow structures

2.1. Notations from probability theory

$\mathcal{P}(X)$ is the space of probability measures on a given base set $X$. We say that a sequence $(\mu_{n})$ of measures in $\mathcal{P}(X)$ converges narrowly to a limit $\mu$ in that space if

[TABLE]

for all bounded and continuous functions $f\in C^{0}_{b}(X)$. The push-forward $T_{\#}\mu$ of a measure $\mu\in\mathcal{P}(X)$ under a measurable map $T\colon X\to Y$ is the uniquely determined measure $\nu\in\mathcal{P}(Y)$ such that, for all $g\in C^{0}_{b}(Y)$,

[TABLE]

With a slight abuse of notation — identifying absolutely continuous measures with their densities — we denote the space of probability densities on ${\mathbb{R}^{d}}$ of finite second moment by

[TABLE]

Clearly, the reference density $\overline{\rho}$, which is supported on the compact set $K\subset\mathbb{R}^{d}$, belongs to $\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$. If $G\colon K\to\mathbb{R}^{d}$ is a diffeomorphism onto its image (which is again compact), then the push-forward of $\overline{\rho}$’s measure produces again a density $G_{\#}\overline{\rho}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$, given by

[TABLE]

2.2. Gradient flow in the Wasserstein metric

Below, some basic facts about the Wasserstein metric and the formulation of (1.1) as gradient flow in that metric are briefly reviewed. For more detailed information, we refer the reader to the monographs of Ambrosio et al. [1] and Villani [34].

One of the many equivalent ways to define the $L^{2}$-Wasserstein distance between $\rho_{0},\rho_{1}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ is as follows:

[TABLE]

The infimum above is in fact a minimum, and the — essentially unique — optimal map $T^{*}$ is characterized by Brenier’s criterion; see, e.g., Villani [34, Section 2.1]. A trivial but essential observation is that if $\overline{\rho}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ is a reference density with support $K\subset{\mathbb{R}^{d}}$, and $\rho_{0}=(G_{0})_{\#}\overline{\rho}$ with a measurable $G_{0}\colon K\to{\mathbb{R}^{d}}$, then (2.2) can be re-written as follows:

[TABLE]

and the essentially unique minimizer $G^{*}$ in (2.3) is related to the optimal map $T^{*}$ in (2.2) via $G^{*}=T^{*}\circ G_{0}$.

$\mathrm{W}_{2}$ is a metric on $\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$; convergence in $\mathrm{W}_{2}$ is equivalent to weak-$\star$ convergence in $L^{1}({\mathbb{R}^{d}})$ and convergence of the second moment. Since $P$ and hence also $h$ are of super-linear growth at infinity, each sublevel set $\mathcal{E}$ is weak-$\star$ closed and thus complete with respect to $\mathrm{W}_{2}$.

As already mentioned above, solutions $\rho$ to (1.1) constitute a gradient flow for the functional $\mathcal{E}$ from (1.7) in the metric space $(\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d});\mathrm{W}_{2})$. In fact, the flow is even $\lambda$-contractive as a semi-group, thanks to the $\lambda$-uniform displacement convexity of $\mathcal{E}$ (see McCann [26], or Daneri and Savaré [15]), which is a strengthened form of $\lambda$-uniform convexity along geodesics. The $\lambda$-contractivity of the flow implies various properties (see Ambrosio et al. [1, Section 11.2]) like global existence, uniqueness and regularity of the flow, monotonicity of $\mathcal{E}$ and its sub-differential, uniform exponential estimates on the convergence (if $\lambda<0$) or divergence (if $\lambda\geq 0$) of trajectories, quantified exponential rates for the approach to equilibrium (if $\lambda<0$) and the like.

An important further consequence is that the unique flow can be obtained as the limit for $\tau\searrow 0$ of the time-discrete minimizing movement scheme (see Ambrosio et al. [1] and Jordan, Kinderlehrer and Otto [21]):

[TABLE]

This time discretization is well-adapted to approximate $\lambda$-contractive gradient flows. All of the properties of mentioned above are already reflected on the level of these time-discrete solutions.

2.3. Gradient flow in $L^{2}$

Equation (1.6) is the gradient flow of $\mathbf{E}$ on the space $L^{2}(K\to{\mathbb{R}^{d}};\overline{\rho})$ of square integrable (with respect to $\overline{\rho}$) maps $G\colon K\to{\mathbb{R}^{d}}$ (see Evans et al. [18] or Jordan et al. [22]). However, the variational structure behind this gradient flow is much weaker than above: most notably, $\mathbf{E}$ is only poly-convex, but not $\lambda$-uniformly convex. Therefore, the abstract machinery for $\lambda$-contractive gradient flows in Ambrosio et al. [1] does not apply here. Clearly, by equivalence of (1.1) and (1.6) at least for sufficiently smooth solutions, certain properties of the primary gradient flow are necessarily inherited by this secondary flow, but for instance $\lambda$-contractivity of the flow in the $L^{2}$-norm seems to fail.

Nevertheless, it can be proven (see Ambrosio, Lisini and Savaré [2]) that the gradient flow is globally well-defined, and it can again be approximated by the minimizing movement scheme:

[TABLE]

In fact, there is an equivalence between (2.5) and (2.4): simply substitute $(G_{\tau}^{n-1})_{\#}\overline{\rho}$ for $\rho_{\tau}^{n-1}$ and $G_{\#}\overline{\rho}$ for $\rho$ in (2.4); notice that any $\rho\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ can be written as $G_{\#}\overline{\rho}$ with a suitable (highly non-unique) choice of $G\in L^{2}(K\to{\mathbb{R}^{d}};\overline{\rho})$. This equivalence was already exploited in Carrillo et al. [14, 12]. Thanks to the equality (2.3), the minimization with respect to $\rho=G_{\#}\overline{\rho}$ can be relaxed to a minimization with respect to $G$. Consequently, if $(G_{\tau}^{0})_{\#}\overline{\rho}=\rho_{\tau}^{0}$, then $(G_{\tau}^{n})_{\#}\overline{\rho}=\rho_{\tau}^{n}$ at all discrete times $n=1,2,\ldots$. However, while the functional $\mathcal{E}_{\tau}(\cdot;\rho_{\tau}^{n-1})$ in (2.4) is $(\lambda+\tau^{-1})$-uniformly convex in $\rho$ along geodesics in $\mathrm{W}_{2}$, the functional $\mathbf{E}_{\tau}(\cdot;G_{\tau}^{n-1})$ in (2.5) has apparently no useful convexity properties in $G$ on $L^{2}(K\to{\mathbb{R}^{d}};\overline{\rho})$.

3. Definition of the numerical scheme

Recall the Lagrangian formulation of (1.1) that has been given in Lemma 1.1. For definiteness, fix a reference density $\overline{\rho}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$, whose support $K\subset{\mathbb{R}^{d}}$ is convex and compact.

3.1. Discretization in space

Our spatial discretization is performed using a finite subspace of linear maps for the Lagrangian maps $G$. More specifically: let $\mathscr{T}$ be some (finite) simplicial decomposition of $K$ with nodes $\omega_{1}$ to $\omega_{L}$ and $n$-simplices $\Delta_{1}$ to $\Delta_{M}$. In the case $d=2$, which is of primary interest here, $\mathscr{T}$ is a triangulation, with triangles $\Delta_{m}$. The reference density $\overline{\rho}$ is approximated by a density $\overline{\rho}_{\mathscr{T}}\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ that is piecewise constant on the simplices of $\mathscr{T}$, with respective values

[TABLE]

The finite dimensional ansatz space $\mathcal{A}_{\mathscr{T}}$ is now defined as the set of maps $G\colon K\to{\mathbb{R}^{d}}$ that are globally continous, affine on each of the simplices $\Delta_{m}\in\mathscr{T}$, and orientation preserving. That is, on each $\Delta_{m}\subset\mathscr{T}$, the map $G\in\mathcal{A}_{\mathscr{T}}$ can be written as follows:

[TABLE]

with a suitable matrix $A_{m}\in\mathbb{R}^{d\times d}$ of positive determinant and a vector $b_{m}\in\mathbb{R}^{d}$.

For the calculations that follow, we shall use a more geometric way to describe the maps $G\in\mathcal{A}_{\mathscr{T}}$, namely by the positions $G_{\ell}=G(\omega_{\ell})$ of the images of each node $\omega_{\ell}$. Denote by $(\mathbb{R}^{d})_{\mathscr{T}}^{L}\subset\mathbb{R}^{L\cdot d}$ the space of $L$-tuples $\vec{G}=(G_{\ell})_{\ell=1}^{L}$ of points $G_{\ell}\in{\mathbb{R}^{d}}$ with the same simplicial combinatorics (including orientation) as the $\omega_{\ell}$ in $\mathscr{T}$. Clearly, any $G\in\mathcal{A}_{\mathscr{T}}$ is uniquely characterized by the $L$-tuple $\vec{G}$ of its values, and moreover, any $\vec{G}\in(\mathbb{R}^{d})_{\mathscr{T}}^{L}$ defines a $G\in\mathcal{A}_{\mathscr{T}}$.

More explicitly, fix a $\Delta_{m}\in\mathscr{T}$, with nodes labelled $\omega_{m,0}$ to $\omega_{m,d}$ in some orientation preserving order, and respective image points $G_{m,0}$ to $G_{m,d}$. With the standard $d$-simplex given by

[TABLE]

introduce the linear interpolation maps $r_{m}\colon\mathord{\bigtriangleup}^{d}\to K$ and $q_{m}\colon\mathord{\bigtriangleup}^{d}\to{\mathbb{R}^{d}}$ by

[TABLE]

Then the affine map (3.2) equals to $q_{m}\circ r_{m}^{-1}$. In particular, we obtain that

[TABLE]

For later reference, we give a more explicit representation for the transformed entropy $\mathbf{E}$ for $G\in\mathcal{A}_{\mathscr{T}}$, and for the $L^{2}$-distance between two maps $G,\hat{G}\in\mathcal{A}_{\mathscr{T}}$. Substitution of the special form (3.2) into (1.8) produces

[TABLE]

with the internal energy (recall the definition of $\widetilde{h}$ from (1.8))

[TABLE]

and the potential energy

[TABLE]

For the $L^{2}$-difference of $G$ and $G^{*}$, we have

[TABLE]

Using Lemma B.1, we obtain on each simplex $\Delta_{m}$:

[TABLE]

3.2. Discretization in time

Let a time step $\tau>0$ be given; in the following, we symbolize the spatio-temporal discretization by $\boxplus$, and we write $\boxplus\to 0$ for the joint limit of $\tau\to 0$ and vanishing mesh size in $\mathscr{T}$.

The discretization in time is performed in accordance with (2.5): we modify $\mathbf{E}_{\tau}$ from (2.5) by restriction to the ansatz space $\mathcal{A}_{\mathscr{T}}$. This leads to the minimization problem

[TABLE]

For a fixed discretization $\boxplus$, the fully discrete scheme is well-posed in the sense that for a given initial map $G_{\boxplus}^{0}\in\mathcal{A}_{\mathscr{T}}$, an associated sequence $(G_{\boxplus}^{n})_{n\geq 0}$ can be determined by successive solution of the minimization problems (3.7). One only needs to verify:

Lemma 3.1**.**

For each given $G^{*}\in\mathcal{A}_{\mathscr{T}}$, there exists at least one global minimizer $G\in\mathcal{A}_{\mathscr{T}}$ of $\mathbf{E}_{\boxplus}(\cdot;G^{*})$.

Remark 3.2**.**

We do not claim uniqueness of the minimizers. Unfortunately, the minimization problem (3.7) inherits the lack of convexity from (2.5), whereas the correspondence between (2.5) and the convex problem (2.4) is lost under spatial discretization. A detailed discussion of $\mathbf{E}_{\boxplus}$’s (non-)convexity is provided in Appendix C.

Proof of Lemma 3.1.

We only sketch the main arguments. For definiteness, let us choose (just for this proof) one of the infinitely many equivalent norm-induced metrics on the $dL$-dimensional vector space $V_{\mathscr{T}}$ of all continuous maps $G\colon K\to\mathbb{R}^{d}$ that are piecewise affine with respect to the fixed simplicial decomposition $\mathscr{T}$: given $G,G^{\prime}\in V_{\mathscr{T}}$ with their respective point locations $\vec{G},\vec{G}^{\prime}\in\mathbb{R}^{dL}$, i.e., $\vec{G}=(G_{\ell})_{\ell=1}^{L}$ for $G_{\ell}=G(\omega_{\ell})$, define the distance between these maps as the maximal $\mathbb{R}^{d}$-distance $\|G_{\ell}-G^{\prime}_{\ell}\|$ of corresponding points $G_{\ell}\in\vec{G}$, $G^{\prime}_{\ell}\in\vec{G}^{\prime}$. Clearly, this metric makes $V_{\mathscr{T}}$ a complete space.

It is easily seen that the subset $\mathcal{A}_{\mathscr{T}}$ — which is singled out by requiring orientation preservation of the $G$’s — is an open subset of $V_{\mathscr{T}}$. It is further obvious that the map $G\mapsto\mathbf{E}_{\boxplus}(G;G^{*})$ is continuous with respect to the metric. The claim of the lemma thus follows if we can show that the sub-level

[TABLE]

is a non-empty compact subset of $V_{\mathscr{T}}$. Clearly, $G^{*}\in S_{c}$, so it suffices to verify compactness.

$S_{c}$ is bounded. We are going to show that there is a radius $R>0$ such that no $G\in S_{c}$ has a distance larger than $R$ to $G^{*}$. From non-negativity of $\mathbf{E}$, and from the representations (3.5) and (3.6), it follows that

[TABLE]

where $\underline{\mu}_{\mathscr{T}}=\min_{\Delta_{m}}\mu_{\mathscr{T}}^{m}$. It is now easy to compute a suitable value for the radius $R$.

$S_{c}$ is a closed subset of $V_{\mathscr{T}}$. It suffices to show that the limit $\overline{G}\in V_{\mathscr{T}}$ of any sequence $(G^{(k)})_{k=1}^{\infty}$ of maps $G^{(k)}\in S_{c}$ belongs to $\mathcal{A}_{\mathscr{T}}$. By definition of our metric on $V_{\mathscr{T}}$, global continuity and piecewise linearity of the $G^{(k)}$ trivially pass to the limit $\overline{G}$. We still need to verify that $\overline{G}$ is orientation-preserving. Fix a simplex $\Delta_{m}$ and consider the corresponding matrices $A_{m}^{(k)}$ and $\overline{A}_{m}$ from (3.2). Since the $G^{(k)}$ converge to $\overline{G}$ in the metric, also $A_{m}^{(k)}\to\overline{A}_{m}$ entry-wise. Now, by non-negativity of $\widetilde{h}$, we have for all $k$ that

[TABLE]

and since $\widetilde{h}(s)\to+\infty$ as $s\downarrow 0$, it follows that $\det A^{(k)}_{m}>0$ is bounded away from zero, uniformly in $k$. But then also $\det\overline{A}_{m}>0$, i.e., the $m$th linear map piece of the limit $\overline{G}$ preserves orientation. ∎

3.3. Fully discrete equations

We shall now derive the Euler-Lagrange equations associated to the minimization problem (3.7), i.e., for each given $G^{*}:=G_{\boxplus}^{n-1}\in\mathcal{A}_{\mathscr{T}}$, we calculate the variations of $\mathbf{E}_{\boxplus}(G;G^{*})$ with respect to the degrees of freedom of $G\in\mathcal{A}_{\mathscr{T}}$. Since that function is a weighted sum over the triangles $\Delta_{m}\in\mathscr{T}$, it suffices to perform the calculations for one fixed triangle $\Delta_{m}$, with respective nodes $\omega_{m,0}$ to $\omega_{m,d}$, in positive orientation. The associated image points are $G_{m,0}$ to $G_{m,d}$. Since we may choose any vertex to be labelled $\omega_{m,0}$, it will suffice to perform the calculations at one fixed image point $G_{m,0}$.

• •

mass term:

[TABLE]

• •

internal energy: observing that — recall (1.2) —

[TABLE]

we obtain

[TABLE]

where

[TABLE]

is the uniquely determined vector in $\mathbb{R}^{d}$ that is orthogonal to the $(d-1)$-simplex with corners $G_{m,1}$ to $G_{m,d}$ (pointing away from $G_{m,0}$) and whose length equals the $(d-1)$-volume of that simplex.

• •

potential energy:

[TABLE]

Now let $\omega_{\ell}$ be a fixed vertex of $\mathscr{T}$. Summing over all simplices $\Delta_{m}$ that have $\omega_{\ell}$ as a vertex, and choosing vertex labels in accordance with above, i.e., such that $\omega_{m,0}=\omega_{\ell}$ in $\Delta_{m}$, produces the following Euler-Lagrange equation:

[TABLE]

3.4. Approximation of the initial condition

For the approximation $\rho^{0}_{\boxplus}=(G^{0}_{\boxplus})_{\#}\overline{\rho}_{\mathscr{T}}$ of the initial datum $\rho^{0}=G^{0}_{\#}\overline{\rho}$, we require:

• •

$\rho^{0}_{\boxplus}$ converges to $\rho^{0}$ narrowly;

• •

$\mathcal{E}(\rho^{0}_{\boxplus})$ is $\boxplus$-uniformly bounded, i.e.,

[TABLE]

In our numerical experiments, we always choose $\overline{\rho}:=\rho^{0}$, in which case $G^{0}\colon K\to\mathbb{R}^{d}$ can be taken as the identity on $K$, and we choose accordingly $G^{0}_{\boxplus}$ as the identity as well. Hence $\rho^{0}_{\boxplus}=\overline{\rho}_{\mathscr{T}}$, which converges to $\rho^{0}=\overline{\rho}$ even strongly in $L^{1}(K)$. Moreover, since $h$ is convex, it easily follows from Jensen’s inequality that

[TABLE]

and therefore,

[TABLE]

In more general situations, in which $G^{0}$ is not the identity, a sequence of approximations $G^{0}_{\boxplus}$ of $G^{0}$ is needed. Pointwise convergence $G^{0}_{\boxplus}\to G^{0}$ is more than sufficient to guarantee narrow convergence of $\rho_{\boxplus}^{0}$ to $\rho^{0}$, but the uniform bound (3.11) might require a well-adapted approximation, especially for non-smooth $G^{0}$’s.

4. Limit trajectory

In this section, we assume that a sequence of vanishing discretizations $\boxplus\to 0$ is given, and we study the respective limit of the fully discrete solutions $(G_{\boxplus}^{n})_{n\geq 0}$ that are produced by the inductive minimization procedure (3.7). For the analysis of that limit trajectory, it is more natural to work with the induced densities and velocities,

[TABLE]

instead of the Lagrangian maps $G_{\boxplus}^{n}$ themselves. Note that $\mathbf{v}_{\boxplus}^{n}$ is only well-defined on the support of $\rho_{\boxplus}^{n}$ — that is, on the image of $G_{\boxplus}^{n}$ — and can be assigned arbitrary values outside. Let us introduce the piecewise constant in time interpolations $\widetilde{\rho}_{\boxplus}\colon[0,T]\times{\mathbb{R}^{d}}\to\mathbb{R}_{\geq 0}$, and $\widetilde{\mathbf{v}}_{\boxplus}\colon[0,T]\times{\mathbb{R}^{d}}\to\mathbb{R}^{d}$ as usual,

[TABLE]

Note that $\widetilde{\rho}(t,\cdot)\in\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ and $\widetilde{\mathbf{v}}_{\boxplus}(t,\cdot)\in L^{2}({\mathbb{R}^{d}}\to\mathbb{R}^{d};\widetilde{\rho}_{\boxplus}(t,\cdot))$ at each $t\geq 0$.

4.1. Energy estimates

We start by proving the classical energy estimates on minimizing movements for our fully discrete scheme.

Lemma 4.1**.**

For each discretization $\boxplus$ and for any indices $\overline{n}>\underline{n}\geq 0$, one has the a priori estimate

[TABLE]

Consequently:

1. (1)

$\mathbf{E}$* is monotonically decreasing, i.e., $\mathcal{E}(\widetilde{\rho}_{\boxplus}(t))\leq\mathcal{E}(\widetilde{\rho}_{\boxplus}(s))$ for all $t\geq s\geq 0$;* 2. (2)

$\widetilde{\rho}_{\boxplus}$* is Hölder-$\frac{1}{2}$-continuous in $\mathrm{W}_{2}$, up to an error $\tau$,*

[TABLE] 3. (3)

$\widetilde{\mathbf{v}}_{\boxplus}$* is square integrable with respect to $\widetilde{\rho}_{\boxplus}$,*

[TABLE]

Proof.

By the definition of $G_{\boxplus}^{n}$ as a minimizer, we know that $\mathbf{E}_{\boxplus}(G_{\boxplus}^{n};G_{\boxplus}^{n-1})\leq\mathbf{E}_{\boxplus}(G;G_{\boxplus}^{n-1})$ for any $G\in\mathcal{A}_{\mathscr{T}}$, and in particular for the choice $G:=G_{\boxplus}^{n-1}$, which yields:

[TABLE]

Summing these inequalies for $n=\underline{n}+1,\ldots,\overline{n}$, recalling that $\mathcal{E}(\rho_{\boxplus}^{n})=\mathbf{E}(G_{\boxplus}^{n}|\overline{\rho}_{\mathscr{T}})$ by (1.8) and that $\mathrm{W}_{2}(\rho_{\boxplus}^{n},\rho_{\boxplus}^{n-1})^{2}\leq\int_{K}|G_{\boxplus}^{n}-G_{\boxplus}^{n-1}|^{2}\overline{\rho}\,\mathrm{d}\omega$ by (2.3), produces (4.1).

Monotonicity of $\mathcal{E}$ in time is obvious.

To prove (4.2), choose $\underline{n}\leq\overline{n}$ such that $s\in((\underline{n}-1)\tau,\underline{n}\tau]$ and $t\in((\overline{n}-1)\tau,\overline{n}\tau]$. Notice that $\tau(\overline{n}-\underline{n})\leq t-s+\tau$. If $\underline{n}=\overline{n}$, the claim (4.2) is obviously true; let $\underline{n}<\overline{n}$ in the following. Combining the triangle inequality for the metric $\mathrm{W}_{2}$, estimate (4.1) above and Hölder’s inequality for sums, we arrive at

[TABLE]

Finally, changing variables using $x=G_{\boxplus}^{n}(\omega)$ in (4.4) yields

[TABLE]

and summing these inequalities from $n=1$ to $n=N_{\tau}$ yields (4.3). ∎

4.2. Compactness of the trajectories and weak formulation

Our main result on the weak limit of $\widetilde{\rho}_{\boxplus}$ is the following.

Theorem 4.2**.**

Along a suitable sequence $\boxplus\to 0$, the curves $\widetilde{\rho}_{\boxplus}\colon\mathbb{R}_{\geq 0}\to\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$ convergence pointwise in time, i.e., $\widetilde{\rho}_{\boxplus}(t)\to\rho_{*}(t)$ narrowly for each $t>0$, towards a Hölder-$\frac{1}{2}$-continuous limit trajectory $\rho_{*}\colon\mathbb{R}_{\geq 0}\to\mathcal{P}_{2}^{\text{ac}}(\mathbb{R}^{d})$.

Moreover, the discrete velocities $\widetilde{\mathbf{v}}_{\boxplus}$ possess a limit $\mathbf{v}_{*}\in L^{2}(\mathbb{R}_{\geq 0}\times{\mathbb{R}^{d}};\rho_{*})$ such that $\widetilde{\mathbf{v}}_{\boxplus}\widetilde{\rho}_{\boxplus}\overset{*}{\rightharpoonup}\mathbf{v}_{*}\rho_{*}$ in $L^{1}(\mathbb{R}_{\geq 0}\times{\mathbb{R}^{d}})$, and the continuity equation

[TABLE]

holds in the sense of distributions.

Remark 4.3**.**

The Hölder continuity of $\rho_{*}$ implies that $\rho_{*}$ satisfies the initial condition (1.1b) in the sense that $\rho_{*}(t)\to\rho^{0}$ narrowly as $t\downarrow 0$.

Proof of Theorem 4.2.

We closely follow an argument that is part of the general convergence proof for the minimizing movement scheme as given in Ambrosio et al. [1, Section 11.1.3]. Below, convergence is shown for some arbitrary but fixed time horizon $T>0$; a standard diagonal argument implies convergence at arbitrary times.

First observe that by estimate (4.2) — applied with $0=s\leq t\leq T$ — it follows that $\mathrm{W}_{2}(\widetilde{\rho}_{\boxplus}(t),\rho_{\boxplus}^{0})$ is bounded, uniformly in $t\in[0,T]$ and in $\boxplus$. Since further $\rho^{0}_{\boxplus}$ converges narrowly to $\rho^{0}$ by our hypotheses on the initial approximation, we conclude that all densities $\widetilde{\rho}_{\boxplus}(t)$ belong to a sequentially compact subset for the narrow convergence. The second observation is that the term on the right hand side of (4.2) simplifies to $(2\overline{\mathcal{E}})^{\frac{1}{2}}|t-s|^{\frac{1}{2}}$ in the limit $\boxplus\to 0$. A straightforward application of the “refined version” of the Ascoli-Arzelà theorem (Proposition 3.3.1 in Ambrosio et al. [1]) yields the first part of the claim, namely the pointwise narrow convergence of $\widetilde{\rho}_{\boxplus}$ towards a Hölder continuous limit curve $\rho_{*}$.

It remains to pass to the limit with the velocity $\widetilde{\mathbf{v}}_{\boxplus}$. Towards that end, we define a probability measure $\widetilde{\gamma}_{\boxplus}\in\mathcal{P}(Z_{T})$ on the set $Z_{T}:=[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ as follows:

[TABLE]

for every bounded and continuous function $\varphi\in C^{0}_{b}(Z_{T})$. For brevity, let $\widetilde{M}_{\boxplus}\in\mathcal{P}([0,T]\times\mathbb{R}^{d})$ be the $(t,x)$-marginals of $\widetilde{\gamma}_{\boxplus}$, that have respective Lebesgue densities $\frac{\rho_{\boxplus}(t,x)}{T}$ on $[0,T]\times\mathbb{R}^{d}$. Thanks to the result from the first part of the proof, $\widetilde{M}_{\boxplus}$ converges narrowly to a limit $M_{*}$, which has density $\frac{\rho_{*}(t,x)}{T}$. On the other hand, the estimate (4.3) implies that

[TABLE]

We are thus in the position to apply Theorem 5.4.4 in Ambrosio et al. [1], which yields the narrow convergence of $\widetilde{\gamma}_{\boxplus}$ towards a limit $\gamma_{*}$. Clearly, the $(t,x)$-marginal of $\gamma_{*}$ is $M_{*}$. Accordingly, we introduce the disintegration $\gamma_{(t,x)}$ of $\gamma_{*}$ with respect to $M_{*}$, which is well-defined $M_{*}$-a.e.. Below, it will turn out that $\gamma_{*}$’s $v$-barycenter,

[TABLE]

is the sought-for weak limit of $\widetilde{\mathbf{v}}_{\boxplus}$. The convergence $\widetilde{\mathbf{v}}_{\boxplus}\widetilde{\rho}_{\boxplus}\overset{*}{\rightharpoonup}\mathbf{v}_{*}\rho_{*}$ and the inheritance of the uniform $L^{2}$-bound (4.3) to the limit $\mathbf{v}_{*}$ are further direct consequences of Theorem 5.4.4 in Ambrosio et al. [1].

The key step to establish the continuity equation for the just-defined $\mathbf{v}_{*}$ is to evaluate the limit as $\boxplus\to 0$ of

[TABLE]

for any given test function $\phi\in C^{\infty}_{c}((0,T)\times\mathbb{R}^{d})$ in two different ways. First, we change variables $t\mapsto t+\tau$ in the second integral, which gives

[TABLE]

For the second evaluation, we write

[TABLE]

and substitute accordingly $x\mapsto x-\tau\widetilde{\mathbf{v}}_{\boxplus}(t,x)$ in the second integral, leading to

[TABLE]

The error term $\mathfrak{e}_{\boxplus}[\phi]$ above is controlled via Taylor expansion of $\phi$ and by using (4.3),

[TABLE]

Equality of the limits for both evaluations of $J_{\boxplus}[\phi]$ for arbitrary test functions $\phi$ shows the continuity equation (4.5). ∎

Unfortunately, the convergence provided by Theorem 4.2 is generally not sufficient to conclude that $\rho_{*}$ is a weak solution to (1.1), since we are not able to identify $\mathbf{v}_{*}$ as $\mathbf{v}[\rho_{*}]$ from (1.3b). The problem is two-fold: first, weak-$\star$ convergence of $\widetilde{\rho}_{\boxplus}$ is insufficient to pass to the limit inside the nonlinear function $P$. Second, even if we would know that, for instance, $P(\widetilde{\rho}_{\boxplus})\overset{*}{\rightharpoonup}P(\rho_{*})$, we would still need a $\boxplus$-independent a priori control on the regularity (e.g., maximal diameter of triangles) of the meshes generated by the $G_{\boxplus}^{n}$ to justify the passage to limit in the weak formulation below.

The main difficulty in the weak formulation that we derive now is that we can only use “test functions” that are piecewise affine with respect to the changing meshes generated by the $G_{\boxplus}^{n}$. For definiteness, we introduce the space

[TABLE]

Lemma 4.4**.**

Assume $S\colon{\mathbb{R}^{d}}\to\mathbb{R}^{d}$ is such that $S\circ G_{\boxplus}^{n}\in\mathcal{D}(\mathscr{T})$. Then:

[TABLE]

Proof.

For all sufficiently small $\varepsilon>0$, let $G_{\varepsilon}=(\mathrm{id}+S)\circ G_{\boxplus}^{n}$. By definition of $G_{\boxplus}^{n}$ as a minimizer, we have that $\mathbf{E}_{\boxplus}(G_{\varepsilon};G_{\boxplus}^{n-1})\geq\mathbf{E}_{\boxplus}(G_{\boxplus}^{n};G_{\boxplus}^{n-1})$. This implies that

[TABLE]

We discuss limits of the three terms under the integral for $\varepsilon\searrow 0$. For the metric term:

[TABLE]

and since $S$ is bounded, the last term vanishes uniformly on $K$ for $\varepsilon\searrow 0$. For the internal energy, since $\mathrm{D}G_{\varepsilon}=\mathrm{D}(\mathrm{id}+\varepsilon S)\circ G_{\boxplus}^{n}\cdot\mathrm{D}G_{\boxplus}^{n}$, and recalling (3.8),

[TABLE]

Since the piecewise constant function $\det\mathrm{D}G_{\boxplus}^{n}$ has a positive lower bound, the convergence as $\varepsilon\searrow 0$ is uniform on $K$. Finally, for the potential energy,

[TABLE]

Again, the convergence is uniform on $K$. Passing to the limit in the integral (4.8) yields

[TABLE]

The same inequality is true with $-S$ in place of $S$, hence this inequality is actually an equality. Since $\rho_{\boxplus}^{n}=(G_{\boxplus}^{n})_{\#}\overline{\rho}_{\mathscr{T}}$, a change of variables $x=S_{\boxplus}^{n}(\omega)$ produces (4.7). ∎

Corollary 4.5**.**

In addition to the hypotheses of Theorem 4.2, assume that

1. (1)

$P(\widetilde{\rho}_{\boxplus})\overset{*}{\rightharpoonup}p_{*}$* in $L^{1}([0,T]\times\Omega)$;* 2. (2)

each $G_{\boxplus}^{n}$ is injective; 3. (3)

as $\boxplus\to 0$, all simplices in the images of $\mathscr{T}$ under $G_{\boxplus}^{n}$ have non-degenerate interior angles and tend to zero in diameter, uniformly w.r.t. $n$.

Then $\rho_{*}$ satisfies the PDE

[TABLE]

in the sense of distributions.

Proof.

Let a smooth test function $\zeta\in C^{\infty}_{c}(\mathbb{R}^{d}\to\mathbb{R}^{d})$ be given. For each $\boxplus$ and each $n$, a $\zeta_{\boxplus}^{n}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ with $\zeta_{\boxplus}^{n}\circ G_{\boxplus}^{n}\in\mathcal{D}(\mathscr{T})$ can be constructed in such a way that

[TABLE]

uniformly on $\mathbb{R}^{d}$, and uniformly in $n$ as $\boxplus\to 0$. This follows from our hypotheses on the $\boxplus$-uniform regularity of the Lagrangian meshes: inside the image of $G_{\boxplus}^{n}$, one can simply choose $\zeta_{\boxplus}^{n}$ as the affine interpolation of the values of $\zeta$ at the points $G_{\boxplus}^{n}(\omega_{\ell})$. Outside, one can take an arbitrary approximation of $\zeta$ that is compatible with the piecewise-affine approximation on the boundary of $G_{\boxplus}^{n}$’s image; one may even choose $\zeta_{\boxplus}^{n}\equiv\zeta$ at sufficient distance to that boundary. The uniform convergences (4.10) then follow by standard finite element analysis.

Further, let $\eta\in C^{\infty}_{c}(0,T)$ be given. For each $t\in((n-1)\tau,n\tau]$, substitute $S(t,x):=\eta(t)\zeta_{\boxplus}^{n}(x)$ into (4.7). Integration of these equalities with respect to $t\in(0,T)$ yields

[TABLE]

We pass to the limit $\boxplus\to 0$ in these integrals. For the first, we use that $P(\widetilde{\rho})\overset{*}{\rightharpoonup}p_{*}$ by hypothesis, for the last, we use Theorem 4.2 above. Since any test function $S\in C^{\infty}_{c}((0,T)\times\Omega)$ can be approximated in $C^{1}$ by linear combinations of products $\eta(t)\zeta(x)$ as above, we thus obtain the weak formulation of

[TABLE]

In combination with the continuity equation (4.5), we arrive at (4.9). ∎

Remark 4.6**.**

In principle, our discretization can also be applied to the linear Fokker-Planck equation with $P(r)=r$ and $h(r)=r\log r$. In that case, one automatically has $P(\widetilde{\rho})\overset{*}{\rightharpoonup}p_{*}\equiv P(\rho_{*})$ thanks to Theorem 4.2. Corollary 4.5 above then provides an a posteriori criterion for convergence: if the Lagrangian mesh does not deform too wildly under the dynamics as the discretization is refined, then the discrete solutions converge to the genuine solution.

5. Consistency in 2D

In this section, we prove consistency of our discretization in the following sense. Under certain conditions on the spatial discretization $\mathscr{T}$, any smooth and positive solution $\rho$ to the initial value problem (1.1) projects to a discrete solution that satisfies the Euler-Lagrange equations up to a controlled error. We restrict ourselves to $d=2$ dimensions.

5.1. Smooth Lagrangian evolution

First, we derive an alternative form of the velocity field $\mathbf{v}$ from (1.3b) in terms of $G$.

Lemma 5.1**.**

For $\rho=G_{\#}\overline{\rho}$ with a smooth diffemorphism $G\colon K\to{\mathbb{R}^{d}}$, we have

[TABLE]

Consequently, the Lagrangian map $G$ — relative to the reference density $\overline{\rho}$ — for a smooth solution $\rho$ to (1.1) satisfies

[TABLE]

Proof.

On the one hand,

[TABLE]

and on the other hand, by definition of the push forward,

[TABLE]

Hence

[TABLE]

Observing that (1.2) implies that $rh^{\prime\prime}(r)=P^{\prime}(r)$, we conclude (5.2) directly from (1.3b). ∎

5.2. Discrete Euler-Lagrange equations in dimension $d=2$

In the planar case $d=2$, the Euler-Lagrange equation (3.10) above can be rewritten in a more convenient way.

In the following, fix some vertex $\omega_{\times}$ of the triangulation, which is indicent to precisely six triangles. For convenience, we assume that these are labelled $\Delta_{0}$ to $\Delta_{5}$ in counter-clockwise order. Similarly, the six neighboring vertices are labeled $\omega_{0}$ to $\omega_{5}$ in counter-clockwise order, so that $\Delta_{k}$ has vertices $\omega_{k}$ and $\omega_{k+1}$, where we set $\omega_{6}:=\omega_{0}$.

Using these conventions and recalling Lemma B.2, the expression for the vector $\nu$ in (3.9) simplifies to

[TABLE]

Summing the Euler-Lagrange equation (3.10) over $\Delta_{0}$ to $\Delta_{5}$, we obtain

[TABLE]

where the momentum term $\mathbf{p}_{\times}$ and the impulse $\mathbf{J}_{\times}$, respectively, are given by

[TABLE]

We shall now prove our main result on consistency. The setup is the following: a sequence of triangulations $\mathscr{T}_{\varepsilon}$ on $K$, parametrized by $\varepsilon>0$, and a sequence of time steps $\tau_{\varepsilon}=\mathcal{O}(\varepsilon)$ are given. We assume that there is an $\varepsilon$-independent region $K^{\prime}\subset K$ on which the $\mathscr{T}_{\varepsilon}$ are almost hexagonal in the following sense: each node $\omega_{\times}\in K^{\prime}$ of $\mathscr{T}_{\varepsilon}$ has precisely six neighbors — labelled $\omega_{0}$ to $\omega_{5}$ in counter-clockwise order — and there exists a rotation $R\in\mathrm{SO}(2)$ such that

[TABLE]

for $k=0,1,\ldots,5$.

Now, let $G\colon[0,T]\times K\to{\mathbb{R}^{d}}$ be a given smooth solution to the Lagrangian evolution equation (5.2), and fix a time $t\in(0,T)$. For all sufficiently small $\varepsilon>0$, we define maps $G_{\varepsilon},G_{\varepsilon}^{*}\in\mathcal{A}_{\mathscr{T}_{\varepsilon}}$ by linear interpolation of the values of $G(t;\cdot)$ and $G(t-\tau;\cdot)$, respectively, on $\mathscr{T}_{\varepsilon}$. That is, $G_{\varepsilon}(\omega_{\ell})=G(t;\omega_{\ell})$ and $G^{*}_{\varepsilon}(\omega_{\ell})=G(t-\tau;\omega_{\ell})$, at all nodes $\omega_{\ell}$ in $\mathscr{T}_{\varepsilon}$. Theorem 5.2 below states that the pair $G_{\varepsilon},G_{\varepsilon}^{*}$ is an approximate solution to the discrete Euler-Lagrange equations (5.3) at all nodes $\omega_{\times}$ of the respective triangulation $\mathscr{T}_{\varepsilon}$ that lie in $K^{\prime}$.

The hexagonality hypothesis on the $\mathscr{T}_{\varepsilon}$ is strong, but some very strong restriction of $\mathcal{A}_{\mathscr{T}_{\varepsilon}}$’s geometry is apparently necessary. See Remark 5.4 following the proof for further discussion.

Theorem 5.2**.**

Under the hypotheses and with the notations introduced above, the Euler-Lagrange equation (5.3) admits the following asymptotic expansion:

[TABLE]

as $\varepsilon\to 0$, uniformly at the nodes $\omega_{\times}\in K^{\prime}$ of the respective $\mathscr{T}_{\varepsilon}$.

Remark 5.3**.**

Up to an error $\mathcal{O}(\varepsilon^{3})$, the geometric pre-factor $\frac{\sqrt{3}}{2}\varepsilon^{2}$ equals to one third of the total area of the hexagon with vertices $\omega_{0}$ to $\omega_{5}$, and is thus equal to the integral of the piecewise affine hat function with peak at $\omega_{\times}$.

Proof of Theorem 5.2.

Throughout the proof, let $\varepsilon>0$ be fixed; we shall omit the $\varepsilon$-index for $\mathscr{T}_{\varepsilon}$ and $\tau_{\varepsilon}$. First, we fix a node $\omega_{\times}$ of $\mathscr{T}\cap K^{\prime}$. Thanks to the equivariance of both (5.2) and (5.3) under rigid motions of the domain, we may assume that $R$ in (5.7) is the identity, and that $\omega_{\times}=0$.

We collect some relations that are helpful for the calculations that follow. Trivially,

[TABLE]

Moreover, we have that

[TABLE]

On the other hand, by definition of $\mu_{\mathscr{T}}^{k}$ in (3.1), it follows that

[TABLE]

Combining (5.10) and (5.11) yields

[TABLE]

In accordance with the definition of $G_{\varepsilon}$ and $G_{\varepsilon}^{*}$ from $G$ detailed above, let $G_{\times}:=G(t,\omega_{\times})$ and $G^{*}_{\times}=G(t-\tau,\omega_{\times})$, and define $G_{k}$, $G_{k}^{*}$ for $k=0,\ldots,5$ in the analogous way. Further, we introduce $\mathrm{D}G_{\times}=\mathrm{D}G(t,\omega_{\times})$, $\mathrm{D}^{2}G_{\times}=\mathrm{D}^{2}G(t,\omega_{\times})$, $\partial_{t}G_{\times}=\partial_{t}G(t,\omega_{\times})$.

To perform an expansion in the momentum term, first observe that

[TABLE]

for each $k=0,1,\ldots,5$, and so, using that $\tau=\mathcal{O}(\varepsilon)$ by hypothesis,

[TABLE]

Using (5.12) and then (5.9) yields

[TABLE]

This is (5.8a).

For the impulse term, we start with a Taylor expansion to second order in space:

[TABLE]

We combine this with the observation that $(\omega_{k}|\omega_{k+1})^{-1}=\mathcal{O}(\varepsilon^{-1})$ to obtain:

[TABLE]

where

[TABLE]

Plugging this in leads to

[TABLE]

where we have use the auxiliary algebraic results from Lemma B.2, Lemma B.3, and Lemma B.4.

For the remaining part of the impulse term, a very rough approximation is sufficient:

[TABLE]

holds for any $g$ that is a convex combination of $G_{\times},G_{0},\ldots,G_{5}$, where the implicit constant is controlled in terms of the supremum of $\mathrm{D}^{2}V$ and $\mathrm{D}G$ on $K^{\prime}$. With that, we simply have, using again (5.12):

[TABLE]

Together, this yields (5.8b). ∎

Remark 5.4**.**

The hypotheses of Theorem (5.2) require that the $\mathscr{T}_{\varepsilon}$ are almost hexagonal on $K^{\prime}$. This seems like a technical hypothesis that simplifies calculations, but apparently, some strong symmetry property of the $\mathscr{T}_{\varepsilon}$ is necessary for the validity of the result.

To illustrate the failure of consistency — at least in the specific form considered here — assume that $V\equiv 0$ and $\overline{\rho}\equiv 1$, and consider a sequence of triangulations $\mathscr{T}_{\varepsilon}$ for which there is a node $\omega_{\times}$ such that (5.7) holds with the $\sigma_{k}$ being replaced by a different six-tuple of vectors $\sigma^{\prime}_{k}$. Repeating the steps of the proof above, it is easily seen that $\mathbf{p}_{\times}=a\varepsilon^{2}\,\partial_{t}G(t;\omega_{\times})+\mathcal{O}(\varepsilon^{3})$, with an $\varepsilon$-independent constant $a>0$ in place of $\sqrt{3}/2$, and that

[TABLE]

with

[TABLE]

If a result of the form (5.8b) — with $\sqrt{3}/2$ replaced by $a$ — was true, then this implies in particular that

[TABLE]

holds with some constant $a^{\prime}>0$ for arbitrary matrices $\mathrm{D}G_{\times}\in\mathbb{R}^{2\times 2}$ of positive determinant and tensors $\mathrm{D}^{2}G_{\times}\in\mathbb{R}^{2\times 2\times 2}$ that are symmetric in the second and third component. A specific example for which (5.13) is not true is given by

[TABLE]

in combination with $\mathrm{D}G_{\times}=\mathds{1}$, and a $\mathrm{D}^{2}G_{\times}$ that is zero except for two ones, at the positions $(1,2,2)$ and $(2,1,1)$. In Lemma B.5, we show that the left-hand side in (5.13) equals to $1\choose 1$; on the other hand, the right-hand side is clearly zero.

Note that this counter-example is significant, insofar as the skew (in fact, degenerate) hexagon described by the $\sigma_{k}^{\prime}$ in (5.14) corresponds to a popular method for triangulation of the plane.

6. Numerical simulations in $d=2$

6.1. Implementation

The Euler-Lagrange equations for the $d=2$-dimensional case have been derived in (5.3). We perfom a small modification in the potential term in order to simplify calculations with presumably minimal loss in accuracy:

[TABLE]

with the short-hand notation

[TABLE]

On the main diagonal, the Hessian amounts to

[TABLE]

Off the main diagonal, the entries of the Hessian are given by

[TABLE]

The scheme consists of an inner (Newton) and an outer (time stepping) iteration. We start from a given initial density $\rho_{0}$ and define the solution at the next time step inductively by applying Newton’s method in the inner iteration. To this end we initialise $G^{(0)}:=G^{n}$ with $G^{n}$, the solution at the $n$th time step, and define inductively

[TABLE]

where the update $\delta G^{(s+1)}$ is the solution to the linear system

[TABLE]

The effort of each inner iteration step is essentially determined by the effort to invert the sparse matrix $\mathbf{H}[G^{(s)}]$. As soon as the norm of $\delta G^{(s+1)}$ drops below a given stopping threshold, define $G^{n+1}:=G^{(s+1)}$ as approximate solution in the $n+1$st time step.

In all experiments the stopping criterion in the Newton iteration is set to $10^{-9}$.

6.2. Numerical experiments

In this section we present results of our numerical experiments for (1.1) with a cubic porous-medium nonlinearity $P(r)=r^{3}$ and different choices for the external potential $V$,

[TABLE]

Numerical experiment 1: unconfined evolution of Barenblatt profile

As a first example, we consider the “free” cubic porous medium equation, that is (6.1) with $V\equiv 0$. It is well-known (see, e.g., Vazquez [33]) that in the long-time limit $t\to\infty$, arbitrary solutions approach a self-similar one,

[TABLE]

where $\mathcal{B}_{3}$ is the associated Barenblatt profile

[TABLE]

where $C_{3}=(2\pi)^{-\frac{2}{3}}\approx 0.29$ is chosen to normalize $\mathcal{B}_{3}$’s mass to unity.

In this experiment, we are only interested in the quality of the numerical approximation for the self-similar solution (6.2). To reduce numerical effort, we impose a four-fold symmetry of the approximation: we use the quarter circle as computational domain $K$, and interprete the discrete function thereon as one of four symmetric pieces of the full discrete solution. To preserve reflection symmetry over time, homogeneous Neumann conditions are imposed on the artificial boundaries. This is implemented by reducing the degrees of freedom of the nodes along the $x$- and $y$-axes to tangential motion.

We initialize our simulation with a piecewise constant approximation of the profile of $\rho^{*}$ from (6.3) at time $t=0.01$. We choose a time step $\tau=0.001$ and the final time $T=2$. In Figure 1, we have collected snapshots of the approximated density at different instances of time. The Barenblatt profile of the solution is very well pertained over time.

Remark 6.1**.**

It takes less than 2 minutes to complete this simulation on standard laptop (Matlab code on a mid-2013 MacBook Air 11” with 1.7 GHz Intel Core i7 processor).

Figure 2 shows surface plots of the discrete solution at different times in comparison with the Barenblatt profile at the respective time. By construction of the scheme, the initial mass is exactly conserved in time as the discrete solution propagates. The left plot in Figure 3 shows the decay in the energy and gives quantitative information about the difference of the discrete solution to the analytical Barenblatt solution. The numerical solution shows good agreement with the analytical energy decay rate $c=2/3$.

We also compute the $l_{1}$-error of the discrete solution to the exact Barenblatt profile and observe that it remains within the order of the fineness of the triangulation. The mass of the discrete solution is perfectly conserved, as guaranteed by the construction of our method.

To estimate the convergence order of our method, we run several experiments with the above initial data on different meshes. We fix the ratio $\tau/h_{\rm max}^{2}=0.4$ and compute the $l_{1}$-error at time $T=0.2$ on triangulations with $h_{\rm max}=0.2,\,0.1,\,0.05,\,0.025.$ We expect the error to decay as a power of $h_{\rm max}$. The double logarithmic plot should reveal a line with its slope indicating the numerical convergence order. The right plot in Figure 3 shows the result, the estimated numerical convergence order which is obtained from a least-squares fitted line through the points is equal to $1.18$. This indicates first order convergence of the scheme with respect to the spatial discretisation parameter $h_{\rm max}$.

Numerical experiment 2: Asymptotic self-similarity

In our second example, we are still concerned with the free cubic porous medium equation, (6.1) with $V\equiv 0$. This time, we wish to give an indication that the discrete approximation of the self-similar solution from (6.2) from the previous experiment might inherit the global attractivity of its continuous counterpart. More specifically, we track the discrete evolution for the initial datum

[TABLE]

until time $T=0.1$ and observe that it appears to approach the self-similar solution from above. Snapshots of the simulation are collected in Figure 4.

Numerical experiment 3: two peaks merging into one under the influence of a confining potential

In this example we consider as initial condition two peaks, connected by a thin layer of mass, given by

[TABLE]

We choose a triangulation of the square $[-1.5,1.5]^{2}$ and initialise the discrete solution piecewise constant in each triangle, with a value corresponding to (6.5), evaluated in the centre of mass of each triangle. We solve the porous medium equation with a confining potential, i.e. (1.1) with $P(r)=r^{m}$ and $V(x,y)=5(x^{2}+y^{2})/2$. The time step is $\tau=0.001$ and the final time is $T=0.2.$

Figure 5 shows the evolution from the initial density. As time increases the peaks smoothly merge into each other. As the thin layer around the peaks is also subject to the potential the triangulated domain shrinks in time. Even if we do not know how to prevent theoretically the intersection of the images of the discrete Lagrangian maps, this seems not to be a problem in practice.

As time evolves, the discrete solution approaches the steady state Barenblatt profile given by

[TABLE]

where $C$ is chosen as the mass of the density. The plot in Figure 6 shows the exponential decay of the $l_{1}$-distance of the discrete solution to the steady state Barenblatt profile (6.6). We observe that the decay agrees very well with the analytically predicted decay $\exp(-5t)$ until $t=0.08$. For larger times, one would monitor triangle quality numerically, and re-mesh, locally coarsening the triangulation where necessary.

Numerical experiment 4: one peak splitting under the influence of a quartic potential

We consider as the initial condition

[TABLE]

We choose a triangulation of the unit circle and initialise the discrete solution piecewise constant in each triangle, with a value corresponding to (6.7), evaluated in the centre of mass of each triangle. We solve the porous medium equation with a quartic potential, i.e. (1.1) with $P(r)=r^{m}$ and $V(x)=5(x^{2}+(1-y^{2})^{2})/2$. The time step is $\tau=0.005$ and the final time is $T=0.02.$

Figure 7 shows the evolution of the initial density. As time increases the initial density is progressively split, until two new maxima emerge which are connected by a thin layer. For larger times, when certain triangles become excessively distorted, one would monitor triangle quality numerically, and re-mesh, locally refining the triangulation where necessary.

Appendix A Proof of the Lagrangian representation

Proof of Lemma 1.1.

We verify that the density function given by $(G_{t}^{-1})_{\#}\rho_{t}$ on $K\subset{\mathbb{R}^{d}}$ is constant with respect to time $t$; the identity (1.5) then follows since

[TABLE]

Firstly, from the definition of the inverse,

[TABLE]

for all $t$, differentiating with respect to time yields

[TABLE]

and so, using (1.4) and (1.3b),

[TABLE]

Now, let $\varphi$ be a smooth test function, and consider

[TABLE]

As $\varphi$ was arbitrary, $(G_{t}^{-1})_{\#}\rho_{t}$ is constant with respect to time. ∎

Appendix B Technical lemmas

Lemma B.1**.**

Given $g_{0},g_{1},\ldots,g_{d}\in\mathbb{R}^{d}$, then

[TABLE]

Proof.

Thanks to the symmetry of the integral with respect to the exchange of the components $\omega_{j}$, the left-hand side of (B.1) equals to

[TABLE]

We calculate the integrals, using Fubini’s theorem. First integral:

[TABLE]

Second integral:

[TABLE]

Third integral:

[TABLE]

Substitute this into (B.2):

[TABLE]

Collecting terms yields the right-hand side of (B.1). ∎

Lemma B.2**.**

For each $A\in\mathbb{R}^{2\times 2}$, we have $\mathbb{J}A\mathbb{J}^{T}=(\det A)\,A^{-T}$.

Proof.

This is verified by direct calculation:

[TABLE]

Lemma B.3**.**

With $\sigma_{k}\in\mathbb{R}^{2}$ defined as in (5.7), we have that

[TABLE]

Proof.

With the abbreviations $\phi_{x}=\frac{\pi}{3}x$ and $\psi=\frac{\pi}{3}$:

[TABLE]

Lemma B.4**.**

Let the scheme $B:=(b_{pqr})_{p,q,r\in\{1,2\}}\in\mathbb{R}^{2\times 2\times 2}$ of eight numbers $b_{pqr}\in\mathbb{R}$ be symmetric in the last two indices, $b_{pqr}=b_{prq}$. With $\sigma_{k}\in\mathbb{R}^{2}$ defined as in (5.7), we have that

[TABLE]

Proof.

In principle, this lemma can be verified by a direct calculation, by writing out the six terms in the sum explicitly and using trigonometric identities. Below, we give a slightly more conceptual proof, in which we use symmetry arguments to reduce the number of expressions significantly.

For the matrix involving $B$, we obtain

[TABLE]

while clearly

[TABLE]

The sum of the diagonal entries of the matrix product are easily calculated,

[TABLE]

with the trigonometric expressions

[TABLE]

To key step is to calculate the sum over $k=0,1,\ldots,5$ of the products of $T_{k}$ with the respective vector

[TABLE]

Several simplifications of this sum can be performed, thanks to the particular form of the $\gamma_{pqr,k}$ and elementary trigonometric identities. First, observe that $\sigma_{k+3}=-\sigma_{k}$, and hence that $\gamma_{pqr,k+3}=-\gamma_{pqr,k}$. Since further $\eta_{k+3}=-\eta_{k}$, it follows that

[TABLE]

Second, $\eta$ can be evaluated explicitly for $k=1,2,3$:

[TABLE]

Third, since $\sigma_{0,1}=-\sigma_{3,1}$ and $\sigma_{1,1}=-\sigma_{2,1}$, as well as $\sigma_{0,2}=\sigma_{3,2}$ and $\sigma_{1,2}=\sigma_{2,2}$, we obtain that

[TABLE]

By putting this together, we arrive at

[TABLE]

where $\mathfrak{e}_{pqr}=1$ if $p+q+r$ is even, and $\mathfrak{e}_{pqr}=0$ if $p+q+r$ is odd. By elementary computations,

[TABLE]

and so the final result is:

[TABLE]

which is (B.3). ∎

Lemma B.5**.**

With $\sigma_{k}^{\prime}\in\mathbb{R}^{2}$ defined as in (5.14), and with $B=(b_{pqr})_{p,q,r\in\{1,2\}}\in\mathbb{R}^{2\times 2\times 2}$ such that $b_{pqr}=0$ except for $b_{122}=b_{211}=1$, we have that

[TABLE]

Proof.

This is a slightly tedious, but straightforward calculation. First, by the choice of $B$,

[TABLE]

and so, by definition of the $\sigma_{k}^{\prime}$ in (5.14),

[TABLE]

For the inverse matrices S_{k}:=\big{(}\sigma_{k}^{\prime}\big{|}\sigma_{k+1}^{\prime}\big{)}^{-1}, we obtain

[TABLE]

For the traces T_{k}:=\operatorname{tr}\big{[}S_{k}\beta_{k}\big{]}, we thus obtain the values:

[TABLE]

In conclusion,

[TABLE]

which is (B.7). ∎

Appendix C Lack of convexity

Below, we discuss why the minimization problem (3.7) is not convex. More precisely, we show that $G\mapsto\mathbf{E}_{\boxplus}(G;\hat{G})$ is not convex as a function of $G$ on the affine ansatz space $\mathcal{A}_{\mathscr{T}}$. Since $\mathbf{E}_{\boxplus}(G;\hat{G})$ is a convex combination of the expressions \mathbb{H}_{m}\big{(}(A_{m}|b_{m});(\hat{A}_{m}|\hat{b}_{m})\big{)}, it clearly suffices to discuss the convexity of the latter.

We consider a curve $s\mapsto(A_{m}+s\alpha_{m}|b_{m}+s\beta_{m})$ and evaluate the second derivatives of the components of the functional at $s=0$. First,

[TABLE]

Second,

[TABLE]

If we assume that $\nabla^{2}V\geq\lambda\mathds{1}$, then we obtain for the sum of these two contributions that

[TABLE]

For the remaining term, however, we obtain — using the abbreviations $\widetilde{g}(s)=s\widetilde{h}^{\prime}(s)$ and $\widetilde{f}(s)=s\widetilde{g}^{\prime}(s)$ — that

[TABLE]

Now observe that $\widetilde{f}(s)=P^{\prime}(1/s)-sP(1/s)$ is a non-negative, and $\widetilde{g}(s)=-sP(1/s)$ is a non-positive function. Thus, from the two terms in the final sum, the first one is generally non-negative whereas the second one is of indefinite sign. Choosing

[TABLE]

the sum is obviously negative.