A Lagrangian scheme for the solution of nonlinear diffusion equations using moving simplex meshes
José A. Carrillo, Bertram Düring, Daniel Matthes, David S. McCormick
AA LAGRANGIAN SCHEME FOR THE SOLUTION OF NONLINEARDIFFUSION EQUATIONS USING MOVING SIMPLEX MESHES
JOS´E A. CARRILLO, BERTRAM D ¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK
Abstract.
A Lagragian numerical scheme for solving nonlinear degenerate Fokker-Planckequations in space dimensions d ≥ d = 2. Avariety of numerical experiments for the porous medium equation indicates that the schemeis well-adapted to track the growth of the solution’s support. Introduction
We study a Lagrangian discretization of the following type of initial value problem for anonlinear Fokker–Planck equation: ∂ t ρ = ∆ P ( ρ ) + ∇ · ( ρ ∇ V ) on R > × R d , (1.1a) ρ ( · ,
0) = ρ on R d . (1.1b)This problem is posed for the time-dependent probability density function ρ : R ≥ × R d → R ≥ ,with initial condition ρ ∈ P ac2 ( R d ). We assume that the pressure P : R ≥ → R ≥ can be writtenin the form P ( r ) = rh (cid:48) ( r ) − h ( r ) for all r ≥ , (1.2)for some non-negative and convex h ∈ C ( R ≥ ) ∩ C ∞ ( R > ) and that V ∈ C ( 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 − , m > , V = 0) and the fastdiffusion equation ( P ( r ) = r m / ( m − , m < , V = 0), and extends to related problems withgiven external potentials V . To motivate our discretization, we first briefly recall the Lagrangianform of the dynamics: rewriting (1.1) as a transport equation, we obtain ∂ t ρ + ∇ · (cid:0) ρ v [ ρ ] (cid:1) = 0 , (1.3a)with a velocity field v that depends on the solution ρ itself, v [ ρ ] = −∇ (cid:0) h (cid:48) ( ρ ) + V (cid:1) . (1.3b) Date : September 11, 2018.DM was supported by the DFG Collaborative Research Center TRR 109, “Discretization in Geometry andDynamics”. BD and DSMcC were supported by the Leverhulme Trust research project grant “Novel discretisationsfor higher-order nonlinear PDE” (RPG-2015-69). a r X i v : . [ m a t h . NA ] F e b JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK
More general systems can be written in this form. For instance, interaction potentials leading toaggregation equations, relativistic heat equations, p -Laplacian equations and Keller-Segel typemodels can also be included; see Carrillo and Moll [12] and the references therein for a goodaccount of models enjoying this structural form. Here, we reduce to models with nonlineardegenerate diffusion, i.e. h (0) = h (cid:48) (0) = 0, and confining potentials in order to explore our newdiscretization.The system (1.3) naturally induces a Lagrangian representation of the dynamics, which canbe summarized as follows. Below, we use the notation G ρ for the push-forward of ρ under amap G : R d → R d ; the definition is recalled in (2.1). Lemma 1.1.
Assume that ρ : [0 , T ] × R d → R ≥ is a smooth positive solution of (1.1) . Let ρ : R d → R ≥ be a given reference density, and let G : R d → R d be a given map such that G ρ = ρ . Further, let G : [0 , T ] × R d → R d be the flow map associated to (1.3b) , satisfying ∂ t G t = v [ ρ t ] ◦ G t , G (0 , · ) = G , (1.4) where ρ t := ρ ( t, · ) and G t := G ( t, · ) : R d → R d . Then ρ t = ( G t ) ρ (1.5) at any t ∈ [0 , T ] . In short, the solution G to (1.4) is a Lagrangian map for the solution ρ to (1.1). This factis an immediate consequence of (1.3a); for convenience of the reader, we recall the proof inAppendix A. Notice that (1.5) can be substituted for ρ in the expression (1.3b) for the velocity,which makes (1.4) an autonomous evolution equation for G : ∂ t G t = −∇ (cid:20) h (cid:48) (cid:18) ρ det D G t (cid:19)(cid:21) ◦ G t − ∇ V ◦ G t . (1.6)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´e [1]) that (1.1) is a gradient flow for the relative Renyi entropyfunctional E ( ρ ) = ˆ R d (cid:2) h ( ρ ( x )) + V ( x ) ρ ( x ) (cid:3) d x, (1.7)with respect to the L -Wasserstein metric on the space P ac2 ( R d ) of probability densities on R d .It appears to be less well known (see Evans et al. [18], Carrillo and Moll [12], or Carrillo andLisini [11]) that also (1.6) is a gradient flow, namely for the functional E ( G | ρ ) := E ( G ρ ) = ˆ K (cid:20)(cid:101) h (cid:18) det D Gρ (cid:19) + V ◦ G (cid:21) ρ d ω, (cid:101) h ( s ) := s h ( s − ) , (1.8)on the Hilbert space L ( K → R d ; ρ ) of maps from K to R d , where ρ is a reference measuresupported on K ⊂ R d . We shall discuss these gradient flow structures in more detail in Section 2below.The particular spatio-temporal discretization of the initial value problem (1.1) that we studyin this paper is based on these facts. Instead of numerically integrating (1.1a) to obtain thedensity ρ directly, we approximate the Lagrangian map G — using a finite subspace of linearmaps in space — which then allows us to recover ρ a posteriori via (1.5). Our approximation for G is constructed by solving a succession of minimization problems that are naturally derived fromthe gradient flow structure behind (1.4). Via variational methods, this provides compactnessestimates on the discrete solutions. AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 3
The use of a finite subspace of linear maps for the Lagrangian maps has a geometric inter-pretation: the induced densities are piecewise constant on triangles whose vertices move in time.In the sequel, we further assume that lim s →∞ sh (cid:48)(cid:48) ( s ) = + ∞ in order to prevent the collapse ofthe images of our Lagrangian maps G t , as proven below. However, we do not yet know how toprevent (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 isobtained by directly solving the PDE (1.6) numerically with finite differences or Galerkin ap-proximation via finite element methods. In fact, the main difference between these two strategiescan be summarized as follows: while Carrillo et al. [14, 12] follows the strategy minimize firstthen discretize , our present approach is to discretize first then minimize . In other words, inCarrillo et al. [14, 12] one minimizes first, obtaining as Euler-Lagrange equations the implicitEuler discretization of (1.6) in [14] approximated by the explicit Euler method in [12], and thendiscretized in space. In the present approach, we discretize first approximating the space ofLagrangian maps by a suitable finite subspace of linear maps, and then we minimize obtaining anonlinear system of equations to find the approximated Lagrangian map within that set of linearmaps.Let us mention that other numerical methods have been developed to conserve particularproperties of solutions of the gradient flow (1.1). Finite volume methods preserving the decayof energy at the semi-discrete level, along with other important properties like non-negativityand mass conservation, were proposed in the papers [4, 9]. Particle methods based on suitableregularisations 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 theenergy E in (1.7) by substituting particles by non-overlapping blobs was proposed and analysed inCarrillo et al. [10, 13]. Moreover, the numerical approximation of the JKO variational scheme hasalready 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 thatgradient-flow-based Lagrangian methods in one dimension for higher-order, drift diffusion andFokker-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 : thescheme that we present builds on the non-obvious “secondary” gradient flow representationof (1.1) in terms of Lagrangian maps. The benefits are monotonicity of the transformed en-tropy functional E and an L -control on the metric velocity for our fully discrete solutions, thateventually lead to weak compactness of the trajectories in the continuous limit. We remarkthat 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 Sec-tion 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 likemetric contractivity of the flow under the discretization, except in the rather special situation ofone space dimension (see Matthes and Osberger [25]). This is mainly due to the non-existenceof finite-dimensional submanifolds of P ac2 ( R d ) that are complete with respect to generalizedgeodesics.The second motivation is that Lagrangian schemes are a natural choice for numerical fronttracking , see, e.g., Budd [7] for first results on the numerical approximation of self-similar solu-tions to the porous medium equation. We recall that due to the assumed degeneracy P (cid:48) (0) = 0of the diffusion in (1.1), solutions that are compactly supported initially remain compactly sup-ported at all times. A numerically accurate calculation of the moving edge of support is chal-lenging, since the solution can have a very complex behavior near that edge, like the waitingtime phenomenon (see Vazquez [33]). Our simulation results for ∂ t ρ = ∆( ρ ) — that possess an JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK analytically known, compactly supported, self-similar Barenblatt solution — indicated that ourdiscretization 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 ingradient flows pertaining our work. Section 3 is devoted to the introduction of the linear set ofLagragian maps and the derivation of the numerical scheme. Section 4 shows the compactnessof the approximated sequences of discretizations and we give conditions leading to the eventualconvergence of the scheme towards (1.1). Section 5 deals with the consistency of the schemein two dimensions while Section 6 gives several numerical tests showing the performance of thisscheme. 2.
Gradient flow structures
Notations from probability theory. P ( X ) is the space of probability measures on agiven base set X . We say that a sequence ( µ n ) of measures in P ( X ) converges narrowly to alimit µ in that space if ˆ X f ( x ) d µ n ( x ) → ˆ X f ( x ) d µ ( x )for all bounded and continuous functions f ∈ C b ( X ). The push-forward T µ of a measure µ ∈ P ( X ) under a measurable map T : X → Y is the uniquely determined measure ν ∈ P ( Y )such that, for all g ∈ C b ( Y ), ˆ X g ◦ T ( x ) d µ ( x ) = ˆ Y g ( y ) d ν ( y ) . 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 P ac2 ( R d ) = (cid:26) ρ ∈ L ( R d ) ; ρ ≥ , ˆ R d ρ ( x ) d x = 1 , ˆ R d (cid:107) x (cid:107) ρ ( x ) d x < ∞ (cid:27) . Clearly, the reference density ρ , which is supported on the compact set K ⊂ R d , belongs to P ac2 ( R d ). If G : K → R d is a diffeomorphism onto its image (which is again compact), then thepush-forward of ρ ’s measure produces again a density G ρ ∈ P ac2 ( R d ), given by G ρ = ρ det D G ◦ G − . (2.1)2.2. Gradient flow in the Wasserstein metric.
Below, some basic facts about the Wasser-stein 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] andVillani [34].One of the many equivalent ways to define the L -Wasserstein distance between ρ , ρ ∈P ac2 ( R d ) is as follows:W ( ρ , ρ ) := inf (cid:26) ˆ R d (cid:107) T ( x ) − x (cid:107) ρ ( x ) d x ; T : R d → R d measurable , T ρ = ρ (cid:27) . (2.2)The infimum above is in fact a minimum, and the — essentially unique — optimal map T ∗ ischaracterized by Brenier’s criterion; see, e.g., Villani [34, Section 2.1]. A trivial but essentialobservation is that if ρ ∈ P ac2 ( R d ) is a reference density with support K ⊂ R d , and ρ = ( G ) ρ with a measurable G : K → R d , then (2.2) can be re-written as follows:W ( ρ , ρ ) = inf (cid:26) ˆ K (cid:107) G ( ω ) − G ( ω ) (cid:107) ρ ( ω ) d ω ; G : K → R d measurable , G ρ = ρ (cid:27) , (2.3) AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 5 and the essentially unique minimizer G ∗ in (2.3) is related to the optimal map T ∗ in (2.2) via G ∗ = T ∗ ◦ G .W is a metric on P ac2 ( R d ); convergence in W is equivalent to weak- (cid:63) convergence in L ( R d )and convergence of the second moment. Since P and hence also h are of super-linear growth atinfinity, each sublevel set E is weak- (cid:63) closed and thus complete with respect to W .As already mentioned above, solutions ρ to (1.1) constitute a gradient flow for the functional E from (1.7) in the metric space ( P ac2 ( R d ); W ). In fact, the flow is even λ -contractive as asemi-group, thanks to the λ -uniform displacement convexity of E (see McCann [26], or Daneriand Savar´e [15]), 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) and thelike.An important further consequence is that the unique flow can be obtained as the limit for τ (cid:38) minimizing movement scheme (see Ambrosio et al. [1] and Jordan,Kinderlehrer and Otto [21]): ρ nτ := argmin ρ ∈P ac2 ( R d ) E τ ( ρ ; ρ n − τ ) , E τ ( ρ, ˆ ρ ) := 12 τ W ( ρ, ˆ ρ ) + E ( ρ ) . (2.4)This time discretization is well-adapted to approximate λ -contractive gradient flows. All of theproperties of mentioned above are already reflected on the level of these time-discrete solutions.2.3. Gradient flow in L . Equation (1.6) is the gradient flow of E on the space L ( K → R d ; ρ )of square integrable (with respect to ρ ) maps G : K → R d (see Evans et al. [18] or Jordanet al. [22]). However, the variational structure behind this gradient flow is much weaker thanabove: most notably, E is only poly-convex, but not λ -uniformly convex . Therefore, the abstractmachinery for λ -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 propertiesof the primary gradient flow are necessarily inherited by this secondary flow, but for instance λ -contractivity of the flow in the L -norm seems to fail.Nevertheless, it can be proven (see Ambrosio, Lisini and Savar´e [2]) that the gradient flow isglobally well-defined, and it can again be approximated by the minimizing movement scheme: G nτ := argmin G ∈ L ( K → R d ; ρ ) E τ (cid:0) G ; G n − τ (cid:1) , E τ ( G ; ˆ G ) = 12 τ ˆ K (cid:107) G − ˆ G (cid:107) d ρ + E ( G | ρ ) . (2.5)In fact, there is an equivalence between (2.5) and (2.4): simply substitute ( G n − τ ) ρ for ρ n − τ and G ρ for ρ in (2.4); notice that any ρ ∈ P ac2 ( R d ) can be written as G ρ with a suitable (highlynon-unique) choice of G ∈ L ( K → R d ; ρ ). This equivalence was already exploited in Carrilloet al. [14, 12]. Thanks to the equality (2.3), the minimization with respect to ρ = G ρ can berelaxed to a minimization with respect to G . Consequently, if ( G τ ) ρ = ρ τ , then ( G nτ ) ρ = ρ nτ at all discrete times n = 1 , , . . . . However, while the functional E τ ( · ; ρ n − τ ) in (2.4) is ( λ + τ − )-uniformly convex in ρ along geodesics in W , the functional E τ ( · ; G n − τ ) in (2.5) has apparentlyno useful convexity properties in G on L ( K → R d ; ρ ).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 ρ ∈ P ac2 ( R d ), whose support K ⊂ R d is convex and compact. JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK
Discretization in space.
Our spatial discretization is performed using a finite subspaceof linear maps for the Lagrangian maps G . More specifically: let T be some (finite) simplicialdecomposition of K with nodes ω to ω L and n -simplices ∆ to ∆ M . In the case d = 2, whichis of primary interest here, T is a triangulation, with triangles ∆ m . The reference density ρ isapproximated by a density ρ T ∈ P ac2 ( R d ) that is piecewise constant on the simplices of T , withrespective values ρ m T := µ m T | ∆ m | for the simplex masses µ m T := ˆ ∆ m ρ ( ω ) d ω. (3.1)The finite dimensional ansatz space A T is now defined as the set of maps G : K → R d that areglobally 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: G ( ω ) = A m ω + b m for all ω ∈ ∆ m , (3.2)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 (cid:96) = G ( ω (cid:96) ) of the images of each node ω (cid:96) . Denote by ( R d ) L T ⊂ R L · d the space of L -tuples (cid:126)G = ( G (cid:96) ) L(cid:96) =1 of points G (cid:96) ∈ R d with the same simplicial combinatorics(including orientation) as the ω (cid:96) in T . Clearly, any G ∈ A T is uniquely characterized by the L -tuple (cid:126)G of its values, and moreover, any (cid:126)G ∈ ( R d ) L T defines a G ∈ A T .More explicitly, fix a ∆ m ∈ T , with nodes labelled ω m, to ω m,d in some orientation preservingorder, and respective image points G m, to G m,d . With the standard d -simplex given by (cid:52) d := ξ = ( ξ , . . . , ξ d ) ∈ R d ≥ ; d (cid:88) j =1 ξ j ≤ , introduce the linear interpolation maps r m : (cid:52) d → K and q m : (cid:52) d → R d by r m ( ξ ) = ω m, + d (cid:88) j =1 ( ω m,j − ω m, ) ξ j ,q m ( ξ ) = G m, + d (cid:88) j =1 ( G m,j − G m, ) ξ j . Then the affine map (3.2) equals to q m ◦ r − m . In particular, we obtain thatdet A m = det D q m det D r m = det Q m T [ G ]2 | ∆ m | where Q m T [ G ] := (cid:0) G m, − G m, (cid:12)(cid:12) · · · (cid:12)(cid:12) G m,d − G m, (cid:1) . (3.3)For later reference, we give a more explicit representation for the transformed entropy E for G ∈ A T , and for the L -distance between two maps G, ˆ G ∈ A T . Substitution of the specialform (3.2) into (1.8) produces E ( G | ρ T ) = (cid:88) ∆ m ∈ T µ m T (cid:2) H m T ( G ) + V m T ( G ) (cid:3) (3.4)with the internal energy (recall the definition of (cid:101) h from (1.8)) H m T ( G ) := (cid:101) h (cid:18) det A m ρ m T (cid:19) = (cid:101) h (cid:18) det Q m T [ G ]2 µ m T (cid:19) AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 7 and the potential energy V m T ( G ) = ∆ m V ( A m ω + b m ) d ω = (cid:52) V (cid:0) r m ( ω ) (cid:1) d ω. For the L -difference of G and G ∗ , we have (cid:107) G − G ∗ (cid:107) L ( K ; ρ T ) = ˆ K (cid:107) G − G ∗ (cid:107) ρ T d ω = (cid:88) ∆ m ∈ T µ m T L m T ( G, G ∗ ) . (3.5)Using Lemma B.1, we obtain on each simplex ∆ m : L m T ( G, G ∗ ) := ∆ m (cid:107) G ( ω ) − G ∗ ( ω ) (cid:107) d ω = (cid:52) (cid:107) r m ( ω ) − r ∗ m ( ω ) (cid:107) d ω = 2( d + 1)( d + 2) (cid:88) ≤ i ≤ j ≤ d ( G m,i − G ∗ m,i ) · ( G m,j − G ∗ m,j ) . (3.6)3.2. Discretization in time.
Let a time step τ > (cid:1) , and we write (cid:1) → τ → T .The discretization in time is performed in accordance with (2.5): we modify E τ from (2.5) byrestriction to the ansatz space A T . This leads to the minimization problem G n (cid:1) := argmin G ∈A T E (cid:1) (cid:0) G ; G n − (cid:1) (cid:1) where E (cid:1) ( G ; G ∗ ) = 12 τ (cid:107) G − G ∗ (cid:107) L ( K ; ρ T ) + E ( G | ρ T ) . (3.7)For a fixed discretization (cid:1) , the fully discrete scheme is well-posed in the sense that for a giveninitial map G (cid:1) ∈ A T , an associated sequence ( G n (cid:1) ) n ≥ can be determined by successive solutionof the minimization problems (3.7). One only needs to verify: Lemma 3.1.
For each given G ∗ ∈ A T , there exists at least one global minimizer G ∈ A T of E (cid:1) ( · ; G ∗ ) . Remark 3.2. We do not claim uniqueness of the minimizers. Unfortunately, the minimizationproblem (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 (cid:1) ’s(non-)convexity is provided in Appendix C.Proof of Lemma 3.1. We only sketch the main arguments. For definiteness, let us choose (justfor this proof) one of the infinitely many equivalent norm-induced metrics on the dL -dimensionalvector space V T of all continuous maps G : K → R d that are piecewise affine with respect tothe fixed simplicial decomposition T : given G, G (cid:48) ∈ V T with their respective point locations (cid:126)G, (cid:126)G (cid:48) ∈ R dL , i.e., (cid:126)G = ( G (cid:96) ) L(cid:96) =1 for G (cid:96) = G ( ω (cid:96) ), define the distance between these maps as themaximal R d -distance (cid:107) G (cid:96) − G (cid:48) (cid:96) (cid:107) of corresponding points G (cid:96) ∈ (cid:126)G , G (cid:48) (cid:96) ∈ (cid:126)G (cid:48) . Clearly, this metricmakes V T a complete space.It is easily seen that the subset A T — which is singled out by requiring orientation preservationof the G ’s — is an open subset of V T . It is further obvious that the map G (cid:55)→ E (cid:1) ( G ; G ∗ ) iscontinuous with respect to the metric. The claim of the lemma thus follows if we can show thatthe sub-level S c := { G ∈ A T ; E (cid:1) ( G ; G ∗ ) ≤ c } with c := E ( G ∗ | ρ T )is a non-empty compact subset of V T . Clearly, G ∗ ∈ S c , so it suffices to verify compactness. JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK S c is bounded. We are going to show that there is a radius
R > G ∈ S c hasa distance larger than R to G ∗ . From non-negativity of E , and from the representations (3.5)and (3.6), it follows that c ≥ τ (cid:107) G − G ∗ (cid:107) L ( K ; ρ T ) ≥ µ T τ (cid:88) ∆ m ∈ T L m T ( G, G ∗ )= µ T ( d + 1)( d + 2) τ (cid:88) ≤ i ≤ j ≤ d ( G m,i − G ∗ m,i ) · ( G m,j − G ∗ m,j ) ≥ µ T d + 1)( d + 2) τ L (cid:88) (cid:96) =1 (cid:107) G (cid:96) − G ∗ (cid:96) (cid:107) , where µ T = min ∆ m µ m T . 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 ( G ( k ) ) ∞ k =1 of maps G ( k ) ∈ S c belongs to A T . By definition of our metric on V T , global continuity andpiecewise linearity of the G ( k ) trivially pass to the limit G . We still need to verify that G isorientation-preserving. Fix a simplex ∆ m and consider the corresponding matrices A ( k ) m and A m from (3.2). Since the G ( k ) converge to G in the metric, also A ( k ) m → A m entry-wise. Now, bynon-negativity of (cid:101) h , we have for all k that c ≥ E ( G ( k ) | ρ T ) ≥ µ m T (cid:101) h (cid:32) det A ( k ) m ρ m T (cid:33) , and since (cid:101) h ( s ) → + ∞ as s ↓
0, it follows that det A ( k ) m > k . But then also det A m >
0, i.e., the m th linear map piece of the limit G preserves orientation. (cid:3) Fully discrete equations.
We shall now derive the Euler-Lagrange equations associatedto the minimization problem (3.7), i.e., for each given G ∗ := G n − (cid:1) ∈ A T , we calculate thevariations of E (cid:1) ( G ; G ∗ ) with respect to the degrees of freedom of G ∈ A T . Since that functionis a weighted sum over the triangles ∆ m ∈ T , it suffices to perform the calculations for one fixedtriangle ∆ m , with respective nodes ω m, to ω m,d , in positive orientation. The associated imagepoints are G m, to G m,d . Since we may choose any vertex to be labelled ω m, , it will suffice toperform the calculations at one fixed image point G m, . • mass term: ∂∂G m, L m T ( G, G ∗ ) = 2( d + 1)( d + 2) ∂∂G m, (cid:88) ≤ i ≤ j ≤ d ( G m,i − G ∗ m,i ) · ( G m,j − G ∗ m,j )= 2( d + 1)( d + 2) G m, − G ∗ m, ) + d (cid:88) j =1 ( G m,j − G ∗ m,j ) • internal energy: observing that — recall (1.2) — (cid:101) h (cid:48) ( s ) = dd s (cid:2) sh ( s − ) (cid:3) = h ( s − ) − s − h (cid:48) ( s − ) = − P ( s − ) , (3.8)we obtain ∂∂G m, H m T ( G ) = ∂∂G m, (cid:101) h (cid:18) det Q m T [ G ]2 µ m T (cid:19) = 12 µ m T P (cid:18) µ m T det Q m T [ G ] (cid:19) ν m T [ G ] , AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 9 where ν m T [ G ] := − ∂∂G m, det Q m T [ G ] = (det Q m T [ G ]) ( Q m T [ G ]) − T d (cid:88) j =1 e j (3.9)is the uniquely determined vector in R d that is orthogonal to the ( d − G m, to G m,d (pointing away from G m, ) and whose length equals the ( d − • potential energy: ∂∂G m, V m T ( G ) = ∂∂G m, (cid:52) V (cid:0) r m ( ξ ) (cid:1) d ξ = (cid:52) ∇ V (cid:0) r m ( ξ ) (cid:1) (1 − ξ − · · · − ξ d ) d ξ. Now let ω (cid:96) be a fixed vertex of T . Summing over all simplices ∆ m that have ω (cid:96) as a vertex, andchoosing vertex labels in accordance with above, i.e., such that ω m, = ω (cid:96) in ∆ m , produces thefollowing Euler-Lagrange equation:0 = (cid:88) ω (cid:96) ∈ ∆ m µ m T (cid:20) d + 1)( d + 2) τ (cid:18) G m, − G ∗ m, ) + d (cid:88) j =1 ( G m,j − G ∗ m,j ) (cid:19) + 12 µ m T P (cid:18) µ m T det Q m T [ G ] (cid:19) ν m T [ G ] + (cid:52) ∇ V (cid:0) r m ( ξ ) (cid:1) (1 − ξ − · · · − ξ d ) d ξ (cid:21) . (3.10)3.4. Approximation of the initial condition.
For the approximation ρ (cid:1) = ( G (cid:1) ) ρ T of theinitial datum ρ = G ρ , we require: • ρ (cid:1) converges to ρ narrowly; • E ( ρ (cid:1) ) is (cid:1) -uniformly bounded, i.e., E := sup E ( ρ (cid:1) ) < ∞ . (3.11)In our numerical experiments, we always choose ρ := ρ , in which case G : K → R d can be takenas the identity on K , and we choose accordingly G (cid:1) as the identity as well. Hence ρ (cid:1) = ρ T ,which converges to ρ = ρ even strongly in L ( K ). Moreover, since h is convex, it easily followsfrom Jensen’s inequality that ˆ ∆ m h (cid:0) ρ ( x ) (cid:1) d x ≥ | ∆ m | h ( ρ m T ) , and therefore, E ( ρ (cid:1) ) ≤ E ( ρ ) . In more general situations, in which G is not the identity, a sequence of approximations G (cid:1) of G is needed. Pointwise convergence G (cid:1) → G is more than sufficient to guarantee narrow con-vergence of ρ (cid:1) to ρ , but the uniform bound (3.11) might require a well-adapted approximation,especially for non-smooth G ’s. 4. Limit trajectory
In this section, we assume that a sequence of vanishing discretizations (cid:1) → G n (cid:1) ) n ≥ that are produced by theinductive minimization procedure (3.7). For the analysis of that limit trajectory, it is morenatural to work with the induced densities and velocities, ρ n (cid:1) := ( G n (cid:1) ) ρ, v n (cid:1) := id − G n − (cid:1) ◦ ( G n (cid:1) ) − τ , instead of the Lagrangian maps G n (cid:1) themselves. Note that v n (cid:1) is only well-defined on the supportof ρ n (cid:1) — that is, on the image of G n (cid:1) — and can be assigned arbitrary values outside. Let usintroduce the piecewise constant in time interpolations (cid:101) ρ (cid:1) : [0 , T ] × R d → R ≥ , and (cid:101) v (cid:1) : [0 , T ] × R d → R d as usual, (cid:101) ρ (cid:1) ( t ) = ρ n (cid:1) , (cid:101) v (cid:1) ( t ) = v n (cid:1) with n such that t ∈ (( n − τ, nτ ] . Note that (cid:101) ρ ( t, · ) ∈ P ac2 ( R d ) and (cid:101) v (cid:1) ( t, · ) ∈ L ( R d → R d ; (cid:101) ρ (cid:1) ( t, · )) at each t ≥ Energy estimates.
We start by proving the classical energy estimates on minimizing move-ments for our fully discrete scheme.
Lemma 4.1.
For each discretization (cid:1) and for any indices n > n ≥ , one has the a prioriestimate E ( ρ n (cid:1) ) + τ n (cid:88) n = n +1 (cid:32) W ( ρ n (cid:1) , ρ n − (cid:1) ) τ (cid:33) ≤ E ( ρ n ) . (4.1) Consequently: (1) E is monotonically decreasing, i.e., E ( (cid:101) ρ (cid:1) ( t )) ≤ E ( (cid:101) ρ (cid:1) ( s )) for all t ≥ s ≥ ; (2) (cid:101) ρ (cid:1) is H¨older- -continuous in W , up to an error τ , W (cid:0)(cid:101) ρ (cid:1) ( t ) , (cid:101) ρ (cid:1) ( s ) (cid:1) ≤ (cid:113) E ( ρ (cid:1) ) (cid:0) | t − s | + τ (cid:1) for all t ≥ s ≥ . (4.2)(3) (cid:101) v (cid:1) is square integrable with respect to (cid:101) ρ (cid:1) , ˆ T ˆ R d (cid:107) (cid:101) v (cid:1) (cid:107) (cid:101) ρ (cid:1) d x d t ≤ E ( ρ (cid:1) ) . (4.3) Proof.
By the definition of G n (cid:1) as a minimizer, we know that E (cid:1) ( G n (cid:1) ; G n − (cid:1) ) ≤ E (cid:1) ( G ; G n − (cid:1) ) forany G ∈ A T , and in particular for the choice G := G n − (cid:1) , which yields:12 τ ˆ K (cid:107) G n (cid:1) − G n − (cid:1) (cid:107) ρ T d ω + E ( G n (cid:1) | ρ T ) ≤ E ( G n − (cid:1) | ρ T ) . (4.4)Summing these inequalies for n = n + 1 , . . . , n , recalling that E ( ρ n (cid:1) ) = E ( G n (cid:1) | ρ T ) by (1.8) andthat W ( ρ n (cid:1) , ρ n − (cid:1) ) ≤ ´ K | G n (cid:1) − G n − (cid:1) | ρ d ω by (2.3), produces (4.1).Monotonicity of E in time is obvious.To prove (4.2), choose n ≤ n such that s ∈ (( n − τ, nτ ] and t ∈ (( n − τ, 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 , estimate (4.1) above and H¨older’s inequalityfor sums, we arrive atW (cid:0)(cid:101) ρ (cid:1) ( t ) , (cid:101) ρ (cid:1) ( s ) (cid:1) = W ( ρ n (cid:1) , ρ n (cid:1) ) ≤ n (cid:88) n = n +1 W ( ρ n (cid:1) , ρ n − (cid:1) ) ≤ n (cid:88) n = n +1 τ n (cid:88) n = n +1 W ( ρ n (cid:1) , ρ n − (cid:1) ) τ = (cid:2) τ ( n − n ) (cid:3) τ n (cid:88) n = n +1 (cid:32) W ( ρ n (cid:1) , ρ n − (cid:1) ) τ (cid:33) ≤ [ t − s + τ ] (cid:2) (cid:0) E ( ρ n (cid:1) ) − E ( ρ n (cid:1) ) (cid:1)(cid:3) ≤ (cid:2) | t − s | + τ (cid:3) E ( ρ (cid:1) ) . AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 11
Finally, changing variables using x = G n (cid:1) ( ω ) in (4.4) yields τ ˆ R d (cid:107) v n (cid:1) (cid:107) ρ n (cid:1) d x + E ( G n (cid:1) ) ≤ E ( G n − (cid:1) ) , and summing these inequalities from n = 1 to n = N τ yields (4.3). (cid:3) Compactness of the trajectories and weak formulation.
Our main result on the weaklimit of (cid:101) ρ (cid:1) is the following. Theorem 4.2.
Along a suitable sequence (cid:1) → , the curves (cid:101) ρ (cid:1) : R ≥ → P ac ( R d ) convergencepointwise in time, i.e., (cid:101) ρ (cid:1) ( t ) → ρ ∗ ( t ) narrowly for each t > , towards a H¨older- -continuouslimit trajectory ρ ∗ : R ≥ → P ac ( R d ) .Moreover, the discrete velocities (cid:101) v (cid:1) possess a limit v ∗ ∈ L ( R ≥ × R d ; ρ ∗ ) such that (cid:101) v (cid:1) (cid:101) ρ (cid:1) ∗ (cid:42) v ∗ ρ ∗ in L ( R ≥ × R d ) , and the continuity equation ∂ t ρ ∗ + ∇ · ( ρ ∗ v ∗ ) = 0 (4.5) holds in the sense of distributions. Remark 4.3.
The H¨older continuity of ρ ∗ implies that ρ ∗ satisfies the initial condition (1.1b) in the sense that ρ ∗ ( t ) → ρ narrowly as t ↓ .Proof of Theorem 4.2. We closely follow an argument that is part of the general convergenceproof 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 diagonalargument implies convergence at arbitrary times.First observe that by estimate (4.2) — applied with 0 = s ≤ t ≤ T — it follows thatW ( (cid:101) ρ (cid:1) ( t ) , ρ (cid:1) ) is bounded, uniformly in t ∈ [0 , T ] and in (cid:1) . Since further ρ (cid:1) converges narrowlyto ρ by our hypotheses on the initial approximation, we conclude that all densities (cid:101) ρ (cid:1) ( t ) belongto a sequentially compact subset for the narrow convergence. The second observation is that theterm on the right hand side of (4.2) simplifies to (2 E ) | t − s | in the limit (cid:1) →
0. A straight-forward application of the “refined version” of the Ascoli-Arzel`a theorem (Proposition 3.3.1 inAmbrosio et al. [1]) yields the first part of the claim, namely the pointwise narrow convergenceof (cid:101) ρ (cid:1) towards a H¨older continuous limit curve ρ ∗ .It remains to pass to the limit with the velocity (cid:101) v (cid:1) . Towards that end, we define a probabilitymeasure (cid:101) γ (cid:1) ∈ P ( Z T ) on the set Z T := [0 , T ] × R d × R d as follows: ˆ Z T ϕ ( t, x, v ) d (cid:101) γ (cid:1) ( t, x, v ) = ˆ T ˆ R d ϕ (cid:0) t, x, (cid:101) v (cid:1) ( t, x ) (cid:1) (cid:101) ρ (cid:1) ( t, x ) d x d tT , for every bounded and continuous function ϕ ∈ C b ( Z T ). For brevity, let (cid:102) M (cid:1) ∈ P ([0 , T ] × R d ) bethe ( t, x )-marginals of (cid:101) γ (cid:1) , that have respective Lebesgue densities ρ (cid:1) ( t,x ) T on [0 , T ] × R d . Thanksto the result from the first part of the proof, (cid:102) M (cid:1) converges narrowly to a limit M ∗ , which hasdensity ρ ∗ ( t,x ) T . On the other hand, the estimate (4.3) implies that ˆ Z T | v | d (cid:101) γ (cid:1) ( t, x, v ) = ˆ [0 ,T ] × R d | (cid:101) v (cid:1) ( t, x ) | d (cid:102) M (cid:1) ( t, x ) ≤ E . We are thus in the position to apply Theorem 5.4.4 in Ambrosio et al. [1], which yields the narrowconvergence of (cid:101) γ (cid:1) towards a limit γ ∗ . Clearly, the ( t, x )-marginal of γ ∗ is M ∗ . Accordingly, weintroduce the disintegration γ ( t,x ) of γ ∗ with respect to M ∗ , which is well-defined M ∗ -a.e.. Below,it will turn out that γ ∗ ’s v -barycenter, v ∗ ( t, x ) := ˆ R d v d γ ( t,x ) ( v ) , (4.6)2 JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK
0. A straight-forward application of the “refined version” of the Ascoli-Arzel`a theorem (Proposition 3.3.1 inAmbrosio et al. [1]) yields the first part of the claim, namely the pointwise narrow convergenceof (cid:101) ρ (cid:1) towards a H¨older continuous limit curve ρ ∗ .It remains to pass to the limit with the velocity (cid:101) v (cid:1) . Towards that end, we define a probabilitymeasure (cid:101) γ (cid:1) ∈ P ( Z T ) on the set Z T := [0 , T ] × R d × R d as follows: ˆ Z T ϕ ( t, x, v ) d (cid:101) γ (cid:1) ( t, x, v ) = ˆ T ˆ R d ϕ (cid:0) t, x, (cid:101) v (cid:1) ( t, x ) (cid:1) (cid:101) ρ (cid:1) ( t, x ) d x d tT , for every bounded and continuous function ϕ ∈ C b ( Z T ). For brevity, let (cid:102) M (cid:1) ∈ P ([0 , T ] × R d ) bethe ( t, x )-marginals of (cid:101) γ (cid:1) , that have respective Lebesgue densities ρ (cid:1) ( t,x ) T on [0 , T ] × R d . Thanksto the result from the first part of the proof, (cid:102) M (cid:1) converges narrowly to a limit M ∗ , which hasdensity ρ ∗ ( t,x ) T . On the other hand, the estimate (4.3) implies that ˆ Z T | v | d (cid:101) γ (cid:1) ( t, x, v ) = ˆ [0 ,T ] × R d | (cid:101) v (cid:1) ( t, x ) | d (cid:102) M (cid:1) ( t, x ) ≤ E . We are thus in the position to apply Theorem 5.4.4 in Ambrosio et al. [1], which yields the narrowconvergence of (cid:101) γ (cid:1) towards a limit γ ∗ . Clearly, the ( t, x )-marginal of γ ∗ is M ∗ . Accordingly, weintroduce the disintegration γ ( t,x ) of γ ∗ with respect to M ∗ , which is well-defined M ∗ -a.e.. Below,it will turn out that γ ∗ ’s v -barycenter, v ∗ ( t, x ) := ˆ R d v d γ ( t,x ) ( v ) , (4.6)2 JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK is the sought-for weak limit of (cid:101) v (cid:1) . The convergence (cid:101) v (cid:1) (cid:101) ρ (cid:1) ∗ (cid:42) v ∗ ρ ∗ and the inheritance of theuniform L -bound (4.3) to the limit v ∗ are further direct consequences of Theorem 5.4.4 inAmbrosio et al. [1].The key step to establish the continuity equation for the just-defined v ∗ is to evaluate thelimit as (cid:1) → J (cid:1) [ φ ] := 1 τ (cid:34) ˆ T ˆ R d φ ( t, x ) (cid:101) ρ (cid:1) ( t, x ) d x d t − ˆ T ˆ R d φ ( t, x ) (cid:101) ρ (cid:1) ( t − τ, x ) d x d t (cid:35) for any given test function φ ∈ C ∞ c ((0 , T ) × R d ) in two different ways. First, we change variables t (cid:55)→ t + τ in the second integral, which gives J (cid:1) [ φ ] = ˆ T ˆ R d φ ( t, x ) − φ ( t + τ, x ) τ (cid:101) ρ (cid:1) ( t, x ) d x d t (cid:1) → −→ − ˆ T ˆ R d ∂ t φ ( t, x ) ρ ∗ ( t, x ) d x d t. For the second evaluation, we write ρ n − (cid:1) = (cid:0) G n − (cid:1) ◦ ( G n (cid:1) ) − (cid:1) ρ n (cid:1) = (cid:0) id − τ v n (cid:1) (cid:1) ρ n (cid:1) , and substitute accordingly x (cid:55)→ x − τ (cid:101) v (cid:1) ( t, x ) in the second integral, leading to J (cid:1) [ φ ] = ˆ T ˆ R d φ ( t, x ) − φ (cid:0) t, x − τ (cid:101) v (cid:1) ( t, x ) (cid:1) τ (cid:101) ρ (cid:1) ( t, x ) d x d t = ˆ T ˆ R d ∇ φ ( t, x ) · (cid:101) v (cid:1) ( t, x ) (cid:101) ρ (cid:1) ( t, x ) d x d t + e (cid:1) [ φ ]= ˆ Z T ∇ φ ( t, x ) · v d (cid:101) γ (cid:1) ( t, x, v ) + e (cid:1) [ φ ] (cid:1) → −→ ˆ Z T ∇ φ ( t, x ) · v d γ ∗ ( t, x, v )= ˆ [0 ,T ] × R d ∇ φ ( t, x ) · (cid:20) ˆ R d v d γ ( t,x ) ( v ) (cid:21) d M ∗ ( t, x )= ˆ T ˆ R d ∇ φ ( t, x ) · v ∗ ( t, x ) ρ ∗ ( t, x ) d x d t. The error term e (cid:1) [ φ ] above is controlled via Taylor expansion of φ and by using (4.3), (cid:12)(cid:12) e (cid:1) [ φ ] (cid:12)(cid:12) ≤ ˆ T ˆ R d τ (cid:107) φ (cid:107) C (cid:13)(cid:13)(cid:101) v (cid:1) ( t, x ) (cid:13)(cid:13) (cid:101) ρ (cid:1) ( t, x ) d x d t ≤ E(cid:107) φ (cid:107) C T τ.
Equality of the limits for both evaluations of J (cid:1) [ φ ] for arbitrary test functions φ shows thecontinuity equation (4.5). (cid:3) Unfortunately, the convergence provided by Theorem 4.2 is generally not sufficient to concludethat ρ ∗ is a weak solution to (1.1), since we are not able to identify v ∗ as v [ ρ ∗ ] from (1.3b). Theproblem is two-fold: first, weak- (cid:63) convergence of (cid:101) ρ (cid:1) is insufficient to pass to the limit insidethe nonlinear function P . Second, even if we would know that, for instance, P ( (cid:101) ρ (cid:1) ) ∗ (cid:42) P ( ρ ∗ ),we would still need a (cid:1) -independent a priori control on the regularity (e.g., maximal diameterof triangles) of the meshes generated by the G n (cid:1) to justify the passage to limit in the weakformulation below.The main difficulty in the weak formulation that we derive now is that we can only use “testfunctions” that are piecewise affine with respect to the changing meshes generated by the G n (cid:1) .For definiteness, we introduce the space D ( T ) := (cid:8) Γ : K → R d ; Γ is globally continuous, and is piecewise affine w.r.t. ∆ m (cid:9) . AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 13
Lemma 4.4.
Assume S : R d → R d is such that S ◦ G n (cid:1) ∈ D ( T ) . Then: ˆ R d P ( ρ n (cid:1) ) ∇ · S d x − ˆ R d ∇ V · S ρ n (cid:1) d x = ˆ R d S · v n (cid:1) ρ n (cid:1) d x. (4.7) Proof.
For all sufficiently small ε >
0, let G ε = (id+ S ) ◦ G n (cid:1) . By definition of G n (cid:1) as a minimizer,we have that E (cid:1) ( G ε ; G n − (cid:1) ) ≥ E (cid:1) ( G n (cid:1) ; G n − (cid:1) ). This implies that0 ≤ ε ˆ K (cid:18) τ (cid:2) (cid:107) G ε − G n − (cid:1) (cid:107) − (cid:107) G n (cid:1) − G n − (cid:1) (cid:107) (cid:3) + (cid:20)(cid:101) h (cid:18) det D G ε ρ T (cid:19) − (cid:101) h (cid:18) det D G n (cid:1) ρ T (cid:19)(cid:21) + (cid:2) V ◦ G ε − V (cid:3)(cid:19) ρ T d ω. (4.8)We discuss limits of the three terms under the integral for ε (cid:38)
0. For the metric term:12 τ ε (cid:2) (cid:107) G ε − G n − (cid:1) (cid:107) − (cid:107) G n (cid:1) − G n − (cid:1) (cid:107) (cid:3) = G n (cid:1) − G n − (cid:1) τ · G ε − G n (cid:1) ε + 12 τ ε (cid:107) G ε − G n (cid:1) (cid:107) = (cid:20)(cid:18) id − T n (cid:1) τ (cid:19) · S (cid:21) ◦ G n (cid:1) + ε τ (cid:107) S (cid:107) ◦ G n (cid:1) , and since S is bounded, the last term vanishes uniformly on K for ε (cid:38)
0. For the internal energy,since D G ε = D(id + εS ) ◦ G n (cid:1) · D G n (cid:1) , and recalling (3.8),1 ε (cid:20)(cid:101) h (cid:18) det D G ε ρ T (cid:19) − (cid:101) h (cid:18) det D G n (cid:1) ρ T (cid:19)(cid:21) = 1 ε (cid:20)(cid:101) h (cid:18) det D G n (cid:1) ρ T det( + ε D S ) ◦ G n (cid:1) (cid:19) − (cid:101) h (cid:18) det D G n (cid:1) ρ T (cid:19)(cid:21) ε (cid:38) −→ det D G n (cid:1) ρ T (cid:101) h (cid:48) (cid:18) det D G n (cid:1) ρ T (cid:19) (cid:18) lim ε (cid:38) det( + ε D S ) ε (cid:19) ◦ G n (cid:1) = − det D G n (cid:1) ρ T P (cid:18) ρ T det D G n (cid:1) (cid:19) tr[D S ] ◦ G n (cid:1) = − det D G n (cid:1) ρ T (cid:2) P ( ρ n ) ∇ · S (cid:3) ◦ G n (cid:1) . Since the piecewise constant function det D G n (cid:1) has a positive lower bound, the convergence as ε (cid:38) K . Finally, for the potential energy,1 ε (cid:2) V ◦ (id + εS ) ◦ G n (cid:1) − V ◦ G n (cid:1) (cid:3) ε (cid:38) −→ (cid:2) ∇ V · S (cid:3) ◦ G n (cid:1) . Again, the convergence is uniform on K . Passing to the limit in the integral (4.8) yields0 ≤ ˆ K (cid:20)(cid:18) id − T n (cid:1) τ (cid:19) · S (cid:21) ◦ G n (cid:1) ρ T d ω − ˆ K (cid:2) P ( ρ n ) ∇ · S (cid:3) ◦ G n (cid:1) det D G n (cid:1) d ω + ˆ K (cid:2) ∇ V · S (cid:3) ◦ G n (cid:1) ρ T d ω. The same inequality is true with − S in place of S , hence this inequality is actually an equality.Since ρ n (cid:1) = ( G n (cid:1) ) ρ T , a change of variables x = S n (cid:1) ( ω ) produces (4.7). (cid:3) Corollary 4.5.
In addition to the hypotheses of Theorem 4.2, assume that (1) P ( (cid:101) ρ (cid:1) ) ∗ (cid:42) p ∗ in L ([0 , T ] × Ω) ; (2) each G n (cid:1) is injective; (3) as (cid:1) → , all simplices in the images of T under G n (cid:1) have non-degenerate interiorangles and tend to zero in diameter, uniformly w.r.t. n . Then ρ ∗ satisfies the PDE ∂ t ρ ∗ = ∆ p ∗ + ∇ · ( ρ ∗ ∇ V ) (4.9) in the sense of distributions.Proof. Let a smooth test function ζ ∈ C ∞ c ( R d → R d ) be given. For each (cid:1) and each n , a ζ n (cid:1) : R d → R d with ζ n (cid:1) ◦ G n (cid:1) ∈ D ( T ) can be constructed in such a way that ζ n (cid:1) → ζ, ∇ · ζ n (cid:1) → ∇ · ζ (4.10)uniformly on R d , and uniformly in n as (cid:1) →
0. This follows from our hypotheses on the (cid:1) -uniform regularity of the Lagrangian meshes: inside the image of G n (cid:1) , one can simply choose ζ n (cid:1) as the affine interpolation of the values of ζ at the points G n (cid:1) ( ω (cid:96) ). Outside, one can take anarbitrary approximation of ζ that is compatible with the piecewise-affine approximation on theboundary of G n (cid:1) ’s image; one may even choose ζ n (cid:1) ≡ ζ 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 − τ, nτ ], substitute S ( t, x ) := η ( t ) ζ n (cid:1) ( x )into (4.7). Integration of these equalities with respect to t ∈ (0 , T ) yields ˆ T ˆ R d P ( (cid:101) ρ (cid:1) ) ∇ · S d x d t − ˆ T ˆ R d ∇ V · S d x d t = ˆ T ˆ R d S · (cid:101) v (cid:1) (cid:101) ρ (cid:1) d x d t. We pass to the limit (cid:1) → P ( (cid:101) ρ ) ∗ (cid:42) p ∗ byhypothesis, for the last, we use Theorem 4.2 above. Since any test function S ∈ C ∞ c ((0 , T ) × Ω)can be approximated in C by linear combinations of products η ( t ) ζ ( x ) as above, we thus obtainthe weak formulation of ρ ∗ v ∗ = ∇ p ∗ + ρ ∗ ∇ V. In combination with the continuity equation (4.5), we arrive at (4.9). (cid:3)
Remark 4.6.
In principle, our discretization can also be applied to the linear
Fokker-Planckequation with P ( r ) = r and h ( r ) = r log r . In that case, one automatically has P ( (cid:101) ρ ) ∗ (cid:42) p ∗ ≡ P ( ρ ∗ ) thanks to Theorem 4.2. Corollary 4.5 above then provides an a posteriori criterion forconvergence: if the Lagrangian mesh does not deform too wildly under the dynamics as thediscretization 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 certainconditions on the spatial discretization T , any smooth and positive solution ρ to the initial valueproblem (1.1) projects to a discrete solution that satisfies the Euler-Lagrange equations up to acontrolled error. We restrict ourselves to d = 2 dimensions.5.1. Smooth Lagrangian evolution.
First, we derive an alternative form of the velocity field v from (1.3b) in terms of G . Lemma 5.1.
For ρ = G ρ with a smooth diffemorphism G : K → R d , we have v [ ρ ] ◦ G = V [ G ] := P (cid:48) (cid:18) ρ det D G (cid:19) (D G ) − T (cid:18) tr (cid:2) (D G ) − D G (cid:3) T − ∇ ρρ (cid:19) − ∇ V ◦ G. (5.1) Consequently, the Lagrangian map G — relative to the reference density ρ — for a smoothsolution ρ to (1.1) satisfies ∂ t G = V [ G ] . (5.2) AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 15
Proof.
On the one hand, D (cid:2) h (cid:48) ( ρ ) ◦ G (cid:3) = (cid:2) D h (cid:48) ( ρ ) (cid:3) ◦ G D G, and on the other hand, by definition of the push forward,D (cid:2) h (cid:48) ( ρ ) ◦ G (cid:3) = D h (cid:48) (cid:18) ρ det D G (cid:19) = h (cid:48)(cid:48) (cid:18) ρ det D G (cid:19) (cid:18) ρ det D G (cid:19) (cid:18) D ρρ − tr (cid:2) (D G ) − D G (cid:3)(cid:19) = (cid:2) ρh (cid:48)(cid:48) ( ρ ) (cid:3) ◦ G (cid:18) D ρρ − tr (cid:2) (D G ) − D G (cid:3)(cid:19) . Hence ∇ h (cid:48) ( ρ ) ◦ G = (cid:2) ρh (cid:48)(cid:48) ( ρ ) (cid:3) ◦ G (D G ) − T (cid:18) ∇ ρρ − tr (cid:2) (D G ) − D G (cid:3) T (cid:19) . Observing that (1.2) implies that rh (cid:48)(cid:48) ( r ) = P (cid:48) ( r ), we conclude (5.2) directly from (1.3b). (cid:3) Discrete Euler-Lagrange equations in dimension d = 2 . In the planar case d = 2, theEuler-Lagrange equation (3.10) above can be rewritten in a more convenient way.In the following, fix some vertex ω × of the triangulation, which is indicent to precisely sixtriangles. For convenience, we assume that these are labelled ∆ to ∆ in counter-clockwiseorder. Similarly, the six neighboring vertices are labeled ω to ω in counter-clockwise order, sothat ∆ k has vertices ω k and ω k +1 , where we set ω := ω .Using these conventions and recalling Lemma B.2, the expression for the vector ν in (3.9)simplifies to ν k T = − J ( G k +1 − G k ) , where J = (cid:18) −
11 0 (cid:19) . Summing the Euler-Lagrange equation (3.10) over ∆ to ∆ , we obtain p × = J × , (5.3)where the momentum term p × and the impulse J × , respectively, are given by p × = 112 (cid:88) k =0 µ k T (cid:20) (cid:18) G × − G ∗× τ (cid:19) + (cid:18) G k − G ∗ k τ (cid:19) + (cid:18) G k +1 − G ∗ k +1 τ (cid:19)(cid:21) (5.4) J × = (cid:88) k =0 µ k T (cid:20) µ k T P (cid:18) µ k T det( G k − G × | G k +1 − G × ) (cid:19) J ( G k +1 − G k ) (5.5) − (cid:52) ∇ V (cid:0) (1 − ξ − ξ ) G × + ξ G k + ξ G k +1 (cid:1) (1 − ξ − ξ ) d ξ (cid:21) . (5.6)We shall now prove our main result on consistency. The setup is the following: a sequence oftriangulations T ε on K , parametrized by ε >
0, and a sequence of time steps τ ε = O ( ε ) are given.We assume that there is an ε -independent region K (cid:48) ⊂ K on which the T ε are almost hexagonal in the following sense: each node ω × ∈ K (cid:48) of T ε has precisely six neighbors — labelled ω to ω in counter-clockwise order — and there exists a rotation R ∈ SO(2) such that R ( ω k − ω × ) = εσ k + O ( ε ) with σ k = (cid:18) cos π k sin π k (cid:19) (5.7)for k = 0 , , . . . , Now, let G : [0 , T ] × K → R d be a given smooth solution to the Lagrangian evolution equa-tion (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 ε ( ω (cid:96) ) = G ( t ; ω (cid:96) ) and G ∗ ε ( ω (cid:96) ) = G ( t − τ ; ω (cid:96) ), at all nodes ω (cid:96) in T ε . Theorem 5.2 below statesthat 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 (cid:48) .The hexagonality hypothesis on the T ε is strong, but some very strong restriction of A T ε ’sgeometry 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-Lagrangeequation (5.3) admits the following asymptotic expansion: p × = √ ε ρ ( ω × ) ∂ t G ( t ; ω × ) + O ( ε ) , (5.8a) J × = √ ε ρ ( ω × ) V [ G ]( t ; ω × ) + O ( ε ) , (5.8b) as ε → , uniformly at the nodes ω × ∈ K (cid:48) of the respective T ε . Remark 5.3.
Up to an error O ( ε ) , the geometric pre-factor √ ε equals to one third of thetotal area of the hexagon with vertices ω to ω , and is thus equal to the integral of the piecewiseaffine hat function with peak at ω × .Proof of Theorem 5.2. Throughout the proof, let ε > ε -index for T ε and τ ε . First, we fix a node ω × of T ∩ K (cid:48) . 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, (cid:88) k =0 σ k = 0 , (cid:88) k =0 ω k = O ( ε ) . (5.9)Moreover, we have that | ∆ k | = det( ω k | ω k +1 ) = ε det( σ k | σ k +1 ) + O ( ε ) = √ ε + O ( ε ) . (5.10)On the other hand, by definition of µ k T in (3.1), it follows that µ k T = | ∆ k | ∆ k ρ d ω = 12 det( ω k | ω k +1 ) (cid:20) ρ (cid:18) ω k + ω k +1 (cid:19) + O ( ε ) (cid:21) = 12 det( ω k | ω k +1 ) (cid:20) ρ × + ε ∇ ρ × · σ k + σ k +1 O ( ε ) (cid:21) . (5.11)Combining (5.10) and (5.11) yields µ k T = ε (cid:32) √ ρ × + O ( ε ) (cid:33) . (5.12)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 , . . . , G × = D G ( t, ω × ), D G × = D G ( t, ω × ), ∂ t G × = ∂ t G ( t, ω × ).To perform an expansion in the momentum term , first observe that G ( t − τ ; ω k ) = G ( t ; ω k ) − τ ∂ t G ( t ; ω k ) + O ( τ ) , AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 17 for each k = 0 , , . . . ,
5, and so, using that τ = O ( ε ) by hypothesis, G k − G ∗ k τ = ∂ t G ( t ; ω k ) + O ( τ ) = ∂ t G × + O ( ε ) + O ( τ ) = ∂ t G × + O ( ε ) . Using (5.12) and then (5.9) yields p × = 112 τ (cid:88) k =0 ε (cid:32) √ ρ × + O ( ε ) (cid:33) (cid:2) ∂ t G × + O ( ε ) (cid:3) = √ ε ρ × ∂ t G × + O ( ε ) . This is (5.8a).For the impulse term , we start with a Taylor expansion to second order in space: G k = G × + D G × ω k + 12 D G × : [ ω k ] + O ( ε ) . We combine this with the observation that ( ω k | ω k +1 ) − = O ( ε − ) to obtain: µ k T det( G k − G × | G k +1 − G × )= det( ω k | ω k +1 )det D G × ρ × + ε ∇ ρ × · σ k + σ k +1 + O ( ε )det (cid:2) ( ω k | ω k +1 ) + (D G × ) − (cid:0) D G × : [ ω k ] (cid:12)(cid:12) D G × : [ ω k +1 ] (cid:1) + O ( ε ) (cid:3) = ρ × det D G × ε ∇ ρ × ρ × · σ k + σ k +1 O ( ε )det (cid:2) + (D G × ) − (cid:0) D G × : [ ω k ] (cid:12)(cid:12) D G × : [ ω k − ] (cid:1) ( ω k | ω k +1 ) − + O ( ε ) (cid:3) = ρ × det D G × (cid:18) ε (cid:26) χ k − ϑ k (cid:27) + O ( ε ) (cid:19) , where χ k = ∇ ρ × ρ × · σ k + σ k +1 ,ϑ k = tr (cid:2)(cid:0) (D G × ) − D G × : [ σ k ] (cid:12)(cid:12) (D G × ) − D G × : [ σ k +1 ] (cid:1) ( σ k | σ k +1 ) − (cid:3) . Plugging this in leads to (cid:88) k =0 (cid:26) P (cid:18) ρ × det D G × (cid:19) + ε P (cid:48) (cid:18) ρ × det D G × (cid:19) (cid:26) χ k − ϑ k (cid:27) + O ( ε ) (cid:27) J D G × ( ω k +1 − ω k )= 12 P (cid:18) ρ det D G × (cid:19) J D G × (cid:32) (cid:88) k =0 ( ω k +1 − ω k ) (cid:33) + ε P (cid:48) (cid:18) ρ × det D G × (cid:19) J D G × J T (cid:32) (cid:88) k =0 { χ k − ϑ k } J ( σ k +1 − σ k ) (cid:33) + O ( ε )= 0 + √ ε P (cid:48) (cid:18) ρ × det D G × (cid:19) (D G × ) − T (cid:26) tr (cid:2) (D G × ) − D G × (cid:3) T − ∇ ρ × ρ × (cid:27) + O ( ε ) , 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: ∇ V ( g ) = ∇ V ( G × ) + O ( ε ) holds for any g that is a convex combination of G × , G , . . . , G , where the implicit constant iscontrolled in terms of the supremum of D V and D G on K (cid:48) . With that, we simply have, usingagain (5.12): (cid:88) k =0 µ k T (cid:52) ∇ V (cid:0) (1 − ξ − ξ ) G × + ξ G k + ξ G k +1 (cid:1) (1 − ξ − ξ ) d ξ = 6 ε (cid:32) √ ρ × + O ( ε ) (cid:33) (cid:0) ∇ V ( G × ) + O ( ε ) (cid:1) = √ ε ρ × ∇ V ( G × ) + O ( ε ) . Together, this yields (5.8b). (cid:3)
Remark 5.4.
The hypotheses of Theorem (5.2) require that the T ε are almost hexagonal on K (cid:48) .This seems like a technical hypothesis that simplifies calculations, but apparently, some strongsymmetry 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 here —assume that V ≡ and ρ ≡ , and consider a sequence of triangulations T ε for which there isa node ω × such that (5.7) holds with the σ k being replaced by a different six-tuple of vectors σ (cid:48) k .Repeating the steps of the proof above, it is easily seen that p × = aε ∂ t G ( t ; ω × ) + O ( ε ) , withan ε -independent constant a > in place of √ / , and that J × = − ε P (cid:48) (cid:18) G × (cid:19) (D G × ) − T (cid:88) k =0 ϑ (cid:48) k J ( σ (cid:48) k +1 − σ (cid:48) k ) + O ( ε ) , with ϑ (cid:48) k = tr (cid:2)(cid:0) (D G × ) − D G × : [ σ (cid:48) k ] (cid:12)(cid:12) (D G × ) − D G × : [ σ (cid:48) k +1 ] (cid:1) ( σ (cid:48) k | σ (cid:48) k +1 ) − (cid:3) . If a result of the form (5.8b) — with √ / replaced by a — was true, then this implies inparticular that (cid:88) k =0 ϑ (cid:48) k J ( σ (cid:48) k +1 − σ (cid:48) k ) = a (cid:48) tr (cid:2) (D G × ) − D G × (cid:3) (5.13) holds with some constant a (cid:48) > for arbitrary matrices D G × ∈ R × of positive determinantand tensors D G × ∈ R × × that are symmetric in the second and third component. A specificexample for which (5.13) is not true is given by σ (cid:48) = (cid:18) (cid:19) = − σ (cid:48) , σ (cid:48) = (cid:18) (cid:19) = − σ (cid:48) , σ (cid:48) = (cid:18) (cid:19) = − σ (cid:48) , (5.14) in combination with D G × = , and a D G × that is zero except for two ones, at the positions (1 , , and (2 , , . In Lemma B.5, we show that the left-hand side in (5.13) equals to (cid:0) (cid:1) ; onthe other hand, the right-hand side is clearly zero.Note that this counter-example is significant, insofar as the skew (in fact, degenerate) hexagondescribed by the σ (cid:48) k in (5.14) corresponds to a popular method for triangulation of the plane. Numerical simulations in d = 26.1. Implementation.
The Euler-Lagrange equations for the d = 2-dimensional case have beenderived in (5.3). We perfom a small modification in the potential term in order to simplify AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 19 calculations with presumably minimal loss in accuracy: Z × [ G ; G ∗ ] = (cid:88) k =0 µ k T (cid:20) (cid:18) G × − G ∗× τ (cid:19) + (cid:18) G k − G ∗ k τ (cid:19) + (cid:18) G k +1 − G ∗ k +1 τ (cid:19)(cid:21) + (cid:88) k =0 (cid:20) (cid:101) h (cid:48) (cid:18) det( G k − G × | G k +1 − G × )2 µ k T (cid:19) J ( G k +1 − G k ) + µ k T ∇ V ( G k + ) (cid:21) , with the short-hand notation G k + = 13 ( G × + G k + G k +1 ) . On the main diagonal, the Hessian amounts to H ×× [ G ] = (cid:32) (cid:88) k =0 µ k T τ (cid:33) + (cid:88) k =0 µ k T (cid:101) h (cid:48)(cid:48) (cid:18) det( G k − G × | G k +1 − G × )2 µ k T (cid:19) (cid:2) J ( G k +1 − G k ) (cid:3)(cid:2) J ( G k +1 − G k ) (cid:3) (cid:62) + (cid:88) k =0 µ k T ∇ V ( G k + )Off the main diagonal, the entries of the Hessian are given by H × k [ G ] = µ k T + µ k − T τ + 14 µ k T (cid:101) h (cid:48)(cid:48) (cid:18) det( G k − G × | G k +1 − G × )2 µ k T (cid:19) (cid:2) J ( G k +1 − G k ) (cid:3)(cid:2) J ( G k +1 − G × ) (cid:3) (cid:62) − µ k − T (cid:101) h (cid:48)(cid:48) (cid:32) det( G k − − G × | G k − G × )2 µ k − T (cid:33) (cid:2) J ( G k − G k − ) (cid:3)(cid:2) J ( G k − − G × ) (cid:3) (cid:62) + µ k T ∇ V ( G k + ) + µ k − T ∇ V ( G k − ) . The scheme consists of an inner (Newton) and an outer (time stepping) iteration. We startfrom a given initial density ρ and define the solution at the next time step inductively byapplying 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 G ( s +1) := G ( s ) + δG ( s +1) , where the update δG ( s +1) is the solution to the linear system H [ G ( s ) ] δG ( s +1) = − Z [ G ( s ) ; G n ] . The effort of each inner iteration step is essentially determined by the effort to invert the sparsematrix 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 − .6.2. Numerical experiments.
In this section we present results of our numerical experimentsfor (1.1) with a cubic porous-medium nonlinearity P ( r ) = r and different choices for the externalpotential V , ∂ t ρ = ∆( u ) + ∇ · ( u ∇ V ) . (6.1) Numerical experiment 1: unconfined evolution of Barenblatt profile.
As a first example, we con-sider the “free” cubic porous medium equation, that is (6.1) with V ≡
0. It is well-known (see,e.g., Vazquez [33]) that in the long-time limit t → ∞ , arbitrary solutions approach a self-similarone, ρ ∗ ( t, x ) = t − dα B (cid:0) t − α x (cid:1) with α = 16 , (6.2)where B is the associated Barenblatt profile B ( z ) = (cid:18) C − (cid:107) z (cid:107) (cid:19) + , (6.3)where C = (2 π ) − ≈ .
29 is chosen to normalize B ’s mass to unity.In this experiment, we are only interested in the quality of the numerical approximationfor the self-similar solution (6.2). To reduce numerical effort, we impose a four-fold symmetryof the approximation: we use the quarter circle as computational domain K , and interpretethe discrete function thereon as one of four symmetric pieces of the full discrete solution. Topreserve reflection symmetry over time, homogeneous Neumann conditions are imposed on theartificial boundaries. This is implemented by reducing the degrees of freedom of the nodes alongthe x - and y -axes to tangential motion. We initialize our simulation with a piecewise constantapproximation of the profile of ρ ∗ from (6.3) at time t = 0 .
01. We choose a time step τ = 0 . T = 2. In Figure 1, we have collected snapshots of the approximated density Figure 1.
Numerical experiment 1: fully discrete evolution of our approxima-tion for the self-similar solution to the free porous medium equation. Snapshotsare taken at times t = 0 . t = 0 . t = 0 .
25, and t = 2 . AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 21
Figure 2.
Numerical experiment 1: comparison of the discrete solution (inter-polated surface plots with triangulation) with the Barenblatt profile (solid anddashed black lines along the identity) at different times. -2 -1 time -1 e n e r g y -2 -1 h max -3 -2 -1 l - e rr o r Figure 3.
Numerical experiment 1: decay of the energy of the discrete solutionin comparison with the analytical decay t − / of the Barenblatt solution (left).Numerical convergence for fixed ratio τ /h = 0 . Remark 6.1.
It takes less than 2 minutes to complete this simulation on standard laptop (
Mat-lab 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 withthe Barenblatt profile at the respective time. By construction of the scheme, the initial mass isexactly conserved in time as the discrete solution propagates. The left plot in Figure 3 showsthe decay in the energy and gives quantitative information about the difference of the discretesolution to the analytical Barenblatt solution. The numerical solution shows good agreementwith the analytical energy decay rate c = 2 /2 JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK
Figure 2 shows surface plots of the discrete solution at different times in comparison withthe Barenblatt profile at the respective time. By construction of the scheme, the initial mass isexactly conserved in time as the discrete solution propagates. The left plot in Figure 3 showsthe decay in the energy and gives quantitative information about the difference of the discretesolution to the analytical Barenblatt solution. The numerical solution shows good agreementwith the analytical energy decay rate c = 2 /2 JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK We also compute the l -error of the discrete solution to the exact Barenblatt profile andobserve that it remains within the order of the fineness of the triangulation. The mass of thediscrete 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 aboveinitial data on different meshes. We fix the ratio τ /h = 0 . l -error attime T = 0 . h max = 0 . , . , . , . . We expect the error to decayas a power of h max . The double logarithmic plot should reveal a line with its slope indicatingthe numerical convergence order. The right plot in Figure 3 shows the result, the estimatednumerical convergence order which is obtained from a least-squares fitted line through the pointsis equal to 1 .
18. This indicates first order convergence of the scheme with respect to the spatialdiscretisation parameter h max . Numerical experiment 2: Asymptotic self-similarity.
In our second example, we are still con-cerned with the free cubic porous medium equation, (6.1) with V ≡
0. This time, we wish togive an indication that the discrete approximation of the self-similar solution from (6.2) from theprevious experiment might inherit the global attractivity of its continuous counterpart. Morespecifically, we track the discrete evolution for the initial datum ρ ( x, y ) = 3000( x + y ) exp[ − | x | + | y | )] + 0 . T = 0 . 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 ρ ( x, y ) = exp[ − x − . + ( y − . )] + exp[ − x + 0 . + ( y + 0 . )] + 0 . . (6.5)We choose a triangulation of the square [ − . , . and initialise the discrete solution piecewiseconstant in each triangle, with a value corresponding to (6.5), evaluated in the centre of mass ofeach 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 + y ) /
2. The time step is τ = 0 .
001 and the final time is T = 0 . . Figure 5 shows the evolution from the initial density. As time increases the peaks smoothlymerge into each other. As the thin layer around the peaks is also subject to the potential thetriangulated domain shrinks in time. Even if we do not know how to prevent theoretically theintersection of the images of the discrete Lagrangian maps, this seems not to be a problem inpractice. As time evolves, the discrete solution approaches the steady state Barenblatt profilegiven by B ( z ) = (cid:18) C − || z || (cid:19) + , (6.6)where C is chosen as the mass of the density. The plot in Figure 6 shows the exponential decayof the l -distance of the discrete solution to the steady state Barenblatt profile (6.6). We observethat the decay agrees very well with the analytically predicted decay exp( − t ) until t = 0 . Numerical experiment 4: one peak splitting under the influence of a quartic potential.
We con-sider as the initial condition ρ ( x, y ) = 1 − ( x + y ) . (6.7)We choose a triangulation of the unit circle and initialise the discrete solution piecewise constantin each triangle, with a value corresponding to (6.7), evaluated in the centre of mass of each AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 23
Figure 4.
Numerical experiment 2: fully discrete evolution for the initial den-sity from (6.4) under the free porous medium equation. Snapshots are taken attimes t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . P ( r ) = r m and V ( x ) = 5( x + (1 − y ) ) /
2. The time step is τ = 0 .
005 and the final time is T = 0 . . Figure 7 shows the evolution of the initial density. As time increases the initial density isprogressively split, until two new maxima emerge which are connected by a thin layer. For largertimes, when certain triangles become excessively distorted, one would monitor triangle qualitynumerically, and re-mesh, locally refining the triangulation where necessary.
Figure 5.
Numerical experiment 3: evolution of two peaks merging under theporous medium equation with a confining potential.
Appendix A. Proof of the Lagrangian representation
Proof of Lemma 1.1.
We verify that the density function given by ( G − t ) ρ t on K ⊂ R d isconstant with respect to time t ; the identity (1.5) then follows since ρ t = ( G t ◦ G − t ) ρ t = ( G t ) (cid:2) ( G − t ) ρ t (cid:3) = ( G t ) (cid:2) ( G − ) ρ (cid:3) = ( G t ) ρ. Firstly, from the definition of the inverse, G − t ◦ G t = id AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 25 time l - d i s t a n ce Figure 6.
Numerical experiment 3: two merging peaks: plot of the l -distanceof the discrete solution to the steady state Barenblatt profile in comparison withthe analytical decay c exp( − t ). Figure 7.
Numerical experiment 4: evolution of the initial density under theporous medium equation with a quartic potential.for all t , differentiating with respect to time yieldsD( G − t ) ◦ G t ∂ t G t + ∂ t ( G − t ) ◦ G t = 0 , and so, using (1.4) and (1.3b), ∂ t ( G − t ) = − D( G − t )( ∂ t G t ◦ G − t ) = − D( G − t ) v [ ρ t ] . (A.1) Now, let ϕ be a smooth test function, and considerdd t ˆ ϕ ( G − t ) ρ t = dd t ˆ ( ϕ ◦ G − t ) ρ t = ˆ ( ϕ ◦ G − t ) ∂ t ρ t + ˆ D ϕ ◦ G − t ∂ t ( G − t ) ρ t = − ˆ ( ϕ ◦ G − t )[ ∇ · ( ρ t v ( ρ t ))] − ˆ (D ϕ ◦ G − t ) D( G − t ) v ( ρ t ) ρ t by (1.1) and (1.4)= ˆ (D ϕ ◦ G − t )D( G − t ) [ v ( ρ t ) − v ( ρ t )] ρ t integrating by parts= 0 . As ϕ was arbitrary, ( G − t ) ρ t is constant with respect to time. (cid:3) Appendix B. Technical lemmas
Lemma B.1.
Given g , g , . . . , g d ∈ R d , then (cid:52) d (cid:13)(cid:13)(cid:13) g + d (cid:88) j =1 ω j ( g j − g ) (cid:13)(cid:13)(cid:13) d ω = 2( d + 1)( d + 2) (cid:88) ≤ i ≤ j ≤ d g i · g j . (B.1) Proof.
Thanks to the symmetry of the integral with respect to the exchange of the components ω j , the left-hand side of (B.1) equals to (cid:107) g (cid:107) + 2 (cid:18) (cid:52) ω d d ω (cid:19) (cid:88) ≤ j ≤ d g · ( g j − g )+ (cid:18) (cid:52) ω d d ω (cid:19) (cid:88) ≤ j ≤ d (cid:107) g j − g (cid:107) + 2 (cid:18) (cid:52) ω d − ω d d ω (cid:19) (cid:88) ≤ i Third integral: (cid:52) ω d − ω d d ω = 1 |(cid:52) d | ˆ (cid:20) ˆ − ω d ω d − ω d (1 − ω d − − ω d ) d − |(cid:52) d − | d ω d − (cid:21) d ω d = |(cid:52) d − ||(cid:52) d | ˆ (cid:20) ˆ z (1 − z )( z − y ) y d − d y (cid:21) d z = d ( d − ˆ (cid:20) d − − d (cid:21) (1 − z ) z d d z = 1 d + 1 − d + 2 = 1( d + 1)( d + 2) . Substitute this into (B.2): (cid:18) − d + 1 + d + d ( d + 1)( d + 2) (cid:19) (cid:107) g (cid:107) + (cid:18) d + 1 − d + 2( d + 1)( d + 2) (cid:19) (cid:88) ≤ j ≤ d g · g j + 2( d + 1)( d + 2) (cid:88) ≤ j ≤ d (cid:107) g j (cid:107) + 2( d + 1)( d + 2) (cid:88) ≤ i Lemma B.2. For each A ∈ R × , we have J A J T = (det A ) A − T .Proof. This is verified by direct calculation: J A J T = (cid:18) − 11 0 (cid:19) (cid:18) a a a a (cid:19) (cid:18) − (cid:19) = (cid:18) a − a − a a (cid:19) = (det A ) A − T . (cid:3) Lemma B.3. With σ k ∈ R defined as in (5.7) , we have that (cid:88) k =0 J ( σ k − σ k +1 ) (cid:18) σ k + σ k +1 (cid:19) T = √ . Proof. With the abbreviations φ x = π x and ψ = π : (cid:88) k =0 J ( σ k − σ k +1 ) (cid:18) σ k + σ k +1 (cid:19) T = 13 (cid:88) k =0 (cid:18) sin φ k +1 − sin φ k cos φ k − cos φ k +1 (cid:19)(cid:18) cos φ k + cos φ k +1 sin φ k + sin φ k +1 (cid:19) T = 13 (cid:88) k =0 (2 sin ψ (cid:18) cos φ k + sin φ k + (cid:19) (2 cos ψ (cid:18) cos φ k + sin φ k + (cid:19) T = sin ψ (cid:88) k =0 (cid:32) φ k + φ k + sin φ k + φ k + sin φ k + φ k + (cid:33) = √ (cid:88) k =0 (cid:20) + (cid:18) cos φ k +1 sin φ k +1 sin φ k +1 − cos φ k +1 (cid:19)(cid:21) = √ . (cid:3) Lemma B.4. Let the scheme B := ( b pqr ) p,q,r ∈{ , } ∈ R × × of eight numbers b pqr ∈ R besymmetric in the last two indices, b pqr = b prq . With σ k ∈ R defined as in (5.7) , we have that (cid:88) k =0 tr (cid:2)(cid:0) σ k (cid:12)(cid:12) σ k +1 (cid:1) − (cid:0) B : [ σ k ] (cid:12)(cid:12) B : [ σ k +1 ] (cid:1)(cid:3) J ( σ k − σ k +1 ) = 2 √ [ B ] T . (B.3) Proof. In principle, this lemma can be verified by a direct calculation, by writing out the sixterms in the sum explicitly and using trigonometric identities. Below, we give a slightly moreconceptual proof, in which we use symmetry arguments to reduce the number of expressionssignificantly.For the matrix involving B , we obtain (cid:0) B : [ σ k ] (cid:12)(cid:12) B : [ σ k +1 ] (cid:1) = (cid:18) b σ k, + b σ k, + 2 b σ k, σ k, b σ k +1 , + b σ k +1 , + 2 b σ k +1 , σ k +1 , b σ k, + b σ k, + 2 b σ k, σ k, b σ k +1 , + b σ k +1 , + 2 b σ k +1 , σ k +1 , (cid:19) , while clearly (cid:0) σ k (cid:12)(cid:12) σ k +1 (cid:1) − = 2 √ (cid:18) σ k +1 , − σ k +1 , − σ k, σ k, (cid:19) . The sum of the diagonal entries of the matrix product are easily calculated, T k := tr (cid:2)(cid:0) σ k (cid:12)(cid:12) σ k +1 (cid:1) − (cid:0) B : [ σ k ] (cid:12)(cid:12) B : [ σ k +1 ] (cid:1)(cid:3) = 2 √ (cid:88) p,q,r =1 b pqr γ pqr,k , with the trigonometric expressions γ ,k = σ k, σ k +1 , − σ k +1 , σ k, , γ ,k = σ k, σ k +1 , − σ k +1 , σ k, ,γ ,k = γ ,k = σ k, σ k, σ k +1 , − σ k +1 , σ k +1 , σ k, ,γ ,k = σ k +1 , σ k, − σ k, σ k +1 , , γ ,k = σ k +1 , σ k, − σ k, σ k +1 , ,γ ,k = γ ,k = σ k +1 , σ k +1 , σ k, − σ k, σ k, σ k +1 , . To key step is to calculate the sum over k = 0 , , . . . , T k with the respectivevector η k = J ( σ k − σ k +1 ) = (cid:18) σ k +1 , − σ k, σ k, − σ k +1 , (cid:19) . Several simplifications of this sum can be performed, thanks to the particular form of the γ pqr,k and elementary trigonometric identities. First, observe that σ k +3 = − σ k , and hencethat γ pqr,k +3 = − γ pqr,k . Since further η k +3 = − η k , it follows that γ pqr,k +3 η k +3 = γ pqr,k η k . (B.4)Second, η can be evaluated explicitly for k = 1 , , η = 12 (cid:18) √ (cid:19) , η = (cid:18) (cid:19) , η = 12 (cid:18) −√ (cid:19) . (B.5)Third, since σ , = − σ , and σ , = − σ , , as well as σ , = σ , and σ , = σ , , we obtain that γ pqr, = 0 if p + q + r is odd , and γ pqr, = ( − p + q + r γ pqr, . (B.6)By putting this together, we arrive at (cid:88) k =0 γ pqr,k η k (B.4) = 2 (cid:88) k =0 γ pqr,k η k (B.5) = (cid:18) √ (cid:0) γ pqr, − γ pqr, (cid:1) γ pqr, + 2 γ pqr, + γ pqr, (cid:19) (B.6) = (cid:18) √ (cid:0) − ( − p + q + r (cid:1) γ pqr, (cid:0) − p + q + r (cid:1)(cid:0) γ pqr, + γ pqr, (cid:1)(cid:19) = (cid:18) √ γ pqr, (1 − e pqr )2 (cid:0) γ pqr, + γ pqr, (cid:1) e pqr (cid:19) , AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 29 where e pqr = 1 if p + q + r is even, and e pqr = 0 if p + q + r is odd. By elementary computations, p + q + r odd, k = 0 : γ , = √ , γ , = 0 , γ , = γ , = √ ; p + q + r even, k = 0 : γ , = − , γ , = , γ , = γ , = 0; p + q + r even, k = 1 : γ , = , γ , = , γ , = γ , = , and so the final result is: (cid:88) k =0 tr (cid:2)(cid:0) σ k (cid:12)(cid:12) σ k +1 (cid:1) − (cid:0) B : [ σ k ] (cid:12)(cid:12) B : [ σ k +1 ] (cid:1)(cid:3) J ( σ k − σ k +1 )= (cid:88) k =0 T k η k = 2 √ (cid:88) p,q,r =1 (cid:32) b pqr (cid:88) k =0 γ pqr,k η k (cid:33) = 2 √ (cid:18) b + b b + b (cid:19) , which is (B.3). (cid:3) Lemma B.5. With σ (cid:48) k ∈ R defined as in (5.14) , and with B = ( b pqr ) p,q,r ∈{ , } ∈ R × × suchthat b pqr = 0 except for b = b = 1 , we have that (cid:88) k =0 tr (cid:2)(cid:0) σ (cid:48) k (cid:12)(cid:12) σ (cid:48) k +1 (cid:1) − (cid:0) B : [ σ (cid:48) k ] (cid:12)(cid:12) B : [ σ (cid:48) k +1 ] (cid:1)(cid:3) J ( σ k − σ k +1 ) = − (cid:18) (cid:19) . (B.7) Proof. This is a slightly tedious, but straightforward calculation. First, by the choice of B , β k := (cid:0) B : [ σ (cid:48) k ] (cid:12)(cid:12) B : [ σ (cid:48) k +1 ] (cid:1) = (cid:18) ( σ (cid:48) k, ) ( σ (cid:48) k +1 , ) ( σ (cid:48) k, ) ( σ (cid:48) k +1 , ) (cid:19) , and so, by definition of the σ (cid:48) k in (5.14), β = β = (cid:18) (cid:19) , β = β = (cid:18) (cid:19) , β = β = (cid:18) (cid:19) . For the inverse matrices S k := (cid:0) σ (cid:48) k (cid:12)(cid:12) σ (cid:48) k +1 (cid:1) − , we obtain S = (cid:18) − 10 2 (cid:19) = − S , S = (cid:18) − (cid:19) = − S , S = (cid:18) − (cid:19) = − S . For the traces T k := tr (cid:2) S k β k (cid:3) , we thus obtain the values: T = T = − , T = T = 12 , T = T = 0 . In conclusion, (cid:88) k =0 T k J ( σ k − σ k +1 ) = J (cid:20) − 12 ( σ − σ ) + 12 ( σ − σ ) (cid:21) = J (cid:18) − (cid:19) = − (cid:18) (cid:19) , which is (B.7). (cid:3) Appendix C. Lack of convexity Below, we discuss why the minimization problem (3.7) is not convex. More precisely, we showthat G (cid:55)→ E (cid:1) ( G ; ˆ G ) is not convex as a function of G on the affine ansatz space A T . Since E (cid:1) ( G ; ˆ G ) is a convex combination of the expressions H m (cid:0) ( A m | b m ); ( ˆ A m | ˆ b m ) (cid:1) , it clearly sufficesto discuss the convexity of the latter. We consider a curve s (cid:55)→ ( A m + sα m | b m + sβ m ) and evaluate the second derivatives of thecomponents of the functional at s = 0. First,I := d d s (cid:12)(cid:12)(cid:12)(cid:12) s =0 (cid:18) τ ∆ m (cid:12)(cid:12) ( A m − ˆ A m + sα m ) ω + ( b m − ˆ b m ) + sβ m (cid:12)(cid:12) d ω (cid:19) = 1 τ ∆ m | α m ω + β m | d ω. Second, II := d d s (cid:12)(cid:12)(cid:12)(cid:12) s =0 ∆ m V (cid:0) ( A m + sα m ) ω + ( b m + sβ m ) (cid:1) d ω = ∆ m ( α m ω + β m ) T · ∇ V ( A m ω + b m ) · ( α m ω + β m ) d ω. If we assume that ∇ V ≥ λ , then we obtain for the sum of these two contributions thatI + II ≥ (cid:18) τ + λ (cid:19) ∆ m | α m ω + β m | d ω. For the remaining term, however, we obtain — using the abbreviations (cid:101) g ( s ) = s (cid:101) h (cid:48) ( s ) and (cid:101) f ( s ) = s (cid:101) g (cid:48) ( s ) — thatd d s (cid:12)(cid:12)(cid:12)(cid:12) s =0 (cid:101) h (cid:18) det( A m + sα m ) ρ m (cid:19) = dd s (cid:12)(cid:12)(cid:12)(cid:12) s =0 (cid:26)(cid:101) g (cid:18) det( A m + sα m ) ρ m (cid:19) tr (cid:2) ( A m + sα m ) − α m (cid:3)(cid:27) = (cid:101) f (cid:18) det A m ρ m (cid:19) (cid:0) tr (cid:2) A − m α m (cid:3)(cid:1) − (cid:101) g (cid:18) det A m ρ m (cid:19) tr (cid:2)(cid:0) A − m α m (cid:1) (cid:3) . Now observe that (cid:101) f ( s ) = P (cid:48) (1 /s ) − sP (1 /s ) is a non-negative, and (cid:101) 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 α m := A m (cid:18) (cid:19) , such that (cid:0) tr (cid:2) A − m α m (cid:3)(cid:1) = 0 , tr (cid:2)(cid:0) A − m α m (cid:1) (cid:3) = 2 , the sum is obviously negative. References [1] L. Ambrosio, N. Gigli, and G. Savar´e. Gradient flows in metric spaces and in the space of probability measures .Lectures in Mathematics ETH Z¨urich. Birkh¨auser Verlag, Basel, second edition, 2008.[2] L. Ambrosio, S. Lisini, and G. Savar´e. Stability of flows associated to gradient vector fields and convergenceof iterated transport maps. manuscripta mathematica , 121(1):1–50, 2006.[3] J.-D. Benamou, G. Carlier, Q. M´erigot, and E. Oudet. Discretization of functionals involving the monge–amp`ere operator. Numerische Mathematik , pages 1–26, 2015.[4] M. Bessemoulin-Chatard and F. Filbet. A finite volume scheme for nonlinear degenerate parabolic equations. SIAM J. Sci. Comput. , 34(5):B559–B583, 2012.[5] A. Blanchet, V. Calvez, and J. A. Carrillo. Convergence of the mass-transport steepest descent scheme forthe sub-critical Patlak-Keller-Segel model. SIAM J. Numer. Anal. , 46(2):691–721, 2008.[6] A. Blanchet, V. Calvez, and J. A. Carrillo. Convergence of the mass-transport steepest descent scheme forthe subcritical Patlak-Keller-Segel model. SIAM J. Numer. Anal. , 46(2):691–721, 2008.[7] C. J. Budd, G. J. Collins, W. Z. Huang, and R. D. Russell. Self-similar numerical solutions of the porous-medium equation using moving mesh methods. R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. ,357(1754):1047–1077, 1999. AGRANGIAN SCHEME FOR NONLINEAR DIFFUSION 31 [8] V. Calvez and T. O. Gallou¨et. Particle approximation of the one dimensional Keller-Segel equation, stabilityand rigidity of the blow-up. Discrete Contin. Dyn. Syst. Ser. A , 36(3):1175–1208, 2015.[9] J. A. Carrillo, A. Chertock, and Y. Huang. A finite-volume method for nonlinear nonlocal equations with agradient flow structure. Communications in Computational Physics , 17(01):233–258, 2015.[10] J. A. Carrillo, Y. Huang, F. S. Patacchini, and G. Wolansky. Numerical study of a particle method forgradient flows. Kinet. Relat. Models , 10(3):613–641, 2017.[11] J. A. Carrillo and S. Lisini. On the asymptotic behavior of the gradient flow of a polyconvex functional.In Nonlinear partial differential equations and hyperbolic wave phenomena , volume 526 of Contemp. Math. ,pages 37–51. Amer. Math. Soc., Providence, RI, 2010.[12] J. A. Carrillo and J. S. Moll. Numerical simulation of diffusive and aggregation phenomena in nonlinearcontinuity equations by evolving diffeomorphisms. SIAM J. Sci. Comput. , 31(6):4305–4329, 2009/10.[13] J. A. Carrillo, F. S. Patacchini, P. Sternberg, and G. Wolansky. Convergence of a particle method for diffusivegradient flows in one dimension. SIAM J. Math. Anal. , 48(6):3708–3741, 2016.[14] J. A. Carrillo, H. Ranetbauer, and M.-T. Wolfram. Numerical simulation of nonlinear continuity equationsby evolving diffeomorphisms. J. Comput. Phys. , 327:186–202, 2016.[15] S. Daneri and G. Savar´e. Eulerian calculus for the displacement convexity in the Wasserstein distance. SIAMJ. Math. Anal. , 40(3):1104–1122, 2008.[16] P. Degond and F.-J. Mustieles. A deterministic approximation of diffusion equations using particles. SIAMJ. Sci. Statist. Comput. , 11(2):293–310, 1990.[17] B. D¨uring, D. Matthes, and J. P. Miliˇsic. A gradient flow scheme for nonlinear fourth order equations. DiscreteContin. Dyn. Syst. Ser. B , 14(3):935–959, 2010.[18] L. Evans, O. Savin, and W. Gangbo. Diffeomorphisms and nonlinear heat flows. SIAM journal on mathe-matical analysis , 37(3):737–751, 2005.[19] L. Gosse and G. Toscani. Identification of asymptotic decay to self-similarity for one-dimensional filtrationequations. SIAM J. Numer. Anal. , 43(6):2590–2606 (electronic), 2006.[20] L. Gosse and G. Toscani. Lagrangian numerical approximations to one-dimensional convolution-diffusionequations. SIAM J. Sci. Comput. , 28(4):1203–1227 (electronic), 2006.[21] R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the Fokker-Planck equation. SIAMJ. Math. Anal. , 29(1):1–17, 1998.[22] O. Junge, D. Matthes, and H. Osberger. A fully discrete variational scheme for solving nonlinear fokker-planckequations in higher space dimensions. arXiv preprint arXiv:1509.07721 , 2015.[23] P.-L. Lions and S. Mas-Gallic. Une m´ethode particulaire d´eterministe pour des ´equations diffusives nonlin´eaires. C. R. Acad. Sci. Paris S´er. I Math. , 332(4):369–376, 2001.[24] S. Mas-Gallic. The diffusion velocity method: a deterministic way of moving the nodes for solving diffusionequations. Transport Theory and Statistical Physics , 31(4-6):595–605, 2002.[25] D. Matthes and H. Osberger. Convergence of a variational Lagrangian scheme for a nonlinear drift diffusionequation. ESAIM Math. Model. Numer. Anal. , 48(3):697–726, 2014.[26] R. J. McCann. A convexity principle for interacting gases. Adv. Math. , 128(1):153–179, 1997.[27] H. Osberger and D. Matthes. Convergence of a variational Lagrangian scheme for a nonlinear drift diffusionequation. ESAIM Math. Model. Numer. Anal. , 48(3):697–726, 2014.[28] H. Osberger and D. Matthes. Convergence of a fully discrete variational scheme for a thin-film equation. Accepted at Radon Ser. Comput. Appl. Math. , 2015.[29] H. Osberger and D. Matthes. A convergent Lagrangian discretization for a nonlinear fourth order equation. Found. Comput. Math. , pages 1–54, 2015.[30] F. Otto. The geometry of dissipative evolution equations: the porous medium equation. Comm. PartialDifferential Equations , 26(1-2):101–174, 2001.[31] G. Peyr´e. Entropic approximation of Wasserstein gradient flows. SIAM J. Imaging Sci. , 8(4):2323–2351,2015.[32] G. Russo. Deterministic diffusion of particles. Comm. Pure Appl. Math. , 43(6):697–733, 1990.[33] J. L. V´azquez. The porous medium equation: mathematical theory . Oxford University Press, 2007.[34] C. Villani. Topics in optimal transportation , volume 58 of Graduate Studies in Mathematics . AmericanMathematical Society, Providence, 2003.[35] M. Westdickenberg and J. Wilkening. Variational particle schemes for the porous medium equation and for thesystem of isentropic euler equations. ESAIM: Mathematical Modelling and Numerical Analysis , 44(1):133–166, 2010.2 JOS´E A. CARRILLO, BERTRAM D¨URING, DANIEL MATTHES, AND DAVID S. MCCORMICK