Optimal convergence rate of the multitype sticky particle approximation of one-dimensional diagonal hyperbolic systems with monotonic initial data

Brenier and Grenier [SIAM J. Numer. Anal., 1998] proved that sticky particle dynamics with a large number of particles allow to approximate the entropy solution to scalar one-dimensional conservation laws with monotonic initial data. In [arXiv:1501.01498], we introduced a multitype version of this dynamics and proved that the associated empirical cumulative distribution functions converge to the viscosity solution, in the sense of Bianchini and Bres-san [Ann. of Math. (2), 2005], of one-dimensional diagonal hyperbolic systems with monotonic initial data of arbitrary finite variation. In the present paper, we analyse the L 1 error of this approximation procedure, by splitting it into the discretisation error of the initial data and the non-entropicity error induced by the evolution of the particle system. We prove that the error at time t is bounded from above by a term of order (1 + t)/n, where n denotes the number of particles, and give an example showing that this rate is optimal. We last analyse the additional error introduced when replacing the multitype sticky particle dynamics by an iterative scheme based on the typewise sticky particle dynamics, and illustrate the convergence of this scheme by numerical simulations.


Introduction
Systems of sticky particles have been known to reproduce the phenomenological behaviour of one-dimensional conservation laws in various physical contexts, in particular in astrophysics or in the study of gas dynamics [18,16].In such systems, finitely many particles evolve on the real line at constant velocity and stick together at collisions, with preservation of mass and momentum but dissipation of energy.The relation between these discrete systems and the equations of continuum physics was formalised by Brenier and Grenier [6], who showed that sticky particle dynamics with a large number of particles allow to approximate the entropy solution to scalar one-dimensional conservation laws with monotonic initial data.We also refer to Bouchut [4], Grenier [11], and E, Rykov and Sinai [8] for previous results in this direction.Based on this idea, we recently introduced a multitype sticky particle dynamics [13] in order to approximate the viscosity solution, in the sense of Bianchini and Bressan [2], of one-dimensional diagonal hyperbolic systems with monotonic initial data of arbitrary finite variation.
These sticky particle dynamics provide natural numerical schemes for the corresponding solutions to scalar conservation laws or diagonal hyperbolic systems.It is thus of interest to control the approximation error due to this procedure.This is the purpose of this article.We shall rely on the remark that sticky particle dynamics generically induce exact weak solutions to the considered equation, but for discrete initial data.Besides, these weak solutions need not satisfy Kružkov's entropy or Bianchini-Bressan's viscosity condition.This leads us to split the total approximation error into a discretisation error of the initial data, and a non-entropicity error induced by the evolution of the particle system.
The discretisation error of the initial data is addressed in Section 2. In particular, if the initial conditions have a compactly supported distributional derivative, this error in L 1 distance for n particles is proved to be bounded from above by a term of order 1/n.The error due to the evolution of the particle system is studied in Section 3 for the case of scalar conservation laws and in Section 4 for the case of diagonal hyperbolic systems.In both cases, this error at time t ≥ 0 is proved to be bounded from above by a term of order t/n.This leads to a global convergence rate of order (1 + t)/n in the number n of particles.The precise statements for scalar conservation laws and diagonal hyperbolic systems are respectively given in Theorems 3.1 and 4.4, which are the main results of this paper.These results are finally illustrated with numerical simulations in Section 5. We emphasise the fact that the sticky particle approach is essentially restricted to the one-dimensional case and mention that (non-)existence and (non-)uniqueness issues related to its multidimensional generalisation were pointed out by Bressan and Nguyen [7].
The remainder of this introduction is dedicated to a detailed presentation of the sticky particle dynamics (SPD) and the multitype sticky particle dynamics (MSPD).
1.1.SPD and scalar conservation laws.This subsection is dedicated to the introduction of the SPD, which allows to approximate the entropy solution to scalar conservation laws in one space dimension.
1.1.1.Scalar conservation laws.Let us consider the scalar conservation law (1.1) for a nonconstant, monotonic and bounded initial condition u 0 .Up to an affine transform of the flux function Λ, one can assume that u 0 is the cumulative distribution function (CDF) of a probability measure m on the real line, which we denote u 0 = H * m where H(x) = 1 {x≥0} is the Heaviside function.The space of probability measures on the real line is denoted by P(R).Then Λ only needs to be defined on the interval [0, 1], and it shall be assumed to have the following regularity.
The following existence and uniqueness result follows from Kružkov's theorem, see [15, Theorem 2.3.5 and Proposition 2.3.6,pp.36-37].In addition, it satisfies the following properties: (i) preservation of total variation: for all t ≥ 0, u(t, •) coincides dx-almost everywhere with the CDF of a probability measure m t on the real line; (ii) finite speed of propagation: if u 0 (a) = 0, then u(t, a − tL C ) = 0 for all t ≥ 0, and if u 0 (b) = 1, then u(t, b + tL C ) = 1 for all t ≥ 0; (iii) stability: if u and v refer to the entropy solutions to the scalar conservation law with respective initial data u 0 and v 0 , then for all t ≥ 0, 1.1.2.Sticky Particle Dynamics.For n ≥ 1, we denote by D n the polyhedron of R n defined by Let x ∈ D n be a vector of initial positions and λ = (λ 1 , . . ., λ n ) ∈ R n a vector of initial velocities.Under the SPD, the particle with index k ∈ {1, . . ., n} has initial position x k , initial velocity λ k and mass 1/n.It evolves at constant velocity on the real line, up to the first collision with another particle.At collisions, the particles stick together and form a cluster: its mass is given by the number of colliding particles over n, and its velocity by the average of the pre-collisional velocities of the particles.More generally, when several clusters collide, they form a single cluster with conservation of total mass and momentum.
For all t ≥ 0, the position of the k-th particle at time t ≥ 0 is denoted by φ k [λ](x; t), and it is easy to check that the process (φ[λ](x; t)) t≥0 defined by φ[λ](x; t) := (φ 1 [λ](x; t), . . ., φ n [λ](x; t)) induces a continuous flow in D n .Its stability with respect to the initial configuration and the vector of initial velocities is detailed in Proposition 1.2 below.Before stating this result, let us define the normalised L 1 distance on D n by Proposition 1.2.[13, Proposition 3.1.9,(i)] Let x, y ∈ D n and λ, µ ∈ R n .For all 0 ≤ s ≤ t, Given an initial configuration x ∈ D n , we define the empirical distribution of the SPD at time t ≥ 0 by and the associated empirical CDF by Given n ≥ 1 and x ∈ D n , it is easily checked that the empirical CDF u n [x] satisfies the properties (i), (ii) and (iii) of Theorem 1.1.The preservation of the total variation is obvious, and the finite speed of propagation is just the transcription of the fact that the modulus of the initial velocities λ k is bounded by L C , uniformly with respect to n.Finally, the L 1 stability follows from (1.2) with µ = λ.
It can also be shown that u n [x] is a weak solution to the scalar conservation law (1.1) with discrete initial condition [13,Proposition 4.2.1].However for fixed n, it does not necessarily satisfy the entropy condition of Theorem 1.1.This property is recovered when taking the limit of an infinite number of particles, as is expressed by the following result, the proof of which is originally due to Brenier and Grenier [6], see also [12] and [13,Lemma 8.2.3] for appropriate generalisations.
Theorem 1.3.Let Λ satisfy Assumption (C) and let m ∈ P(R).Let (x(n)) n≥1 be a sequence of initial configurations such that, for all n ≥ 1, x(n) ∈ D n and the empirical distribution converges weakly to m.For all t ≥ 0, the empirical distribution µ t [x(n)] converges weakly to the probability measure m t ∈ P(R) such that u(t, x) := H * m t (x) is the unique entropy solution of the scalar conservation law (1.1) with initial condition u 0 = H * m.Equivalently, for all t ≥ 0, the empirical CDF u n [x(n)](t, •) converges dx-almost everywhere to u(t, •).
Using Theorem 1.3 to pass to the limit n → +∞ in (1.2), the stability inequality of (iii) can be extended in order to take the dependence of the entropy solution on the flux function into account.This is done in the next proposition, which is of independent interest, and the proof of which is postponed to Appendix A. Proposition 1.4.Let Λ, M : [0, 1] → R satisfying Assumption (C), and u 0 , v 0 be CDFs on the real line.Denote λ := Λ and µ := M , and call u and v the entropy solutions of the scalar conservation law with respective flux function Λ and M, and respective initial data u 0 and v 0 .Then, for all 0 ≤ s ≤ t, 1.1.4.Rate of convergence.Given a sequence of initial configurations (x(n)) n≥1 satisfying the assumptions of Theorem 1.3, our first purpose in this article is to estimate the error when approximating u with u n [x(n)].On account of the stability property stated in Theorem 1.1, a fairly natural distance to measure this error is the L 1 distance on R. Indeed, this stability property allows us to write, for all t ≥ 0, where we have introduced the entropy solution The two terms in the right-hand side above are of a very different nature, and can be estimated separately: the first term corresponds to the discretisation error of the measure m, while the second term only measures the non-entropicity induced by the evolution of the particle system for a given initial condition The discretisation error of the measure m is addressed in Section 2. There, we use the fact that, given two probability measures m and m on the real line, the L 1 distance between H * m and H * m is the Wasserstein distance of order 1 between m and m , defined by where the infimum runs over all the probability measures m ∈ P(R 2 ) such that, for all Borel sets A, A ⊂ R, This is due to the fact that, on the real line, the measure where U refers to the Lebesgue measure on [0, 1], realises the infimum in (1.4).In this definition, the pseudo-inverse F −1 of a CDF is defined by We deduce that As a consequence, the first term in the right-hand side of (1.3) rewrites Precise bounds on W 1 (m, µ n ) in terms of n for the optimal discretisation of m are derived in Lemma 2.1.They depend heavily on the tail of m.In particular, an important remark to be done at this point is the following.Assume that m has an infinite first order moment.Then, since by (1.4) and the triangle inequality, any approximation of m by a measure µ n with finite first order moment necessarily satisfies W 1 (m, µ n ) = +∞.As a consequence, there is no purpose in trying to compute a rate of convergence in this case.Therefore, although our results hold true without any assumption on m, they only have a nontrivial content when m has a finite first order moment.The non-entropicity error is then addressed in Section 3, where given arbitrary n ≥ 1 and x ∈ D n , an estimation is first derived on the L 1 distance between u n [x](t, •) and u ∞ [x](t, •) in Proposition 3.2.This result holds under the following strengthening of Assumption (C).
Combining the results of Section 2 with Proposition 3.2 yields complete rates of convergence of the SPD, as is stated in Theorem 3.1.
(C) For all γ ∈ {1, . . ., d}, the function λ γ is continuous on [0, 1] d .Under Assumption (C), we denote We also require the system to be uniformly strictly hyperbolic, in the sense of the following assumption.
(USH) There exists L USH ∈ (0, +∞) such that ∀γ ∈ {1, . . ., d − 1}, inf On account of the fact that the system (1. where ω γ γ:k (x) denotes the (scaled) rank of the particle γ : k within the system of type γ , formally defined by where particles of different types sharing the same location are counted according to the postcollisional rank imposed by Assumption (USH).
The resulting dynamics in D d n is the MSPD, denoted by (Φ(x; t)) t≥0 .More details on its construction are given in [13,Subsection 3.2], where it is proved that it defines a continuous flow.Its stability with respect to the initial configuration is described by the next result, which forms the core of the article [13].
Proposition 1.5.[13, Theorem 2.5.2]Under Assumptions (LC) and (USH), there exists L 1 ∈ [1, +∞) depending only on d and the ratio L LC /L USH such that, for all n ≥ 1, for all x, y ∈ D In contrast with Proposition 1.2, we note that no stability result with respect to the characteristic fields λ 1 , . . ., λ d is available.1.2.3.Approximation of the diagonal hyperbolic system.Similarly to the scalar case, we define the empirical distribution of the system of type γ in the MSPD at time t ≥ 0 as the probability measure on the real line given by and the vector of empirical CDFs is a weak solution to the system (1.7).Given a sequence of initial configurations (x(n)) n≥1 such that, for all n ≥ 1, x ∈ D d n and, for all γ ∈ {1, . . ., d}, the empirical distribution converges weakly to some m γ ∈ P(R), one could by analogy with Theorem 1.3 expect u n [x(n)] to converge to a weak solution of the system (1.7) with initial data u γ 0 = H * m γ , satisfying some specific entropy-like condition making it unique and physically meaningful.Although it is true that, up to extracting a subsequence, u n [x(n)] actually converges to a weak solution of the system [13, Theorem 2.4.5],following Bianchini and Bressan's construction we were only able in [13] to identify the limit if it satisfies some semigroup and stability estimate with respect to the initial data.For this purpose, we introduced [13, Definition 2.6.4] a discretisation operator For all γ ∈ {1, . . ., d}, the empirical measure converges weakly to m γ , see [13, Lemma 8.1.5],and we obtained the following convergence theorem for the MSPD, which follows from the more general statement of [13, Theorem 2.6.5].
Theorem 1.6.Under Assumptions (LC) and For all t ≥ 0, for all γ ∈ {1, . . ., d}, the empirical distribution µ γ t [χ n m] converges weakly to the probability measure is the unique semigroup solution to the diagonal hyperbolic system (1.7) with initial data u γ 0 = H * m γ .Equivalently, for all t ≥ 0, for all γ ∈ {1, . . ., d}, the empirical CDF Each coordinate u γ of the semigroup solution u satisfies the preservation of total variation and finite speed of propagation properties stated in Theorem 1.1.The stability property writes as follows: for all m, m ∈ P(R) d , the semigroup solutions u and v to the system (1.7) with respective initial data u 0 , v 0 defined by The notion of semigroup solution is adapted from [2] and detailed in [13, Definition 8.2.5].
In Corollary 4.6, we shall generalise this result and prove the convergence of the empirical cumulative distribution functions of the MSPD to the semigroup solution for a large class of sequences of initial configurations (x(n)) n≥1 approximating a vector of initial measures m = (m 1 , . . ., m d ) such that, for all γ ∈ {1, . . ., d}, m γ has a finite first-order moment.
1.2.4.Rate of convergence and approximation by the iterated TSPD.The second purpose of this paper is to supplement Theorem 1.6 with a rate of convergence of u n [χ n m](t, •), or more generally u n [x(n)](t, •) for a sequence of initial configurations (x(n)) n≥1 approximating m, to u(t, •).As in the scalar case, we split the approximation error by writing where u ∞ [x(n)] is the semigroup solution, in the sense of Theorem 1.6, of the system (1.7 In the last line above, we used the stability estimate of Theorem 1.6Â on semigroup solutions.
The estimation of the discretisation error of the initial data follows from the analysis of the scalar case carried out in Section 2. We shall use the distance The analysis of the error due to the evolution of the MSPD turns out to be more delicate than in the scalar case, and we shall resort to an approximation of this dynamics by the following scheme.Fix a time step ∆ > 0 and define the process ( Φ∆ (x; t)) t≥0 in D d n by letting, on each time interval [(L − 1)∆, L∆], L ≥ 1, the particles evolve according to the TSPD with initial velocity vector λ( Φ∆ (x; (L − 1)∆)).This amounts to neglecting the collisions between clusters of different types on [(L − 1)∆, L∆] and only updating the velocities with respect to the new ordering of the system at the end of each such interval.From a computational point of view, this approximated dynamics is expected to be easier to simulate than the MSPD, as one does not have to keep track of all the collisions.The approximation error induced by this scheme is addressed in Section 4, where we also use it as a theoretical tool to obtain rates of convergence in Theorem 4.4.We finally present a numerical implementation of this scheme in Section 5.

W 1 discretisation error of initial conditions
In this section, we want to approximate the probability measure m on the real line, with CDF F = H * m, by the empirical measure On account of (1.6), one has where we recall the definition (1.5) of the pseudo-inverse F −1 and complement it with the convention The choice (2.1) the median of the image of the uniform law on [(k − 1)/n, k/n] by F −1 , for all k ∈ {1, . . ., n}, minimises this Wasserstein distance.Let us now analyse this optimal choice.First, When that is to say m has a finite first order moment, one deduces that where the right-hand side goes to 0 when n → +∞ by Lebesgue's theorem.With the discussion of the case where the first order moment of m is infinite at the end of §1.1.4,we deduce the first statement in the next lemma.
Lemma 2.1.Let µ n be defined as in the discussion above.The behaviour of W 1 (m, µ n ) depends on the tail of m as is described below.
(i) One has lim n→+∞ W 1 (m, µ n ) = 0 if and only if m has a finite first order moment.
(ii) For all n ≥ 1, where the right-hand side may be infinite.
Proof.The assertion (ii) follows from the comparison of W 1 (m, µ n ) with the expected Wasserstein distance between m and the empirical measure 1 n n k=1 δ X k of independent random variables X 1 , . . ., X n with identical distribution m, a detailed study of which was carried out by Bobkov and Ledoux [3].By the optimality of the choice of µ n , we first have According to the proof of [3, Theorem 3.2], the right-hand side is not greater than If m([a, b]) = 1, then the mass of the measure ν on (0, 1) such that ν((0, u) Let us finally assume that m(dx Since m is the image of the Lebesgue measure on [0, 1] by F −1 , one has for all u ∈ (1/(2n), 1).The continuity of translations in L 1 ensures that it is enough to check that .
This follows from the weak convergence of 2 n k=1 1 {u∈((2k−1)/(2n),k/n)} du to 1 {u∈(0,1)} du and the density of continuous and bounded functions in Similarly, the condition is an assumption on the decay of the tails of m.Since F On the other hand, if for some α > 2, We finally discuss some other cases where W 1 (m, µ n ) can be estimated.If m has a positive density f on (0, ∞), one has On the other hand, if u → 1/f (F −1 (u)) is nondecreasing on (u, 1) with u ∈ (0, 1), then We apply these estimates on the following two examples.

Rate of convergence of the SPD to scalar conservation laws
Under the assumptions of Theorem 1.3, the purpose of this section is to estimate the L 1 distance between the entropy solution u of the scalar conservation law (1.1) with initial condition u 0 = H * m for some m ∈ P(R), and the empirical CDF u n [x(n)] associated with the SPD started at some configuration x(n) = (x 1 (n), . . ., x n (n)) ∈ D n , over finite time horizons.Theorem 3.1.Let Λ satisfy Assumption (LC).Then for all n ≥ 1, for all Proof.Following the approach described in the introduction of the article, we first write, for all t ≥ 0, where u ∞ [x(n)] refers to the entropy solution of the scalar conservation law (1.1) The L 1 stability property of Theorem 1.1 for the scalar conservation law (1.1)yields Proposition 3.2.Under Assumption (LC), let n ≥ 1 and let x ∈ D n .For all t ≥ 0, Proof.Instead of comparing u n [x] directly with the entropy solution u ∞ [x], we first compare it with the empirical CDF of the SPD with 2n particles, started in the configuration where the positions in x are duplicated, and then iterate this comparison.More precisely, for all n ≥ 1, let us define the operator for all k ∈ {1, . . ., n}, and denote by xm ∈ D 2 m n the m-th iteration of this operator.We first prove that, for all n ≥ 1, for all x ∈ D n , In this purpose, we interpret the function u n [x] as the empirical CDF of the SPD with 2n particles, with initial configuration x, and with modified velocity coefficients λ1 , . . ., λ2n defined by for all k ∈ {1, . . ., n}.Then (1.2) yields We now fix n ≥ 1 and use (3.3) to write that, for all M ≥ 1, To complete the proof, we have to check that To this aim, we note that, for all M ≥ 1, the empirical measure µ 0 [x M ] associated with xM does not depend on M and coincides with µ 0 [x].As a consequence, by Theorem 1.3, Using the finite speed of propagation for both •) has a compact support which does not depend on M .As a consequence, (3.4) follows from Lebesgue's theorem and the proof is completed.
for all n ≥ 1, x ∈ D n and t ≥ 0. Note that this remark holds true even without Assumption (LC).
On the contrary, one can construct examples of convex flux functions for which the bound of Proposition 3.2 provides the right order of magnitude for u Let us indeed consider the Burgers equation Λ(u) = u 2 /2 with the initial condition m = δ 0 , the Dirac mass at 0. Fix n ≥ 1 and let x = (0, . . ., 0) ∈ D n .On the one hand, the entropy solution u On the other hand, the particles drift away from each other in the SPD, so that which is of the same order as the bound given in Proposition 3.2.where we recall that Φ(x; t) refers to the MSPD while Φ[ λ(x)](x; t) is the TSPD with initial velocity vector λ(x) given by (1.8).By the flow property of the MSPD, the L-th iteration Φ L ∆ (x) is nothing but Φ(x; L∆).On the other hand, the L-th iteration Φ∆ (x) is easier to compute as one does not have to take interactions between particles of different types into account.The purpose of this subsection is to estimate the error of the approximation of the MSPD by this iterated TSPD scheme.To this aim, given x ∈ D d n , we first extend the iterated TSPD into a continuous process ( Φ∆ (x; t)) t≥0 in D d n by interpolating between the points of the uniform grid with step ∆ thanks to the TSPD.More precisely, for all L ≥ 1, for all t ∈ [(L − 1)∆, L∆], we define Our purpose is to prove the following result.The value of C is given in (4.6) below.Of course, for t restricted to the uniform grid (L∆) L≥0 , the estimate (4.1) implies sup We begin by proving this restricted estimation before deducing (4.1).We first use the flow property of both the TSPD and the MSPD [13, Proposition 3.2.8] to write (4.2) where the last inequality follows from the discrete stability estimate of Proposition 1.5.The next lemma allows to estimate each term of the sum in the right-hand side.We recall from [13] that, for all y ∈ D d n , the set R(y) contains the pairs of particles (α : i, β : j) ∈ (P d n ) 2 such that α < β and y α i < y β j , so that in the MSPD started at y, these particles undergo a collision at a finite and positive time τ coll α:i,β:j (y).We also define the first collision time in the MSPD started at y by t * (y) := min{τ coll α:i,β:j (y), (α : i, β : j) ∈ R(y)} ∈ (0, +∞], where we take the convention that t * (y) = +∞ if R(y) is empty.Proof of Lemma 4.2.Let t 0 := 0 < t 1 < • • • < t R < ∆ =: t R+1 refer to the successive instants of collisions (i.e.t * (y), t * (Φ(y; t * (y))), etc.) in the MSPD started at y, on the time interval [0, ∆].For all r ∈ {0, . . ., R + 1}, let us denote y r := Φ(y; t r ) and ỹr := Φ[ λ(y)](y; t r ), so that y 0 = ỹ0 = y while y R+1 = Φ ∆ (y) and ỹR+1 = Φ∆ (y).Besides, the definition of the MSPD and flow property of both the TSPD and the MSPD yield, for all r ∈ {1, . . ., R + 1}, As a consequence, where the last line is obtained by applying the stability property of the SPD (1.2) typewise.Using the triangle inequality ) and summing by parts, we obtain Our purpose is to exhibit a bound of this order (with respect to n) for the whole sum, uniformly in L. We shall rely on the following remark, which is a consequence of Assumption (USH) and of the boundedness of the velocities.We deduce from the first part of the remark that, for all ∈ {1, . . ., L}, R( Φ −1 ∆ (x)) ⊂ R(x), which allows us to rewrite Now for all (α : i, β : j) ∈ R(x), let α:i,β:j be the lowest index in {1, . . ., L} such that (α : i, β : If there is no such index then the contribution of (α : i, β : j) in the sum above is null.Otherwise, the second part of the remark implies that after at most m := 2L C /L USH iterations, the pair (α : i, β : j) can no longer belong to R( Φ α:i,β:j −1+m ∆ (x)).We deduce that (4.5) and therefore conclude that Injecting this bound in (4.4), we complete the proof of (4.1) with the constant (4.6) Now, for t ∈ [(L − 1)∆, L∆], by reasoning like in the derivation of (4.2), one obtains (α:i,β:j)∈R(y) (α:i,β:j)∈R(y) Then, for all t ≥ 0, where the constant C is given in (4.6).
Proof.The estimations involving ũn,∆ are an immediate consequence of Proposition 4.1 and the estimations involving u n [x(n)].To prove those estimations, following the approach described in the introduction of the article, we first write, for all t ≥ 0, where u ∞ [x(n)] refers to semigroup solution to the diagonal hyperbolic system with initial condition ).The L 1 stability property of Theorem 1.6 yields which, according to (iii) in Lemma 2.1, is smaller than 1 for all γ, and x(n) is given by (4.7).The following proposition allows to control the second term in the right-hand side of (4.8).Proposition 4.5.Under Assumptions (LC) and (USH), let n ≥ 1 and x ∈ D d n .For all t ≥ 0, Proof.For all n ≥ 1, let us extend the duplication operator introduced in the proof of Proposition 3.2 by setting where xγ 2k−1 = xγ 2k = x γ k for all k ∈ {1, . . ., n} and γ ∈ {1, . . ., d}.Notice that (4.9) Like in the proof of Proposition 3.2, we are going to estimate The direct analysis of the MSPD being delicate, we introduce a ficticious step ∆ > 0 to transform this analysis into the comparison between the TSPD with n particles and 2n duplicated particles on each time-step.Let t > 0 and L = t/∆ .One has Hence, by the triangle inequality and the stability of the MSPD, (4.10) . Now, by the triangle inequality and (4.9) The second term in the right-hand side is a comparison at time ∆ between the TSPD with n particles starting from Φ −1 ∆ (x) and the TSPD with 2n particles starting from the duplicated vector Φ −1 ∆ (x).Reasoning like in the derivation of (3.3), we bound it from above by ∆dL LC /(2n).By Lemma 4.2, where N ∆ has been defined in (4.3).Plugging these estimations in (4.10) and dealing in the same way with the last term in the right-hand side of (4.10), we deduce that Clearly, for all α : i, β : We now use the same arguments as in the derivation of (4.5) to obtain that, for all α : i, β : To this aim we only have to replace Remark 4.3 with the following: let y ∈ D d n and α : i, β : j ∈ P d n such that α < β, and denote z := Φ ∆ (y).
• If (α : i, β : j) ∈ R(y), then for all (i , j The remainder of the argument is the same and leads to (4.11).As a consequence, Taking the limit ∆ → 0, we deduce that With this estimation replacing (3.3), we shall conclude like in the proof of Proposition 3.2.We therefore need to show that To this aim, we denote m n := µ 0 [x] and recall the definition (1.9) of the discretisation operator to write Using the same arguments as in the proof of Proposition 3.2 but with Theorem 1.6 replacing Theorem 1.3, in particular the finite speed of propagation for both the MSPD and u ∞ [x], we obtain that the second term in the right-hand side vanishes.On the other hand, the discrete stability estimate of Proposition 1.5 yields As a consequence, the corresponding empirical CDF converge dx-almost everywhere.Since the support of each measure by (1.10), and thereby completes the proof.
As an immediate corollary of Theorem 4.4, we obtain a convergence result for the MSPD to the semigroup solution of the system (1.7) which holds under less stringent conditions on the sequence of initial configurations.Corollary 4.6.Under Assumptions (LC) and (USH), let m ∈ P(R) d and let (x(n)) n≥1 be such that for all n ≥ 1, x(n) ∈ D d n and (4.12) For all t ≥ 0, the empirical CDF u n [x(n)](t, •) converges in L 1 (R) d to the semigroup solution u(t, •) of the system (1.7) with initial data u 0 = (H * m 1 , . . ., H * m d ).
Under the condition (4.12), this corollary extends both [13, Theorem 2.4.5] and [13, Theorem 2.6.5]: in the former it allows to identify the limit, in the latter it relaxes the assumption that x(n) = χ n m.We however underline that a necessary condition for (4.12) to hold is that, for all γ ∈ {1, . . ., d}, m γ have a finite first-order moment.

Numerical implementation
5.1.Scalar conservation laws.This subsection is dedicated to the numerical implementation of the SPD in order to approximate the entropy solution to the scalar conservation law (1.1).The algorithm we use is described in §5.1.1.Then two case studies are presented: §5.1.2addresses the Burgers equation (5.1) with the CDF of the two-atom measure as initial datum, while §5.1.3deals with the conservation law with concave flux function with the CDF of the two-sided exponential measure (or Laplace distribution) (5.4) m(dx) = 1 2 exp(−|x|)dx as initial datum.
5.1.1.Numerical computation of the SPD.Given t ≥ 0, a vector of (ranked) initial positions x = (x 1 , . . ., x n ) ∈ D n , and a vector λ = (λ 1 , . . ., λ n ) ∈ R n of initial velocities, we use a remark due to Brenier and Grenier [6,Section 4] to devise an algorithm computing the vector φ[λ](x; t) of positions at time t in O(n) elementary operations, without following the detailed trajectory of each particle on the time interval [0, t].
Let us define the free transport flow ψ[λ](x; t) ∈ R n by, for all k ∈ {1, . . ., n}, and introduce the functions P t and Q t on [0, 1] such that P t (0) = Q t (0) = 0 and, for all k ∈ {1, . . ., n}, with linear interpolation on [(k − 1)/n, k/n].Brenier and Grenier pointed out that P t is the convex hull of Q t .Thus, our algorithm consists in computing the vector {Q t (k/n), k = 1, . . ., n} first, which obviously requires O(n) operations, and then deducing P t from the algorithm described by the pseudo-code below, which follows Andrew's monotone chain algorithm [1], based on the Graham scan [10].Note that the parallelisation of the computation of the vector {Q t (k/n), k = 1, . . ., n} is straightforward, while a parallelisable method to determine the convex hull of a sorted point set in the plane was devised in [9].
The input is the vector Q of size n+1, with elements Q(k) = Q t (k/n) indexed by k ∈ {0, . . ., n}.The algorithm constructs a list L of integers k (ranked in decreasing order) containing the successive points at which P t (k/n) = Q t (k/n).Once L is computed, P t is reconstructed by linear interpolation, which finally yields φ[λ](x; t).The length of the list L is denoted by |L| and its elements are indexed starting from 1.It is assumed that its first and second elements can be accessed and removed in constant time.remove the first element from L end while insert k at the beginning of L end for It is easily checked, by induction on k ∈ {1, . . ., n}, that after the (k − 1)-th iteration of the for loop, L contains the indices of the points at which Q t coincides with the convex hull of {Q t (k /n), k = 0, . . ., k}.Indeed, the while loop consists in removing from L the points L(1) at which the piecewise linear function interpolating between the values of Q t at the points L(2), L(1) and k is concave, see Figure 1.From the fact that the list L is browsed at each iteration, one could think that the algorithm uses O(n 2 ) operations.This is actually not the case, as elements of L for which the test in the while loop returns true are removed from L and will therefore not appear again in the next iterations.As a consequence, the actual complexity of this algorithm is O(n).
Notice that an algorithm computing the explicit space-time points of collisions in the SPD would also take O(n) operations, as there are at most n − 1 such points.As a consequence, the method presented here has the same computational efficiency as the explicit simulation of the SPD.However, the Brenier-Grenier trick allows for a significative simplification of the implementation.5.1.2.Burgers equation with two-atom initial measure.We consider the Burgers equation (5.1) with the CDF of the two-atom measure (5.2) as initial datum.In the SPD, two fans of particles are created, respectively originating from the points −1 and 1, see Figure 2 (a).These fans correspond to the fact that the entropy solution is the superposition of two rarefaction waves, respectively located at time t on [−1, −1 + t/2] and [1 + t/2, 1 + t], see Figures 2 (b) and (c).
The L 1 error between the particle approximation and the solution is plotted as a function of n, for several given terminal times t, on Figure 2 (d).In accordance with Proposition 3.2 and Remark 3.3, and since there is no discretisation error of the initial condition here (for even n), it is observed that this error is of the order of magnitude 1/n.(5.3) and the CDF of the two-sided exponential measure (5.4) as initial datum.As is made clear on Figure 3 (a), the particles progressively aggregate at 0. It results in the formation of a shock wave in the solution, see Figures 3 (b) and (c).The L 1 error is displayed on Figure 3 (d) and exhibits the following behaviour: given t ≥ 0, there exists a critical number of particles such that: • below this number, the error does not vary with n, • above this number, the error decreases when n increases at the same rate as for the discretisation of the initial measure.
This behaviour is explained by the fact that, for n small, all the particles have arrived at 0 at time t, so that the approximate solution is the Heaviside function whatever n.But as soon as n is large enough to allow some particles to have an initial position far enough from 0 so that they have not reached 0 at time t yet, then the contribution of these particles in the approximate solution allows the latter to fit better the part of u outside of the shock wave, at the same rate as for the initial discretisation since the shape of u outside of the shock wave is merely an affine transformation of the initial profile.Following the conclusions of Section 2, this discretisation error is of the order log(n)/n.
Of course, the larger t is, the larger the magnitude of the shock wave is, therefore the better the Heaviside function approximates u and the more particles it takes to reach the critical number, which explains the ordering of the different curves on the picture.Figure 3. Numerical results for a concave flux function with two-sided exponential initial measure: (a) trajectories of 30 particles in the space-time plane, (b) value of the solution u(t, x) in the space-time plane, (c) profile of the solution u(t, x) at successive times t = 0, 0.5, 1, . . ., 5, (d) logarithmic plot of the L 1 error between the approximation obtained with 2 p particles and the solution, as a function of p = 1, . . ., 12.The different lines correspond to different values of t, namely t = 0, 1, . . ., 17, and the higher curves correspond to the smaller times.5.2.Diagonal hyperbolic systems.We now turn to the numerical resolution of the diagonal hyperbolic system (1.7) thanks to the MSPD.
A first method to simulate the trajectory of the MSPD obviously consists in computing the exact space-time position of each collision (between particles of the same type, or between particles of different types).The number of such collisions is at most of order n 2 : indeed, because of Assumption (USH), there are at most n 2 d(d − 1)/2 collisions between particles of different types, and whenever particles of the same type collide, the space-time point of the next collision with a particle of another type is the same for all the particles in the current cluster.As a consequence, an algorithm computing the exact trajectory of each particle is expected to perform O(n 2 ) elementary operations.We however believe that such an algorithm with optimal complexity would require a rather technical implementation.
We therefore suggest to use a second method, which consists in approximating the MSPD by the iterated TSPD scheme described in Subsection 4.1.Thanks to the Brenier-Grenier algorithm introduced in the scalar case, each iteration of the SPD is easily implemented and requires O(n) elementary operations.Then we shall show below that updating the velocities after each step also requires O(n) operations.As a consequence, computing the iterated TSPD on L steps requires O(nL) elementary operations.On the other hand, the error between the solution to the system (1.7) and the approximated solution provided by the iterated TSPD scheme was proved in Theorem 4.4 to be of order O(t/n + ∆) at time t.For the terms t/n and ∆ to be of the same order, one therefore has to run the iterated TSPD scheme on L t/∆ n iterations, which leads to a total number of elementary operations in O(n 2 ).As a conclusion, this method has the same cost as the exact simulation of the MSPD, but it seems easier to implement.
It follows from this discussion that to reach an approximation error of order at time t ≥ 0, the iterated TSPD scheme requires O(t 2 / 2 ) elementary operations.In comparison, standard upwind schemes for the hyperbolic system (1.7), with a time step ∆t and a mesh size ∆x satisfying the CFL condition L C ∆t ≤ ∆x, are generally expected first-order accurate [14], so that the approximation error at time t writes C(t)∆x, with a constant C(t) that depends neither on ∆t nor on ∆x, and grows at least linearly with t.Besides, at each iteration of such a scheme, O(1/∆x) elementary operations are necessary to compute the values of the fluxes and of the solution on the grid, so that after L t/∆t iterations, O(t/(∆t∆x)) elementary operations have been performed.As a consequence, the minimal number of elementary operations to reach a precision of order at time t is obtained when the CFL condition is saturated, and it is in O(tC(t) 2 / 2 ), which has the same dependence on as the iterated TSPD scheme.5.2.1.Description of the algorithm.In order to simulate the MSPD, we use the approximation of the latter by the TSPD on small time steps ∆, as is described in Subsection 4.1.Given x ∈ D d n , we thus compute ΦL ∆ (x) instead of Φ(x; L∆).To this aim, we use an elementary iterative algorithm which will therefore perform L steps.At each iteration ∈ {1, . . ., L}, it is necessary to: (i) compute the vector of velocities for the TSPD started from Φ −1 ∆ (x), (ii) compute the evolution of each subsystem of particles according to the TSPD.
Of course, the second step uses the algorithm described in Subsection 5.1 and therefore makes O(nd) elementary operations.The first step is realised by the following pseudo-code, the input of which is an array x of size d × n, such that x(gamma,k) contains the initial position x γ k of the particle γ : k.We recall that, for fixed γ, x γ k ≤ x γ k+1 for all k ∈ {1, . . ., n − 1}.current_indices = vector [0, ..., 0] of size d while min(current_indices) < n gamma = max( argmin( x(g,current_indices(g)) ; g such that current_indices(g)<n ) ) k = current_indices(gamma) current_indices(gamma) = k+1 velocity(gamma,k) = lambda(gamma,current_indices) end while In this pseudo-code, lambda(gamma,[k_1, ..., k_d]) returns the velocity of the particle γ : k, so that at the end of the algorithm, the vector velocity(gamma,:) contains the initial velocities of the particles of type gamma for the SPD.
There are nd iterations of the while loop, and to select gamma it is necessary to scan the vector current_indices, which costs d operations.As a consequence, the computation of the velocities requires O(nd 2 ) operations.Since d is a physical parameter, we only consider the complexity with respect to the numerical parameter n and therefore conclude that the computation of ΦL ∆ (x) is made in O(nL) operations.

5.2.2.
Case study: p-system.The p-system is a simple model for isentropic gas dynamics in one space dimension, where u is the specific volume of the gas and v is its velocity.The function p(u) determines the pressure in terms of the specific volume, and must generically satisfy p (u) < 0 for all u ≥ 0. In the sequel, we shall furthermore assume that there exists ν > 0 such that (5.5) ν u=0 −p (u)du = 1, which is the appropriate condition to develop our probabilistic approach, and in which case the specific volume will take its values in [0, ν].
Let us define the Riemann invariants w − and w + for this 2 × 2 system [15] by where, for all u ∈ [0, ν], g(u) := The assumption (5.5) ensures that, if w − , w + ∈ [0, 1], then (w + − w − )/2 belongs to the image of g, and it is immediately checked that u and v are recovered from the formulas (5.6) Furthermore, the Riemann invariants satisfy ∂ t w ± = ±c(u)∂ x w ± , which rewrites under the form of the 2 × 2 diagonal system so that, for initial conditions w − 0 and w + 0 given by the CDFs of probability measures, the system satisfies Assumption (USH) with constant L USH = 2 .
We now present numerical approximations of u and v for the choice of pressure function where ν > 0 is a given reference specific volume and κ > 0 is a dimensionless shape parameter.This choice implies We first address the case where the initial conditions w − 0 and w + 0 are the respective CDFs of the shifted two-sided exponential distributions m − and m + defined by m ± (dx) := 1 2 exp(−|x ± x 0 |), for some x 0 ≥ 0. In this case, w + 0 (x) ≥ w − 0 (x) for all x ∈ R, which implies that w + (t, •) − w − (t, •) remains nonnegative at all times t ≥ 0. This is indeed easily checked at the level of the MSPD, a trajectory of which is plotted on Figure 4: if w − 0 and w + 0 are the empirical CDFs of two vectors x ± = (x ± 1 , . . ., x ± n ) ∈ D n , then w + 0 ≥ w − 0 if and only if, for all k ∈ {1, . . ., n}, x + k ≤ x − k .By (5.8), for all t ≥ 0 we have Φ + k (x; t) ≤ x + k − t, Φ − k (x; t) ≥ x − k + t, with the obvious notation x = (x − , x + ) ∈ D 2 n , so that the corresponding empirical CDF satisfies w + (t, •) ≥ w − (t, •).That this inequality still holds true in the limit n → +∞ can be checked using the notion of trajectories introduced in [13, Section 5], see in particular [13,Corollary 5.1.2].
One can observe on Figure 4 that particles of the same type never collide with each other.This is due to the fact that, for fixed w − ∈ [0, 1], the mapping w + → λ + (w − , w + ) is increasing on [w − , 1].As a consequence, two consecutive particles of type + with no particle of type − between them can only have velocities taking them away from each other.The same phenomenon occurs for particles of type −.However, collisions between particles of different types modify the velocities of these particles.Thus, the Riemann invariants w − and w + , respectively plotted on Figures 5  (a  As a sticky particle dynamics where particles never stick together may seem a little disappointing, we now choose initial conditions that do not satisfy the condition that w + 0 (x) ≥ w − 0 (x) for all x ∈ R. To this aim, we still assume w − 0 to be given by the CDF of m − (dx) = exp(−|x − x 0 |)/2 with x 0 ≤ 0, and take w + 0 (x) = 1 {x≥0} .Particles of both type can now aggregate into clusters, as is depicted on Figure 6.But on account of Assumption (USH), after a finite time (that generally depends on the number of particles), the property that w + ≥ w − is recovered and the particles start drifting away from each other again.This is also observed on Figure 6, where two clusters blow up under the effect of a collision.As a result, the functions w − , w + , u and v exhibit shocks on short times, and are essentially given by interacting rarefaction waves on longer times, see Figure 7.

1. 2 . 1 .
Diagonal hyperbolic systems.Let us fix an integer d ≥ 2, and consider the diagonal hyperbolic system (1.7) which, according to (iii) in Lemma 2.1, is smaller than (b − a)/(2n) if m([a, b]) = 1 and x(n) is given by (3.1).The following proposition allows to control the second term in the right-hand side of (3.2).

Remark 3 . 3 .
If the flux function Λ is concave, then u n [x] is actually the entropy solution to the scalar conservation law (1.1) with discrete initial condition[5, Lemma 3.3], so that

Figure 1 .
Figure 1.An iteration of the algorithm computing L. The left-hand figure displays the composition of L at the beginning of the (k-1)-th iteration.On the central figure, the point k-1 is removed from L. The right-hand figure displays the composition of L at the end of the k-th iteration.

Figure 2 .
Figure 2. Numerical results for the Burgers equation with two-atom initial measure: (a) trajectories of 50 particles in the space-time plane, (b) value of the solution u(t, x) in the space-time plane, (c) profile of the solution u(t, x) at successive times t = 0, 8, 16, . . ., 80 and (d) logarithmic plot of the L 1 error between the approximation obtained with 2 p particles and the solution, as a function of p = 1, . . ., 9. The different lines correspond to different values of t, namely t = 10, 20, . . ., 50, with the higher curves corresponding to the larger times.The slope of each line is −0.693 − log 2, which expresses the order 1/n of the error.

Figure 4 .
Figure 4. Trajectory of the MSPD (obtained with the iterated TSPD scheme with ∆ = 0.03) with 20 particles per type associated with the p-system for w + 0 ≥ w − 0 .Blue rays correspond to particles of type −, red rays correspond to particles of type +.Here x 0 = 0.1 and ν = 0.5, κ = 5.

Figure 6 .
Figure 6.Trajectory of the MSPD (obtained with the iterated TSPD scheme with ∆ = 0.02) with 20 particles per type associated with the p-system for initial conditions with shocks.Blue rays correspond to particles of type −, red rays correspond to particles of type +.Red particles first remain aggregated into a single cluster up to the collision with the median blue particle, which makes it blow up then.Here x 0 = −4 and ν = 0.5, κ = 5.

Figure 7 .
Figure 7. Representation in the space-time plane of the Riemann invariants w − (a) and w + (b), of the specific volume u (c) and of the velocity v (d) for initial conditions with shocks.The simulation uses the iterated TSPD algorithm with parameters n = 200 and ∆ = 0.02, and ν = 0.5, κ = 5.
The first step towards the definition of the MSPD is the introduction of the Typewise Sticky Particle Dynamics (TSPD).Given an array of initial velocities λ = (λ ) .(y) = t r appears twice in the sum in the right-hand side above, we deduce that [13,P d n λγ k (y r ) − λγ k (y r−1 ) ≤ L LC n γ:k∈P d n γ =γ n k =1 1 {τ coll α:i,β:j (y)=tr} ,where (α : i, β : j) in the last term refers to (γ : k, γ : k ) if γ < γ and to (γ : k , γ : k) if γ > γ .We note that the argument here is similar to the one employed in the proof of Fact 2 in[13, Lemma 7.2.4].Since each pair (α : i, β : j) ∈ R(y) such that τ coll α:i,β:j