Denise Aregba-Driollet

Frédéric Krantz
email: krantz@math.u-bordeaux1.fr

DISCRETE KINETIC FINITE VOLUME SCHEMES FOR MULTIDIMENSIONAL SYSTEMS OF CONSERVATION LAWS ON UNSTRUCTURED GRID

Keywords: Conservation Laws, BGK Models, Kinetic Schemes, Finite Volume Schemes, Unstructured Meshes

We construct a family of Riemann solver free finite volume schemes for systems of conservation laws. The method is based on a relaxation approximation which takes the form of a BGK system. We prove convergence for unstructured grids in the scalar case. We present 2D computations for Euler equations.

INTRODUCTION

In this paper we construct some finite volume discretisations for hyperbolic systems of conservation laws in two space dimensions, in an open subset Ω ⊂ R 2 with C 2 boundary ∂Ω :

∂ t U + ∂ x 1 A 1 (U ) + ∂ x 2 A 2 (U ) = 0 in ]0, +∞[×Ω, (1)
where U (x, y, t) ∈ R p and A 1 , A 2 are smooth functions with values in R p . We denote u 0 an initial value for this system:

U (x 1 , x 2 , 0) = U 0 (x 1 , x 2 ). (2)
D. Aregba-Driollet and R. Natalini 2 designed a class of numerical schemes based on a discrete kinetic approximation of the Cauchy problem for [START_REF] Aregba-Driollet | Kinetic approximation of a boundary value problem for conservation laws[END_REF]. These schemes were extended to the initial-boundary value problem by D. Aregba-Driollet and V. Milišić 1 . All these schemes were designed on cartesian grids, here we extend these works to unstructured grids.

To approximate the solutions of (1), we choose a semilinear hyperbolic system of BGK type :

   ∂ t f ǫ l + λ l1 ∂ x 1 f ǫ l + λ l2 ∂ x 2 f ǫ l = 1 ǫ (M l (U ǫ ) -f ǫ l ) l ∈ {1, • • • , N } f ǫ l (0, x) = M l (U 0 (x)) (3)
where (λ ld ), 1 ≤ l ≤ N , d = 1, 2 are fixed real numbers, U ǫ is defined by U ǫ = N l=1 f ǫ l , and ǫ is a positive parameter. The link between (1) and ( 3) is done via the following compatibility conditions: Definition 1.1. Let I ⊆ R p be a fixed domain. A Lipschitz continuous function M = (M l ) 1≤l≤N : I → (R p ) N is a (local) Maxwellian Function for [START_REF] Aregba-Driollet | Kinetic approximation of a boundary value problem for conservation laws[END_REF] and with respect to I if the following (compatibility) conditions are verified for all U ∈ I:

N l=1 M l (U ) = U (4) N l=1 λ ld M l (U ) = A d (U ) d = 1, 2. ( 5
)
It is easy to see formally that if conditions (4)(5) are satisfied and if U ǫ tends to a limit function U as ǫ → 0 then U is a solution of (1). In the following we shall assume that the function M is a Maxwellian Function.

Actually one needs for some stability conditions to ensure this convergence, namely:

Definition 1.2. The functions M = (M l ) 1≤l≤N is Monotone Non Decreasing (MND) on I if all the eigenvalues of M ′ l (U ) are real and non negative for all U ∈ I and all l = 1, . . . , N . In particular in the scalar case p = 1 it means that M ′ l (U ) ≥ 0 on I for all l.

It can be shown that the MND property is closely related to Liu's subcharacteristic condition interlacing the eigenvalues of system (1) and the characteristic velocities λ ld , see Aregba-Driollet and Milišić 1 . In the scalar case, suppose that the maxwellian function is MND on [-U 0 ∞ , U 0 ∞ ]. Then U ǫ → U , the unique entropy solution of the Cauchy problem (1-2), see Natalini 8 . This result has been extended to the initial-boundary problem in Milišić 7 . In the case of systems with a strictly convex entropy, see Bouchut 3 for the relationship between the MND condition and the entropy properties of the BGK approximation.

In this paper, for the numerical experiments we focus our attention on the following 4 orthogonal velocities model because the needed information for applying the MND property is minimal: one has just to know the spectral ray of the jacobian of the fluxes. The velocities vectors are

-→ λ 1 = λ x 1 0 -→ λ 2 = λ y 0 -1 -→ λ 3 = λ x -1 0 -→ λ 4 = λ y 0 1 ,
and the Maxwellian Functions are

M 1 (U ) = 1 4 [U + 2A 1 (U ) λ x ] M 2 (U ) = 1 4 [U - 2A 2 (U ) λ y ] M 3 (U ) = 1 4 [U - 2A 1 (U ) λ x ] M 4 (U ) = 1 4 [U + 2A 2 (U ) λ y ]
Here the stability condition reads

λ x > 2 sup U ∈I |σ(A ′ 1 (U ))| λ y > 2 sup U ∈I |σ(A ′ 2 (U ))|,
where σ(•) is the spectrum of the jacobian matrices. The plan of the paper is the following: in section 2 we construct the schemes and give a convergence result for the scalar case. In the last section numerical experiments are performed on Euler equations of gas dynamics.

KINETIC FINITE VOLUME SCHEMES

In this part we design a numerical scheme for the relaxing semilinear system (3) associated with (1), [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. We present here a cell-center type finite volume method on unstructured grid. Our aim is to construct a scheme of the following form :

         U n+1 i = U n i - ∆t |C i | e=C i ∩C j |e|g ij (U n i , U n j , n e ) U 0 i = 1 |C i | C i U 0 (x 1 , x 2 ) dx 1 dx 2 . (6)

Geometry and notations

We shall collect here the assumptions concerning geometry and the notations used in this work. Definition 2.1. Let a k-polygon be a closed, convex polygon with k vertices. The set T h := {C i /C i is a k-polygon for i ∈ I ⊆ N} (where I ⊆ N is an index set) is called an unstructured grid of Ω ⊂ R n if the two following properties are satisfied :

1. Ω = i∈I C i 2. For two different C i and C j , we have C i ∩ C j = ∅, or C i ∩ C j = a common vertex of C i , C j or C i ∩ C j = a common edge of C i , C j .
Definition 2.2. Let (C i ) i∈I denote an unstructured grid of R 2 . We shall use the following notation :

λ l (n e ) = λ l • n e . h = sup i∈I h i , where h i is the exterior diameter of C i .
U n i : approximation of the exact solution U on C i at time n∆t.

Let us make the following assumptions on the geometry :

Assumption 2.1. We assume that there are two constants c 1 et c 2 such that

0 < c 1 ≤ ∆t h ≤ c 2 (7)
if ∆t, h → 0. Moreover, we assume that there exists a constant c V > 0 such that

sup i h 2 |C i | ≤ c V (8)

Design of the schemes

The timestep is denoted by ∆t and discrete times are t 0 = 0 and t n+1 = t n + ∆t. Each U n i is an approximation of the mean value of U (•, t n ) on the cell C i . At time t 0 = 0, if U 0 is the initial value, we take

U 0 i = 1 |C i | C i U 0 (x 1 , x 2 )dx 1 dx 2 .
For the microscopic (kinetic) variables, we use as initial conditions :

f 0 l,i = M l (U 0 i ).
We use a fractional step method : first we solve the homogeneous hyperbolic problem on [t n , t n+1 ]:

∂ t f ǫ l + λ l1 ∂ x 1 f ǫ l + λ l2 ∂ x 2 f ǫ l = 0 l ∈ {1, • • • , N } f ǫ l (x 1 , x 2 , t n ) = f ǫ,n (x 1 , x 2 ) (9)
As the system is diagonal, we may consider each equation separately. That is, we choose a conservative scheme of flux ϕ l for the l th -equation. Then we obtain a scheme under conservative form :

f ǫ,n+1/2 l,i = f ǫ,n l,i - ∆t |C i | e=C i ∩C j |e|ϕ l (f ǫ,n l,i , f ǫ,n l,j , n e ).
Next, taking f ǫ,n+1/2 l,i

as initial condition at time t n , we solve the ordinary differential system :

∂ t f ǫ l,i = 1 ǫ (M l (U ǫ,n i ) -f ǫ l,i ) f ǫ l,i (t n ) = f ǫ,n+1/2 l,i
Thanks to the compatibility condition (4), we may solve explicitly this system. Then we write the solution :

f ǫ,n+1 l,i = (1 -e -∆t ǫ )M l (U ǫ,n+1/2 i ) + e -∆t ǫ f ǫ,n+1/2 l,i
, where

U ǫ,n+1/2 i = N l=1 f ǫ,n+1/2 l,i . Note that U ǫ,n+1 = U ǫ,n+1/2 .
We have constructed a family of numerical scheme for the semilinear system (3), that differ by the choice of the homogeneous scheme.

When ǫ → 0, we obtain the relaxed limit of the scheme :

           f n l,i = M l (U n i ) f n+1/2 i = H ∆ (∆t)f n i U n+1 i = N l=1 f n+1/2 l,i
where H ∆ (∆t) represents the numerical scheme applied to the kinetic equations [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV: A new approach to numerical convection[END_REF].

In order to obtain a finite volume formulation (6), we have to choose for each kinetic equation a numerical scheme, which can be defined by its numerical flux, called ϕ l for the l th equation. According to the form of the relaxed scheme, the numerical scheme is the following :

           U n+1 i = U n i - ∆t |C i | e=C i ∩C j N l=1 |e|ϕ l (M l (U n i ), M l (U n j ), n e ) U 0 i = 1 |C i | C i U 0 (x 1 , x 2 )dx 1 dx 2 . (KFVS)
This scheme is of the form (6) with

g ij (U i , U j , n e ) = N l=1 ϕ l (M l (U i ), M l (U j ), n e ). ( 10
)
3 Convergence in the scalar case

In the scalar case p = 1 we are able to prove convergence of schemes (KFVS) following the lines of Kröner and Rokyta's work 6 . A preliminary result is the existence of kinetic entropies: Theorem 3.1 (Bouchut 3 ). Let η a strictly C 1 -convex entropy, the condition M ′ l (u) ≥ 0, ∀u ∈ I, ∀l ∈ {1, . . . , N } is equivalent to existence of kinetic entropies S l,η defined in the following way :

S l,η (f ) = R 1 2 (|f -M l (k)| -|M l (k)|)η ′′ (k)dk + 1 2 f (η ′ (-∞) + η ′ (∞)).
They own the following properties :

Functions S l,η are C 2 -convex in [M l (-u 0 ∞ ), M l ( u 0 ∞ )]. ( 11
) ∀w ∈ I, S ′ l,η (M l (w)) = η ′ (w), ( 12
) N l=1 λ ld S l,η (M l (w)) = A η,d (w). ( 13
) N l=1 S l,η (M l (w)) = η(w) -η(0). (14)
The identity (12) holds in the sense of subdifferentials : for all f ∈ [M l (-u 0 ∞ ), M l ( u 0 ∞ )], and for all w ∈ I,

S l,η (M l (w)) -S l,η (f ) -η ′ (w)(M l (w) -f ) ≤ 0. ( 15
)
It ensures, for each kinetic equation, the existence of an entropy pair (S l,η , F l ), where F ld (s) = λ ld S l,η (s).

We make the following assumptions on the kinetic numerical fluxes: Assumption 3.1. For all l ∈ {1, . . . , N }, let ϕ l be a numerical flux of linear three-points scheme in conservation form, consistent with λ l f l • n e i.e.

ϕ l (u, u, n e ) = λ l u • n e (16)
We assume that ϕ l is Lipschitz-continuous. In particular, suppose that for all M > 0 there is a constant

C ϕ l (M ) such that for all u, u ′ , v, v ′ ∈ [-M, M ] |ϕ l (u, v) -ϕ l (u ′ , v ′ )| ≤ C ϕ l (M )(|u -u ′ | + |v -v ′ |) (17)
and that ϕ l is conservative, i.e.

ϕ l (u, v, n e ) = -ϕ l (v, u, -n e ) (18)
Moreover, we assume that ϕ l is monotone :

∂ ∂u ϕ l (u, v) ≥ 0 ≥ ∂ ∂v ϕ l (u, v) (19)
Assumption 3.2. Let η : R → R be Lipschitz-continuous and convex, and let (S l,η , F l ) be an entropy pair for [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV: A new approach to numerical convection[END_REF]. We assume that for all l ∈ [1, N ] there exists a numerical entropy flux G l which is :

(i) consistent with F l (u) • n e = λ l • n e S l,η (u) i.e G l (u, u, n e ) = λ l • n e S l,η (u) (20)
(ii) Lipschitz-continuous, such that for all M > 0 there is a constant

C G l (M ) such that for all u, u ′ , v, v ′ ∈ [-M, M ] ⊂ R |G l (u, v) -G l (u ′ , v ′ )| ≤ C G l (M )(|u -u ′ | + |v -v ′ |) (21) (iii) conservative, i.e. for all j ∈ N i G l (u, v, n e ) = -G l (v, u, -n e ) ( 22
)
and (iv) satisfies the kinetic compatibility condition

∂G l ∂v (p, q, n e ) = S ′ l,η (q) ∂ϕ l ∂q (p, q, n e ) ( 23
)
Let us note that if we choose the (linear) Godunov numerical flux for each ϕ l then assumptions 3.1 and 3.2 are verified. On another hand this choice is the optimal one for first order, as it is the less diffusive one. Lemma 3.1. For all i ∈ I, we have :

e⊂C i e=C i ∩C j |e| n e = 0
A key result to prove convergence is the L ∞ stability of the scheme, which implies monotony and entropy dissipation.

Proposition 3.1 (L ∞ stability). Let u 0 ∈ L ∞ (R 2 ) ∩ L 1 (R 2 ). Let u n
i defined by the kfv scheme (KFVS) such that the fluxes ϕ l satisfy the assumptions 3.1. We denote

M = u 0 L ∞ (R 2 ) and C(M ) = l C l (M l ( u 0 ∞ )).
We assume that the CFL-type condition

∆t h ≤ 1 k • C(M )c V , (24)
where c V is defined as in [START_REF] Natalini | A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws[END_REF], is fullfilled. We assume that the maxwellian functions M l are MND on [-M, M ]. Then we have, for all n ∈ N,

u n L ∞ (R 2 ) ≤ u 0 L ∞ (R 2 ) . (25)
Proof. The KFVS method is given by :

u n+1 i = u n i - ∆t |C i | e⊂C i e=Γ ij |e| g ij (u n i , u n j , n e ) = u n i - ∆t |C i | e⊂C i e=Γ ij |e| (g ij (u n i , u n j , n e ) -g ij (u n i , u n i , n e )),
by lemma 3.1. Then

u n+1 i = u n i - ∆t |C i | e⊂C i e=Γ ij D n ij (u n i -u n j ).
where

D n ij = |e| g ij (u n i , u n j , n e ) -g ij (u n i , u n i , n e ) u n i -u n j for u n j = u n i (26)
and

D n ij = 0 for u n j = u n i . (27)
Next we can write

u n+1 i = u n i     1 - ∆t |C i | e⊂C i e=Γ ij D n ij     + ∆t |C i | e⊂C i e=Γ ij D n ij u n j .
Let us assume that for a fixed n ≥ 0, u n L ∞ (R 2 ) ≤ M . By assumption 3.1 and monotony of the maxwellian functions 0 ≤ D n ij ≤ C(M )h, and using the CFL-type condition (24) and assumption 2.1

∆t |C i | e⊂C i e=Γ ij D n ij ≤ k ∆ t min i∈I |C i | C(M )h ≤ 1.
Finally we have, due to the convex combination

|u n+1 i | ≤ |u n i |     1 - ∆t |C i | e⊂C i e=Γ ij D n ij     + ∆t |C i | e⊂C i e=Γ ij D n ij |u n j | ≤ u n L ∞ (R 2 )     1 - ∆t |C i | e⊂C i e=Γ ij D n ij     + ∆t |C i | e⊂C i e=Γ ij D n ij u n L ∞ (R 2 ) ≤ u n L ∞ (R 2 ) ≤ M.
Then we get u n+1 L ∞ (R 2 ) ≤ M and we can complete the proof by induction.

Kröner and Rokyta 6 proved convergence for upwind finite volume schemes. The above stability result allows us to prove monotony and to define an entropy numerical flux, which satisfies the compatiblity condition (23). We are then able to prove the convergence theorem: Theorem 3.2. We consider the Cauchy problem (1)-( 2) and the approximate solution {u h } defined by the kinetic finite volume (first-order) method (KFVS). We suppose that the maxwellian functions are MND on the interval [-u 0 ∞ , u 0 ∞ ] and we take M > u 0 ∞ . We suppose that the conditions [START_REF] Milišić | Stability and convergence of discrete kinetic approximations to an initialboundary value problem for conservation laws[END_REF], [START_REF] Natalini | A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws[END_REF] for the triangulation hold, and that the CFL condition (24)

∆t h ≤ 1 k • C(M )C V , (28)
is satisfied. Then the method (KFVS) converges weak-⋆ to the unique entropy solution in the sense of Kružkhov of the Cauchy problem (1)- [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF].

u h → u weak -⋆, u ∈ L ∞ (R 2 × R + ).

Extension to second order in space

We focus now on the extension to second order of these family of kinetic schemes. The semi-discrete version of these scheme is :

|C i | dU i dt + e=C i ∩C j |e|g(U i , U j , n e ) = 0, C i ∈ T h . (29)
For internal edges, the fonction g(., ., n e ) is equal to the kinetic flux as shown above. The extension to second order accury consists in replacing the states U i and U j in relation (29) by nonlinear extrapolations U i,e and U j,e on each side of the boundary, constructed from the previous states as explained in the sequel. We follow here the method given by where α i (z) is a limiting coefficient satisfying the following conditions :

0 ≤ α i (z) ≤ 1, z ∈ {ρ, ρu, ρv, p}, C i ∈ T h k(z i -m i (z)) ≤ α i (z)∇z| i • (y i,e -x i ) ≤ k(M i (z) -z i ), ∀e ⊂ ∂C i , C i ∈ T h . (38)
Then α i (z) is chosen as large as possible and less than or equal to 1 in order to satisfy the constraints :

α i (z) = min 1, min(M i (z) -z i , z i -m i (z)) max{|∇z| i • (y i,e -x i )|, e ⊂ ∂C i } . ( 39
)
For k = 1, we recover the initial limiter proposed by Van Leer 9 . When k = 1/2, we obtain the minmod limiter proposed by Harten 5 . The intermediate value k = 3/4, which was named STS by Dubois and Michaux is a good choice between the two others.

Second order in time

When all values z i,e are known for all control voumes, all faces and all fields, extrapolated states U i,e are naturally defined by going back to the conservatives variables. Then we introduce these states as arguments of the flux function g.. and obtain by this way a new system of ordinary differential equations :

|C i | dU i dt + e⊂∂C i
|e|g(U i,e , U i,j(e) , n e ) = 0, C i ∈ T h .

To have a second order accuracy in time, we use a two-steps Runge-Kutta scheme :

|C i | ∆t (U (1) i
-U n i ) + e⊂∂C i |e|g(U n i,e , U n i,j(e) , n e ) = 0, C i ∈ T h , (41)

|C i | ∆t (U (2) i -U (1)
i ) +

e⊂∂C i |e|g(U (1)
i,e , U

i,j(e) , n e ) = 0,

C i ∈ T h , (42)
U n+1 i = 1 2 (U (2) i + U n i ), C i ∈ T h . ( 43
)

NUMERICAL EXPERIMENTS

We approximate a weak solution U of Euler's system of gas dynamics in two dimensions

U =     ρ ρu ρv E     , A 1 (U ) =     ρu ρu 2 + p ρuv (E + p)u     , A 2 (U ) =     ρv ρuv ρv 2 + p (E + p)v     , (44)
where we have to add the equation of state for a γ-law gas

E = p γ -1 + 1 2 ρ(u 2 + v 2 ). (45)

Denise Aregba-Driollet, Frédéric Krantz F. Dubois and O. Michaux 4 , which is a generalization of Van Leer's 10 MUSCL scheme for unstructured meshes.

Let us introduce the set V i of neighbouring cells of a given cell C i :

We also introduce the point y i,e on the interface e that links the barycenters x i and x i,j(e) :

y i,e =(1θ i,e )x i + θ i,e x i,j(e) , y i,e ∈ e, e ⊂ ∂C i , C i ∈ T h .

Next, for z equal to one scalar variable of the family :

we evaluate a mean value z i,e on the interface e :

and the gradient ∇z| i of field z(•) in volume C i with a Green formula :

An extrapolation of field z(•) is given by

but the variation ∇z| i • (y i,ex i ) has to be limited. To this end we define the minimum m i (z) and the maximum M i (z) of field z in the neighbouring cells :

If the value z i is extremum among the neighbouring ones, we impose that the extrapolated value z i,e is equal to the cell value z i :

We introduce a nonlinear extrapolation of the field z(•) between center x i and boundary face y i,e (e ⊂ ∂C i ) :

Denise Aregba-Driollet, Frédéric Krantz

The jacobian matrix A ′ 1 (U ) has the four eigenvalues

As in one dimension, the sound speed is c = γp/ρ. The eigenvalues in the y-direction are similar, with the roles of u and v reversed.

A mach 3 wind tunnel with a step

This example is a classical test example, which has been proven to be useful for a large number of numerical schemes. A very complete presentation can be found in the paper of Woodward and Colella 11 .

The test case describes a Mach 3 flow in a wind tunnel. The tunnel is 1 length unit high and 3 length units long. The step is 0.2 length units high and is located 0.6 length units from the left-hand end of the tunnel. The walls are reflective.

At the left, we impose a supersonic inflow boundary condition, while at the right side we impose an outgoing condition. Initially the wind tunnel is filled with a gamma-law gas, with γ = 1.4, which everywhere has density 1.4, pressure 1.0, x-velocity 3.0 and y-velocity 0.0.

On the boundary we impose reflecting boundary condition that is

We use four different schemes, all constructed on the four velocities model and the upwind (Godunov) flux. Every solution is displayed with 30 isolines at time t=3.5. The first-order scheme (figure 1) gives naturally the least accuracy, but the Mach stem is present.

With the the second-order scheme (figure 2), the shocks are thinner but a numerical instability is evident near the bottom wall, and behind the Mach stem.

In order to improve the previous scheme, where the velocities are the same in each cell, we fixed these velocities by cell. The result is clearly better (figure 3), This scheme gives a more accurate reprensation of the general shape and position of the shocks .

The last scheme (figure 4) is the same as the previous combined with a second-order in time Runge -Kutta scheme. The results are quite similar.

Double mach reflection of a strong shock

This test describes the reflection of a planar Mach shock in air from a wedge. The setup is of a Mach 10 shock which initially makes a 60 degrees angle with a reflecting wall. The undisturbed air ahead of the shock has a density of 1.4 and pressure of 1.0. The computational domain is [0, 4] × [0, 1] and the reflecting wall lies along the bottom of the domain, beginning at x = 1/6. See Woodward and Colella 11 for more precision about this case.

For this test-case we use three different schemes constructed with the four velocities and the upwind flux.

First, the first-order scheme (figure 5) gives the least accurate results. The jet formed by the double Mach reflection is unresolved. The strong shocks are too large.

The second scheme is the second-order scheme with velocities fixed by cells (figure 6). The results are superior to those of first-order. The jet is resolved and the weak shock is better decribed. However there are some unphysical structures.

The third scheme is the same as the previous combined with a second-order in time Runge-Kutta scheme (figure 7). The results are better, the shocks are thinner and the jet is quite well described. However the weak shock generated at the kink in the main reflected shock is quite broad.

CONCLUSIONS

The kinetic models (3) offer an alternative to the Godunov approach for the construction of numerical approximations of hyperbolic systems of conservation laws. The higher order extension in space and in time is easy and this approach is very flexible with regard to the choice of the underlying kinetic model as well as to the one of the discretisation parameters.

In the scalar case, the convergence is proved. We have performed numerical experiments on 2D Euler system with unstructured grids. The chosen kinetic model is as simple as possible, and the implementation needs for minimal information on the spectral properties of the jacobian of the fluxes, namely their spectral ray. We observe nevertheless that the numerical results are qualitatively very satisfying: all the essential structures of the solutions are present.

Other choices of the maxwellian function may need for more precise information and lead also to more accurate results but in all cases one has just to solve linear transport equations.


Denise Aregba-Driollet
email: aregba@math.u-bordeaux1.fr

THE RIEMANN PROBLEM FOR KERR EQUATIONS AND NON-UNIQUENESS OF SELFSIMILAR ENTROPY SOLUTIONS

Keywords: 2000 Mathematics Subject Classification. Primary: 35L65 ; Secondary: 35A02, 78A60 Systems of conservation laws, Riemann problem, Lax entropy solutions, Liu's criteria, nonlinear Maxwell's equations

We solve the Riemann problem for a nonlinear full wave Maxwell system arising in nonlinear optics. This system is hyperbolic, some eigenvalues have non-constant multiplicity and are neither genuinely nonlinear, nor linearly degenerate. In a particular 2×2 reduced case, we are able to exhibit two distinct selfsimilar entropy solutions. We compute the amounts of entropy dissipation and compare them.

1. Introduction. In nonlinear optics, the propagation of electromagnetic waves in a crystal can be modelized by the so-called Kerr model, which consists of Maxwell's equations ∂ t D -curlH = 0, ∂ t B + curlE = 0, with divD = divB = 0 and the constitutive relations

B = µ 0 H D = D(E) = ε 0 (1 + ε r |E| 2 )E.
Here µ 0 , ε 0 are the free space permeability and permittivity and ε r is the relative permittivity, see [START_REF] Shen | The Principles of Nonlinear Optics[END_REF] for further details. The model is a 6 × 6 system of conservation laws in the unknown u = (D, H):

∂ t D -curlH = 0, ∂ t H + µ -1 0 curl(P(D)) = 0 ( 1
)
where P is the reciprocal function of D. Denoting q(e) = ε 0 (e + ε r e 3 ), e ∈ R, p = q -1 ,

we have

E = P(D) = D ε 0 (1 + ε r p 2 (|D|))
.

As proposed in [START_REF] Carbou | Relaxation approximation of some nonlinear Maxwell initialboundary value problem[END_REF] we also introduce the one dimensional model satisfied by solutions D(x, t) = (0, d(x, t), 0), H(x, t) = (0, 0, h(x, t)) and x = x 1 ∈ R. In that framework the solutions of Kerr model 1 satisfy the following p-system:

∂ t d + ∂ x h = 0, ∂ t h + µ -1 0 ∂ x p(d) = 0. (2)
As p ′ > 0 it is strictly hyperbolic but the properties of the function p differ from the ones which appear in the general framework of gas dynamics or viscoelasticity.

Here: p(0) = 0, p ′ > 0, and p is strictly convex on ] -∞, 0], strictly concave on [0, +∞[. Known existence results for system 1 are related to strong solutions, see [START_REF] Racke | Lectures on nonlinear evolution equations. Initial value problems[END_REF], [START_REF] Kanso | Sur le modèle de Kerr-Debye pour la propagation des ondes électromagnétiques[END_REF] and references therein. A first insight into weak solutions, which is also useful in the aim of designing numerical schemes, is the study of the Riemann problem: a direction ω ∈ R 3 , |ω| = 1, and u ± ∈ R 6 being fixed, one looks for the solution of system 1 with initial data

u(x, 0) = u -if x • ω < 0, u + if x • ω > 0. ( 3
)
This work is devoted to the resolution of this problem and to the link between the solutions for the full model 1 and the ones for the reduced system 2. We point out the fact that we do not suppose that the initial data are divergence free (ie that (D + -D -) • ω = 0 and (H + -H -) • ω = 0) because in numerical applications, this condition is not exactly satisfied in general.

The electromagnetic energy is a mathematical entropy [START_REF] Carbou | Relaxation approximation of some nonlinear Maxwell initialboundary value problem[END_REF]. In the reduced 2 × 2 case it reduces to the classical entropy of the p-system.

The characteristic fields of system 1 have been described in [START_REF] Aregba-Driollet | Kerr-Debye relaxation shock profiles for Kerr equations[END_REF], the following proposition summarizes the results: Proposition 1. [START_REF] Aregba-Driollet | Kerr-Debye relaxation shock profiles for Kerr equations[END_REF] The Kerr system 1 is hyperbolic diagonalizable: for all ω ∈ R 3 , |ω| = 1, the eigenvalues are given by

λ 1 ≤ λ 2 = -λ < λ 3 = λ 4 = 0 < λ 5 = λ ≤ λ 6 = -λ 1 (4)
where, denoting c = (ε 0 µ 0 ) -1 2 the light velocity:

λ 2 1 = c 2 1 + ε r |E| 2 , λ 2 = c 2 1 + ε r (|E| 2 + 2(E • ω) 2 ) (1 + ε r |E| 2 )(1 + 3ε r |E| 2 ) . (5
Ω(ω) = {(D, H) ∈ R 6 ; ω × D = 0}.
The system is not strictly hyperbolic and the second and fifth characteristic fields are neither genuinely nonlinear, nor linearly degenerate. Hence, Lax' wellknown result [START_REF] Lax | Hyperbolic systems of conservation laws[END_REF] does not apply. Here for |u +u -| small enough, we construct a "Lax' solution" and prove that such a construction is unique.

In the 2×2 case, the system 2 is strictly hyperbolic but similarly to the 6×6 case, the characteristic fields are genuinely nonlinear only in the domain Ω 1 = {(d, h) ∈ R 2 ; d = 0}. Nevertheless in that case we can construct the solution by using multiple waves, following Wendroff [START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF] and Liu [START_REF] Liu | The Riemann problem for general 2 × 2 conservation laws[END_REF].

System 2 being a particular case of 1, we find Lax' solutions of 1 which are also weak solutions of 2, but they are different from the "Liu's solutions". Moreover, we shall prove that the electromagnetic energy is dissipated by both solutions, so that there exists (at least) two selfsimilar entropy solutions of the Riemann problem for system 2.

A first study of the Riemann problem can be found in [START_REF] De | High-order scheme for a nonlinear Maxwell system modelling Kerr effect[END_REF]: the problem is solved for a reduced 4 × 4 system and it is assumed that D • ω is identically zero. Related numerical schemes are constructed in 1D and 2D transverse electric configurations.

In [START_REF] Aregba-Driollet | Kerr-Debye relaxation shock profiles for Kerr equations[END_REF] we studied Kerr shocks and related shock profiles provided by the Kerr-Debye model, which is a hyperbolic quasilinear relaxation approximation of Kerr model. We studied Lax and Liu's admissibility criteria and proved that only Lax shocks give rise to Kerr-Debye shock profiles.

In [START_REF] Kanso | Sur le modèle de Kerr-Debye pour la propagation des ondes électromagnétiques[END_REF], Godunov's scheme, which requires the solution of the Riemann problem, was implemented for a a two-dimensionnal transverse electric configuration. Actually this case reduces to the one of the p-system 2. Liu's solution was implemented. The results were found to coincide with the ones obtained by a Kerr-Debye relaxation scheme.

The plan of the paper is the following. In section 2 we determine the simple waves and the wave functions. In section 3 we construct the solution of the Riemann problem. Section 4 is devoted to the 2 × 2 case.

2. Wave functions. If u -= u + are connected by a k-Lax shock or a k-rarefaction wave or a k-contact discontinuity, u -and u + are said to be connected by a k-wave.

A plane discontinuity σ, u + , u -is a weak solution u of 1 such that u(x, t) = u -if x • ω < σt, u(x, t) = u + else.
All the plane discontinuities have been studied in [START_REF] Aregba-Driollet | Kerr-Debye relaxation shock profiles for Kerr equations[END_REF]. The centred rarefaction waves are computed in [START_REF] Aregba-Driollet | Godunov schemes for Kerr equations[END_REF]. The results are the following.

Proposition 2. Contact discontinuities. Stationary contact discontinuities are characterized by ω × [H] = 0, ω × [E] = 0. ( 6
)
The divergence free ones are constant. A discontinuity σ, u + , u -is a contact discontinuity associated to λ 1 or λ 6 if and only if

|E + | = |E -| σ 2 = c 2 (1 + ǫ r |E + | 2 ) -1 = c 2 (1 + ǫ r |E -| 2 ) -1 and ω • [D] = 0 [H] = σω × [D].

Moreover the only discontinuities satisfying Rankine-Hugoniot conditions and such that |E

-| = |E + | are the above contact discontinuities.
The shocks and rarefactions are related to the second and sixth characteristic fields. We recall that a discontinuity σ, u -, u + is a Lax' k-shock if

λ k (u + ) < σ < λ k+1 (u + ), λ k-1 (u -) < σ < λ k (u -).
Those inequalities are entropy conditions and also ensure that one can construct the solution of the Riemann problem as a superposition of simple waves. In our case, defining

f (d, d 0 ) = c 2 d 1 + ǫ r p 2 d 2 0 + d 2 , (d, d 0 ) ∈ R 2 ,
we express the shocks by using the function S defined as

S(d 1 , d 2 , d 0 ) = ((f (d 2 , d 0 ) -f (d 1 , d 0 ))(d 2 -d 1 )) 1 2 .
The rarefaction waves are obtained via the integral curves of the eigenvectors of the system. ζ being a unitary vector orthogonal to ω, we use the function R ζ defined by

R ζ (d 1 , d 2 , d 0 ) = d2 d1 λ(d 0 ω + sζ)ds, d 1 ≤ d 2 , d 0 ∈ R.
Finally let φ ζ be the function defined for d 1 ≥ 0, d 2 ≥ 0 and d 0 ∈ R by

φ ζ (d 1 , d 2 , d 0 ) = S(d 2 , d 1 , d 0 ) if d 1 ≥ d 2 , -R ζ (d 1 , d 2 , d 0 ) if d 1 < d 2 .
Proposition 3. φ ζ is a decreasing C 1 function with respect to d 2 and for all d ≥ 0, d 0 ∈ R:

φ ζ (d, 0, d 0 ) = cd 1 + ǫ r p 2 ( d 2 0 + d 2 )
, lim

d2→+∞ φ ζ (d, d 2 , d 0 ) = -∞.
The 2 and 5 waves are characterized as follows, see [START_REF] Aregba-Driollet | Godunov schemes for Kerr equations[END_REF], [START_REF] Aregba-Driollet | Kerr-Debye relaxation shock profiles for Kerr equations[END_REF] for the proof:

D ± = d 0 ω + d ± ζ, H + = H -+ φ ζ (d -, d + , d 0 )ω × ζ .
u -and u + are connected by a 5-wave if there exist two distinct nonnegative real numbers d -, d + such that

D ± = d 0 ω + d ± ζ, H + = H -+ φ ζ (d + , d -, d 0 )ω × ζ .
3. Solution of the full wave Riemann problem. Suppose that u ± = (D ± , H ± ) and ω ∈ R 3 , |ω| = 1, are given. We look for intermediate states u 1 , u * , u * * , u 2 such that:

• u -and u 1 are connected by a 1-contact discontinuity,

• u 1 and u * are connected by a 2-wave,

• u * and u * * are connected by a stationary contact discontinuity,

• u * * and u 2 are connected by a 5-wave,

• u 2 and u + are connected by a 6-contact discontinuity. In the following we shall denote d ± 0 = D ± • ω. Suppose that a solution exists. For the 1 and 6 contact discontinuities, the following conditions have to be fulfilled:

D 1 • ω = D -• ω = d - 0 , |D 1 | = |D -|, D 2 • ω = D + • ω = d + 0 , |D 2 | = |D + |, (7)
H 1 -H -= σ -ω × (D 1 -D -), H + -H 2 = σ + ω × (D + -D 2 ). ( 8
)
where

σ ± = ±c (1 + ǫ r |E ± |) -1 2 .
For the 2 and 5 waves we know that D In the second case we have e * = e * * , which also reads as

f (d * , d - 0 ) = f (d * * , d + 0 ). ( 9
)
First case:

ω × D * = ω × D * * = 0. In that case, D * = d - 0 ω, D * * = d + 0 ω.
σ 2 = - f (d 1 , d - 0 ) -f (0, d - 0 ) d 1 = σ -.
In the same way, u * * and u 2 are the left and right states of a 5-shock propagating with speed σ + . Consequently the contact discontinuities merge with the shocks. Let us denote

V = ω × (H + -H --ω × (σ + D + -σ -D -)) .
(10) Conditions 6, 8 and Rankine-Hugoniot conditions on the shocks imply that V = 0 and

H * = H --ω × σ -D -, H * * = H + -ω × σ -D + . (11
D 1 = d - 0 ω + d 1 ζ, D 2 = d + 0 ω + d 2 ζ, (12)
D * = d - 0 ω + d * ζ, D * * = d + 0 ω + d * * ζ, ( 13
) with d 1 ≥ 0, d * > 0, d * * > 0, d 2 ≥ 0. Let us denote d = D • ζ, h = H • (ω × ζ). By 7-8:    d 1 = |ω × (ω × D -)|, d 2 = |ω × (ω × D + )|, h 1 = h -+ σ -(d 1 -d -) h 2 = h + + σ + (d 2 -d + ). (14)
As u * , u * * ∈ Ω(ω), one can define the 2 and 5-wave curves:

H * -H 1 = φ ζ (d 1 , d * , d - 0 )ω × ζ , H 2 -H * * = φ ζ (d 2 , d * * , d + 0 )ω × ζ . (15)
By 6, h * = h * * and

h * = h 1 + φ ζ (d 1 , d * , d - 0 ) = h 2 -φ ζ (d 2 , d * * , d + 0 ).
Therefore, using 9, we see that d * and d * * are solution of the two by two system

   f (d * , d - 0 ) = f (d * * , d + 0 ) , h 1 + φ ζ (d 1 , d * , d - 0 ) = h 2 -φ ζ (d 2 , d * * , d + 0 ). (16)
As φ ζ is decreasing and d * , d * * are positive:

φ ζ (d 1 , d * , d 0 ) + φ ζ (d 2 , d * * , d 0 ) < φ ζ (d 1 , 0, d - 0 ) + φ ζ (d 2 , 0, d + 0 ) = σ + d 2 -σ -d 1 .
This inequality is useful to determine ζ. As a matter of fact, using 8, we have also

H * = H -+ σ -ω × (D 1 -D -) + φ ζ (d 1 , d * , d - 0 )ω × ζ, H * * = H + -σ + ω × (D + -D 2 ) -φ ζ (d 2 , d * * , d + 0 )ω × ζ . Again by 6: V = σ + d 2 -σ -d 1 -φ ζ (d 1 , d * , d - 0 ) -φ ζ (d 2 , d * * , d + 0 ) ζ.
Therefore V = 0 and

ζ = V |V | . (17)
Those results can be summarized as follows.

Proposition 5. Consider u -, u + such that the Riemann problem for system 1 has a solution which is a superposition of simple waves. Let V be the vector defined in 10. Only the following two cases occur: 1) V = 0, u -and u * are connected by a 2-Lax shock propagating with velocity σ -, u + and u * * are connected by a 5-Lax shock propagating with velocity σ + , D * = (D -• ω)ω, D * * = (D + • ω)ω, H * and H * * are given by 11.

2) V = 0, ζ is defined by 17, u 1 and u 2 are determined by conditions 12, 14 and u * , u * * are determined by [START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF]15 and the solution of system 16.

Sufficient conditions are as follows.

Theorem 3.1. Let u -, u + be a Riemann data for system 1 in the direction ω. There exists η > 0 such that if |(D --D + ) • ω| < η then the Riemann problem has a unique solution in the class of the functions which are superpositions of simple waves. Let V be the vector defined in 10.

If V = 0, then the solution is the superposition of a 2-Lax' shock, a stationary contact discontinuity and a 5-Lax' shock.

If V = 0, then the solution is the superposition of a 1-contact discontinuity, a 2wave (Lax' shock or rarefaction), a stationary contact discontinuity, a 5-wave (Lax' shock or rarefaction) and a 6-contact discontinuity.

Moreover we can construct the solution in every case.

To prove this theorem, we remark that if d + 0 = d - 0 , then system 16 can be solved. For the general case we use the implicit function theorem. We refer to [START_REF] Aregba-Driollet | Godunov schemes for Kerr equations[END_REF] for the detailed proof.

4. Non-uniqueness of selfsimilar entropy solutions. Here we consider the Riemann problem 1-3 with ω = (1, 0, 0) T , D ± = (0, D 2,± , 0), H ± = (0, 0, H 3,± ). As (D + -D -) • ω = (H + -H -) • ω = 0, the stationary contact discontinuity is trivial. We have V = (0, v, 0) and if v = 0 then ζ = (0, sgn(v), 0). Hence Lax' solution is of the form (0, D 2 , 0), (0, 0, H 3 ) and (D 2 , H 3 ) is a weak solution of the Riemann problem for the p-system 2.

Reciprocally, weak solutions (d, h) of 2 give solutions (0, d, 0), (0, 0, h) of 1. Following [START_REF] Liu | The Riemann problem for general 2 × 2 conservation laws[END_REF] (see also [START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF]) one can compute a weak solution of the Riemann problem for the p-system by using multiple waves. In that case the shocks satisfy Liu's (E) condition, which generalize Lax' conditions when p has inflexion points. This solution differs from the first one.

As an example, we computed both solutions in a particular case. The results for a fixed time t are depicted in Figure 1. The Lax' solution of 1-3 consists of a 1-contact discontinuity, a 2-shock, a 5-shock and a 6-contact discontinuity. Liu's solution consists of a shock propagating to the left and a shock propagating to the right. Observe that in that case the sign of d can change through the shock, while Lax' shock condition applied to 1 imply that the sign of d = D 2 is constant.

To be more precise, we denote -µ 1 = µ 2 = µ At this point, we have two distinct solutions of the Riemann problem for 2, or equivalently for 1. Moreover, we can prove that both dissipate the physical entropy, namely the electromagnetic energy. For those solutions, this entropy is just the wellknown entropy of the p-system:

Proposition 4 .

 4 If u -= u + are connected by a 2 or a 5 wave, then D -= D + and D -• ω = D + • ω. Moreover ω × (ω × D -) and ω × (ω × D + ) are colinear. Reciprocally, let us consider u -and u + such that D -= D + and D -•ω = D + •ω = d 0 . If ω × D + = 0, we set D = D + . Else, ω × D -= 0 and we set D = D -. We define ζ by ζ = -ω × (ω × D) |ω × (ω × D)| . u -and u + are connected by a 2-wave if there exist two distinct nonnegative real numbers d -, d + such that

1 ,

 1 D * , ω are coplanar and D 2 , D * * , ω are coplanar. Moreover [D] • ω = 0. There exist unitary vectors ζ 1 , ζ 2 , orthogonal to ω such thatD 1 = d - 0 ω + d 1 ζ 1 , D * = d - 0 ω + d * ζ 1 and D 2 = d + 0 ω + d 2 ζ 2 , D * * = d + 0 ω + d * * ζ 2 and d 1 , d * , d * * , d 2 are non negative.The stationary contact discontinuity is defined by conditions 6. One hasE * = e * 0 ω + e * ζ 1 , E * * = e * * 0 ω + e * * ζ 2, wheree * = d * ǫ 0 (1 + ǫ r p 2 (|D * |)) , e * * = d * * ǫ 0 (1 + ǫ r p 2 (|D * * |)). Therefore e * ω × ζ 1 = e * * ω × ζ 2 . Hence either e * = e * * = 0 or those quantities are both positive and ζ 1 = ζ 2 . The first case occurs if and only if ω × D * = ω ×D * * = 0.

  u 1 and u * are the left and right states of a 2-shock propagating with speed

  ) If D -× ω = 0 then u -= u * . Else u -and u * are connected by a 2-Lax shock. In the same way, if D + × ω = 0 then u + = u * * , else u + and u * * are connected by a 5-Lax shock. Second case: D * × ω = 0 and D * * × ω = 0. In this case, ζ 1 = ζ 2 = ζ and

Figure 1 .

 1 Figure 1. Two solutions of the p-system for the same initial data. Left: d, right: h.

  η(d, h) = E(d) + 1 2 µ 0 h 2 , E(d) = 1 2 ǫ 0 (e 2 + 3ǫ r 2 e 4), e = p(d)

with entropy flux Q(d, h) = eh. Contact discontinuities and rarefactions conserve the entropy. For the shocks a straightforward calculation leads to the following, see [START_REF] Aregba-Driollet | Godunov schemes for Kerr equations[END_REF] for details. Proposition 6. Entropy dissipation for the 2 × 2 system

• Lax' and Liu's shocks satisfy the entropy dissipation property: We have solved the Riemann problem for the full wave Kerr system and we have proved the non-uniqueness of selfsimilar entropy solutions. Which is the physical solution? Liu's solution is more dissipative than the other one. 1D and 2D numerical experiments with a physically relevant relaxation model, the Kerr-Debye system, lead to Liu's solutions see [START_REF] Aregba-Driollet | Numerical approximation of Kerr-Debye equations[END_REF], [START_REF] Kanso | Sur le modèle de Kerr-Debye pour la propagation des ondes électromagnétiques[END_REF], but this is possibly due to numerical viscosity: it is wellknown that contact discontinuities are not easy to catch numerically. On another hand the full wave model should be the more realistic one. Further investigations are in progress.


Hynek Kova

¸oise Truc
email: francoise.truc@ujf-grenoble.fr

DICATAM Hynek Kovařík
email: hynek.kovarik@ing.unibs.it

SCHR ÖDINGER OPERATORS ON A HALF-LINE WITH INVERSE SQUARE POTENTIALS

We consider Schrödinger operators Hα given by equation (1.1) below. We study the asymptotic behavior of the spectral density E(Hα, λ) for λ → 0 and the L 1 → L ∞ dispersive estimates associated to the evolution operator e -itHα . In particular we prove that for positive values of α, the spectral density E(Hα, λ) tends to zero as λ → 0 with higher speed compared to the spectral density of Schrödinger operators with a short-range potential V . We then show how the long time behavior of e -itHα depends on α. More precisely we show that the decay rate of e -itHα for t → ∞ can be made arbitrarily large provided we choose α large enough and consider a suitable operator norm.

Introduction

This paper is concerned with Schrödinger operators

H α = - d 2 dx 2 + α x 2 , α ≥ - 1 4 , (1.1)
in L 2 (R + ) with Dirichlet condition at x = 0. In particular, we are interested in the dependence of various spectral properties of H α on the parameter α. Note that potentials of the type α/x 2 have a special role, since the resulting operator H α is scaling invariant. Moreover, it is known that the potentials which satisfy V (x) ∼ x -2 as x → ∞ represent a borderline case for certain important spectral inequalities such as dispersive or Strichartz estimates, see [GVV].

It is therefore not surprising that Schrödinger operators with inverse square potentials have recently attracted certain attention; we might mention for example the heat kernel bounds obtained in [START_REF] Milman | Heat kernel bounds and desingularizing weights[END_REF][START_REF] Milman | Global heat kernel bounds via desingularizing weights[END_REF], or Strichartz estimates in dimension three studied in [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF][START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with potentials of critical decay[END_REF]. Inverse square potentials appear naturally also in connection with two-dimensional Schrödinger operators with Aharonov-Bohm-type magnetic field, see [FFFP, GK].

Main objects of our interest here are the spectral density

E(H α , λ) = 1 2πi lim ε→0+ (H α -λ -iε) -1 -(H α -λ + iε) -1 , λ > 0 (1.2)
of H α , and the unitary group e -itHα . In particular, we are going to study the asymptotic behavior of E(H α , λ) for λ → 0 and the L 1 → L ∞ dispersive estimates associated to the evolution operator e -itHα . It is very well-known that the asymptotic behavior of E(H α , λ) for small λ is closely related to the asymptotic behavior of e -itHα for large t.

There is a huge amount of literature on this subject, see e.g. [START_REF] Erdogan | A weighted dispersive estimate for Schrödinger operators in dimension two[END_REF][START_REF] Goldberg | Dispersive estimates for Schrödinger operators in dimensions one and three[END_REF][START_REF] Jensen | Spectral properties of Schrödinger operators and time-decay of the wave functions[END_REF][START_REF] Murata | Asymptotic expansions in time for solutions of Schrödinger-type equations[END_REF][START_REF] Schlag | Dispersive estimates for Schrödinger operators in dimension two[END_REF][START_REF] Schlag | Dispersive estimates for Schrödinger operators: a survey[END_REF][START_REF] Weder | L p -L p ′ Estimates for the Schrödinger Equation on the Line and Inverse Scattering for the Nonlinear Schrödinger Equation with a Potential[END_REF][START_REF] Weder | The L p -L p estimate for the Schrödinger equation on the half-line[END_REF] and references therein. We are not going to discuss this connection any further since it will not be used in our proofs.

For general one-dimensional Schrödinger operators of the type H V = -d 2 dx 2 + V the behavior of both E(H V , λ) and e -itH V is known provided the potential V decays fast enough at infinity. In particular, if zero is a regular point of H V , (which is the generic case), then

E(H V , λ) ∼ λ 1 2 , λ → 0, (1.3)
in a suitable operator topology, see [START_REF] Goldberg | Transport in the one-dimensional Schrödinger equation[END_REF][START_REF] Murata | Asymptotic expansions in time for solutions of Schrödinger-type equations[END_REF][START_REF] Schlag | Dispersive estimates for Schrödinger operators: a survey[END_REF][START_REF] Weder | The L p -L p estimate for the Schrödinger equation on the half-line[END_REF]. Accordingly, for such short range potentials, under certain regularity conditions, Murata [Mu] proved

w -1 e -itH V w -1 L 2 (R)→L 2 (R) ≤ C t -3 2 ∀ t > 2, (1.4)
where w is a weight function with a sufficient growth at infinity. The corresponding L 1 → L ∞ was established by Schlag

ρ -1 e -itH V ρ -1 L 1 (R)→L ∞ (R) ≤ C t -3 2 ∀ t > 2, (1.5)
with ρ(x) = (1 + |x|), see [START_REF] Schlag | Dispersive estimates for Schrödinger operators: a survey[END_REF]. It is important to mention that the decay conditions on V , under which all the above results were obtained, imply that

V (x) = o(x -2 ) as |x| → ∞.
The goal of the present note is to show that if V is of type α x -2 with α > 0, then the asymptotic relation (1.3) is no longer valid and has to replaced by a new one, and, on the other hand, the estimates (1.4) and (1.5) can be improved. In particular E(H α , λ) decays faster to zero than in (1.3), see Theorem 2.1. Accordingly the decay in the dispersive estimate (1.5) can be improved provided the weight function ρ grows fast enough at infinity, see Theorem 2.4. Although our results regard a family of Schrödinger operators with explicit potentials, it can be expected that similar results should hold also if H α is perturbed by a sufficiently short-range perturbation.

It should be finally mentioned that our main results, i.e. Theorems 2.1 and 2.4, fail in the case of Schrödinger operators on the whole line due to the presence of the zero resonance.

Main results

2.1. Notation. We set ρ(x) = 1 + x on R + . For any s ∈ R we denote

L 2 s (R + ) = {u : ρ s u L 2 (R + ) < ∞}, u 0,s := ρ s u L 2 (R + ) .
Let B(s, s ′ ) be the space of bounded linear operators from L 2 s (R + ) to L 2 s ′ (R + ) and let • B(s,s ′ ) denote the corresponding operator norm. Finally, we put ν = 1/4 + α.

(2.1)

We have

Theorem 2.1. Let α > -1/4. Then for any ε > 0 and any

s ≥ ν + 1 + ε it holds E(H α , λ) = E 0 λ ν + O(λ ν+ε ) λ → 0+ (2.2) in B(s, -s)
, where E 0 is the integral operator in L 2 (R + ) with the kernel

E 0 (x, y) = (x y) ν+ 1 2 2 ν Γ 2 (ν + 1) .
Remark 2.2. Equation (2.2) shows that for positive values of α the density E(H α , λ) is of lesser order than in the case of a short-range potential, see equation (1.3).

Remark 2.3. For a throughout discussion of threshold expansion of resolvents of onedimensional operators with short-range potentials we refer to [JN]. Asymptotic behaviour of Schrödinger groups generated by operators with inverse square decay on conical manifolds was studied in [Wa] in the setting of weighted L 2 -spaces.

Theorem 2.4. Let α ≥ -1/4. Then for any s ∈ [0, ν + 1/2] there exists a constant C(α, s) such that

ρ -s e -itHα ρ -s L 1 (R + )→L ∞ (R + ) ≤ C(α, s) t -1 2 -s ∀ t > 0.
(2.3)

Remark 2.5. For -1/4 < α ≤ 0 the dispersive estimate (2.3) can be derived from [START_REF] Fanelli | Time decay of scaling critical electromagnetic Schrödinger flows[END_REF]Thm.1.11] by considering the restriction of inequality [START_REF] Fanelli | Time decay of scaling critical electromagnetic Schrödinger flows[END_REF]Eq.(1.29)] to radial functions. On the other hand, the result for α > 0, namely the faster decay of e -itHα in t is new. The maximal decay rate t -1-ν , achieved by the choice s = ν + 1/2, should be compared with the t -3 2 decay rate in the estimate (1.5).

Remark 2.6. Note also that in the border-line case α = -1/4, which means ν = 0, Theorem 2.4 with the choice s = 1/2 gives the decay rate t -1 , which is the decay rate of the free evolution e it∆ in dimension two. This is not surprising since the operator

-d 2 dx 2 -1 4x 2 in L 2 (R + )
with Dirichlet boundary condition at zero is unitarily equivalent, by means of the unitary mapping f

(x) → √ x f (x), to the Laplacian -∆ in L 2 (R 2 )
restricted to radial functions.

An immediate consequence of Theorem 2.4 is the following Corollary 2.7. Let α ≥ -1/4. Then for any s ∈ [0, ν + 1/2] and any β > s + 1/2 there exists a constant C 2 , depending only on α, β and s, such that

ρ -β e -itHα ρ -β L 2 (R + )→L 2 (R + ) ≤ C 2 t -1 2 -s ∀ t > 0.
(2.4)

Proof. Let u ∈ L 2 β (R + ) and let f = ρ s u.
Then by the Cauchy-Schwarz inequality f ∈ L 1 (R + ) and

f 2 L 1 (R + ) = ∞ 0 ρ s ρ -β ρ β u dx 2 ≤ C 1 ρ β u 2 L 2 (R + ) , (2.5)
where we have used the fact that β > s + 1/2. Hence from Theorem 2.4 and (2.5) we obtain

ρ -β e -itHα u 2 L 2 (R + ) = ∞ 0 ρ(x) 2s-2β | ρ -s e -itHα ρ -s f | 2 dx ≤ C 2 (α, s) t -1-2s f 2 L 1 (R + ) ≤ C 2 (α, s) C 1 t -1-2s ρ β u 2 L 2 (R + )
. This proves (2.4). Inequality (2.4) should be compared with the estimate (1.4) valid for short-range potentials.

3. Proofs 3.1. Proof of Theorem 2.1. For simplicity we shall drop the index α in the sequel and write E(λ) instead of E(α, λ). We also use the notation

R(λ, x, y) = lim ε→0+ (H α -λ -iε) -1 (x, y).
We first study the solutions u ∈ L 2 (R + ) of the generalized eigenvalue equation

-u ′′ - α x 2 = λu. (3.1) After setting u(x) = √ x ψ( √ λ x), equation (3.1) writes z 2 ψ ′′ -x + z ψ ′ + (z 2 -ν 2 ) ψ = 0 , (3.2) with z = √ λ x.
The latter is a Bessel equation of the first kind, see [AS,Sec.9.1]. We now find two solutions u 1 , u 2 of (3.1) which satisfy u 1 (0) = 0 and

u 2 ∈ L 2 (1, ∞) for Im λ > 0. Since | J ν (z) + i Y ν (z) | ∼ 2 zπ | e iz | , |z| → ∞, Im z > 0,
see [AS,Eqs.9.1.3,9.2.3], and

J ν (z) = (z/2) ν Γ(ν + 1) + o(z ν ) z → 0. (3.3)
by [AS, Eq.9.1.7], the sought solutions u 1 and u 2 take the form

u 1 (x) = √ x J ν ( √ λ x) (3.4) u 2 (x) = √ x (J ν ( √ λ x) + i Y ν ( √ λ x)) . (3.5)
Hence by the theory of Sturm-Liouville problems we obtain the resolvent kernel

R(λ, x, y) = iπ 2 √ xy J ν ( √ λ x) (J ν ( √ λ y) + i Y ν ( √ λ y)) (x ≤ y) (3.6) R(λ, x, y) = iπ 2 √ xy J ν ( √ λ y) (J ν ( √ λ x) + i Y ν ( √ λ x)) (x ≥ y) (3.7)
The Stone formula (1.2) then implies that

E(λ, x, y) = 1 π Im R(λ, x, y) = 1 2 √ xy J ν ( √ λ x) J ν ( √ λ y) . (3.8)
From (3.3) we now easily verify that lim

λ→0+ λ -ν E(λ, x, y) = (x y) ν+ 1 2 2 ν Γ 2 (ν + 1) = E 0 (x, y). (3.9)
Let us define the rest term E 1 (λ) as the integral operator in L 2 (R + ) with the kernel given by

E 1 (λ, x, y) = E(λ, x, y) -E 0 (x, y) λ ν . (3.10)
To prove Theorem 2.1 we need the following Lemma 3.1. For any ε > 0 and any s > ν + 1 + ε we have

E 1 (λ) B(s,-s) = O(λ ν+ε ) λ → 0 + . (3.11)
Proof.-

We will use the fact that

E 1 (λ) B(s,-s) = ρ -s E 1 (λ)ρ -s L 2 (R + )→L 2 (R + ) .
(3.12)

From (3.10) we get that

ρ -s λ -ν E(λ) ρ -s = ρ -s E 0 ρ -s + ρ -s λ -ν E 1 (λ) ρ -s . (3.13)
Note that the operator E 0 is Hilbert-Schmidt in B(s, -s). This follows from the identity (3.12) applied to E 0 . Hence by applying the Taylor formula to the operator ρ -s λ -ν E(λ) ρ -s at λ = 0 we find that the claim of the Lemma will follow if we show that

ρ -s ∂ λ (λ -ν E(λ)) ρ -s HS(R + ) = ρ -s ∂ λ (λ -ν E 1 (λ)) ρ -s HS(R + ) = O(λ -1+ε ) λ → 0, (3.14)
where • HS(R + ) denotes the Hilbert-Schmidt norm in L 2 (R + ). Using the recurrence relations for the derivatives of J ν :

J ′ ν (z) = -J ν+1 (z) + ν z J ν (z) J ′ ν (z) = J ν+1 (z) - ν z J ν (z),
see [AS,Eq.9.1.27], we get from (3.8)

∂ λ (λ -ν E(λ, x, y)) = - λ -ν-1/2 4 √ x y (x J ν+1 ( √ λ x) J ν ( √ λ y) + y J ν+1 ( √ λ y) J ν ( √ λ x) . (3.15)
Hence by the Cauchy-Schwarz inequality

ρ -s ∂ λ (λ -ν E(λ)) ρ -s 2 HS(R + ) = = ∞ 0 ∞ 0 |∂ λ (λ -ν E(λ, x, y))| 2 ρ(x) -2s ρ(y) -2s dxdy ≤ C λ -1-2ν I(λ) J (λ) (3.16)
where

I(λ) = ∞ 0 x 3 J 2 ν+1 ( √ λ x) (1 + x) -2s dx J (λ) = ∞ 0 y J 2 ν ( √ λ y) (1 + y) -2s dy
and C is a constant independent of λ. To estimate the last two integrals we will need a point-wise estimate on the Bessel function J ν . From the integral representation

J ν (z) = 2 ( z 2 ) ν √ π Γ(ν + 1 2 ) 1 0 (1 -t 2 ) ν-1 2 cos(zt) dt,
see [AS,Eq.9.1.20], it follows that |J ν (z)| ≤ C ν z ν for all z > 0 and ν > 0.

On the other hand, by [AS,Eq.9.1.20] we have |J ν (z)| ≤ 1 for all z > 0 and ν ≥ 0. A combination of these two upper bounds then implies that for any -1/2 ≤ µ ≤ ν there exists a constant C(µ, ν) such that

|J ν (z)| ≤ C(µ, ν) z µ ∀ z > 0, ∀ µ ∈ - 1 2 , ν .
(3.17) Using (3.17) with µ = ν -1 + 2ε in I(λ) and with µ = ν in J (λ) together with the fact that s ≥ ν + 1 + ε, we find

I(λ) = O(λ ν-1+2ε ), J (λ) = O(λ ν ) λ → 0.
In view of (3.16) this implies (3.14) and therefore completes the proof.

Theorem 2.1 now follows from (3.10) and (3.11).

3.2. Proof of Theorem 2.4. We will prove Theorem 2.4 by estimating the integral kernel of the operator e -itHα . To provide a formula for the integral kernel, we will follow [START_REF] Kovařík | Heat kernels of two-dimensional magnetic Schrödinger and Pauli operators[END_REF]Sec.5], where the formula for the integral of the heat semi-group e -tHα was established, see also [GK]. Equation (3.8) in combination with the Weyl-Titchmarsh-Kodaira Theorem, cf. [START_REF] Dunford | Linear Operators, Part II[END_REF]Chap.13], shows that the operator H α is unitarily equivalent to a multiplication operator, namely we have

U ν H α U -1 ν f (p) = p f (p), f ∈ U ν (D(H α )), (3.18)
where D(H α ) denotes the operator domain of H α and the mappings U

ν , U -1 ν : L 2 (R + ) → L 2 (R + ) are given by (U ν g)(p) = ∞ 0 g(x) √ x J ν (x √ p) dx (U -1 ν f )(x) = 1 2 ∞ 0 f (p) √ x J ν (x √ p) dp (3.19)
The mapping U ν and U -1 ν define unitary operators on L 2 (R + ). Let g ∈ C ∞ 0 (R + ). By [START_REF] Teschl | Mathematical Methods in Quantum Mechanics With Applications to Schrödinger Operators[END_REF]Thm.3.1] e -itHα g = lim ε→0+ e -(ε+it) Hα g.

(3.20)

In view of (3.18) we thus get where we have used [START_REF] Erdelyi | Tables of integral transforms[END_REF]Eq.4.14(39)] to calculate the integral with respect to p. Moreover, from [AS,Eq.9.6.18] it follows that the function is bounded on every compact interval uniformly with respect to ε > 0. Since the support of g is compact, we can use the dominated theorem and interchange the limit and integration in (3.21). Taking the limit ε → 0 and using the identity I ν (iz) = e -iνπ/2 J ν (z), see [AS,Eq.9.6.3], we obtain The last equations now imply that

lim ε→0+ e -(ε+it)Hα g (r) = lim ε→0+ U -1 ν e -(ε+it) p U ν g (x) = lim ε→0+ 1 2 ∞ 0 √ xy ∞ 0 e -(ε+it) p J ν (x √
ρ -s e -itHα ρ -s f L ∞ (R + ) ≤ C(α, s) t -1 2 -s f L 1 (R + )
for all f ∈ L 1 (R + ). This proves inequality (2.3).

  the upper bound (3.17) with µ = s -1/2 ∈ [-1/2, ν] and z = xy 2t .

Acknowledgements

H.K. would like to thank the Institut Fourier in Grenoble for the warm hospitality extended to him during his stay. The work of H.K. has been partially supported by the MIUR-PRIN'2010-11 grant for the project "Calcolo della variation". Both authors thank the referee whose remarks and comments helped them improve the original version of the text.


Maya Briani
email: m.briani@iac.cnr.it

Denise Aregba-Driollet

Roberto Natalini
email: r.natalini@iac.cnr.it

How to improve the decay of the numerical error for large times: the case of dissipative BGK systems

Keywords: 2000 Mathematics Subject Classification. Primary: 65M12; Secondary: 35L65 Finite differences methods, dissipative hyperbolic problems, BGK systems, asymptotic behavior, asymptotic high order schemes

We introduce new finite differences schemes to approximate one dimensional dissipative semilinear hyperbolic systems with a BGK structure. Using accurate analytical time-decay properties of the local truncation error, it is possible to design schemes based on standard upwinding schemes, which are increasingly accurate for large times when computing small perturbations of constants asymptotic states.

1. Introduction. Consider the following class of one dimensional BGK systems:

∂ t f i + λ i ∂ x f i = M i (u) -f i , i = 1, ..., m. (1)
Here f i ∈ R k , u := m i=1 f i , x ∈ R and t > 0, and the functions M i = M i (u) ∈ R k are smooth functions of u such that:

m i=1 M i (u) = u.
To obtain the time decay rates of these solutions, we need to rewrite the problem in more suitable coordinates. Following [START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF], we rewrite the BGK system in its conservative-dissipative form for the new unknowns

Z = (u, Z) T .
It is proved in [START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF] that, under some dissipativity conditions and for initial data which are small and smooth in some suitable norms, the time decay of the global solutions, for large times and in the L ∞ -norm, is given by

∂ l x u = O(t -1/2-l/2 ), ∂ l x Z = O(t -1-l/2
), and similar estimates are available for their time derivatives. Notice that the improved estimate for Z can only be obtained in these new coordinates.

The aim of this paper is to give a brief overview of the way it is possible to take advantage of these precise decay estimates to build up more accurate numerical schemes. Actually, we can see that, for standard numerical schemes, like for instance the upwind scheme with the source term approximated pointwise by the standard Euler scheme, the truncation error has the following decay as t → +∞:

T u (x, t) = O(∆x t -3/2 ) + O(∆t t -3/2 ), T Z (x, t) = O(∆x t -3/2 ) + O(∆t t -3/2 ).
(

) 2
It can be seen numerically that the corresponding absolute error, for a fixed space step, decays as

e u (t) = O(t -1/2 ),e z (t) = O(t -1 ),
which implies that the relative error is essentially constant in time.

Here, our main goal is to improve the decay estimates on the truncation order to achieve an effective decay in time of the relative error, both in u and Z. To obtain this result, we perform a detailed analysis of the behavior of the truncation error for a general class of schemes, called "Time Asymptotically High Order" (TAHO) schemes, which generalize those introduced in [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2X2 dissipative hyperbolic systems[END_REF]. Thanks to this analysis, we are able to select some schemes such that the truncation order behaves as

T u (x, t) = O(∆x t -2 ), T Z (x, t)) = O(∆x t -2 ), (3)
for a fixed CFL ratio and such that the numerical error observed in the practical tests improves of t -1/2 on other schemes. The plan of the paper is the following. In Section 2, we introduce our analytical framework. The main schemes are derived in Section 3, where we show how to improve the time decay of their local truncation error. Section 4 presents some numerical tests which show the nice behavior of our new schemes in two test cases.

2. The analytical framework. Following [START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF], we rewrite system 1 in its conservativedissipative form. This means that we assume that there exists an invertible matrix

D = D 11 D 12 D 21 D 22 , (4)
such that, setting m 1 = k, m 2 = k(m -1), the new unknown

Z = Df = (u, Z) T ∈ R k × R m2 , (5)
solves the system

∂ t u + A 11 ∂ x u + A 12 ∂ x Z = 0, ∂ t Z + A 21 ∂ x u + A 22 ∂ x Z = Q(u) -Z, (6)
where A is symmetric and Q(u) is quadratic in u, i.e.: Q(0) = 0 and Q′ (0) = 0. Observe that, after this transformation, which a priori is not unique, the source term is zero in the first component and the second one is the sum of a quadratic term and of the dissipative term -Z. Moreover, when transforming system 1 in system 7, we can always assume that blocks D 11 and D 12 have the special form

D 11 = I k , D 12 = (I k I k • • • I k ) ∈ R k×m2 ,
and, setting Λ = diag(λ 1 I k , ..., λ m I k ), we have that

A = A 11 A 12 A 21 A 22 = DΛD -1 .
Therefore, we can rewrite our system in a more compact form:

∂ t Z + A∂ x Z = -Z + DM (u). (7)
To guarantee the existence of the matrix D in 5, we can assume that our system is strictly entropy dissipative in the sense of [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF] and verifies the Shizuta-Kawashima condition [START_REF] Shizuta | Systems of equations of hyperbolic-parabolic type with applications to the discrete Boltzmann equation[END_REF][START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF][START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF]. For instance, using Bouchut's Entropy dissipation condition [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation law[END_REF], it is possible to prove the existence of a matrix D in 5, with all the above properties.

2.1. The simplest example: the Jin-Xin 2 × 2 relaxation system. Consider the following system

∂ t u + ∂ x v = 0, ∂ t v + λ 2 ∂ x u = F (u) -v, (8)
with λ > 0. The unknowns u and v are scalar and the function F = F (u) is smooth, with F (0) = 0. This case is obtained from 1 for k = 1, m = 2, and λ 2 = -λ 1 = λ, by setting

u = f 1 + f 2 , v = λ(f 2 -f 1 ), F (u) = λ(M 2 -M 1 ).
Under the condition

λ > |F ′ (0)|, ( 9
)
the problem is dissipative, at least in a small neighborhood of the origin, in the sense of [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF] and the Shizuta-Kawashima condition is verified.

In this case the conservative-dissipative form is obtained by using

D = 1 1 -µa + µa - ,
where a = F ′ (0), µ = (λ 2 -a 2 ) -1/2 is real and positive, a ± = λ ± a > 0, from assumption 10.

2.2.

A 3 × 3 BGK example. Let us now compute the conservative-dissipative form for the following 3 × 3 BGK model,

   ∂ t f 1 -λ∂ x f 1 = M 1 (u) -f 1 , ∂ t f 2 = M 2 (u) -f 2 , ∂ t f 3 + λ∂ x f 3 = M 3 (u) -f 3 .
Let F = F (u) be a smooth scalar function such that F (0) = 0 and let γ be such that γ ′ (u) = |F ′ (u)|, with γ(0) = 0. We choose our three maxwellian functions as follows, for β ∈]0, 1[ and λ > 0

M 1 (u) = 1 2 γ(u) -F (u) λ + βu ,M 3 (u) = 1 2 γ(u) + F (u) λ + βu , M 2 (u) = u -M 1 (u) -M 3 (u) = (1 -β)u - γ(u) λ .
The functions M i , i = 1, 2, 3, are strictly increasing if for any u under consideration

λ > |F ′ (u)| 1 -β ,
and so the entropy dissipation condition [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation law[END_REF] is verified. Let a = F ′ (0) and α = |a| + βλ, following the results in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation law[END_REF][START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF][START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF], the matrix D for the transformation in the conservative-dissipative form 6 is given by

D =         1 1 1 α+a α-a λ(α-a) α(α+a) 0 -λ(α-a) α(α+a) -λ-α α -λ-α α + λ √ α(λ-α) -λ-α α         .
3. The numerical approximation. In this section we first introduce general finite difference approximations for system 1. Then, we compute the local truncation error of these schemes and we discuss its decays properties. The main result is given in Theorem 3.1, where a class of Time Asymptotically High Order (TAHO) schemes is fully characterized. First, we approximate the differential part following the direction of the characteristic velocities, so we study the methods for the system in diagonal form 1.

We denote by f = (f 1 , ..., f m ) the exact solution. Let ∆x the uniform meshlength and x j = j ∆x the spatial grid points for all j ∈ Z. The time levels t n , with t 0 = 0, are also spaced uniformly with mesh-length ∆t = t n+1 -t n for n ∈ N. We denote by ρ the CFL ratio ρ = ∆t/∆x, which is taken constant through all the paper.

We consider the Cauchy problem for system 1 possibly subjected to some stability conditions. The initial data f 0 is supposed to be smooth and approximated by its node values. The approximate solution (f 1 j,n , ..., f m j,n ) T , f i j,n ∈ R k , i = 1, ..., m, for j ∈ Z and n ∈ N, is given by

f i j,n+1 -f i j,n ∆t + λ i 2∆x f i j+1,n -f i j-1,n - q i 2∆x δ 2 x f i j,n = l=-1,0,1 B i l (u j+l,n ) -β i l f i j+l,n , (10)
with f i j,0 = f i 0 (x j ) and δ 2 x f j,n = (f j+1,n -2f j,n + f j+1,n ), for all i = 1, ..., m. The artificial diffusion terms q i are diagonal matrices in R k×k + . The source term approximation is defined, for l = -1, 0, 1, by the diagonal matrices β i l ∈ R k×k and by the vectors of functions B i l (•) ∈ R k . We assume the scheme 11 is consistent with system 1, i.e, for all i = 1, ..., m

β i -1 + β i 0 + β i 1 = I k + ∆xC i , B i -1 (u) + B i 0 (u) + B i 1 (u) = M i (u) + ∆xC i (u), where C i = diag(c i 1 , ..., c i k ) ∈ R k×k and C i (u)
are k functions to be defined.

Decay properties of the local truncation error.

In this section we focus on the local truncation error for the general scheme 11. By applying the time decay properties given in [START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF], we will show how it is possible to build up numerical schemes which are more accurate for large times. Set, for i = 1, ..., m,

C = diag(C i ), C = DCD -1 , C(u) = (C i (u)) T , γ i = (β i 1 -β i -1
). Scheme 11 is clearly consistent. Now, using the time decay estimates in [START_REF] Bianchini | Asymptotic behavior of smooth solutions for partially Dissipative hyperbolic systems with a convex entropy[END_REF] and similar estimates for their time derivatives, we obtain for a general approximation the following estimates for the local truncation error, as t → +∞:

T u (x, t) = O(∆x t -3/2 ) + O(∆t t -3/2 ),T z (x, t) = O(∆x t -3/2 ) + O(∆t t -3/2 ).
We would like to improve the decay property of this local truncation error to build up more accurate numerical schemes. The main idea is to chose the free parameters of the scheme to delete the terms that decay more slowly in the Taylor expansion of the local truncation error (see [START_REF] Aregba-Driollet | Time Asymptotic High Order schemes for dissipative BGK hyperbolic systems[END_REF] for details).

Let g i = diag(γ (i-1)k+1 , ..., γ ik ) for i = 1, ..., m and G = diag(g 1 , ..., g m ). If we make the following choice for the coefficients of the scheme 11,

C = - ρ 2 I km , C = CM (u) = - ρ 2 M (u), ( 11
)
g i = - 1 2 q i h -1 i - ρ 2 h -1 i P -(λ i I k -A 11 ) 2 (λ i I k -A 11 ) -1 h i ( 12
)
and

Γ ′ i (u) = g i M ′ i (u) + ρ 2   (h -1 i -M ′ i (u))A 11 + λ i M ′ i (u) -h -1 i m j=1 λ j M ′ j (u)   , (13)
both for i = 1, ..., m, then the scheme 11 is TAHO and the local truncation error reads

T u (x, t) = O ∆x t -2 + O ∆x 2 t -3/2 ,T z (x, t) = O ∆x t -2 + O ∆x 2 t -3/2 .
For the proof and further considerations in case A 11 = 0 we refer to [START_REF] Aregba-Driollet | Time Asymptotic High Order schemes for dissipative BGK hyperbolic systems[END_REF].

4. Numerical tests. In this Section we show how, for large time simulations, TAHO schemes give better numerical results than standard approximations for both examples considered in Sections 2.1 and 2.2. Specifically, we shall compare our TAHO scheme with two numerical approximations: i) a source pointwise approximation, denoted by STD and defined by 11 with Γ i = 0, γ i = 0, C i = 0 and C i = 0, for i = 1 =, ..., m; ii) a source upwinding approximation, denoted by ROE and defined by

f i j,n+1 -f i j,n ∆t + λ i 2∆x f i j+1,n -f i j-1,n - |λ i | 2∆x δ 2 x f i j,n = M i (u n j-1 ) + 2M i (u n j ) + M i (u n j+1 ) 4 + sgn(λ i ) 4 (M i (u n j-1 ) -M i (u n j+1 )) - f i j-1,n + 2f i j,n + f i j+1,n 4 - sgn(λ i ) 4 (f i j-1,n -f i j+1,n ). ( 14
)
To complete the definition of scheme 11 coupled with conditions 12-14 it is still necessary to choose some free parameters, such as for the 2 × 2 case B 1,2 0 (•) and β 1,2 0 . For both cases considered, such parameters can be defined by applying monotonicity conditions to the scheme. We refer to [START_REF] Aregba-Driollet | Time Asymptotic High Order schemes for dissipative BGK hyperbolic systems[END_REF] for more details.

For all tests, we focus our attention on the numerical error as a function of time: we plot the error e(t) = (u H -U h )(t) L ∞ as the time t = n∆t increases, where u H is the reference solution obtained by the ROE scheme 15, with ∆x = O(10 -4 ). Then, given different numerical approximations U h , we look for a constant C and γ which best fit the equality

e(t) = (u H -U h )(t) L ∞ = Ct -γ . ( 15
)
Given N data points (t i , e(t i )) i=1,N , we shall define γ and C as the solution of the following least squares problem, min

C,γ N i=1 | ln(e(t i )) -ln(Ct -γ )| 2 .
For all schemes, we fix the steps ratio ρ to verify all the CFL conditions. Since all schemes are of first order approximation, to emphasize the good behaviour of TAHO compared to the others schemes, we compute the numerical solutions U h by using a quite big grid step ∆x = O(10 -1 ).

All numerical results we present show that for standard approximations, such as STD and ROE, the absolute error e(t), for a fixed space step, decays as e u (t) = O(t -1/2 ), e z (t) = O(t -1 ), while for the TAHO scheme, it improves of t -1/2 on the previous schemes.

4.1.

Results for the Jin-Xin 2 × 2 system. We fix q = λ and we compare for the 2×2 case the TAHO scheme coupled with monotonicity assumptions, with ROE and STD scheme. We shall consider as initial datum the function

u 0 = χ [-1,1] -x 2 + 1 , z 0 = 1 λ F (u 0 (x)),
and we fix F (u) = a u -u 2 . The numerical results 1 show a better performance of the TAHO scheme; for both conservative and dissipative variable, the numerical solution obtained by TAHO fit better the benchmark curve. Again, the decay of the errors e u (t) and e z (t), Figure 1-(c)-(d), goes faster for the TAHO scheme, as confirmed by Table 1. There the decay parameter γ is numerically computed for all three schemes. The value obtained for the TAHO scheme improve of t -1/2 on the others. We stress on that the numerical solutions are computed with quite big step ∆x = 0.1.

4.2.

Results for the 3 × 3 system. As initial data, we take the smooth function u 0 defined by

u 0 (x) = χ [-1,1] exp 1 - 1 1 -x 2 .
Then we set f 0 (x) = M (u 0 (x)). We choose a = 1, λ = 2.1, β = (α -a)/λ = 0.1. The discretization parameters are ∆x = 0.1, ρ = 1 2λ , which satisfy all monotonicity requirements, see [START_REF] Aregba-Driollet | Time Asymptotic High Order schemes for dissipative BGK hyperbolic systems[END_REF].

The numerical results show as in the 2×2 case a better performance of the TAHO scheme. In Figure 2-(a)-(b), we plot the time evolution of the l ∞ errors e u (t) and e z (t). They show how for the TAHO scheme both errors decay as time increases more quickly than other. This result is also confirmed by Table 2, where the values of γ and C are computed. Looking at the different values of γ, it is clear that for the TAHO approximation the decay velocity of the absolute error improves of t -1/2 on the previous schemes. As expected by our asymptotic analysis, for the TAHO scheme the absolute errors e u,z (t) decay faster as the time increases. This result is confirmed in Table 1, where we compute the decay parameters γ of absolute errors previously plotted.

Theorem 3 . 1 (

 31 Local Truncation Error). Let ∆t/∆x = ρ be fixed and let H = diag(h 1 , . . . , h m ) be the block diagonal matrix given by H = D T D and set P = m i=1 λ 2 i h -1 i . Assume A 11 = 0 and that the following condition holds: the matrix (λ i I k -A 11 ) is invertible for i = 1, ..., m.

Figure 1 .

 1 Figure 1. The 2 × 2 Test, see subsection 4.1. (a)-(b) Zoom on the solutions u and z respectively obtained by the different schemes at final time T . The plot show that TAHO scheme gives better results than others with a quite big step ∆x = 0.1.(c)-(d) Time evolution of the l ∞ errors e u (t) and e z (t) defined in 16 for the different schemes. As expected by our asymptotic analysis, for the TAHO scheme the absolute errors e u,z (t) decay faster as the time increases. This result is confirmed in Table1, where we compute the decay parameters γ of absolute errors previously plotted.

Table 1 .

 1 The 2 × 2 Test, see subsection 4.1. Evaluation of constants γ and C for e u (t) = C u t -γu and e z (t) = C z t -γz defined in 16. For standard approximation STD and ROE, the absolute error decays as e u,z (t) ≈ O(t -1/2 ); while, for the TAHO scheme it improves of t -1/2 .

		scheme		C u		γ u		C z		γ z
		STD	0.013797 0.374708 0.010744 0.341554
		ROE 0.004874 0.333634 0.007850 0.439996
		TAHO 0.111380 1.151517 0.495480 1.451030
				U								Z
	0.048						ROE		0.046						ROE
							STD								STD
	0.046					TAHO ex		0.044					TAHO ex
	0.044								0.042
	0.042								0.04
	0.04								0.038
	0.038								0.036
	0.036								0.034
	178	180	182	184	186	188	190	192	178	180	182	184	186	188	190	192

scheme

C u γ u C z γ z STD 0.0052 0.54 0.0064 1.1 ROE 0.0027 0.66 0.0036 1.2 TAHO 0.006 1 0.012 1.62 Table 2. The 3 × 3 Test, see section 4.2. Evaluation of constants γ and C for e u (t) = C u t -γu and e z (t) = C z t -γz defined in 16. For STD and ROE approximations, the numerical results show that the absolute error decays as e u (t) = O(t -1/2 ) and e z (t) = O(t -1 ); while, for TAHO's it improves of t -1/2 .


Joaquín Rodriguez

Scheduling theory and constraint programming applied to rail traffic management

. The last step is a solution algorithm that is many times based a complete or partial enumeration procedure (EP) of feasible solutions where large subsets are discarded. This article will focus on the path ST → CP → EP to solve railway management problems. I will first briefly introduce each field, and after applications to the railway traffic management.

Scheduling theory

Scheduling theory has been studied since 1950s, the first problems were taken on the field of management of operations in workshops. Scheduling has been used in many other areas but a specific scheduling problem were defined for each situation. This result a huge number of scheduling algorithms for a large number of specific problems. A widely-used category of problem is the job-shop scheduling (JSS) problem. A JSS problem consists of a set of n jobs that can run on m machines. A job consists of operations which may also be referred to as tasks or activities. Each activity requires a machine (also named a resource) during its execution. There are two types of constraints in a JSS problem:

• precedence constraints between activities;

• resource constraints that assert that machines are unary resources, i.e. no two activities requiring the same resource may execute at the same time. The preemption is not permitted, i.e. the processing of an operation of a job cannot be interrupted. Each job may follow one or more machine sequences. If each job follows one machine sequence, operations are pre-assigned to particular resources, as a result scheduling means the decision when to execute operations. If each job may follow different machine sequences, there may be several alternative resources to which the operation can be allocated. In that case resource allocation is part of the scheduling task. Some widely used objectives are : the minimization of the maximal completion time (so called makespan), the minimization of the maximal lateness, the minimization of the sum of all the flow the times that a jobs spend in the workshop (so called flow times), the minimization of the sum of all the tardiness of the jobs. As many JJS problems have been shown to be NP-hard, researchers often develop heuristic methods to find good solutions rather than exact methods to find optimal solutions.

Scheduling theory in railway traffic management

One the first work which use the scheduling theory for rail traffic management modeling go back to the 70's with the work by Spzigel [Spzigel, 1973]. The problem considered was the timetable construction of a single-track line. In a single-track line, each pair of stations are linked by a track segment where at most one train can be at any time. Spzigel pointed out the similarity of the time-table construction and the JSS problem : each train is an ordered sequence of activities. In his model, track segments are unary resources used to perform an activity. In the case of stations, Spzigel's model assumes that stations can accept any number of trains, i.e. that their capacity is infinite. The objective function is the weighted sum of the journey times which is equivalent to the weighted sum of the times each train enters the last segment. Resolution uses a Branch and Bound method, with the values being chosen in priority from the disjunctive constraints between the trains that are in conflict on a track segment. The instances used relate to a line with 5 track segments, with half the trains traveling in each direction. The journey time on the segments is the same and the number of trains varies between 4 and 10.

Other authors [START_REF] Higgins | Optimal scheduling of trains on a single line track[END_REF], Kreuger et al., 1997, Isaai and Singh, 2000] follow the idea of modeling segments of single-track lines as unary resources or by posting disjunctive constraints between activities which require the segment.

When trains are running on the same direction a more detailed model is needed to take into account the safety constraints to protect against rear collision. To avoid rear collisions, the block signals enforce train headway such that at most one train can be in a block at any time. Therefore the idea was to model each block segment with a unary resource. It can be notice that this model include the model of Spzigel and it is the most widely used approach [START_REF] Brännlund | Railway timetabling using lagrangian relaxation[END_REF], Oliveira, 2001, Kreuger et al., 2001, Mascis and Pacciarelli, 2002, Barber et al., 2004,Geske, 2005,Mazzarello and Ottaviani, 2005,Schlenker, 2006[START_REF] D'ariano | A branch and bound algorithm for scheduling trains in a railway network[END_REF], Tornquist and Persson, 2007, Corman et al., 2010] In [Mascis and Pacciarelli, 2002] Mascis et al. studied specific JSS model with blocking and no-wait constraints. This category of JJS model was used in [Mascis et al., 2002] and Figure 2: Deadlock example forbidden by no swap blocking constraint other works are based on this approach [START_REF] Mazzarello | A traffic management system for real-time traffic optimisation in railways[END_REF][START_REF] D'ariano | A branch and bound algorithm for scheduling trains in a railway network[END_REF], Corman et al., 2010]. The no-wait constraints allow to model that each activity will start immediately after its direct predecessor in a job. The blocking constraints allow to model that each activity of a job after complete processing on a resource will remain in it until the resource for the next activity become available for processing. Due attention must be paid to use no swap blocking constraints or any other equivalent technique to not missing some deadlocks. The figures 2 and 2 illustrate a deadlock example which could be missed without the no swap blocking constraints. This kind of possible infeasible train movements of the JSS block model was recently pointed out in [Harrod, 2011] and solved with a directed hypergraph model.

Our first work using a JSS model for rail traffic management was in [START_REF] Rodriguez | Constraint programming for real-time train circulation management problem in railway nodes[END_REF]]. The problem was to schedule and route trains to minimize secondary delays during operation on complex junctions or stations. We defined a more accurate model than the JSS block model because a block is too imprecise to represent the incompatibilities correctly i.e. potential train conflicts or deadlocks. Track circuits are more appropriate for specifying the resource used to determine the incompatibilities between train times. Therefore, the idea was to model each track circuit for train detection with a unary resource. This is illustrated in the figure 4. This model take into account the temporal position of the head and tail of train, i.e. a train is not considered as a point but it has a length. The start and end variables of the activities are an expression of the time when the head enters or the tail clears the track circuits. This modeling feature allows to easily formulate temporal constraints for the signal watching time, several aspects signaling, the clearing time, the switching time and the sectional route release implemented by the interlocking system. One can notice that all these constraints are complying with the so-called «blocking time theory» defined in the UIC leaflet 406 [UIC, 2004]. The clearing time constraint between consecutive activities can prevent infeasible swap movement when there are deadlocks such as in figure .

Constraint programming

Constraint programming is a technique used to represent constraint satisfaction problems. In very general terms, a constraint satisfaction problem (CSP) is defined by the following ordered triple P = (X, D, C) such that:

• X = {x 1 , . . . , x m } is a set of variables, • D = {d 1 , . . . , d m }
is a set of finite domains where a domain d i is a set of possible values for the variable x i ,

• C = {c 1 , . . . , c e } is a set of constraints where each constraint c i limits the combinations of possible values for a subset of variables X.

•

For a given CSP instance, it is possible to either look for one solution, all the solutions, or an optimum solution with respect to an objective function. In the last case, a

Constraint programming in railway traffic management

To our knowledge, Fukumori [Fukumori, 1980, Fukumori and[START_REF] Fukumori | [END_REF] was among the first to propose the use of constraint programming in order to model a rail traffic management problem. This research dealt with constructing a timetable on the basis of specifications of a set of constraints that applied to the trains. The decision variables related to trains overtaking or crossing in stations. Consequently, the main variables in the model were departure and arrival times at stations. The criterion to be optimized was each train's journey time in order to limit excessive train standing time in stations. The model's constraints expressed:

(Ct. 1) the minimum headways which trains must observe in the sections of track between stations, (Ct. 2) the minimum durations for stopping, crossing and overtaking in stations;

(Ct. 3) the maximum number of trains which each station can take.

A Depth First Search (DFS) resolution algorithm was used. This makes trade-offs with regard to track use conflicts on the basis of the priorities assigned to the trains. The instantiation order of the variables is the same as the journey order of the trains. When this research was carried out, CP tools were not in common use and a constraint propagation algorithm was therefore developed specifically. This algorithm seeks consistency with regard to the bounds of the variable domain, which is similar to arc-B-consistency [Lhomme, 1993]. Fukumori's model was validated on problem instances from one day's operation of the Osaka and Nara lines [START_REF] Fukumori | Fundamental algorithm for train scheduling based on artificial intelligence[END_REF]. The examples given in the article relate to the instances for 14 stations and 40 trains.

A number of subsequent studies have dealt with the timetable construction problem using Fukumori's model as a basis. The CAPRES model developed by Hachemane [Hachemane, 1997], uses the same constraints with the addition of some constraints which are specific to the applications in question. The most striking extension is that stations may be modeled with greater precision. For example, this model can include incompatibilities between routes into and out of stations. However, this option is only used to prove the feasibility of a solution. The constraint consistency algorithm is of the type used by Fukumori and in this case too it is implemented without a CP library. The CAPRES model has been used in many capacity studies (see a survey in [START_REF] Rivier | Analysing the capacity of railway networks: Summing up the experience[END_REF]), and was subsequently included in the Viriato software package [START_REF] Sma Und Partner | Future module : CAPRES -Elaboration And Saturation Of Railway Network[END_REF].

The PRaCoSy project model described by [START_REF] Chiu | A constraint-based interactive train rescheduling tool[END_REF]] also makes use of the constraints defined by Fukumori. It involves a different application, that of assisting an operator and designing a timetable in order to cope with disruptions. Two new criteria have been considered in this case:

• Minimizing the modifications to the theoretical timetable,

• Minimizing the largest delay in comparison with the theoretical timetable.

Two resolution heuristics have been tested. The first uses an instantiation order based on the smallest value in the variable domain. The second uses an instantiation order which as far as possible conserves the values in the initial timetable. Thus, the variables for the conflicting trains and for the trains which share part of their route with the conflicting trains are the last to be instantiated. The second heuristic gave better results for the problem instances from the Nanjingxi-Shanghai line. Implementing this instance with the CP Ilog Solver library gave a problem with 598 variables and 3005 constraints. However, the paper quotes no indicators to give an idea of the difficulty of the problems dealt with (e.g. the number of conflicting trains).

Ingolitto et al. [START_REF] Ingolitto | An efficient method to schedule new trains on a heavily loaded railway network[END_REF]] also make use of Fukumori's model and CSPs. As with CAPRES, the problem involved timetable saturation. A resolution heuristic is used which divides the line in two in order to make a scheduling decision between the saturating trains. It should be noted that the algorithm only verifies the constraints for the use of single tracks. If the verification fails, the times of the trains are changed.

Isaai and Singh [START_REF] Isaai | An intelligent constraint-based search method for a single-line passenger-train scheduling problems[END_REF] propose a model of the same type as Fukumori in order to solve a timetable construction problem. Their resolution method combines the CP heuristic in order to obtain a solution that is feasible and a taboo heuristic to improve this solution in a neighborhood. The 5 data sets used for the tests were obtained from two lines with single and double tracks. These lines had a total of about 50 stations. The number of conflicts was estimated at between 25 and 62. The paper's results show that the taboo method improves the initial solution given by the CP heuristic.

The decision variables used in the above research principally relate to train scheduling with shared resources (station platforms, track sections). Constraint programming has also been applied to problems with additional variables for train routing,i.e;the resources allocation. Adding of these variables greatly increases the size of the search space which makes it very difficult to find a solution which is both acceptable and optimum.

To our knowledge, Gosselin in [Gosselin, 1993] describes the first application which uses constraint programming to deal simultaneously with the problems of train scheduling and train routing. In view of the number of variables and constraints, the problem was limited to the search for an acceptable solution. To be able to do this, it is essential to define the constraints that propagate between the times of station arrival and departure events and alternative routes. The resolution method is guided by a route evaluation function which enables it to find solutions more rapidly. In order to restrict the size of the search space, the number of possible routes between stations was limited to 4. The model was applied to 10 instances on the Melbourne suburban network, with the size of the instances varying between 200 and 2000 trains. Resolution times were less than one hour which satisfied customer requirements.

Poujade et al. [START_REF] Poujade | Sagitaire: Allocation des quais de la gare Paris-Nord[END_REF] also used a CP model to deal with a routing problem. The application involved allocating platform lines at the Gare du Nord in Paris. In this application, the train times at the reference points were fixed. The difficulty of the problem lay in the number of routes and the number and type of constraints to be satisfied. The initial problem in the application was over-constrained, i.e. there was no possible solution. The constraints were then divided into three categories: those that absolutely had to be satisfied, those which could be modified and those which had to be satisfied as far as possible. In the resolution method, the possible number of routes for a train was so great (several thousand) that these routes were only calculated after the platform lines had been assigned (terminus and/or intermediate stops). Then, route selection focused on those which were as direct as possible. In the event of failure, new platform line assignments were tried out. If the search ended in failure, the sets of conflicting trains were identified in order for the user to be able to choose which constraint to relax and then initiate a new search. This application is operational and used on a daily basis. The routes cover a zone with a radius of 16 km around the Gare du Nord which includes 10 stations, 100 platforms for 1700 trains and the model has a resolution time of 15 minutes. The «track circuit» JSS model [START_REF] Rodriguez | Constraint programming for real-time train circulation management problem in railway nodes[END_REF]] mentioned in section 2 was formulated as a Constraint Satisfaction Optimization Problem and implemented with a CP library (IBM Ilog Solver/Scheduler). It was successively improved to deal with more complex instance problems [Rodriguez, 2007a]. In [START_REF] Rodriguez | A study of the use of state resources in a constraint-based model for routing and scheduling trains[END_REF] we introduce the use of state resources to early prune scheduling decision leading to opposite direction conflicts. Thereafter [START_REF] Rodriguez | A study of an incremental texture-based heuristic for the train routing and scheduling problem[END_REF] we define a new heuristic which performs incrementally decisions at very local level i.e. at track circuit level. At each step, the algorithm performs decisions of allocation of a track circuit or of scheduling a pair of train runs on this track circuit. For this later solution method, the experimental study was a set of 18 instance problems with increasing difficulty. The difficulty of the instance problems is related to the number of trains (6 to 40 trains) and conflicts between trains (1 to 26). The instances were generated by selecting train services from the services of a peak hour of the station of Lille-Flandres, the main station of Lille in the North of France. The instances considered correspond to instances with 2000 to 22000 variables and 3000 to 26000 constraints. This later solution method have shown that good quality solutions can be obtained with processing times (less than 30 seconds) which are compatible with the operational constraints.

Prospect

This paper has briefly presented ST and CP models for the railway traffic management problems. We focused on one path ST→ CP→EP of the graph of figure 1. Despite contributions to improve the methods of solving problems, this model is still only predictive but not explanatory as it is argued in [Hooker, 2007]. Indeed all models are generally confined to a set of computational techniques but they do not help to understand the phenomenon they are studying : the railway traffic management. To our knowledge, there is still no modeling contribution that highlight the structure of a problem of railway traffic management and can be used to understand why a solution is good. Moreover, there is no reason to consider a single path in the graph of figure 1. On the contrary, all paths, i.e. all techniques, should be available to exploit the structure of the problem and solve it. Therefore, future prospects are to integrate different optimization techniques (ST, GT, IP) to obtain efficient predictive models but also explanatory models of the railway traffic management.

Figure 1 :

 1 Figure 1: Steps to model railway traffic management with scheduling theory

Figure 3 :

 3 Figure 3: Time-Distance diagram of the deadlock example of fig. 2

Figure 4 :

 4 Figure 4: Time-Distance diagram of a JSS model with track circuit resources