A Lagrangian Scheme for the Solution of Nonlinear Diffusion Equations Using Moving Simplex Meshes

A Lagrangian numerical scheme for solving nonlinear degenerate Fokker–Planck equations in space dimensions d≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\ge 2$$\end{document} 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=2$$\end{document}. 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.


Nonlinear Fokker-Planck Equations
We study a variational Lagrangian discretization of the following type of initial value problem: This problem is posed for the time-dependent probability density function ρ : R ≥0 × R d → R ≥0 , with a given initial density ρ 0 . We assume that the pressure P : R ≥0 → R ≥0 can be written in the form for some non-negative and convex h ∈ C 1 (R ≥0 ) ∩ C ∞ (R >0 ), and that V ∈ C 2 (R d ) is a non-negative potential without loss of generality. Problem (1.1) encompasses a large class of diffusion equations, such as-for power-type nonlinearities P(r ) = r m and vanishing potential V ≡ 0-the heat equation (m = 1), porous medium equations (m > 1) and fast diffusion equations (m < 1). By a slight abuse of notation, we refer to (1.1) with more general P and non-vanishing V as nonlinear Fokker-Planck equations. In this paper, we assume a degenerate diffusion, that is h(0) = h (0) = 0, and a confining potential, that is V is convex, not necessarily strict. For technical reasons, we further need to assume that lim s→∞ sh (s) = +∞. (1.3) Since our particular spatio-temporal discretization of the initial value problem (1.1) is based on the Lagrangian representation of its dynamics, and on its variational formulation, we briefly recall both of them now.

Lagrangian Formulation
Equation (1.1) can be written as a transport equation, with a velocity field v that depends on the solution ρ itself, Various further evolution equations can be written in the form (1.4a), such as non-local aggregation equations (see, e.g., Ambrosio et al. [1]); Keller-Segel type models (see, e.g., Blanchet et al. [5]); and also fourth order thin film equations (see, e.g., Otto [34]) or quantum equations (see, e.g., Gianazza et al. [21]). To simplify the presentation, we stick to equations of nonlinear Fokker-Planck type (1.1a).
The system (1.4) naturally induces a Lagrangian representation of the dynamics, which can be summarized as follows. Below, the reference density ρ is a probability density supported on some compact set K ⊂ R d , and we use the notation G # ρ for the push-forward of ρ under a map G : K → R d ; the definition is recalled in (2.1). Lemma 1.1 Assume that ρ : [0, T ] × R d → R ≥0 is a smooth positive solution of (1.1). Let G 0 : K → R d be a given map such that G 0 # ρ = ρ 0 . Further, let G : [0, T ] × R d → R d be the flow map associated to (1.4b), satisfying where ρ t := ρ(t, ·) and G t := G(t, ·) : R d → R d . Then, at any t ∈ [0, T ], In short, the solution G to (1.5) is a Lagrangian map for the solution ρ to (1.1). This fact is an immediate consequence of (1.4a); for convenience of the reader, we recall the proof in "Appendix A". Subsequently, (1.6) can be substituted for ρ in the expression (1.4b) for the velocity, which makes (1.5) an autonomous evolution equation for G: (1.7) A more explicit form of (1.7) is derived in (5.2).

Variational Structure
It is well-known (see Otto [35] or Ambrosio et al. [1]) that (1.1) is a gradient flow for the relative Renyi entropy functional (1.8) with respect to the L 2 -Wasserstein metric on the space P ac 2 (R d ) of probability densities on R d with finite second moment. It appears to be less well known (see Evans et al. [20], Carrillo and Moll [13], or Carrillo and Lisini [12]) that also (1.7) is a gradient flow, namely for the functional on the Hilbert space L 2 (K → R d ; ρ) of square integrable maps from K to R d . We shall discuss these gradient flow structures in more detail in Sect. 2 below.

Discretization and Approximation Results
Our discretization in space is based on the Lagrangian formulation. Instead of numerically integrating (1.1a) to obtain the density ρ directly, we approximate the associated Lagrangian maps G that satisfy (1.7): specifically, we assume that a simplicial decomposition T of K is given, and we restrict G to the finite dimensional subspace A T of continuous maps from K to R d that are piecewise linear with respect to T . A posteriori, we recover an approximation of ρ via (1.6). That ansatz for the Lagrangian maps corresponds to a simple geometric picture: the induced densities are piecewise constant on triangles whose vertices move in time.
For the discretization in time, we exploit the aforementioned variational structure of (1.7): namely, we adopt the celebrated minimizing movement scheme that is known to provide a robust approximation of gradient flows. In the context at hand, this scheme reads as follows: let a time step τ > 0 and an initial condition G 0 ∈ A T be given. (Here and below, symbolizes the space-time mesh generated by T on K and τ on R >0 .) Then the nth time iterate G n ∈ A T -that serves as our approximation of G(nτ ; ·)-is chosen inductively for n = 1, 2, . . . as the minimizer in the respective problem 1 2τ G − G n 2 L 2 (K →R d ;ρ) + E(G|ρ) −→ min, (1.10) where the minimization is carried out over the finite dimensional space A T . With the sequence (G n ) n=0,1,... of approximating Lagrangian maps at hand, we define piecewiseconstant-in-time interpolations for the derived density ρ and velocity v as usual via Our analytical results on the scheme can be summarized as follows.
• The sequence of fully discrete minimization problems (1.10) is well-posed: see Lemma 3.1. We thus obtain a sequence (G n ) n=0,1,... for each sufficiently fine discretization . • The G n are entropy-diminishing and are -uniformly Hölder continuous: see Lemma 4.1.
• Consequently, the induced densities ρ converge weakly to an absolutely continuous limit trajectory ρ, and the fluxes ρ v converge weakly to a limit of the form ρv: see Theorem 4.2. The identification of the limit velocity v, however, is only possible under strong additional hypotheses: see Corollary 4.5. • In d = 2 dimensions, we prove numerical consistency in the sense that, if G is a smooth solution to (1.7), then its restriction to the mesh satisfies the fully discrete Euler-Lagrange equations associated to (1.10), with a quantifiable error that vanishes in a suitable continuous limit: see Theorem 5.2.
• Our previously mentioned consistency results requires that the triangulation T of K is almost ideally hexagonal: see Eq. (5.7). We discuss why consistency cannot be expected if that condition is violated: see Remark 5.4.

Comparison with Results in the Literature
The approach presented in this paper is an alternative to the one developed by Carrillo et al. [13,15], where G is obtained by directly solving the PDE (1.7) numerically with finite differences or Galerkin approximation via finite element methods. In other words, while Carrillo et al. [13,15] follows the strategy minimize first then discretize, our present approach is to discretize first then minimize. In the former approach, the minimization (1.10) is performed on the spatially continuous level, yielding Euler-Lagrange equations that are then discretized in space; in the present approach, the space of Lagrangian maps is approximated by the finite dimensional subspace A T , and the minimization problem (1.10) on A T yields a nonlinear system of Euler-Lagrange equations that are directly solvable numerically. 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 nonnegativity and mass conservation, were proposed in the papers [4,8,10]. Particle methods based on suitable regularizations of the flux of the continuity Eq. (1.1) have been proposed in the papers [18,27,28,37]. A particle method based on the steepest descent of a regularized internal part of the energy E in (1.8) by substituting particles by non-overlapping blobs was proposed and analysed in Carrillo et al. [11,14]. Deterministic particle methods for diffusions have been recently explored, see [9] and the references therein. High-order relaxation schemes for nonlinear diffusion problems have been proposed in Cavalli et al. [16], while high-resolution schemes for nonlinear convection-diffusion problems are introduced in Kurganov et al. [26]. 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,7,23,40]) or solving for the optimal map in a JKO step (see [3,25]). Finally, note that gradient-flow-based Lagrangian methods in one dimension for higherorder, drift diffusion and Fokker-Planck equations have recently been proposed in the papers [19,[31][32][33].
There are two main arguments in favour of our taking this indirect approach of solving (1.7) 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 include monotonicity of the transformed entropy functional E and 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, which is reviewed in Sect. 2 below. However, despite significant effort in the recent past-see, e.g., the references [3,5,14,15,19,22,25,29,36,40]-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 [29]). This is mainly due to the non-existence of finite-dimensional submanifolds of P ac 2 (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 [6] 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 (0) = 0 of the diffusion in (1.1), solutions that are compactly supported initially remain compactly supported for 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 [38]). Our simulation results for ∂ t ρ = (ρ 3 )which possesses an analytically known, compactly supported, self-similar Barenblatt solution -indicate that our discretization is indeed able to track the edge of support quite accurately.
The expected convergence of our scheme, with implicit Euler stepping in time and piecewise linear approximation of the Lagrangian maps, is of first order in both space and time. This is confirmed in our experiments. For an improved approximation, particularly of the moving fronts, numerical schemes with a higher order of consistency would be desirable. In principle, such schemes could be constructed along the same lines, for example, by replacing the implicit Euler method by a Runge-Kutta method in time, and the piecewise constant ansatz space A T by finite elements with functions of higher global regularity in space. However, it is unclear if a similar degree of structure preservation can be achieved for these schemes, and their analysis would be very different from the one presented here.

Structure of the Paper
This work is organized as follows. In Sect. 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 Sect. 6 gives several numerical tests showing the performance of this scheme.

Notations from Probability Theory
P(X ) is the space of probability measures on a given base set X . We say that a sequence (μ n ) of measures in P(X ) converges narrowly to a limit μ in that space if for all bounded and continuous functions f With a slight abuse of notation -identifying absolutely continuous measures with their densities-we denote the space of probability densities on R d of finite second moment by Clearly, the reference density ρ, which is supported on the compact set K ⊂ R d , belongs to P ac 2 (R d ). If G : K → R d is a diffeomorphism onto its image (which is again compact), then the push-forward of ρ's measure produces again a density G # ρ ∈ P ac 2 (R d ), given by

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 [39]. One of the many equivalent ways to define the L 2 -Wasserstein distance between ρ 0 , ρ 1 ∈ P ac 2 (R d ) is as follows: (2. 2) 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 [39, Section 2.1]. A trivial but essential observation is that if ρ ∈ P ac 2 (R d ) is a reference density with support K ⊂ R d , and ρ 0 = (G 0 ) # ρ with a measurable G 0 : K → R d , then (2.2) can be re-written as follows: and the essentially unique minimizer G * in (2.3) is related to the optimal map T * in (2.2) via G * = T * • G 0 . W 2 is a metric on P ac 2 (R d ); convergence in W 2 is equivalent to weak-convergence in L 1 (R d ) and convergence of the second moment. Since P and hence also h are of super-linear growth at infinity, each sublevel set E is weak-closed and thus complete with respect to W 2 .
As already mentioned above, solutions ρ to (1.1) constitute a gradient flow for the functional E from (1.8) in the metric space (P ac 2 (R d ); W 2 ). In fact, assuming that the potential V is λ-convex (i.e., ∇ 2 V ≥ λ1), the flow is even λ-contractive as a semi-group, thanks to the λ-uniform displacement convexity of E (see McCann [30], or Daneri and Savaré [17]), which is a strengthened form of λ-uniform convexity along geodesics. The λ-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 E and its sub-differential, uniform exponential estimates on the convergence (if λ > 0) or divergence (if λ ≤ 0) of trajectories, quantified exponential rates for the approach to equilibrium (if λ > 0) and the like.
An important further consequence is that the unique flow can be obtained as the limit for τ 0 of the time-discrete minimizing movement scheme (see Ambrosio et al. [1] and Jordan, Kinderlehrer and Otto [24]): This time discretization is well-adapted to approximate λ-contractive gradient flows. All of the properties of mentioned above are already reflected on the level of these time-discrete solutions. Evans et al. [20] or Jordan et al. [25]). However, the variational structure behind this gradient flow is much weaker than above: most notably, E is only poly-convex, but not λ-uniformly convex. Therefore, the abstract machinery for λcontractive gradient flows in Ambrosio et al. [1] does not apply here. Clearly, by equivalence of (1.1) and (1.7) at least for sufficiently smooth solutions, certain properties of the primary gradient flow are necessarily inherited by this secondary flow, but for instance λ-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:

Gradient Flow in L 2
(2.5) In fact, there is an equivalence between (2.5) and (2.4): simply substitute (G n−1 τ ) # ρ for ρ n−1 τ and G # ρ for ρ in (2.4); notice that any ρ ∈ P ac 2 (R d ) can be written as G # ρ with a suitable (highly non-unique) choice of G ∈ L 2 (K → R d ; ρ). This equivalence was already exploited in Carrillo et al. [13,15]. Thanks to the equality (2.3), the minimization with respect to ρ = G # ρ can be relaxed to a minimization with respect to G.

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 ρ ∈ P ac 2 (R d ), whose support K ⊂ R d is a compact, convex polytope.

Discretization in Space
Our spatial discretization is performed using a finite subspace of linear maps for the Lagrangian maps G. More specifically: let T be some (finite) simplicial decomposition of K with nodes ω 1 to ω L and n-simplices 1 to M . In the case d = 2, which is of primary interest here, T is a triangulation, with triangles m . The reference density ρ is approximated by a density ρ T ∈ P ac 2 (R d ) that is piecewise constant on the simplices of T , with respective values The finite dimensional ansatz space A T is now defined as the set of maps G : K → R d that are globally continous, affine on each of the simplices m ∈ T , and orientation preserving. That is, on each m ⊂ T , the map G ∈ A T can be written as follows: with a suitable matrix A m ∈ R d×d of positive determinant and a vector b m ∈ R d . For the calculations that follow, we shall use a more geometric way to describe the maps G ∈ A T , namely by the positions G = G(ω ) of the images of each node ω . Denote by the same simplicial combinatorics (including orientation) as the ω in T . Clearly, any G ∈ A T is uniquely characterized by the L-tuple G of its values, and moreover, any More explicitly, fix a m ∈ T , with nodes labelled ω m,0 to ω 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 Then the affine map (3.2) equals to q m • r −1 m ; this is shown schematically in Fig. 1 for the case d = 2. In particular, we obtain that For later reference, we give a more explicit representation for the transformed entropy E for G ∈ A T , and for the L 2 -distance between two maps G,Ĝ ∈ A T . Substitution of the special form (3.2) into (1.9) produces For the L 2 -difference of G and G * , we have Using Lemma B.1, we obtain on each simplex m :

Discretization in Time
Let a time step τ > 0 be given; in the following, we symbolize the spatio-temporal discretization by , and we write → 0 for the joint limit of τ → 0 and vanishing mesh size in T . The discretization in time is performed in accordance with (2.5): we modify E τ from (2.5) by restriction to the ansatz space A T . This leads to the minimization problem For a fixed discretization , the fully discrete scheme is well-posed in the sense that for a given initial map G 0 ∈ A T , an associated sequence (G n ) n≥0 can be determined by successive solution of the minimization problems (3.7). One only needs to verify: 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 E '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 d Ldimensional vector space V T of all continuous maps G : K → R d that are piecewise affine with respect to the fixed simplicial decomposition T : given , define the distance between these maps as the maximal R d -distance G −G of corresponding points G ∈ G, G ∈ G . Clearly, this metric makes V T a complete space.
It is easily seen that the subset A T -which is singled out by requiring orientation preservation of the G's-is an open subset of V T . It is further obvious that the map G → E (G; G * ) is continuous with respect to the metric. The claim of the lemma thus follows if we can show that the sub-level is a non-empty compact subset of V T . Clearly, G * ∈ 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 ∈ S c has a distance larger than R to G * . From non-negativity of E, and from the representations (3.5) and (3.6), it follows that It is now easy to compute a suitable value for the radius R. S c is a closed subset of V T . It suffices to show that the limit G ∈ V T of any sequence By definition of our metric on V T , global continuity and piecewise linearity of the G (k) trivially pass to the limit G. We still need to verify that G is orientation-preserving. Fix a simplex m and consider the corresponding and since h(s) → +∞ as s ↓ 0, it follows that det A (k) m > 0 is bounded away from zero, uniformly in k. But then also det A m > 0, i.e., the mth linear map piece of the limit G preserves orientation.

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 n−1 ∈ A T , we calculate the variations of E (G; G * ) with respect to the degrees of freedom of G ∈ A T . Since that function is a weighted sum over the triangles m ∈ T , it suffices to perform the calculations for one fixed triangle m , with respective nodes ω m,0 to ω 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 ω m,0 , it will suffice to perform the calculations at one fixed image point G m,0 .
• mass term: is the uniquely determined vector in 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: Now let ω be a fixed vertex of T . Summing over all simplices m that have ω as a vertex, and choosing vertex labels in accordance with above, i.e., such that ω m,0 = ω in m , produces the following Euler-Lagrange equation: (3.10)

Approximation of the Initial Condition
For the approximation ρ 0 = (G 0 ) # ρ T of the initial datum ρ 0 = G 0 # ρ, we require: • ρ 0 converges to ρ 0 narrowly; • E(ρ 0 ) is -uniformly bounded, i.e., (3.11) In our numerical experiments, we always choose ρ := ρ 0 , in which case G 0 : K → R d can be taken as the identity on K , and we choose accordingly G 0 as the identity as well. Hence ρ 0 = ρ T , which converges to ρ 0 = ρ even strongly in L 1 (K ). Moreover, since h is convex, it easily follows from Jensen's inequality that and therefore, In more general situations, in which G 0 is not the identity, a sequence of approximations G 0 of G 0 is needed. Pointwise convergence G 0 → G 0 is more than sufficient to guarantee narrow convergence of ρ 0 to ρ 0 , but the uniform bound (3.11) might require a well-adapted approximation, especially for non-smooth G 0 's.

Limit Trajectory
In this section, we assume that a sequence of vanishing discretizations → 0 is given, and we study the respective limit of the fully discrete solutions (G n ) n≥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, instead of the Lagrangian maps G n themselves. Note that v n is only well-defined on the support of ρ n -that is, on the image of G n -and can be assigned arbitrary values outside. Let us introduce the piecewise constant in time interpolations ρ : Note that ρ(t, ·) ∈ P ac 2 (R d ) and v (t, ·) ∈ L 2 (R d → R d ; ρ (t, ·)) at each t ≥ 0.

Energy Estimates
We start by proving the classical energy estimates on minimizing movements for our fully discrete scheme.
Monotonicity of E in time is obvious. To prove (4.2), choose n ≤ n such that s ∈ ((n − 1)τ, nτ ] and t ∈ ((n − 1)τ, nτ ]. Notice that τ (n − n) ≤ t − s + τ . If n = n, the claim (4.2) is obviously true; let n < n in the following. Combining the triangle inequality for the metric W 2 , estimate (4.1) above and Hölder's inequality for sums, we arrive at Finally, changing variables using x = G n (ω) in (4.4) yields and summing these inequalities from n = 1 to n = N τ yields (4.3).

Compactness of the Trajectories and Weak Formulation
Our main result on the weak limit of ρ is the following.
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 ≤ t ≤ T -it follows that W 2 ( ρ (t), ρ 0 ) is bounded, uniformly in t ∈ [0, T ] and in . Since further ρ 0 converges narrowly to ρ 0 by our hypotheses on the initial approximation, we conclude that all densities ρ (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 (2E) 1 2 |t − s| 1 2 in the limit → 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 ρ towards a Hölder continuous limit curve ρ * .
It remains to pass to the limit with the velocity v . Towards that end, we define a probability measure γ ∈ P(Z T ) on the set Z T := [0, T ] × R d × R d as follows: for every bounded and continuous function ϕ Thanks to the result from the first part of the proof, M converges narrowly to a limit M * , which has density ρ * (t,x) T . On the other hand, the estimate (4.3) implies that We are thus in the position to apply Theorem 5.4.4 in Ambrosio et al. [1], which yields the narrow convergence of γ towards a limit γ * . Clearly, the (t, x)-marginal of γ * is M * . Accordingly, we introduce the disintegration γ (t,x) of γ * with respect to M * , which is welldefined M * -a.e.. Below, it will turn out that γ is the sought-for weak limit of v . The convergence v ρ * v * ρ * and the inheritance of the uniform L 2 -bound (4.3) to the limit 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 v * is to evaluate the limit as → 0 of for any given test function φ ∈ C ∞ c ((0, T ) × R d ) in two different ways. First, we change variables t → t + τ in the second integral, which gives For the second evaluation, we write The error term e [φ] above is controlled via Taylor expansion of φ and by using (4.3),

Equality of the limits for both evaluations of J [φ]
for arbitrary test functions φ shows the continuity Eq. (4.5).
Unfortunately, the convergence provided by Theorem 4.2 is generally not sufficient to conclude that ρ * is a weak solution to (1.1), since we are not able to identify v * as v[ρ * ] from (1.4b). The problem is two-fold: first, weak-convergence of ρ is insufficient to pass to the limit inside the nonlinear function P. Second, even if we would know that, for instance, P( ρ ) * P(ρ * ), we would still need a -independent a priori control on the regularity (e.g., maximal diameter of triangles) of the meshes generated by the G 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 n . For definiteness, we introduce the space D(T ) := : K → R d ; is globally continuous, and is piecewise affine w.r.t. m .

Lemma 4.4 Assume S : R d → R d is such that S • G n ∈ D(T ). Then:
Proof For all sufficiently small ε > 0, let G ε = (id + S) • G n . By definition of G n as a minimizer, we have that E (G ε ; G n−1 ) ≥ E (G n ; G n−1 ). This implies that We discuss limits of the three terms under the integral for ε 0. For the metric term: and since S is bounded, the last term vanishes uniformly on K for ε 0. For the internal energy, since DG ε = D(id + εS) • G n · DG n , and recalling (3.8), Since the piecewise constant function det DG n has a positive lower bound, the convergence as ε 0 is uniform on K . Finally, for the potential energy, Again, the convergence is uniform on K . Passing to the limit in the integral (4.8) yields The same inequality is true with −S in place of S, hence this inequality is actually an equality. Since ρ n = (G n ) # ρ T , a change of variables x = S n (ω) produces (4.7). Then ρ * satisfies the PDE in the sense of distributions.
Proof Let a smooth test function ζ ∈ C ∞ c (R d → R d ) be given. For each and each n, a ζ n : R d → R d with ζ n • G n ∈ D(T ) can be constructed in such a way that ζ n → ζ, ∇ · ζ n → ∇ · ζ (4.10) uniformly on R d , and uniformly in n as → 0. This follows from our hypotheses on theuniform regularity of the Lagrangian meshes: inside the image of G n , one can simply choose ζ n as the affine interpolation of the values of ζ at the points G n (ω ). Outside, one can take an arbitrary approximation of ζ that is compatible with the piecewise-affine approximation on the boundary of G n 's image; one may even choose ζ n ≡ ζ at sufficient distance to that boundary. The uniform convergences (4.10) then follow by standard finite element analysis.
Further, let η ∈ C ∞ c (0, T ) be given. For each t ∈ ((n − 1)τ, nτ ], substitute S(t, x) := η(t)ζ n (x) into (4.7). Integration of these equalities with respect to t ∈ (0, T ) yieldŝ We pass to the limit → 0 in these integrals. For the first, we use that P( ρ) * p * by hypothesis, for the last, we use Theorem 4.2 above. Since any test function S ∈ C ∞ c ((0, T ) × ) can be approximated in C 1 by linear combinations of products η(t)ζ (x) as above, we thus obtain the weak formulation of In combination with the continuity Eq. (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( ρ) * p * ≡ P(ρ * ) 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.

Consistency in 2D
In this section, we prove consistency of our discretization in the following sense. Under certain conditions on the spatial discretization T , any smooth and positive solution ρ 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.

Smooth Lagrangian Evolution
First, we derive an alternative form of the velocity field v from (1.4b) in terms of G.
Consequently, the Lagrangian map G-relative to the reference density ρ -for a smooth solution ρ to (1.1) satisfies Proof On the one hand, and on the other hand, by definition of the push forward, Observing that (1.2) implies that rh (r ) = P (r ), we conclude (5.2) directly from (1.4b).

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 ω × of the triangulation, which is incident to precisely six triangles. For convenience, we assume that these are labelled 0 to 5 in counter-clockwise order. Similarly, the six neighboring vertices are labeled ω 0 to ω 5 in counter-clockwise order, so that k has vertices ω k and ω k+1 , where we set ω 6 := ω 0 .
Using these conventions and recalling Lemma B.2, the expression for the vector ν in (3.9) simplifies to Summing the Euler-Lagrange equation (3.10) over 0 to 5 , we obtain where the momentum term p × and the impulse J × , respectively, are given by We shall now prove our main result on consistency. The setup is the following: a sequence of triangulations T ε on K , parametrized by ε > 0, and a sequence of time steps τ ε = O(ε) are given. We assume that there is an ε-independent region K ⊂ K on which the T ε are almost hexagonal in the following sense: each node ω × ∈ K of T ε has precisely six neighborslabelled ω 0 to ω 5 in counter-clockwise order-and there exists a rotation R ∈ SO(2) such that for k = 0, 1, . . . , 5. Now, let G : [0, T ] × K → R d be a given smooth solution to the Lagrangian evolution Eq. (5.2), and fix a time t ∈ (0, T ). For all sufficiently small ε > 0, we define maps G ε , G * ε ∈ A T ε by linear interpolation of the values of G(t; ·) and G(t − τ ; ·), respectively, on T ε . That is, G ε (ω ) = G(t; ω ) and G * ε (ω ) = G(t − τ ; ω ), at all nodes ω in T ε . Theorem 5.2 below states that the pair G ε , G * ε is an approximate solution to the discrete Euler-Lagrange equations (5.3) at all nodes ω × of the respective triangulation T ε that lie in K .
The hexagonality hypothesis on the T ε is strong, but some very strong restriction of A T ε '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:
as ε → 0, uniformly at the nodes ω × ∈ K of the respective T ε .

Remark 5.3
Up to an error O(ε 3 ), the geometric pre-factor √ 3 2 ε 2 equals to one third of the total area of the hexagon with vertices ω 0 to ω 5 , and is thus equal to the integral of the piecewise affine hat function with peak at ω × .
Proof of Theorem 5.2 Throughout the proof, let ε > 0 be fixed; we shall omit the ε-index for T ε and τ ε . First, we fix a node ω × of T ∩ K . 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 ω × = 0.
We collect some relations that are helpful for the calculations that follow. Trivially, Moreover, we have that On the other hand, by definition of μ k T in (3.1), it follows that Combining (5.10) and (5.11) yields In accordance with the definition of G ε and G * ε from G detailed above, let G × := G(t, ω × ) and G * × = G(t − τ, ω × ), and define G k , G * k for k = 0, . . . , 5 in the analogous way. Further, we introduce DG × = DG(t, ω × ), To perform an expansion in the momentum term, first observe that for each k = 0, 1, . . . , 5, and so, using that τ = O(ε) by hypothesis, Using (5.12) and then (5.9) yields This is (5.8a).
For the impulse term, we start with a Taylor expansion to second order in space: We combine this with the observation that (ω k |ω k+1 ) −1 = O(ε −1 ) to obtain: Plugging this in leads to where we have use the auxiliary algebraic results from Lemmas B.2, B.3 and B.4.
For the remaining part of the impulse term, a very rough approximation is sufficient: holds for any g that is a convex combination of G × , G 0 , . . . , G 5 , where the implicit constant is controlled in terms of the supremum of D 2 V and DG on K . With that, we simply have, using again (5.12): Together, this yields (5.8b).

Remark 5.4
The hypotheses of Theorem (5.2) require that the T ε are almost hexagonal on K . This seems like a technical hypothesis that simplifies calculations, but apparently, some strong symmetry property of the T ε is necessary for the validity of the result.
To illustrate the failure of consistency-at least in the specific form considered hereassume that V ≡ 0 and ρ ≡ 1, and consider a sequence of triangulations T ε for which there is a node ω × such that (5.7) holds with the σ k being replaced by a different six-tuple of vectors σ k . Repeating the steps of the proof above, it is easily seen that p × = aε 2 ∂ t G(t; ω × )+O(ε 3 ), with an ε-independent constant a > 0 in place of √ 3/2, and that If a result of the form (5.8b)-with √ 3/2 replaced by a-was true, then this implies in particular that holds with some constant a > 0 for arbitrary matrices DG × ∈ R 2×2 of positive determinant and tensors D 2 G × ∈ R 2×2×2 that are symmetric in the second and third component. A specific example for which (5.13) is not true is given by 14) in combination with DG × = 1, and a D 2 G × 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 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 σ k in (5.14) corresponds to a popular method for triangulation of the plane. 6 Numerical Simulations in d = 2

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: with the short-hand notation On the main diagonal, the Hessian amounts to Off the main diagonal, the entries of the Hessian are given by ).
The scheme consists of an inner (Newton) and an outer (time stepping) iteration. We start from a given initial density ρ 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 nth time step, and define inductively where the update δG (s+1) is the solution to the linear system The effort of each inner iteration step is essentially determined by the effort to invert the sparse matrix H[G (s) ]. As soon as the norm of δG (s+1) drops below a given stopping threshold, define G n+1 := G (s+1) as approximate solution in the n + 1st time step.
In all experiments the stopping criterion in the Newton iteration is set to 10 −9 .

Numerical Experiments
In this section we present results of our numerical experiments for (1.1) with a cubic porousmedium nonlinearity P(r ) = r 3 and different choices for the external potential V ,

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 ≡ 0. It is well-known (see, e.g., Vazquez [38]) that in the long-time limit t → ∞, arbitrary solutions approach a self-similar one, where B 3 is the associated Barenblatt profile where C 3 = (2π) − 2 3 ≈ 0.29 is chosen to normalize 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 ρ * from (6.3) at time t = 0.01. We choose a time step τ = 0.001 and the final time T = 2. In Fig. 2, 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.   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 τ/ h 2 max = 0.4 and compute the l 1error at time T = 0.2 on triangulations with h max = 0.2, 0.1, 0.05, 0.025. We expect the error to decay as a power of h max . The double logarithmic plot should reveal a line with its slope indicating the numerical convergence order. The right plot in Fig. 4 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 max .

Numerical experiment 2: Asymptotic self-similarity
In our second example, we are still concerned with the free cubic porous medium Eq. (6.1) with V ≡ 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 ρ 0 (x, y) = 3000(x 2 + y 2 ) exp[−5(|x| + |y|)] + 0.1 ( 6 . 4 ) 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 Fig. 5.

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 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 τ = 0.001 and the final time is T = 0.2. Figure 6 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 where C is chosen as the mass of the density. The plot in Fig. 7 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.  Figure 8 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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.

Lemma B.5
With σ k ∈ R 2 defined as in (5.14), and with B = (b pqr ) p,q,r ∈{1,2} ∈ R 2×2×2 such that b pqr = 0 except for b 122 = b 211 = 1, we have that Proof This is a slightly tedious, but straightforward calculation. First, by the choice of B, and so, by definition of the σ k in (5.14), For the traces T k := tr S k β k , we thus obtain the values: In conclusion, 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 → E (G;Ĝ) is not convex as a function of G on the affine ansatz space A T . Since E (G;Ĝ) is a convex combination of the expressions H m (A m |b m ); (Â m |b m ) , it clearly suffices to discuss the convexity of the latter. We consider a curve s → (A m + sα m |b m + sβ m ) and evaluate the second derivatives of the components of the functional at s = 0. First, I :=