A discrete Bakry-Emery method and its application to the porous-medium equation

The exponential decay of the relative entropy associated to a fully discrete porous-medium equation in one space dimension is shown by means of a discrete Bakry-Emery approach. The first ingredient of the proof is an abstract discrete Bakry-Emery method, which states conditions on a sequence under which the exponential decay of the discrete entropy follows. The second ingredient is a new nonlinear summation-by-parts formula which is inspired by systematic integration by parts developed by Matthes and the first author. Numerical simulations illustrate the exponential decay of the entropy for various time and space step sizes.


Introduction
The Bakry-Emery method allows one to establish convex Sobolev inequalities and to compute exponential decay rates towards equilibrium for solutions to diffusion equations [2,3]. The key idea of Bakry and Emery is to differentiate a so-called entropy functional twice with respect to time and to relate the second-order derivative to the entropy production. Our aim is to develop a discrete version of this technique, and in this paper, we present a step forward in this direction.
The study of discrete Bakry-Emery methods and related topics is rather recent. Caputo et al. [5] computed exponential decay rates for time-continuous Markov processes, using the Bochner-Bakry-Emery method. Given a stochastic process with density u(t) and the entropy functional H c (u(t)), the core of the Bakry-Emery approach is to find a constant λ > 0 such that the inequality d 2 H c /dt 2 ≥ −λdH c /dt holds for all time. Integrating this inequality, one may show that dH c /dt ≤ −λH c which implies that H c (u(t)) ≤ e −λt H c (u(0)) for all t > 0, i.e., the entropy decays exponentially fast along u(t). The relation between d 2 H c /dt 2 and dH c /dt is achieved in [5] by employing a discrete Bochner-type identity which replaces the Bochner identity of the continuous case. The Bochner-Bakry-Emery method was extended by Fathi and Maas in [13] in the context of Ricci curvature bounds and used by the authors of [19] to derive discrete Beckner inequalities.
Another approach has been suggested by Mielke [21]. He investigated geodesic convexity properties of nonlocal transportation distances on probability spaces such that continuous-time Markov chains can be formulated as gradient flows. Very related results have been obtained independently by Chow et al. [10] and Maas [20]. The geodesic convexity property implies exponential decay rates [1]. Mielke showed that the inequality d 2 H c /dt 2 ≥ −λdH c /dt is equivalent to the positive semi-definiteness of a certain matrix such that matrix algebra can be applied. This idea was extended recently to certain nonlinear Fokker-Planck equations [7].
All these examples involve spatial semi-discretizations of diffusion equations. Temporal semi-discretizations often employ the implicit Euler scheme since it gives entropy dissipation, dH c /dt ≤ 0, under rather general conditions; see, e.g., the implicit Euler finite-volume approximations in [8,14]. Entropy-dissipating higher-order semi-discretizations have been analyzed in [11,17,18]. However, there seem to be no results for fully discrete schemes using the Bakry-Emery approach. In this paper, we make a first step to fill this gap.
In order to understand the mathematical difficulty in fully discrete schemes, consider the abstract Cauchy problem (1) ∂ t u + A(u) = 0, t > 0, u(0) = u 0 , where A : D(A) → X ′ is a (nonlinear) operator defined on its domain D(A) ⊂ X of the Banach space X with dual X ′ . If the dual product A(u), H ′ c (u) is nonnegative, where H ′ c (u) is the (Fréchet) derivative of the entropy and u(t) a solution to (1), then dH c dt = ∂ t u, H ′ c (u) = − A(u), H ′ c (u) ≤ 0, showing entropy dissipation. Next, consider the implicit Euler scheme where u k is an approximation of u(kτ ) and A h is an approximation of A still satisfying A h (u k ), H ′ (u k ) ≥ 0. Here, H(u k ) is the discrete entropy, which is supposed to be convex. Then entropy dissipation is preserved by the scheme since The problem is to estimate the discrete analog of d 2 H c /dt 2 . It turns out that the inequality in (2) is too weak, we need an equation; see Section 2.1 for details. We overcome this difficulty by developing two ideas. The first idea is to identify the elements which are necessary to build an abstract discrete Bakry-Emery method. Unlike in the continuous case, we distinguish between the discrete entropy production P := −τ −1 (H(u k ) − H(u k−1 )) and the Fisher information F := A h (u k ), H ′ (u k ) . We explain this difference in Section 2.
The Bakry-Emery method relies on an estimate of τ −1 (F (u k ) − F (u k−1 )), which approximates d 2 H c /dt 2 . For this estimate, discrete versions of suitable integrations by parts and chain rules are necessary. Our second idea is to "translate" a nonlinear integration-by-parts formula to the discrete case, using the systematic integration by parts method of [16]. This leads to a new inequality for numerical three-point schemes as explained next.
Again, consider first the continuous case. We show in Lemma 7 that for all (A, B) ∈ R c := {(A, B) ∈ R 2 : (2A − B − 1)(A + B − 2) < 0} and all smooth positive functions w, where the constant κ c > 0 depends on A and B; see (16) below. The proof is based on systematic integration by parts [16]. The discrete counterpart is the following inequality: For any 0 < ε ≤ 1, there exists a region R of admissible values (A, B), containing the line where κ = εA; see Lemma 8. Interestingly, the inequality is not true for each term but only with the sum. The admissible set R for (3) is generally smaller than R c ; see Section 3. We conjecture that R = R c for κ = 0.
Inequality (3) is the first nonlinear summation-by-parts formula derived from a systematic method. We believe that this idea will lead to a whole family of new finite-difference inequalities useful in numerical analysis, and we will explore this in a future work.
To overcome the difficulty with the entropy production inequality, we introduce the new variable v = u α and write the porous-medium equation in the form The advantage of this formulation is that the entropy becomes linear in the variable v. thus avoiding inequality (2). We discretize (5) by an implicit Euler finite-difference scheme. Let τ > 0 be the time step, h > 0 the space step, and let v k i = (u k i ) α be an approximation of (h −1 ih (i−1)h u(x, kτ )dx) α , i = 1, . . . , N. The iterative scheme reads as We show in Lemma 10 the existence of solutions to (6) as well as the preservation of nonnegativity. However, the total mass which is the price that we have to pay for the estimation of the entropy production. We discuss this point in Section 5. Our main result reads as follows.
be a nonnegative solution to (6) and set u k i = (v k i ) 1/α . Let 0 < ε < 1. Then there exist a region S ⊂ (0, ∞) 2 , containing the line α − β = 1, and a number U > 0 such that all (α, β) ∈ S with α > 1 and β ≥ 1, it holds that is the discrete (relative) entropy, and C p = h 2 /(4 sin 2 (hπ)) ≥ 1/(4π 2 ) the the discrete Poincaré constant. Moreover, the total mass h N i=1 u k i is increasing in k and converges to U as k → ∞. Remark 2 (Exponential versus algebraic decay). The exponential decay rate depends on the minimum of the solution, which is not surprising. Indeed, because of the degeneracy, we cannot generally expect exponential decay; an example is the Barenblatt solution. Algebraic decay rates for implicit Euler finite-volume schemes have been derived in, e.g., [8]. When the minimum is positive, the equation is no longer degenerate, and exponential decay follows.
Remark 3 (Shannon entropy). Unfortunately, the theorem does not apply to the Shannon entropy h i u i log u i , corresponding to α → 1, since λ → 0 as α → 1. The reason is that for α → 1, the entropy production P cannot be bounded from above by the Fisher information F and so, Assumption A1 of our abstract Bakry-Emery method does not hold; see Section 2.2.
Remark 4 (Discrete gradient flow). Erbar and Maas [12] showed that the gradient flow of the Shannon entropy with respect to a nonlocal transportation measure equals the discrete porous-medium equation in one space dimension. The porous-medium equation in several space dimensions was solved by Benamou et al. [4] by providing a spatial discretization of this equation as a convex optimization problem. In both references, no decay rates have been derived.
The set S is illustrated in Figure 1 for two different values of ε. Numerical computations indicate that S converges (in the set theoretical sense) to the set S c defined in (4) if ε → 0 but for fixed ε > 0, S is strictly contained in S c . The set S c , defined by −1 < α −β < 2, is shown in light blue for comparison; it contains the dark blue region S.
The paper is organized as follows. The abstract Bakry-Emery result is presented in Section 2, and in Section 3, inequality (3) is verified. Theorem 1 is proved in Section 4. Numerical examples are presented in Section 5, and some auxiliary inequalities are recalled in the Appendix.

An abstract Bakry-Emery method
In this section, we present our abstract result. In order to identify the key ingredients of the Bakry-Emery method, we recall the basic ideas for continuous evolution equations.
2.1. The continuous Bakry-Emery method. Let us first consider the abstract Cauchy problem The nonlinear operator A is defined on some domain D(A) of a Banach space X. We do not specify the properties of A nor its domain since they are not needed in the following. As mentioned in the introduction, the idea of the Bakry-Emery method is to differentiate the entropy functional H c : D(A) → [0, ∞) twice with respect to time along solutions to (7). We define the entropy production P c (u(t)) : , the entropy production is nonnegative and the entropy is nonincreasing along solutions to (7). We call F c (u) := A(u), H ′ c (u) the generalized Fisher information since if A(u) = −∆u on T d and H c (u) = T d u(log u − 1)dx, we obtain the Fisher information functional F c (u) = 4 T d |∇ √ u| 2 dx. Clearly, P c (u(t)) = F c (u(t)) along solutions u(t) to (7).
and we conclude exponential decay of t → F c (u(t)) with rate λ c > 0. In particular, lim t→∞ F c (u(t)) = 0. Then, integrating the previous inequality over (t, ∞), it follows that Assuming that also and by Gronwall's lemma, t → H c (u(t)) converges exponentially fast to zero with rate λ c . We see that two assumptions are essential: the functional inequality (8) and the limit (10). On the discrete level, we need to distinguish between the (discrete) entropy production and the (discrete) Fisher information since dH c /dt and A(u), H ′ c (u) may differ on the discrete level. We assume that both functionals can be estimated by each other. Instead of the functional inequality (8) we assume a discrete version of inequality (9). Finally, a discrete version of (10) is required.

A discrete Bakry-Emery method. Let two functions
where v, w ∈ R N and τ > 0. We call H an entropy, F the Fisher information, and P the entropy production. The following result does not need any reference to the solution of a discrete problem. Then where λ = (C m /C M )κ and η = log(1 + τ λ)/(τ λ) ∈ (0, 1).
The discrete decay rate λ is generally smaller than the decay rate κ of the Fisher information, since η < 1 and we may have C m < C M . If the entropy production and the Fisher information coincide, i.e. C m = C M = 1, then λ = κ.
Proof. By Assumption A2, it follows that lim k→∞ F (v k ) = 0. By Assumption A2 again and the second inequality in Assumption A1, we have . Taking the sum from k = ℓ + 1 to k = m > ℓ + 1, we find that ). Passing to the limit m → ∞, observing that lim m→∞ F (v m ) = 0 and, by Assumption A3, , which holds for all ℓ ∈ N. It remains to use the first inequality in Assumption A1 to conclude that

A nonlinear summation-by-parts formula
To apply the abstract Bakry-Emery method to the porous-medium equation, we need to verify the assumptions of Proposition 5. The key condition is Assumption A2. To verify it, we "translate" some integrations by parts to the discrete level. It is convenient to investigate the continuous situation first in order to formulate the discrete formula that is needed to show Assumption A2.
Consider the nonlinear diffusion equation where β > 0, and introduce the (relative) entropy Here, u = T u 0 (x)dx is the constant steady state. (Recall that meas(T) = 1.) Then, for any positive smooth solution to (11), Proof. Integrating by parts, the time derivatives of H c (u(t)) become We wish to estimate the second time derivative. To this end, we set w = u (α+β−1)/2 , . Then the derivatives can be written as In Lemma 7 below we show that there exists κ c > 0 such that and we infer that Furthermore, by the Poincaré inequality applied to w x (see Lemma 11), and it follows that Denoting by P c = −dH c /dt the entropy production, this inequality can be formulated as dP c /dt ≤ −λ c P c . Gronwall's lemma then implies that P c (u(t)) ≤ P c (u 0 )e −λct for t > 0 and in particular lim t→∞ P c (u(t)) = 0. Integrating (13) over (t, s) with t < s and passing to the limit s → ∞, we see that Thus, there exists a sequence t j → ∞ such that w x (t j ) L 2 (T) → 0. Following the arguments of [6, Prop. 1ii], 1 it follows that lim t j →∞ H c (u(t j )) = 0, and since t → H c (u(t)) is nonincreasing, any sequence converges, lim t→∞ H c (u(t)) = 0. We integrate (13) over (t, ∞) and use lim t→∞ (dH c /dt)(u(t)) = lim t→∞ H c (u(t)) = 0: Thus, another application of Gronwall's lemma gives the conclusion.
It remains to prove Lemma 7. Set Proof. The idea of the proof is to employ systematic integration by parts [16]. Since we can formulate (15) as the following problem: Find c ∈ R and κ c > 0 such that for all smooth positive functions w, Calculating the derivatives and setting ξ 1 = w x /w, ξ 2 = w xx /w, this inequality is equivalent to The idea is to interpret the integrand as a polynomial in the variables ξ 1 , ξ 2 and to solve the following polynomial decision problem: Find c ∈ R and κ c > 0 such that for all . This problem can be solved explicitly. Clearly, it must hold that A ≥ κ c > 0. We distinguish two cases: κ c = A and κ c < A.
First let 0 < κ c < A. Then (18) is valid if the discriminant is nonpositive, Choosing the minimizing value the discriminant is nonpositive if and only if Set κ c = εA for 0 < ε < 1. Then the previous inequality is true if and only if We infer that if We state now a discrete version of inequality (15). Lemma 8. Let w 0 , . . . , w N +1 ∈ R satisfy w N = w 0 , w N +1 = w 1 and let 0 < ε ≤ 1. There exists a region R ⊂ R 2 , containing the line A = 1, such that for all (A, B) ∈ R, where κ = εA > 0.
The lemma is trivial as stated since (21) clearly holds for R = {(A, B) : A = 1} with κ = 1. Figure 2 ilustrates the numerical admissible regions for (A, B) for two different values of ε. The admissible region R is smaller than the region R c for the continuous case but it approaches the latter region as κ → 0. The idea of the proof of (21) is to add the following discrete version of the integrationby-parts formula (17), namely The sum vanishes because of the periodic boundary conditions. Here ρ > 0 is a free parameter, and the function M(x, y) is a symmetric mean value, i.e., it satisfies for all x, y, λ ≥ 0. For the numerical simulations below, we will choose ρ = (A + B + 1)/3 such that the mean function does not need to be specified. Then (21) holds if we can show the following inequality for all admissible (A, B) and w i = 0: We verify this inequality pointwise, i.e. setting X = w i+1 /w i and Y = w i−1 /w i , we wish to find c ∈ R, κ > 0 such that for all X, Y > 0, The first term (X A + Y A − 2)(X + Y − 2) becomes negative in certain regions; see Figure  3. It is compensated by the second term (shift term) on the right-hand side of (23) if we choose the constant c according to (19) with κ = κ c as in (16).  (19).
Unfortunately, it seems to be difficult to prove (23) analytically in full generality. Note that polynomial quantifier elimination does not apply if A and B are not integers, since the function T (X, Y ) generally is not a polynomial. Instead, we verify (23) analytically for all (A, B) ∈ R c and all (X, Y ) in some neighborhood of (1, 1). Lemma 9. Let T be given by (23) and let (A, B) ∈ R c , where R c is defined in (14). Then there exists a neighborhood W of (1, 1) such that for all (X, Y ) ∈ W , holds for c as in (19) and with κ c = κ as in (16).
If the step size h > 0 is small enough, we expect that the quotients w i+1 /w i are close to one for all i = 0, . . . , N − 1. This means that (X, Y ) lies in a neighborhood of (1, 1), and the lemma applies.
Proof. We use the local coordinates u = (X + Y − 2)/h 2 and v = (X − Y )/(2h), which correspond to (central) second-order and first-order derivatives. Then X = 1 + hv + h 2 u/2 and Y = 1 − hv + h 2 u/2. We develop T as a function of h at h = 0. For this, we observe that M(1, 1) = 1 and M X (1, 1) = M Y (1, 1) = 1/2. Indeed, we infer from the properties (22) that which implies that M X (1, 1) = 1/2. Calculating the Taylor series of the terms in T with respect to h at h = 0 leads to In particular, as expected, the explicit choices of both ρ and M(x, y) do not change the behavior of the shift term locally around the equilibrium w i−1 = w i = w i+1 or h = 0. Moreover, min{1, X A+B−1 , Y A+B−1 } = 1 + O(h) and (X + Y − 2) 2 = u 2 h 4 . Combining these expressions gives The polynomial is the same as in (18). The proof of Lemma 7 shows that it is nonnegative for all (A, B) ∈ R c with c as in (19) and κ c as in (16). We deduce that T (X, Y ) ≥ 0 holds for all (A, B) ∈ R c if h ∈ R is sufficiently small. This proves the lemma.

Discrete porous-medium equation
We apply the abstract Bakry-Emery method to a finite-difference approximation of the porous-medium equation, i.e., we choose A(u) = −(u β ) xx on T for suitable functions u. Let τ > 0 be the time step and h > 0 the space step. A natural scheme would be for all i = 1, . . . , N, k ∈ N, and u k i is the constant steady state. We choose α > 1 and β > 1.
Unfortunately, the abstract Bakry-Emery method cannot be applied to this scheme. The problem is the second inequality in Assumption A1. Indeed, using the inequality y α − z α ≥ αz α−1 (y − z) for all y, z ≥ 0, which follows from the convexity of z → z α for α > 1, inserting the numerical scheme and then summing by parts, we find that This expression cannot be estimated further; it may even have the wrong sign. We need a scheme that avoids the use of the inequality y α − z α ≥ αz α−1 (y − z). We stress the fact that this problem does not occur in the semi-discrete scheme and this expression is nonpositive (since α > 1). Our idea is to make the entropy production linear in its argument. For this, we introduce the new variable v k i = (u k i ) α . In the (continuous) variable v = u α , the evolution equation The discrete entropy and Fisher information become where V > 0 has to be determined and γ = (α + β − 1)/(2α). The entropy production can be estimated, using summation by parts, as According to Lemma 13, the entropy production can be estimated from below and above in terms of the Fisher information.
After this motivation, we prove the existence of solutions to (24). Proof. We give only a sketch of the proof since the existence of solutions follows from a standard fixed-point theorem. We only provide the a priori estimates needed for this argument. First multiply (24) by (v k i ) − = min{v k i , 0} and sum over i = 1, . . . , N.

By summation by parts, this becomes
and this is the desired a priori estimate.
Next, we turn to the proof of our main result.
The idea of the proof is to expand G(τ ) around zero: for some ξ ∈ (0, τ ). We show that the right-hand side can be bounded from above by −τ KF (v) for some K > 0, which verifies Assumption A2. This idea has been first employed in [8]. Clearly, we have G(0) = 0. The first derivatives of G equal First, we claim that G ′′ (τ ) ≤ 0 for any τ > 0. Indeed, we replace v i − τ a i by v i and obtain It holds that c 1 ≥ 0, since this inequality is equivalent to (2γ − 1)v γ i+1 ≥ (γ − 1)v γ i , and this is true for 1/2 ≤ γ ≤ 1 (which is equivalent to β ≥ 1 and α − β ≥ −1). Moreover, the discriminant 4c 1 c 3 − c 2 2 ≥ 0 is equivalent to , which also holds true for 1/2 ≤ γ ≤ 1. This shows that G ′′ (τ ) ≤ 0 and consequently, It remains to compute G ′ (0). Inserting the definition of a i , we find that By the discrete Poincaré-Wirtinger inequality (Lemma 12), applied with and hence, with This shows Assumption A2 of Proposition 5 and, in particular, after applying Gronwall's lemma, lim k→∞ F (v k ) = 0. It remains to prove that Assumption A3, i.e. lim k→∞ H(v k ) = 0, holds. We know that so, for any fixed i = 1, . . . , N, (v k i ) is bounded. Therefore, there exists a sequence k j → ∞ such that v k j i → y i for some y i ≥ 0. By the discrete Poincaré-Wirtinger inequality (Lemma 12), applied to Since lim k→∞ F (v k ) = 0, we deduce that (v k j i ) and V k j have the same limit, say y := y i . Set U := y 1/α . This defines the entropy According to Proposition 5, the discrete entropy converges exponentially with decay rate Next, we claim that the total mass h N i=1 u k i is nondecreasing in k. Indeed, by the concavity of z → z 1/α (recall that α > 1), we have y 1/α − z 1/α ≥ (1/α)y (1−α)/α (y − z) for all y, z ≥ 0 and hence, Inserting scheme (6), we find that since v k i satisfies periodic boundary conditions. This shows the claim. The monotonicity of the total mass and the convergence property , and the convergence is monotone. This finishes the proof.
The profile will hit the boundary of Ω = (0, 1) at time t end . For the fast diffusion case β = 0.5, we take x 0 = 0.5, t 0 = 10 −2 , and C = t (β−1)/(β+1) 0 such that the maximum of the initial profile equals 1. Figure 4 illustrates the evolution of the total mass for α = 2, β = 0.5 (left) and α = 3, β = 4 (right). As predicted in Theorem 1, the total mass is indeed increasing in time. The mass defect scales well with both the time step τ and the grid size h, where the influence of τ is more prevalent.  The time decay of the (relative) entropy H is shown in Figure 5 for various space and time steps. We observe that the decay is indeed exponential. Here, the steady state u ∞ (which is needed to define the relative entropy) is given by u ∞ = h N i=0 u kmax i , where k max is the final time step. This choice clearly depends on the scheme since the mass is not conserved. The relative entropy converges exponentially even when (α, β) is chosen outside of the admissible region; see Figure 6.
where C p = h 2 /(4 sin 2 (hπ)) ≥ 1/(4π 2 ). This constant is sharp.  Proof. The second inequality is proven in [8,Lemma A.3]. For the proof of the first inequality, we divide it by y a+b and set z = x/y. Then the inequality is equivalent to which after expansion can be equivalently written as (z a/2 −z b/2 ) 2 ≥ 0, and this is true.