Monge-Kantorovich distance for PDEs: the coupling method

TL;DR

This paper reviews the use of the coupling method to analyze the decay of Monge-Kantorovich distances between solutions of various PDEs, highlighting its effectiveness and limitations across different equations.

Contribution

It introduces a unified coupling approach for PDEs to study solution distances and discusses alternative dual methods, expanding the analytical toolkit for PDE analysis.

Findings

Coupling method effectively shows decay of solution distances for several PDEs.

The method's limitations are identified for certain nonlinear PDEs like porous media.

Alternative dual Monge-Kantorovich approaches are proposed for challenging cases.

Abstract

We informally review a few PDEs for which the Monge-Kantorovich distance between pairs of solutions, possibly with some judicious cost function, decays: heat equation, Fokker-Planck equation, heat equation with varying coefficients, fractional heat equation with varying coefficients, homogeneous Boltzmann equation for Maxwell molecules, and some nonlinear integro-differential equations arising in neurosciences. We always use the same method, that consists in building a coupling between two solutions. This amounts to solve a well-chosen PDE posed on the Euclidian square of the physical space, i.e. doubling the variables. Finally, although the above method fails, we recall a simple idea to treat the case of the porous media equation. We also introduce another method based on the dual Monge-Kantorovich problem.

Equations

$$\left\{\begin{array}[]{l}d_{\varrho}(u_{1},u_{2})=\displaystyle\inf_{v\in K(u_{1},u_{2})}\int\hskip-8.0pt\int\varrho(x,y)v(x,y)dx\,dy,\\[10.0pt] K(u_{1},u_{2})=\{v:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}_{+}\quad\hbox{such that}\quad\int v(x,y)dy=u_{1}(x),\;\int v(x,y)dx=u_{2}(y)\}.\end{array}\right.$$

$${\varrho }_p(x,y)=\frac{|x-y{|}^p}{p},$$

$${\partial }_{t}u-\Delta u=0,x\in {\mathbb{R}}^d,t\ge 0.$$

$$d_{\varrho }(u_1(t),u_2(t))\le d_{\varrho }(u_1^0,u_2^0).$$

$$\frac{\partial v}{\partial t}-{\Delta }_{x}v-{\Delta }_{y}v-2{\nabla }_x\cdot {\nabla }_{y}v=0,x,y\in {\mathbb{R}}^d,t\ge 0$$

$$v_1(x,t)=\int v(x,y,t)dy,v_2(y,t)=\int v(x,y,t)dx$$

$$\left\{\begin{array}[]{l}\frac{\partial v_{1}(x,t)}{\partial t}-\Delta_{x}v_{1}(x,t)=0,\quad x\in\mathbb{R}^{d},\;t\geq 0,\\[5.0pt] v_{1}(x)=u^{0}_{1}(x),\quad x\in\mathbb{R}^{d}\end{array}\right.$$

$$d_{\varrho }(u_1(t),u_2(t))\le \int \int \varrho (x,y)v(x,y,t)dxdy=\int \int r(|x-y|)v(x,y,t)dxdy.$$

$$\frac{d}{dt}\int\hskip-8.0pt\int r(|x-y|)v(x,y,t)dxdy=\!\int\hskip-8.0pt\int v(x,y,t)\Big{(}\Delta_{x}[r(|x-y|)]+\Delta_{y}[r(|x-y|)]+2\nabla_{x}\nabla_{y}[r(|x-y|)]\Big{)}dxdy=\!0.$$

$$d_{\varrho }(u_1(t),u_2(t))\le \int \int r(|x-y|)v^0(x,y)dxdy=\int \int \varrho (x,y)v^0(x,y)dxdy$$

$$(V(x,t)-V(y,t))\cdot (x-y)\le \alpha |x-y{|}^2,x,y\in {\mathbb{R}}^d,t\ge 0.$$

$${\partial }_{t}u-\Delta u+\mathrm{div}(V(x,t)u)=0,x\in {\mathbb{R}}^d,t\ge 0.$$

$$d_p(u_1(t),u_2(t))\le d_p(u_1^0,u_2^0)\exp (\alpha pt).$$

$${\partial }_{t}v-{\Delta }_{x}v-{\Delta }_{y}v-2{\nabla }_x\cdot {\nabla }_{y}v+{\mathrm{div}}_x(V(x,t)v)+{\mathrm{div}}_y(V(y,t)v)=0,x,y\in {\mathbb{R}}^d,t\ge 0$$

$$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int\hskip-8.0pt\int\frac{|x-y|^{p}}{p}v(x,y,t)dxdy&=0+\displaystyle\int\hskip-8.0pt\int v(x,y,t)|x-y|^{p-2}(x-y)\cdot(V(x,t)-V(y,t))dxdy\\[10.0pt] &\leq\displaystyle\alpha\int\hskip-8.0pt\int|x-y|^{p}v(x,y,t)dxdy\end{array}$$

$$d_{p}(u_{1}(t),u_{2}(t))\leq\int\hskip-8.0pt\int\frac{|x-y|^{p}}{p}v(x,y,t)dxdy\leq\Big{(}\int\hskip-8.0pt\int\frac{|x-y|^{p}}{p}v^{0}(x,y)dxdy\Big{)}e^{\alpha pt}$$

$${\partial }_{t}u+\mathrm{div}[V(x,I(t))u]=0,x\in {\mathbb{R}}^d,t\ge 0,$$

$$I(t)={\int }_{{\mathbb{R}}^d}\psi (x)u(x,t)dx.$$

$$(x-y)\cdot\big{(}V(x,I)-V(y,I)\big{)}\leq-\alpha|x-y|^{2},\qquad\forall x,y\in\mathbb{R}^{d},\;I\geq 0.$$

$$d_2(u(t),{\delta }_{⟨x(t)⟩})=\int \frac{|x-⟨x(t)⟩{|}^2}{2}u(x,t)dx\le e^{-2\alpha t}\int \frac{|x-⟨x(0)⟩{|}^2}{2}{u}^0(x,t)dx=e^{-2\alpha t}{d}_2(u(0),{\delta }_{⟨x(0)⟩}).$$

$$\beta =\Vert D_{I}V{\Vert }_{\infty }\Vert D\psi {\Vert }_{\infty }<\alpha$$

$$d_2(u(t),{\delta }_{X(t)})=\int \frac{|x-X(t){|}^2}{2}u(x,t)dx\le e^{2(\beta -\alpha )t}\int \frac{|x-X^0{|}^2}{2}{u}^0(x)dx=d_2(u(0),{\delta }_{X^0}).$$

$${\partial }_{t}v+{\mathrm{div}}_x[V(x,I_1(t))v]+{\mathrm{div}}_y[V(y,I_2(t))v]=0.$$

$$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int\hskip-8.0pt\int\frac{|x-y|^{2}}{2}&u_{1}(x,t)u_{2}(y,t)dxdy=\displaystyle\int\hskip-8.0pt\int(x-y)\cdot\big{(}V(x,I_{1}(t))-V(y,I_{2}(t))\big{)}u_{1}(x,t)u_{2}(y,t)dxdy\\[10.0pt] =&\displaystyle\int\hskip-8.0pt\int(x-y)\cdot\big{(}V(x,I_{1}(t))-V(y,I_{1}(t))\big{)}u_{1}(x,t)u_{2}(y,t)dxdy\\[10.0pt] &+\displaystyle\int\hskip-8.0pt\int(x-y)\cdot\big{(}V(y,I_{1}(t)))-V(y,I_{2}(t))\big{)}u_{1}(x,t)u_{2}(y,t)dxdy\\[10.0pt] \leq&-\displaystyle\alpha\int\hskip-8.0pt\int\!|x-y|^{2}u_{1}(x,t)u_{2}(y,t)dxdy\\ &+\displaystyle\|D_{I}V\|_{\infty}|I_{1}(t)-I_{2}(t)|\Big{(}\int\hskip-8.0pt\int|x-y|^{2}u_{1}(x,t)u_{2}(y,t)dxdy\Big{)}^{1/2}.\end{array}$$

$$\int \int \frac{|x-y{|}^2}{2}u(x,t)u(y,t)dxdy\le e^{-2\alpha t}\int \int \frac{|x-y{|}^2}{2}{u}^0(x,t)u^0(y,t)dxdy.$$

$$\int |x-⟨x(t)⟩{|}^{2}u(x,t)dx\le e^{-2\alpha t}\int |x-⟨x(0)⟩{|}^2{u}^0(x,t)dx.$$

$$I_1(t)-I_2(t)=\int \int [\psi (x)-\psi (y)]u_1(x,t)u_2(y,t)dxdy,$$

$$\big{|}I_{1}(t)-I_{2}(t)\big{|}\leq\|D\psi\|_{\infty}\int\hskip-8.0pt\int|x-y|u_{1}(x,t)u_{2}(y,t)dxdy\leq\|D\psi\|_{\infty}\Big{(}\int\hskip-8.0pt\int|x-y|^{2}u_{1}(x,t)u_{2}(y,t)dxdy\Big{)}^{1/2}.$$

$$\frac{d}{dt}\int \int \frac{|x-y{|}^2}{2}{u}_1(x,t)u_2(y,t)dxdy\le (\beta -\alpha )\int \int |x-y{|}^2{u}_1(x,t)u_2(y,t)dxdy.$$

$$\frac{\partial u}{\partial t}-\frac{{\partial }^2}{\partial x^2}[a(x)u]=0,x\in \mathbb{R},t\ge 0.$$

Full text

2010 Mathematics Subject Classification: 35A05; 35K55; 60J99; 28A33

Keywords and phrases: Monge-Kantorovich distance; Coupling; Fokker-Planck equation; Fractional Laplacian; Homogeneous Boltzmann equation; Porous media equation; Integro-differential equations.

Introduction

It is usual to study the well-posedness, stability and large-time behavior of stochastic processes (e.g. solutions to Stochastic Differential Equations) by using coupling methods: we consider two such processes, with different initial conditions, driven by suitably correlated randomness, and we measure the $\varrho$-Monge-Kantorovich distance $d_{\varrho}$ between their distributions.

We work in $\mathbb{R}^{d}$ and we always assume that the cost function $\varrho:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}$ satisfies $\varrho(x,x)=0$ and $\varrho(x,y)=\varrho(y,x)>0$ for $x\neq y$. We recall that for two probability densities $u_{1},u_{2}$ on $\mathbb{R}^{d}$,

[TABLE]

Observe that $d_{\varrho}$ is not always really a distance because it does not automatically satisfy the triangular inequality. However, this is the case when, for some $p\geq 1$,

[TABLE]

and we put $d_{p}=d_{\varrho_{p}}$.

The probabilistic coupling method can easily be written in terms of the Kolmogorov equation of the coupled process. The goal of the present survey paper is to describe, in an informal way, this method, using only arguments based on Partial Differential Equations (PDEs in short). The difficulties and novelties rely on the choice of the cost function and on the choice of coupling between two solutions by solving a well-chosen PDE posed on the Euclidian square of the physical space, $\mathbb{R}^{2d}$ in general. Each time, we try to emphasize the main technical difficulties that would allow one to justify the computations.

For example, considering the Brownian motion leads to the heat equation. We first give a simple proof that the heat equation is non-expansive (weak contraction) for any smooth cost function of the form $\varrho(x,y)=r(|x-y|)$. This is standard but the PDE literature seems to ignore this simple approach. The method can be extended to various cases. The Fokker-Planck equation is the simplest extension. The case of the heat equation with variable coefficients, of the form $\partial_{t}u-\Delta(a(x)u)=0$, is more involved: in one dimension, the distance $d_{1}$ plays a central role and is always non-expansive (under technical conditions); we illustrate the general structure in higher dimension and show that if the cost function $\varrho$ satisfies some elliptic PDE, which does seem to enter a class with generic existence results, then $d_{\varrho}$ is non-expansive along solutions. The method also applies to some jump processes: fractional heat equation with variable coefficients in dimension one, scattering equations, kinetic scattering equations, Boltzmann equation for Maxwell molecules.

For the porous media equation, the situation is more intricate and the above method does not seem to apply. However, we recall from [5] another, somehow related and rather simple, path to treat this equation.

Concerning piecewise deterministic jump processes and (inhomogeneous) kinetic scattering equations, we present a new result showing that the 1-Monge-Kantorovich distance is non-expansive.

Finally, concerning jump processes, in particular those related to the discretized heat equation, we present another approach, based on the dual formulation of the Monge-Kantorovich distance.

Recently , the topic of Monge-Kantorovich distance has developed quickly for PDEs and integro-differential equations (IDEs) after new understanding of optimal transportation and of the Brenier-Kantorovich map by [9, 2]. There are several approaches to use the Monge-Kantorovich distance in PDEs. A geometrical approach based on gradient flow structures has been introduced in [24] and extended in [12, 7], in particular for the porous media equation, for interacting particle systems and for granular flows. Also, many results on PDEs have been derived from the splitting algorithm named JKO after [20]. See the book [28] for a complete presentation of these results. Let us also mention that the special structure associated with dimension 1 has been used to prove strict contraction for the porous media equation [11] for the distance $d_{2}$, and to treat other equations as scalar conservation laws [3] or the Keller-Segel system [10]. Methods based on optimal transportation have also been recently used to treat singular congestion (incompressible) equations arising in crowd modeling, see for instance [8, 23, 14].

Most of the recent papers using the Monge-Kantorovich distance for PDEs have been using the gradient flow structure which is closely related to a variational formulation of the fluxes. Here, with several examples of conservative equations, which do not necessarily have a gradient flow structure, we control the Monge-Kantorovich distance using the coupling method. We often borrow our examples from the stochastic processes which represent the PDEs thanks to their Kolmogorov equation. The cases of variable coefficients are particularly interesting because they often require some special choice of the cost function.

We organize our examples as follows. We begin with three simple examples: heat equation, Fokker-Planck equation, and a class of nonlinear transport equations. We show directly that the Monge-Kantorovich distances are non-expansive along these equations. Then we turn, in Section 2, to the heat equation with variable coefficients. In Section 3 we consider some IDEs: scattering equations, including kinetic scattering and inhomogeneous fractional heat equation. Zero-th order terms, describing absorption and re-emission, as they appear in models of neural networks, can also be treated by adapting the method; this is explained in Section 4. The famous Tanaka theorem for the homogeneous Boltzmann equation can be included in our framework and this is done in Section 5. We treat the porous media equation in Section 6. Finally, we exemplify in Section 7 how the same results can be proved using the dual formulation of the Monge-Kantorovich distance.

1 Heat, Fokker-Planck and transport equations

In order to explain the coupling method in a very simple, but still relevant, framework, we begin with the heat equation. Then we turn to drift and transport terms.

1.1 Heat equation

Here is the well-known result, see e.g. [28], we want to quickly recall.

Theorem 1

Consider any increasing function $r:[0,\infty)\mapsto[0,\infty)$ such that the cost function $\varrho:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}_{+}$ defined by $\varrho(x,y)=r(|x-y|)$ is of class $C^{2}$. Consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}^{d}$, and the corresponding solutions $u_{1},u_{2}$ to the heat equation

[TABLE]

For any $t\geq 0$, one has

[TABLE]

Proof. We consider an initial density $v^{0}:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}_{+}$ with marginals $u^{0}_{1}$ and $u^{0}_{2}$, that is such that $v^{0}\in K(u_{1}^{0},u_{2}^{0})$. We next consider the solution $v(x,y,t)$ to the degenerate heat equation

[TABLE]

starting from $v^{0}$. Clearly, it holds that $v(x,y,t)\geq 0$, because of the non-negativity of the operator in (3), which can be written in the variables $(x+y,x-y)$ as $-\Delta_{x+y}$.

We then define the marginals

[TABLE]

and show that $v_{1}=u_{1}$ and $v_{2}=u_{2}$: for instance, integrating (3) with respect to $y$, one finds

[TABLE]

and uniqueness of the solution of the heat equation gives us $v_{1}=u_{1}$.

Recalling (1), we conclude that

[TABLE]

Finally, we may also compute, using (3) and integrating by parts,

[TABLE]

Therefore, for any initial data $v^{0}\in K(u_{1}^{0},u_{2}^{0})$,

[TABLE]

and minimizing among such $v^{0}$ completes the proof. $\square$

The only technical question is to justify the integration by parts, which is immediate if we assume enough moments initially (otherwise there is nothing to prove), at least when we restrict ourselves to power cost functions $\varrho_{p}(x,y)=|x-y|^{p}/p$ with $p\geq 2$. Notice that the well-posedness for (3) follows from the observation that we actually deal with $-\Delta_{x+y}$. It is also possible, under some conditions, to treat the case of some non smooth cost functions, e.g. $\varrho_{p}$ for some $p\in[1,2)$: this issue is discussed in Section 2.

1.2 Fokker-Planck equation

The coupling method can be extended to the Fokker-Planck equation, see [6] for some more elaborate consequences. The result can be stated as follows

Theorem 2

Consider some function $V:\mathbb{R}^{d}\times\mathbb{R}_{+}\mapsto\mathbb{R}^{d}$ such that, for some $\alpha\in\mathbb{R}$,

[TABLE]

Consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}^{d}$ and the corresponding solutions $u_{1},u_{2}$ to the Fokker-Planck equation

[TABLE]

For any $t\geq 0$, any $p\geq 1$, one has

[TABLE]

This inequality is well-known, see for example [28] §9.1.5 and the references therein. One can also find relations to several deep and recent functional analysis tools. This goes far beyond our present purpose.

Proof. This is the same proof as for the heat equation, with longer expressions. We consider any $v^{0}$, with marginals $u^{0}_{1}$ and $u^{0}_{2}$, and the solution $v$ to the equation

[TABLE]

starting from $v^{0}$. One easily checks that $v(x,y,t)\geq 0$ and, integrating (6) with respect to $y$, that $v_{1}(x,t):=\int v(x,y,t)dy$ solves (5) and starts from $u^{0}_{1}$, whence $v_{1}=u_{1}$. The second marginal is treated similarly, and we conclude that $d_{p}(u_{1}(t),u_{2}(t))\leq p^{-1}\int\hskip-4.0pt\int|x-y|^{p}v(x,y,t)dxdy$. Finally, using the same computation as for the heat equation, with some additional terms, we see that

[TABLE]

by assumption (4). The result follows using the Gronwall lemma

[TABLE]

and minimizing in $v^{0}$. $\square$

1.3 A nonlinear transport equation

We next consider a fully deterministic problem which arises in several types of modeling, such as polymers, cell division, neuron networks, etc :

[TABLE]

where the nonlinearity stems from the quantity $I(t)$ defined, with a given weight $\psi:\mathbb{R}^{d}\mapsto\mathbb{R}$,

[TABLE]

We again complement this equation with an initial condition $u^{0}\geq 0$ with mass $\int u^{0}=1$.

Theorem 3

Assume that $V:\mathbb{R}^{d}\times\mathbb{R}\mapsto\mathbb{R}^{d}$ and $\psi:\mathbb{R}^{d}\mapsto\mathbb{R}$ are of class $C^{1}$, and that for some $\alpha>0$

[TABLE]

Setting $\langle x(t)\rangle=\int xu(x,t)dx$, we have

[TABLE]

Assume additionally that

[TABLE]

and fix any initial point $X^{0}\in\mathbb{R}^{d}$. Consider the solution $X$ to $X^{\prime}(t)=-V(X(t),\psi(X(t))$ starting from $X^{0}$. For all $t\geq 0$, one has

[TABLE]

It holds that $(\delta_{X(t)})_{t\geq 0}$ solves (7) in a weak sense. A more general result, involving any pair of solutions, can be found in [28].

Proof. We consider two solutions $u_{1}$ and $u_{2}$ to (7), and denote by $I_{1}(t)$ and $I_{2}(t)$ the corresponding functions, see (8). As we are interested in the case where one of the two solutions is a Dirac mass (for each $t\geq 0$, $u_{2}(t)=\delta_{X(t)}$), we can consider the trivial coupling $v(x,y,t)=u_{1}(x,t)u_{2}(y,t)$, which of course has the correct marginals, and satisfies

[TABLE]

Therefore, we may compute

[TABLE]

We first apply this in the case of single solution $u:=u_{1}=u_{2}$, whence $I_{1}=I_{2}$, and we directly conclude by Gronwall’s lemma that

[TABLE]

This classically rewrites as

[TABLE]

as desired. Next, when considering two solutions, we notice that

[TABLE]

whence

[TABLE]

Therefore,

[TABLE]

Applying this to the case where $u_{2}(t)=\delta_{X(t)}$ concludes the proof. $\square$

For consistency with the other presentations in this section, we have written this result for an $L^{1}$ density $v$ with a finite second moment, but the extension to a probability measure is immediate.

2 Heat equation with variable coefficients

We consider the heat equation with variable coefficient. This is much more intricate than the previous examples. In 1 dimension, we use the $d_{1}$ distance and recover a result implicitly included in [22]. In higher dimension, we indicate a general way to construct cost functions. This leads to a poorly explored degenerate elliptic PDE, see however [25] and the references therein.

2.1 One-dimensional case

We consider some $a:\mathbb{R}\mapsto\mathbb{R}_{+}$ and the following heat equation.

[TABLE]

Theorem 4

Assume that $d=1$ and that $a=\sigma^{2}$ for some $\sigma\in C^{1/2}(\mathbb{R})$. Consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}$ and the corresponding solutions $u_{1},u_{2}$ to (10). For all $t\geq 0$, one has

[TABLE]

Proof. We give a proof for $\sigma\in C^{\alpha}(\mathbb{R})$, with $\alpha>1/2$, the remark below explains how to treat $\alpha=1/2$. We consider any probability density $v^{0}(x,y)$ with marginals $u_{1}^{0}$ and $u_{2}^{0}$ and consider the coupling equation

[TABLE]

starting from $v^{0}$. This equation preserves non-negativity. A simple way to see this is the following computation: mutliplying (11) by $-v_{-}$, integrating on $\mathbb{R}^{2}$ and using some integrations by parts, one can check that

[TABLE]

Since $\int\hskip-4.0pt\int v_{-}^{2}(x,y,0)dx\,dy=0$, the result follows from the Gronwall lemma if $\sigma$ is smooth. Otherwise, one can work by approximation.

Integrating (11) with respect to $y$, we see that $v_{1}(x,t):=\int v(x,y,t)dy$ solves (10) and starts from $u^{0}_{1}$, whence $v_{1}=u_{1}$. The second marginal is treated similarly, and we conclude that $d_{1}(u_{1}(t),u_{2}(t))\leq\int\hskip-4.0pt\int|x-y|v(x,y,t)dxdy$. Because of its singularity, we need to regularize the absolute value as a $W^{2,\infty}$ function and define

[TABLE]

Using the Hölder constant $C_{\sigma}$ of $\sigma(\cdot)$, we see that

[TABLE]

because $v(t)$ is a probability measure. Since now $2\alpha-1>0$, we may let $\varepsilon\to 0$ and we find that

[TABLE]

We conclude, as usual, by minimizing in $v^{0}$. $\square$

Remark 5

The condition $\sigma\in C^{1/2}(\mathbb{R})$ is enough. To treat this exponent, a better construction of the regularization is required, using the so-called Yamada function:

[TABLE]

There are other technical issues here. For example, the well-posedness of (11), which is necessary to identify the marginals of the solution $v$ to the coupling equation, is not so easy. A possible direction is to use results established in [16], in the spirit of [15].

2.2 A general construction of the weight

In order to unravel the algebraic structure behind the choice of the weight $\varrho$, we now consider the general case of dimension $d$. We assume that $a:\mathbb{R}^{d}\mapsto{\mathcal{M}}_{d\times d}(\mathbb{R})$ is everywhere symmetric and nonnegative, of the form

[TABLE]

for some $\sigma:\mathbb{R}^{d}\mapsto{\mathcal{M}}_{d\times K}(\mathbb{R})$, and we consider the heat equation

[TABLE]

completed with an initial probability density $u^{0}$ on $\mathbb{R}^{d}$.

Proposition 6

Assume that $\sigma$ is regular enough and consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}^{d}$ and the corresponding solutions $u_{1},u_{2}$ to (13). For all $t\geq 0$, one has

[TABLE]

for any smooth cost $\varrho:\mathbb{R}^{d}\mapsto\mathbb{R}_{+}$ satisfying

[TABLE]

When $a$ is constant, we recover that any $C^{2}$ cost function of the form $\varrho(x,y)=r(|x-y|)$ works. In dimension $1$, $\varrho(x,y)=|x-y|$ is indeed a (weak) solution to (14). We do not know of a theory to solve (14), in dimension $d\geq 2$, for a general coefficient $a$, so that we do not know if this result is useful. Notice that equation (14) should be completed by the boundary value $\varrho(x,x)=0$ with some growth condition to mimic $|x-y|^{p}$.

Proof. We consider any probability density $v^{0}(x,y)$ with marginals $u_{1}^{0}$ and $u_{2}^{0}$ and consider the coupling equation

[TABLE]

starting from $v^{0}$. We show as usual that $\int v(x,y,t)dy=u_{1}(x,t)$ and that $\int v(x,y,t)dx=u_{2}(y,t)$. Moreover, we have $v(x,y,t)\geq 0$: we multiply the coupling equation by $-v_{-}$ and integrate, finding

[TABLE]

with

[TABLE]

which can also be written

[TABLE]

The other term is

[TABLE]

which can also be written after integration by parts

[TABLE]

Assuming that the entries $\sigma_{ik}$ are bounded with two bounded derivatives, we conclude by Gronwall’s lemma that $v_{-}\equiv 0$, since we initially have $\int\hskip-4.0pt\int v_{-}^{2}(x,y,0)dx\,dy=0$.

Recalling (1), we conclude that $d_{\varrho}(u_{1}(t),u_{2}(t))\leq\int\hskip-4.0pt\int\varrho(x,y)v(x,y,t)dxdy$. Since finally

[TABLE]

by assumption, we conclude as usual. $\square$

We leave open the question to formalize this approach rigorously, in particular for degenerate coefficients $\sigma$, and to build other examples where one can prove the existence of a weight $\varrho$.

3 Scattering and integral kernels

We now turn to equations that describe the probability law of various jump processes. These are well-known results except the case of kinetic scattering in Subsection 3.2 which seems to be new.

3.1 Simple scattering

For $x\in\mathbb{R}^{d}$, we parameterize the pre-jump location $X=\Phi(x,h)$ by $h\in\mathbb{R}^{d}$, distributed according to a bounded measure $\mu$. We assume that for all fixed $h\in\mathbb{R}^{d}$,

[TABLE]

and we use the notation $X\mapsto x=\Phi^{-1}(X,h)$ for the inverse in $x$ (with $h$ fixed).

We consider the scattering problem

[TABLE]

with initial condition $u^{0}$, a probability density on $\mathbb{R}^{d}$. Actually, this equation is to be understood in the weak sense: integrating the right hand side against a test function $\varphi(x)$, we see that

[TABLE]

which shows that the determinant $\det(D_{x}\Phi(x,h))$ is only used informally. We briefly prove the following result, which is classical, see for instance [1].

Theorem 7

Assume (15), fix $p\in[1,\infty)$ and suppose theres is $\delta\in\mathbb{R}$ such that for all $X,Y\in\mathbb{R}^{d}$,

[TABLE]

Consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}^{d}$ and the corresponding solutions $u_{1},u_{2}$ to (16). For all $t\geq 0$, one has

[TABLE]

The homogeneous scattering corresponds to $\Phi(x,h)=x+h$ and obviously fulfills the above assumptions.

Proof. For a probability density $v^{0}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $u^{0}_{1}$ and $u^{0}_{2}$, we consider the solution $v$ of the coupled equation built in such a way that the jumps parameter $h$ is common to the two variables. Namely, we choose

[TABLE]

starting from $v^{0}$. We clearly have $v\geq 0$, and integrating in $y$ and using the change of variable $y\mapsto\Phi(y,h)$, we find that $v_{1}(x,t)=\int v(x,y,t)dy$ satisfies (16). Since it starts from $u^{0}_{1}$, we conclude that $v_{1}=u_{1}$. The second marginal is treated similarly, and we conclude as usual that $d_{p}(u_{1}(t),u_{2}(t))\leq p^{-1}\int\hskip-4.0pt\int|x-y|^{p}v(x,y,t)dxdy$. Next, we compute, using (18):

[TABLE]

We used the changes of variables $X=\Phi(x,h)$ and $Y=\Phi(y,h)$ (with $h$ fixed). Recalling (17), we conclude that

[TABLE]

Using the Gronwall lemma, we thus find that

[TABLE]

and we conclude as usual, minimizing in $v^{0}$. $\square$

The most general scattering equation reads

[TABLE]

and equation (16) corresponds to the homogeneous cases when $\int\pi(x_{*},x)dx_{*}=1$, and the above method can easily be adapted. For the inhomogeneous case, see Section 4.

3.2 Kinetic scattering

We next consider some kinetic scattering models, that means we work in the phase space. We consider some finite measure $\mu$ on $\mathbb{R}^{d}$, some application $V:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ such that, for all $h\in\mathbb{R}^{d}$,

[TABLE]

and the kinetic scattering equation

[TABLE]

completed with an initial data $f^{0}(x,v)\geq 0$ with $\int f^{0}dxdv=1$.

Theorem 8

Assume (20). Set $K=\mu(\mathbb{R}^{d})$ and suppose that for some $L\in\mathbb{R}_{+}$, for all $v,w\in\mathbb{R}^{d}$,

[TABLE]

Suppose that $K\geq KL+1$. Consider two probability densities $f^{0}_{1},f^{0}_{2}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ and the corresponding solutions $f_{1},f_{2}$ to (16). It holds that for all $t\geq 0$, (here $d_{1}$ is associated to the cost function $\varrho((x,v),(y,w))=|x-y|+|v-w|$)

[TABLE]

Proof. As usual, we consider any probability density $F^{0}((x,v),(y,w))$ on $(\mathbb{R}^{d}\times\mathbb{R}^{d})^{2}$ with marginals $f^{0}_{1}$ and $f^{0}_{2}$, and we consider $F((x,v),(y,w),t)$ starting from $F^{0}$ and solving

[TABLE]

This function is clearly nonnegative and has the correct marginals. For example, with $F_{1}(x,v,t)=\int F(x,y,v,w,t)dydw$, we see that

[TABLE]

because $\int F((x,\Phi(v,h)),(y,\Phi(w,h)),t)\det(D_{v}\Phi(w,h))dydw=F_{1}(x,\Phi(v,h),t)$: use the substitution $V=\Phi(w,h)$ (with $h$ fixed). Since $F_{1}(0)=f_{1}(0)$, we conclude that $F_{1}(t)=f_{1}(t)$. Hence we conclude that $d_{1}(f_{1}(t),f_{2}(t))\leq\int\hskip-4.0pt\int(|x-y|+|v-w|)F((x,v),(y,w),t)dxdydvdw$.

Next, using the equation for $F$, we find with $V=\Phi(v,h)$ and $W=\Phi(w,h)$,

[TABLE]

by (22). Since now $K\geq 1+KL$ by assumption, we deduce that

[TABLE]

and complete the proof as usual, minimizing in $F^{0}$. $\square$

Remark 9

Fix $a>0$. If using the Monge-Kantorovich distance with weight $\varrho=a|x-y|+|v-w|$, the condition $K\geq 1+KL$ is replaced by the condition $K>a+KL$.

3.3 Fractional heat equation with variable coefficients

Informally, the fractional Laplacian is a variant of the integral equation treated in Subsection 3.1. However there is a particular interest when the coefficients depend on space, an example we borrow from [21, 17]. Consider the parabolic equation with derivatives of order $\alpha\in(0,2)$

[TABLE]

Theorem 10

Assume that $\alpha\in(1,2)$ and that $\sigma\in C^{1/\alpha}$ and consider two initial probability densities $u^{0}_{1}$ and $u^{0}_{2}$ on $\mathbb{R}$ and the corresponding solutions $u_{1}$ and $u_{2}$ to (23). For all $t\geq 0$,

[TABLE]

Proof. We consider an initial probability density $v^{0}$ on $\mathbb{R}^{2}$ with marginals $u_{1}^{0}$ and $u_{2}^{0}$ and the solution $v$ to the problem (written in weak form): for all smooth $\varphi:\mathbb{R}^{2}\mapsto\mathbb{R}$,

[TABLE]

starting from $v^{0}$. The solution is clearly nonnegative and one checks as usual that for each $t\geq 0$, the marginals of $v(t)$ are $u_{1}(t)$ and $u_{2}(t)$: for example, we apply the above formula with $\varphi$ depending only on $x$ and deduce that $v_{1}(x,t)=\int v(x,y,t)dy$ solves the weak form of (23), whence $v_{1}=u_{1}$ since $v_{1}(0)=u_{1}(0)$. Consequently, we have $d_{\alpha-1}(u_{1}(t),u_{2}(t))\leq(\alpha-1)^{-1}\int\hskip-4.0pt\int|x-y|^{\alpha-1}v(x,y,t)dxdy$, and, using the same arguments as usual, it suffices to show that

[TABLE]

This follows from the fact that for all $x,y\in\mathbb{R}$, setting $u=\frac{\sigma(x)-\sigma(y)}{x-y}$,

[TABLE]

The proof of this last equality can be found in [17, Lemma 9-(ii)], case $a_{+}=a_{-}$ and $\beta=\alpha-1$. Observe that

[TABLE]

so that (3.3) makes sense with $\varphi(x,y)=|x-y|^{\alpha}$, thanks to our regularity assumption on $\sigma$. $\square$

Here again, as in Section 2.1, the main technical difficulty is to prove the well-posedness of (23), in particular when $\sigma$ may degenerate. This is useful to check that the solution $v$ to the coupled equation has the correct marginals.

4 Inhomogeneous integral equations

Our next purpose is to give an example on the way to take into account $x$-dependency in IPDE models, for instance when considering a measure $\mu(x,h)$ in the scattering equation (16). We exemplify this issue with a simple equation we borrow from [18]. Consider an interval $I$ of $\mathbb{R}$, a rate function $d\geq 0$ defined on $I$ and some probability density $b$ on $I$. We consider the conservative equation

[TABLE]

starting from an initial probability density $u^{0}$ on $I$. We notice at once that this equation makes sense for probability measures $u(dx,t)$ (for each $t\geq 0$, $u(dx,t)$ is a probability measure on $I$) in the following weak sense: for all smooth $\varphi:I\mapsto\mathbb{R}$,

[TABLE]

Theorem 11

Consider two probability densities $u^{0}_{1},u^{0}_{2}$ on $\mathbb{R}^{d}$ and the corresponding solutions $u_{1},u_{2}$ to (26). Under one of the two conditions (a) or (b) below, for all $t\geq 0$,

[TABLE]

(a) $I=\mathbb{R}_{+}$, $d(0)=0$, $d$ is increasing, $b=\delta_{0}$, and $\varrho(x,y)=|d^{p}(x)-d^{p}(y)|$ for some $p\geq 1$.

(b) $I=\mathbb{R}_{+}$, $d(x)=\alpha x^{p}+\beta$ for some $\alpha,\beta\geq 0$ and $p\geq 1$, with $\varrho(x,y)=|x^{p}-y^{p}|$, under the condition that $\beta\geq\alpha\int_{0}^{\infty}z^{p}b(z)dz$.

Other assumptions on $I$, $b$, $d$ are possible: it suffices that $\varrho$, $b$ and $d$ satisfy the dual inequality (28) below, which corresponds to (14) for the heat equation with variable coefficients.

Proof. We consider some probability density $v^{0}$ on $I^{2}$ with marginals $u_{1}^{0}$ and $u_{2}^{0}$ and define the probability measure $v(dx,dy,t)$ as solving, for all smooth $\varphi:I^{2}\mapsto\mathbb{R}$,

[TABLE]

It holds true that $v(t)$ is a probability measure on $I^{2}$ for each $t\geq 0$, and that its marginals are $u_{1}(t)$ and $u_{2}(t)$. For example, applying the coupling equation with $\varphi$ depending only on $x$ and using that

[TABLE]

one verifies that $\int_{y\in I}v(dx,dy,t)$ solves (26), whence $\int_{y\in I}v(dx,dy,t)=u_{1}(dx,t)$ by uniqueness. Hence for any cost function $\varrho:I^{2}\mapsto\mathbb{R}_{+}$, we have $d_{\varrho}(u_{1}(t),u_{2}(t))\leq\int\hskip-4.0pt\int\varrho(x,y)v(dx,dy,t)$. Furthermore, we easily compute, using that $\varrho(z,z)=0$ for all $z\in I$, that $b$ is a probability density, and that $\min(r,s)+(r-s)_{+}+(s-r)_{+}=\max(r,s)$,

[TABLE]

Therefore, using the same arguments as usual, the result will follows from the fact that for all $x,y\in I$,

[TABLE]

(a) Assuming that $I=\mathbb{R}_{+}$, that $d(0)=0$, that $d$ is increasing, that $\rho(x,y)=|d(x)-d(y)|^{p}$ for some $p\geq 1$ and that $b=\delta_{0}$, we check that, when e.g. $x\geq y\geq 0$,

[TABLE]

which holds true since indeed, for any $s\geq t\geq 0$, $(s^{p}-t^{p})s\geq t^{p}(s-t)$ because $p\geq 1$.

(b) Assume next that $I=\mathbb{R}_{+}$, $d(x)=\alpha x^{p}+\beta$ for some $\alpha,\,\beta\geq 0$ and $p\geq 1$ and that $\varrho(x,y)=|x^{p}-y^{p}|$. We have to verify that, for all $x\geq y\geq 0$,

[TABLE]

Setting $m=\int_{0}^{\infty}z^{p}b(z)dz$, it suffices to check that $\alpha x^{p}+\beta\geq\alpha(m+y^{p})$. This of course holds true if $\beta\geq\alpha m$. $\square$

Observe that the strong equation corresponding to the weak form (27) is nothing but

[TABLE]

5 Homogeneous Boltzmann equation

In his seminal paper [27], Tanaka observed that the homogeneous Boltzmann equation for Maxwell molecules is non-expansive for the 2-Monge-Kantorovich distance. This result was extended to inelastic collisions in [4] and a survey of results concerning homogeneous kinetic equations can be found in [13]. Also, [26] managed to study the corresponding dissipation in order to quantify the convergence to equilibrium of the solutions and, even more interesting, to prove the convergence to equilibrium of Kac particle system, with a rate of convergence not depending on the number of particles.

The homogeneous Boltzmann writes

[TABLE]

The collision kernel $B$ is assumed to satisfy $\int_{0}^{\pi}B(\theta)d\theta=1$. As is well-known, this equation writes, in weak form, for all mapping $\varphi:\mathbb{R}^{3}\mapsto\mathbb{R}$,

[TABLE]

Theorem 12

Consider two initial probability densities $f^{0}_{1},f^{0}_{2}$ on $\mathbb{R}^{3}$ with a finite moment of order $2$ and the corresponding solutions $f_{1},f_{2}$ to (29). Then, for all $t\geq 0$, one has

[TABLE]

Proof. We fix a probability density $F^{0}$ on $(\mathbb{R}^{3})^{2}$ with marginals $f^{0}_{1}$ and $f^{0}_{2}$ and build a coupled equation with the same principle as for scattering, that is the jump parameters are taken in common to the two variables, in such a way that the post-collisional velocities are as close as possible. We consider the solution $F(v,w,t)$, starting from $F^{0}$, to the following coupling equation written in weak form: for all mapping $\Psi:\mathbb{R}^{3}\times\mathbb{R}^{3}\mapsto\mathbb{R}$,

[TABLE]

where, for $v,w,v_{*},w_{*}\in\mathbb{R}^{3}$, $\theta\in(0,\pi)$ and $\varphi\in(0,2\pi)$, we have set

[TABLE]

where $I=\frac{(v-v_{*})\wedge(w-w_{*})}{|(v-v_{*})\wedge(w-w_{*})|}$, and $I_{1},\,I_{2}$ are chosen so that $(\frac{v-v_{*}}{|v-v_{*}|},I,I_{1})$ and $(\frac{w-w_{*}}{|w-w_{*}|},I,I_{2})$ are two direct orthormal bases, and where

[TABLE]

It is clear that $F$ remains nonnegative for all times. Also, it holds that $\int_{\mathbb{R}^{3}}F(v,w,t)dw=f_{1}(v,t)$ and $\int_{\mathbb{R}^{3}}F(v,w,t)dv=f_{2}(w,t)$. For example concerning $f_{1}$, we apply the weak coupling equation to some $\Psi$ depending only on $v$ and we show that $\int_{\mathbb{R}^{3}}F(v,w,t)dw$ solves (30). This follows from the fact that, when fixing $(v,w)$ and $(v_{*},w_{*})$, the expression between brackets in (31) only depends on $\sigma$, so that for any function $H:\mathbb{S}^{2}\mapsto\mathbb{R}$, we may write

[TABLE]

We conclude that

[TABLE]

Consequently, it holds that $d_{2}(f_{1}(t),f_{2}(t))\leq\int\hskip-4.0pt\int|v-w|^{2}F(v,w,t)dvdw=:h(t)$, and it suffices, as usual, to show that $h^{\prime}(t)\leq 0$. For this, it suffices to verify that for all fixed $v,w,v_{*},w_{*}\in\mathbb{R}^{3}$, all $\theta\in(0,\pi)$,

[TABLE]

A simple computation, using that $\int_{0}^{2\pi}\sigma d\varphi=\int_{0}^{2\pi}\omega d\varphi=0$, shows that

[TABLE]

We used that $|v-v_{*}||w-w_{*}|I_{1}\cdot I_{2}=(v-w)\cdot(v_{*}-w_{*})$. All in all, we arrive at

[TABLE]

and the proof is complete. $\square$

6 Porous media equation

We now consider the generalized porous media equation written, with $A:\mathbb{R}_{+}\mapsto\mathbb{R}$ of class $C^{2}$, as

[TABLE]

It was discovered in [24], see also [12], using a gradient flow approach, that this equation is non-expansive for $d_{2}$, under a few conditions on $A$ including convexity. The coupling method used in the whole present paper does not seem to apply directly. However, using Brenier’s map, this property follows as proved in [5] by an argument closely related to the coupling method. We present this argument, staying at an informal level.

Theorem 13

Consider some $C^{2}$ function $A:\mathbb{R}_{+}\mapsto\mathbb{R}$ such that $B(r)=\int_{0}^{r}wA^{\prime\prime}(w)dw\geq 0$ for all $r\geq 0$ and such that $r\mapsto r^{1/d-1}B(r)$ is non-decreasing. Consider two probability densities $u_{1}^{0},u_{2}^{0}$ on $\mathbb{R}^{d}$ and the corresponding solutions $u_{1},u_{2}$ to (32). Then for all $t\geq 0$, one has

[TABLE]

This applies to the porous media equation, i.e. with $A(u)=m^{-1}u^{m}$, as soon as $m\geq 1$. The justification of the computation requires at least that $\int_{\mathbb{R}^{d}}B(u_{1}(x,t))<\infty$. See [5] for the rigorous proof, which assumes that the solutions are smooth and positive.

Proof. We consider Brenier’s map [9] for $u_{1}^{0}$ and $u_{2}^{0}$, i.e. a convex function $\Phi:\mathbb{R}^{d}\mapsto\mathbb{R}$ such that $d_{2}(u_{1}^{0},u_{2}^{0})=\frac{1}{2}\int_{\mathbb{R}^{d}}|x-\nabla\Phi(x)|^{2}u_{1}^{0}(x)dx$ and $\nabla\Phi\#u_{1}^{0}=u_{2}^{0}$. We next consider the probability measure $v(dx,dy,t)$ (for each $t\geq 0$, $v(t)\in\mathcal{P}((\mathbb{R}^{d})^{2})$) solving the coupling equation

[TABLE]

and starting from $v^{0}(dx,dy)=u_{1}^{0}(dx)\delta_{\nabla\Phi(x)}(dy)\in K(u^{0}_{1},u^{0}_{2})$, defined by the formula $v^{0}(\mathcal{A})=\int_{R^{d}}{\bf 1}_{\{(x,\nabla\Phi(x))\in\mathcal{A}\}}u_{1}^{0}(x)dx$ for all Borel set $\mathcal{A}\subset\mathbb{R}^{d}$. Again, we only use the weak form

[TABLE]

for all smooth $\varphi:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}$.

One easily verifies, as usual, that for each $t\geq 0$, $v(t)$ has $u_{1}(t)$ and $u_{2}(t)$ for marginals: for example, applying the weak equation of $v$ to some $\varphi$ depending only on $x$ shows that $\int_{y\in\mathbb{R}^{d}}v(dx,dy,t)$ is a (weak) solution to (32) and since it starts from $u_{1}^{0}$, we conclude by uniqueness. As a conclusion, for all $t\geq 0$, $d_{2}(u_{1}(t),u_{2}(t))\leq I(t)$, where $I(t)=\frac{1}{2}\int\hskip-4.0pt\int|x-y|^{2}v(dx,dy,t)$. Next, we observe that

[TABLE]

where

[TABLE]

In particular, by definition of $v^{0}$,

[TABLE]

where we recall that $B(r)=\int_{0}^{r}wA^{\prime\prime}(w)dw$. Integrating by parts, we thus find

[TABLE]

This uses that for any convex function $\Phi:\mathbb{R}^{d}\mapsto\mathbb{R}$, we have $d^{-1}\Delta\Phi(x)\geq[\det(D^{2}\Phi(x))]^{1/d}$. But since $\nabla\Phi\#u_{1}^{0}=u_{2}^{0}$, we have, for any $\varphi:\mathbb{R}^{d}\mapsto\mathbb{R}$,

[TABLE]

so that $\det D^{2}\Phi(x)=u_{1}^{0}(x)/u_{2}(\nabla\Phi(x))$ (see [19] for an account on this Monge-Ampère equation). All in all, we have checked that

[TABLE]

Proceeding similarly, we see that

[TABLE]

Performing the substitution $x=(\nabla\Phi)^{-1}(y)$, we end with

[TABLE]

We thus find, with the notation $y=\nabla\Phi(x)$

[TABLE]

Since $r\mapsto r^{1/d-1}B(r)$ is non-decreasing by assumption, we conclude that $I^{\prime}(0)\leq 0$. The above considerations hold true at any time, and not only at $t=0$. In other words, for all $t\geq 0$, we can find a function $I_{t}:[t,\infty)\mapsto\mathbb{R}$ such that $I_{t}(t)=d_{2}(u_{1}(t),u_{2}(t))$, $I_{t}^{\prime}(t)\leq 0$ and $d_{2}(u_{1}(s),u_{2}(s))\leq I_{t}(s)$ for all $0\leq t\leq s$.One immediately concludes that for all $t\geq 0$,

[TABLE]

so that $t\mapsto d_{2}(u_{1}(t),u_{2}(t))$ is non-increasing. $\square$

7 An approach by duality

In order to complete the presentation, we quickly mention another possible and original approach, based on duality. We consider the simplest model, i.e. the heat equation in dimension $1$, but all the models treated in the present paper, except the porous media equation, may be treated similarly, with more complicated discretization procedures and more involved computations.

Proof of Theorem 1 when $d=1$ for $\rho(x,y)=|x-y|^{p}$ with $p\geq 1$. We consider two solutions $u_{1},u_{2}$ to (2), starting from probability measures $u_{1}^{0},u_{2}^{0}$ with finite $p$-moment. For $h>0$, we consider the solutions $u_{1,h},u_{2,h}$, starting from $u_{1}^{0},u_{2}^{0}$, to the discrete heat equation

[TABLE]

It can be written in weak form

[TABLE]

It is standard that $u_{1,h}\to u_{1}$ and $u_{2,h}\to u_{2}$ as $h\to 0$ (in the weak topology of measures for instance). We will verify that for each $h>0$, it holds that $d_{p}(u_{1,h}(t),u_{2,h}(t))\leq d_{p}(u_{1}^{0},u_{2}^{0})$, for any $p\geq 1$, and this will complete the proof.

We fix $p\geq 1$ and introduce the set $Q_{p}$ of pairs $(\varphi,\psi)$ of functions from $\mathbb{R}^{d}$ to $\mathbb{R}$ such that for all $x,y\in\mathbb{R}^{d}$, $\varphi(x)+\psi(y)\leq|x-y|^{p}$. For any pair of probability densities $f,g$ on $\mathbb{R}^{d}$, the Monge-Kantorovich distance can also be expressed by duality, see [28], as

[TABLE]

For $(\varphi,\psi)\in Q_{p}$, we set $\Delta_{\varphi,\psi}(t)=\int\varphi(x)u_{1,h}(x,t)dx+\int\psi(y)u_{2,h}(y,t)dy$. Using that $(\varphi(\cdot+h),\psi(\cdot+h))$ and $(\varphi(\cdot-h),\psi(\cdot-h))$ both belong to $Q_{p}$, we find

[TABLE]

This implies that

[TABLE]

Taking the supremum over all pairs $(\varphi,\psi)$ in $Q_{p}$ and dividing by $p$, we conclude that

[TABLE]

By the Gronwall lemma, we conclude that $e^{2h^{-2}t}d_{p}(u_{1,h}(t),u_{2,h}(t))\leq e^{2h^{-2}t}d_{p}(u_{1}^{0},u_{2}^{0})$ as desired. $\square$

Unfortunately, we are not able to use a similar procedure directly on the (non discretized) heat equation.