Aekta Aggarwal

Rinaldo M Colombo

AND Paola Goatin

NONLOCAL SYSTEMS OF CONSERVATION LAWS IN SEVERAL SPACE DIMENSIONS

Keywords: Nonlocal Conservation Laws, Lax-Friedrichs Scheme AMS subject classifications. 35L65, 65M12

We present a Lax-Friedrichs type algorithm to numerically integrate a class of nonlocal and nonlinear systems of conservation laws in several space dimensions. The convergence of the approximate solutions is proved, also providing the existence of solution in a slightly more general setting than in other results in the current literature. An application to a crowd dynamics model is considered. This numerical algorithm is then used to test the conjecture that as the convolution kernels converge to a Dirac δ, the nonlocal problem converges to its non-nonlocal analogue.

1. Introduction. We are concerned with the following class of nonlocal systems of N conservation laws in d space dimensions:

∂ t U + div x F (t, x, U, η * U ) = 0 (1.1)
with a nonlocal flow of the type

F : R + × R n × R N × R m → R N ×d t x U η * U → F (t, x, U, η * U )
and where the coupling among the equations is due only to the nonlocal terms η * U , given by, for ℓ = 1, . . . , m,

(η * U ) ℓ (t, x) = R n N k=1
η ℓ,k (x -ξ) U k (t, ξ) dξ , and η : R n → R m×N .

Below we consider a Lax-Friedrichs type algorithm that yields a sequence of approximate solutions to (1.1) that, up to a subsequence, converges to a weak entropy solution of (1.1). As usual, t is time, x the space variable, U the vector of the unknown densities, F the matrix valued flow and η is a matrix of smooth convolution kernels. Systems of the type (1.1), with the particular coupling considered below, comprise, for instance, sedimentation models [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF] and various crowd dynamics models [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]. The nonlocal nature of (1.1) is particularly suitable in describing the behavior of crowds, where each member moves according to her/his evaluation of the crowd density and its variations within her/his horizon. We refer to the several models recently considered in the literature, e.g. [START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]. In one space dimension, vehicular traffic offers entirely similar situations, see [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF]. In this connection, we recall that nonlocal systems of conservation laws are widely used in biological applications, e.g., in structured population dynamics, see [START_REF] Perthame | Transport equations in biology[END_REF]. The algorithm presented below can easily be adapted to the typical case of the initial boundary value problem for a renewal equation, once the flow at the boundary is pointing into the domain and the source terms are dealt through a standard operator splitting technique, see [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Chapter 17]. Further applications of nonlocal conservation laws are found in [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] for supply chain models, in [START_REF] Amadori | An integro-differential conservation law arising in a model of granular flow[END_REF] for granular material dynamics and the case of a conveyor belt is considered in [START_REF] Göttlich | Modeling, simulation and validation of material flow on conveyor belts[END_REF].

The key point in the convergence proof lies in obtaining a uniform bound on the total variation of the approximate solutions. However, a bound on the L ∞ norm of the solution is provided below, as well as the L 1 Lipschitz continuous dependence of the solution on time. As a consequence, the present result also extends to a wider class of equations the existence part of the results in [START_REF] Amorim | On the numerical integration of scalar nonlocal conservation laws[END_REF][START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF][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]. More precisely, in the scalar 1D case, a first convergence proof for a Lax-Friedrichs type algorithm was obtained in [START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF], see also [START_REF] Amorim | On the numerical integration of scalar nonlocal conservation laws[END_REF][START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF]. Analytical well posedness results for nonlocal conservation laws in the scalar case were obtained in [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] and in [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF] in the case of systems. In these references, as in the present work, the actual number of space dimensions is not relevant, if not for the computational cost.

We remark that in the construction below, no divergence free assumption is required, see [START_REF] Chainais-Hillairet | Finite volume schemes for a nonlinear hyperbolic equation. Convergence towards the entropy solution and error estimate[END_REF] and the convolution kernel η is not required to have compact support.

From an analytical point of view, the limit where η converges to a Dirac measure δ is particularly appealing. Indeed, in that limit system (1.1) tends to a general system of conservation laws in several space dimensions and (at the time of this writing) its well posedness is a well known open problem. Nevertheless, using the algorithm below, in Section 3 we numerically investigate this limit, although the present estimates clearly do not allow to rigorously justify any conclusion in this direction.

The next section presents the main results, while Section 3 is devote to sample numerical integrations of (1.1). The technical details are deferred to Section 4.

Main Results. Throughout, we denote

R + = [0, +∞[ and R N + = [0, +∞[ N .
Aiming at a formal simplification, we consider throughout the case of d = 2 space dimensions and denote the space variables by (x, y) ∈ R 2 . Moreover, we fix our attention to the k-th equation in (1.1), for k =, . . . , N , namely

∂ t U k + ∂ x f k (t, x, y, U k , η * U ) + ∂ y g k (t, x, y, U k , ϑ * U ) = 0 (2.1)
where

U (t, x, y) = (U 1 , . . . , U N )(t, x, y) ∈ R N , U k (t, x, y) ∈ R , η(x, y) =     η 11 . . . η 1N . . . • • • . . . η m1 . . . η mN     (x, y) ∈ R m×N , η ℓ,k (x, y) ∈ R , ϑ(x, y) =     ϑ 11 . . . ϑ 1N . . . • • • . . . ϑ m1 . . . ϑ mN     (x, y) ∈ R m×N , ϑ ℓ,k (x, y) ∈ R , and (η * U )(t, x, y) ∈ R m , (ϑ * U )(t, x, y) ∈ R m , where for ℓ = 1, . . . , m , (η * U ) ℓ (t, x, y) = R 2 N k=1 η ℓ,k (x -x ′ , y -y ′ ) U k (t, x ′ , y ′ ) dx ′ dy ′ ∈ R , (ϑ * U ) ℓ (t, x, y) = R 2 N k=1 ϑ ℓ,k (x -x ′ , y -y ′ ) U k (t, x ′ , y ′ ) dx ′ dy ′ ∈ R .
Throughout, the notation |ξ| is used for ξ ∈ R, whereas for ξ ∈ R ν , ξ stands for the Euclidean norm. The following assumptions on the maps f k = f k (t, x, y, u, A), g k = g k (t, x, y, u, B) and η, ϑ in (2.1), for k = 1, . . . , N are used below:

(H0) f k , g k ∈ C 2 (R + × R 2 × R× R m ; R); ∂ u f k , ∂ u g k ∈ L ∞ (R + × R 2 × R× R m ; R); for all t ∈ R + , (x, y) ∈ R 2 and A, B ∈ R m , f k (t, x, y, 0, A) = g k (t, x, y, 0, B) = 0. (H1) There exists a constant M > 0 such that max ∂ x f k (t, x, y, u, A) , ∂ y g k (t, x, y, u, B) ≤ M |u|, max ∇ A f k (t, x, y, u, A) , ∇ B g k (t, x, y, u, B) ≤ M |u|, max ∂ 2 xx f k (t, x, y, u, A) , ∂ 2 yy g k (t, x, y, u, B) ≤ M |u|, max ∂ 2 xy f k (t, x, y, u, A) , ∂ 2 xy g k (t, x, y, u, B) ≤ M |u|, max ∂ x ∇ A f k (t, x, y, u, A) , ∂ x ∇ B g k (t, x, y, u, B) ≤ M |u|, max ∂ y ∇ A f k (t, x, y, u, A) , ∂ y ∇ B g k (t, x, y, u, B) ≤ M |u|, max ∇ 2 AA f k (t, x, y, u, A) , ∇ 2 BB g k (t, x, y, u, B) ≤ M |u|. (H2) ∂ u f k , ∂ u g k ∈ W 1,∞ (R + × R 2 × R × R m ; R). (H3) η, ϑ ∈ (C 2 ∩ W 2,∞ )(R 2 ; R m×N ). Remark 1.
Standard situations of interest lead to compactly supported initial data yielding solutions that are also compactly supported at all finite times. In these cases, it is sufficient that the above uniform boundedness requirements are satisfied on any compact subset of R 2 . Equivalently, on f k and g k it is enough to assume only conditions (H0) and (H3). The proof is immediate and follows the same lines, for instance, of that in [START_REF] Borsche | Differential equations modeling crowd interactions[END_REF]Corollary 2.3].

Recall the definition [START_REF] Amorim | On the numerical integration of scalar nonlocal conservation laws[END_REF][START_REF] Betancourt | On nonlocal conservation laws modelling sedimentation[END_REF][START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][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] of solution to the Cauchy problem for (2.1). Definition 2.1.

Let U o ∈ L ∞ (R 2 ; R N ). A map U : [0, T ] → L ∞ (R 2 ; R N
) is a solution to (2.1) with initial datum U o if, for k = 1, . . . , N , setting for all w ∈ R f k (t, x, y, w) = f k t, x, y, w, (η * U )(t, x, y) , gk (t, x, y, w) = g k t, x, y, w, (ϑ * U )(t, x, y) , the map U k is a Kružkov solution to the conservation law

∂ t U k + ∂ x f k (t, x, y, U k ) + ∂ y gk (t, x, y, U k ) = 0 U k (0, x, y) = U k o (x, y) . (2.2)
Above, by Kružkov solution we refer to [17, Definition 1]. It is immediate to observe that (H0) ensures that U (t, x, y) ≡ 0 solves (2.1) in the sense of Definition 2.1. Fix a rectangular grid with sizes ∆x and ∆y in R 2 and choose a time step ∆t. For later use, we also introduce the usual notation

t n = n ∆t, n ∈ Z ; x i = i ∆x, x i+1/2 = i + 1 2 ∆x, i ∈ Z ; y j = j ∆y, y j+1/2 = j + 1 2 ∆y, j ∈ Z ; λ x = ∆t/∆x , λ y = ∆t/∆y . Throughout, an initial datum U o ∈ (L ∞ ∩ BV)(R 2 ; R N
) is fixed and we denote

u k,0 ij = 1 ∆x ∆y x i+1/2 x i-1/2 y i+1/2 y i-1/2 U k o (x, y) dx dy for i, j ∈ Z .
We define a piecewise constant approximate solution u ∆ ≡ u 1 ∆ , . . . , u N ∆ to (2.1) by

u k ∆ (t, x, y) = u k,n ij for      t ∈ [t n , t n+1 [ , x ∈ [x i-1/2 , x i+1/2 [ , y ∈ [y j-1/2 , y j+1/2 [ , where n ∈ N i ∈ Z j ∈ Z k ∈ {1, . . . , n}
through the 5-points algorithm based on dimensional splitting, see [START_REF] Crandall | The method of fractional steps for conservation laws[END_REF]Section 3]:

u k,n+1/2 ij = u k,n ij -λ x F k,n i+1/2,j (u k,n ij , u k,n i+1,j ) - F k,n i-1/2,j (u k,n i-1,j , u k,n ij ) u k,n ij = u k,n+1/2 ij -λ y G k,n i,j+1/2 (u k,n+1/2 ij , u k,n+1/2 i,j+1 ) -G k,n i,j-1/2 (u k,n+1/2 i,j-1 , u k,n+1/2 ij ) (2.
3) where for fixed α, β in 0, 2/3 , we define

F k,n i+1/2,j (u, v) = f k,n i+1/2,j (u) + f k,n i+1/2,j (v) 2 - α(v -u) 2 λ x , G k,n i,j+1/2 (u, v) = g k,n i,j+1/2 (u) + g k,n i,j+1/2 (v) 2 - β(v -u) 2 λ y , f k,n i+1/2,j (u) = f k (t n , x i+1/2 , y j , u, A n i+1/2,j ) , g k,n i,j+1/2 (u) = g k (t n , x i , y j+1/2 , u, B n i,j+1/2 ) .
The convolution terms are computed through quadrature formulae, i.e.,

A n i+1/2,j = ∆x ∆y     l,p∈Z N k=1 η 1,k i+1/2-l,j-p u k,n l+1/2,p , . . . , l,p∈Z N k=1 η m,k i+1/2-l,j-p u k,n l+1/2,p     , B n i,j+1/2 = ∆x ∆y     l,p∈Z N k=1 ϑ 1,k i+1/2-l,j-p u k,n l+1/2,p , . . . , l,p∈Z N k=1 ϑ m,k i+1/2-l,j-p u k,n l+1/2,p     , (2.4
) where, for instance, u k,n l+1/2,p is a convex combination of u k,n l,p and u k,n l+1,p , moreover

η ℓ,k i+1/2,j = η ℓ,k (x i+1/2 , y j ) and ϑ ℓ,k i+1/2,j = ϑ ℓ,k (x i+1/2 , y j ) . (2.5)
Throughout, we require that ∆t is chosen in order to satisfy the CFL condition

λ x ≤ 2 min{3α, 2 -3α} 1 + 6 max k ∂ u f k L ∞ , λ y ≤ 2 min{3β, 2 -3β} 1 + 6 max k ∂ u g k L ∞ , (2.6)
and we assume that the mesh sizes are sufficiently small, in the sense that ∆x, ∆y ≤ 1/(3M ), where M is as in (H0).

The main result of this paper is the following: Theorem 2.2. Let assumptions (H0), (H1), (H3) and conditions (2.6) hold. Fix an initial datum

U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, the algorithm (2.3) defines a sequence of approximate solutions which converges, up to a subsequence, to

a solution U ∈ C 0 R + ; L 1 (R 2 ; R N + ) to (2.1
) in the sense of Definition 2.1. Moreover, U is locally Lipschitz continuous with respect to the L 1 norm and for all k ∈ {1, . . . , N } and for all t ∈ R + , the following bounds hold:

U (t) L ∞ (R 2 ;R N ) ≤ e Ct U o L ∞ (R 2 ;R N ) with C as in (4.7), U k (t) L 1 (R 2 ;R) = U k o L 1 (R 2 ;R) , TV(U k (t)) ≤ e K1 t TV(U k o ) + K 2 e K1 t -1 with K 1 , K 2 as in (4.19) , U (t + τ ) -U (t) L 1 (R 2 ;R N ) ≤ C(t) τ with C(t) as in (4.22) .
The proof is deferred to Section 4 and is based on the lemmas below. Observe that Theorem 2.2 also ensures the existence of solutions to (2.1) with a flow more general than that considered in [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]. The solution is proved to be a Lipschitz continuous function of time with respect to the L 1 topology. Moreover, in view of Remark 1, when the initial datum is compactly supported, (H0) and (H3) are sufficient to guarantee the existence of a solution and the bounds above on any finite time interval. As first step, we prove the positivity of the approximate solutions defined by (2.3). Lemma 2.3 (Positivity). Let (H0), (H1), (H3) and

(2.6) hold. Fix a U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, the approximate solution u ∆ defined by the algorithm

(2.3) satisfies, for all k ∈ {1, . . . , N }, t ∈ R + and (x, y) ∈ R 2 . u k ∆ (t, x, y) ≥ 0 .
A consequence of positivity, due to the conservative nature of the algorithm, is the L 1 bound.

Lemma 2.4 (L 1 bound). Let (H0), (H1), (H3) and (2.6) hold. Fix an initial datum

U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, the approximate solution u ∆ defined by the algorithm (2.3) satisfies, for all k ∈ {1, . . . , N } and t ∈ R + .

u k ∆ (t) L 1 = u k ∆ (0) L 1 . The proof is omitted, since it is standard, see for instance [3, Lemma 2.4]. Lemma 2.4 allows to prove the L ∞ bound. Lemma 2.5 (L ∞ bound). Let (H0), (H1), (H3) and (2.6) hold. Fix an initial datum U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, there exists a C depending only on η, ϑ, f 1 , . . . , f N , g 1 , . . . , g N , see (4.7), such that the approximate solution u ∆ defined by the algorithm (2.3) satisfies for all t ∈ R + ,

u ∆ (t) L ∞ ≤ e C t(1+ Uo L 1 ) U o L ∞ .
(2.7)

The previous results allow to get the following key bound on the total variation. Lemma 2.6 (BV bound). Let (H0), (H1), (H2), (H3) and (2.6) hold. Fix an initial datum

U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, there exist constants K 1 and K 2 depending only on U o L 1 , η, ϑ, f 1 , . . . , f N , g 1 , . . . , g N , , see (4.19), such that the approximate solution u ∆ defined by the algorithm (2.3) satisfies for all n, for all t ∈ t n , t n+1 and for all k = 1, . . . ,

N ij u k,n i+1,j -u k,n ij ∆y + u k,n i,j+1 -u k,n ij ∆x ≤ e K1 t ij u k,0 i+1,j -u k,0 ij ∆y + u k,0 i,j+1 -u k,0 ij ∆x + K 2 e K1 t -1 .
Remark that K 2 depends on all the components of U o . Hence an initial datum with a constant component may well lead to a solution with no constant components.

To prove the discrete entropy inequality, we introduce the following Kružkov numerical entropy flows, see [

Φ k,n,κ i+1/2,j (u 1 , u 2 ) = f k,n i+1/2,j (u 1 ∨ κ, u 2 ∨ κ) -f k,n i+1/2,j (u 1 ∧ κ, u 2 ∧ κ) , Γ k,n,κ i,j+1/2 (u 1 , u 2 ) = g k,n i,j+1/2 (u 1 ∨ κ, u 2 ∨ κ) -g k,n i,j+1/2 (u 1 ∧ κ, u 2 ∧ κ) .
Lemma 2.7 (Discrete entropy condition). Let assumptions (H0), (H1), (H2), (H3) and conditions (2.6) hold. Fix an initial datum

U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, the approximate solution u ∆ defined by the algorithm (2.3) satisfies for all i, j ∈ Z and for all k = 1, . . . , N and for all κ ∈ R the discrete entropy inequality

u k,n+1 ij -κ -u k,n ij -κ + λ x Φ k,n,κ i+1/2,j (u k,n ij , u k,n i+1,j ) -Φ k,n,κ i-1/2,j (u k,n i-1,j , u k,n ij ) + λ x sgn(u k,n+1/2 ij -κ) f k,n i+1/2,j (κ) -f k,n i-1/2,j (κ) + λ y Γ k,n,κ i,j+1/2 (u k,n+1/2 ij , u k,n+1/2 i,j+1 ) -Γ k,n,κ i,j-1/2 (u k,n+1/2 i,j-1 , u k,n+1/2 ij ) + λ y sgn(u k,n+1 ij -κ) g k,n i,j+1/2 (κ) -g k,n i,j-1/2 (κ) ≤ 0 .
The proof is entirely analogous to that of [3, Proposition 2.8]. Lemma 2.8 (L 1 Lipschitz Continuity in time). Let assumptions (H0), (H1), (H2), (H3) and conditions (2.6) hold. Fix an initial datum

U o ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R N + ).
Then, for any n ∈ N, there exists a constant C depending on U o L 1 , TV(U o ), t n , λ x , λ y α, β and on the functions η, ϑ, f k and g k , for k ∈ {1, . . . , N } such that for n = 1, . . . , n the approximate solution u ∆ defined by (2.3)

satisfies u ∆ (t n+1 ) -u ∆ (t n ) L 1 ≤ C ∆t .
The proof is a straightforward generalization of [3, Lemma 2.7].

3. Numerical Integrations.

3.1.

A Crowd Dynamics Sample Integration. We now show that the crowd dynamics model introduced in [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] fits in the framework of (2.1). Assume that a group of pedestrians can be described through the density U = U (t, x, y) that satisfies the nonlocal continuity equation

∂ t U + ∇ • U (1 -U ) (1 -U * µ) v = 0 . (3.1)
Here, the given vector field v = v(x, y) ≡ v 1 (x, y), v 2 (x, y) describes the path followed by the pedestrians. The smooth, non negative and compactly supported function µ models the way in which each individual averages the density around her/his position to adjust her/his speed. Lemma 3.1. System (3.1) fits into (2.1) setting

N = 1 , m = 1 , f (t, x, y, U, A) = U (1 -U ) (1 -A) v 1 (x, y) g(t, x, y, U, B) = U (1 -U ) (1 -B) v 2 (x, y) and η = µ ϑ = µ . (3.2) Moreover, if v ∈ (C 2 ∩W 2,∞ )(R 2 , R 2 ) and µ ∈ (C 2 ∩W 2,∞ )(R 2 ; R), then Theorem 2.2 applies to any initial datum in (L 1 ∩ L ∞ ∩ BV)(R 2 ; [0, 1]).
The proof in Section 4 essentially relies on the invariance of the interval [0, 1] for u.

We now consider a specific situation. Two groups of people move from different places along a corridor [0, 10] × [-1, 1] towards one of its ends and stop in front of {9.5} × [-1, 1]. Therefore, in (3.1) we choose: and with λ x = λ y = 0.15. The convergence rate γ is then calculated at time t = 1.0 by computing the L 1 distance between the numerical solutions u h and u h 2 obtained for the grid size h and h 2 , for each grid size h in (3.5). The results recorded in Table 3.1 show that the numerical scheme (2.3) has a convergence rate γ ≈ 0.5..

v(x, y) = (1 -y 2 ) 3 exp -1/(x -9.5) 2 χ ]-∞,9.5]×[-1,1] (x, y) -2 y exp(1 -1/y 2 ) ρ(y) , (3.3
h u h -u h 2 L 1 γ = log 2 u h -u h 2 L 1 u h 2 -u h 4

The Nonlocal Keyfitz-Kranzer

System. We consider the example of the classical Keyfitz-Kranzer system [START_REF] Keyfitz | A system of nonstrictly hyperbolic conservation laws arising in elasticity theory[END_REF] which is given by

∂ t U 1 + ∂ x U 1 ϕ(U ) = 0 ∂ t U 2 + ∂ x U 2 ϕ(U ) = 0 (3.6)
is a 2 × 2 system of conservation laws in one space dimension in the unknown U = (U 1 , U 2 ). System (3.6) has been extended and studied from many different points of view. For instance, the existence of solutions in the case of several space dimensions was obtained in [2, Theorem 2.6]. The convergence of approximate solutions was proved in [START_REF] Koley | Finite difference schemes for the symmetric keyfitzkranzer system[END_REF] in the case of one dimensional symmetric systems. At the time of this writing, an open problem in the analytical theory of (1.1) is whether the formal convergence to its non-nonlocal analogue ∂ t U +div x F (t, x, U, U ) = 0 as η tends to Dirac delta can be made rigorous. The present numerical procedures can be helpful in having some insight in this issue. To this end, we consider the classical Keyfitz-Kranzer system in 2 dimensions

∂ t U 1 + ∂ x U 1 ϕ 1 (U ) + ∂ y U 1 ϕ 2 (U ) = 0 ∂ t U 2 + ∂ x U 2 ϕ 1 (U ) + ∂ y U 2 ϕ 2 (U ) = 0 (3.7)
with its non-local generalization given by

∂ t U 1 + ∂ x U 1 ϕ 1 (U * µ) + ∂ y U 1 ϕ 2 (U * µ) = 0 ∂ t U 2 + ∂ x U 2 ϕ 1 (U * µ) + ∂ y U 2 ϕ 2 (U * µ) = 0 (3.8)
which is a system of 2 nonlocal conservation laws in 2 space dimensions. This system fits into the framework of Theorem 2.2. Lemma 3.2. System (3.8) fits into (2.1) setting:

N = 2 , m = 2 , f (t, x, y, U, A) = U ϕ(A) g(t, x, y, U, B) = U ϕ(B) and η = ϑ = µ 0 0 µ . Moreover, if ϕ 1 , ϕ 2 ∈ (C 2 ∩ W 2,∞ )(R 2 ; R) and µ ∈ (C 2 ∩ W 2,∞ )(R 2 ; R), then for any compactly supported initial datum (U 1 o , U 2 o ) ∈ (L 1 ∩ L ∞ ∩ BV)(R 2 ; R 2 + )
, problem (3.8) admits a solution defined for all t ∈ R + and satisfying the estimates in Theorem 2.2. The proof of this Lemma, i.e., that (H0) and (H3) hold in the case of (3.8), is immediate and, hence, omitted.

The following numerical integration of (3. so that R 2 µ(x, y) dx dy = 1. As initial datum, we choose:

ϕ 1 (A 1 , A 2 ) = sin(A 1 2 + A 2 2 ) , ϕ 2 (B 1 , B 2 ) = cos(B 1 2 + B 2 2 ) (3.
U o (x, y) =              (1, √ 3) (x, y) ∈ ]0, 0.4] × ]0, 0.4] ( √ 2, 1) (x, y) ∈ ]-0.4, 0] × ]0, 0.4] ( 1 2 , 1 3 ) (x, y) ∈ ]-0.4, 0] × ]-0.4, 0] ( √ 3, √ 2
u n = (1 + H) n U o + (1 + H) n -1 H K ≤ (1 + H) n (U o + n K) .
Proof. The first equality is immediate. Then, using the binomial theorem, we get

(1 + H) n -1 H = n h=1 n h H h-1 ≤ n-1 h=0 n n -1 h H h ≤ n(1 + H) n .
Lemma 4.2. Under conditions (H3), the algorithm (2.3) yields the estimates

A n i+1/2,j -A n i-1/2,j ≤ ∆x ∂ x η L ∞ u ∆ (t n ) L 1 , (4.1)
A n i+1/2,j+1 -A n i+1/2,j ≤ ∆y ∂ y η L ∞ u ∆ (t n ) L 1 , B n i,j+1/2 -B n i,j-1/2 ≤ ∆y ∂ y ϑ L ∞ u ∆ (t n ) L 1 , B n i+1,j+1/2 -B n i,j+1/2 ≤ ∆x ∂ x ϑ L ∞ u ∆ (t n ) L 1 , A n i+3/2,j -2A n i+1/2,j + A n i-1/2,j ≤ 2 ∆x 2 ∂ 2 xx η L ∞ u ∆ (t n ) L 1 ,(4.2) B n i,j+3/2 -2B n i,j+1/2 + B n i,j-1/2 ≤ 2 ∆y 2 ∂ 2 yy ϑ L ∞ u ∆ (t n ) L 1 , A n i+1/2,j -A n i+1/2,j+1 -A n i-1/2,j +A n i-1/2,j+1 ≤ ∆x ∆y ∂ 2 xy η L ∞ u ∆ (t n ) L 1 , B n i,j+1/2 -B n i+1,j+1/2 -B n i,j-1/2 +B n i+1,j+1/2 ≤ ∆x ∆y ∂ 2 xy ϑ L ∞ u ∆ (t n ) L 1 .
Proof. Consider the estimates on the A terms, the others being entirely similar:

A n i+1/2,j -A n i-1/2,j ≤ ∆x ∆y l,m∈Z (η i+1/2-l,j-m -η i-1/2-l,j-m ) u ∆ (t n , x l+1/2 , y m ) ≤ ∆x ∆y l,m∈Z η i+1/2-l,j-m -η i-1/2-l,j-m u ∆ (t n , x l+1/2 , y m ) ≤ ∆x ∆y l,m∈Z u ∆ (t n , x l+1/2 , y m ) x i+1/2-l x i-1/2-l ∂ x η(x, y j-m ) dx ≤ ∆x ∂ x η k L ∞ ∆x ∆y l,m∈Z u ∆ (t n , x l+1/2 , y m ) ≤ ∆x ∂ x η k L ∞ u ∆ (t n ) L 1
, proving the first estimate, the second one is entirely analogous. Similarly,

A n i+3/2,j -2A n i+1/2,j + A n i-1/2,j ≤ ∆x ∆y l,m∈Z u ∆ (t n , x l+1/2 , y m ) η i+3/2-l,j-m -2η i+1/2-l,j-m + η i-1/2-l,j-m ≤ ∆x ∆y l,m∈Z u ∆ (t n , x l+1/2 , y m ) ∆x ∂ x η(τ i+1-l , y j-m ) -∂ x η(τ i-l , y j-m ) ≤ ∆x ∆y l,m∈Z u ∆ (t n , x l+1/2 , y m ) ∆x τ i+1-l τ i-l ∂ 2 xx η(s, y j-m ) ds ≤ 2∆x 2 ∂ 2 xx η L ∞ u ∆ (t n ) L 1 where τ i-l ∈ x i-1/2-l , x i+1/2-l .
Proof of Lemma 2.3. In the spirit of [3, Lemma 2.5], note that (2.3) can be rewritten as follows:

u k,n+1/2 ij = u k,n ij -a k,n i-1/2,j (u k,n ij -u k,n i-1,j ) + b k,n i+1/2,j (u k,n i+1,j -u k,n ij ) -λ x F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i-1/2,j (u k,n ij , u k,n ij ) , (4.3)
where

a k,n i-1/2,j = λ x F k,n i-1/2,j (u k,n ij , u k,n ij ) -F k,n i-1/2,j (u k,n i-1,j , u k,n ij ) u k,n ij -u k,n i-1,j , b k,n i+1/2,j = λ x F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i+1/2,j (u k,n ij , u k,n i+1,j ) u k,n i+1,j -u k,n ij .
We now prove that a k,n i-1/2,j and b k,n i+1/2,j are non negative. Indeed,

a k,n i-1/2,j = λ x F k,n i-1/2,j (u k,n ij , u k,n ij ) -F k,n i-1/2,j (u k,n i-1,j , u k,n ij ) u k,n ij -u k,n i-1,j = λ x u k,n ij -u k,n i-1,j   f k,n i-1/2,j (u k,n ij ) - f k,n i-1/2,j (u k,n i-1,j ) + f k,n i-1/2,j (u k,n ij ) 2   + α 2 = λ x u k,n ij -u k,n i-1,j   f k,n i-1/2,j (u k,n ij ) -f k,n i-1/2,j (u k,n i-1,j ) 2   + α 2 = λ x ∂ u f k,n i-1/2,j (ζ k,n i-1/2,j ) + α 2 for a suitable ζ k,n i-1/2,j lying between u k,n i-1,j and u k,n ij . Hence a k,n i-1/2,j ≥ -λ x ∂ u f k L ∞ + α 2 ≥ 0 , a k,n i-1/2,j ≤ λ x ∂ u f k L ∞ + α 2 ≤ 1 3 (4.4) by (2.6)
. Entirely similar computations hold for the term b k,n i-1/2,j , so that 1-a k,n i-1/2,jb k,n i-1/2,j , a k,n i-1/2,j and b k,n i-1/2,j are the coefficients of a convex combination. Using (H0), (H1) and the CFL condition (2.6) we now estimate:

F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i-1/2,j (u k,n ij , u k,n ij ) = f k (t n , x i+1/2 , y j , u k,n ij , A n i+1/2,j ) -f k (t n , x i+1/2 , y j , u k,n ij , A n i-1/2,j ) +f k (t n , x i+1/2 , y j , u k,n ij , A n i-1/2,j ) -f k (t n , x i-1/2 , y j , u k,n ij , A n i-1/2,j ) ≤ f k (t n , x i+1/2 , y j , u k,n ij , A n i+1/2,j ) + f k (t n , x i+1/2 , y j , u k,n ij , A n i-1/2,j ) f k (t n , x i+1/2 , y j , u k,n ij , A n i-1/2,j ) -f k (t n , x i-1/2 , y j , u k,n ij , A n i-1/2,j ) ≤ 2 ∂ u f k L ∞ u k,n ij + ∂ x f k (t n , xi , y j , u k,n ij , A n i-1/2,j ) ∆x ≤ 2 ∂ u f k L ∞ + M ∆x u k,n ij
for some xi ∈ x i-1/2 , x i+1/2 . Now using condition (2.6) we have:

F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i-1/2,j (u k,n ij , u k,n ij ) ≤ 1 3λ x u k,n ij , (4.5)
Using (4.4) and (4.5) in (4.3) we get:

u k,n+1/2 ij = (1 -a k,n i-1/2,j -b k,n i+1/2,j )u k,n ij + a k,n i-1/2,j u k,n i-1,j + b k,n i+1/2,j u k,n i+1,j
Proof of Lemma 2.6. We follow the proof of [

-u k,n+1/2 ij = C n ij -λ x D n ij ,
where

C n ij = (u k,n i+1,j -u k,n ij ) -λ x F k,n i+3/2,j (u k,n i+1,j , u k,n i+2,j ) -F k,n i+1/2,j (u k,n ij , u k,n i+1,j ) -λ x -F k,n i+3/2,j (u k,n ij , u k,n i+1,j ) + F k,n i+1/2,j (u k,n i-1,j , u k,n ij ) D n ij = F k,n i+3/2,j (u k,n ij , u k,n i+1,j ) -F k,n i+1/2,j (u k,n ij , u k,n i+1,j ) +F k,n i-1/2,j (u k,n i-1,j , u k,n ij ) -F k,n i+1/2,j (u k,n i-1,j , u k,n ij )
and we now estimate the various terms above separately. Consider the term C n ij .

C n ij = (u k,n i+1,j -u k,n ij )   1 + λ x F k,n i+1/2,j (u k,n ij , u k,n i+1,j ) -F k,n i+1/2,j (u k,n ij , u k,n ij ) u k,n i+1,j -u k,n ij   +(u k,n i+1,j -u k,n ij )   -λ x F k,n i+3/2,j (u k,n i+1,j , u k,n i+1,j ) -F k,n i+3/2,j (u k,n ij , u k,n i+1,j ) u k,n i+1,j -u k,n ij   +(u k,n i+2,j -u k,n i+1,j )   -λ x F k,n i+3/2,j (u k,n i+1,j , u k,n i+2,j ) -F k,n i+3/2,j (u k,n i+1,j , u k,n i+1,j ) u k,n i+2,j -u k,n i+1,j   +(u k,n ij -u k,n i-1,j )   λ x F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i+1/2,j (u k,n i-1,j , u k,n ij ) u k,n ij -u k,n i-1,j   = (u k,n i+1,j -u k,n ij )(1 -b k,n i+ 1 2 ,j -ãk,n i+ 1 2 ,j ) + (u k,n i+2,j -u k,n i+1,j )b k,n i+ 3 2 ,j +(u k,n ij -u k,n i-1,j )ã k,n i-1 2 ,j
where ãk,n

i-1 2 ,j = λ x F k,n i+1/2,j (u k,n ij , u k,n ij ) -F k,n i+1/2,j (u k,n i-1,j , u k,n ij ) (u k,n ij -u k,n i-1,j ) and ãk,n i-1
2 ,j ∈ [0, 1/3] can be proved exactly as was done for a k,n i-1 2 ,j . By convexity, we have

ij C n ij ≤ ij u k,n i+1,j -u k,n ij . (4.8) Concerning D n ij , D n ij = 1 2 f k,n i+3/2,j (u k,n ij ) -2f k,n i+1/2,j (u k,n ij ) + f k,n i-1/2,j (u k,n ij ) (4.9) + 1 2 f k,n i+3/2,j (u k,n i+1,j ) -f k,n i+1/2,j (u k,n i+1,j ) (4.10) - 1 2 f k,n i+1/2,j (u k,n i-1,j ) -f k,n i-1/2,j (u k,n i-1,j ) . ( 4
= 1 2 f k (t n , x i+3/2 , y j , u k,n i+1,j , A n i+3/2,j ) -f k (t n , x i+1/2 , y j , u k,n i+1,j , A n i+3/2,j ) + 1 2 f k (t n , x i+1/2 , y j , u k,n i+1,j , A n i+3/2,j ) -f k (t n , x i+1/2 , y j , u k,n i+1,j , A n i+1/2,j ) - 1 2 f k (t n , x i+1/2 , y j , u k,n i-1,j , A n i+1/2,j ) -f k (t n , x i-1/2 , y j , u k,n i-1,j , A n i+1/2,j ) - 1 2 f k (t n , x i-1/2 , y j , u k,n i-1,j , A n i+1/2,j ) -f k (t n , x i-1/2 , y j , u k,n i-1,j , A n i-1/2,j ) = 1 2 ∆x ∂ x f k (t n , x i+1 , y j , u k,n i+1,j , A n i+3/2,j ) -∂ x f k (t n , x i , y j , u k,n i-1,j , A n i+1/2,j ) + 1 2 ∇ A f k (t n , x i+1/2 , y j , u k,n i+1,j , A n i+1,j )(A n i+3/2,j -A n i+1/2,j ) - 1 2 ∇ A f k (t n , x i-1/2 , y j , u k,n i-1,j , A n i,j )(A n i+1/2,j -A n i-1/2,j ) for suitable x i ∈ x i-1/2 , x i+1/2 and A n i,j ∈ I(A n i-1/2,j , A n i+1/2,j ). Introduce xi+1/2 ∈ ]x i , x i+1 [, xi ∈ x i-1/2 , x i+1/2 , u k,n ij ∈ I(u k,n i-1,j , u k,n i+1,j ), ǔk,n ij ∈ I(u k,n i-1,j , u k,n i+1,j ), Ân i,j ∈ I(A n i-1/2,j , A n i+1/2,j ), Ǎn i+1/2,j ∈ I(A n i,j , A n i+1,j ) and Ãn i+1/2,j ∈ I(A n i+1/2,j - A n i-1/2,j , A n i+3/2,j -A n i+1/2,j
) so that we have (4.10) + (4.11)

≤ ∆x 2 ∆x ∂ 2 xx f k (t n , xi+1/2 , y j , u k,n i,j , Ân i+1,j ) + ∂ x ∂ u f k (t n , xi+1/2 , y j , u k,n i,j , Ân i+1,j ) u k,n i+1,j -u k,n i-1,j + ∂ x ∇ A f k (t n , xi+1/2 , y j , u k,n i,j , Ân i+1,j ) A n i+3/2,j -A n i+1/2,j + 1 2 ∆x ∂ x ∇ A f k (t n , xi , y j , ǔk,n i,j , Ǎn i+1/2,j ) Ãn i+ 1 2 ,j + ∂ u ∇ A f k (t n , xi , y j , ǔk,n i,j , Ǎn i+1/2,j ) Ãn i+ 1 2 ,j u k,n i+1,j -u k,n i-1,j + ∇ 2 AA f k (t n , xi , y j , ǔk,n i,j , Ǎn i+1/2,j ) Ãn i+ 1 2 ,j A n i+1,j -A n i,j + ∇ A f k (t n , xi , y j , ǔk,n i,j , Ǎn i+1/2,j ) A n i+3/2,j -2A n i+3/2,j + A n i-1/2,j
By Ãn

i+ 1 2 ,j ∈ I(A n i+1/2,j -A n i-1/2,j , A n i+3/2,j -A n i+1/2,j ) and A n i,j ∈ I(A n i-1/2,j , A n i+1/2,j ), using Lemma 4.2, we have Ãn i+ 1 2 ,j ≤ ∆x ∂ x η L ∞ u ∆ (t n ) L 1 and A n i+1,j -A n i,j ≤ ∆x ∂ x η L ∞ u ∆ (t n ) L 1 . Now using Lemma 2.4, we have (4.10) + (4.11) ≤ 1 2 M ∆x 2 u k,n i,j + ∆x ∂ x ∂ u f L ∞ u k,n i+1,j -u k,n i-1,j + M ∆x 2 u k,n i,j ∂ x η L ∞ U o L 1 + 1 2 M ∆x 2 ǔk,n i,j ∂ x η L ∞ U o L 1 + ∆x ∂ u ∇ A f L ∞ ∂ x η L ∞ U o L 1 u k,n i+1,j -u k,n i-1,j +M ∆x 2 ǔk,n i,j ∂ x η 2 L ∞ U o 2 L 1 + 2M ∆x 2 ǔk,n i,j ∂ 2 xx η L ∞ U o L 1 = 1 2 M ∆x 2 u k,n i,j 1 + ∂ x η L ∞ U o L 1 + 1 2 M ∆x 2 ǔk,n i,j ∂ x η 2 L ∞ U o 2 L 1 + 2 ∂ 2 xx η L ∞ U o L 1 + ∂ x η L ∞ U o L 1 + 1 2 ∆x u k,n i+1,j -u k,n i-1,j ∂ u ∇ A f L ∞ ∂ x η L ∞ U o L 1 + ∂ x ∂ u f L ∞ .
Now, let us estimate (4.9). where

1 2 f k,n i+3/2,j (u k,n ij ) -2f k,n i+1/2,j (u k,n ij ) + f k,n i-1/2,j (u k,n ij ) = 1 2 f k (t n , x i+3/2 , y j , u k,n ij , a k,n i+3/2,j ) -2f k (t n , x i+1/2 , y j , u k,n ij , a k,n i+1/2,j ) +f k (t n , x i-1/2 , y j , u k,n ij , a k,n i-1/2,j ) = 1 2 ∆x ∂ x f k (t n , xi+1 , y j , u k,n ij , Ȃn i+1,j ) -∂ x f k (t n , xi , y j , u k,n ij , Ȃn ij ) + 1 2 ∇ A f k (t n , xi+1 , y j , u k,n ij , Ȃn i+1,j )(A n i+3/2,j -A n i+1/2 ) - 1 2 ∇ A f k (t n ,
+ ∆x 2 M u k,n ij ∂ x η L ∞ U o L 1 + 1 2 ∆x 2 M u k,n ij ∂ x η L ∞ U o L 1 2 + ∆x 2 M u k,n ij ∂ 2 xx η L ∞ U o
∂ x η L ∞ U o L 1 +∆y ∂ u ∇ A f L ∞ ∂ y η L ∞ U o L 1 u k,n i,j+1 -u k,n i-1,j +M ∆x∆y ǔk,n i-1/2,j+1/2 ∂ x η L ∞ ∂ y η L ∞ U o 2 L 1 +M ∆x∆y ǔk,n i-1/2,j+1/2 ∂ 2 xy η L ∞ U o L 1 = M ∆x∆y u k,n i-1/2,j+1/2 1 + ∂ x η L ∞ U o L 1 +M ∆x∆y ǔk,n i-1/2,j+1/2 × ∂ x η L ∞ ∂ y η L ∞ U o 2 L 1 + ∂ 2 xy η L ∞ U o L 1 + ∂ x η L ∞ U o L 1 +∆y u k,n i,j+1 -u k,n i-1,j ∂ u ∇ A f L ∞ ∂ y η L ∞ U o L 1 + ∂ y ∂ u f L ∞ ,
K 2 = M 2 U o L 1 1 + ∂ x η L ∞ U o L 1 2 + ∂ y η L ∞ U o L 1 + ∂ 2 xy η L ∞ U o L 1 K 3 = 1 2 ∂ u ∇ A f L ∞ ∂ y η L ∞ U o L 1 + ∂ y ∂ u f L ∞ .
(4.17 Then, Theorem 2.2 applies and yields a solution to ∂ t U +∇• p(U ) q(U * µ) v = 0 in the sense of Definition 2.1. By [START_REF] Nikolaevich | First order quasilinear equations with several independent variables[END_REF]Theorem 3], u attains values inside [0, 1], hence it also solves (3.1).

  ) μ(x, y) = (0.16 -x 2 -y 2 ) 3 χ {(x,y) : x 2 +y 2 ≤0.16} (x, y) , µ(x, y) = 1 R 2 μ μ(x, y) , so that R 2 µ(x, y) dx dy = 1. Here, ρ ∈ C 2 c (R; R) is chosen so that ρ(y) = 1 for all y ∈ [-1, 1].The role of ρ is to allow conditions (H0), (H1) and (H2) to hold and does

Fig. 3 . 1 .

 31 Fig. 3.1. Solution to (3.1)-(3.4) at times t = 0.00, 4.02, 8.04, 12.06, 15.98, 20.00 with space mesh ∆x = ∆y = 0.0125, α = β = 0.3333 in (2.6) and ∆t ≈ 9.62E -4.

  not modify the solutions due to the invariance of the space domain [0, 10] × [-1, 1]. As initial datum we choose: U o (x, y) = χ [1,4]×[0.1,0.8] (x, y) + χ [2,5]×[-0.8,-0.1] (x, y) . (3.4) The numerical integration of (3.1)-(3.4) obtained with the algorithm described in Section 2 is shown in Figure 3.1. The two groups are seen moving to the right and accumulating in front of x = 9.5. By Definition 2.1, [17, Theorem 3] and the choice (3.4) of the initial datum, the maximal density of solution can not exceed 1, coherently with the result of the integration. The numerical domain [0, 10] × [-1, 1] is chosen invariant with respect to the flow specified in (3.3). To compute the convergence rate of the scheme (2.3), we apply the algorithm to problem (3.1)-(3.4) on the domain [0, 10] × [-1, 1] on the time interval [0, 1] with different grid sizes, namely h = ∆x = ∆y = 0.05, 0.025, 0.0125, 0.00625 . (3.5)

( 2 . 3 )

 23 on the domain [0, 10] × [-1, 1] at time t = 1.0 for the approximate solutions to the problem (3.1)-(3.4) with λx = λy = 0.15 .

  8) is referred to the domain [-1.1, 1.1]× [-1.1, 1.1] and to the time interval [0, 0.1] with

  9) and μ(x, y) = r 2 -(x 2 + y 2 ) 3 χ {(x,y) : x 2 +y 2 ≤r 2 } (x, y) , µ(x, y) = μ(x, y) R 2 μ (3.10)

Fig. 3 . 2 .Fig. 3 . 3 .

 3233 Fig. 3.2. Integration of (3.8)-(3.10) with initial datum (3.11) at times t = 0.03, 0.07, 0.1. with r = .0125 is depicted in Figure 3.2. The solution of the non-nonlocal analogue of (3.8)-(3.10) is shown in Figure 3.3. Both the simulations were obtained with a space mesh ∆x = ∆y = 1.25E -3, λ x = λ y = 0.1 and time step ∆t = 1.25E -4. The color scales in figures 3.2 and 3.3 range from the minimum to the maximum of the computed solutions U 1 and U 2 for (3.8)-(3.10) over the time interval [0, 0.1].We consider the following radii of the convolution kernel µ, see (3.10): r = 0.05, 0.0375, 0.025, 0.0125(3.12)

7 )

 7 and solutions to(3.8) with different values of r as in (3.12), all with the same initial datum(3.11), with the same space meshes ∆x = ∆y = 1.25E -3 and with a time step ∆t = 1.25E -4.non-nonlocal problem, nevertheless they suggest that a positive answer is possible. 4. Technical Details. Throughout, for any pair a, b ∈ R D , by I(a, b) we denote the segment {a + t(b -a) : t ∈ [0, 1]}. We make extensive use of the following lemmas. Lemma 4.1. Let H, K and U o be positive. Define u n+1 = (1 + H)u n + K. Then,

. 11 )

 11 Let us first consider the estimates (4.10) and (4.11): (4.10) +(4.11)

which implies that ij ( 4 . 15 )

 415 ≤ K 2 + K 3 ∆y ij u k,n i,j+1 -u k,n i-1,j ,

≤

  2K 2 ∆t + K 3 r ∆t ij u k,n i+1,j -u k,n ij ∆y + u k,n i,j+1 -u k

  3, Proposition 2.8], [11, Section 3], [12, Section 4]:

Table 3 . 1

 31 Convergence rate γ for the numerical scheme

Table 3 .

 3 2.The present numerical integrations can not prove the convergence of the nonlocal to

			u 1		u 2
	❍ ❍ ❍ ❍ ❍ t r	0.05	0.0375 0.025 0.0125	.05	0.0375 0.025 0.0125
	0.03	0.2246 0.2236 0.2222 0.2207 0.2183 0.2177 0.2168 0.2158
	0.07	0.3589 0.3577 0.3565 0.3543 0.3476 0.3461 0.3454 0.3447
	0.10	0.4479 0.4453 0.4424 0.4384 0.4135 0.4116 0.4100 0.4077

Table 3 . 2

 32 Difference of the L 1 norm of the solution (u 1 , u 2 ) to (3.

  3, Proposition 2.6]. Consider first the term ij u

	k,n+1/2 i+1,j	-u	k,n+1/2 ij	∆y, setting u	k,n+1/2 i+1,j

  xi , y j , u k,n ij , Ȃn xx f k (t n , xi+1/2 , y j , u k,n ij , Án i+1/2,j ) +∂ x ∇ A f k (t n , xi+1/2 , y j , u k,n ij , Án i+1/2,j )( Ȃn i+1,j -Ȃn ∇ A f k (t n , xi+1/2 , y j , u k,n ij , Àn i+1/2,j ) Ãn

				ij )(A n i+1/2,j -A n i-1/2 )
	=	1 2	∆x ∆x∂ 2
				ij )
		+	1 2	∆x∂ x i+1/2,j
		+	1 2	∇ 2 AA f k (t n , xi+1/2 , y j , u k,n ij , Àn i+1/2,j ) Ãn i+1/2,j ( Ȃn i+1,j -Ȃn ij )
		+	1 2	∇ A f k (t n , xi+1/2 , y j , u k,n ij , Àn i+1/2,j )(A n i+3/2,j -2A n i+1/2,j + A n i-1/2,j ) ,
	for suitable xi ∈ x i-1/2 , x i+1/2 , Ȃn i,j ∈ I(A n i-1/2,j , A n i+1/2,j ), xi , xi ∈ ]x i , xi+1 [, Án i,j , Àn i,j ∈ I( Ȃn i,j , Ȃn i+1,j ), and Ãn i+1/2,j ∈ I(A n i+1/2,j -A n i-1/2 , A n i+3/2,j -A n i+1/2 ). Passing to the modulus,
	(4.9) ≤	1 2	∆x 2 M u k,n ij

  L 1 ∈ I(A n i-1/2,j , A n i-1/2,j+1 ), which implies that Ãn ij ≤ ∆y ∂ y η L ∞ u ∆ (t n ) L 1 , ≤ ∆x ∂ x η L ∞ u ∆ (t n ) L 1 , A n i+1/2,j -A n i+1/2,j+1 -A n i-1/2,j + A n i-1/2,j+1 ≤ ∆x ∆y ∂ 2 xy η L ∞ u ∆ (t n ) L 1 . Therefore 2 (4.15) ≤ M ∆x∆y u k,n i-1/2,j+1/2 +∆y ∂ y ∂ u f k L ∞ u k,n i,j+1 -u k,n i-1,j +M ∆x∆y u k,n i-1/2,j+1/2 ∂ x η L ∞ U o L 1 +M ∆x∆y ǔk,n

	A n i+1/2,j+1 ) and A n i-1/2,j+1/2 A n i+1/2,j+1/2 -A n i-1/2,j+1/2 i-1/2,j+1/2
	so that
			ij	λ x D n ij ∆y ≤	ij	λ x (4.9) + (4.10) + (4.11) ∆y
								
						≤ ∆t	 K + K 1	ij	u k,n i+1,j -u k,n ij	∆y	
	and, adding up (4.8) to the latter expression above
	ij	u	k,n+1/2 i+1,j	-u	k,n+1/2 ij

∆y ≤ (1 + K 1 ∆t) ij u k,n i+1,j -u k,n ij ∆y + K ∆t . (4.12)

Acknowledgment: This work was partially supported by the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) ERC Grant Agreement n. 257661. The second author was partially supported by the INDAM-GNAMPA project Conservation Laws: Theory and Applications.

-a k,n i-1/2,j -b k,n i+1/2,j u k,n ij + a k,n i-1/2,j u k,n i-1,j + b k,n i+1/2,j u k,n i+1,j ≥ 0 , completing the proof.

Proof of Lemma 2.5. Refer to (4.3). Using (4.1) and (H1) in (4.3), we obtain:

Now using Lemma 2.4, (4.4) and (4.6) in (4.3), we get

Applying the same procedure for the y component gives

, and by Lemma 4.1 we have the

so that (2.7) holds with

completing the proof.

where

We now deal with ij u k,n+1/2 i,j+1

where

Following the above treatment of C n ij , we obtain:

As in the estimate above for

by convexity. Passing to Dn ij , we have:

We first estimate (4.15) using Lagrange's mean value theorem.

for suitable y j+1/2 ∈ y j , y j+1 and

.

We now use Lemma 4.2 since it holds that Ãn ij ∈ I(

Collecting the various terms above, we get

and, finally,

where, using the definitions (4.13), (4.17) and (4.18),

Entirely analogous estimates lead to

completing the proof of the lemma.

Proof of Lemma 2.8. Following [3, Lemma 2.7] and [14, Section 3], we have

Consider (4.20), since (4.21) can be treated in the same way. Using Lemma 4.2

which, together with the analogous estimate for

∆x ∆y completes the proof of the Lemma with Finally, the uniformity of the bounds proved in the Lemmas allows to pass each of these bounds to the solution, completing the proof.

Proof of Lemma 3.1. It is sufficient to verify that (H0), (H1), (H2) and (H3) hold. The latter condition is immediate. Clearly, f and g as defined in (3.2) are both of class C 2 in all their arguments.

Let now p ∈ C 3 c (R; R) be such that p(U ) = U (1 -U ) for all U ∈ [0, 1] and q ∈ C 3 c (R; R) be such that q(A) = 1 -A for all A ∈ [0, 1]. Define f (t, x, y, U, A) = p(U ) q(A) v 1 (x, y) and ĝ(t, x, y, U, B) = p(U ) q(B) v 2 (x, y) .

Clearly, by the above choices and by (3.3), f and ĝ both are of class C 2 ∩ W 3,∞ .


Paula Craciun

Josiane Zerubia

TOWARDS EFFICIENT SIMULATION OF MARKED POINT PROCESS MODELS FOR BOAT EXTRACTION FROM HIGH RESOLUTION OPTICAL REMOTELY SENSED IMAGES

Keywords: remote sensing, object extraction, marked point processes, high performance computing

Marked point processes have proven their efficiency in solving object extraction problems from high resolution optical images. However, these complex mathematical models are difficult to simulate which usually results in high computation times. A new parallel sampler has been recently developed. Nevertheless, this sampler does not yield good results in the presence of large objects. We propose modifications to the original parallel sampler to cope with such situations. Furthermore, we implement an additional mask to increase the speedup in the case of boat extraction.

INTRODUCTION

Marked point process (MPP) models have been successfully applied to object extraction in high resolution remotely sensed optical images. Applications range from flamingos and boats counting to road or building extraction [START_REF] Descombes | Stochastic Geometry for Image Analysis[END_REF]. The configuration of objects in an image can be viewed as a realization of a MPP. A complex probability density function is defined to describe both the interactions between the objects and the likelihood of such a configuration with respect to the considered image. To find the best configuration, an optimization problem has to be solved by means of sampling algorithms such as Reversible Jump MCMC (RJMCMC) [START_REF] Green | Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF]. The main ideea of RJMCMC resides in the current configuration being iteratively perturbed by a proposition kernel while an acceptance ratio is being computed for the new configuration. Nevertheless, RJMCMC poses a computational burden due to its sequential nature which results in high computation times to reach convergence. Several attempts have been made to improve the optimization procedure. A first improvement came with the implementation of Multiple Birth and Death (MBD) [START_REF] Descamps | Automatic flamingo detection using multiple birth and death process[END_REF], an optimization method which allows for multiple perturbations to be performed simultaneously. This algorithm has been further refined to Multiple Birth and Cut (MBC) [START_REF] Gamal-Eldin | Multiple birth and cut algorithm for point process optimization[END_REF] which combines the advantages of MBD with those of the Graph-Cut algorithm. The most important advantage of MBC was the reduction of the number of parameters, thus resulting in an increased convergence speed. A major breakthough came with the implementation of the first parallel version of RJM-CMC [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF]. In spite of its sequential nature, RJMCMC has been successfully parallelized by splitting the search space in smaller cells and performing the optimization locally and concurrently in each of them. This approach resulted in an efficient sampler that yielded small computation times and brought back the focus on such optimization methods which were previously considered too slow for real applications. However, the parallel implementation on GPU presented in [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF] does not take into account important considerations such as large objects located at the intersection of neighboring cells, which results in the need of inter-process communication or access to shared memory. In this paper we propose modifications to the initial sampler to account for such situations. Furthermore, we develop a fast and efficient method to discriminate between water and land areas. This latter attempt is motivated by our interest in boat detection. This paper is organized as follows: we start with a brief introduction of the general MPP framework in section 2. Section 3 describes the water/land discrimination method. Section 4 presents the weaknesses of the current parallel sampler described in [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF] and contains the modifications brought to it. Results are presented in section 5 and finally conclusions are drawn in section 6.

MARKED POINT PROCESS MODEL FOR BOAT EXTRACTION

We consider a MPP of ellipses. When seen from above, i.e. satellite images at nadir, the boats tend to have an elliptical shape. The object space W, is a bounded set in R 5 :

W = [0, X M ] × [0, Y M ] × [a m , a M ] × [b m , b M ] × [0, π]
where X M and Y M are the width and height of the input image, respectively, a m and a M are the minimum and the maximum of the semi-major axis of the ellipse, b m and b M are the minimum and the maximum of the semi-minor axis of the ellipse and ω ∈ [0, π] is the orientation of the ellipse. We are interested in a particular family of MPPs, namely the Gibbs processes [START_REF] Stoyan | Stochastic Geometry and its Applications[END_REF]. Denoting the observed image with y, the density of the considered MPP is given by:

f θ (X = x|y) = 1 c(θ|y) exp -U θ (x,y)
where c(θ|y) is a normalizing function given by:

c(θ|y) = Ω exp -U θ (x,y)µ(dx) .
Here, θ is a parameter vector which allows the model to be flexible and fit different images. It has to be adjusted according to the given image. µ(•) is the intensity measure of the reference Poisson process, U θ (x, y) is called the energy, Ω is the configuration space and x denotes the object configuration.

Using the Maximum A Posteriori (MAP) criterion, the most likely configuration which allows the extraction of objects corresponds to the global minimum of the total energy U θ (x, y):

x ∈ Argmax x∈Ω f θ (X = x|y) = Argmin x∈Ω [U θ (x|y)].
Parameter estimation techniques are presented in detail in [START_REF] Chatelain | Parameter estimation for a marked point process. application to object extraction from remote sensing images[END_REF][START_REF] Ben Hadj | Parameter estimation for a marked point process within a framework of multidimensional shape extraction from remote sensing images[END_REF]. Once the parameter vector θ is determined, a solution to the resulting optimization problem has to be found. For the extraction of boats outside the harbor area we use the same model presented in [START_REF] Descamps | Automatic flamingo detection using multiple birth and death process[END_REF] for flamingo extraction. In the case of harbor images, details on the model used are presented in [START_REF] Craciun | Unsupervised marked point process model for boat extraction in harbors from high resolution optical remotely sensed image[END_REF].

WATER/LAND DISCRIMINATION

A first step to improve the performance of this framework is to limit the search for objects to the area where they are expected to appear. Since our main interest is to extract boats, it is meaningful to search for them only in the areas where water is present. This motivates our attempt to discriminate between water and land area. The water area can be clearly characteristics, as shown in Figure 1 (left). The difficulty lies in finding appropriate features to distinguish between these two area types.

In this paper we focus on a simple approach to make the separation between water and shadows. We are interested in extracting the water area under the condition that the number of misdetection should be zero, whereas the number of false alarms should be minimal. Given the low radiometric values of water, a threshold is determined by means of bimodal histogram splitting. This approach is motivated by our interest to distinguish between water and non-water area. 1 shows the difference in feature values for the two types of areas. Based on the results presented in Table 1, we can devise an algorithm to extract the water area, similar to the one presented in [START_REF] Dare | Shadow analysis in high-resolution satellite imagery of urban areas[END_REF]. The main difference is our emphasis on size. Thus, after we compute the size and the radiometric variance for each connected component identified in the threshold image, we retain only those components that have both large size and low radiometric variance. Figure 1 (right) presents the results obtained. We can observe that the water area has been correctly and entirely identified.

For the purpose of boat extraction, we erode the obtained image with a circular structuring element, twice the size of the boats we are interested in extracting.

HIGH PERFORMANCE COMPUTING

The first major step towards a time-efficient approach to simulate MPPs came with the development of the parallel sampler devised in [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF]. The main idea in the 2D case is to divide the search space using a quadtree K. The cells of the quadtree are divided into 4 independent sets called mic-sets and denoted S mic . At each iteration, one mic-set is selected and the optimization step is performed in parallel within all the cells, c, contained in it. The general proposition kernel Q is formulated as a mixture of uniform sub-kernels Q c,γ , where γ ∈ Γ = {birth and death, translation, rotation, scale}. The probability of each proposition kernel is computed as q c,t = P r(γ) # cells in K . The computational efficiency of the sam- pler has been proven for a large number of applications. The sampler in [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF] makes use of GPU computing to perform the optimization. GPU computing has confirmed its usefulness in the case of data-parallel processing. In other words, if there is no need for individual processors to communicate and share data between each other, nor to access shared memory, GPU computing outstands multi-core parallelization approaches. 2. Ellipses are allowed to cross the boundary of a cell but this approach leads to several detections of the same object (see Figure 2 (c))

(a) (b) (c) (d) (e) (f) (g) (h)
To overcome this problem, we implement a multi-core version of the sampler with a shared memory between the cores.

Algorithm 1 Parallel sampler 1. Initialize X 0 = x 0 and t = 0;

2. Compute the water mask;

3. Compute the data-driven space partitioning tree K and truncate it based on the water mask;

4. At iteration t with X t = x:

• Choose a mic-set S mic ∈ K and a kernel type γ ∈ Γ according to the probability c∈Smic q c,γ

• For each cell c in S mic :

-Perturb x in cell c to a configuration x ′ according to Q c,γ (x → •); -Retrieve the configuration z from the neighboring cells; -Compute the Green ratio:

R = Q c,γ (x ′ ∪ z → x ∪ z) Q c,γ (x ∪ z → x ′ ∪ z) exp U (x ∪ z) -U (x ′ ∪ z) T t
-Choose X t+1 = x ′ with probability min(1, R) and X t+1 = x otherwise;

• Update T t+1 = αT t (in our tests α = 0.95).

At each iteration, for each cell in a chosen mic-set S mic , the configurations within the neighboring cells are taken into con-# CPU's Computation times for Figure 3 A multi-core computer with 18 Intel Xeon(R) 2.30 GHz physical CPU's and 6 virtual CPU's was used for the computations. The computation time needed for each of the two images with respect to the number of CPU's and the use of the water mask are presented in Table 2. A birth map was considered for each extraction, independent of the existence of the water mask. As expected, the extraction times drop considerably with the increase in the number of CPU's. Results show that an increase in the number of objects present in the scene translates into a considerable increase in computation time. Finally, the usefulness of the water mask is proven by an additional decrease in computation time in both cases.

RESULTS AND DISCUSSIONS

CONCLUSIONS AND FUTURE WORK

In this paper, we have presented two approaches for increasing the efficiency of simulating marked point process models for boat extraction. We discriminated between water and land area to reduce the search space. Furthermore, we have modified the parallel sampler devised in [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF] to cope with large objects and implemented it on a multi-core computer with shared memory access between the processors. This approach was motivated by the necessity to access the configurations in the neighborng cells. This resulted in improved extraction results and low computation times. Future work will include an extension of this framework to time series.

Fig. 1 .

 1 Fig. 1. left: Harbor image c Astrium-EADS; right: Identified water area identified as a large area of low radiometric values. Nevertheless, the shadows of tall buildings within cities exhibit similar

Fig. 3 .

 3 Fig. 3. (a), (e) Image of boats in and outside a harbor c Astrium-EADS; (b),(f) Water/Land discrimination results; (c),(g) Extraction results; (d),(h) Close-ups on extraction results.

Fig. 2 .

 2 Fig. 2. (a) A large boat split in two due to the space partitioning; (b), (c) Errors which arise when processing cells independently of each other However, this is not our case. One object can be contained in two neighboring cells, as shown in Figure 2 (a). Imposing data-parallel processing in each cell independently results in two possible types of detection errors: 1. The object is split at cell boundary and thus, two smaller objects are detected (see Figure 2 (b));

Figure 3

 3 Figure 3 presents the extraction results.Figure 3 (a) of size 1620 × 1450 represents an image of a harbor with a large number of objects, while Figure 3 (e) of size 2862 × 2676 represents a large area outside of a harbor with only a few boats. Figures 3 (b) and (f) show the area used for boat extraction after identifying the water area.A multi-core computer with 18 Intel Xeon(R) 2.30 GHz physical CPU's and 6 virtual CPU's was used for the computations. The computation time needed for each of the two images with respect to the number of CPU's and the use of the water mask are presented in Table2. A birth map was considered for each extraction, independent of the existence of the water mask. As expected, the extraction times drop considerably with the increase in the number of CPU's. Results show that an increase in the number of objects present in the scene translates into a considerable increase in computation time. Finally, the usefulness of the water mask is proven by an additional decrease in computation time in both cases.

Table 1 .

 1 This will identify areas of low radiometric values, which include both water and shadows. Two main features can be used to differentiate Comparison

	Component type	Size	Mean	Variance
	Shadow 1	8641	12.8962	1.03389
	Shadow 2	8211	13.0898	1.22846
	Shadow 3	9986	13.0068	1.00696
	Water	969675	13.2784	0.377285

of characteristics of water and shadow components between water and shadows: size of the area and radiometric variance within it. Table

Table 2 .

 2 Comparison of characteristics of water and shadow components sideration. Note that this does not give rise to synchronization problems when accessing the shared memory since processors only read information in the neighboring cells and do not modify it. The modified parallel sampler is detailed in Algorithm 1.

			(a)	Computation times for Figure 3 (d)
		Water mask	No water mask	Water mask	No water mask
	1	1h 38min 57sec	2h 4min 5sec	1h 27min 48sec	2h 9min 35sec
	2	1h 20min 23sec	1h 45min 41sec	33min 43sec	40min 21sec
	4	53min 37sec	1h 01min 12sec	11min 29sec	15min 26sec
	8	44min 7sec	49min 23sec	4min 59sec	5min 46sec
	16	20min 19sec	23min 1sec	3min 53sec	4min 31sec
	24	9min 41sec	12min 54sec	2min 32sec	3min 8sec

The authors thank Astrium-EADS for providing the data and for the partial funding of this research.