Daniel Barlet
email: daniel.barlet@univ-lorraine.fr

Barlet Daniel

Algebraic differential equations associated to some polynomials (seconde version)

Keywords:

published or not. The documents may come

Algebraic differential equations associated to some polynomials

1 Introduction and statement of the theorem.

1.1 Introduction.

Since the works of J. Milnor [Mi.68] and E. Brieskorn [Br.70] it is a classical problem to compute the Gauss-Manin connection of a holomorphic function. The first examples are the Brieskorn-Pham singularities and then the quasi-homogeneous isolated singularities. We shall give here a simple calculation which shows that in the next case, that is to say for a non quasi-homogeneous polynomial in n + 1 variables with n + 2 monomials, it is possible to obtain for any given integral of the type γs µ.dx df , where µ is a monomial and γ s a horizontal family of n-cycles in the fibers of f , an explicit differential (global) equation. This differential equation is defined over Q [λ] where λ = 0 is the parameter of such a polynomial (defined modulo a diagonal change of variables) and has a very special form described in our theorem below. This motivate the study of the kind of differential equation which appears in these examples. We prove two interesting factorization theorems which apply to such an equation in section 3.

Statement of the theorem.

We consider a polynomial f ∈ C[x 0 , . . . , x n ] which is the sum of n + 2 monomials f = n+2 j=1 m j where m j := σ j .x α j , with σ j ∈ C * and α j ∈ N n+1 . Define the matrix with (n + 1) lines and (n + 2) columns M = (α i,j ) and let M be the square matrix obtained from M by adding a first line equal to (1, . . . , 1). We shall assume the following condition:

(C) The rank of M is n + 2.

Remark that this condition is equivalent to the fact that such a f is not quasihomogeneous.

Example. Let p 0 , . . . , p n be integers bigger or equal to 2 and β ∈ (N * ) n+1 such that n 0

β i p i = 1, then f (x 0 , . . . , x n ) := n 0 σ i .x p i i + ρ.x β
satisfies the condition (C), assuming that σ 0 . . . σ n .ρ = 0.

Consider now σ 1 , . . . , σ n+2 ∈ C * and define as before

f (x) := n+2 j=1 σ j .x α j .
A diagonal linear change of variables allows to reduce the study to the case where

f (x) = n+1 j=1
x α j + λ.x α n+2 (@) for some λ ∈ C * . So in what follows, we shall assume that m j = x α j for j ∈ [1, n+1] and m n+2 = λ.x n+2 where λ ∈ C * is a fixed parameter.

In our setting it will be convenient to construct differential equations using, instead of the Weyl algebra C < s, ∂ ∂s >, the algebra A := C < a, b >≃ C < s, ∂ ∂s -1 > defined by the commutation relation a.b -b.a = b 2 which is the translation of the Leibnitz rule, as we shall make a acts on holomorphic functions as multiplication by s and b as ( ∂ ∂s ) -1 := s 0 . The traduction of our result in usual terms will be immediate (see the remark following the theorem below).

For such a polynomial f we define the unitary sub-algebra A f ⊂ C[x 0 , . . . , x n ] generated by the monomials of f . It is the image of R := C[m 1 , . . . , m n+2 ] by the algebra morphism defined by m j → x α j , j ∈ [1, n + 1] and m n+2 → λ.x α n+2 , where λ = 0 is a fixed parameter. We define also

E f := Ω n+1 df ∧ dΩ n-1
where Ω where γ s is the horizontal family of cycles in the fibers of f associated to γ, satisfies the global algebraic differential equation

k := C[x 0 , . . . , x n ] ⊗ C Λ k ((C n+1 ) * ) where (C n+1 ) * := ⊕ n i=0 C .dx i . Then E f is a left A-
P (µ)[F µ,γ ] = 0.
2. Each of the elements P d and P d+h are products in the algebra A of elements of the form a -r.b where r is a rational number.

3. The elements P d and P d+h are computable by an easy linear algebra algorithm.

Remark. To reach the algebraic differential equation in the usual way just look at b -d-h P (µ) which is the sum of a monic polynomial of degree d + h (with rational roots) in s. ∂ ∂s , (see the property 2.) minus c.λ r .( ∂ ∂s ) h composed with a monic polynomial of degree d in s. ∂ ∂s (also with rational roots). The order of this differential equation is d + h and the coefficient of ( ∂ ∂s ) d+h is equal to s d+h -c.λ r .s d , showing that the singular points of this differential equation outside the origine are solutions of s h -c.λ r = 0.

Example. Assume that p ∈ f -1 (0) is a singular point of a polynomial f satisfying the condition (C). Then choose for U a Milnor ball for f at p and for S the corresponding punctured disc with center 0. Then they fulfill the hypothesis in the theorem, and so we obtain, for each given monomial µ, an algebraic global differential equation satisfied by the function F µ,γ for any choice of γ in the n-th homology group of the Milnor fiber of f at the point p.

2 Proof of the theorem.

2.1 Some A-modules.

We consider a polynomial f ∈ C[x 0 , . . . , x n ] with n + 2 monomials m 1 , . . . , m n+2 . So f = n+2 j=1 m j . We write m j := σ j .x α j with σ j ∈ C * , and α j ∈ N n+1 . Using a linear diagonal change of variables and reordering the variables if necessary we shall assume the following conditions :

i) We have σ j = 1 for all j ∈ [1, n + 1] and σ n+2 = λ ∈ C * .
ii) We assume that α 1 , . . . , α n+1 are linearly independent in Q n+1 .

Then we shall write

α n+2 = n+1 j=1 ρ j .α j with ρ j ∈ Q. (@)
We define

H := {j ∈ [1, n + 1] / ρ j = 0}; J + := {j ∈ [1, n + 1] / ρ j > 0}; J -:= {j ∈ [1, n + 1] / ρ j < 0}.
Let |r| be the smallest positive integer such that |r|.ρ j := p j is an integer for each j ∈ [1, n + 1]. Write now the relation above as

|r|.α n+2 + j∈J - (-p j ).α j = j∈J + p j .α j .
Now define d + h and d as respectively the supremum and infimum of the two positive integers |r| + j∈J -(-p j ) and j∈J + p j . Here d and h are positive because the positivity of d is consequence of the fact that |r| ≥ 1 and that at least one p j is positive. The positivity of h is consequence of the fact that the equality of these two integers would imply that the first line in M satisfies the same linear relation (@) than all the other lines in M , contradicting our hypothesis (C). Now the relation above gives the relation between the monomials (m j ) j∈ [1,n+2] :

m |r| n+2 . j∈J - m -p j j = λ |r| . j∈J + m p j j
and we shall write it

m ∆ = λ r .m δ (@@)
where ∆ and δ are in N n+2 of respective weight d + h and d, and with zero component for each h ∈ H. Note that we have r = ±|r| and so that r is in Z * .

We shall also use the following observation later on :

Observation. For j ∈ [1, n + 1] the j-th element of the first column of the matrix M-1 is zero if and only if j is in H.

Consider the free polynomial algebra R := ⊕ ∞ q=0 R q = C[m 1 , . . . , m n+2 ] graded by the total degree in the variables m 1 , . . . , m n+2 and consider the algebra morphism ϕ : R → C[x 0 , . . . , x n ] sending the variable m j in R to the monomial x α j in C[x 0 , . . . , x n ] for j ∈ [1, n + 1] and m n+2 to λ.x α n+2 where λ ∈ C * is a fixed parameter. We shall need an extra variable µ and we associated to it a multi-index β ∈ N n+1 . Then we consider the free, rank 1 R-module R.µ and the R-linear map ϕ µ : R.µ → C[x 0 , . . . , x n ] by sending µ to x β . We construct now on the two vector spaces R.µ and C[x 0 , . . . , x n ] the C -linear operators b i , i ∈ [0, n] and a as follows : i) For γ ∈ N n+1 and i ∈ [0, n] let Γ i (γ, µ) := β i + n+k j=1 α i,j .γ j and define in the vector space R.µ

(Γ i (γ, µ) + 1).b i (m γ .µ) = n+2 j=1 α i,j .m j .m γ .µ and a(m γ .µ) = n+2 j=1 m j .m γ .µ
ii) For ζ ∈ N n+1 and i ∈ [0, n] define on the vector space C[x 0 , . . . , x n ] :

(ζ i +1).b i (x ζ ) = n+1 j=1 α i,j .x α j +ζ + λ.x α n+2 +ζ = x ζ .(x i . ∂f ∂x i ) and a(x ζ ) := f.x ζ
Note that the operator b i on C[x 0 , . . . , x n ] may also be written as :

P → ∂f ∂x i . x i 0 P.dt = x i . ∂f ∂x i . 1 x i . x i 0 P.dt = n+2 j=1
α i,j .m j . 1 x i .

x i 0 P.dt .

Remark that C[x 0 , . . . , x n ].x β ⊂ C[x 0 , . . . ,
x n ] is stable by the operators b i and a for any choice of β.

The commutation relations a.b i -b i .a = b 2 i are easy to verify on R.µ and C[x 0 , . . . , x n ]. It is also easy to see that ϕ µ commutes with these operators. We shall denote V p (µ) the vector subspace of R p .µ defined inductively as follows:

1. V 0 (µ) = {0}. 2. V p+1 (µ) := i =0 (b i -b 0 ).R p .µ + b 0 .V p (µ).
Remark that, by definition, we have

b i .V p (µ) ⊂ (b i -b 0 ).V p (µ) + b 0 .V p (µ) ⊂ V p+1 (µ)
for all i ∈ [0, n] and for all p ≥ 0. We have also a.

V p (µ) ⊂ V p+1 (µ) because of the identity a.(b i -b 0 ) = (b i -b 0 ).a + (b i -b 0 ).b i + b 0 .(b i -b 0 ).
Then we define the graded vector space S(µ) := ⊕ p≥0 S p (µ) where we put

S p (µ) := R p .µ V p (µ).
The quotient map π µ : R.µ → S(µ) is graded.

Note that the algebra A has a natural grading A = ⊕ ∞ q=0 A q where A q is the vector subspace of homogeneous elements in a, b of degree q.

Remark. An element P ∈ A q which is monic in a can be written

P = (a -r 1 .b) . . . (a -r q .b)
for some suitable choice of complex numbers r 1 , . . . , r q ; see [B.09] or write b -q .P as a polynomial in b -1 .a.

On S(µ) we have a structure of left graded A-module deduced from the action of the operators a and b := b 0 = • • • = b n which are equal on S(µ). Note that if we consider R.µ as a left A-module with the actions of a and b 0 , the quotient map π µ : R.µ → S(µ) becomes a surjective map of graded left A-modules.

Proposition 2.1.1 Fix β ∈ N n+1 and q ∈ N. For each γ ∈ N n+2 such that |γ| = q, there exists an element χ q (µ)(m γ ) ∈ A q , which is a product of homogeneous elements in (a,b) of degree 1 of the form η.a + θ.b with η and θ in Q, such that we have in S(µ) the equality :

π µ (m γ .µ) = χ q (µ)(m γ ).µ.
Moreover, for every γ such that γ h = 0 for each h ∈ H, the corresponding element χ q (µ)(m γ ) ∈ A q is monic in a up to a non zero rational number.

Proof. Consider the following equalities in S(µ) :

a[π µ (m γ .µ)] = f.m γ .µ.dx = n+2 j=1 π µ (m j .m γ .µ) (Γ i (γ, µ) + 1).b[π µ (m γ .µ)] = n+2 j=1 α i,j .π µ (m j .m γ .µ) ∀i ∈ [0, n].
(@@@) as a linear system with matrix M and unknown the column vector

π µ (m j .m γ .µ) , j ∈ [1, n + 2] in S(µ) n+2 .
The invertibility of M (on the field Q) gives that for each j ∈ [1, n+2] we may write π µ (m j .m γ .µ) = (η j .a + θ j .b).π µ (m γ .µ) for some rational numbers η j and θ j . Assuming that, for each γ with |γ| = q the element χ q (µ)(m γ ) is a product of q homogeneous elements in (a,b) of degree 1 of the form η.a + θ.b with η and θ in Q, we conclude that for each γ with |γ| = q + 1 the element χ q+1 (µ)(m γ ) also satisfies this condition. Moreover, as the first column of M-1 has a non zero j-th component for j ∈ H (see the observation above), we have η j = 0 for each j ∈ H. So, assuming inductively that, for any γ such that |γ| = q, with zero j-th component for j ∈ H, we have already constructed a suitable χ q (µ)(m γ ) ∈ A q this implies that for any γ ∈ N n+2 such that |γ| = q + 1 and with no component in H, we may write m γ = m j .m γ for some j ∈ H and some γ ∈ N n+2 such that |γ| = q has zero components in H, and define χ q+1 (µ)(m γ ) := (η j .a + θ j .b).χ q (µ)(m γ ) where η j is in Q * . So we conclude that for such a γ we may choose χ q+1 (µ)(m γ ) as a monic element in a belonging to A q+1 , up to a non zero rational number. Now using only the actions of a and b 0 gives a structure of left A-module on R.µ and a graded A-linear map .µ : A → R.µ sending 1 to µ. Composed with the quotient map π µ we obtain a graded A-linear map π 0 : A → S µ . The following corollary is an obvious consequence of the proposition above.

Corollary 2.1.2 The graded A-linear map π 0 : A → S µ is surjective.

End of the proof of the theorem 1.2.1

The following lemma is now easy.

Lemma 2.2.1 The map ϕ µ induces a commutative diagram

A π 0 B B B B B B B B B .µ / / R.µ π ϕµ / / C[x 0 , . . . , x n ] π ′ S µ ψµ / / E f where the map ψ µ is A-linear. The image of ψ µ is equal to A.[x β .dx] ⊂ E f .
proof. The A-linearity is a direct consequence of the compatibility of the map ϕ µ with the operators a and b i , i ∈ [0, n] and the fact that these operators induce a and b on S(µ) and on E f . The identification of the image of ψ µ is then a consequence of the corollary of the proposition 2.1.1 and of the commutativity of the diagram above.

As the relation (@@) holds in C[x 0 , . . . , x n ], we have also m ∆ .x β .dx = λ r .m δ .x β .dx in E f . As ∆ and δ have no components in H, there exists elements P d+h and P d respectively in A d+h and A d which are equal, modulo some invertible rational numbers, to

χ d+h (µ)[m ∆ ] and χ d (µ)[m δ ].
Then we obtain in E f the equality

P d+h -cλ r .P d [x β .dx] = 0
for some non zero rational number c, thanks to the A-linearity of the map ψ µ . Now assume the existence of S and U as in the property 1. of the statement of the theorem. Note that we may assume S simply connected. Then define O b (S) as the vector space of holomorphic functions on S such that |s| 1+ε .|f (s)| → 0 when s → 0, for some ε > 0. We want to show that the C -linear map

I : E f -→ O b (S) ⊗ C H n (F, C) * given by µ → γ → γs µ.dx df is A-linear. Here the actions on O b (S)⊗ C H n (F, C) * are given by the rules a(g ⊗ v) := s.g ⊗ v and b(g ⊗ v) := ( s 0 g) ⊗ v.
The commutation with a of the map I is clear. The commutation of the map I with b is consequence of the derivation formula for the integral of a holomorphic n-form :

∂ ∂s γs ξ = γs dξ df .
Note that for a holomorphic n-form the function γs ξ has only bounded terms in its asymptotic expansion when s → 0 (see [M.75]), so its derivative is in the space O b (S) defined above because such a function admits a "standard" convergent asymptotic (multivalued) expansion when s → 0. proof. Let F be the vector with components f, x 0 . ∂f ∂x 0 , . . . , x n . ∂f ∂xn and W the vector with components m 1 , . . . , m n+2 ; then we have F = M.W and this allows to compute each m j as a linear combination of the components of F as M is invertible. Moreover, for j ∈ H, the coefficient of f in this linear combination does not vanishes. Replacing in the equation (@@) :

m ∆ = c.λ r .m δ
gives a monic degree d + h = |∆| polynomial in f with coefficients in the ring C[x 0 . ∂f ∂x 0 , . . . , x n . ∂f ∂xn ] which vanishes identically.

For instance, let f := x 2 +y 3 +z 4 +λ.x.y.z. Then let

X := f -1 2 x.f ′ x -1 3 y.f ′ y -1 4 z.f ′ z . The computation gives 2 6 .3 4 .4 3 .X 12 = λ 12 .(X + 13x.f ′ x ) 6 (X + 13y.f ′ y ) 4 (X + 13z.f ′ z ) 3
which gives an integral dependance relation of degree 13 of f on the ring

Q[x.f ′ x , y.f ′ y , z.f ′ z ], for λ fixed in Q * .
Remark. In fact we have shown that, in general, for λ given in Q * , f is integrally dependant on the ring Q[x 0 . ∂f ∂x 0 , . . . , x n . ∂f ∂xn ].

Examples.

First let me give an example with 3 variables, small (total) degree = 6 such that d = 61 and h = 15. Let f = x 3 .y + y 4 .z + z 5 .x + λ.x 2 .y 2 .z 2 . Let α j , j ∈ [1, 4], the exponents of the monomials of f . Note that the singularity of f at the origin is not isolated. 20) 20 . This is not simple to see by a direct computation.

We shall give now a family of examples with n = 3 (so 4 variables) with some symetry. Let x 0 , x 1 , x 2 , x 3 be the coordinates in C 4 , and note σ the circular permutation of the coordinates defined by σ(x i ) = x i+1 , i ∈ Z 4.Z. Now consider a multi-index α ∈ N 4 with the following conditions : i) |α| ≥ 5.

ii) The matrix M 0 given by M 0 (i, j) := α i+j with (i, j) ∈ Z 4.Z has rank 4 1 .

Then we shall consider the polynomial

f (x) := 3 j=0 σ j (x α ) + λ.x 1
where λ is a non zero complex parameter and where we use the following conventions : for a monomial µ in the variables x 0 , x 1 , x 2 , x 3 we define σ(µ) to be the monomial obtained from µ by applying the circular permutation σ to the variables, and we note 1 the multi-index (1, 1, 1, 1).

The relation between the five monomials in f is simply

λ |α| . 3 j=0 σ j (x α ) = (λ.x 1 ) |α| .
So we will compute the images of the monomials of both sides of this relation by the map ψ 1 in E f as P d+h [dx] and P d [dx] where P m are in A m and monic in a up to a non zero constant. Here we have d = 4 and h = |α| -4.

Computation of

3 j=0 σ j x α . It will be convenient to use the notation B(x γ ) for the vector with components (γ i + 1).b i (x γ ), i ∈ [0, 3] ≃ Z 4.Z, and Σ(x γ ) for the vector with components σ j (x γ ), j ∈ [0, 3] ≃ Z 4.Z, when x γ is a monomial. Then we have, for each monomial m ∈ C[x 0 , . . . , x n ] the equation

B(m) = m. M 0 .Σ(x α ) + λ.x 1 .v 0 1 Note that this is true if and only if α 0 + α 2 = α 1 + α 3 .
where v 0 is the column vector v 0 := (1, 1, 1, 1). We shall use successively the values 1,

x α , x α .σ(x α ), x α .σ(x α ).σ 2 (x α ) for m in this formula. But remark that, because we have b i (m) = b(m) for each i ∈ [0, 3], we may write in S(1) = C[m 1 , . . . , m 5 ] V (1):

B(1) = b(1).v 0 , B(x α ) = b(x α ).(v 0 + M 0 .v 1 ) B(x α .σ(x α )) = b(x α .σ(x α )).(v 0 + M 0 .v 2 ) + λ.x 1 .x α .σ(x α ).v 0 B(x α .σ(x α ).σ 2 (x α )) = b(x α .σ(x α ).σ 2 (x α )).(v 0 + M 0 .v 3 ) + λ.x 1 .x α .σ(x α ).σ 2 (x α ).v 0
with the following values for the (column) vectors

v 0 := t (1, 1, 1, 1), v 1 := t (1, 0, 0, 0), v 2 = t (1, 1, 0, 0) and v 3 = t (1, 1, 1, 0).
Now remark that M 0 .v 0 = |α|.v 0 and recall that M 0 is invertible. This allows to deduce the following equalities in E f :

|α|.Σ(x α )) = (b(1) -λ.x 1 ).v 0 |α|.x α .Σ(x α ) = (b(x α ) -λ.x 1 .x α ).v 0 + |α|.b(x α ).v 1 |α|.x α .σ(x α ).Σ(x α ) = b(x α .σ(x α )) -λ.x 1 .x α .σ(x α ) .v 0 + |α|.b(x α .σ(x α )).v 2 |α|.x α .σ(x α ).σ 2 (x α ).Σ(x α ) = b(x α .σ(x α .σ 2 (x α )) -λ.x 1 .x α .σ(x α .σ 2 (x α ) .v 0 + + |α|.b(x α .σ(x α )).σ 2 (x α ).v 3
The first equality implies the equality of σ j (x α ) for all j ∈ [0, 3] and then we have a(1) = 4.x α + λ.x 1 and b(1) = |α|.x α + λ.x 1 and then (a -b)(1) = (4 -|α|).x α . Using the other equalities gives after some easy computations

(a -2b)(x α ) = (4 -|α|).x α .σ(x α ) (a -3b) x α .σ(x α ) = (4 -|α|).x α .σ(x α ).σ 2 (x α ) (a -4b) x α .σ(x α ).σ 2 (x α )) = (4 -|α|).x α .σ(x α ).σ 2 (x α ).σ 3 (x α ).
So we have

(a -4b)(a -3b)(a -2b)(a -b)[1] = (|α| -4) 4 .x α .σ(x α ).σ 2 (x α ).σ 3 (x α ). (@)
Now we compute the action of b on the monomial m p := (x 0 .x 1 .x 2 .x 3 ) p . We have

B(m p ) = (p + 1).b(m p ).v 0 = M 0 .m p .Σ(x α ) + λ.m p+1 .v 0
and so m p .σ j (x α ) is independent on j ∈ [0, 3] and equal to So we conclude that the element 2 (a -4b).

|α|-2 p=0 a -

4(p + 1) |α| .b -λ |α| .(|α| -4) |α|-4 .(a -3b)(a -2b)(a -b)
annihilated the class of dx in E f .

Note that this element depends on α only by the number |α|. So the differential equation for the periods of [dx] for such an f depends only on the total degree of the monomial x α and not of the precise choice of such a monomial x α . Remark also that for |α| = 5 the element of A under the brackets is invariant by the anti-automorphism θ 4 of A defined by the relations

θ 4 (x.y) = θ 4 (y).θ 4 (x), θ 4 (1) = 1, θ 4 (a) = a -4b, θ 4 (b) = -b.
So it is remarkable that we find the same differential equation with this symetry property for the period of dx ∧ dy ∧ dz ∧ dt for the following polynomials :

x 5 + y 5 + z 5 + t 5 + λ.x.y.z. 3 The factorization theorems.

3.1 The decomposition theorem.

Here we shall work with the completion à of the algebra A relative to its b-adic filtration given by the two-sided ideals b ν .A for ν ∈ N.

Recall that a (a,b)-module is a left Ã-module which is free and of finite rank on the commutative sub-algebra C We shall say that a (a,b)-module is monogenic when it may be generated by one element as a left Ã-module.

Recall also that a (a,b)-module E is local, respectively simple pole, if there exists an integer N ≥ 1, respectively N = 1, such that a N .E ⊂ b.E.

The purpose of this section is to prove a decomposition theorem for any (a,b)module E. In the case where E := à Ã.P where P ∈ A is monic in a, as in the differential equations constructed in the previous section, it will give a factorisation of P in à P = P 1 . . . . P k

where each P i is monic in a and with a class modulo b. à of the form (a -v i ) d i where v i is a complex number and d i a positive integer. So, in the decomposition theorem which corresponds to the "decomposition" of the differential equation in its different singular points, we find in this case local factors of the same type (i.e of the form à Ã.P i ).

G v of E such that G v b.G v ≃ F v , and such that (a -v) d(v) .G v ⊂ b.G v , where d(v)
is the multiplicity of the root v in the minimal polynomial of the action of a on F . Moreover, G v is of rank dim C F v and we have the following decompositon as a direct sum of (a,b)-modules

E = ⊕ v G v .
Proof. First we shall construct G v for a given eigenvalue v of a acting on F . We may assume without loss of generality that v = 0, and we note d(0) = d for short. Define

G := {x ∈ E / a n.d .x ∈ b n .E ∀n ∈ N}.
This vector subspace is clearly stable by a. It is also stable by b using the following identities which are easy to prove by induction on ν ∈ N

a ν .b = b.(a + b) ν (a + b) ν = a ν + ν.a ν-1 .b a ν .b = b.a ν + ν.b.a ν-1 .b If a (n+1).d .x ∈ b n+1 .E and a n.d .b.x ∈ b n+1 .E then a (n+1).d .b.x = b.a (n+1).d .x + (n + 1).d.b.a d-1 .(a n.d .b.x) ∈ b n+2 .E.
So we obtain that b.G ⊂ H where H ⊂ G is defined by

H := {y ∈ E / a n.d .y ∈ b n+1 .E}.
We shall show now that for any x ∈ E such that x lies in F 0 + b.E, there exists an element ξ ∈ G such that ξ -x ∈ b.E. Fix a vector space decomposition E := F ⊕ b.E and denote Fv the pull back of F v by the bijection F → F . Define then

L := ⊕ v =0 Fv .
Assume we have already found y 1 , . . . , y n-1 in L such that

a n.d .(x - n-1 j=1 b j .y j ) = b n .x n (@ n )
for some x n ∈ E. For n = 1 this is clear from our assumption that x ∈ F 0 + b.E and a d .F 0 ⊂ b.E. Assume that (@ n ) is true for n ≥ 1. We shall prove (@ n+1 ). Write

x n = t n + z n + b.w n where t n is in F0 , z n in L and w n ∈ E. As a is bijective on L = L + b.E b.E it is also bijective on b n .L = b n .L + b n+1 .E b n+1 .E, because the identity a.b n = b n .a + nb n+1
gives the commutativity of the diagram

E b.E a / / b n E b.E b n b n .E b n+1 .E a / / b n .E b n+1 .E So we may find y n ∈ Fv such that a n.d .b n .y n = b n .(z n + b.ζ n ) where ζ n is in E. Then we obtain a (n+1).d .b n .y n = a d .b n .(z n + b.ζ n )
and so Remark. As each G v is a direct factor (as (a,b)-module) of E, each of these (a,b)-modules is monogenic when E is monogenic.

a (n+1).d .(x - n j=1 b j .y j ) = a d .b n .x n -a (n+1).d .b n .y n = a d .b n .t n + a d .b n .z n + a d .b n+1 .w n -a d .b n .(z n + b.ζ n ). But a d .b n .t n = b n .a d .t n + b n+1 .E and it belongs to b n+1 .E as t n ∈ F0 implies that a d .t n is in b.E. So we conclude that a (n+1).d .(x -n j=1 b j .y j ) is in b n+1 .E. This proves the existence of y := +∞ j=1 b j .y j ∈ b. C[[b]].L such that for any n ∈ N we have a n.d .(x -y) = a n.d .(x - n-1 j=0 b j .y j ) + a n.d .b n . ∞ j=1 b j .y n+j ∈ b n .E.
The following corollary is now immediate. Of course, if we decompose the generator e of E according to the direct decomposition

E = ⊕ v i ,i∈[1,l] G v i : e = l i=1
e v i then P l generates the annihilator of e v l in E (or in G v l ). So we have an isomorphism à Ã.P l ≃ G v l by sending 1 to e v l . Again the polynomial P l-1 generates the annihilator of P l .e v l-1 . As P l .e v l-1 is a generator of G v l-1 because in a local (a,b)-module the topology defined by the a-filtration is complete, so the action of an element in à with a non zero constant term is invertible, we have an isomorphism à Ã.P l-1 ≃ G v l-1 obtained by sending 1 to P l .e v l-1 etc... Remark that changing the order of the v i changes the polynomials P i in the decomposition above but only change the order of the set of the isomorphism classes of the à Ã.P i .

The following consequence of the theorem 1.2.1 and the results above is obvious.

Corollary 3.1.4 In the situation of the theorem 1.2.1 the element P d is a left multiple of the Bernstein element of the left Ã-module generated by

[µ.dx] in E f ⊗ C[b] ⊗ C[[b]] b -torsion.
This is, of course, a non trivial information on the monodromies at the origin for the periods integrals considered in the beginning of this paper.

Irregularity and the second factorization theorem.

In this paragraph we examine what happens in the previous factorization when the initial part of P is not monic in a. Proof. First remark that if I and J are normal sub-modules such that E I and E J are regular, then E I ∩ J is again regular because E I ∩ J is a sub-module of the direct sum of the regular modules E I and E J. Then let I be the intersection of all normal sub-modules J of E such that E J is regular. As any descending chain of normal sub-modules is finite, it is clear that I is the smallest normal sub-module of E such that E I is regular. Let f : I → F be a Ã-linear map, where F has a simple pole, and let K be the kernel of f . Then in the exact sequence of (a,b)-modules 0 → I K → E K → E I → 0 the (a,b)-modules I K5 and E I are regular. So is E K. Then E K is a sub-module of a simple pole (a,b)-module, and then I ⊂ K. This implies f = 0, and so I is totally irregular. Lemma 3.2.3 Let P ∈ Ã be a monic polynomial in a of degree d. Assume that the initial form of P in (a,b) is equal to b q , for some integer q < d. Then the (a,b)-module E := Ã Ã.P is totally irregular.

Proof. Let f : E → F be a Ã-linear map such that F is a simple pole (a,b)-module. Let x := f (1) ∈ F . We have P.x = 0 in F , and so b q .x is in b q+1 .F , because of our assumption on the initial form of P and on the simple pole for F . So x ∈ b.F . Then the image of f lies in b.F which is again a simple pole (a,b)-module. Iteration of this gives that the image of f is in b n .F for any integer n. But ∩ n∈N b n .F = 0, and so f = 0. Theorem 3.2.4 Let P be a element in à with the following properties :

(i) P is a monic polynomial in a of degree d + h with h ≥ 1.

(ii) The initial form in (a,b) of P has degree d ≥ q ≥ 0, and is of the form ρ.b q .P d-q with ρ ∈ C * where P d-q is monic in a.

(iii) Modulo b. Ã, we have P = a d+h .

Then there exists elements Z and Q in à which are polynomials in a of respective degrees q + h and at most d -q -1, with respective valuations in (a,b) at least q + 1 and d -q + 1, such that we have in à the equality (ρ.b q + Z).(P d-q + Q) = P.

Defining the (a,b) modules E := Ã Ã.P, F := Ã Ã.(P d-q + Q) and F := Ã Ã.(ρ.b q + Z),

we have the exact sequence of (a,b)-modules of respective ranks q + h, d + h, d -q :

0 → I → E → F → 0
and I is the totally irregular part of E.

Remark. Assuming that the polynomial B associated to P d-q by the relation (-b) d-q .B(-b -1 .a) = P d-q in A, has negative rational roots, the (a,b)-module F is a geometric (a,b)-module generated by one element as left Ã-module ( see [B.09]) and B is the Bernstein polynomial of F .

Proof. We shall construct by induction on n ≥ d + 1 homogeneous elements in (a,b) ξ n-d+q and η n-q of respective degrees n -d + q and n -q, which are polynomial in a of respective degrees less or equal to q + h and d -q -1, such that for each N ≥ d + 1 we have

(b q + N n=d+1
ξ n-d+q ).(P d-q + N n=d+1 η n-q ) = P + Y N +1 (@ N )

where Y N +1 has valuation in (a,b) ≥ N + 1 and is a polynomial in a of degree ≤ d + h.

For N = d + 1 denote (P ) d+1 the homogeneous part in (a,b) of degree d + 1 of P , and make the right division by P d-q in A :

(P ) d+1 = ξ q+1 .P d-q + R d+1

where R d+1 has degree ≤ d -q -1 in a and is homogeneous in (a,b) of degree d + 1. So we may write R d+1 = ρ.b q .η d-q+1 . Let Y 0 d+2 := P -ρ.b q .P d-q -(P ) d+1 . Then we have (ρ.b q + ξ q+1 ).(P d-q + η d-q+1 ) = ρ.b q .P d-q + ρ.b q .η d-q+1 + ξ q+1 .P d-q + ξ q+1 .η d-q+1

= ρ.b q .P d-q + (P ) d+1 + ξ q+1 .η where ξ N +1-d+q is homogeneous in (a,b) of degree N + 1 -d + q, where R N +1 is homogeneous in (a,b) of degree N + 1, and they have respectively degrees less than q + h and d -q -1 in a. Then we may write R N +1 = ρ.b q .η N -q+1 . Now we have where we defined

Z N +2 := ξ n-d+q ). N +1 d+1 η n-q + ( N d+1
ξ n-d+q ).η N +1-q and Y N +2 := Y 0 N +2 + Z N +2 . This complete the proof of our induction step. Now the series ∞ n=d+1 ξ n-d+q and ∞ n=d+1 η n-q converge in à to Z and Q which are polynomials in a of degree6 respectively equal to q + h and ≤ d -q -1. So we obtain the desired decomposition of P by defining Z := ∞ n=d+1 ξ n-d+q and Q := ∞ n=d+1 η n-q . Now the lemma 3.2.3 implies that I is totally irregular and the lemma 3.2.2 implies that E I is regular. To complete the proof is then easy.

F

  module where a acts by multiplication by f and where b acts as df ∧ d -1 .We shall fix in the sequel a monomial µ := x β in C[x 0 , . . . , x n ] and we shall denote ϕ µ : R → E f the map given by m j → [m j .µ.dx]. Theorem 1.2.1 In the situation above there exist positive integers d and h such that for any choice of the monomial µ ∈ C[x 0 , . . . , x n ] we may construct in a natural way (explained below) homogeneous elements in (a, b) of degree d and d + h in the algebra A, noted P d (µ) and P d+h (µ), which are monic in a, and such that the element P (µ) := P d+h (µ) -c.λ r .P d (µ) ∈ A with some r ∈ Z * and c ∈ Q * , has the following properties : 1. Assume that on a connected open set S in C * with the origin on its boundary, we have an open set U ⊂ f -1 (S) such that f : U → S is a C ∞ -fibration with fiber F ; then for any compact n-cycle γ ∈ H n (F, C), the holomorphic function on S defined by (see [M.74])

(

  |α|) -1 . (p + 1).b(m p ) -λ.m p+1 . m p ) = λ.(|α| -4).m p+1 .

  [[b]] of Ã. For more on (a,b)-modules (generalized Brieskorn modules) see [B.06].

For a proof

  of the following proposition see [B.09]. Proposition 3.1.1 Let I be a left ideal in à such that E := à I is a rank d (a,b)-module and such that x d is the minimal polynomial of the action of a on E b.E. Then there exists a monic degree d polynomial in a with coefficients in C[[b]] such that I = Ã.P . Moreover, such a P is unique, and E is regular if and only if the homogeneous initial part of P in (a,b) is of degree d. When E is regular, the Bernstein element 3 of E is the initial part of P . Conversely, if E is a monogenic local (a,b)-module, for any generator e of E, its annihilator in à is of the type I = Ã.P where P is a monic polynomial in a of degree d with coefficients in C[[b]]. Theorem 3.1.2 (Decomposition) Let E be an (a,b)-module. For each spectral subspace F v for the eigenvalue v of the action of a on F := E b.E, there exists a normal 4 submodule

  So we have found an element z := x -y ∈ G such that x -z ∈ b.E. Note that this implies that the map G b.G → E b.E induced by the inclusion G ⊂ E has exactely F 0 as image, because a d .x ∈ b.E implies that x ∈ F 0 + b.E. Let us show now that H = b.G. As the inclusion b.G ⊂ H is already proved, take x ∈ H. We have x = b.z for some z ∈ E (using n = 0 in the definition of H), and now a n.d .b.z = b.a n.d .z + n.d.b.a n.d-1 .b.z = b.a n.d .z + n.d.b.a d-1 .a (n-1).d .b.z and as a(n-1).d .b.z ∈ b n .E and a n.d .b.z ∈ b n+1 .E ( x ∈ H) we get b.a n.d .z ∈ b n+1 .E and so a n.d .z ∈ b n .E. So z is in G. This implies that a d .G ⊂ H = b.G. We shall prove now that G is normal in E. We shall first show that if x lies in b k .E ∩ G for some k ≥ 1 there exists y ∈ G such that x -b k .y ∈ b k+1 .E. Assume that x ∈ G satisfies x = b k .z for some z ∈ E and k ≥ 1. Write z = u+v +b.ξ with u ∈ F0 , v ∈ L and ξ ∈ E. Then x = b k .z = b k .u+b k .v +b k+1.ξ is the spectral decomposition relative to the action of a on x in b k .E b k+1 .E. This implies that v = 0 and so x = b k .u + b k+1 .ξ and the previous step of our proof gives an element y ∈ G such that y -u ∈ b.E. So we have x -b k .y ∈ b k+1 .E proving our assertion. So, for such an x ∈ b k .E ∩ G we can construct inductively y k , . . . , y n , . . . in G such that x = b k . +∞ j=0 b j .y k+j . So x is in b k .G, and G is normal. As a consequence, the map G b.G → E b.E is injective and is an isomorphism of G b.G on F 0 . So G has rank dim C (F 0 ), and satisfies a d .G ⊂ b.G.To prove that the sum of the normal sub-modules G v is direct is an easy exercice left to the reader. Then to prove the direct sum decomposition of E it enough to prove that E b.E is the direct sum of the images of the G v b.G v . But this is the spectral decomposition of E b.E for the action induced by a from what we proved above.

  Corollary 3.1.3 Let P ∈ Ã be a monic polynomial in a with coefficients in C[[b]], and define E := Ã Ã.P . Let v 1 , . . . , v l be the distinct eigenvalues of a acting on E b.E, and let d 1 , . . . , d l be the dimension of the corresponding spectral subspaces. Then there exists monic polynomials P 1 , . . . , P l respectiveley in a -v 1 , . . . , a -v l , of respective degrees d 1 , . . . , d l , with coefficients in C[[b]], such that we have in A the equality P = P 1 . . . P l .

  Definition 3.2.1 A (a,b)-module E will be called totally irregular (at 0) if any Ã-linear map f : E → F in a simple pole (a,b)-module F is the zero map. This terminology is compatible with the usual notion of irregularity of a differential equation at the point 0. It concerns in fact only G 0 the local part at 0 in the previous decomposition theorem 3.1.2 for an (a,b)-module. A quotient of a totally irregular (a,b)-module is again totally irregular. Lemma 3.2.2 Let E be a local (a,b)-module. Then there exists a smallest normal sub-module I ⊂ E such that the quotient E I is regular. This sub-module I is totally irregular.

  d-q+1 = P + Y 0 d+2 + ξ q+1 .η d-q+1 = P + Y d+2where Y d+2 has a valuation in (a,b) ≥ d + 2 and degree ≤ d + h in a.Assume now that we have proved (@ N ) for someN ≥ d + 1. Let Y 0 N +1 the homogeneous part of Y N +1 of degree N +1 in (a,b), and write Y 0 N +2 := Y N +1 -Y 0 N +1 . The right division by P d-q in à gives -Y 0 N +1 = ξ N +1-d+q .P d-q + R N +1

  Y N +1 + ξ N +1-d+q .P d-q + ρ.b q .η N -q+1 + Z N +2 = P + Y 0 N +2 + Z N +2 = P + Y N +2

  In the situation of the theorem 1.2.1 the element f is integrally dependent on the subring C[x 0 . ∂f ∂x 0 , . . . , x n . ∂f ∂xn ] ⊂ C[x 0 , . . . , x n ].

2.3 Integral dependance of f . Lemma 2.3.1

  The equality 61.α 4 = 34.α 1 + 22.α 2 + 20.α 3 gives the relation (λ.x 2 .y 2 .z 2 ) 61 = λ 61 .(x 3 .y) 34 .(y 4 .z) 22 .(z 5 .x) 20 . (R) Of course to compute explicitely P : P 76 -c.λ -61 .P 61 is a little tedious and without real interest. But computing the monomials in the relation (R) modulo b.A gives easily the value of the rational number c. And so the non zero critical values of f are solutions of the equation s 15 = c.λ -61 with c = -(61) 61 .(15) 15 (34) 34 .(22) 22 .(

The Bernstein element Q of E is the homogeneous element in à monic in a of degree d, such that the Bernstein polynomial of E is given by (-b) d .B(-b -1 .a) = Q.

A sub-module G ⊂ E is normal when it satisfies G ∩ b.E = b.G.

As f induces an injection of I K in F which has a simple pole, it is regular.

Remark that ξ q+h has degree exactly q + h and that ξ m has degree ≤ q + h -1 for m = q + h.

Example. Let P d+h and P d-q be homogeneous elements in A with respective degree d + h and d. Assume that P d+h and P d-q are monic in a, then P := P d+h + ρ.b q .P d-q satisfies the hypothesis of the theorem.

In the case of the previous example, let P d+h = ξ h+q .P d-q + ρ.b h+q+1 .η d-q-1 be the right division of P d+h by P d-q in A. Then we have Z = ρ.b q + ξ h+q + Y with Y a degree ≤ h + q -1 polynomial in a with valuation in (a,b) at least equal to 2h + q, and Q = P d-q + b h+1 .η d-q-1 + T where ηd-q-1 and T are degree ≤ d -q -1 polynomials in a with valuation in (a,b) respectively equal to d -q -1 and at least equal to d + 2h -q.

Consequence. If an element x of an (a,b)-module E is killed by an element in à which is monic in a of degree d + h, h ≥ 1, equal to a d+h modulo b. Ã, with initial form equal to ρ.b q .P d-q , where ρ is in C * , where q ≥ 1 and where P d-q is homogeneous of degree d -q in (a,b) and monic in a, then the left Ã-module generated by x in E is a regular (a,b)-module with rank at most d -q and its Bernstein element is a right divisor of P d-q in A.

References.

•


Nicolas Royackkers

Michel Desvignes

Houssam Fawal

Marinette Revenu

Detection and Statistical Analysis of Human Cortical Sulci

HAL is


G Balarac

G.-H Cottet

J.-M Etancelin

J.-B Lagaert

F Perignon

C Picard

Multi-scale problems, high performance computing and hybrid numerical methods

The turbulent transport of a passive scalar is an important and challenging problem in many applications in fluid mechanics. It involves different range of scales in the fluid and in the scalar and requires important computational resources. In this work we show how hybrid numerical methods, combining Eulerian and Lagrangian schemes, are natural tools to address this mutli-scale problem. One in particular shows that in homogeneous turbulence experiments at various Schmidt numbers these methods allow to recover the theoretical predictions of universal scaling at a minimal cost. We also outline how hybrid methods can take advantage of heterogeneous platforms combining CPU and GPU processors.

Introduction

Numerical simulations have become a routine tool to develop, prototype and/or validate products and processes in industry. Applications encompass virtually all sectors of activity from Aeronautics, Automotive industry and Oil exploration to Circuit design, Biomechanics and Animations studios, to name a few. With the need to perform more and more realistic simulations and the advent of supercomputers, available in national or regional centers, the field of High Performance Computing (HPC) is not anymore restricted to academia and scientific grand challenges but starts to reach SMEs.

HPC requires easy and flexible access to HPC facilities, obviously, and to master the appropriate programming language, but also to question the numerical meth- ods and algorithms that are used in the simulations. These methods and algorithms should be adapted both to the physics of the problems to solve and to the architecture of the simulation platforms. Moreover since these platforms are often of an hybrid nature, that is combine different type of processors, typically CPU and GPU processors, one may also wish to develop or use methods which couple different types of algorithms that can be optimally distributed to different types of processors. This is particularly desirable if the problem to solve is multi-level by nature. In that case the different scales that are to be represented can also be resolved on different types of processors. This is hybrid computing. In some sense the nature of the problem and of the hardware inspires the type of mathematical and numerical models that should be used for optimal efficiency.

The purpose of this paper is to describe ongoing work in our group towards hybrid computing for applications in turbulent transport of passive scalar. In the next section we briefly describe the physical context of this work. In section 3 we describe a hybrid method coupling a semi-Lagrangian method for the scalar transport and a spectral method for incompressible flows and we show some results obtained with this method. Section 4 is devoted to the implementation of scalar transport on GPU processors.

Universal scaling in turbulent transport

The prediction of the dynamics of a scalar advected by a turbulent flow is an important challenge in many applications. Some of these applications are illustrated in Figures 1 to 3. Figure 1 shows the dynamics of a pollutant ejected by a sewer in the Los Angeles bay at two different times. Figure 2 shows the atomization of a jet. In that case the transported quantity is the interface water-air interface [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]. Figure 3 shows a similar experiment but in the context of combustion. In this case the transported quantities are concentrations of chemical species [START_REF] Lodato | Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows[END_REF]. All these illustrations share a common feature, namely that very small scales spontaneously appear and need to be captured if accurate predictions are needed for the location of the pollutant, the size of the droplets or the combustion efficiency, respectively, are sought.

The production of small scales in an advected scalar indeed reflects some fundamental turbulence properties and is driven by the value of the Schmidt number, Sc, the ratio between the viscosity of the fluid and the diffusivity of the scalar. if Sc > 1, the so-called Batchelor scale η B which measures the size of the smallest scalar fluctuations is smaller than the smallest length scales of the turbulent flow (the Kolmogorov scale η K ). These scales are related by η B = η K / √ Sc. More precisely, for Sc > 1, Batchelor [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF] reports that the classical Corrsin-Obukhov cascade associated with a k -5/3 law (where k is the wave number) for the scalar variance spectrum [START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF][START_REF] Obukhov | The structure of the temperature field in a turbulent flow[END_REF] is followed by a viscous-convective range with a k -1 power law. This viscous-convective range is followed by the dissipation range, where various theoretical scalings have been proposed for the spectrum [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF][START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF]. A direct conse-quence of this fact is that, for Sc > 1, in numerical simulations the scalar is more demanding, in terms of grid resolution, than the flow itself. It is therefore natural to envision numerical approaches which use different grid resolutions for the scalar and the momentum.

Hybrid particle-spectral method

We consider in the following scalar equation

∂ θ ∂t

+ u • ∇θ = ∇ • (κ∇θ ) (1)
coupled with the incompressible Navier-Stokes equation

∂ u ∂t + u • ∇u = ∇ • (ν∇u) -∇p , ∇ • u = 0, (2)
in a periodic box. κ is the molecular scalar diffusivity, ν the flow viscosity and u the flow velocity field. Using different grid resolutions for the scalar and the flow has already been considered for instance in [START_REF] Cottet | Subgrid particle resolution for the turbulent transport of a passive scalar[END_REF][START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF]. In the latter reference a compact finite-difference method was used for the scalar and a pseudo-spectral method for the flow. A significant speed-up over a pure spectral solver with high resolution for both the momentum and the scalar was obtained. Our choice here is to combine a particle method for the scalar and a spectral method for the flow.

Our motivation to choose a particle method for the scalar advection comes form the fact that, for large Schmidt numbers, the scalar dynamics is essentially driven by advection, a regime for which Lagrangian or semi-Lagrangian methods are well suited. Moreover in such method, the stability limits for the time-step are governed by the amount of strain in the flow, and not by the grid size. In the present context where high resolutions of the scalar are desired, this is definitely a feature that is expected lead to an important speed up.

More precisely, the method we use for the scalar is a semi-Lagrangian (or remeshed) particle method, where at every time step particles carrying the scalar values are moved along the streamlines of the velocity then remeshed on a regular cartesian grid. Remeshing particles on a regular grid is a way to guarantee the accuracy of particle methods. This approach has been systematically used and validated in a number of simulation of vortex flows [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | A comparison of spectral and vortex methods in three-dimensional incompressible flows[END_REF][START_REF] Cottet | Advances in direct numerical simulations of 3d wall-bounded flows by vortex-in-cell methods[END_REF][START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Application to the sphere at Re=300, 500 and 1000[END_REF][START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF] or in combination with level set methods for interface capturing [START_REF] Hieber | A lagrangian particle level set method[END_REF][START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF][START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF]. Remeshing particles at every time-step also allows to easily couple the method to grid based methods, in particular when velocity values are computed on a grid. These methods can be summarized by the following formula

θ n+1 i = ∑ j θ n j Γ x n+1 j -x i ∆ x θ . ( 3
)
where θ n j denotes the value of the scalar at the grid point x j and at time t n = n∆t θ , x n+1 j is the location after one advection step of the particle initialized at time t n on the grid point x j , and ∆ x θ and ∆t θ denote the grid size and the time-step. In the above formula Γ is an interpolating kernel, the smoothness and the moment properties of which govern the spatial overall accuracy of the method [START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF]. In this work we chose the following kernel second order kernel

Γ (x) =            1 12 (1 -|x|) 25|x| 4 -38|x| 3 -3|x| 2 + 12|x| + 12 if 0 ≤ |x| < 1 1 24 (|x| -1) (|x| -2) 25|x| 3 -114|x| 2 + 153|x| -48 if 1 ≤ |x| < 2 1 24 (3 -|x|) 3 (5|x| -8) (|x| -2) if 2 ≤ |x| < 3 0 if 3 ≤ |x|.
The scalar time-step is given by ∆t θ = (|∇u| max ) -1 . As already mentioned this value does not depend on the scalar grid size.

For the momentum equation we use a classical pseudo-spectral method, with the 3/2 rule to de-alias inertial terms and a second-order Runge-Kutta scheme is used both for the time-stepping of the spectral method and to advect particles. Precise descriptions of the methods and of the experimental set up are given in [START_REF] Lagaert | Hybrid spectral-particle method for the turbulent transport of a passive scalar[END_REF][START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF].

Figure 4 shows a comparison of the scalar spectra obtained by the present coupling method and pure spectral method in an experiment of decaying homogeneous turbulence. In the hybrid method two different resolutions where used for the scalar. This experiment shows that, provided the particle method is used with slightly more grid points than needed by the spectral method, the scalar values are well recovered all the way to the dissipation scale. In this Direct Numerical Simulation, the Schmidt number was equal to 50 and the momentum equation was solved with 256 modes in each direction. To evaluate the efficiency of the hybrid method, we show in Table 1 CPU times for the full spectral method and the hybrid method for Sc = 50. All runs correspond to fully resolved simulations for the Navier-Stokes equations. One can see that, because it can use much larger time-steps, the hybrid method, even when it uses slightly more points to accurately resolve the finest scales, leads to significant savings over the pure spectral method. Additional validation and diagnostics are give in [START_REF] Lagaert | Hybrid spectral-particle method for the turbulent transport of a passive scalar[END_REF].

The computational efficiency of the hybrid method allows to address more challenging cases and to investigate in a systematic fashion the universal scaling laws in the case of forced homogeneous turbulence. Table 2 summarizes the simulation set up corresponding to two values of the Reynolds number and several Schmidt numbers. For the highest Schmidt number, Sc = 128, the simulation used 3064 3 computational elements for the scalar equation, on a IBM Blue Gene supercomputer. The ratio between the time-step used in the present simulation and that which would have been required in a comparable spectral simulation is almost equal to 100. Figure 5 shows the compensated spectra of the scalar for a Reynolds number based on the Taylor mico-scale R λ [START_REF] Lesieur | Turbulence in fluids, Fluid mechanics and its applications[END_REF] equal to 130. These spectra do exhibit a k -1 decay on a range which increases with the Schmidt number. Beyond this viscous-convective range, the spectra follow an exponential decay coinciding with the scaling law proposed by Kraichnan [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF].

R
210 512 3 3e -3 0.7 770 3 2e -2 2e -3 4 1024 3 1.5e -3
Table 2 Setup of simulations performed in forced homogeneous turbulence. ∆t u is the time step used to solve the Navier-Stokes equation with a pseudo-spectral solver. ∆t θ is the time step used to solve the scalar transport equation with the particle method. ∆t θ spec is the time step which would be needed if a pseudo-spectral method was used for the same number of scalar grid points [START_REF] Lagaert | Hybrid spectral-particle method for the turbulent transport of a passive scalar[END_REF] Figure 6 shows vorticity and scalar contours in a cross section of the computational box for the simulation corresponding R λ ≈ 130 and Sc = 128. It illustrates the scale separation between the flow and the scalar for these parameters. The extension of the hybrid method to the coupling of particle methods with finite-volume methods to address engineering configurations is under way.

Towards hybrid computing

As already mentioned the multi-scale nature of turbulent transport makes natural the idea of hybrid computing methodologies where different part of the problems are distributed to different types of hardware. To be able to implement hybrid algorithms on hybrid architectures, one needs to develop frameworks and libraries with a high level description which allows to distribute different solvers and grids to different parts of the clusters in a seamless fashion. Both particle advection and particle remeshing are local operations. This limits the communications between computational elements and makes semi-Lagrangian particle methods well suited to parallel implementations [START_REF] Sbalzarini | PPM -A highly efficient parallel particle-mesh library for the simulation of continuum systems[END_REF]. Such an implementation is described in [START_REF] Rossinelli | GPU accelerated simulations of bluff body flows using vortex methods[END_REF] for the 2D Navier-Stokes equations and in [START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF] for 3D linear transport equations. In [START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF], to achieve good portability, the computational frameworks are written using OpenCL.

The efficiency of GPU algorithms is very much conditioned by memory access strategies. To minimize the resulting computational overhead, we use a directional splitting where particles are pushed and remeshed successively along each direction. This allows to send a given number of independent particle lines on a single workgroup. This strategy requires to transpose data after each direction has been processes. However on modern GPU cards, transpositions can be achieved at a cost close to that of a simple copy operation. Figure 5 shows the computational cost of our GPU implementations [START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF] in double precision arithmetics for different remeshing kernels, for 2D and 3D experiments using about 16 million points. In these experiments a second order splitting was used to alternate one dimensional particle advection and remeshing. The number of points in the kernel stencils in each direction varied from 4 to 8 [START_REF] Cottet | High-order semi-lagrangian particle methods for transport equation: numerical analysis and implementation issues[END_REF]. These calculations were done on a NVIDIA Tesla K20m. These performances reached between 20 and 50 % of the peak performance of the GPU, depending on the size of the stencil and represented a speed up of about 25 over a multi-threaded MPI implementation running on 8 Xeon E5-2640 cores Hybrid computing would consist of combining the above implementation of scalar transport at high resolution with flow calculations on CPU processors. Based on timing obtained in our CPU and GPU implementations, in the case when the full scalar grid fits on a single GPU, up to resolutions of 512 3 , a target toy configuration where computational times on CPU and GPU would be similar, consist of a 128 3 flow resolution running on 8 CPUs together with a 512 3 scalar resolution running on the GPU, for an overall computational time of about 1s per iteration. To obtain such performance it is essential that communications between velocity data processed on the CPUs and the GPU are processed in an optimal way. This is the object of ongoing research.

Conclusion

Combining high order semi-Lagrangian and Eulerian methods is an efficient strategy to address turbulent transport problems. It allows to describe accurately the viscous-convective range and dissipation scales of the scalar at a minimal cost. This is due to the fact that semi-Lagrangian methods are not subject to CFL conditions. The local nature of particle methods naturally opens the way to hybrid computations using heterogeneous hardware. putations reported in section 3 were performed using HPC resources from GENCI-IDRIS (Grant 2012-020611).

  G. Balarac LEGI, CNRS and Université de GrenobleG.-H. Cottet, J.-M. Etancelin, F. Perignon and C. Picard Laboratoire Jean Kuntzmann, CNRS and Université de Grenoble J.-B. Lagaert Laboratoire de Mathématiques, Université Paris 11 1

Fig. 1 Fig. 2

 12 Fig. 1 Transport of a pollutant in the bay of Los Angeles at two different times. Courtesy of E. Blayo (Université Joseph Fourier, Grenoble)

Fig. 3

 3 Fig. 3 Reacting jet . Courtesy of L. Vervisch, INSA Rouen.

Fig. 4

 4 Fig. 4 Spectra of the scalar variance E θ (k,t) at two two different times for Sc = 50. Right picture is a zoom of left picture on the smallest scales.

Fig. 5

 5 Fig. 5 Compensated spectra for the scalar variance at R λ ≈ 130. The arrow shows the direction of increasing Schmidt numbers. In the dissipative region, the circles show the law proposed by Kraichnan and the squares show the law proposed by Batchelor in the dissipative scales. The vertical axis shows the spectra compensated by the Batchelor law predicting a k -1 decay in the intermediate scale.

Fig. 6

 6 Fig. 6 Cross-section colored by the vorticity magnitude (left, blue regions are for the lowest vorticity values and red regions are for the highest vorticity values) and by the passive scalar (middle, blue regions are for the lowest scalar values and red regions are for the highest scalar values) for R λ ≈ 130 and Sc = 128. The zooms (right) for the vorticity magnitude (top) and the scalar (bottom) correspond to the white box with a length approximately equal to the Kolmogorov scale.

Fig. 7

 7 Fig. 7 Code profiling of one time step for different problem sizes in double precision.

Table 1

 1 Numerical efficiency of the different methods on a decaying turbulence experiment-Runs are performed on 2048 cores of a Blue Gene Q. N u , N θ denotes the spatial resolution for velocity and scalar and ∆t u , ∆t θ are the numerical time steps for momentum and scalar equations. CPU times correspond to the simulation time t = 6.

	Method	N u	N θ	∆t u (×10 -4 )	∆t θ (×10 -4 )	total CPU time
	Spectral	1024 3	1024 3	2.5	2.5	43 590s
	Spectral	256 3	1024 3	2.5	2.5	16 671s
	Hybrid	256 3	1024 3	10	100	1 139s
	Hybrid	256 3	1280 3	10	100	1 328s

Acknowledgments

This work was partially supported by the Agence Nationale pour la Recherche (ANR) under Contracts No. ANR-2010-JCJC-091601 and ANR-2010-COSI-0009. G.-H. C. is also grateful for the support from Institut Universitaire de France. Com-


Sebastien Blandin
email: sblandin@sg.ibm.com

Paola Goatin
email: paola.goatin@inria.fr

Well-posedness of a conservation law with non-local flux arising in traffic flow modeling

Keywords: Scalar conservation laws, Non-local flux, Lax-Friedrichs scheme, Macroscopic Traffic Models

We prove the existence and stability of entropy weak solutions of a scalar conservation law with non-local flux arising in traffic flow modeling. The result is obtained providing accurate L ∞ , BV and L 1 estimates for the sequence of approximate solutions constructed by an adapted Lax-Friedrichs scheme.

Introduction

Conservation laws with non-local terms arise in a variety of physical applications. Space-integral terms are considered for example in models for granular flows [START_REF] Amadori | An integro-differential conservation law arising in a model of granular flow[END_REF], sedimentation [6], crowd motion [START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF], or more general problems like gradient constrained equations [START_REF] Amorim | On a nonlocal hyperbolic conservation law arising from a gradient constraint problem[END_REF]. Also, non-local in time terms arise in conservation laws with memory, starting from [START_REF] Dafermos | Solutions in L ∞ for a conservation law with memory[END_REF].

Macroscopic traffic flow models usually consist of one or two first-order hyperbolic partial differential equations accounting for the conservation of the number of cars, as in the celebrated Lighthill-Whitham-Richards model [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF], and for a momentum conservation or balance equation, see for example the Aw-Rascle model [4], the general phase transition model [7], or the generalized second order model [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF]. The resulting (system of) classical conservation laws usually turns out to be well-posed.

Equations with non-local flux have been recently introduced in traffic flow modeling to account for the reaction of drivers or pedestrians to the surrounding density of other individuals, see [3,[START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF][START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF][START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF][START_REF] Sopasakis | Stochastic modeling and simulation of traffic flow: asymmetric single exclusion process with Arrhenius look-ahead dynamics[END_REF]. While pedestrians are likely to react to the presence of people all around them, drivers will mainly adapt their velocity to the downstream traffic, assigning a greater importance to closer vehicles. In this paper, we consider the following mass conservation equation for traffic flow with non-local mean velocity:

∂ t ρ(t, x) + ∂ x ρ(t, x)v x+η x ρ(t, y)w η (y -x) dy = 0 , (1)
defined for t ∈ R + and x ∈ R, η > 0. Above, the convolution kernel w η ∈ C 1 ([0, η]; R + ) is a non-increasing function such that η 0 w η (x)dx = 1 (for example, w η (x) ≡ 1/η or w η (x) = 2(1 -x/η)/η). With slight abuse of notation, we define the downstream convolution product as ρ * d w η (t, x) := x+η x ρ(t, y)w η (y -x) dy .

(

) 2
This choice is intended to model the behavior of drivers reacting to what happens in front of them, thus adapting their velocity with respect to the downstream density.

For the sake of simplicity, and aiming to get sharp estimates, in this paper we take the mean velocity function to be v(r) = 1 -r. The same approach could be applied to more general continuous decreasing functions.

Setting V (t, x) = v(ρ * d w η (t, x)), we rewrite (1) as ∂ t ρ(t, x) + ∂ x (ρ(t, x)V (t, x)) = 0 , and we couple it with an initial datum ρ(0, x) = ρ 0 (x) ∈ BV(R; [0, 1]).

(3)

We will consider solutions ρ = ρ(t, x) satisfying the following definition.

Definition 1 A function ρ ∈ (L 1 ∩L ∞ ∩BV)(R + ×R; [0, 1]
) is a weak solution of (1), (3), if

+∞ 0 +∞ -∞ (ρϕ t + ρ(t, x)v(ρ * d w η )ϕ x ) dxdt + +∞ -∞
ρ 0 (x)ϕ(0, x) dx = 0 (4) for all ϕ ∈ C 1 c (R 2 ; R).

Following [16, Definition 1], [6, Definition 4.1] and [10, Definition 2.1], we can also give an entropy criterion to select admissible solutions.

Definition 2 A function ρ ∈ (L 1 ∩ L ∞ ∩ BV)(R + × R; [0, 1]) is an entropy weak solution if +∞ 0 +∞ -∞ (|ρ -κ|ϕ t + |ρ -κ| V ϕ x -sgn(ρ -κ) k V x ϕ) (t, x)dxdt + +∞ -∞ |ρ 0 (x) -κ|ϕ(0, x) dx ≥ 0 (5) for all ϕ ∈ C 1 c (R 2 ; R + ) and κ ∈ R.
We remark that, unlike [6, 10], in our case the assumptions of Kružkov are not fully satisfied: the convolution kernel w η is not continuous on R. Therefore, uniqueness of the above entropy weak solutions is not guaranteed.

The main results of this paper are collected in the following theorem.

Theorem 1 Let ρ 0 ∈ BV(R; [0, 1]) and w η ∈ C 1 ([0, η]; R + ) be a non-increasing function such that η 0 w η (x)dx = 1.
Then the Cauchy problem

∂ t ρ + ∂ x (ρv(ρ * d w η )) = 0, x ∈ R, t > 0, ρ(0, x) = ρ 0 (x), x ∈ R,
admits a weak entropy solution in the sense of Definitions 1 and 2, such that

min R {ρ 0 } ≤ ρ(t, x) ≤ max R {ρ 0 }, for a.e. x ∈ R, t > 0. (6)
The existence of a weak entropy solution is proved by constructing a converging sequence of finite volume approximate solutions, defined using an adapted Lax-Friedrichs scheme. Compared to the previous studies [3, 6, 9-11], our result does not require the kernel w η to be smooth, and provides very accurate L ∞ bounds on the solutions (which are given by the maximum principle (6) above). Moreover, equation ( 1) provides a meaningful model of road traffic evolution. Yet, we can not provide a uniqueness result. Nevertheless, stability estimates issued from the numerical approximation guarantee that the solutions obtained as limit of our finite volume scheme are stable with respect to the initial data in the L 1 norm, see Section 2.4.

The paper is organized as follows. Section 2 contains the description of the finite volume scheme used to construct approximate solutions, and the proofs of its fine properties: maximum principle, bounded total variation, discrete entropy inequalities and L 1 stability estimates. Relying on these results, the proof of Theorem 1 is detailed in Section 3. Some numerical simulations illustrating the properties of the solutions are collected in Section 4.

A Lax-Friedrichs numerical scheme

We take a space step ∆x such that η = N ∆x, for some N ∈ N, and a time step ∆t subject to a CFL condition which will be specified later. For j ∈ Z and n ∈ N, let x j+1/2 = j∆x be the cells interfaces, x j = (j -1/2)∆x the cells centers and t n = n∆t the time mesh. We want to construct a finite volume approximate solution ρ ∆x (t, x) = ρ n j for (t, x) ∈ C n j = [t n , t n+1 [×]x j-1/2 , x j+1/2 [. Let us denote w k η := w η (k∆x) for k = 0, . . . , N -1 and set

V j := 1 -∆x N -1 k=0 w k η ρ j+k , (7)
which involves a quadrature formula to approximate the convolution term. Remark that the above discretization choice for w η implies ∆x

N -1 k=0 w k η ≤ 1 + w η (0)∆x. ( 8
)
We consider the following modified Lax-Friedrichs flux adapted to (1):

F n j+1/2 :=g(ρ n j , . . . , ρ n j+N ) = 1 2 ρ n j V n j + 1 2 ρ n j+1 V n j+1 + α 2 (ρ n j -ρ n j+1 ) , (9)
where α ≥ 1 is the viscosity coefficient. This gives the following N + 2 points finite volume scheme

ρ n+1 j = H(ρ n j-1 , . . . , ρ n j+N ) , (10)
where H(ρ j-1 , . . . , ρ j+N ) := ρ j -λ (g(ρ j , . . . , ρ j+N ) -g(ρ j-1 , . . . , ρ j+N -1 ))

=

ρ j + λα 2 (ρ j-1 -2ρ j + ρ j+1 ) + λ 2 (ρ j-1 V j-1 -ρ j+1 V j+1 ) , with λ = ∆t/∆x.
It is straightforward to show that H(ρ, . . . , ρ) = ρ for all ρ ∈ [0, 1] and the numerical flux [START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF] satisfies the minimal assumptions:

-Consistency:

g(ρ, . . . , ρ) = ρ 1 -ρ ∆x N -1 k=0 w k η . (12)
In particular, g(ρ, . . . , ρ) → ρ(1 -ρ * d w η ) as ∆x ց 0 for ρ constant. -Lipschitz continuity: there exists K > 0 such that

|g(ρ j , . . . , ρ j+N ) -g(ρ, . . . , ρ)| ≤ K max 0≤k≤N |ρ j+k -ρ| for ρ, ρ j+k ∈ [0, 1], k = 0, . . . , N . (8) guarantees that this is is verified for K = 2 + α + w η (0)∆x.
Assuming ρ i ∈ [0, 1] for i = j -1, . . . , j + N , straightforward computations give:

∂H ∂ρ j-1 = λ 2 α + 1 -∆x N -1 k=0 w k η ρ j+k-1 -∆x w 0 η ρ j-1 (13a) ∂H ∂ρ j = 1 -λ α + 1 2 ∆x w 1 η ρ j-1 ≥ 1 -λ α + ∆x w η (0) 2 (13b) ∂H ∂ρ j+1 = λ 2 α -1 + ∆x N -1 k=0 w k η ρ j+k+1 -∆x w 2 η ρ j-1 + ∆x w 0 η ρ j+1 (13c) ∂H ∂ρ j+k = λ 2 ∆x w k-1 η ρ j+1 -w k+1 η ρ j-1 k = 2, . . . , N -2 (13d) ∂H ∂ρ j+N -1 = λ 2 ∆x w N -2 η ρ j+1 (13e) ∂H ∂ρ j+N = λ 2 ∆x w N -1 η ρ j+1 (13f)
Observe that (13e) and (13f) are non-negative. Moreover, the CFL condition

∆t ≤ 2 2α + ∆x w η (0) ∆x (14)
ensures the positivity of (13b) and the assumption

α ≥ 1 + ∆x w η (0) (15)
guarantees the increasing monotonicity w.r.t. ρ j+1 in (13c), and combined with [START_REF] Bretti | Numerical algorithms for simulations of a traffic model on road networks[END_REF], guarantees the non-negativity of (13a). To obtain ( 14) and ( 15) we used the fact that w k η ≤ w η (0) for all k = 0, . . . , N -1, by non-increasing monotonicity assumption. On the contrary, the sign of (13d) cannot be a-priori determined. (Nevertheless, we observe that if ρ j-1 = ρ j+1 the derivatives in (13d) are non-negative due to the decreasing monotonicity of w η .) Therefore, the numerical scheme (9), [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] is not monotone, and classical convergence results do not apply. Nevertheless, we are able to recover the necessary L ∞ and BV bounds.

Maximum principle and L ∞ estimates

The sought L ∞ bound is a direct consequence of a maximum principle property.

Proposition 1 For any initial datum ρ

0 j , j ∈ Z, let ρ m = min j∈Z {ρ 0 j } ∈ [0, 1] and ρ M = max j∈Z {ρ 0 j } ∈ [0, 1].
Then the finite volume approximation ρ n j , j ∈ Z and n ∈ N, constructed using scheme (9), [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] satisfies the bounds

ρ m ≤ ρ n j ≤ ρ M
for all j ∈ Z and n ∈ N, under the CFL condition [START_REF] Eymard | Finite volume methods[END_REF].

The proof is based on the following lemma.

Lemma 1 Let 0 ≤ ρ m ≤ ρ n j ≤ ρ M ≤ 1 for all j ∈ Z. Then H(ρ m , ρ m , ρ m , ρ j+2 , . . . , ρ j+N -2 , ρ m , ρ m ) ≥ ρ m , (16)
H(ρ M , ρ M , ρ M , ρ j+2 , . . . , ρ j+N -2 , ρ M , ρ M ) ≤ ρ M . ( 17
)
Proof From [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF] we get

H(ρ m , ρ m , ρ m , ρ j+2 , . . . , ρ j+N -2 , ρ m , ρ m ) = ρ m + ∆t 2 ρ m N -1 k=0 w k η (ρ j+k+1 -ρ j+k-1 ) ,
and we observe that

N -1 k=0 w k η (ρ j+k+1 -ρ j+k-1 ) = ρ m w N -2 η + w N -1 η -w 0 η -w 1 η + N -2 k=1 ρ j+k w k-1 η -w k+1 η ≥ ρ m w N -2 η + w N -1 η -w 0 η -w 1 η + ρ m N -2 k=1 w k-1 η -w k+1 η = ρ m N k=1 w k-1 η - N -2 k=-1 w k+1 η = 0 ,
where the inequality is due to the non-increasing monotonicity of w η . Inequality ( 17) can be recovered following the same procedure.

Proof of Proposition 1

We apply the mean value theorem between the points R n j = (ρ n j-1 , . . . , ρ n j+N ) and R n m = (ρ m , ρ m , ρ m , ρ n j+2 , . . . , ρ n j+N -2 , ρ m , ρ m ), which by ( 16) tells us

ρ n+1 j = H(R n j ) = H(R n m ) + ∇H(R ξ ), R n j -R n m ≥ ρ m + ∇H(R ξ ), R n j -R n m , (18)
for

R ξ = (1 -ξ)R n m + ξR n j , for some ξ ∈ [0, 1]. It is now enough to observe that ∂H ∂ρ j+k (R ξ )(R n j -R n m ) k = 0 k = 2, . . . , N -2, since (R n j -R n m ) k = 0 for k = 2, . . . , N -2.
Therefore, under the assumptions ( 14) and ( 15), we can conclude that ∇H(R ξ ), R n j -R n m ≥ 0 and therefore by [START_REF] Leveque | Numerical methods for conservation laws[END_REF] we have proved that

ρ n+1 j ≥ ρ m .
The upper bound ρ n+1 j ≤ ρ M is recovered smilarly by considering

R n M = (ρ M , ρ M , ρ M , ρ n j+2 , . . . , ρ n j+N -2 , ρ M , ρ M ) in place of R n
m and using (17).

BV estimates

BV estimates cannot be derived here using standard general approaches. Indeed, scheme [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF] does not fit the classical assumption of TVD schemes, see [15, Lemma 2.2], and [5] for the generalization of these conditions to Npoints schemes. Nevertheless, accurate estimates show that the approximate solutions constructed using our numerical scheme have bounded total variation and preserve monotonicity.

Proposition 2 Let ρ 0 ∈ BV(R; [0, 1]
), and let ρ ∆x be given by [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]. If α ≥ 1 + 2 ∆x w η (0) and the CFL condition ∆t ≤ 2∆x/(2α + 3∆x w η (0)) holds, then for every T > 0 the following discrete space BV estimate is satisfied

TV(ρ ∆x (T, •)) ≤ C(w η , ρ 0 , T ) := e 2wη(0)T TV(ρ 0 ). ( 19
)
In particular, the numerical scheme [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], ( 11) is monotonicity preserving.

Proof Scheme ( 10), [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], can be rewritten as

ρ n+1 j = ρ n j - λ 2 α + 1 -∆x N -1 k=0 w k η ρ n j+k+1 -∆x w 0 η ρ n j-1 ∆ n j-1/2 + λ 2 α -1 + ∆x N -1 k=0 w k η ρ n j+k+1 + ∆x (w 0 η + w 1 η ) ρ n j-1 ∆ n j+1/2 + λ 2 ρ n j-1 ∆x N -1 k=2 (w k-1 η + w k η ) ∆ n j+k-1/2 + λ 2 ∆x w N -1 η ρ n j-1 ∆ n j+N -1/2 ,
where we have set ∆ n j+k-1/2 = ρ n j+k -ρ n j+k-1 for k = 0, . . . , N . In the same way we get

ρ n+1 j+1 = ρ n j+1 - λ 2 α + 1 -∆x N -1 k=0 w k η ρ n j+k+2 -∆x w 0 η ρ n j ∆ n j+1/2 + λ 2 α -1 + ∆x N -1 k=0 w k η ρ n j+k+2 + ∆x (w 0 η + w 1 η ) ρ n j ∆ n j+3/2 + λ 2 ρ n j ∆x N -1 k=2 (w k-1 η + w k η ) ∆ n j+k+1/2 + λ 2 ∆x w N -1 η ρ n j ∆ n j+N +1/2 .
Therefore computing the difference gives

∆ n+1 j+1/2 = = λ 2 α + 1 -∆x N -1 k=0 w k η ρ n j+k+1 -∆x w 0 η ρ n j-1 ∆ n j-1/2 + 1 - λ 2 2α -∆x w 0 η ρ n j + ∆x (w 0 η + w 1 η )ρ n j-1 -∆x N -1 k=0 w k η ∆ j+k+3/2 ∆ n j+1/2 + λ 2 α -1 + ∆x (w 0 η + w 1 η ) ρ n j -∆x (w 1 η + w 2 η ) ρ n j-1 + ∆x N -1 k=0 w k η ρ n j+k+2 ∆ n j+3/2 + λ 2 ∆x N -2 k=2 [(w k-1 η + w k η )ρ n j -(w k η + w k+1 η )ρ n j-1 ] ∆ n j+k+1/2 + λ 2 ∆x [(w N -2 η + w N -1 η )ρ n j -w N -1 η ρ n j-1 ]∆ n j+N -1/2 + λ 2 ∆x w N -1 η ρ n j ∆ n j+N +1/2 .
Adding and subtracting (w k η + w k+1 η )ρ n j in the fourth term of the above righthand side, and noting that

[(w N -2 η +w N -1 η )ρ n j -w N -1 η ρ n j-1 ]∆ n j+N -1/2 = [w N -2 η ρ n j +w N -1 η ∆ n j-1/2 ]∆ n j+N -1/2 , we get ∆ n+1 j+1/2 = = λ 2 α + 1 -∆x w 0 η (ρ n j-1 + ρ n j+1 ) -∆x (w 1 η + w 2 η ) ρ n j+2 -∆x N -2 k=2 w k+1 η ρ n j+k ∆ n j-1/2 + 1 - λ 2 2α -∆x w 0 η ρ n j + ∆x (w 0 η + w 1 η ) ρ n j-1 -∆x N -1 k=0 w k η ∆ j+k+3/2 ∆ n j+1/2 + λ 2 α -1 + ∆x (w 0 η + w 1 η ) ρ n j -∆x (w 1 η + w 2 η ) ρ n j-1 + ∆x N -1 k=0 w k η ρ n j+k+2 ∆ n j+3/2 + λ 2 ρ n j ∆x N -2 k=2 (w k-1 η -w k+1 η ) ∆ n j+k+1/2 + λ 2 ∆x w N -2 η ρ n j ∆ n j+N -1/2 + λ 2 ∆x w N -1 η ρ n j ∆ n j+N +1/2 .
Observe that the first coefficient in the summation is non-negative for ∆x sufficiently small such that α ≥ 2w η (0)∆x, since

∆x w 0 η (ρ n j-1 + ρ n j+1 ) + ∆x(w 1 η + w 2 η )ρ n j+2 + ∆x N -2 k=2 w k+1 η ρ n j+k ≤ ∆x w 0 η + ∆x N -1 k=0 w k η ≤ 1 + 2w η (0)∆x
by [START_REF] Bretti | Numerical algorithms for simulations of a traffic model on road networks[END_REF]. The second coefficient in the summation is non-negative under the slightly stronger CFL assumption

∆t ≤ 2 2α + 3∆x w η (0) ∆x,
and the third term is non-negative if α ≥ 1 + 2 ∆x w η (0). Hence all the above coefficients are non-negative, and the above formula guarantees in particular that the scheme (10), ( 11) is monotonicity preserving. Taking the absolute values in the above expression and summing over j ∈ Z we get

j ∆ n+1 j+1/2 = = j λ 2 α + 1 -∆x w 0 η (ρ n j-1 + ρ n j+1 ) -∆x (w 1 η + w 2 η ) ρ n j+2 -∆x N -2 k=2 w k+1 η ρ n j+k ∆ n j-1/2 + j 1 - λ 2 2α -∆x w 0 η ρ n j + ∆x (w 0 η + w 1 η ) ρ n j-1 -∆x N -1 k=0 w k η ∆ j+k+3/2 ∆ n j+1/2 + j λ 2 α -1 + ∆x (w 0 η + w 1 η ) ρ n j -∆x (w 1 η + w 2 η ) ρ n j-1 + ∆x N -1 k=0 w k η ρ n j+k+2 ∆ n j+3/2 + j λ 2 ρ n j ∆x N -2 k=2 (w k-1 η -w k+1 η ) ∆ n j+k+1/2 + j λ 2 ∆x w N -2 η ρ n j ∆ n j+N -1/2 + j λ 2 ∆x w N -1 η ρ n j ∆ n j+N +1/2 .
Rearranging the indices we obtain

j ∆ n+1 j+1/2 = j ∆ n j+1/2 1 + ∆t 2 N -2 k=2 (w k-1 η -w k+1 η )ρ n j+k+1 -(w 1 η + w 2 η )ρ n j-2 + N -2 k=2 (w k-1 η -w k+1 η )ρ n j-k -(w 1 η + w 2 η )ρ n j+3 + w N -2 η ρ n j+N + w N -1 η ρ n j+N +1 + w N -2 η ρ n j-N -1 + w N -1 η ρ n j-N ≤ 1 + ∆t 2 2 N -2 k=2 (w k-1 η -w k+1 η ) + 2w N -2 η + 2w N -1 η j ∆ n j+1/2 ≤ (1 + 2w η (0)∆t) j ∆ n j+1/2 .
Therefore we recover the following estimate for the total variation TV(ρ ∆x (T, •)) ≤ (1 + 2w η (0)∆t) T /∆t TV(ρ ∆x (0, •)) ≤ e 2wη(0)T TV(ρ 0 ) .

From Proposition 2, the following space-time BV estimate can be derived as in [START_REF] Eymard | Finite volume methods[END_REF]Corollary 5.1], relying on the Lipschitz continuity of g .

Corollary 1 Let ρ 0 ∈ BV(R; [0, 1]), and let ρ ∆x be given by [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]. If α ≥ 1 + 2 ∆x w η (0) and ∆t ≤ 2∆x/(2α + 3∆x w η (0)), then for every T > 0 there exists C = C(w η , ρ 0 , T, α) such that

TV(ρ ∆x ; [0, T ] × R) ≤ C. ( 20
)

Discrete entropy inequalities

Following [3, Proposition 2.8], we derive a discrete entropy inequality for the approximate solution generate by [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]. Let us denote by Proposition 3 Let ρ n j , j ∈ Z, n ∈ N, be given by ( 9), ( 10) Then, if α ≥ 1 and the CFL condition [START_REF] Eymard | Finite volume methods[END_REF] holds, for all j ∈ Z, n ∈ N we have

G j+1/2 (u, v) = 1 2 u V n j + 1 2 v V n j+1 + α 2 (u -v) , F κ j+1/2 (u, v) = G j+1/2 (u ∧ κ, v ∧ κ) -G j+1/2 (u ∨ κ, v ∨ κ),
ρ n+1 j -κ -ρ n j -κ + λ F κ j+1/2 (ρ n j , ρ n j+1 ) -F κ j-1/2 (ρ n j-1 , ρ n j ) (21)
+ λ 2 sgn(ρ n+1 j -κ) κ V n j+1 -V n j-1 ≤ 0 for all κ ∈ R. Proof Setting Hj (u, v, z) = v -λ G j+1/2 (v, z) -G j-1/2 (u, v) ,
the function Hj is monotone non-decreasing in its first variable, monotone non-decreasing in its second variable for αλ ≤ 1, which is guaranteed by the CFL condition ( 14), and monotone non-decreasing in its third variable for α ≥ 1, which is guaranteed by [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]. Moreover, we have the identity

Hj (ρ n j-1 ∧ κ, ρ n j ∧ κ, ρ n j+1 ∧ κ) -Hj (ρ n j-1 ∨ κ, ρ n j ∨ κ, ρ n j+1 ∨ κ) = ρ n j -κ -λ F κ j+1/2 (ρ n j , ρ n j+1 ) -F κ j-1/2 (ρ n j-1 , ρ n j ) .
By monotonicity,

Hj (ρ n j-1 ∧ κ, ρ n j ∧ κ, ρ n j+1 ∧ κ) -Hj (ρ n j-1 ∨ κ, ρ n j ∨ κ, ρ n j+1 ∨ κ) = Hj (ρ n j-1 , ρ n j , ρ n j+1 ) ∧ Hj (κ, κ, κ) -Hj (ρ n j-1 , ρ n j , ρ n j+1 ) ∨ Hj (κ, κ, κ) = Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -Hj (κ, κ, κ) = sgn Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -Hj (κ, κ, κ) × Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -Hj (κ, κ, κ) = sgn Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ + λ 2 κ V n j+1 -V n j-1 × Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ + λ 2 κ V n j+1 -V n j-1 ≥ sgn Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ × Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ + λ 2 κ V n j+1 -V n j-1 = Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ + λ 2 sgn Hj (ρ n j-1 , ρ n j , ρ n j+1 ) -κ κ V n j+1 -V n j-1 = ρ n+1 j -κ + λ 2 sgn ρ n+1 j -κ κ V n j+1 -V n j-1 ,
by definition of the scheme ( 10), [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], which gives (21).

L 1 stability estimates

Since Kružkov theory cannot be applied as it is to our case, because of lack of regularity of the convolution kernel w η , we prove explicit L 1 estimates that ensure the stability of the scheme (10), [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF].

Proposition 4 Let ρ 0 , ρ0 ∈ BV(R; [0, 1]) be two initial data, and denote by ρ ∆x , ρ∆x the corresponding approximate solutions constructed applying the modified Lax-Friedrichs scheme [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]:

ρ n+1 j = ρ n j + λα 2 ρ n j-1 -2ρ n j + ρ n j+1 + λ 2 ρ n j-1 V n j-1 -ρ j+1 V n j+1 , (22) ρn+1 j
= ρn j + λα 2 ρn j-1 -2ρ n j + ρn j+1 + λ 2 ρn j-1 V n j-1 -ρj+1 V n j+1 , ( 23
)
where we have set

V n j = 1-∆x N -1 k=0 w k η ρ n j+k and V n j = 1-∆x N -1 k=0 w k η ρn j+k .
Then, under assumptions [START_REF] Eymard | Finite volume methods[END_REF] and [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF], the following estimate holds:

ρ ∆x (T, •) -ρ∆x (T, •) L 1 ≤ K(w η , ρ 0 , ρ0 , T ) ρ 0 -ρ0 L 1 (24)
with K(w η , ρ 0 , ρ0 , T ) := exp(T w η (0)(1+0.5 min{C(w η , ρ 0 , T ), C(w η , ρ0 , T )})).

Proof Subtracting (23) from ( 22) we get

ρ n+1 j -ρn+1 j = = (1 -λα)(ρ n j -ρn j ) + λα 2 (ρ n j-1 -ρn j-1 ) + λα 2 (ρ n j+1 -ρn j+1 ) + λ 2 ρ n j-1 (V n j-1 -V n j-1 ) + (ρ n j-1 -ρn j-1 ) V n j-1 - λ 2 (ρ n j+1 -ρn j+1 ) V n j+1 + ρ n j+1 (V n j+1 -V n j+1 ) = 1 -λα - λ 2 ∆x w 1 η ρ n j-1 (ρ n j -ρn j ) + λ 2 (α + V n j-1 -∆x w 0 η ρ n j-1 )(ρ n j-1 -ρn j-1 ) + λ 2 (α -V n j+1 + ∆x w 0 η ρ n j+1 -∆x w 2 η ρ n j-1 )(ρ n j+1 -ρn j+1 ) + λ 2 ∆x N -2 k=2 w k-1 η ρ n j+1 -w k+1 η ρ n j-1 (ρ n j+k -ρn j+k ) + λ 2 ∆x w N -2 η ρ n j+1 (ρ n j+N -1 -ρn j+N -1 ) + λ 2 ∆x w N -1 η ρ n j+1 (ρ n j+N -ρn j+N ).
Observe that the coefficient of the first term is positive thanks to ( 14) and the coefficients of the second and third terms are positive thanks to [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]. Therefore, taking the absolute values in the above equality we get

ρ n+1 j -ρn+1 j ≤ ≤ 1 -λα - λ 2 ∆x w 1 η ρ n j-1 ρ n j -ρn j + λ 2 (α + V n j-1 -∆x w 0 η ρ n j-1 ) ρ n j-1 -ρn j-1 + λ 2 (α -V n j+1 + ∆x w 0 η ρ n j+1 -∆x w 2 η ρ n j-1 ) ρ n j+1 -ρn j+1 + λ 2 ∆x N -2 k=2 w k-1 η ρ n j+1 -w k+1 η ρ n j-1 ρ n j+k -ρn j+k + λ 2 ∆x w N -2 η ρ n j+1 ρ n j+N -1 -ρn j+N -1 + λ 2 ∆x w N -1 η ρ n j+1 ρ n j+N -ρn j+N .
Summing over j ∈ Z, rearranging the indexes and observing that by monotonicity of w η and using the triangular inequality

w k-1 η ρ n j+1 -w k+1 η ρ n j-1 ≤ (w k-1 η -w k+1 η )ρ n j+1 + w k+1 η ρ n j+1 -ρ n j-1 we get j ρ n+1 j -ρn+1 j ≤ 1 + ∆t 2 N -2 k=2 (w k-1 η -w k+1 η ) + w N -2 η + w N -1 η + ∆t 2 N -2 k=2 w k+1 η ρ n j-k+1 -ρ n j-k-1 j ρ n j -ρn j .
Therefore,

ρ ∆x (T, •) -ρ∆x (T, •) L 1 ≤ 1 + ∆t 2 + TV(ρ ∆x (T )) 2 w η (0)) T /∆t ρ 0 -ρ0 L 1 ,
which gives the desired estimate (24).

Proof of Theorem 1

We apply the classical procedure of Lax-Wendroff theorem, see [START_REF] Eymard | Finite volume methods[END_REF]Theorem 5.3] and [START_REF] Leveque | Numerical methods for conservation laws[END_REF]Theorem 12.1], to show that the scheme ( 10), [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], converges to a weak solution of (1), in the sense of Definition 1.

By Proposition 1 and Corollary 1, we can apply Helly's theorem, stating that there exists a subsequence, still denoted by ρ ∆x , that converges to some

ρ ∈ (L 1 ∩ L ∞ ∩ BV)(R + × R; [0, 1]) in the L 1 loc -norm. Let ϕ ∈ C 1 c (R 2
) and multiply [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] by ϕ(t n , x j ). Summing over j ∈ Z and n ∈ N we get

n j ϕ(t n , x j ) ρ n+1 j -ρ n j = -λ n j
ϕ(t n , x j ) g(ρ n j , . . . , ρ n j+N ) -g(ρ n j-1 , . . . , ρ n j+N -1 ) .

Summing by parts we obtain

j ϕ(0, x j )ρ 0 j + n j ϕ(t n , x j ) -ϕ(t n-1 , x j ) ρ n j (25) + λ n j (ϕ(t n , x j+1 ) -ϕ(t n , x j )) g(ρ n j , . . . , ρ n j+N ) = 0.
Then we multiply (25) by ∆x getting

∆x j ϕ(0, x j )ρ 0 j + ∆x∆t n j ϕ(t n , x j ) -ϕ(t n-1 , x j ) ∆t ρ n j (26) + ∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x g(ρ n j , . . . , ρ n j+N ) = 0.
By strong L 1 loc convergence of ρ ∆x → ρ, it is straightforward to see that the first two terms in (26) converge to

+∞ -∞ ρ 0 (x)ϕ(0, x) dx + +∞ 0 +∞ -∞ ρ(t, x)ϕ t (t, x) dxdt
as ∆x ց 0. Concerning the last term, since ρ n j ∈ [0, 1] we observe that

g(ρ n j , . . . , ρ n j+N ) -ρ n j V n j ≤ α 2 ρ n j+1 -ρ n j + 1 2 ρ n j+1 V n j+1 -ρ n j V n j ≤ α 2 ρ n j+1 -ρ n j + 1 2 (ρ n j+1 -ρ n j )V n j+1 + ρ n j (V n j+1 -V n j ) ≤ 1 + α 2 ρ n j+1 -ρ n j + 1 2 ρ n j ∆x N -1 k=0 w k η ρ n j+k+1 -ρ n j+k ≤ 1 + α 2 ρ n j+1 -ρ n j + 1 2 w η (0) TV(ρ ∆x (t n , •)) ∆x ≤ 1 + α 2 ρ n j+1 -ρ n j + C ′ (w η , ρ 0 , T ) ∆x , (27)
where we have set C ′ (w η , ρ 0 , T ) = w η (0)C(w η , ρ 0 , T )/2 for T ≥ t n . Therefore, the last term in (26) can be rewritten as

∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x g(ρ n j , . . . , ρ n j+N ) = ∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x ρ n j V n j + ∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x × g(ρ n j , . . . , ρ n j+N ) -ρ n j V n j .
By consistency, the first term in the above decomposition converges to

+∞ 0 +∞ -∞ ρ(t, x)v(ρ * d w η (t, x))ϕ x (t, x) dxdt ,
while the second term can be bounded using (27): set T > 0 and R > 0 such that ϕ(t, x) = 0 for t > T and |x| > R, and let n T ∈ N and j 0 , j

1 ∈ Z such that T ∈ ]n T ∆t, (n T + 1)∆t], -R ∈ ]x j0-1/2 , x j0+1/2 ] and R ∈ ]x j1-1/2 , x j1+1/2 ], then ∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x g(ρ n j , . . . , ρ n j+N ) -ρ n j V n j ≤ ∆x∆t ϕ x ∞ n T n=0 j1 j=j0 1 + α 2 ρ n j+1 -ρ n j + C ′ (w η , ρ 0 , T ) ∆x = 1 + α 2 ϕ x ∞ ∆x∆t n T n=0 j1 j=j0 ρ n j+1 -ρ n j + ϕ x ∞ C ′ (w η , ρ 0 , T )2R T ∆x ≤ 1 + α 2 ϕ x ∞ T 0 R -R |ρ ∆x (t, x + ∆x) -ρ ∆x (t, x)| dxdt + ϕ x ∞ C ′ (w η , ρ 0 , T )2R T ∆x ≤ 1 + α 2 ϕ x ∞ C(w η , ρ 0 , T )∆x + ϕ x ∞ C ′ (w η , ρ 0 , T )2R T ∆x ,
which clearly goes to zero when ∆x ց 0.

Concerning the entropy condition, we proceed as above to show that (21) converges to (5). Multiplying ( 21) by ∆x ϕ(t n , x j ) ≥ 0 and then summing by parts we get

0 ≤ ∆x j ϕ(0, x j ) ρ 0 j -κ + ∆x∆t n j ϕ(t n , x j ) -ϕ(t n-1 , x j ) ∆t ρ n j -κ + ∆x∆t n j ϕ(t n , x j+1 ) -ϕ(t n , x j ) ∆x F κ j+1/2 (ρ n j , ρ n j+1 ) -∆x∆t n j sgn ρ n+1 j -κ κ V n j+1 -V n j-1 2∆x ϕ(t n , x j ) .
Following the same steps as above, the first three terms in the sum clearly converge to

+∞ 0 +∞ -∞ (|ρ(t, x) -k|ϕ t + |ρ(t, x) -k|v(ρ * d w η (t, x)) ϕ x ) dxdt + +∞ -∞ |ρ 0 (x) -κ|ϕ(0, x) dx
as ∆x ց 0. The third term can be decomposed as

n j sgn ρ n+1 j -κ κ V n j+1 -V n j-1 2∆x ϕ(t n , x j ) = sgn ρ n+1 j -κ -sgn ρ n j -κ κ V n j+1 -V n j-1 2∆x ϕ(t n , x j ) (28) + sgn ρ n j -κ κ V n j+1 -V n j-1 2∆x ϕ(t n , x j ) .
The first term in (28) can be controlled by C(w η , ρ 0 )∆x, and the second clearly converges to

- +∞ 0 +∞ -∞ sgn(ρ(t, x) -k) k V x (t, x) ϕ(t, x)dxdt ,
providing the entropy inequality (5).

Numerical tests

In this section, we perform numerical simulations to show evidence of some properties of equation [START_REF] Amadori | An integro-differential conservation law arising in a model of granular flow[END_REF]. In particular, we compare solutions depending on the kernel choice and its support.

For the tests presented below, the space domain is given by the interval [-1, 1], the space discretization mesh ∆x = 0.002 and η = 0.1, where not specified otherwise. Absorbing conditions are imposed at the boundaries.

Kernel support location

We aim at comparing the behavior of solutions of (1) with centered and upstream convolutions, obtained taking the velocity v as

1 -∆x N/2 k=-N/2 w k η ρ j+k ≃ v x+η/2 x-η/2 ρ(t, y)w η (y -x) dy , (29)
which models drivers looking both ahead and behind, and

1 -∆x 1 k=-N w k η ρ j+k ≃ v x x-η ρ(t, y)w η (y -x) dy , (30)
corresponding to drivers taking into account traffic densities behind them.

Figure 1 shows similarly the density profiles at time t = 0.5 corresponding to the initial condition

ρ 0 (x) = 0.6 if x < 0, 0.2 if x > 0. ( 31
)
The solution to the classical scalar conservation law

∂ t ρ + ∂ x (ρ(1 -ρ)) = 0 (32)
would be a rarefaction wave, as displayed in Figure 1,(a).

Figure 2 shows the density profiles at time t = 0.2 corresponding to the initial condition The solution to the classical scalar conservation law would be a shock traveling with speed σ = -0.3. Cases (b), (c) and (d) in Figures 1,2 show the density profiles corresponding to the constant kernel w η (x) = 1/η with downstream (2), central (29) and upstream (30) supports respectively. (The approximate solutions do not change significantly using other non-increasing kernels.) We observe that the downstream convolution gives a smooth profile with constant total variation and monotonicity (as expected from Proposition 2) close to the solution of the corresponding classical problem (32) showed in Figures 1,2,(a). On the contrary, central and upstream convolutions cause oscillations formation and total variation blow-up (for upstream convolution).

ρ 0 (x) = 0.4 if x < 0, 0.9 if x > 0. ( 33
) (a) (b) (c) (d)

Kernel monotonicity

We are interested in studying the effect of the monotonicity of the kernel on the solution characteristics. To this aim, we consider the downstream convolution (2) and we compare the solution corresponding to the constant kernel w η (x) = 1/η, to the linear decreasing kernel w η (x) = 2(η -x)/η 2 and to the linear increasing kernel w η (x) = 2x/η 2 , for the Riemann-type initial datum (33). Figure 3 shows the density profiles at time t = 0.5 (left) and the total variation TV(ρ(t, •); [-1, 1]) for t ∈ [0, 1] (right) corresponding to the different kernel choices. Numerical simulations confirm that non-increasing monotonicity is necessary for the scheme to be monotonicity preserving. In particular, for Riemann type initial data, the total variation remains constant, while it increases for monotonically increasing kernels. We also compute numerical convergence orders. Following [START_REF] Bretti | Numerical algorithms for simulations of a traffic model on road networks[END_REF], we define

γ(∆x) = log 2 e(∆x) e(∆x/2) , (34)
where the L 1 -error is computed at final time T = 0.3 as where x j = (j -1/2)∆x/2 for ∆x = 0.01, 0.005, 0.0025, 0.00125, 0.000625. Table 1 show that non-increasing kernels give good convergence rates (the rate is far better for the linear decreasing kernel), while in the case of linear increasing kernel convergence is not clearly established. The L 1 -errors reported in the table represent the L 1 -distances to the reference solutions corresponding to ∆x = 0.00015625. 1 Convergence orders (34) and L 1 -errors to the reference solutions corresponding to ∆x = 0.00015625 for constant, linear decreasing and linear increasing kernels at final time T = 0.5 corresponding to a Riemann-like initial datum with ρ L = 0.4, ρ R = 0.9.

e(∆x) = ρ ∆x (T, •) -ρ ∆x/2 (T, •) L 1 = ∆x 2 j ρ ∆x (T, x j ) -ρ ∆x/2 (T, x j ) , (35)
wη(x) = 1/η wη(x) = 2(η -x)/η 2 wη(x) = 2x/η 2 ∆x γ(∆x) L 1 -error γ(∆x) L 1 -error γ(∆x) L 1 -

Limit η ց 0

In this section, we investigate the convergence to the solutions of the classical conservation law (32) as η ց 0. To this end, we compare the solutions to [START_REF] Amadori | An integro-differential conservation law arising in a model of granular flow[END_REF] with different values of η = 0.1, 0.01, 0.001, for ∆x = 0.001 and initial datum (33). As expected, the solution computed with η = 0.001 coincides with the numerical solution of (1) computed using the Lax-Friedrichs type scheme

ρ n+1 j = ρ n j + λα 2 ρ n j-1 -2ρ n j + ρ n j+1 + λ 2 f (ρ n j-1 ) -f (ρ n j+1 ) . (36)
In fact, when η = ∆x, N = 1 in (7) and the scheme [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF], [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF] coincides with (36).

Highly oscillating initial datum

We consider the following initial datum ρ 0 (x) = 0.5 (1 + sin(10πx)) if x ∈ ]-0.5, 0.5[ , 0.5 otherwise, (37) see Figure 5. We are interested in exploring the effect of the location of the support of w η (x) = 1/η on the total variation of the solution. Figure 6 shows the density profiles at t = 0.5 (left) and the total variation TV(ρ(t, •); [-1, 1]) for t ∈ [0, 0.5] (right). We observe that the downstream convolution (2) smooths down the solution towards the constant profile ρ(t, x) = 0.5 for t → ∞. The central convolution (29) still reduces oscillations, but the upstream convolution induces blow up both in the L ∞ and in the BV norms.

  with a ∧ b = max(a, b) and a ∨ b = min(a, b).

Fig. 1

 1 Fig. 1 Density profiles at t = 0.5 corresponding to a Riemann-like initial datum with ρ L = 0.6, ρ R = 0.2. Case (a) corresponds to the rarefaction profile that solves the classical equation (32). Cases (b), (c) and (d) correspond to the kernel wη(x) = 1/η with downstream (2), central (29) and upstream (30) supports respectively.

Fig. 2

 2 Fig. 2 Density profiles at t = 0.2 corresponding to a Riemann-like initial datum with ρ L = 0.4, ρ R = 0.9. Case (a) corresponds to the shock profile that solves the classical equation (32). Cases (b), (c) and (d) correspond to the kernel wη(x) = 1/η with downstream (2), central (29) and upstream (30) supports respectively.

Fig. 3

 3 Fig. 3 Density profiles at t = 0.5 (left) and the total variation TV(ρ(t, •); [-1, 1]) for t ∈ [0, 1] (right) corresponding to the Riemann-like initial datum with ρ L = 0.4, ρ R = 0.9 for wη(x) = 1/η (top), wη(x) = 2(ηx)/η 2 (middle) and wη(x) = 2x/η 2 (bottom).

Fig. 4

 4 Fig. 4 Density profiles at t = 0.5 (top) and and t = 1 (bottom) corresponding to the Riemann-like initial datum with ρ L = 0.4, ρ R = 0.9 for wη(x) = 1/η and η = 0.1, 0.01, 0.001.

Fig. 5 Fig. 6

 56 Fig.5Oscillating initial datum ρ 0 given by (37)

Acknowledgements

The authors wish to thank Rinaldo M. Colombo and Frédéric Lagoutière for valuable discussions.

This research was partially supported by the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) / ERC Grant Agreement n. 257661.