Numerical schemes for the aggregation equation with pointy potentials

The aggregation equation is a nonlocal and nonlinear conservation law commonly used to describe the collective motion of individuals interacting together. When interacting potentials are pointy, it is now well established that solutions may blow up in ﬁnite time but global in time weak measure valued solutions exist. In this paper we focus on the convergence of particle schemes and ﬁnite volume schemes towards these weak measure valued solutions of the aggregation equation


Introduction
This work is devoted to the numerical approximation of measure valued solutions to the so-called aggregation equation in R d . This is a nonlinear and nonlocal conservation law that is commonly used to model the dynamics of (a density of) individuals interacting together through an interaction potential. Denoting by W the interaction potential, its gradient ∇ x W (x − y) measures the relative force exerted by a unit mass located at a point y onto a unit mass located at a point x, and the aggregation equation reads We complement this system with the initial condition ρ(0, ·) = ρ ini .
When considering fully attractive potentials, ∇W may have some discontinuities. This is the case, for instance, for potentials of the type W (x) = w(|x|) whose gradient may have a singularity in 0 even for smooth w. It is known that in this situation, L p weak solutions may blow up in finite time, in the sense that the L p norm, for p > 1, blows up in finite time (see e.g. [4][5][6]). Since the aggregation equation is conservative, a notion of solutions in the sense of measures has been developped, for which global-in-time existence may be obtained. Two different approaches may be found in the literature. On one side, weak measure valued solutions have been defined in [11] thanks to the theory of gradient flows for a Wasserstein metric [2]. On the other side, based on the theory of Filippov [24], weak measure valued solutions have been defined as a pushforward by a flow in [13] (see also [27] for the one dimensional case). In this latter approach, the aggregation equation is seen as a transport equation with a velocity field ∇ x W * ρ which satisfies a one-sided Lipschitz-continuity condition (see (5) below for the definition of one-sided Lipschitz-continuity condition). It has been also proved in [13] that both approaches are equivalent (and there is also an equivalence with entropy solution to conservation law in the one dimensional case as stated in [8,27]).
In this work, we are interested in the numerical treatment of the aggregation equation (1) with pointy potential. Numerical investigations of regular solutions of the aggregation equation, i.e. before the blow-up time, may be found for instance in [10,16]. However, there are no convergence results towards weak measure valued solutions after blow-up time. Based on the approach using Filippov flows, convergence towards weak measure valued solutions has been proven in [28] for one-dimensional numerical solutions constructed by a finite volume scheme. It has been extended in higher dimension for some finite volume schemes on structured meshes [13]. Moreover, a precise error estimate has been obtained recently by some of the authors in [19] showing that the convergence of an upwind numerical scheme is of order 1/2 in Wasserstein distance. The obtention of this convergence order relies on the construction of a stochastic characteristics thanks to a probabilistic interpretation of the numerical scheme in the spirit of [17,18].
One aim of this work is to provide some numerical illustrations of this latter convergence result and to investigate convergence on unstructured meshes. In particular, we will observe that on unstructured meshes, numerical solutions to (1), when computed by some standard finite volume schemes, do not converge towards weak measure valued solutions of (1). Note that this is in contrast with the result proven in [19]; there, it is shown that forward semi-lagrangian schemes (of upwind type) do converge at the order 1/2; of course, the latter schemes are conceptually very different from the finite volume typed schemes that we are handling here. This paper is organized as follows. Section 1 is devoted to a brief summary of the existence and uniqueness result of weak measure valued solutions to the aggregation equation, and recalls and summarizes some material already presented in [13,19]. A consequence of this existence result is the convergence of a particle scheme in which a finite number of particles is considered and its dynamics is discretized thanks to a Euler scheme. A proof of the convergence at the order 1/2 is provided in section 2 (a weak convergence proof, without rate, was already provided in [11,13]). In section 3, we provide some illustration of the convergence order for some finite volume schemes in dimension 1. The results show a convergence order that is better than expected when the potential is pointy. Finally, section 4 provides finite volume numerical results in higher dimension on unstructured meshes.

Existence result of weak measure valued solutions
In this section, we summarize the existence and uniqueness result for the aggregation equation (1) that may be found in [13] (see also [19], and [29], in which some slight generalizations are proposed). We first start by setting our assumptions on the interaction potential W and some useful notations.

Assumptions and notations
We may always assume, up to a rescaling, that the total mass of the system is 1. Then the inital data is assumed to be a probability measure. Moreover, we may assume, up to a translation, that the center of mass is 0, i.e. R d xρ ini (dx) = 0. The interaction potential W : R d → R is assumed to satisfy the following properties: Such a potential will be referred to as a pointy potential. Typical examples are the fully attractive potentials W (x) = 1 − e −|x| , which is −1-convex, and W (x) = |x|, which is 0-convex. Notice that the Lipschitz-continuity of the potential allows to bound the velocity field: there exists a nonnegative constant a ∞ such that for all Remark also that (A3) implies that λ ≤ 0 in (A1). In the following, we may avoid assumption (A3) to allow λ > 0 and get better estimates, but in this case we make the further assumption that the initial datum of the Cauchy problem has a compact support. In this case, a ∞ will be a bound for |∇W (x)| on the support of ρ ini . We denote by M b (R d ) the space of Borel signed measures whose total variation is finite. For ρ a measure in M b (R d ) and Z a measurable map, we denote Z # ρ the pushforward measure of ρ by Z; it satisfies, for any continuous function φ, We call P(R d ) the subset of M b (R d ) made of probability measures. For p ≥ 1, we define the space of probability measures with finite pth order moment by Here and in the following, | · | stands for the Euclidean norm in R d , and ·, · for the Euclidean inner product. The space P p (R d ) is equipped with the Wasserstein distance d p defined by (see e.g. [2,33,38]) where Γ(µ, ν) is the set of measures on R d × R d with marginals µ and ν, i.e.
By a weak compactness argument, we know that the infimum in the definition of d p is actually a minimum. A measure that realizes the minimum in the definition (3) of d p is called an optimal plan, the set of which is denoted by Γ 0 (µ, ν). Then, for all γ 0 ∈ Γ 0 (µ, ν), we have

Filippov flow for linear transport equation
Let us first consider the linear conservative transport equation We assume that the velocity field has a weak regularity, more precisely for α ∈ L 1 loc ([0, +∞)). It has been established in [24] that a Filippov characteristic flow could be defined. For s ≥ 0 and x ∈ R d , a Filippov characteristic starting from x at time s is defined as a continuous function In this definition, Convess(E) denotes the essential convex hull of the set E. We remind briefly the definition for the sake of completeness (see [1,24] for more details). We denote by Conv(E) the classical convex hull of E, i.e., the smallest closed convex set containing E. Given the vector field b(t, ·) : R d → R d , its essential convex hull at point x is defined as where N 0 is the set of zero Lebesgue measure sets and B(x, r) is the ball of center x and radius r > 0. Moreover, we have the semi-group property: for any t, τ, s ∈ [0, +∞) such that t ≥ τ ≥ s and x ∈ R d , From now on, we will make use of the notation Y (t, x) = Y (t; 0, x), for a Filippov characteristic. Since characteristics may be constructed, then solutions to the conservative transport equation (4) with a given bounded and one-sided Lipschitz-continuous velocity field could be defined as the pushforward of the initial condition by the Filippov characteristic flow, i.e. ρ(t) = Y (t) # ρ 0 . The well-posedness of this solution has been established in [32]. Moreover stability properties have been recently established in [7].

Existence and uniqueness of a Filippov flow
We are now in position to state an existence result of a Filippov flow for the aggregation equation (1). For ρ ∈ C([0, T ], P 2 (R d )), we define the velocity field a ρ by where we have used the notation Due to the λ-convexity of W , see (A2), we deduce that, for all x, y in R d \ {0}, Moreover, since W is even, ∇W is odd and by taking y = −x in (8), we deduce that inequality (8) still holds for ∇W , even when x or y vanishes: This latter inequality provides an one-sided Lipschitz-continuity (OSL) estimate for the velocity field a ρ defined in (7), i.e. we have As a consequence, we are in the framework to construct a Filippov flow for this velocity field. Such construction has been established in [13]. More precisely the statement reads: [13, Theorem 2.5 and 2.9] [19, Theorem 2.1] Let W satisfy assumptions (A0)-(A3) and let ρ ini be given in P 2 (R d ). Then, there exists a unique solution ρ ∈ C([0, +∞); P 2 (R d )) satisfying, in the sense of distributions, the aggregation equation where a ρ is defined by (7). This solution may be represented as the family of pushforward measures (ρ(t) : is the unique Filippov characteristic flow associated to the velocity field a ρ . Moreover, the flow Z ρ is Lipschitz-continuous and we have At last, if ρ and ρ are the respective solutions of (10) with ρ ini and ρ ini, as initial conditions in P 2 (R d ), then This existence result also holds true when we assume that W only satisfies assumptions (A0)-(A2) and ρ ini ∈ P 2 (R d ) is compactly supported (see [19,

Upwind discretization
We denote by ∆t the time step and consider a Cartesian grid with step ∆x i in the ith direction, i = 1, . . . , d; we then let ∆x := max i ∆x i . We also introduce the following notations. For a multi-index  = (a 1 , . . . , a d ).
For a given nonnegative measure ρ ini ∈ P 2 (R d ), we put, for any J ∈ Z d , The notation (a) + = max{0, a} stands for the positive part of the real a and respectively (a) − = max{0, −a} for the negative part. The macroscopic velocity is defined by Since W is even, we also have: Based on a probabilistic approach, it has been proved in [19] that the above upwind scheme converges at order 1/2 in the Wasserstein distance d 2 . More precisely the convergence result reads: with a ∞ as in (2). For ρ ini ∈ P 2 (R d ), let ρ = (ρ(t)) t≥0 be the unique measure solution to the aggregation equation with initial datum ρ ini , as given by Theorem 1.1. Define ((ρ n J ) J∈Z d ) n∈N as in (12)-(13)- (14) and let Then, there exists a nonnegative constant C, only depending on λ, a ∞ and d, such that, for all n ∈ N * , Note that one has J ρ 0 J = ρ(R d ) = 1 because the datum is a probability measure, but of course all the results of this paper remain true when the initial mass is any finite positive real number (in that case it suffices to rescale the datum by dividing it by the initial mass, then perform the computation, and at last multiply the result by the initial mass).

Particle scheme
In this section, we address another scheme than the upwind discretization introduced in Subsection 1.4. Convergence at order 1/2 is proven in Theorem 2.1 below. Numerical examples are given in Section 3.

Definition of the scheme
Let ρ ini ∈ P 2 (R d ) be the initial datum for (1), and (ρ 0 J ) J∈Z d be a discrete version of ρ 0 defined for every where M J = Π d i=1 [(J i − 1/2)∆x, (J i + 1/2)∆x[, so that we can define an approximation for ρ, with x J = ∆x × J. The solution to (1) with this discrete datum, denoted by ρ ∆x (t), satisfies ρ ∆x (t) = Z(t)#ρ 0 ∆x , where Z is the associated characteristic flow. From the stability property (see Theorem 1.1), one has d 2 (ρ(t), ρ ∆x (t)) ≤ e |λ|t d 2 (ρ ini , ρ 0 ∆x ), and, as d 2 (ρ ini , ρ 0 ∆x ) ≤ C∆x for a certain constant C, In the flow Z(t), each characteristic starts from a point x J , J ∈ Z d , and transports a particle of mass ρ 0 J . Denoting by (Y J (t)) t≥0 the trajectory of the characteristic starting from x J , the Filippov flow reduces to the sticky particles dynamics (see [9]) This is a direct consequence of the fact that (10) is satisfied in the sense of distributions. Since we are in the aggregative case, two particles may collide in finite time. If for instance the particle I (that is to say, the one associated with ρ 0 I ) collides with the particle K, then they form a bigger particle with mass m I + m K and the dynamics continues with one particle less.
In order to approximate (1) in a discrete way, we propose the explicit Euler type scheme for (20): for every J ∈ Z d . The corresponding approximation of ρ(t n ) is then defined as Theorem 2.1. Let ρ ini ∈ P 2 (R d ). Let ρ ∆x be defined as in (22) thanks to the particle scheme (21), (18).
Proof. We first notice the fact that the function | ∇W | is bounded on the support of the solution by a ∞ in both cases (i) and (ii). Indeed, it is obvious for (i), and for (ii) we use Lemma 2.2 below which states that the numerical solution is compactly supported. Thanks to (19) and the triangle inequality, in order to prove the estimate, we only have to estimate the term d 2 (ρ n ∆x , Z(t n , ·)#ρ 0 ∆x ) = d 2 (X n #ρ 0 ∆x , Z(t n , ·)#ρ 0 ∆x ).
We know that

Thus (25) becomes
which writes, taking ε = ∆t, As assumed in the theorem, we choose ∆t small enough to ensure ∆t < 1 and ∆t < 1/(2λ). In this way, one has 2λ(1 − ∆t)∆t < 1, and we get by induction In the above proof we have used the following Lemma to guarantee the boundedness of the velocity: We consider also, up to a translation, that the center of mass is 0, i.e.

Numerical examples
Implementing the scheme (21), we can observe the trajectories of the particles for typical potentials. Figures  1 and 2 show the trajectories of 20 particles where the initial datum is a standard gaussian law truncated and renormalized over the interval [−3, 3]. On figure 1, we observe the positions of the particles with respect to time, with the smooth potential W (x) = x 2 (this illustrates the fact that no collisions occur). Figure 2 show the results with the pointy potential W (x) = |x|. We can note that despite the previously established convergence of the scheme, the numerical aggregation does not correspond to a "proper" gluing of the particles. Indeed, the time discretization makes possible the crossing of the different trajectories. This explains the oscillations with small amplitude observed in this latter figure.

Numerical illustration of the convergence order results in dimension 1
In this section, we numerically illustrate the convergence order results obtained in Theorems 2.1 and 1.2 on one dimensional test cases.
where m is computed in such a way that the integral of the initial condition over (−1, 1) is equal to 1. Actually, the numerical initial condition is computed in a slightly different way than announced in (18): one takes where ρ ini is identified to its density function, and, moreover, the normalizing coefficient m is computed in such a way that J ρ 0 J = 1. Following the particle scheme (21), we run simulations with 2 k particles, for k in the set {6, . . . , 12} and for two different potentials: the pointy potential W (x) = |x|, and the smooth potential W (x) = |x| 2 /2.
The errors in the Wasserstein distances d 1 and d 2 are computed at time T = 1 for each solution, relatively to the next one. That is, the error for a solution with 2 k particles at time T = 1 is computed using the solution with 2 k+1 particles.
In a one dimensional setting, the Wasserstein distances have an explicit expression. Here we choose to compute it with the help of the Python package POT [39], since we are also using it later in two dimensions.
On Figure 3, one can see the convergence curves obtained for the two potentials.

Wasserstein errors
Space step size Convergence of the particle scheme for a pointy potential W1 W2 Slope 1 Figure 3. Order of convergence of the particle scheme for the smooth potential W (x) = x 2 /2 (left) and for the pointy potential W (x) = |x| (right). Theorem 2.1 states a convergence order for the particle scheme of 1/2 in time and 1 in space for the Wasserstein distance d 2 . For the simulations, we choose a time step of the same order as the space step, meaning that an order 1/2 is expected numerically. The results clearly show a better order. For both potentials and both distances, the order is 1. This suggests that our estimate in Theorem 2.1 might not be optimal, at least for smooth initial conditions.

Finite volume scheme: illustration of Theorem 1.2
The one dimensional problem presented here consists in two Dirac masses of weight 0.5 at a distance of 1 from each other at initial time. We are again considering the pointy and the smooth potential of Section 3.1: (1) for the pointy potential W (x) = |x|, given the initial condition, the exact solution can be computed. The two Dirac masses move toward each other, both at constant velocity 0.5, to merge at time 1 and form a single Dirac masse with a weight of 1. (2) for the smooth potential W (x) = x 2 2 , the exact solution is also known. The two Dirac masses move toward each other but will never merge. The distance between the two Dirac masses is indeed e −t in this case. We consider the domain [−0.75, 0.75] and set the two initial masses at the points −0.5 and 0.5. Following the notation of Section 1.4, the values ρ 0 J are thus all zeros, except for the two cells J − and J + containing the points −0.5 and 0.5 respectively. For these two cells, we set ρ 0 J− = ρ 0 J+ = 0.5/∆x. Using the expression of the velocity given by (14) and the upwind scheme of (13), we compute both the Wasserstein distances d 1 and d 2 between the numerical and the exact solution at time t = 0.75.
To study the convergence order of the upwind scheme, we run, for the two potentials, several simulations with a number of discretization points of 2 k for k in the set {6, . . . , 12}. Because the exact solution is known, we compute the Wasserstein errors with respect to this exact solution. Figure 4 shows the results obtained for the two Wasserstein distances d 1 and d 2 .  Figure 4. Order of convergence for the smooth potential W (x) = x 2 /2 (left) and for the pointy potential W (x) = |x| (right).
The numerical order of convergence in the case of the smooth potential is rather clearly 0.5 for both d 1 and d 2 . Recalling estimate (17) of Theorem 1.2, this is the expected result in distance d 2 . However, for the non-smooth pointy potential, it seems that we recover a first order of convergence. This difference can be explained by how two close Dirac masses interact with each other through the potential. Indeed, the velocity depends on the gradient of the potential. For the smooth potential, when two masses are close to each other, this gradient is close to zero. Compared to the numerical diffusion of the scheme, the agregation phenomenon is less important and the order 0.5 is obtained. On the contrary, the pointy potential has a constant gradient that does not depend on the distance between the two masses. Even at close range, the agregation phenomenon is important in this case. Moreover, it is acting as the opposite of the numerical diffusion and counter balanced it. This seems to lead to a better order of convergence (what would need a rigorous proof). Note that the link with the Burgers equation with decreasing initial datum when the potential is |x|, made in [27], makes it clear that this superconvergence phenomenon is linked with the same superconvergence for the Burgers equation with standard finite volume schemes, see [20], and also [34] for the viscous Burgers equation.

Numerical results on unstructured meshes in dimension 2
It has already been proved and numerically checked that the upwind scheme (13) converges on Cartesian meshes for these aggregation equations (see [18]). It is easy to extend this result to another diffusive scheme such as the Rusanov scheme. We are interested here in their behaviors when used with non-Cartesian meshes.
We first recall briefly the definition of the upwind and finite volume schemes on general meshes, for more details we refer to [25]. Let us consider an admissible finite volume mesh of R d denoted T (see Definition 6.1 in [25]). For an element K ∈ T , we denote |K| the measure of K and V(K) the set of its neighbours. For L ∈ V(K), we denote L ∩ K the common interface between K and L and by ν KL the unit normal oriented from K to L.
We introduce the following explicit in time finite volume scheme This scheme is initiated by the condition ρ 0 K = 1 |K| K ρ ini (dx). The function g allows to define the numerical flux. In this part we will consider the two following numerical methods : • Rusanov g(ρ n K , ρ n L , ν KL ) = 1 2 (ρ n K a n K · ν KL + ρ n L a n L · ν KL + a ∞ (ρ n L − ρ n K )) .
• Upwind g(ρ n K , ρ n L , ν KL ) = (a n K · ν KL ) + ρ n K − (a n L · ν KL ) − ρ n L . We use a two dimensional version of the toy problem of Section 3 to evaluate the numerical order of convergence. The triangular mesh is made of one line of squares, all cut in two along the same diagonal. Thus, the domain is the set [−0.75, 0.75] × [−∆x, ∆x], where ∆x is the mesh step size. The initial condition is the same as in Section 3. We identify the two cells J − and J + whose center is the closest to the points (−0.5, 0) and (0.5, 0) respectively. For these two cells, we set ρ 0 J− = ρ 0 J+ = 0.5/A, where A is the area of a cell. Everywhere else, the initial condition is zero. The exact solution remains the same as this problem is essentially a one dimensional problem.
Considering again the pointy and the smooth potential, we run the simulations for several mesh step size and compute both the Wasserstein errors d 1 and d 2 . In order to compare these results, we run the same simulations with the Rusanov scheme. Figure 5 shows the results obtained with the two schemes. Convergence of upwind and Rusanov schemes for a pointy potential upwind -W1 upwind -W2 rusanov -W1 rusanov -W2 slope 1 Figure 5. Order of convergence for the smooth potential W (x) = x 2 /2 (left) and for the pointy potential W (x) = |x| (right).
In the case of the smooth potential, the schemes behave nicely and we recover the order 0.5 that we had in the one dimensional setting. Concerning the pointy potential, however, some remarks are necessary: (1) The upwind scheme seems not to converge to the correct solution, while the Rusanov scheme does. The solution with the upwind scheme is not blowing up but the velocity at wich the two Dirac masses are getting close to each other is higher than the exact one. Some tests have been run to understand precisely the reason of this non-convergence but no convincing results have been reached. (2) The order at which the Rusanov scheme is converging is 1. The same explanation as the one in Section 3 about the numerical diffusion being counter balanced by the non-smooth gradient might apply here.