Franc ¸ois Cuvelier Nadia Belaribi

Francesco Russo

Probabilistic and deterministic algorithms for space multidimensional irregular porous media equation

Keywords: Stochastic particle algorithm, porous media equation, monotonicity, stochastic differential equations, non-parametric density estimation, kernel estimator AMS classification. MSC 2010: 65C05, 65C35, 82C22, 35K55, 35K65, 35R05, 60H10, 60J60, 62G07, 65M06

The object of this paper is a multi-dimensional generalized porous media equation (PDE) with not smooth and possibly discontinuous coefficient β, which is well-posed as an evolution problem in L 1 (R d ). This work continues the study related to the one-dimensional case by the same authors. One expects that a solution of the mentioned PDE can be represented through the solution (in law) of a nonlinear stochastic differential equation (NLSDE). A classical tool for doing this is a uniqueness argument for some Fokker-Planck type equations with measurable coefficients. When β is possibly discontinuous, this is often possible in dimension d = 1. If d > 1, this problem is more complex than for d = 1. However, it is possible to exhibit natural candidates for the probabilistic representation and to use them for approximating the solution of (PDE) through a stochastic particle algorithm. We compare it with some numerical deterministic techniques that we have implemented adapting the method of a paper of Cavalli et al. whose convergence was established when β is Lipschitz. Special emphasis is also devoted to the case when the initial condition is radially symmetric. On the other hand assuming that β is continuous (even though not smooth), one provides existence results for a mollified version of the (NLSDE) and a related partial integro-differential equation, even if the initial condition is a general probability measure.

Introduction

The main target of this work is to construct and implement a stochastic algorithm which approximates the solution of a multidimensional porous media type equation with monotone possibly irregular coefficient.

In the whole paper T will be a strictly positive real number and d a strictly positive integer. We consider the parabolic problem on R d given by

   ∂ t u(t, x) ∈ 1 2 ∆β(u(t, x)), t ∈]0, T ], u(0, x) = u 0 (dx), x ∈ R d , (1.1)
in the sense of distributions, where u 0 is an initial probability measure. If u 0 has a density we will still denote it by the same letter. We look for the solutions of (1.1) with time evolution in L 1 (R d ), i.e., u :]0, T ] × R d → R such that u(t, •) ∈ L 1 (R d ), ∀t ∈]0, T ]. We formulate the following assumptions; they will be in particular valid in this Introduction.

Assumption A.

(i) β : R → R such that β is monotone.

(ii) β(0) = 0 and β continuous at zero.

(iii) We assume the existence of λ > 0 such that (β + λid)(R) = (R) and id(x) ≡ x. There is ℓ ≥ 1 such that the following is verified.

Assumption B(ℓ). There exists a constant

C β > 0 such that |β(u)| ≤ C β |u| ℓ , ℓ ≥ 1.
The analysis of (1.1) is generally done in the framework of monotone partial differential equations, including the case when β is discontinuous. In that case, by filling the gaps β is considered as multi-valued. In this sequel of this introduction, for the sake of simplicity, we will almost always use a single-valued formulation.

We define Φ : R → R + , setting

Φ(u) =        β(u) u if u = 0, C if u = 0, (1.2)
where C ∈ [lim inf u→0 + Φ(u), lim sup

u→0 + Φ(u)].
Note that when β(u) = u.|u| m-1 , m > 1, the partial differential equation (PDE) in (1.1) is nothing else but the classical porous media equation. In this case Φ(u) = |u| m-1 2 . We are particularly interested in the case when β is continuous excepted for a possible jump at one positive point, say u c > 0. A typical example is:

β(u) = H(u -u c ).u, (1.3)
H being the Heaviside function and u c will be called critical value or critical threshold.

Definition 1.1. (i) We will say that the PDE in (1.1), or β is non-degenerate if there is a constant c 0 > 0 such that Φ ≥ c 0 , on each compact of R + .

(ii) We will say that the PDE in (1.1), or β is degenerate if lim u→0 + Φ(u) = 0.

Remark 1.2. (i) We observe that β may be neither degenerate nor non-degenerate.

(ii) β defined in (1.3) is degenerate. β(u) = (H(uu c ) + ǫ)u is non-degenerate, for any ǫ > 0.

Equation (1.3) constitutes a model intervening in some self-organized criticality (often called SOC) phenomena, see [START_REF] Bak | How Nature Works: The science of Self-Organized Criticality[END_REF] for a significant monography on the subject and [START_REF] Barbu | Stochastic porous media equations and self-organized criticality[END_REF][START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF] for recent related references. Sand piles are typical related models, which were first introduced in the discrete setting: for instance the BTW (Bak-Tang-Wiesenfeld) model, see [START_REF] Bak | Self-organized criticality[END_REF] and a refined version, the so-called Zhang model. Inspired from the latter model, [START_REF] Bantay | Avalanche dynamics from anomalous diffusion[END_REF] introduced continuous sand pile models, in which appear a porous media equation of the type (1.1) with β defined in (1.3). Two different effects appear: the avalanche and the regular arrival of sand. A natural description of the global phenomenon is a stochastic perturbation by noise of the mentioned equation, i.e. a generalized stochastic porous media equation. The two effects appearing in very different scales, there is sense to analyze them separately. The deterministic PDE (1.1) is a natural description of the avalanche effect and in this paper we concentrate on that one. Recent work related to self-organized criticality and SPDEs was done by [START_REF] Barbu | Self-organized criticality via stochastic partial differential equations, Potential theory and stochastics in Albac[END_REF][START_REF]Stochastic porous media equations and self-organized criticality[END_REF].

In the one-dimensional case, under Assumption B [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF]Proposition 3.4] proved existence and uniqueness of solutions (in the sense of distributions) for (1.1) when the initial condition u 0 is a bounded integrable function. Indeed that result was essentially a clarification of older and celebrated results quoted in [START_REF] Benilan | The continuous dependence on ϕ of solutions of u t -∆ϕ(u) = 0[END_REF]. In Proposition 3.1, we extend that result when the dimension d is greater than 1, under the validity of Assumptions A and B(ℓ) for some ℓ ≥ 1.

As we mentioned, the paper focuses on the possible probabilistic representation for (1.1) when d ≥ 2, in the sense expressed below. We are interested in finding a process (Y t ) such that the marginal laws u(t, •) admit a density, that we will still denote by the same letter, and u is a solution of (1.1). In fact, we look for it in the form of a stochastic non-linear diffusion, i.e., a solution of a non-linear stochastic differential equation (NLSDE) of the type

           Y t =Y 0 + t 0 Φ d (u(s, Y s ))dW s , u(t, •) =Law density of Y t , ∀t ∈]0, T ], u(0, •) =u 0 , (1.4)
where Φ d (u) = Φ(u).I d and I d is the unit matrix on R d . Y 0 is an R d -valued, u 0 -distributed random variable and W is a d-dimensional classical Brownian motion independent of Y 0 .

To the best of our knowledge the first author who considered a probabilistic representation for the solutions of non-linear deterministic partial differential equations was McKean [START_REF] Jr | Propagation of chaos for a class of non-linear parabolic equations[END_REF]. However, in his case, the coefficients were smooth. In the one-dimensional case, a probabilistic interpretation of (1.1) when β(u) = u.|u| m-1 , m > 1, was provided in [START_REF] Benachour | Processus associés à l'équation des milieux poreux[END_REF]. Since the original article of McKean, many papers were produced and it is impossible to list them all: [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF] provides a reasonable list. If β(u) = u.|u| m-1 , m ∈]0, 1[, the partial differential equation in (1.1) is in fact the so-called fast diffusion equation. In the case when d = 1, [START_REF] Belaribi | About Fokker-Planck equation with measurable coefficients: application to the fast diffusion equation[END_REF] provides a probabilistic representation for the Barenblatt type solutions of (1.1).

Under Assumptions A and B [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], supposing that u 0 has a bounded density, [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF] (resp. [13]) proves existence and uniqueness of the probabilistic representation (in law) when β is nondegenerate (resp. degenerate). Besides, a theoretical probabilistic representation of the PDE perturbed by a multiplicative noise is given in [START_REF] Barbu | A stochastic Fokker-Planck equation and double probabilistic representation for the stochastic porous media type equation[END_REF].

Earlier, in the multi-dimensional case [START_REF] Jourdain | Propagation of chaos and fluctuations for a moderate model with smooth initial data[END_REF] concentrated on the case when β is non-degenerate and Φ is Lipschitz, continuously differentiable at least up to order 3, and with some further regularity assumptions on u 0 . The authors established existence and uniqueness of the probabilistic representation and the so-called propagation of chaos, see [START_REF] Sznitman | Topics in propagation of chaos[END_REF] for a rigorous formulation of this concept. When β is not smooth, the probabilistic representation remains an open problem in the multidimensional case.

The connection between the solutions of (1.4) and (1.1) is given by the following result.

Proposition 1.3. Let us assume the existence of a solution Y for (1.4). Let u :]0, T ] × R d → R + be a Borel function such that u(t, •) is the law density of Y t , t ∈]0, T ]. Then u provides a solution in the sense of distributions of (1.1) with u 0 = u(0, •).

The proof of previous result is well-known, but we recall here the basic argument.

Proof of Proposition 1.3. Let t ∈ [0, T ] and ϕ ∈ D(R d ), Y be a solution of (1.4). We apply Itô's formula to ϕ(Y t ) to obtain

ϕ(Y t ) = ϕ(Y 0 ) + d i=1 t 0 ∂ y i ϕ(Y s )Φ(u(s, Y s ))dW i s + 1 2 t 0 ∆ϕ(Y s )Φ 2 (u(s, Y s ))ds.
Taking the expectation we get Using then integration by parts and the fact that, according to (1.2) β(u) = Φ 2 (u)u, the expected results follows. ✷ Proposition 1.3 constitutes the easy step in the analysis of the probabilistic representation. The most difficult and interesting case consists in writing the converse implication. When β is non-degenerate, one idea for this consists in taking the solution u of (1.1). By [27, Section 6, Theorem 1] there is a solution (Y t ) in law of the stochastic differential equation

Y t = ξ + t 0 Φ d (u(s, Y s ))dW s , (1.5)
where ξ is a random element, distributed according to u 0 , independent of an underlying ddimensional Brownian motion W . The remaining difficult part consists in identifying the marginal laws of (Y t ), t ∈]0, T ] with u(t, x)dx. A general tool for this is a uniqueness theorem for Fokker-Planck type equations with measurable coefficients, of the type

∂ t u(t, x) = ∆(a(t, x)u(t, x)), t ∈ [0, T ], u(0, •) = u 0 (dx), (1.6)
with a(t, x) = 1 2 Φ 2 (u(t, x)). When d = 1 and a is bounded this was the object of [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF]Theorem 3.8]. Extensions where considered in [START_REF] Belaribi | About Fokker-Planck equation with measurable coefficients: application to the fast diffusion equation[END_REF]Theorem 3.1] when a is possibly degenerate and is allowed to be unbounded under some technical conditions. In fact [START_REF] Belaribi | About Fokker-Planck equation with measurable coefficients: application to the fast diffusion equation[END_REF]Theorem 3.1] also deals with the multidimensional case. When a is bounded and given two measure-valued solutions z 1 and z 2 , the Fokker-Planck uniqueness theorem applies saying that z 1 = z 2 whenever (z 1z 2 )(t, •) has a density z(t, •) for almost all t and z belongs to [START_REF] Krylov | Controlled diffusion processes[END_REF]Section 2.3,Theorem 4], show that whenever Y is a solution of (1.5), the measure f → E T 0 f (t, Y t )dt admits a density which belongs to L p ([0, T ]×R d ) for p ≤ d+1 d . When d = 1, p can be chosen equal to 2 and the Fokker-Planck uniqueness theorem applies, which allows [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF] to prove the probabilistic representation when β is non-degenerate. If d > 1, Alexandrov-Pucci-Krylov estimates are not enough to fulfill the assumptions of [START_REF] Belaribi | About Fokker-Planck equation with measurable coefficients: application to the fast diffusion equation[END_REF]Theorem 3.1].

L 2 ([0, T ] × R d ). If d ≥ 1, the Alexandrov-Krylov-Pucci estimates, see
In the present paper, we do not solve the general problem of the probabilistic representation for (1.1), however, the solutions to (1.5) are natural candidates and we establish some theoretical related results. A mollified version of (1.4) will be stated in (4.1). It consists in replacing in the first line of (1.4), Φ d (u) with Φ d (K H * v H ), where v H are the marginal laws of the solution Y of (4.1) and K H is a mollifying kernel. In Proposition 4.3, at least when Φ is bounded, nondegenerate and continuous, we establish the equivalence between the existence for (4.1) and the existence for the corresponding integro-differential deterministic PDE (4.11). In Proposition 4.2, under the same assumptions on Φ, we prove existence in law for (4.1). This also provides existence for (4.11).

Since a basic obstacle arises in dimension d > 1, a natural simplified problem to study is the probabilistic representation of (1.1) when the initial condition u 0 only depends on the radius (radially symmetric). In fact, in that case the problem can be reduced to dimension one, studying the stochastic differential equation fulfilled by the norm of the process at power d: this is the object of section 5. The reduction to dimension 1 has also the interest to build a new class of non-linear diffusions with singular coefficients. Unfortunately also that type of equation is difficult to handle at the theoretical level since, even when Φ = 1 (and so the unique solution Y to (1.4) is a classical d-dimensional Brownian motion), some non-Lipschitz functions at zero naturally appear. In that case the norm of Y is a d-dimensional Bessel process. When β is non-degenerate, the mentioned norm behaves then similarly to a d-dimensional Bessel process. Since that Bessel type process is recurrent for d ≤ 2 and transient for d > 2, it is expected that the zero (which is a singularity) point is much less visited when d > 2.

As mentioned in the beginning of the introduction one purpose of the present paper is to exploit the probabilistic representation in order to simulate the solutions of (1.1). For this we will implement an Euler scheme for stochastic differential equations and a non-parametric density estimation method using Gaussian kernel estimators, see [START_REF] Wand | Kernel smoothing[END_REF]. Since we expect our methods to be robust when the coefficient Φ is irregular, it is not reasonable to make use of higher order discretization schemes involving derivatives of Φ. Concerning the choice of the smoothing parameter ε for the density estimate, we extend to the multidimensional case the techniques used in [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF].

Besides, we have carried out a deterministic numerical method in dimension d = 2 of space, based on a sophisticated procedure developed in [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF], which is one of the most recent references in the subject. In fact, Cavalli et al. [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF] coupled WENO (weighted essentially non-oscillatory) interpolation methods for space discretization, see [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Advanced numerical approximation of nonlinear hyperbolic equations[END_REF], in order to prevent the onset of spurious oscillations, with IMEX (implicit explicit) Runge-Kutta schemes for time advancement, see [START_REF] Pareschi | Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], to obtain a high order method.

The general stochastic particle algorithm is empirically investigated in dimension d = 2. This is done in the following cases. When β(u) = u 3 comparing with the Barenblatt exact solutions and when β is defined by (1.3) comparing with the deterministic numerical technique; indeed in that case exact explicit solutions for (1.1) are unknown.

In the radially symmetric case, we implement the stochastic particle algorithm to the onedimensional reduced non-linear stochastic differential equation. As mentioned earlier when d = 2 (resp. d > 2), in most cases, the solutions are expected to be recurrent (resp. transient); so we suspect that the simulations are more performing when d > 2. However, our experiments show that the error of the approximations with respect to the exact solutions (derived by the classical porous media equation) stably reasonably low for all the values of d including d = 2. In sections 8.1 and 8.2 we compare the radial reduction method with the deterministic approach when d = 2.

The paper is organized as follows. After this introduction and some preliminaries, in section 3 we discuss the well-posedness of the deterministic problem (1.1). In section 4 we discuss the existence of a mollified version of the non-linear stochastic differential equation and its equivalence to the existence of an integro-differential PDE which is a regularized version of (1.1). Section 5 handles the case of a radially symmetric initial condition. In section 6 we describe the general particle algorithm for d = 2. Section 7 summarizes the deterministic technique developed in [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF] and finally section 8 is devoted to numerical experiments.

Preliminaries

In the whole paper, we will denote by We denote by M(R d ) the set of finite measures. We define a multivariate mollifier K H setting,

K H (x) = |H| -1 2 K(H -1 2 x), x ∈ R d , (2.1)
where K is a fixed d-variate smooth C ∞ probability kernel, typically a Gaussian kernel, and H is a symmetric strictly definite positive d × d matrix.

Definition 2.1. Let u 0 be a finite Borel measure on R d . We say that u :]0, T ] × R d → R is a solution in the sense of distributions of (1.1) with initial condition u 0 , if there is

η u :]0, T ]×R d → R, u(t, •), η u (t, •) ∈ L 1 loc (R d )
for almost all t, and

R d u(t, x)ϕ(x)dx = R d u 0 (dx)ϕ(x) + 1 2 t 0 ds R d η u (s, x)∆ϕ(x)dx (2.2)
for all ϕ ∈ D(R d ), and

η u (t, x) ∈ β(u(t, x)) for dt ⊗ dx-a.e. (t, x) ∈ [0, T ] × R d (2.3)
Remark 2.2. (i) By an obvious identification we can also consider

u :]0, T ] → L 1 loc (R d ) ⊂ D ′ (R d ).
(ii) If u is a solution of (1.1) with initial condition u 0 , then u :]0, T ] → D ′ (R d ) extends to a weakly continuous function [0, T ] → D ′ (R d ) still denoted by u such that u(0) = u 0 .

Estimates for the solution of the deterministic equation

Proposition 3.1. Let u 0 ∈ L 1 L ∞ (R d ), u 0 ≥ 0.
We suppose the validity of Assumptions A and B(ℓ) for some ℓ ≥ 1. Then there is a unique solution in the sense of distributions u ∈

(L 1 L ∞ )([0, T ] × R d ) of ∂ t u ∈ 1 2 ∆β(u), u(0, x) = u 0 (x), (3.1
)

with corresponding η u (t, •) ∈ L ∞ ([0, T ] × R d ). Furthermore, ||u(t, .)|| ∞ ≤ ||u 0 || ∞ for every t ∈ [0, T ] and there is a unique version of u such that u ∈ C([0, T ] ; L 1 (R d )) (⊂ L 1 ([0, T ] × R d )). Remark 3.2. (i) An immediate consequence of previous result is that u ∈ L p ([0, T ] × R d )
for every p ≥ 1.

(ii) Assumption B on β is more general than the case of [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF]13] stated for d = 1. In that case we had Assumption B(1).

(iii) Indeed, most of the arguments of the proof of Proposition 3.1 appear implicitly in [START_REF] Benilan | The continuous dependence on ϕ of solutions of u t -∆ϕ(u) = 0[END_REF] and related references. For the comfort of the reader we decided to give an independent complete proof of Proposition 3.1.

Proof of Proposition 3.1. See Appendix 9.1. ✷

The mollified non-linear stochastic differential equation

We suppose again that u 0 (dx) is a probability measure. We consider the following mollified non-linear diffusion equation (in law):

           Y H t =Y 0 + t 0 Φ d ((K H * v H )(s, Y H s ))dW s v H (t, •) =Law density of Y H t , ∀t ∈]0, T ] v H (0, •) =u 0 = Law density of Y 0 , (4.1)
where, Φ d ((K H * v H )) = Φ((K H * v H )).I d and I d is the unit matrix on R d and W is an underlying classical Brownian motion. Remark 4.1. Suppose β non-degenerate. If Φ is supposed to be Lipschitz and continuously differentiable at least up to order 3 and u 0 is absolutely continuous with density in H 2+α for some 0 < α < 1, [25, Proposition 2.2] states existence (even strong existence) and uniqueness of solutions to (4.1).

Existence of solutions for the mollified NLSDE

The result below affirms, in the non-degenerate case, the existence of solutions to (4.1) in the case when Φ is bounded and continuous. Proposition 4.2. Suppose that β is non-degenerate and Assumption B(1) holds. Furthermore, assume that Φ is continuous. Then, the problem (4.1) admits existence in law.

Proof of Proposition 4.2. The assumptions imply the existence of constants c 0 , c 1 > 0 such that c 0 ≤ Φ ≤ c 1 .

Let X = (X t ) t∈[0,T ] be the canonical process on the canonical space C([0, T ]) equipped with its Borel σ-field. We define (ρ δ ) δ>0 to be a family of Gaussian mollifiers converging to the Dirac delta measure and we set Φ δ = Φ * ρ δ . We define u δ

0 := u 0 1 [-1 δ , 1 δ ] * ρ δ,d , where ρ δ,d (y) = d i=1 ρ δ (y i ) if y = (y 1 , . . . , y d ).
We observe that u δ 0 belongs to C ∞ b , for each δ > 0 . Similarly, we consider the d × d matrix Φ δ,d = (Φ * ρ δ )I d , where I d is the unit matrix on R d . We remark that u δ 0 belongs to the space H 2+α (0 < α < 1) considered in [START_REF] Jourdain | Propagation of chaos and fluctuations for a moderate model with smooth initial data[END_REF]Notations]. Since Φ δ and all its derivatives are bounded, Remark 4.1 implies existence (even strong existence) for

the problem      X t =X 0 + t 0 Φ δ,d ((K H * P s )(X s ))dW s , t ∈ [0, T ], P : Law of X, X 0 ∼ u δ 0 (4.2)
where P s denotes the (marginal) law of X s under P . We denote by P := P δ the corresponding probability solving (4.2). In particular, X 0 is square integrable under P δ . Taking into account that u δ 0 → u 0 in law, since (Φ δ ) δ>0 are bounded by Φ ∞ , using Kolmogorov lemma with caution it is possible to show that the family (P δ ) is tight, see [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]Section 2.4,Problem 4.1]. Consequently, by relative compactness, there is a sequence (P δn ), that we will denote (P n ), which converges weakly to some probability P .

It remains to show that P solves the martingale problem related to (4.1). In particular, we will prove that the process

(MP) f (X t ) -f (X 0 ) - 1 2 t 0 ∆f (X s )Φ 2 ((K H * P s )(X s )) ds, t ∈ [0, T ],
is an (F s )-martingale, where (F s ) is the canonical filtration associated with X.

Let E (resp. E n ) be the expectation operator with respect to P (resp. P n ). Let 0 ≤ s < t ≤ T and R = R(X r , r ≤ s) be an (F s )-measurable, bounded and continuous random variable with respect to C([0, s]). In order to show the martingale property (MP) of X, we have to prove that

E f (X t ) -f (X s ) - 1 2 t s ∆f (X r )Φ 2 ((K H * P r )(X r ))dr R = 0, (4.3) for every f ∈ C 2 0 (R d ). Let f ∈ C 2 0 (R d ).
Since (X t ) t∈[0,T ] under P n is a solution of (4.2) with δ = δ n , we have

E n f (X t ) -f (X s ) - 1 2 t s ∆f (X r )Φ 2 δn ((K H * P n r )(X r ))dr R = 0. (4.4)
In order to take the limit in (4.4) we will only have to show that

lim n→+∞ E n [F n (X)] -E [F (X)] = 0, (4.5)
where

F n (ℓ) = t s ∆f (ℓ(r))Φ 2 δn ((K H * P n r )(ℓ(r)))drR(ℓ(ξ), ξ ≤ s), F (ℓ) = t s ∆f (ℓ(r))Φ 2 ((K H * P r )(ℓ(r)))drR(ℓ(ξ), ξ ≤ s).
Since the family of laws (P n ) converges to P , then the sequence of time-marginal laws (P n r ) converges to P r , for every r ≥ 0 and

lim n→+∞ (K H * P n r )(x) = lim n→+∞ R d K H (x -y)P n r (dy), =(K H * P r )(x), ∀x ∈ R d . (4.6)
Now , we split the left-hand side of (4.5) into I 1 (n) + I 2 (n), where

I 1 (n) =E n t s dr∆f (X r ){Φ 2 δn ((K H * P n r )(X r )) -Φ 2 ((K H * P r )(X r ))}R(X ξ , ξ ≤ s) , I 2 (n) =E n t s ∆f (X r )Φ 2 ((K H * P r )(X r ))drR(X ξ , ξ ≤ s) (4.7) -E t s ∆f (X r )Φ 2 ((K H * P r )(X r ))drR(X ξ , ξ ≤ s) .
Since c 0 ≤ Φ δ ≤ c 1 , according to [27, Section 2.3, Theorem 3] and taking into account the notations at the beginning of [27, Section 2.2] then, for every ϕ ∈ D([0, T ] × R d ) there is a constant A = A(c 0 , c 1 ), such that

E n T 0 ϕ(t, X t )dt ≤ A ϕ L d+1 ([0,T ]×R d )) . (4.8)
This implies that the measure ϕ → T 0 P n s (dy)ϕ(s, y) admits a density in

L p ′ ([0, T ] × R d ))
, where p ′ = d+1 d . We denote them by (s, y) → q n (s, y).

Moreover, (4.8) implies that sup

n≥1 [0,T ]×R d |q n (s, y)| p ′ dsdy < +∞. (4.9)
Consequently, there is a subsequence converging weakly in L p ′ to some (s, y) → q(s, y) which belongs to L p ′ . Since the sequence (P n ) converges weakly to P , it follows that the measure ϕ → E T 0 ϕ(t, X t )dt admits q as density. We are now able to show that lim

n→+∞ I 1 (n) = 0. In fact, I 1 (n) is bounded by R ∞ E n t s ∆f (X r ) Φ 2 δn ((K H * P n r )(X r )) -Φ 2 ((K H * P r )(X r )) dr , = R ∞ t s dr R d q n (r, y) ∆f (y) Φ 2 δn ((K H * P n r )(y)) -Φ 2 ((K H * P r )(y)) dy.
Previous expression is bounded by

R ∞ ∆f ∞      [s,t]×R d |q n (r, y)| p ′ drdy      1 p ′ ×      [s,t]×suppf drdy Φ 2 δn ((K H * P n r )(y)) -Φ 2 ((K H * P r )(y)) p      1 p , (4.10)
where p = d + 1. The first integral in (4.10) is bounded by (4.9). By (4.6), since Φ is continuous, Lebesgue dominated convergence theorem implies that the second integral in (4.10) converges to zero as n goes to infinity. this shows that lim n→+∞ I 1 (n) = 0. On the other hand, lim n→+∞ I 2 (n) = 0, because the family of laws (P n ) converges to P and ∆f , R, Φ(K H * P r ) are continuous and bounded for fixed r.

This concludes the proof of Proposition 4.2. ✷

Some complements concerning the mollified equations

We recall that u 0 is a general real probability measure. If Φ is uniformly continuous, we will prove in the sequel that the existence for the mollified NLSDE (4.1), is equivalent to the existence for the following non-

linear integro differential PDE    ∂ t v H (t, •) = 1 2 ∆ Φ 2 ((K H * v H )(t, x))v H (t, •) , t ∈ [0, T ] v H (0, •) = u 0 (dx). ( 4
Y t = Y 0 + t 0 A(s, Y s )dW s , (4.12)
where

A(t, y) = Φ((K H * v)(t, y))I d .
Again by Itô's formula, then the family of marginal laws of Y , z(t, dy), solve

       ∂ t z(t, •) = 1 2 d i=1 ∂ 2 x 2 i Φ 2 ((K H * v)(t, x))z(t, •) , t ∈ [0, T ] z(0, •) = u 0 (dx), (4.13)
in the sense of distributions. Another obvious solution of (4.13), is provided by v which is a solution of (4.11).

In order to identify v with z it will be helpful to prove uniqueness for the solutions of (4.13) in the class of weakly continuous solutions [0, T ] → M(R d ). This will imply that z ≡ v and this will allow to conclude the proof. The key result for doing this is [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF]Lemma 2.3]. Indeed, since the coefficients in (4.13) are continuous bounded and non-degenerate, that lemma leads to the result if we prove uniqueness in law for the family of

Y t = x + t 0 A(s, Y s )dW s , (4.14)
for any x ∈ R d . The validity of previous uniqueness follows by [START_REF] Stroock | Multidimensional diffusion processes[END_REF]Theorem 7.2.1], with γ = AA t and b = 0.

In our case γ(s, x) = Φ 2 (K H * v(s, x))I d . β being non-degenerate, we obviously get condition (2.1) of [35, Theorem 7.2.1], i.e., inf 0≤s≤T inf θ∈R d < θ, γ(s, x)θ > / θ 2 > 0.
It remains to check the corresponding condition (2.2), i.e., lim

y→x sup s∈[0,T ] γ(s, y) -γ(s, x) ∞ = 0, (4.15)
where • ∞ is the supremum of matrix components. Actually, (4.15) is equivalent to

lim y→x sup s∈[0,T ] |Φ 2 ((K H * v)(s, y)) -Φ 2 ((K H * v)(s, x))| = 0. (4.16)
Let k > 0. Since Φ is uniformly continuous then Φ 2 is also uniformly continuous; so there exists δ > 0 such that

|Φ 2 (w 1 ) -Φ 2 (w 2 )| ≤ k, if |w 1 -w 2 | < δ.
Besides, we have

|(K H * v)(s, y) -(K H * v)(s, x)| = R d |(K H (y -z) -K H (x -z)) v(s, dz)| . (4.17)
Therefore, (4.17) is bounded by

x -y ∇K H ∞ T 0 ds v(s, •) var ,
where

T 0 ds v(s, •) var < +∞ because v : [0, T ] → M(R d ) is weakly continuous. Thus, choosing x -y < δ ∇K H ∞ T 0 ds v(s,•) var , gives |(K H * v)(s, y) -(K H * v)(s, x)| < δ. Consequently, sup s∈[0,T ] |Φ((K H * v)(s, y)) -Φ((K H * v)(s, x))| < k.
This concludes the proof of (4.16). Finally, equation (4.14) admits uniqueness in law and the result follows. ✷

The next step should be to prove the convergence of the solution v H of (4.11) to the solution u of (1.1). At this stage we are note able to prove it without assuming that Φ has some smoothness.

Reduction to dimension 1 when the initial condition is radially homogeneous

From now on, without restriction of generality, d will be greater or equal to 2.

Some mathematical aspects of the reduction

In this section we are interested in the solutions of (1.1) whose initial condition is radially symmetric.

From now on R t will denote the transpose of a generic matrix R.

An orthogonal matrix R ∈ R d ⊗ R d is a matrix such that RR t = R t R = I d , where I d is the identity matrix on R d . We denote by O(d) the set of d × d orthogonal matrices. Given a σ-finite Borel measure µ on (R d \{0}), R ∈ O(d), we define µ R as the σ-finite Borel measure such that R d µ R (dx)ϕ(x) = R d µ(dx)ϕ(R -1 x).
If µ is absolutely continuous with density f then µ R is absolutely continuous with density f R :

R d → R, where f R (x) = f (Rx). If u :]0, T ] × R d → R, we set u R (t, x) = u(t, Rx), t ∈]0, T ], x ∈ R d . Definition 5.1. (i) µ is said radially symmetric if for any R ∈ O(d) we have µ R = µ. (ii) A function u 0 : R d → R is said radially symmetric if there is u 0 :]0, +∞[→ R such that u 0 (x) = ū0 ( x ), ∀x = 0. Remark 5.2. If u 0 ∈ L 1 loc (R d \{0}
) then u 0 is radially symmetric if and only if the σ-finite measure u 0 (x)dx is radially symmetric.

Let J be a class of finite Borel measures on R d which is invariant through the action of every orthogonal matrix. Let U be a class of weakly continuous u : [0, T ] → M(R), t → u(t, •) such that for almost all t ∈]0, T ] u(t, •) admits a density, still denoted by u(t, x), x ∈ R d and such that u R ∈ U for any R ∈ O(d). We suppose that (1.1) is well-posed in U for every u 0 ∈ J.

Remark 5.3. Suppose that Assumptions A and B(ℓ), for some ℓ ≥ 1, are fulfilled. A classical choice of J (resp. U) is the cone of bounded non-negative integrable functions on R d (resp.

L 1 L ∞ ([0, T ] × R d )).
We first observe that whenever the initial condition of (1.1) is radially symmetric then the solution conserves this property.

Proposition 5.4. Let u 0 be a finite Borel measure on R d . Let u :]0, T ] × R d → R such that u(t, •) ∈ L 1 (R d ), ∀t ∈]0, T ]. Let u be a solution in the sense of distributions of (1.1) with u 0 as initial condition. (i) Let R ∈ O(d). Then u R is again a solution in the sense of distributions of (1.1), with initial condition u R 0 . (ii) If u 0 ∈ J and u ∈ U then there is ū :]0, T ]×]0, +∞[→ R such that u(t, x) = ū(t, x ), ∀t ∈]0, T ], x ∈ R d \{0}.
Proof of Proposition 5.4. (i) Let ϕ : R d → R be a smooth function with compact support. Since |det(R)| = 1, taking into account that u solves (1.1) we get,

R d u R (t, x)ϕ(x)dx = R d u(t, x)ϕ R -1 (x)dx = R d ϕ R -1 (x)u 0 (dx) + 1 2 t 0 ds R d η u (s, x)∆ϕ R -1 (x)dx,
where η u (s, x) ∈ β(u(s, x)) dsdx a.e. Previous sum is equal to

R d ϕ R -1 (x)u 0 (dx) + 1 2 t 0 ds R d η u (s, x)(∆ϕ) R -1 (x)dx, (5.1) since
∆ ϕ S (x) = (∆ϕ) (Sx), (5.2)
for an orthogonal matrix S; here S = R -1 . We shortly prove (5.2). We recall that

D 2 ϕ is a bounded bilinear form on R d . If e, f ∈ R d , we have D 2 ϕ S (x)(e, f ) = (D 2 ϕ)(Sx)(Se, Sf ).
Let (e i ) 1≤i≤d be an orthonormal basis of R d . We write

∆ϕ S (x) = Tr D 2 ϕ S (x) = d i=1 D 2 ϕ S (x)(e i , e i ) = d i=1 (D 2 ϕ)(Sx)(Se i , Se i ) = Tr (D 2 ϕ)(Sx) , since (Se i ) 1≤i≤d is still an orthonormal basis of R d . Finally (5.2) is established. Then, (5.1) gives R d u R (t, x)ϕ(x)dx = R d ϕ(x)u R 0 (dx) + 1 2 t 0 ds R d η u (s, Rx)(∆ϕ)(x)dx. Since η u (s, Rx) ∈ β(u(s, Rx)) = β(u R (s, x)) dsdx a.e., this establishes (i). (ii) According to Remark 5.2, it is enough to show that u(t, x) = u(t, Rx), ∀(t, x) ∈]0, T ] × R d , for every R ∈ O(d). Since u 0 is radially symmetric, item (i) implies that for any R ∈ O(d),
u R is a solution of (1.1) with u 0 as initial condition. Since u, u R ∈ U, we get u = u R and so, item (ii) follows.

✷

From now on, we will suppose the validity of Assumptions A and B(ℓ), for some ℓ ≥ 1. Let u 0 ∈ J, u ∈ U solution of (1.1) in the sense of distributions.

Remark 5.5. By Proposition 5.4(ii), there is

ũ :]0, T ] × R + → R such that u(t, x) = ũ(t, x d ), t ∈]0, T ], x ∈ R d .
We are now interested in the stochastic differential equation solved by the process (S t ), being defined as follows:

∀t ∈ [0, T ], S t = Y t d = d i=1 (Y i t ) 2 d 2 , (5.3)
where (Y t ) is a given solution of the d-dimensional problem (1.4), which in this section is supposed to exist. We will denote by ν(t, .) the law of S t . We first state a result concerning the relation between ν and ũ.

Lemma 5.6. For almost all t ∈]0, T ], ν(t, •) admits a density ρ → ν(t, ρ) verifying

ν(t, ρ) = C d ũ(t, ρ), ∀ρ > 0. (5.4)
where,

C = 2(π) d 2 Γ( d 2 ) , (5.5)
and Γ is the usual Gamma function. In particular, this gives

C =            2(π) d 2 ( d 2 -1)! , if d is an even number 2 d+1 2 π d-1 2 1 × 3 × 5 × . . . × (d -2) , otherwise.
Remark 5.7. The statement of Proposition 5.4 could allow to define ũ such that u(t, x) = ũ(t, x γ ) for a generic γ > 0, which could also be equal to 1 or 2. The choice of taking γ = d is justified by Lemma 5.6 above. If we take a different γ the quotient ν ũ in (5.4) would be proportional to some power of ρ, producing a bad numerical conditioning.

Proof of Lemma 5.6. Let f be a continuous and bounded function on R + . Since ν(t, •), t ∈ ]0, T ], is the law density of S t , we have

E(f (S t )) = R+ f (ρ)ν(t, ρ)dρ.
(5.6)

Since S t is defined by (5.

3), we have

E(f (S t )) = E(f ( Y t d )) = R d f ( y d )u(t, y)dy.
By Remark 5.5, we obtain

E(f (S t )) = R d f ( y d ) ũ(t, y d )dy.
Using the change of variables with hyperspherical coordinates, we get

E(f (S t )) = R+ Cr d-1 f (r d ) ũ(t, r d )dr,
where, C is given in (5.5). Finally, setting the change of variables ρ = r d , and identifying the result with (5.6), we obtain formula (5.4), for ν. ✷ (iii) If µ 0 is a radially symmetric σ-finite measure on R d \{0} we denote μ0 the σ-finite measure on ]0, +∞[ defined by

+∞ 0 C d μ0 (dr) ψ(r) = R d µ 0 (dx)ψ(x)dx, (5.7) if ψ ∈ D(]0, +∞[), ψ(x) = ψ( x d ). (iv) In particular ν 0 is the law of Y 0 d , i.e., ν 0 = C d u 0 .
From now on we will suppose that β is single-valued. ν defined in (5.4) is a solution of a partial differential equation that we determine below. Proposition 5.9. Let u 0 ∈ J and u ∈ U be the solution in the sense of distributions of the ddimensional problem (1.1) and assume that v(t, •) is defined by (5.4). Then t → v(t, •) verifies in the sense of distributions the following PDE in C(]0, T ]):

∂ t v(t, ρ) = C(1 -d)∂ ρ ρ 1-2 d β d C v(t, ρ) + Cd 2 ∂ 2 ρρ ρ 2-2 d β d C v(t, ρ) , (5.8
)

with initial condition v 0 = C d ũ0
, where ũ0 is defined in (5.7). C is the constant defined in (5.5).

This means in particular that for every

ψ ∈ D(]0, +∞[), +∞ 0 v(t, ρ)ψ(ρ)dρ = +∞ 0 v 0 (ρ)ψ(ρ)dρ -C(1 -d) +∞ 0 ψ ′ (ρ)ρ 1-2 d β d C v(t, ρ) dρ + Cd 2 +∞ 0 ψ ′′ (ρ)ρ 2-2 d β d C v(t, ρ) dρ.
Proof of Proposition 5.9. Let φ ∈ D(]0, +∞[). For t ∈]0, T ], taking into account Remark 5.5 and Lemma 5.6, we have

R+ v(t, ρ) φ(ρ)dρ = R+ C d ũ(t, ρ) φ(ρ)dρ = R d ũ(t, x d ) φ( x d )dx = R d u(t, x)ϕ(x)dx, where ϕ(x) = φ( x d ). By Remark 2.2(ii), t → u(t, •), t ∈ [0, T ] is weakly continuous in D ′ (R d \{0}) so t → v(t, •), t ∈]0, T ] is weakly continuous in D ′ (]0, +∞[),
and it admits a weakly continuous extension on ]0, +∞[. Therefore

< v(0, •), ϕ > = lim t 0 →0 +∞ 0 v(t 0 , ρ) ϕ(ρ)dρ = lim t 0 →0 +∞ 0 C d ũ(t 0 , ρ) ϕ(ρ)dρ = +∞ 0 C d ũ0 (dρ) ϕ(ρ) = +∞ 0 v 0 (dρ) ϕ(ρ).
This shows the initial condition property, i.e.,

< v(0, •), ϕ >= +∞ 0 v 0 (dρ) ϕ(ρ) = R d u 0 (dx)ϕ(x) .
(5.9)

Now, since u is a solution in the sense of distributions of problem (1.1), we have

R+ v(t, ρ) φ(ρ)dρ = 1 2 t 0 ds R d β(u(s, x))∆ϕ(x)dx + R d u 0 (dx)ϕ(x).
Again, by Remark 5.5, we obtain

R+ v(t, ρ) φ(ρ)dρ = 1 2 t 0 ds R d β( ũ(s, x d ))∆ φ( x d )dx + R d u 0 (dx)ϕ(x).
Expressing the Laplacian in terms of the radius i.e.,

∆ φ( x 2 ) = 2(d 2 -d) x d-2 φ′ ( x d ) + d 2 φ′′ ( x d ) x 2(d-1) , x ∈ R d \ {0},
and using again hyperspherical change of variables, lead to

R+ v(t, ρ) φ(ρ)dρ = C(d 2 -d) t 0 ds R+ r 2d-3 β( ũ(s, r d )) φ′ (r d )dr + Cd 2 2 t 0 ds R+ r 3d-3 β( ũ(s, r d )) φ′′ (r d )dr + R d u 0 (dx)ϕ(x).
Then, setting ρ = r d , using integration by parts and taking (5.9) into account, the result follows. ✷ Proposition 5.10 below is related to the probabilistic representation of (5.8).

Proposition 5.10. Let u 0 ∈ J radially symmetric and ũ0 defined through (5.7). For y = 0, ρ > 0, we set

         Ψ 1 (ρ, y) = (d 2 -d)ρ 1-2 d Φ 2 d C y , Ψ 2 (ρ, y) = dρ 1-1 d Φ d C y , (5.10)
where, C is defined by (5.5).

(i) Suppose that (Z t ) is a non-negative process solving the non-linear SDE defined by

       Z t = Z 0 + t 0 Ψ 2 (Z s , p(s, Z s ))dB s + t 0 Ψ 1 (Z s , p(s, Z s ))ds p(t, •) = Law density of Z t , ∀t ∈]0, T ], Z 0 ∼ C d u 0 .
(5.11)

Then, p is a solution, in the sense of distributions, of the PDE (5.8) with initial condition

       Z t = Z 0 + d t 0 Z 1-1 d s Φ d C p(s, Z s ) dB s + (d 2 -d) t 0 Z 1-2 d s Φ 2 d C p(s, Z s ) ds p(t, •) = Law density of Z t , ∀t ∈]0, T ], Z 0 ∼ C d u 0 .
(5.12)

Proof of Proposition 5.10. (i) Let g ∈ D(]0, +∞[) and Z be a solution of problem (5.11). Itô's formula gives

g(Z t ) = g(Z 0 ) + t 0 g ′ (Z s )Ψ 2 (Z s , p(s, Y s ))dB s + t 0 g ′ (Z s )Ψ 1 (Z s , p(s, Y s ))ds + 1 2 t 0 g ′′ (Z s )Ψ 2 2 (Z s , p(s, Y s ))ds.
Taking the expectation we get

R g(y)p(t, ρ)dρ = R g(ρ) C d u 0 (dρ) + R g ′ (ρ)Ψ 1 (ρ, p(s, ρ))p(s, ρ)dρ + 1 2 t 0 ds R g ′′ (ρ)Ψ 2 2 (ρ, p(s, ρ))p(s, ρ)dρ,
where, p(t, •) is the density law of Z t . This implies the result.

(ii) Let Y be a solution of (1.4), u(t, •) the law density of Y t , t ∈]0, T ] and

S t = ϕ(Y t ) = Y t d , t ∈ [0, T ].
We apply Itô's formula to ϕ(Y t ) to obtain

Y t d = Y 0 d + M t + (d 2 -d) t 0 Y s d-2 Φ 2 ( ũ(s, Y s d ))ds, (5.13)
where

M t = d d i=1 t 0 Y s d-2 Φ( ũ(s, Y s d ))Y i s dW i s ,
is a continuous local (F t )-martingale and (F t ) is the canonical filtration of Y . We observe that

[M ] t = d 2 t 0 Y s 2(d-1) Φ 2 ( ũ(s, Y s d ))ds = d 2 t 0 S 2-2 d s Φ 2 ( ũ(s, S s ))ds.
Enlarging the probability space if necessary, we consider a Brownian motion ϒ independent of F. Let (G t ) be the canonical filtration generated by (F t ) and ϒ. We set

B t = B 1 t + t 0 1 {SsΦ( ũ(s,Ss))=0} dϒ s ,
and B 1 is the (F t ) and (G t )-local martingale defined by

B 1 t = 1 d t 0 1 {SsΦ( ũ(s,Ss))>0} S 1-1 d s Φ( ũ(s, S s )) dM s .
(5.14)

The quadratic variation of B 1 is given by t 0 1 {SsΦ( ũ(s,Ss))>0} ds.

Consequently,

[B] t = [B 1 ] t + t 0 1 {SsΦ( ũ(s,Ss))=0} ds = t. Since B is a (G t )-local martingale, by Lévy's characterization theorem we obtain that B is a (G t )-Brownian motion. Now, for t ∈ [0, T ] M t = t 0 1 {SsΦ( ũ(s,Ss))>0} dM s , (5.15
M t = d t 0 S 1-1 d s Φ( ũ(s, S s ))dB 1 s = d t 0 S 1-1 d s Φ( ũ(s, S s ))dB s , t ∈ [0, T ], since d t 0 S 1-1 d s Φ( ũ(s, S s ))1 {SsΦ( ũ(s,Ss))=0} dϒ s = 0.
On the other hand, ũ(s, •) = d C ν(s, •) by Lemma 5.6, where ν(s, •) is the family of marginal laws of S. This concludes the proof. ✷

A toy model: the heat equation via Bessel processes

Suppose that Y t solves the equation

Y t = Y 0 + W t , t ∈ [0, T ],
where Y 0 is a random variable distributed according to a uniform distribution on the d-dimensional sphere S d-1 centered at 0 with radius ℓ 0 > 0, i.e.,

S d-1 = {x ∈ R d , x = ℓ 0 }. Then, S t = Y t d = R d 2
t , where R is the square of a d-dimensional Bessel process starting at ℓ 2 0 . According to [30, Chapter XI, Section 1, Corollary 1.4], the law density of R t is characterized by

r → q d t (ℓ 2 0 , r) = 1 2t r ℓ 2 0 d-2 4 exp - ℓ 2 0 + r 2t I d 2 -1 ℓ 0 √ r t , t ∈]0, T ],
where I d 2 -1 is the so-called modified Bessel function of the first kind and of index d 2 -1, see e.g. [1, p. 375]. Therefore, the law density of S t at time t, which starts at ℓ d 0 is given by

ν ℓ d 0 (t, ρ) = ρ 2-d 2d tℓ d-2 2 0 d exp - ℓ 2 0 + ρ 2 d 2t I d 2 -1 ℓ 0 ρ 1 d t , t ∈]0, T ], ρ > 0. (5.16)
By Proposition 5.10 (ii), replacing Φ ≡ 1 in (5.10) and (5.11) it follows that the process (S t ) is a solution of the equation

S t = ℓ d 0 + d t 0 S 1-1 d s dB s + (d 2 -d) t 0 S 1-2 d s ds.
Remark 5.12. Take J as the family of all probability measures on R d and U as the family of weakly continuous u : [0, T ] → M(R), t → u(t, •) such that for almost all t ∈]0, T ] u(t, •) admits a density, still denoted by u(t, x), x ∈ R d . For t ∈]0, T ], let H t = tI d . It is well-known that given a probability measure u 0 (dy

) on R d , u characterized by u(t, x) = R d K Ht (x - y)u 0 (dy), t ∈]0, T ], x ∈ R d ,
is the unique solution of the heat equation with initial condition u 0 . By Lemma 5.6 and Remark 5.5, the function

u(t, x) = d C ν ℓ d 0 (t, x d ), t ∈]0, T ], x ∈ R d solves the PDE ∂ t u = 1
2 ∆u in the sense of distributions, with initial condition u(0, •) = u 0 (dx), where u 0 is the distribution of a uniform random variable on the d-sphere S d-1 .

Probabilistic numerical implementation

We adopt here the same notations as in section 5.1. In particular u and u 0 were introduced in the lines before Remark 5.5, the process S was defined in (5.3) with ν as marginal laws. One of our aims is to approximate ũ for d ≥ 2 which coincides up to the constant C d , with ν. We remind in particular that ν is a solution of (5.8) with initial condition ν 0 = C d ũ0 . Our program consists in implementing the one-dimensional probabilistic method developed in [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF]: we introduce in this subsection a stochastic particle algorithm based upon the time discretization of (5.12), which will allow us to simulate the solutions ν of (5.8). From now on we fix n, the number of particles.

Let ε > 0. Similarly to [15, Section 3], we first replace (5.12) by the following mollified version:

                 Z ε t = Z 0 + d t 0 (Z ε s ) 1-1 d + Φ d C (φ ε * p ε )(s, Z ε s ) dB s + (d 2 -d) t 0 (Z ε s ) 1-2 d + Φ 2 d C (φ ε * p ε )(s, Z ε s ) ds p ε (t, •) = Law density of Z ε t , ∀t > 0, Z 0 ∼ C d u 0 , (5.17)
where (x) + = max(x, 0) and φ ε is a mollifier obtained from a fixed probability density function φ by the scaling

φ ε (y) = 1 ε φ y ε , y ∈ R. (5.18)
For the numerical experiments, we assume that φ is a Gaussian probability density function with mean 0 and unit standard deviation. Now, we introduce a particles system given by

Z i,ε,n t = Z i 0 + d t 0 (Z i,ε,n s ) 1-1 d + Φ   d Cn n j=1 φ ε (Z i,ε,n s -Z j,ε,n s )   dB i s (5.19) + (d 2 -d) t 0 (Z i,ε,n s ) 1-2 d + Φ 2   d Cn n j=1 φ ε (Z i,ε,n s -Z j,ε,n s )   ds,
where i = 1, . . . , n.

To simulate a trajectory of each (Z i,ε,n t

), i = 1, . . . , n, we discretize in time: we choose a time step ∆t > 0 and N ∈ N, such that T = N ∆t. We denote by t k = k∆t, the discretization times for k = 0, . . . , N .

The Euler explicit scheme of order one, leads then to the following discrete time system, i.e., for every i = 1, . . . , n:

S i t k+1 = S i t k + d S i t k 1-1 d + Φ d C ν(t k , (S i t k ) + ) N i (0, ∆t) + (d 2 -d) S i t k 1-2 d + Φ d C ν(t k , (S i t k ) + ) ∆t, (5.20)
where N i (0, ∆t), i = 1, . . . , n, are i.i.d Gaussian random variables with variance ∆t. At each time step t k , ν(t k , .) is defined by

ν(t k , y) = 1 n n j=1 φ ε y -(S j t k ) + + 1 n n j=1 φ ε y + (S j t k ) + , y ∈ R + . (5.21)
In fact, ν is a density estimator which is convenient for simulating a density on R + . This is based on a symmetrization technique proposed by Silverman [34, Section 2.10]. Indeed, the classical kernel estimator

ν(t k , y) = 1 n n j=1 φ ε y -S j t k , y ∈ R +
gives an over-smoothed approximation of the target density ν killing the natural discontinuity at zero. We have decided to use Silverman's method because it is simple and easy to implement. However, the derivative of ν with respect to y vanishes when y = 0. This constitutes a limit for the method, even though in our case it gives suitable numerical results. To avoid that feature one could apply the technique explained in the Appendix 9.2. We emphasize that the choice of the smoothing parameter ε, intervening in (5.21), is done according to the bandwidth selection procedure that had been described in [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF]Section 4].

Note that when Φ ≡ 1 and d = 2, previous scheme corresponds to the one of [START_REF] Deelstra | Convergence of discretized stochastic (interest rate) processes with stochastic drift term[END_REF]. Further work on this subject was performed by [START_REF] Diop | Sur la discrétisation et le comportement à petit bruit d'EDS multidimensionnelles dont les coefficients sont à dérivées singulières[END_REF] and more recently by [START_REF] Alfonsi | On the discretization schemes for the CIR (and Bessel squared) processes[END_REF].

The multidimensional probabilistic algorithm

In this section, we want to extend the stochastic particles algorithm introduced in [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF] to the multidimensional case. We will determine a numerical solution of (1.1) by simulating a multidimensional interacting particles system. Again, the solution of the non-linear problem (1.1) is approximated through the smoothing of the empirical measure of the particles. For any n ∈ N, we consider a family of n particles propagating in R d , whose positions at time t ≥ 0 are denoted by Y i,H,n t , i = 1, . . . , n, which evolve according to the system

Y i,H,n t = Y i 0 + t 0 Φ d   1 n n j=1 K H (Y i,H,n s -Y j,H,n s )   dW i s , i = 1, . . . , n, (6.1)
where, (W i ) 1≤i≤n are n d-dimensional standard Brownian motions. (Y i 0 ) 1≤i≤n is a family of independent d-dimensional random variables with law density u 0 and independent of the Brownian motions. K H is the mollifier defined in section 2.

Assuming that the propagation of chaos holds, one expects that the regularized empirical measure 1

n n j=1 K H (• -Y j,H,n t )
approaches the solution u of (1.1).

Remark 6.1. In the case where Φ is Lipschitz, continuously differentiable at least up to order 3, with some further regularity assumptions on u 0 , the authors in [25, Theorem 2.7] established the propagation of chaos. At the best of our knowledge there are no such results when Φ is irregular.

Probabilistic numerical implementation

For fixed T > 0, we choose ∆t > 0 and N ∈ N such that T = N ∆t. We introduce the following numerical Euler scheme which provides us a discrete time approximation of the particles positions (Y i,H,n t ) denoted by (X i t k ),

X i,ℓ t k +1 = X i,ℓ t k + Φ   1 n n j=1 K H (X i t k -X j t k )   N i,ℓ ∆t , 1 ≤ i ≤ n, 1 ≤ ℓ ≤ d, (6.2)
where (N i,ℓ ∆t ) 1≤i≤n,ℓ=1,...,d , is a family of independent Gaussian random variables with mean 0 and variance ∆t.

At each time step t k = k∆t, k = 0, . . . , N , we approximate the function u(t k , •) by the smoothed empirical measure of the particles

u H,n (t k , x) = 1 n n j=1 K H (x -X j t k ), x ∈ R d . (6.3)
From now on, we will suppose that K, as defined in (2.1) is a d-dimensional standard normal density. In particular, we have K(x) = d ℓ=1 φ(x ℓ ). Therefore, the function u H,n (t k , •) becomes the so-called multivariate kernel density estimator of u(t k , •) for every time step t k . The only unknown parameter in (6.3) is the symmetric definite positive d × d matrix H; we refer to it as the bandwidth matrix.

Just as in the univariate case, the optimal choice of H crucially determines the performance of the density estimator u H,n . In fact, a large amount of research was done in this area, we refer to [START_REF] Wand | Kernel smoothing[END_REF] and [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF] for a survey of the subject.

First of all, one has to decide about the particular form of H. A full bandwidth matrix allows for more flexibility; however it also introduces more complexity into the estimator since more parameters have to be selected. A simplification of (6.3) can be obtained by imposing the restriction H ∈ D, where D denotes the subclass of diagonal positive definite d × d matrices. Then, for H ∈ D, we have

H = diag(ε 2 1 , . . . , ε 2 d ), so we have K H (x) = d ℓ=1 φ ε ℓ (x ℓ ).
Besides, a further simplification can be done by considering H = ε 2 I d , where I d is the unit matrix on R d . This restriction has the advantage that one has only to deal with a single smoothing parameter, but the considerable disadvantage is that the amount of smoothing is the same in each coordinate direction. Accordingly, we will suppose from now on that H ∈ D, so that we could have more flexibility to smooth by different amounts in each of the coordinate directions.

It remains to choose the components (ε ℓ ) 1≤ℓ≤d of the bandwidth matrix H itself. For this, we will need some methodology for the mathematical quantification of the performance of the kernel density estimator u H,n . In order to balance between the complexity and the efficiency of the bandwidth selection procedure to be used, we proceed as follows.

The ideal criterion of performance for the estimator q of the density q of some random variable Z, consists in minimizing (asymptotically) the quantity

E q -q 2 L 2 (R d ) , (6.4)
where

q(x) = 1 n n j=1 d ℓ=1 φ ε ℓ x ℓ -Z j,ℓ , x ∈ R d ,
where Z j,. , 1 ≤ j ≤ n are R d -valued random elements, i.i.d according to q.

We have chosen instead to minimize (asymptotically) the quantity

E q -q 2 m , (6.5)
where given f : R d → R, f m is defined as follows

f 2 m = d ℓ=1 R dx ℓ   R d-1 f (x) d k =ℓ dx k   2 . ( 6.6)
In fact, • m is a semi-norm on the linear space f :

R d → R such that f ∈ L 1 (R d ) and f ℓ ∈ L 2 (R), where f ℓ (x ℓ ) = R d-1 f (x) d k =ℓ dx k , ℓ = 1, . . . , d.
It is one generalization of the L 2 -norm for d = 1 to the multidimensional case. • m is indeed a semi-norm since it is non-negative and it verifies the pseudo-homogeneity property and the triangle inequality.

Obviously, the minimal quantity over (ε ℓ ) 1≤ℓ≤d of (6.5) equals the sum over ℓ ∈ {1, . . . , d} of the minimal quantities over each ε ℓ of

MISE {u ε ℓ ,n (t, •)} = E u R u ℓ (t, x) -u ε ℓ ,n (t, x) 2 dx, (6.7)
where, u ε ℓ ,n is the univariate kernel density estimator of the marginal law density u ℓ (of the coordinate X ℓ ). Consequently, each bandwidth ε ℓ , ℓ = 1, . . . , d, will be computed according to the procedure developed by [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] and described with details in [15, Section 6].

The deterministic numerical method

The main aim of our work is to approximate solutions of the d-dimensional non-linear problem given by

∂ t u(t, x) ∈ 1 2 ∆β (u(t, x)) , t ∈ [0, +∞[ , u(0, x) = u 0 (x), x ∈ R d , (7.1)
where u 0 is an integrable function and β is given by (1.3). Despite the fact that, up to now at our knowledge, there are no analytical approaches dealing such issues, we got interested into a recent method, proposed by Cavalli et al. [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF], when β is Lipschitz. Actually, we are heavily inspired by [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF] to implement a deterministic procedure simulating solutions of (7.1) which will be compared to the probabilistic ones.

In our numerical simulations, we will consider the case where d = 2. The operational aspects of that method, in the one-dimensional case, were explained in details in [15, Section 5].

Numerical experiments

The probabilistic and deterministic algorithms were both carried out using Matlab. In order to speed up our probabilistic procedure, we have implemented, using the Matlab Parallel Computing Toolbox (PCT), a GPU version of the kernel density estimator in dimension 1 and 2 of space. Using 10 5 particles, this has 500 times reduced the CPU time on our reference computer. As mentioned in section 7, the deterministic numerical solutions are performed via the method provided in [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF]. In fact, we use the WENO spatial reconstruction of order 5 and a third order explicit Runge-Kutta IMEX scheme for time stepping. From now on, we will denote the related time step by ∆t det and the deterministic numerical solution by ûdet .

The general stochastic particle algorithm for d = 2

We have proceeded to the validation of our algorithm in three main situations: the classical porous media equation, the fast diffusion equation and the Heaviside case.

The porous media equation case

In the case where β(u) = u|u| m-1 , m > 1, we recall that the PDE in (1.1) is nothing else but the classical porous media equation (PME). If u 0 = δ 0 , an exact solution is provided by [START_REF] Barenblatt | On some unsteady motions of a liquid and gas in a porous medium[END_REF](see also [START_REF]The porous medium equation[END_REF]Section 17.5]) known as the density of Barenblatt-Pattle:

E(t, x) = t -α D -κ x 2 t -2β 1 m-1 + , x ∈ R d , t > 0, (8.1)
where 1) and

α = d (m -1)d + 2 , β = α d , κ = m -1 m β, D = κ -d 2 IC 2(1-m) 2+d(m-
I = Γ d 2 Γ m m-1 Γ d 2 + m m-1
, Γ being the usual Gamma function.

We would now compare the exact solution (8.1) to an approximated probabilistic solution. Up to now, we are not able to perform an efficient bandwidth selection procedure in the case when the initial condition of PME is a Dirac probability measure, i.e. the law of a deterministic random variable. Since we are nevertheless interested in exploiting (8.1), we have considered a time translation of the exact solution E defined as

U (t, x) = E(t + 2, x), x ∈ R d , t ∈ [0, T ]. (8.2)
Note that one can immediately deduce from (8.2) that U still solves the PME but now with a smooth initial condition given by

u 0 (x) = E(2, x), x ∈ R d . ( 8

The fast diffusion equation (FDE) case

Now, we suppose that β(u) = u|u| m-1 , m ∈]0, 1[. In that case the PDE in (1.1) corresponds to the so-called fast diffusion equation. Similarly as for the porous media equation, there also exists a Barenblatt type solution for the mentioned β when the initial condition u 0 is a delta Dirac measure at zero. Indeed, it is given by the following expression:

E(t, x) = t -α D + κ x 2 t -2β -1 1-m , x ∈ R d , t > 0, (8.4)
where

α = d (m -1)d + 2 , β = α d , κ = 1 -m m β, D = κ -d 2 IC 2(m-1)
d(1-m)-2 and

I = Γ d 2 Γ 1 1-m -d 2 Γ 1 1-m .
Again, we consider a time shifted version of the explicit solution (8.4) for the numerical experiments. Indeed, we define Figure 4(a)-(d) displays the numerical probabilistic solutions at times t = 0, t = 0.5, t = 1 and t = T = 1.5, respectively. Figure 5(a)-(d) shows the deterministic solution values at t = 0, t = 0.5, t = 1 and t = T = 1.5, respectively. Furthermore, Figure 6 describes the time evolution of the L 1 errors on the time interval [0, 1.5], related to both the probabilistic and deterministic algorithms. Using the analytical expression (8.4) of E(t, •), we get sup

U (t, x) = E(t + 1, x), x ∈ R d , t ∈ [0, T ].
x∈R d E(t, x) ≤ D -1 1-m t -α ,
where the right hand side of previous expression goes to zero as t → +∞. This convergence clearly appears in Figures 4 and5 when d = 2. However, the convergence in L 1 does not hold since for every t > 0, Remark 8.1. In previous cases, we had some exact expressions of the solution of (1.1) at our disposal that we could compare with the approximations issued from the deterministic and probabilistic algorithms. The committed error using the deterministic approach is definitely lower than using the probabilistic one. Below we treat the Heaviside case: by default of exact expressions, the deterministic solutions will be used for evaluating the error related to the probabilistic method.

R d E(t, x)dx = 1 if m > m c ,

The Heaviside case

In this part, we will discuss the numerical experiments for a coefficient β given by (1.3). We recall that in this case we do not know an exact solution for the problem (1.1). Consequently, we will compare our probabilistic solution to the numerical deterministic solution obtained using the method developed in [START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF], see also section 7. In fact, we will simulate both numerical solutions according to several initial data u 0 and with different values of the critical threshold u c . Empirically, after various experiments, similarly to the one-dimensional case investigated in [15, Section 6], it appears that for a fixed threshold u c , the numerical solution approaches some limit function which seems to belong to the "attracting" set

J = {f ∈ L 1 (R 2 )| f (x)dx = 1, 0 ≤ f ≤ u c }; (8.5) in fact J is the closure in L 1 of J 0 = {f : R 2 → R + | β(f ) = 0}.
Again, the following theoretical questions arise.

(1) Does indeed u(t, •) have a limit u ∞ when t → ∞?

(2) If yes, does u ∞ belong to J ? (3) If (2) holds, do we have u(t, •) = u ∞ for t larger than a finite time τ ?

a) Gaussian initial condition

For the mentioned β we consider an initial condition u 0 being a Gaussian density with mean µ and invertible covariance matrix Σ, i.e.,

u 0 (x) = p(x, µ, Σ),
where,

p(x, µ, Σ) = 1 (2π) d 2 |Σ| 1 2 exp - 1 2 (x -µ) t Σ -1 (x -µ) , x ∈ R d . (8.6)
Simulation experiments. Test case 1. We set d = 2, u c = 0.07, µ = (0, 0) and Σ = I 2 , where I 2 is the unit matrix on R 2 . We compute both deterministic and probabilistic solutions over the time-space grid [0, 0.9] × [-4, 4] × [-4, 4] with a uniform space step ∆x = 0.05. For the deterministic approximation we set ∆t det = 2 × 10 -4 while for the probabilistic one we use n = 200000 particles and a time step ∆t = 2 × 10 -4 .

Figures 7,8,9, show the deterministic and probabilistic numerical solutions at times t = 0, t = 0.3 and t = T = 0.9, respectively. Furthermore, Figure 10 describes the time evolution of the L 1 -norm of the difference of the two solutions over the time interval [0, 0.9].

b) Bimodal initial condition

Now, we suppose that the initial condition is a mixture of two Gaussian densities with separated modes, i.e.,

u 0 (x) = 1 2 (p(x, µ 1 , Σ 1 ) + p(x, µ 2 , Σ 2 )) ,
where p is defined in (8.

c) Trimodal initial condition

For the β given by (1.3), we consider an initial condition being a mixture of three Gaussian densities with three modes at some distance from each other, i.e.,

u 0 (x) = 1 3 (p(x, µ 1 , Σ 1 ) + p(x, µ 2 , Σ 2 ) + p(x, µ 3 , Σ 3 )) , x ∈ R d , (8.7)
where p is defined in (1.3). Simulation experiments. We fix again d = 2. For this specific type of initial condition u 0 , we consider two test cases depending on the value taken by the critical threshold u c . We set, for instance,

µ 1 = (-2, 2), µ 2 = (2, -2), µ 3 = (0, 0), Σ 1 = (0.1) 2 I 2 , Σ 2 = (0.2) 2 I 2 and Σ 3 = Σ 2 .
Test case 3. We start with u c = 0.15. We consider a time-space grid [0, 0.4]×[-5, 5]×[-5, 5], with a uniform space step ∆x = 0.05. For the deterministic approximation, we set ∆t det = 2×10 -4 . The probabilistic simulation uses n = 200000 particles and a time step ∆t = 2×10 -4 . Figures 15,16, 17 display both the deterministic and probabilistic numerical solutions at times t = 0, t = 0.14 and t = T = 0.4, respectively. Besides, the time evolution of the L 1 -norm of the difference between the two numerical solutions is depicted in Figure 18.

Test case 4. We choose now as critical value u c = 0.035 and a time-space grid [0, 2] × [-5, 5] × [-5, 5], with a uniform space step ∆x = 0.05. We set ∆t det = 3 × 10 -4 and the probabilistic approximation is performed using n = 200000 particles and a time step ∆t = 4 × 10 -4 . Figures [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF][START_REF] Brezis | Uniqueness of solutions of the initial-value problem for u t -∆ϕ(u) = 0[END_REF][START_REF] Cavalli | High-order relaxation schemes for nonlinear degenerate diffusion problems[END_REF] show the numerical (probabilistic and deterministic) solutions at times t = 0, t = 0.66 and t = T = 2. In addition, Figure 22 describes the L 1 -norm of the difference between the two.

d) Uniform and normal densities mixture initial condition

We proceed again with β defined in (1.3). We are now interested in an initial condition u 0 , being a mixture of a normal and an uniform density, i.e.,

u 0 (x) = 1 2 p(x, µ, Σ) + 1 [0,1]×[-1,0] (x) , x ∈ R 2 ,
where p is defined in (8.6). Simulation experiments. Test case 5. We fix u c = 0.15, µ = (0, -1) and Σ = (0.076) 2 I 2 . We perform both the approximated deterministic and probabilistic solutions on the time-space grid [0, 0.5]×[-3, 3]× [-3, 3], with a space step ∆x = 0.05. We use n = 200000 particles and a time step ∆t = 2 × 10 -4 for the probabilistic simulation. Moreover, we set ∆t det = 2 × 10 -4 . Figures 23,24, 25 illustrate those approximated solutions at times t = 0, t = 0.2 and t = T = 0.6. Furthermore, we compute the L 1 -norm of the difference between the numerical deterministic solution and the probabilistic one. That error is displayed in Figure 26.

Long-time behavior of the solutions

As it was mentioned previously, we are interested in the empirical behavior of solutions to (1.1) in the Heaviside case. For this, we first provide Figure 27, which displays the time evolution of the L 1 -norm of the difference between two successive time evaluations of the numerical solutions. That quantity was computed for both deterministic and probabilistic numerical solutions and in the different test cases 1 to 5. In fact, Figure 27, shows that the numerical solutions approach some limit function û∞ . Indeed, they seem to reach û∞ after a finite time τ .

In addition, according to Figure 28, this suggests the existence of a limit function u ∞ , which of course depends on the initial condition u 0 such that u(t, •) = u ∞ , for t ≥ τ . Moreover, u ∞ is expected to belong to the "attracting set" J , since β( û(t, •)) L 1 (R 2 ) equals zero when t is larger than τ , at least when û is the deterministic numerical solution.

. Now, we inquire about the long time behavior of the probabilistic particle algorithm. In fact, we are interested in the dependence of the error over the time. For this, we have simulated the solution of the PME with n = 200000, T = 50, ∆t = 0.02 and m = 3. Figure 29 displays the time evolution of the L 1 -norm of the error. In particular, Figure 29 shows (in the PME case) that the probabilistic algorithm seems to remain stable for a large time T .

Time evolution of u H,n (t i+1 , •) -u H,n (t i , •) L 1 (R 2 ) (dashed lines) and ûdet (t i+1 , •) -ûdet (t i , •) L 1 (R 2 ) (solid

The radially symmetric case

. Validation on exact solutions in hyperspherical coordinates

We make the same conventions as in section 5. that u 0 is radially symmetric. u is the solution of (1.1) and ũ is given in Remark 5.5. ν constitutes the marginal laws of S which is defined in (5.3) which is a solution of (5.11). ν equals ũ up to a constant, see (5.4).

The Fokker-Planck equation for Bessel processes

When β(u) = u, the PDE in (1.1) corresponds to the heat equation. The linked radial equation (5.8) with Φ = 1 (Bessel equation) and with initial condition ℓ d 0 , admits (5.16) as exact solution. We would like to compare the probabilistic approximation ν ε,n defined as the right-hand side of (5.21) to (5.16). In order to avoid computation difficulties for the bandwidth ε in the case when the initial condition ν is the law of a deterministic random variable, we will consider a time shifted version of (5.16), given by

v(t, ρ) = ρ 2-d 2d (t + t 0 )ℓ d-2 2 0 d exp - ℓ 2 0 + ρ 2 d 2(t + t 0 ) I d 2 -1 ℓ 0 ρ 1 d t + t 0 , t ∈ [0, T ], ρ > 0, t 0 > 0. (8.8)
In fact, with this reposition, v(t, •) still solves (5.8) but now with a smooth initial data v(0, ρ) = ν ℓ d 0 (t 0 , ρ).

Simulation experiments.

We set ℓ 0 = 1 and t 0 = 10 -3 . We compute the probabilistic numerical solutions of (5.8) over a time-space grid [0, 0.01]×]0, L], L > 0, for different values of the dimension d = 2, 5, 10 and with a space step ∆x = 0.01. We use a time step ∆t = 10 -4 and n = 200000 particles. Figures 30,[START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] We point out that the performance of our algorithm is satisfying for all values of d ≥ 2 even though when d = 2 the solution is recurrent. In that case the process often attains zero, which is a non regular point of the diffusion term in (5.12).

The radial transformation of the classical porous media equation

When β(u) = u.|u| m-1 , m > 1 and u 0 = δ 0 , we consider again the explicit Barenblatt type solution of (1.1), denoted by E and given in (8.1). Again, we will shift in time this exact solution, in order to avoid simulation problems in the case of a delta Dirac measure as initial condition. In fact, we set U (t, x) = E(t + 1, x), t ∈ [0, T ], x ∈ R d . Then, U still solves (1.1) for the mentioned β and with u 0 (x) = E(1, x) as initial condition. Since u 0 is radially symmetric, we deduce

U (t, x) = U (t, x d ) = (t + 1) -α D -κ x d 2 d (t + 1) -2β 1 m-1 + , x ∈ R d , t ≥ 0.
Then, using (5.4), we get

ν ex (t, ρ) =          C d (t + 1) -α D -κρ 2 d (t + 1) -2β 1 m-1 if ρ ∈ 0, D κ d 2 (t + 1) α , 0 otherwise.
(8.9)

Simulation experiments. We compute the probabilistic numerical solutions of (5.8) when β(u) = u 3 (radial PME), over the time-space grid [0, 1] × [0, 2.5] for different values of the dimension d = 2, 5, 10 and with a space step ∆x = 0.01. We consider a time step ∆t = 10 -2 and n = 200000 particles. Figures 33,[START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]

Comparison between the the radial stochastic algorithm and the 2-dimensional deterministic approach in the Heaviside case

In this part of the work we exploit, at least empirically, the existing relation between the solutions of the multidimensional problem (1.1) and the solutions of the one-dimensional PDE (5.8), in the case when the initial condition u 0 of (1.1) is a radially symmetric function. For this we will simulate the d-power norm of the solution Y of (1.4) via the one-dimensional non-linear diffusion introduced in (5.12). Let u : [0, T ] × R d → R + be a solution of (1.1) with initial condition u 0 which is radially symmetric. According to Remark 5.5 there exists ũ :

[0, T ]×]0, +∞[→ R + such that u(t, x) = ũ(t, x d ), ∀(t, x) ∈ [0, T ] × R d \ {0}. If Y
u(t, x) = d C ν(t, x d ), ∀(t, x) ∈ [0, T ] × R d \ {0}. (8.10)
Consequently, previous expression allows us to simulate the solution u of (1.1) using the solution ν of (5.8). Indeed, it will be enough to replace in (8.10) ν by its kernel density estimator ν ε,n given in (5.21). Those approximations will be compared to the ones obtained via the ddimensional general stochastic algorithm.

From now on we will fix d = 2 and β given by (1.3).

Test case (A)

Let u 0 be a Gaussian density, i.e.,

u 0 (x) = p(x, µ, Σ), x ∈ R 2 ,
where p is defined in (8.6), µ = (0, 0) and Σ = I 2 . Thus, the initial condition ν 0 of (5.8) is the probability density of an exponential distribution of parameter λ = 1 2 , i.e., ν 0 (ρ) = λ exp(-λρ), ∀ρ ∈ R + .

Simulation experiments. We fix u c = 0.07. We compute the probabilistic solution obtained through the one-dimensional radial algorithm using n = 200000 particles and a time step ∆t = 2 × 10 -4 , and we compare it to the 2-dimensional deterministic approximation presented in section 7. We fix ∆t det = 4 × 10 -4 . We have represented both the 2-dimensional deterministic and probabilistic solutions, on the time-space grid [0, 0.9] × [-4, 4] × [-4, 4], with a uniform space step ∆x = 0.05. Figures 36, 37, 38 illustrate those approximated solutions at times t = 0, t = 0.3 and t = T = 0.9. Furthermore, we compute the L 1 -norm of the difference between the numerical deterministic solution and the probabilistic one. Values of that error are displayed in Figure 39.

Test case (B)

Proceeding with the same β, we consider now an initial condition u 0 given by

u 0 (x) = 1 π 2 π exp(- x 4 2 ), x ∈ R 2 .
The corresponding ũ0 is related via (5.4) to the law density of an absolute value of a standard Gaussian r.v. Indeed, the initial condition ν 0 of (5.8) is a probability density defined by

ν 0 (ρ) = 2 π exp(- ρ 2

Test case (C)

We assume that the initial condition u 0 is defined by

u 0 (x) = 1 2π g( x 2 , m 1 , σ 1 ) + g( x 2 , m 2 , σ 2 ) , x ∈ R 2 ,
where g(ρ, m, σ) = f (ρ, m, σ) + f (-ρ, m, σ), ρ ≥ 0, and f is the density function of a onedimensional Gaussian distribution with mean m and standard deviation σ. Therefore, ν 0 of (5.8) is given by

ν 0 (ρ) = 1 2 (g(ρ, m 1 , σ 1 ) + g(ρ, m 2 , σ 2 )) , ∀ρ ∈ R + .
Simulation experiments. We fix for instance u c = 0.07, m 1 = 0, m 2 = 6, σ 1 = 0.

Remark 8.2. (i)

The probabilistic algorithm can be parallelized on a graphical processor unit (GPU) such that we can speed-up its time machine execution; on the other hand, for the deterministic algorithm, this operation is far from being obvious.

(ii) At this point, even though it provides reliable approximation of the solutions, the implementation of the deterministic algorithm in dimension 2 is not optimal. Indeed, it costs a huge amount of CPU time comparing to the deterministic one-dimensional procedure and to the probabilistic algorithm in dimension 1 and 2. Investigations in that direction will be the object of a future work.

(iii) In general, empirically, the different errors committed by the probabilistic algorithm seem to be reasonable, even though not very good.

a) The general two-dimensional probabilistic algorithm behaves well in the case of an unimodal initial condition. Some difficulties arise in the multimodal case; on the other hand we obtain satisfying results in Figure 26 which represents an evolution in the Heaviside case of a bimodal and irregular initial condition. b) The probabilistic algorithm in the radial case behaves quite well for all d ≥ 2, if the initial condition is unimodal and the coefficient β is smooth. If β is of Heaviside type, the error becomes more important when d = 2. Unfortunately we have no mean to validate the algorithm for larger values than d = 2, in which case we could expect a better performance. In conclusion the use of probabilistic methods in higher dimension is justified. Contrarily to the one-dimensional case, treated in [START_REF] Belaribi | A probabilistic algorithm approximating solutions of a singular PDE of porous media type[END_REF], the probabilistic techniques are much simpler to formulate than the analytical method.

Appendix

Proof of proposition 3.1

We start with some preliminary considerations. Remark 9.1. Referring to Proposition 2 and Remarks following Theorem 1 in [START_REF] Benilan | The continuous dependence on ϕ of solutions of u t -∆ϕ(u) = 0[END_REF], we know the following.

(i) We consider the problem u

-∆w = f, in D ′ (R d ), (9.1)
where, w(

x) ∈ β(u(x)) a.e. x ∈ R d , f ∈ L 1 (R d ). (ii) Let f ∈ L 1 (R d ).
There is a unique u ∈ L 1 (R d ) for which there exists w ∈ L 1 loc (R d ) such that (u, w) solves (9.1).

(iii) It is then possible to define a (multi-valued) operator A := A β : E → E, where D(A) is the set of u ∈ L 1 (R d ) for which there is w ∈ L 1 loc (R d ) such that w(x) ∈ A is a m-accretive operator. (iv) Since A is m-accretive, J λ := (I + λA) -1 is a contraction in the sense that J λ f -

J λ g L 1 (R) ≤ f -g L 1 (R d ) . In addition, J λ is single valued. If f ∈ L ∞ (R d ) then J λ f ∞ ≤ f ∞ and J λ f ∈ L 1 L ∞ .
In particular, for every positive n, J n λ f ∞ ≤ f ∞ .

(v) If f ≥ g, a.e, J λ f ≥ J λ g, a.e.

(vi) Since β(0) = 0, J λ f ≥ 0, if f ≥ 0, a.e. Let u ∈ C + 0 (R d ). For λ > 0, we set u λ := (I + λA) -1 u. Similarly to item 2. of the proof of [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF]Proposition 3.4], it is enough to show that u λ → u, weakly in E.

(

In fact, it is easy to show that D + (A) (and so also its closure D + (A)) is convex. By Satz 6.13 of [START_REF] Alt | Lineare Funktionalanalysis: Eine anwendungsorientierte Einführung[END_REF], D + (A) is weakly sequentially closed. So (9.3) would imply that u ∈ D + (A), which would conclude the proof of this item. It remains to show (9.3). Since (I+λA) -1 is a contraction on E, the sequence (u λ ) is bounded in E. Since u λ ∈ D(A), by definition there exists w λ ∈ L 1 loc (R d ) such that w λ ∈ 1 2 β(u λ (x)) for dx-a.e. x ∈ R d , ∆w λ ∈ L 1 (R d ). Moreover, u = u λ -λ∆w λ .

(9.4)

We recall that (I + λA) -1 (0) = 0 so that, by items (v) and (vi) of Remark 9.1 u λ = (I + λA) -1 (u) ≥ 0, a.e., for every λ > 0. We show the existence of λ 0 > 0 such that the sequence (u λ ) 0<λ<λ 0 is weakly relatively compact. Since the sequence (u λ ) is bounded in E, it is enough to prove the existence of λ 0 > 0 such that (u λ ) 0<λ<λ 0 is Lebesgue uniformly absolutely continuous. Then, by Dunford-Pettis theorem, see [START_REF] Barbu | Analysis and control of nonlinear infinite-dimensional systems[END_REF]Theorem 1.15], the result would follow. Since u is bounded, by Remark 9.1(iv), the sequence ( u λ ∞ ) λ>0 is bounded by u ∞ .

Let F be a Borel subset of R d and K > 1. For λ < λ 0 , we have Consequently, the result will be established if we prove the existence of λ 0 such that the sequence (u λ ) 0<λ<λ 0 is Lebesgue equi-integrable, i. Note that I 2 ≤ ε because of (9.5). Concerning I 1 (λ), we have Since ϕ M → 1 pointwise, as M → +∞, ∆ϕ M ∞ does not depend on M and ∆ϕ M → 0, then equality (9.8) follows by Lebesgue dominated convergence theorem. Together with (9.7), this gives

I 1 (λ) =
I 1 (λ) = R d
(I + λA) -1 u(x)u(x) ϕ(x)dx.

On the other hand

w λ ∞ ≤ β ( u λ ∞ ) -β (-u λ ∞ ) ≤ β ( u ∞ ) -β (-u ∞ ) .
Therefore, (w λ ) is uniformly bounded and so λw λ → 0 in L ∞ , as λ → 0, which implies that λ∆w λ → 0 in the sense of distributions. By (9.4) also u λ → u, in the sense of distributions. Consequently, I 1 (λ) → 0 as λ → 0 and therefore there is λ 0 such that 0 < λ < λ 0 and I 1 (λ) < ε. Consequently, if K > K 0 , 0 < λ < λ 0 and taking into account (9.6), we obtain Consequently, (u λ ) is weakly relatively compact in E and (9.2) is established.

c) The next step consists in showing that a C 0 -solution is a solution in the sense of distributions. We proceed as in item 3. of the proof of Proposition 3.4 in [START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF]. We consider a family of approximating ε-solutions u ε using the fact that A is m-accretive. By item (iv) of Remark 9.1, one can easily show that sup t≤T u ε (t, •) ∞ ≤ u 0 ∞ . (9.9)

We also consider w ε : [0, T ] × R d → R such that w ε (t, x) ∈ β(u ε (t, x))dtdx a.e. Since β is monotone, previous inequality implies that Let K > 0. (9.10) implies the existence of a subsequence (ε n ) such that w εn converges weakly in L 2 ([0, T ] × B(0, K)) to some η ∈ L 2 ([0, T ] × B(0, K)). It remains to see that η(t, x) ∈ β(u(t, x)) a.e. dt ⊗ dx, in order to prove that u solves (2.2). Since u ε → u in L 1 (R d ), uniformly in t, Lebesgue dominated convergence theorem and (9.9) imply that for any K > 0, the sequence (u εn ) converges to u in L 2 ([0, T ]; B(0, K)). The map u → 1 2 β(u) on L 2 ([0, T ] × B(0, K)) is an m-accretive multi-valued map, see [32, p.164 The following reasoning can be used to simulate density laws having a jump at some fixed value a. This happens in particular when we approximate a density with support in [a, +∞[. Without restriction of generality we suppose a = 0 so that the support is R + (instead of the whole line). Let p : R + → R + be a density function, so R+ p(x)dx = 1. Let ε > 0, n ∈ N. Let X 1 , . . . , X n be a sequence of non-negative independent r.v. distributed as p(x)dx. A classical non-parametric estimator is given by p(x) = 1 n n i=1 K ε (x -X i ). (9.12)

|w ε (t, x)| ≤ β ( u 0 ∞ ) -β (-u 0 ∞ ) a.e. ( 9

  )Φ 2 (u(s, y))u(s, y)dy.

  • the Euclidian norm on R d . Let O be an open subset of R d , by D(R d ) (resp. D(O)) we denote the space of infinitely differentiable functions with compact support (resp. compact support included in O) ϕ : R d → R (resp. ϕ : O → R). D ′ (R d ) stands for the dual of D(R d ), i.e., the linear space of Schwartz. If f : R d → R is a bounded function we will denote f ∞ = sup x∈R d |f (x)|. L 1 loc (O) will denote the space of real functions defined on O whose restriction to each compact of O is integrable.

Remark 5 . 8 .

 58 (i) If u 0 : R d → R is integrable and radially symmetric, u 0 (x) = ũ0 ( x d ), for some ũ0 :]0, +∞[→ R. (ii) If ψ ∈ D(]0, +∞[) then ψ(x) := ψ( x d ) belongs to D(R d \{0}). By a change of variables with hyperspherical coordinates, we get R d u 0 (x)ψ(x)dx = +∞ 0 C d ũ0 (r) ψ(r)dr.

CRemark 5 . 11 .

 511 d ũ0 . (ii) If S is defined by (5.3), with marginal laws denoted by ν, then S verifies (5.11) with Z = S and p = ν. For clarification we rewrite explicitly (5.11)

) since t 0 1

 1 {SsΦ( ũ(s,Ss))=0} d[M ] s = 0. By (5.15) and (5.14) we finally get

. 3 )

 3 Simulation experiments. We set d = 2 and m = 3. We compute both the deterministic and probabilistic numerical solutions over the time-space grid [0, 3] × [-2.5, 2.5] × [-2.5, 2.5], with a uniform space step ∆x = 0.0167. We set ∆t det = 7.5 × 10 -4 , while, we use n = 200000 particles and a time step ∆t = 10 -2 for the probabilistic simulation. Figure1(a)-(d) (resp. 2(a)-(d)) displays the numerical probabilistic (resp. deterministic) solutions at times t = 0, t = 1, t = 2 and t = T = 3, respectively. Besides, Figure 3 describes the time evolution of the L 1 probabilistic and deterministic errors on the time interval [0, 3].

Figure 1 .

 1 Figure 1. PME: Probabilistic numerical solution values at t = 0 (a), t = 1 (b), t = 2 (c) and t = T = 3 (d).

Figure 2 .Figure 3 .

 23 Figure 2. PME: Deterministic numerical solution values at t = 0 (a), t = 1 (b), t = 2 (c) and t = 3 (d).

  Obviously, U solves the FDE with u 0 = E(1, •) as initial condition. Simulation experiments. We set d = 2 and m = 1 2 . We consider the time-space grid [0, 1.5] × [-15, 15] × [-15, 15] over which the probabilistic, the deterministic and the exact solutions are computed. We fix ∆t det = 1.5 × 10 -3 . We use a uniform space step ∆x = 0.4. For the probabilistic simulation we set n = 200000 and ∆t = 10 -2 .

Figure 4 .

 4 Figure 4. FDE: Probabilistic numerical solution values at t = 0 (a), t = 0.5 (b), t = 1 (c) and t = 1.5 (d).

Figure 5 .

 5 Figure 5. FDE: Deterministic numerical solution values at t = 0 (a), t = 0.5 (b), t = 1 (c) and t = 1.5 (d).

Figure 6 .

 6 Figure 6. FDE: Evolution of the L 1 probabilistic (solid line) deterministic (dashed line) errors over the time interval [0, 1.5].

2 d

 2 where m c = 0 for d = 1, 2 and m c = d-if d ≥ 3, see [37, Section 5.6].

Figure 7 .

 7 Figure 7. Test case 1: Probabilistic (left) and deterministic (right) solution values at t = 0.

Figure 8 .

 8 Figure 8. Test case 1: Probabilistic (left) and deterministic (right) solution values at t = 0.3.

Figure 9 .Figure 10 .

 910 Figure 9. Test case 1: Probabilistic (left) and deterministic (right) solution values at t = 0.9.

  6). Simulation experiments. Test case 2. We set d = 2, u c = 0.1. We fix µ 1 = (1, 0), µ 2 = (-1, 0), Σ 1 = (0.1)I 2 and Σ 2 = (0.2)I 2 . The deterministic and probabilistic solutions are simulated over the time-space grid [0, 0.8] × [-3.5, 3.5] × [-3.5, 3.5] with a uniform space step ∆x = 0.05. We set ∆t det = 2 × 10 -4 , while we use n = 200000 particles and a time step ∆t = 2 × 10 -4 , for the probabilistic approximation. Figures 11, 12, 13, show the deterministic and probabilistic numerical solutions at times t = 0, t = 0.27 and t = T = 0.8, respectively. Furthermore, Figure 14 displays the time evolution of the L 1 -norm of the difference over the time interval [0, 0.8].

Figure 11 .

 11 Figure 11. Test case 2: Probabilistic (left) and deterministic (right) solution values at t = 0.

Figure 12 .

 12 Figure 12. Test case 2: Probabilistic (left) and deterministic (right) solution values at t = 0.27.

Figure 13 .Figure 14 .

 1314 Figure 13. Test case 2: Probabilistic (left) and deterministic (right) solution values at t = 0.8.

Figure 15 .

 15 Figure 15. Test case 3: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 16 .

 16 Figure 16. Test case 3: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.14.

Figure 17 .Figure 18 .

 1718 Figure 17. Test case 3: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.4.

Figure 19 .

 19 Figure 19. Test case 4: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 20 .

 20 Figure 20. Test case 4: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.66.

Figure 21 .Figure 22 .

 2122 Figure 21. Test case 4: Probabilistic (left) and deterministic (right) numerical solution values at t = 2.

Figure 23 .

 23 Figure 23. Test case 5: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 24 .

 24 Figure 24. Test case 5: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.2.

Figure 25 .Figure 26 .

 2526 Figure 25. Test case 5: Probabilistic (left) and deterministic (right) numerical solution values at t = 0.6.

Figure 28 .

 28 Figure 27.Timeevolution of u H,n (t i+1 , •)u H,n (t i , •) L 1 (R 2 ) (dashed lines) and ûdet (t i+1 , •)ûdet (t i , •) L 1 (R 2 ) (solidlines) for the Test case 1 (a), Test case 2 (b), Test case 3(c), Test case 4 (d) and Test case 5 (e), respectively .

Figure 29 .

 29 Figure 29. Evolution of the L 1 error over the time interval [0, 50], in the PME case.

  ,32 (a)-(c), show the exact and the numerical solutions at times t = 0, t = 0.005 and t = T = 0.01, for d = 2, 5, 10, respectively. The exact solution defined in(8.8), is depicted by solid lines. Besides, Figures30 (e), 31 (e), 32 (e) describe the time evolution of the L 1 error on the interval [0, 0.01], for d = 2, 5, 10, respectively.

Figure 30 .

 30 Figure 30. Bessel equation, d = 2: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.005 (b), t=0.01 (c). The evolution of the L 1 error over the time interval [0, 0.01] (d).

Figure 31 .

 31 Figure 31. Bessel equation, d = 5: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.005 (b), t=0.01 (c). The evolution of the L 1 error over the time interval [0, 0.01] (d).

Figure 32 .

 32 Figure 32. Bessel equation, d = 10: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.005 (b), t=0.01 (c). The evolution of the L 1 error over the time interval [0, 0.01] (d).

  , 35 (a)-(c) display the exact and the numerical solutions at times t = 0, t = 0.5 and t = T = 1 respectively, for d = 2, 5, 10, respectively.The exact solution, defined in (8.9), is depicted by solid lines. Besides, Figures33 (e), 34 (e), 35 (e) describe the evolution of the L 1 norm of the error on the time interval [0, 1], for d = 2, 5, 10, respectively.

Figure 33 .

 33 Figure 33. Radial PME, d = 2: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.5 (b), t=1 (c). The evolution of the L 1 error over the time interval [0, 1] (d).

Figure 34 .

 34 Figure 34. Radial PME, d = 5: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.5 (b), t=1 (c). The evolution of the L 1 error over the time interval [0, 1] (d).

Figure 35 .

 35 Figure 35. Radial PME, d = 10: Exact (solid line) and probabilistic (doted line) solution values at t=0 (a), t=0.5 (b), t=1 (c). The evolution of the L 1 error over the time interval [0, 1] (d).

Figure 36 .

 36 Figure 36. Test case (A): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 37 .

 37 Figure 37. Test case (A): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.3.

Figure 38 .Figure 39 .

 3839 Figure 38. Test case (A): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.9.

2 ),

 2 Figures 40, 41, 42 illustrate those approximated solutions at times t = 0, t = 0.46 and t = T = 1.4. Besides, Figure43displays the L 1 -norm of the difference between the numerical deterministic solution and the probabilistic one.

Figure 40 .

 40 Figure 40. Test case (B): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 41 .

 41 Figure 41. Test case (B): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.46.

Figure 42 .

 42 Figure 42. Test case (B): Probabilistic (left) and deterministic (right) numerical solution values at t = 1.4.

Figure 43 .

 43 Figure 43. Test case (B): Evolution of the L 1 -norm of the difference between the two solutions over the time interval [0, 1.4].

2 and σ 2

 2 = 0.3. We consider n = 200000 particles and a time step ∆t = 2 × 10 -4 , in order to perform the probabilistic numerical solution via the one-dimensional radial algorithm. We compare it to the 2-dimensional deterministic approximation. We set ∆t det = 2 × 10 -4 . We represent both the 2-dimensional deterministic and probabilistic solutions, on the time-space grid [0, 1] × [-4, 4] × [-4, 4], with a uniform space step ∆x = 0.05. Figures44, 45, 46 display those approximated solutions at times t = 0, t = 0.33 and t = T = 1. Besides, Figure47shows the L 1 -norm of the difference between the two solutions.

Figure 44 .

 44 Figure 44. Test case (C): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.

Figure 45 .

 45 Figure 45. Test case (C): Probabilistic (left) and deterministic (right) numerical solution values at t = 0.33.

Figure 46 .Figure 47 .

 4647 Figure 46. Test case (C): Probabilistic (left) and deterministic (right) numerical solution values at t = 1.

1 2 β

 2 (u(x)) for a.e. x ∈ R d and ∆w ∈ L 1 (R d ). For u ∈ D(A), we set A(u) = -12 w | w as in the definition of D(A) . In particular, if β is single valued then Au = -1 2 ∆β(u).

(

  vii) If f ∈ L ∞ , β being monotone, it easily follows that w ∈ L ∞ . In addition, under Assumption B(ℓ), ℓ ≥ 1, then w ∈ L 1 L ∞ .Proof of Proposition 3.1. The proof follows the same lines as[START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF] Proposition 3.4] even though here the assumption on β is slightly more general and d ≥ 1.Problem (3.1) is related with the non-linear evolution equation onE = L 1 (R d ) 0 ∈ u ′ (t) + A(u(t)), 0 ≤ t ≤ T,(9.2)with Au = A β u = -∆(β(u)). We denote by L 1 + (R d ) the space of non-negative functions in L 1 and D + (A) the space of non-negative elements of D(A).a) Suppose thatD + (A) = L 1 + (R d ). By [32, Corollary IV.8.4], there is a unique C 0 -solution of (9.2). It appears as a limit in C([0, T ]; L 1 (R d )) of approximating ε-solutions. For the notion of C 0 -solution u ∈ C([0, T ]; L 1 (R d )) and related ε-solutions, we refer to[START_REF] Showalter | Monotone operators in Banach space and nonlinear partial differential equations[END_REF] Chapter IV.8]. In particular, u ∈ L 1 ([0, T ] × R d ). b) We now prove that D + (A) is dense in L 1 + (R d ).We denote by C + 0 (R d ), the space of nonnegative continuous functions with compact support.Let u ∈ C + 0 (R d ). Since C + 0 (R d ) is dense in L 1 + (R d) it is enough to show that C + 0 (R d ) ⊂ D + (A), where D + (A) is the L 1 -closure of D + (A).

F

  dx|u λ (x)| ≤ F ∩{|x|<K} dx|u λ (x)| + {|x|≥K} dx|u λ (x)| ≤ u ∞ Leb(F ) + sup 0<λ<λ 0 {|x|≥K} dx|u λ (x)|.

  dx|u λ (x)| = 0. Let ε > 0. Since u ∈ L 1 , we consider K 0 > 0 such that {|x|≥K 0 -1} u(x)dx < ε. (9.5)Let ϕ : R d → R be a smooth function with compact support in the closed ball B(0, K 0 ) such that ϕ(x) = 1 if |x| ≤ K 0 -1, 0 ≤ ϕ ≤ 1, elsewhere. For λ > 0, {|x|≥K 0 } (I + λA) -1 u(x) dx = {|x|≥K 0 } (I + λA) -1 u(x)dx ≤ R d (I + λA) -1 u(x)(1ϕ(x))dx ≤ I 1 (λ) + I 2 ,(9.6)whereI 1 (λ) = R d (I + λA) -1 u(x)u(x) (1ϕ(x))dx, I 2 = R d u(x)(1ϕ(x))dx.

  λA) -1 u(x)u(x) ϕ(x)dx. ϕ M : R d → R be a smooth function with compact support on B(0, M + 1), 0 ≤ ϕ M ≤ 1 and ϕ M (x) = 1 if x ∈ B(0, M ). By (9.4) we getR d u λ (x)ϕ M (x)dxλ R d w λ (x)∆ϕ M (x)dx = R d u(x)ϕ M (x)dx.

  {|x|≥K} (I + λA) -1 u(x) dx ≤ {|x|≥K 0 } (I + λA) -1 u(x) dx ≤ 2ε. So, sup λ<λ 0 {|x|≥K} (I + λA) -1 u(x) dx ≤ 2ε.

. 10 )

 10 The proof is then very close, replacing R with R d , to the one of[START_REF] Blanchard | Probabilistic representation for solutions of an irregular porous media type equation[END_REF] Proposition 3.4], until equation(3.8), that saysR d u(t, x)α(x)dx = R d u 0 (x)α(x)dx + lim ε→0 t 0 R d w ε (s, x).α ′′ (x)dxds. (9.11)

9 . 2 .

 92 , Example 2c ]. So it is weakly strongly closed because of [8, p.37, Proposition 1.1 (i) and (ii)]. This conclude the proof of item c). d) The obtained solution belongs to L ∞ ([0, T ] × R d ) by the same arguments as in item 4. of the proof of [19, Proposition 3.4]. i) Finally, uniqueness of the equation in D ′ ([0, T ] × R d ) follows from Theorem 1 and remark 1.20 of [20]. ✷ Elementary tricks to simulate a r.v with values in [a, +∞[, for some a ∈ R

  Proof of Theorem 4.3. (i) Let Y H be a solution of (4.1). As for the proof of Proposition 1.3, Itô's formula implies that the family of marginal laws of Y H t , denoted by t → v H (t, •), is a solution in the sense of distributions of (4.11). v H is weakly continuous because Y H is a continuous process. (ii) Let v = v H be a solution of (4.11). Since Φ is a bounded non-degenerate Borel function, by [27, Section 2.6, Theorem 1] there exists a process Y = Y H being a solution in law of the SDE

.11) Indeed, we state a result which has an interest by itself since it generalizes the one obtained in [15, Theorem 3.2], to the case when d ≥ 1. Theorem 4.3. Suppose that β fulfills Assumption B(ℓ) and it is non-degenerate. Moreover, we assume that Φ is uniformly continuous. (i) If Y H is a solution to (4.1) then, the law of Y H t , t → v H (t, •), is a solution to (4.11). (ii) If v H : [0, T ] → M(R d ) is weakly continuous and is a solution of (4.11) then, the problem (4.1) admits at least one solution in law. Corollary 4.4. The problem (4.11) admits existence of one solution

v = v H : [0, T ] → M(R d ).

Remark 4.5. We do not know any uniqueness results for the problem (4.11).

  solves(1.4), the second item of Proposition 5.10 says that S = Y d verifies the non-linear diffusion equation (5.12). Ifν : [0, T ] × R + → R + such that ν(t, •) is the law of S t , t ∈ [0, T ],then by Lemma 5.6, we have ũ = d C ν and therefore

Acknowledgments. The work of the first and third named authors was supported by the ANR Project MASTERIE 2010 BLAN-0121-01. Part of the work was done during their stay at the Bernoulli Center at the EPFL Lausanne and during the stay of the third named author at Bielefeld University, SFB 701 (Mathematik).

Proposition 9.2. Let X distributed as p(x)dx, U a Bernoulli r.v of parameter 1/2 and independent of X. We set Z = U X + (U -1)X.

Then Z is distributed according to q(x)dx, where

Proof of Proposition 9.2. Let f be a bounded and continuous function on R + . We have

Then, setting y = -x in the second integral of the previous expression, we obtain

So, the result follows. ✷

Let U 1 , . . . , U n be a sample of random variables distributed according to a Bernoulli law of parameter 1/2. A classical estimator of q being the density of the random variable Z defined in the previous proposition, is then given by q

Therefore, as a consequence of Proposition 9.2, one can construct on R + , the following estimator, p, for the density p

Authors information