R Natalini
email: roberto.natalini@cnr.it

M Ribot
email: ribot@unice.fr

ASYMPTOTIC HIGH ORDER MASS-PRESERVING SCHEMES FOR A HYPERBOLIC MODEL OF CHEMOTAXIS

Keywords: hyperbolic systems with source, initial-boundary value problem, asymptotic behavior, finite difference schemes, chemotaxis AMS subject classifications. Primary: 65M06; Secondary: 35L60, 35L50, 92B05, 92C17

We introduce a new class of finite difference schemes for approximating the solutions to an initial-boundary value problem on a bounded interval for a one dimensional dissipative hyperbolic system with an external source term, which arises as a simple model of chemotaxis. Since the solutions to this problem may converge to non constant asymptotic states for large times, standard schemes usually fail to yield a good approximation. Therefore, we propose a new class of schemes, which use an asymptotic higher order correction, second and third order in our examples, to balance the effects of the source term and the influence of the asymptotic solutions. A special care is needed to deal with boundary conditions, to avoid harmful loss of mass. Convergence results are proven for these new schemes, and several numerical tests are presented and discussed to verify the effectiveness of their behavior.

Introduction.

It is often quite difficult to find an effective numerical approximation to hyperbolic equations with a source term. There are many different problems which could arise, for instance: stiffness of the source term, instability of the solutions, incorrect approximation of stationary solutions. Many ideas were introduced in the last thirty years to face these problems, and many of them are just working for a specific class of problems. Let us mention some families of schemes, sometimes overlapping: well balanced [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF]A steady-state capturing method for hyperbolic systems with geometrical source terms[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF], Runge-Kutta IMEX [START_REF] Pareschi | Implicit-explicit runge-kutta schemes and applications to hyperbolic system with relaxation[END_REF], upwinding source [START_REF] Roe | Upwind differencing schemes for hyperbolic conservation laws with source term[END_REF][START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF][START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF], and asymptotic preserving [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | A micro-macro decomposition-based asymptotic-preserving scheme for the multispecies Boltzmann equation[END_REF]. One of the main ideas, which can be found as a starting point for most of these schemes, is to plug the knowledge of the analytical behavior of the solutions into the scheme, to guarantee a better approximation, at least around some relevant asymptotic states of the problem.

In this paper, we want to study a specific problem, the numerical approximation of a one space dimensional hyperbolic system, which arises as a simple model for cell movement driven by chemotaxis:

   u t + v x = 0, v t + λ 2 u x = uφ x -v, φ t -D φ xx = au -bφ, (1.1)
where λ, a, b, D are all positive constants. Here, the function u denotes the density of cells in a given medium (at this level, they could be bacteria, eukaryotic cells, etc...), v their averaged flux, and φ a chemotactic stimulus, produced by the cells, biasing the movement of the population itself. Such kind of models were originally considered in [START_REF] Segel | A theoretical study of receptor mechanisms in bacterial chemotaxis[END_REF], and later reconsidered in [START_REF] Greenberg | Stability results for a diffusion equation with functional drift approximating a chemotaxis model[END_REF]. Recently, some generalizations of these models have been studied analytically in [START_REF] Hillen | Hyperbolic models for chemotaxis in 1-D, Nonlinear Anal[END_REF][START_REF] Hillen | Existence of weak solutions for a hyperbolic model of chemosensitive movement[END_REF] and later in [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], where the analytical features were almost completely worked out, at least around constant equilibrium states. Multidimensional extensions and more realistic generalizations have been considered for instance in [START_REF] Dolak | Cattaneo models for chemosensitive movement. Numerical solution and pattern formation[END_REF][START_REF] Gamba | Percolation, morphogenesis, and Burgers dynamics in blood vessels formation[END_REF][START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Cavalli | 3D simulations of early blood vessel formation[END_REF][START_REF] Russo | Analysis and numerical approximation of hydrodynamical models of biological movements[END_REF], and the present study is also aimed to shed a light on the approximation of these other models.

Clearly, this problem has to be complemented by initial conditions at t = 0, u(x, .) = u 0 (x), v(x, .) = v 0 (x), φ(x, .) = φ 0 (x), and, if we cast the problem in a bounded interval (0, L), by a suitable choice of boundary conditions. Here, we just deal with no-flux conditions, which are mainly relevant for biological purposes and we use the boundary conditions v(0, .) = v(L, .) = 0, φ x (0, .) = φ x (L, .) = 0, (1.2) but it is possible to use the present framework to extend the present results to more general conditions Let us now explain the specific numerical pathology we have to face when dealing with this problem. If we approximate the two first equations of (1.1) by using a standard method as an upwind explicit method with the source term being just approximated by the Euler scheme, and the third one using the standard Crank-Nicolson scheme, we obtain a scheme, see (3.6) below, with all the suitable properties for this problem, which gives some coherent results around small perturbations of stable constant states. However, it is possible to see that most of the stationary states of this problem are non constant. For instance, if the total mass of bacteria is large enough, with respect to the size of the domain and to the parameters of the system, a time asymptotically stable stationary state for u is given by the S-shaped distribution, here on the left in Figure 1 (more details will be given in Section 2 below). The correspond-Fig. 1. Numerical results for the upwind scheme on problem (1.1): asymptotic function u is displayed on the left, function φ in the middle and function v on the right. Exact asymptotic states are displayed in black stars and numerical results for the upwind scheme in blue diamonds. We can notice that the function v is far from vanishing as it should be. However, results in red are obtained with our alternative scheme described later on and are very accurate. The initial datum is a non-symmetric perturbation of constant state equal to 1135.

ing flux function v needs to be constant, as for all stationary states, and so equal to 0 under the no-flux conditions. However, using this basic scheme, both functions u and φ are well approximated, while there is a quite large error in the approximation of the function v, as seen in the right picture in Figure 1. This is not surprising at all, since the standard scheme reads as follows on the first equation:

u n+1 j = u n j - k 2h (v n j+1 -v n j-1 ) + λ 2 k h (u n j+1 -2u n j + u n j-1 ),
where the last term on the right hand side is just the numerical viscosity induced by the upwind scheme. So, when the function u is near to the equilibrium, and the second derivative of u is large, the space step h needs to be very small, actually hu xx = o(1) is needed, in order to make the function v constant. Therefore, if the asymptotic state is constant, the upwind scheme will be perfectly adapted, whereas for a non-constant asymptotic state, with a large second derivative, a new scheme has to be found. This is a typical situation when dealing with a problem which tends toward an equilibrium state where the flux in the conservation equation has to vanish in the time asymptotic limit. For example, this is the case of kinetic problems (radiative transfer [START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptoticpreserving schemes built on Case's elementary solutions[END_REF] or Boltzmann equation near the equilibrium state). Notice that here the problem does not arise from the stiffness of the source term, so implicit schemes are not expected to give better performances. Moreover, higher order schemes (TVD, WENO) give some limited results, but to radically improve the behavior of the schemes it is necessary to take into account the qualitative behavior of the solutions.

In this paper we propose a solution to this problem by introducing suitable modifications of upwind schemes, which are still globally first order accurate, but such that their truncation error computed on every stationary solution is order two or more. This approach was successfully used by the Asymptotic High Order schemes (AHO) introduced in [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF] to deal with dissipative hyperbolic systems, which are based on a suitable modification of the approximation of the source terms, to take into account for the behavior of the stationary solutions. However, in the present case, there are two new difficulties. First of all, there are infinitely many non constant solutions to this problem, and so it is impossible to design schemes which are exact when computed on all these solutions. Actually in our framework, we cannot move beyond the order 3 on asymptotic states. Besides, we are in a bounded domain and we have to approximate the boundary conditions in a consistent way with the AHO modifications. It turns out that to obtain a good approximation, it is also crucial to have an exact conservation of the numerical mass of the function u, which can be obtained only via a further modification of the natural (upwind) boundary conditions.

Our class of schemes verifies all the above properties. Two of these schemes are asymptotically second and third order respectively. While the second order AHO scheme uses a standard discretization of the diffusive equation, the third order AHO is third order on the solutions of the hyperbolic part with source, and for the diffusion equation uses a second order three points scheme which is fourth order on stationary solutions of the parabolic equation. Then, the coupling with the hyperbolic part of the third order AHO is made using a fourth order discretization of Φ x and the final scheme is really third order in our numerical tests. Please notice that the computational costs of all these schemes do not increase with respect of the standard ones.

Let us now mention that our approach is somewhat related to the well-balanced philosophy, see for instance [START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF] or for the specific case of hyperbolic chemotaxis problems, see [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF]. It is interesting to compare the Asymptotic High Order schemes we obtain in this paper with the results of some very recent works on Well-Balanced schemes [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF][START_REF] Gosse | Well-balanced numerical approximations display asymptotic decay toward Maxwellian distributions for a model of chemotaxis in a bounded interval[END_REF] for the same system considered here, which actually appeared only after the first submission of the present paper. First, both methods start with the problem of balancing the flux term for the conservation equation. The AHO schemes use the upwinding approach to approximate the stationary solutions by using a Taylor expansion in the truncation error, whereas the well-balanced approach for this special problem uses an approximate reconstruction of φ x which is inserted in the "exact" solution of the stationary problem for the hyperbolic part, given by the non-conservative solution to the Riemann problem. Both approaches have good results, however, concerning the accuracy, Well-Balanced schemes are only asymptotically second order, while our schemes are respectively second and third order, always near the equilibrium. An advantage of the AHO schemes is that they are really easy to be computed and that they can be extended in a multi-dimensional setting, since they are not related to the solutions to the Riemann problem. However, adequate boundary conditions and mass conservation are automatically obtained using the well-balanced approach, whereas a particular care for them is needed for our method, even if this can be achieved by simple modifications of the original upwinding conditions.

Finally, we want to remark that if in equation (1.1) we replace the term uφ x with a general function f = f (x, t), we obtain just the dissipative 1D wave equation with an external forcing term, see equation (4.2) below. If f = 0, both AHO and Well-Balanced schemes are exact on the stationary solutions, since these solutions are just linear functions of x. When a primitive function of the source term f is known, it is easy to find AHO and Well-Balanced schemes which are both exact on the stationary solutions. That is the reason why the main focus of our paper is in setting up the right framework for the upwinding approach in the approximation of the external source term f when no analytical primitive function is known, as in the case of chemotaxis, where it is impossible to select a unique stationary solution.

The plan of this paper is as follows. In the next section we recall some analytical features of problem (1.1) and in Section 3 its approximation is proposed using a basic upwind/Crank-Nicolson scheme. A general class of AHO schemes is introduced in Section 4 and their main properties are established. Next we study the problem of boundary conditions and determine a new set of conditions to be imposed to preserve the numerical mass of our approximated solutions. Specific instances of AHO schemes, second and third order, are introduced in Section 6, and their convergence is proved in the following section. Finally, Section 8 is devoted to some numerical experiments, which illustrate the behavior of our schemes and their actual accuracy.

2. Some analytical backgrounds. In [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], some results about global existence and stability of solutions to problem (1.1) were obtained on the real line R and on a bounded interval [0, L], and in this case with the no-flux boundary conditions (1.2). Let us recall the main features of this problem, which are going to be useful in what follows. The first remark is that the total mass of u is conserved in time, namely

I(t) = [0,L] u(x, t) dx = I(0) = [0,L] u(x, 0) dx. (2.1)
Symmetry is also conserved, which is, if the inital datum satisfies (u, v, φ)(L -x, 0) = (u, -v, φ)(x, 0), then for all times t > 0, we have

(u, v, φ)(L -x, t) = (u, -v, φ)(x, t).
Another remark is that for all stationary solutions, the function v vanishes. In [START_REF] Schaaf | Stationary solutions of chemotaxis systems[END_REF], many results about these stationary solutions were obtained and were recalled in [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]. First, there exists a class of constant stationary solutions, namely (u, v, φ) = (U, 0, a b U ), for some given constant U > 0. For bounded intervals, the time asymptotic stability of constant stationary states for the hyperbolic case of system (1.1) was proven in [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]. Namely we showed that (u, v, φ) = (U, 0,

a b U ) is a stable stationary solution for U ≤ U 1 = λ 2 a (D π 2 L 2 + b)
, with an exponential decay rate, whereas the constant states are in general unstable for U > U 1 .

Besides, there exist also infinitely many non-constant stationary solutions. They can be displayed as points on a graph [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] where a stationary solution (u, 0, φ) is represented by the couple (I(0), φ(0)) of the total mass of u defined at equation (2.1) and the value of φ at point x = 0. On this graph, non-constant stationary solutions are organized as branches, beginning at some bifurcation points on the line of constant stationary solutions By numerical experiments we obtained information about the stability of non constant stationary solutions and we found that the first branch, beginning at

U 1 = λ 2 a (D π 2 L 2 + b)
, appears to be also stable. However, if the initial datum is symmetric, the asymptotic solution cannot be the non constant solution on the first branch, which is not symmetric. In that case, it was proven in [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] that the asymptotic behavior of the solutions will be given by the constant stationary state if the mass is lower than

U 2 = λ 2 a (D 4π 2 L 2 + b).
Numerical simulations indicate that, if the mass of u is greater than U 2 , the non constant stationary solution of the second branch will be stable for symmetric data.

3. A basic upwind scheme. To introduce our schemes, first we diagonalize system (1.1), by setting w

= 1 2 (u - v λ ), z = 1 2 (u + v λ
). This yields

           w t -λw x = 1 2λ (-(φ x + λ)w + (λ -φ x )z) , z t + λz x = 1 2λ ((φ x + λ)w + (φ x -λ)z) , φ t -Dφ xx = a(w + z) -bφ, (3.1)
and the boundary conditions now read w(0, .) = z(0, .), w(L, .) = z(L, .), φ x (0, .) = φ x (L, .) = 0.

(3.2)

To introduce the numerical approximation, let us denote by h the space step and by k the time step, linked by the usual hyperbolic CFL condition λk h ≤ 1. We consider the discretization points x j = j h, 0 ≤ j ≤ M + 1, with x 0 = 0 and x M +1 = L. The discretization times will be given by t n = n k, n ∈ N.

Let us denote by w n j (resp. z n j and φ n j ) the approximation of w(x j , t n ) (resp. z(x j , t n ) and φ(x j , t n )). Therefore, the discretization vectors at time t n will be denoted by

W n = (w n 1 , • • • , w n M ) T , Z n = (z n 1 , • • • , z n M ) T and Φ n = (φ n 1 , • • • , φ n M )
T . We approximate the two first equations of system (3.1) by using an upwind explicit scheme and the third one using the Crank-Nicolson scheme in time and a finite difference method of order two in space. The boundary conditions are treated as follows. From conditions (3.2), we have z n 0 = w n 0 and w n M +1 = z n M +1 . Since w n+1 0 and z n+1 M +1 are directly computed using the upwind method, we simply compute the missing values of w and z on the boundaries by

z n+1 0 = w n+1 0 = 1 -λ k h w n 0 + λ k h w n 1 , w n+1 M +1 = z n+1 M +1 = 1 -λ k h z n M +1 + λ k h z n M . (3.3)
The function φ satisfies homogeneous Neumann boundary conditions and we use a standard second order rule to approximate its derivative, which yields:

φ n 0 = 4φ n 1 -φ n 2 3 , φ n M +1 = 4φ n M -φ n M -1 3 . (3.4)
This enables us to compute Φ n x by using the following approximation:

Φ n x = 1 2 h 4 3 (φ n 2 -φ n 1 ), φ n 3 -φ n 1 , • • • , φ n M -φ n M -2 , 4 3 (φ n M -φ n M -1 ) T , (3.5)
with (Φ x ) n 0 = (Φ x ) n M +1 = 0 on the boundaries. This vector is needed for the approximation of the two first equations of system (3.1).

Therefore, setting M and M 2 as the M × M matrices

M =        2/3 -2/3 (0) -1 2 -1 . . . . . . . . . -1 2 -1 (0) -2/3 2/3        and M 2 = 1 + b k 2 I + k 2h 2 DM,
we consider the following scheme, with 1

≤ j ≤ M                w n+1 j = w n j + λ k h (w n j+1 -w n j ) - k 2λ (Φ x ) n j (w n j + z n j ) - k 2 (w n j -z n j ), z n+1 j = z n j -λ k h (z n j -z n j-1 ) + k 2λ (Φ x ) n j (w n j + z n j ) + k 2 (w n j -z n j ), Φ n+1 = M -1 2 (1 -b k 2 )Φ n - k 2h 2 DMΦ n + a k 2 (W n + W n+1 + Z n + Z n+1 ) . (3.6) Since the spectrum of M is contained in the disk D(2, 2) = {µ ∈ C : |µ -2| ≤ 2}, the matrix M 2 is

invertible without any condition on k and h.

The consistency of this scheme is well-known. Moreover, it preserves some properties of the original system; namely, constant stationary solutions, conservation of mass, symmetry with respect to the transformation (x, t, u, v, φ) → (L-x, t, u, -v, φ).

In order to prove it, we need to rewrite scheme (3.6) and the boundary conditions (3.3) in the initial variables u and v using the relations u = z + w and v = λ(z -w).

Conservation of constant steady states. System (3.1) has constant station-

ary solutions of the form (u, v, φ) = (U, 0, a b U ), which read (w, z, φ) = ( U 2 , U 2 , a b U ) in
the diagonal variables. We can check with a simple computation that scheme (3.6) conserves them exactly.

Mass conservation. We have seen that an important feature of system (3.1) is that the mass I(t) of the function u is constant in time. Let us define a discrete version of the mass of the density function u at time t n on the interval [0, L] using the trapezoidal rule,

I n 1 = h u n 0 2 + M i=1 u n i + u n M +1 2
. Using scheme (3.6) and boundary conditions (3.3) in the initial variables u and v, we obtain that I n+1

1 = I n 1
, that is to say that the mass is preserved by the scheme.

Conservation of symmetry. Finally, we can also prove that the symmetry properties are satisfied, namely that if (SU 0 , SV 0 , SΦ

0 ) = (U 0 , -V 0 , Φ 0 ), with S =   1 . . . 1   , then for all n ∈ N, (SU n , SV n , SΦ n ) = (U n , -V n , Φ n ).
Proposition 3.1. The scheme (3.6), supplemented by the numerical boundary conditions (3.3) and (3.4), is a consistent scheme for system (3.1). It is also an exact scheme on constant equilibrium solutions, it conserves the total numerical mass and the symmetry of solutions for symmetric initial data.

4.

A more general class of schemes. In view of the bad behavior of the basic scheme (3.6) near non-constant stationary solutions, as mentioned in the Introduction (but see Section 8, for more details and numerical tests), we are going to deal with a more general class of schemes, which are based on the asymptotically high order schemes (AHO) introduced in [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF]. These schemes are more precise than standard schemes near non-constant asymptotic states and so they can avoid the failure of standard schemes. So, let us outline briefly the core ideas which shape AHO schemes for general hyperbolic systems.

Let U be a solution to the general linear hyperbolic system

U t + AU x = BU + F, (4.1)
and let U n+1 = H(U n , F n ) be a numerical scheme which is linear in U and consistent with system (4.1). Let T H be its local truncation error. Let us assume that the local truncation error is only first order on a smooth solution of (4.1), that is to say T H (U ) = O(h + k).We also consider generic stationary solutions Û to the same problem, namely such that A Ûx = B Û + F. Definition 4.1. We say that the scheme U n+1 = H(U n , F n )is (locally) Asymptotic High Order of order p, which will be denoted by (AHO) p , for system (4.1), if the scheme is of order p on every stationary solution Û , i.e. T H ( Û ) = O(h p ). Now, assume that Û is a stationary time asymptotic state for U , namely U -Û = τ (1/t), as t → ∞,in some appropriate functional space topology, where τ is a modulus of continuity. Then

T H (U ) = T H (U -Û ) + T H ( Û ) ≈ τ (1/t) + O(h p ),
for some modulus of continuity τ , possibly depending on h; i.e.: for large times the effective (local) truncation error for the scheme is higher than the basic order of the scheme. As we shall see later on, this time asymptotic improvement of the scheme is sufficient to balance the scheme around stationary solutions.

In the following, we are going to deal with a general right hand side term f , instead of the term φ x u in system (1.1), and so we shall omit the parabolic equation for φ, which will be treated as before. Therefore, we consider the following linear system, for x ∈ [0, L],

u t + v x = 0, v t + λ 2 u x = f -v (4.2)
with the boundary conditions v(0, .) = v(L, .) = f (0, .) = f (L, .) = 0. Thus, we will consider a more general form of the source term and we will propose a class of schemes which are independent of the form of f and which consequently work for every source term f . In the following, some numerical tests will be presented in the case of an explicitly given function f . We reconsider the system in diagonal variables, that is to say

     w t -λw x = 1 2 (z -w) - 1 2λ f z t + λz x = 1 2 (w -z) + 1 2λ f. (4.3)
Let us denote by ω = w z and rewrite the system in the following form

ω t + Λω x = Bω + F, (4.4
)

with Λ = -λ 0 0 λ , B = 1 2 -1 1 1 -1 , F = 1 2λ
-f f . We also denote by ω n i and F n i two approximations of ω(x i , t n ) and F (x i , t n ), and consider schemes of the following more general form:

ω n+1 i -ω n i k + Λ 2h ω n i+1 -ω n i-1 - q 2h (ω n i+1 -2ω n i +ω n i-1 ) = ℓ=-1,0,1 B ℓ ω n i+ℓ + D ℓ F n i+ℓ , (4.5
) where B ℓ and D ℓ are 2 × 2 matrices acting on the vectors ω n i+ℓ and F n i+ℓ respectively at the points x i+ℓ , and q is the artificial viscosity.

Notice that the main feature of the scheme (4.5) is to distribute the approximation of the source terms on all the points of the stencil, to take into account the upwinding influences acting on the solutions.

We will now find some conditions on the matrices B ℓ and D ℓ in order to have a consistent, monotone scheme, that is of higher order on the stationary solutions. An explicit (AHO) 2 scheme of this type will be given in detail by equations (6.2)-( 6.3)and an explicit (AHO) 3 scheme by equations (6.7)-(6.8).

Consistency, monotonicity and second-order accuracy of the AHOschemes. The matrices B ℓ and D ℓ of the scheme (4.5) will be consistent with the equation (4.4) under the following conditions: ℓ=-1,0,1

B ℓ = B, ℓ=-1,0,1 D ℓ = I 2,2 .
(4.6)

As in [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF], monotonicity conditions, taken for a fixed f , read:

B ℓ i,j ≥ 0, ℓ = -1, 0, 1, i = j, (4.7a) 1 - k h q + kB 0 i,i ≥ 0, ∓(-1) i k h λ 2 + k 2h q + kB ±1 i,i ≥ 0, i = 1, 2. (4.7b)
Now, let us compute the local truncation error of the scheme (4.5) for a stationary solution, that is to say for a function ω such that ωt = 0 and Λω x = B ω + F . Let us set ωn i = ω(x i , t n ). We shall use that (Λ -1 B) 2 = 0 to obtain, using a standard Taylor expansion:

T ωn i = Λ 2h ωn i+1 -ωn i-1 - q 2h (ω n i+1 -2ω n i + ωn i-1 ) - ℓ=-1,0,1 B ℓ ωn i+ℓ - ℓ=-1,0,1 D ℓ F (x i+ℓ , t n ) = B -B -1 + B 0 + B 1 ωn i + I 2,2 -D -1 + D 0 + D 1 F (x i , t n ) - h 2 2 B 1 -B -1 Λ -1 B ωn i + qΛ -1 BΛ -1 + 2 B 1 -B -1 Λ -1 F (x i , t n ) + qΛ -1 + 2 D 1 -D -1 F x (x i , t n ) + O(h 2 ).
We therefore impose the following conditions to have second-order accuracy on every stationary solution:

B 1 -B -1 = - q 2 Λ -1 B = q 4λ -1 1 -1 1 , D 1 -D -1 = - q 2 Λ -1 = q 2λ 1 0 0 -1 . (4.8)
Rewriting the scheme in the (u, v) unknowns. Set

X n = (w n 1 , • • • , w n M , z n 1 , • • • , z n M ) T , F n = 1 2λ (-f n 1 , • • • , -f n M , f n 1 , • • • , f n M ) T .
Let us rewrite scheme (4.5) as

X n+1 = X n - k 2h L δ 1 X n -δ -1 X n + kq 2h ∆ 2 X n +k ℓ=-1,0,1 B ℓ δ ℓ X n +D ℓ δ ℓ F n , (4.9)
with the following notations:

L = Λ ⊗ I M,M , B ℓ = B ℓ ⊗ I M,M , D ℓ = D ℓ ⊗ I M,M , (δ ℓ V ) i = V i+ℓ , ℓ = -1, 0, 1, ∆ 2 V = δ 1 V -2V + δ -1 V. Let R = 1 1 -λ λ be the diagonalizing matrix such that u v = R w z
and denote by R -1 = 1 2

1 -1/λ 1 1/λ its inverse.

We rewrite the scheme (4.9) in terms of the initial variables u, v and f as:

U n+1 = U n + k 2 - 1 h δ 1 V n -δ -1 V n + q h ∆ 2 U n + ℓ=-1,0,1 β ℓ u,u δ ℓ U n + 1 λ ℓ=-1,0,1 β ℓ u,v δ ℓ V n + 1 λ ℓ=-1,0,1 γ ℓ u δ ℓ f n , V n+1 = V n + k 2 - λ 2 h δ 1 U n -δ -1 U n + q h ∆ 2 V n + λ ℓ=-1,0,1 β ℓ v,u δ ℓ U n + ℓ=-1,0,1 β ℓ v,v δ ℓ V n + ℓ=-1,0,1 γ ℓ v δ ℓ f n , (4.10)
where

U n = (u n 1 , • • • , u n M ) T , V n = (v n 1 , • • • , v n M ) T , f n = (f n 1 , • • • , f n M )
T , and

RB ℓ R -1 = 1 2 β ℓ u,u β ℓ u,v /λ λβ ℓ v,u β ℓ v,v , RD ℓ R -1 = 1 2 η ℓ u γ ℓ u /λ η ℓ v γ ℓ v . (4.11)
Let us notice that the two parameters η ℓ u and η ℓ v will not appear in the equations we consider later and therefore can be chosen at our convenience.

In these last variables, consistency relations (4.6) and relations (4.8) are equivalent to the following conditions:

                 β 1 u,u β 1 u,v β 1 v,u β 1 v,v + β 0 u,u β 0 u,v β 0 v,u β 0 v,v + β -1 u,u β -1 u,v β -1 v,u β -1 v,v = 0 0 0 -2 , β 1 u,u β 1 u,v β 1 v,u β 1 v,v - β -1 u,u β -1 u,v β -1 v,u β -1 v,v = q λ 0 1 0 0 , γ 1 u γ 1 v + γ 0 u γ 0 v + γ -1 u γ -1 v = 0 2 , γ 1 u γ 1 v - γ -1 u γ -1 v = q λ -1 0 .
(4.12)

Conservation of symmetry. As explained before, the system possesses a symmetry invariance, that we want to keep at the discrete level. Following computations of Section 3 according to equations (4.10) and (4.12), we find that the following conditions for the coefficients of our scheme:

β 0 u,v = 0, β 0 v,u = 0, γ 0 u = 0, (4.13)
lead to SU n+1 -U n+1 = 0 and SV n+1 + V n+1 = 0, where S is defined in section 3.

We state the final conditions on these new schemes. Proposition 4.2. A scheme in the form (4.10) is a consistent scheme with system (4.2) and of order 2 on stationary solutions, i.e. (AHO) 2 , if conditions (4.12) are verified. Moreover, if it also verifies conditions (4.13), it conserves the symmetry of the solutions for symmetric initial data. Monotonicity, for a given f , holds under conditions (4.7).

5. Mass conservation: how to impose boundary conditions. The boundary conditions for v and for f are clear: v n 0 = v n M +1 = 0 and f n 0 = f n M +1 = 0, whereas those for u must be carefully chosen to obtain an effective approximation. Actually, at an analytical level, the function u verifies a Neumann condition u x = 0 at the boundary, but it is difficult to use this remark. The first simple idea is to use the standard upwind conditions (3.3). Therefore we compute the evolution in time of solutions to problem (4.2) using a scheme of the form (4.10), which verifies all the conditions of Proposition 4.2, with the boundary conditions given by (3.3). In Figure 2 we consider as initial datum a non symmetric perturbation of the constant equal to 1500. We display on the left the exact (non constant) asymptotic state, and the computed one using the AHO-scheme with upwind boundary conditions at convergent time T = 20. We used space and time steps h = 1/90 and k = 5 × 10 -4 . We can see that using the AHO scheme, with the upwind boundary conditions, yields a solution which is very different from the right one. Now, we have compensated the effect of the artificial viscosity, but we have a different problem. Let us recall that one of the analytical feature of this problem is the conservation of the total mass of the function u, namely the quantity

I(t) = [0,L]
u(x, t) dx has to be preserved along the time. However, if we plot the total mass of our numerical approximation vs. time, see the right part of Figure 2, we observe a clear and unphysical loss of mass.

Conservation of mass. Now, we try to modify our boundary conditions to guarantee the conservation of mass. The aim is to define the values of u in 0 and L in order to preserve exactly a discrete mass in time, that is to say such that for all n ∈ N, I n 1 = I 0 1 . Our first attempt was to take the natural discrete mass defined at section 3 and to look for an expression of u n 0 and of u n M +1 as a function of u n i , v n i and f n i Fig. 2. On the left the exact (strongly skewed on the left) and the computed (wrong) asymptotic states, obtained with boundary conditions (3.3) are displayed. On the right, we can see the evolution of mass with respect to the iteration number. We can notice that the mass is strongly decreasing, whereas it should be constant.

for 1 ≤ i ≤ M (static boundary condition). However, it can be proven that it is not possible to find such an expression which would satisfy the condition I n+1 1 -I n 1 = 0. Also, the attempt to conserve mass choosing carefully the coefficients of matrices B ℓ and D ℓ defined at equation (4.5) also fails. We therefore need dynamic boundary conditions if we want to preserve such a mass, i.e. to express u n+1 0 as a function of

u n i , v n i and f n i for 1 ≤ i ≤ M . Let us compute I n+1 1 -I n
1 , using definition of I n 1 of Section 3 and scheme (4.10). Thanks to v n 0 = v n M +1 = f n 0 = f n M +1 = 0 and equations (4.12) and (4.13), we find

I n+1 1 -I n 1 = hk 2 1 k (u n+1 0 -u n 0 ) + 1 h v n 1 + q h (u n 0 -u n 1 ) + β 1 u,u (u n 0 -u n 1 ) - 1 λ β 1 u,v v n 1 - 1 λ γ 1 u f n 1 + hk 2 1 k (u n+1 M +1 -u n M +1 ) - 1 h v n M + q h (u n M +1 -u n M ) + β 1 u,u (u n M +1 -u n M ) + 1 λ β 1 u,v v n M + 1 λ γ 1 u f n M .
If we take a generic consistent boundary condition, we will have that I n+1 1 -I n 1 = O(h + k) and therefore after N time iterations, we will have

I N +1 1 -I 0 1 = O(1)
, which can be very large.

Therefore, we propose the following boundary conditions that ensure the exact conservation of mass, i.e. I n+1 1 -I n 1 = 0:

u n+1 0 = 1 -q k h -kβ 1 u,u u n 0 + k q h + β 1 u,u u n 1 -k 1 h - β 1 u,v λ v n 1 + k λ γ 1 u f n 1 , u n+1 M +1 = 1 -q k h -kβ 1 u,u u n M +1 + k q h + β 1 u,u u n M + k 1 h - β 1 u,v λ v n M - k λ γ 1 u f n M .
(5.1) Let us notice that, if we try to conserve the mass in a weaker sense, that is to say such that I n+1 1 -I n 1 = O(hk), with q = λ, we find the same boundary conditions (3.3) we used with the upwind scheme, which do not work, as previously seen.

Remark 1. Considering a different discrete approximation of the total mass of u, for example

I n 2 = h 3 2 u n 1 + M -1 i=2 u n i + 3 2
u n M , we can find other sets of boundary conditions which enable to conserve the corresponding discrete mass exactly.

Consistency of boundary conditions (5.1). Thanks to asymptotic expansions and to equation u t + v x = 0, we find that:

1 k u n+1 0 -1 -q k h -kβ 1 u,u u n 0 -k q h + β 1 u,u u n 1 + k 1 h - β 1 u,v λ v n 1 - k λ γ 1 u f n 1 = -(q + hβ 1 u,u )u x (x 0 , t n ) + k 2 u tt (x 0 , t n ) -q h k u xx (x 0 , t n ) - 2h λk β 1 u,v v x (x 0 , t n ) - 2h λk γ 1 u f x (x 0 , t n ) + h 2 v xx (x 0 , t n ) + O(k 2 + h 2 ),
which is an approximation of -qu x (x 0 , t n ) of order one. Now, we consider the case of stationary solutions, which satisfy ût = vt = 0 and also vx = 0, λ 2 ûx = f . In that case, the expansion above becomes :

1 k ûn+1 0 -1 -q k h -kβ 1 u,u ûn 0 -k q h + β 1 u,u ûn 1 + k 1 h - β 1 u,v λ vn 1 - k λ γ 1 u f n 1 = -(q + hβ 1 u,u )û x (x 0 , t n ) -q h 2λ 2 f x (x 0 , t n ) - h λ γ 1 u f x (x 0 , t n ) + O(k 2 + h 2 ) = -(q + hβ 1 u,u )û x (x 0 , t n ) + O(k 2 + h 2 ),
since for any (AHO) p scheme, with p ≥ 2 which preserves the symmetry property, we have γ 1 u = -q 2λ thanks to equations (4.12) and (4.13). The following term is equal to

h 2 12λ 2 qf xx (x 0 , t n ) -6β 1 u,u f x (x 0 , t n
) and therefore cannot be cancelled. Consequently, we can say that the boundary condition has the same behavior as the AHO schemes we propose in this article, namely it is generally of order one and of order two on stationary solutions.

Remark 2. We notice that, if f = 0, then u =constant and v = 0 are stationary constant solutions of equation (4.2). Considering the scheme (4.10) and following the same computations as in Section 3, we prove that stationary constant states are still exactly preserved.

Specific examples of (AHO) p schemes. Let us give examples of schemes

which satisfy the conditions of Proposition 4.2. We will present a concrete (AHO) 2 scheme at equation (6.2), complemented with boundary conditions (6.3) and a concrete (AHO) 3 scheme at equation (6.7), complemented with boundary conditions (6.8).

A compatible choice for matrices B ℓ and D ℓ with equations (4.6)-(4.8) leading to an (AHO) 2 scheme will be

B 0 = 1 4 -1 1 1 -1 , B ±1 = 1 8 -1 ∓ q λ 1 ± q λ 1 ∓ q λ -1 ± q λ , D 0 = 1 2 1 0 0 1 , D ±1 = 1 4 1 ± q λ 0 0 1 ∓ q λ . (6.1)
Note that, according to definition (4.11), we have therefore the following coefficients for this example: and q ≥ λ and h ≤ 2(q + λ). Finally conditions (4.13) are also clearly satisfied. In order to have q as small as possible, we choose q = λ and we therefore have the following restrictions on the time and space steps: h ≤ 4λ and k ≤ 4h h + 4λ .

                 β 1 u,u β 1 u,v β 1 v,u β 1 v,v = 0 q/2λ 0 -1/2 , β -1 u,u β -1 u,v β -1 v,u β -1 v,v = 0 -q/2λ 0 -1/2 , β 0 u,u β 0 u,v β 0 v,u β 0 v,v = 0 0 0 -1 , γ 1 u γ 1 v = -q/2λ 1/2 , γ 0 u γ 0 v = 0 1 , γ -1 u γ -1 v = q/2λ 1/
In that case, the scheme can be rewritten as:

       w n+1 i -w n i k -λ w n i+1 -w n i h = 1 4 (z n i -w n i ) + 1 4 (z n i+1 -w n i+1 ) - 1 4λ f n i - 1 4λ f n i+1 , z n+1 i -z n i k + λ z n i -z n i-1 h = 1 4 (w n i -z n i ) + 1 4 (w n i-1 -z n i-1 ) + 1 4λ f n i + 1 4λ f n i-1 , (6.2
) which for this particular example turns out to be exactly the scheme first proposed by Roe [START_REF] Roe | Upwind differencing schemes for hyperbolic conservation laws with source term[END_REF].

However, the major concern of the present work is to propose adapted boundary conditions for Roe's scheme, namely in the case of coefficients (6.1) with q = λ, for the boundary conditions (5.1), we obtain :

u n+1 0 = 1 -λ k h u n 0 + λ k h u n 1 -k 1 h - 1 2λ v n 1 - k 2λ f n 1 , u n+1 M +1 = 1 -λ k h u n M +1 + λ k h u n M + k 1 h - 1 2λ v n M + k 2λ f n M . (6.3)
An extension to third-order accuracy. Let us compute the third order truncation error as in section 4. We obtain

T ωn i = (B -(B -1 + B 0 + B 1 ))ω(x i , t n ) + (I 2,2 -(D -1 + D 0 + D 1 ))F (x i , t n ) - h 2 2(B 1 -B -1 )Λ -1 B ω(x i , t n ) + (qΛ -1 B + 2(B 1 -B -1 ))Λ -1 F (x i , t n ) +(qΛ -1 + 2(D 1 -D -1 ))F x (x i , t n ) - h 2 6 3(B 1 + B -1 )Λ -1 BΛ -1 F (x i , t n ) +(3(B 1 + B -1 ) -B)Λ -1 F x (x i , t n ) + (3(D 1 + D -1 ) -I 2,2 )F xx (x i , t n ) + O(h 3 ).
In order to have third-order accuracy on the truncation error for stationary solutions, we generalize consistency conditions to:

ℓ=-1,0,1

B ℓ = B + hC, ℓ=-1,0,1 D ℓ = I 2,2 + hE, (6.4)
for some suitable matrices C and E. We therefore impose the following five conditions:

C = EB, D 1 -D -1 = - q 2 Λ -1 , D 1 + D -1 = I 2,2 3 , B 1 -B -1 = -EΛ - q 2 Λ -1 B, B 1 + B -1 = B 3 .
Therefore, once E is chosen, we should take

C = EB, D 0 = 2 3 I 2,2 + hE , D ±1 = 1 2 I 2,2 3 ∓ q 2 Λ -1 B 0 = 2 3 I 2,2 + hE B, B ±1 = 1 2 B 3 ∓ EΛ ∓ q 2 Λ -1 B . (6.5)
Let us find some conditions on E to have the same properties as before. First, in order to conserve relations (4.12), we impose that

E 1,1 + E 2,1 = E 1,2 + E 2,2 and E 1,1 + E 1,2 = E 2,1 + E 2,2 , which leads to E 1,1 = E 2,2 and E 1,2 = E 2,1 .
Under these conditions, all relations of (4.12) hold true, but also symmetry conservation conditions (4.13).

In order to have a diffusion as small as possible, we inforce q = λ. To enhance monotonicity conditions (4.7) and relations

E 1,1 = E 2,2 and E 1,2 = E 2,1 , we should have E = 1 12λ 1 1 1 1
with the following conditions on the time and space steps:

h ≤ 6λ, k ≤ h λ + h/3 .
In that case, using equation (4.11), we can compute coefficients to be used in boundary conditions (5.1) as

                 β 1 u,u β 1 u,v β 1 v,u β 1 v,v = 0 -1/6 + q/λ 0 -1/3 , β -1 u,u β -1 u,v β -1 v,u β -1 v,v = 0 1/6 -q/λ 0 -1/3 , β 0 u,u β 0 u,v β 0 v,u β 0 v,v = 0 0 0 -4/3 , γ 1 u γ 1 v = -q/2λ 1/3 , γ 0 u γ 0 v = 0 4/3 , γ -1 u γ -1 v = q/2λ 1/3 .
(6.6) To sum up, we have handled an (AHO) 3 scheme, using coefficients (6.5) in scheme (4.5) and coefficients (6.6) in boundary conditions (5.1).To be more precise, the scheme can be written in the case q = λ as :

                       w n+1 i -w n i k -λ w n i+1 -w n i h = 1 3 (z n i -w n i ) + 1 6 (z n i+1 -w n i+1 ) + 1 24λ f n i-1 - 1 3λ f n i - 5 24λ f n i+1 , z n+1 i -z n i k + λ z n i -z n i-1 h = 1 3 (w n i -z n i ) + 1 6 (w n i-1 -z n i-1 ) + 5 24λ f n i-1 + 1 3λ f n i - 1 24λ f n i+1 , (6.7)
with boundary conditions

u n+1 0 = 1 -λ k h u n 0 + λ k h u n 1 -k 1 h - 5 6λ v n 1 - k 2λ f n 1 , u n+1 M +1 = 1 -λ k h u n M +1 + λ k h u n M + k 1 h - 5 6λ v n M + k 2λ f n M . (6.8)
Let us remark that it is not possible to find a fourth order scheme under the form (4.5), since it leads to incompatible relations between the coefficients of the matrices B ℓ and D ℓ . It may be possible to find a higher order scheme by increasing the stencil of discretisation of the source, but in that case, it would be more difficult to deal with boundary conditions and stability properties.

7. Consistency, stability and convergence. Let w, z be the solutions to system (4.3) and set ω = w z . We define the truncation error for the scheme (4.5) by

T ω(x i , t n ) = ω(x i , t n+1 ) -ω(x i , t n ) k + Λ 2h (ω(x i+1 , t n ) -ω(x i-1 , t n )) - q 2h (ω(x i+1 , t n ) -2ω(x i , t n ) + ω(x i-1 , t n )) - ℓ=-1,0,1 B ℓ ω(x i+ℓ , t n ) - ℓ=-1,0,1 D ℓ F (x i+ℓ , t n )
Theorem 7.1 (Consistency). For every AHO scheme, which satisfies conditions (6.4), we have that, for all x i , 1 ≤ i ≤ M and t n , n > 0,

|T ω(x i , t n )| ≤ C(k + h),
where C depends on the coefficients of the scheme, the second derivatives of the solution ω and the first derivatives of the function F .

Proof. Using a Taylor expansion, equation (4.3) and consistency relations (6.4), we find that

T ω(x i , t n ) = k 2 ∂ tt ω(x i , t n ) + h 2 6 Λ∂ xxx ω(x i , t n ) - h 2 Q∂ xx ω(x i , t n ) -hC ω(x i , t n ) -h(B 1 -B -1 )∂ x ω(x i , t n ) -hE F (x i , t n ) -h(D 1 -D -1 )∂ x F (x i , t n ) + O(k 2 + h 2 ),
which gives the proof.

Stability for the boundary conditions (5.1). To establish stability results for the set of boundary conditions (5.1), let us rewrite them in diagonal variables. The conditions (5.1) become now, for q = λ, and setting δ = λk h :

z n+1 0 = w n+1 0 = (1 -δ)w n 0 + δw n 1 -kβ 1 u,u w n 0 + k 2 (β 1 u,u -β 1 u,v )w n 1 + k 2 (β 1 u,u + β 1 u,v )z n 1 + k 2λ γ 1 u f n 1 , w n+1 M +1 = z n+1 M +1 = (1 -δ)z n M +1 + δz n M -kβ 1 u,u z n M +1 + k 2 (β 1 u,u -β 1 u,v )z n M + k 2 (β 1 u,u + β 1 u,v )w n M - k 2λ γ 1 u f n M . (7.1)
Therefore, in order to prove the L ∞ and the L 1 stability in the case q = λ, we rewrite the scheme (4.9) under a matrix form using previous conditions (7.1). Setting

Y n = w n 0 , • • • , w n M , z n 1 , • • • , z n M +1 T and F n = 1 2λ (0, -f n 1 , • • • , -f n M , f n 1 , • • • , f n M , 0)
T , we have:

Y n+1 = A(δ)Y n + kBY n + kDF n , (7.2)
where

A(δ) = (1 -δ)1 2M,2M + δ               0 1 (0) . . . . . . (0) 0 1 (0) 0 1 1 0 (0) 1 0 (0) . . .      0 • • • 0 0 • • • 0 b -1 1,1 b 0 1,1 b 1 1,1 (0) b 0 1,2 b 1 1,2 (0) . . . . . . . . . b -1 1,2 . . . . . . . . . (0) b -1 1,1 b 0 1,1 (0) b -1 1,2 b 0 1,2 b 1 1,2 b -1 2,1 b 0 2,1 b 1 2,1 (0) b 0 2,2 b 1 2,2 (0)
. . . . . . . . . b -1 2,2 . . . . . . b 1 2,1 . . . (0) b -1 2,1 b 0 2,1 (0) b -1 2,2 b 0 2,2 b 1 2,2 0 • • • 0 0 • • • 0                       +         -β 1 u,u 1 2 (β 1 u,u -β 1 u,v ) 1 2 β 1 u,u + β 1 u,v b -1 1,2 b 1 1,1 b -1 2,2 b 1 2,1 1 2 β 1 u,u + β 1 u,v 1 2 (β 1 u,u -β 1 u,v ) -β 1 u,u         and D =                    0 • • • γ 1 u • • • 0 . . . d 0 1,1 d 1 1,1 (0) d 0 1,2 d 1 1,2 (0)
. . .

d -1 1,1 . . . . . . d -1 1,2 . . . . . . (0) . . . . . . (0) . . . . . . . . . d 0 2,1 d 1 2,1 (0) d 0 2,2 d 1 2,2 (0)
d -1 2,1 . . . . . . d -1 2,2 . . . . . . (0) . . . . . . (0)
. . . . . . Remark 3. We first notice that, since 0 ≤ δ ≤ 1, we have ||A(δ)|| 1 = 1 and that ||A(δ)|| ∞ = 1, whatever the size of the matrix. We also remark that ||B|| 1 , ||B|| ∞ , ||D|| 1 , ||D|| ∞ are independent of the size of the matrix. Theorem 7.2 (L ∞ and L 1 Stability). Let 0 ≤ δ ≤ 1 and t n = nk. The solution Y n to the numerical scheme (7.2), which is just scheme (4.9) with boundary conditions (5.1), can be bounded by

||Y n || 1 ≤ e tn||B||1 ||Y 0 || 1 + t n ||D|| 1 max 0≤j≤n-1 ||F j || 1 .
The same bound holds for the L ∞ norm.

Proof. We obtain

Y n = A(δ) n Y 0 +k n j=1 A(δ) j-1 BY n-j + DF n-j , which leads to ||Y n || 1 ≤ ||Y 0 || 1 + k||B|| 1 n-1 j=0 ||Y j || 1 + nk||D|| 1 max 0≤j≤n-1 ||F j || 1 .
We can also prove by induction that

||Y n || 1 ≤ (1 + k||B|| 1 ) n ||Y 0 || 1 + (1 + k||B|| 1 ) n t n ||D|| 1 max 0≤j≤n-1 ||F j || 1 ,
which achieves the proof. Now that consistency and stability are proven, let us state a convergence result. We define the convergence error at time t n as

e n (h, k) = max 1≤i≤2M +2 |Y n i -T n i |,
where T n is the vector

T n = (w(x 0 , t n ), • • • , w(x M , t n ), z(x 1 , t n ), • • • , z(x M +1 , t n )) T
given by the solution of system (4.3) with initial value T 0 . The vector Y n is given by the scheme (7.2) with initial vector Y 0 such that Y 0 = T 0 . By using the Lax equivalence Theorem, we can prove the following convergence result. Theorem 7.3 (Convergence). For every smooth solution to problem (4.2), with a smooth function f , the approximation given by the scheme (4.9), with the boundary conditions (5.1), converges and we have the following estimate of the numerical error:

e n (h, k) ≤ Ce tn||B||∞ t n (k + h). (7.3)
8. Numerical simulations.

8.1. Comparison upwind vs AHO schemes. Let us first compare the standard upwind scheme (3.6) with the second and third order AHO schemes with boundary conditions (5.1) for the system (4.2). To be more precise, we will consider 1. the (AHO) 2 defined by matrices (6.1), namely scheme (6.2), 2. the (AHO) 3 defined by matrices (6.5) and E = 1 12λ 1 1 1 1

, namely scheme (6.7). Here, we will choose an explicitly given f vanishing at the boundaries, for which we can compute explicitly the stationary solution of system (4.2).

Let us consider f (x) = Cx(1 -x) on the domain [0, L]. In that case, the explicit

stationary solution is u(x) = M L + C λ 2 x 2 2 - x 3 3 - L 2 12 (2 -L) , v(x) = 0, x ∈ [0, 1],
where M = L 0 u 0 (x)dx is the mass of the initial datum u 0 . We display our numerical simulations in Figure 3, with L = 1, λ = 10, C = 40000, u 0 the constant function equal to 35 and v 0 the null function such that the stationary solution is equal to

u(x) = 35+400 x 2 2 - x 3 3 - 1 12
, v(x) = 0, x ∈ [0, 1]. We notice that for the function u the three schemes, namely upwind, (AHO) 2 and (AHO) 3 give equivalent results, obviously with a different accuracy. However, for the function v, the upwind scheme does not vanish as it should be, whereas AHO schemes do.

Fig. 3. The source term f is such that f (x) = Cx(1 -x), the initial datum for u is a constant equal to 35 and the initial datum for v is the null function. On the left, the asymptotic state for u is displayed for different schemes and on the right, the same for v. We notice that results are equivalent, but with different orders of convergence, for the function u and that AHO schemes perform clearly better than upwind scheme for the function v. We used as space and time steps h = 1/90, k = 5 × 10 -4 and final convergent time T = 30.

In Figure 4 we present some plots of the errors in the L 1 norm at a logarithmic scale for the different schemes using as an exact solution the explicit stationary solution and the same initial data as above. Considering the error in the L ∞ norm gives exactly the same results. We can see that the three schemes have the expected orders for the function u, namely one for upwind, two for (AHO) 2 and three for (AHO) 3 . As for the function v, upwind and (AHO) 3 are of the right order whereas the error for (AHO) 2 is so small that the numerical order estimate cannot be handled properly. In !" Log-log plot of the error in L 1 norm as a function of the space step. On the left, the error for the function u is displayed for three different schemes, namely upwind, (AHO) 2 and (AHO) 3 ; we can see on the right the error for v. The initial data are the same as in Figure 3. We used different space steps h betwen 0.02 and 10 -3 , time steps equal to k = h/2λ and final convergent time T = 100.

(# !" ($ !" (! !" (!! !" (!" !" (% !" (& !" (' !" (( !" () !" (* !" (# !" ($ !" (! + + ,-./01 234$ 234# 56176+! 56176+$ 56176+# !" (# !" ($ !" (! !" (!% !" (!$ !" (!" !" (& !"
the following Table 1, we display the slopes of the straight lines shown at Figure 4. 4 Upwind scheme (AHO) 2 scheme (AHO) 3 For the (AHO) 2 scheme, we will use scheme (6.2) for the hyperbolic part coupled with the Crank-Nicolson discretization for the parabolic part, written as the third equation of system (3.6). The coupling in the hyperbolic system is made through the discretization of Φ x which is given in equation (3.5).

However, when we use the (AHO) 3 scheme (6.7) for the hyperbolic part, we also need to be more accurate in the discretization of the parabolic part. Thus, we will use a parabolic AHO scheme of order 4, namely a scheme which is generally of order two, but it is order four on stationary solutions. Notice that, in principle, it should possible to approximate the parabolic equation using the schemes proposed in [START_REF] Scharfetter | Large signal analysis of a silicon Read diode[END_REF], which share in some sense the spirit of the Well-Balanced schemes, namely they use an approximation of the exact local stationary solutions at the interface of the cells to integrate the different contributions of the schemes. However, in the present problem, it is difficult to isolate a specific stationary solution, since in general we could have an infinity of such solutions. So, to keep our scheme as simple as possible, we decided to apply our AHO approach also to the parabolic term.

Let us write more precisely the discretization for the equation φ t -D φ xx = au-bφ with the same notations used in Section 3. Using a Taylor expansion, as we did for instance at section 4, we obtain the following scheme

φ n+1 i = φ n i +D k h 2 (φ n i+1 -2φ n i +φ n i-1 )+ ak 12 (u n i+1 +10u n i +u n i-1 )- bk 12 (φ n i+1 +10φ n i +φ n i-1 ),
which has an order 4 truncation order when we insert the stationary solutions of the equation -D φ xx = au -bφ. Let us notice that the discretization is still done keeping our stencil on 3 points only. A crucial step in making our scheme higher order, is to use a fourth order discretization of Φ x on a 4 points stencil. We use here the following formula in (AHO) 3 :

Φ n x = 1 h 1 22 (-5φ n 3 + 28φ n 2 -23φ n 1 ), 1 132 (-11φ n 4 + 90φ n 3 -9φ n 2 -70φ n 1 ), • • • , 1 12 (-φ n j+2 + 8φ n j+1 -8φ n j-1 + φ n j-2 ) • • • , 1 132 (11φ n M -3 -90φ n M -2 + 9φ n M -1 + 70φ n M ), 1 22 (5φ n M -2 -28φ n M -1 + 23φ n M ) ,
with (Φ x ) n 0 = (Φ x ) n M +1 = 0 on the boundaries. 8.3. Comparison upwind vs AHO schemes for the chemotaxis system (1.1). In this section, we present a comparison between the standard upwind scheme presented at section 3 and the AHO schemes of section 6, coupled with the parabolic part using the ideas of the previous subsection. In problem (1.1), let us choose the coefficients a = b = 1 and the velocity λ = 10. In the following, we will always take at t = 0, φ = u and v = 0. The initial condition for u will be specified later on according to the different cases we consider. We will also always use as space and time steps h = 1/90 and k = 5 × 10 -4 unless differently specified. Let us finally notice that we do not impose a fixed final time to our simulations, which are stopped when the residues decrease below the value 10 -11 . That is why the final times mentioned below are all different.

Let us consider solutions with a mass equal to 1135, that is to say solutions after the first bifurcation point. In the case of a symmetric perturbation of the constant state, the asymptotic state for the function u should be the constant state. In that case, the upwind scheme (3.6) behaves very well and gives similar results to the more complex AHO schemes, see [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]. Now, in the case the asymptotic state is a nonconstant one, namely with in the case of an initial datum being a non-symmetric perturbation of constant state, we can see in Figure 1 that although all schemes give equivalent results for functions u and φ, the function v is not null for the upwind scheme, whereas it is for the AHO schemes. Functions are displayed at final time T = 10. Now, we consider solutions with a mass equal to 4100, that is to say solutions after the second bifurcation point. We begin with a symmetric perturbation of 4100. According to [START_REF] Schaaf | Stationary solutions of chemotaxis systems[END_REF] (see also [START_REF] Guarguaglini | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]), the asymptotic solutions u and φ should be symmetric non-constant functions with derivatives that vanish three times, two at the boundaries and one in the middle of the interval. The function φ should have a minimum equal to 4062 and a maximum equal to nearly 4140. In that case, the same phenomenon as before occurs, which is: upwind scheme gives good results, except for the function v. We also compare the two AHO schemes and we can see that the order two scheme gives better results for v than the order three one for a fixed space step, as shown in Figure 5. Finally, we consider a non-symmetric perturbation of the constant function equal to 4100. In that case, the asymptotic state for u is a function highly concentrated near the origin. We can see, on the left and in the middle of Figure 6, that the standard upwind scheme gives wrong results not only for the function v, but also for the function u for a time equal to T = 0.5, which is the time just before the stop in the computation due to the explosion of the numerical solution. However, AHO schemes converge and we can see, on the right of Figure 6, the asymptotic solutions obtained for the same initial data using both AHO schemes. We can also notice that an accurate result is obtained with boundary condition (5.1), whereas boundary condition (3.3) gives clearly a wrong computation. Fig. 6. The mass of the solution is equal to 4100 and the initial datum is a non-symmetric perturbation of constant state equal to 4100. We used as space and time steps h = 1/90, k = 5 × 10 -4 and time T = 0.5 (on the left and in the middle). On the left, the function u is displayed for some different schemes and in the middle, the function v. We notice that upwind scheme gives wrong results for the function v, but also for the function u, since upwind scheme explodes within a short time, i.e. the time at which solutions are displayed here. Asymptotic states are however reached with AHO schemes and are displayed on the right in red stars. Therefore the apparent lack of accuracy of the AHO schemes on the two first pictures is due to the fact that, at that time, we are still far from the asymptotic state. On the right, we can also notice the good behavior of boundary condition (5.1), unlike boundary condition (3.3).

8.4. Order of the AHO schemes in the case of the system (1.1) of chemotaxis. We investigate now the actual numerical order of our schemes. In that case, since it is difficult to compute an approximation of the exact asymptotic solution with a sufficient accuracy for small values of h, we plot the quantity ||u h -u h/2 || as a function of h, where u h is the approximation obtained with space step h. Let us notice that in that case the slope of the straight lines still gives the order of the considered schemes. In Figure 7, the error for the function u is shown on the left and the error for the function v on the right in L 1 norm. The initial condition is a non-symmetric perturbation of the constant state equal to 1135.

First, we can notice that for the function v, the standard upwind scheme is just order one, (AHO) 3 scheme is indeed of order three, but the (AHO) 2 scheme gives so small errors that the order computation is not reliable. These results perfectly correspond to the ones of Figure 4 for an hyperbolic system with an explicitly given source term. Besides, the results concerning the function u are fully coherent: upwind scheme is order one, (AHO) 2 scheme is order two and (AHO) 3 scheme is order three as expected thanks to the appropriate discretization of the parabolic part described in Subsection 8.2. The behavior of the approximation of φ has exactly the same behavior as the one for the function u. In Table 2, we display the slopes of the straight lines shown at Figure 7. 7 Upwind scheme (AHO) 2 scheme (AHO) 3

  Fig.4. Log-log plot of the error in L 1 norm as a function of the space step. On the left, the error for the function u is displayed for three different schemes, namely upwind, (AHO) 2 and (AHO)3 ; we can see on the right the error for v. The initial data are the same as in Figure3. We used different space steps h betwen 0.02 and 10 -3 , time steps equal to k = h/2λ and final convergent time T = 100.

Fig. 5 .

 5 Fig.5. The mass of the solution is equal to 4100 and the initial datum is a symmetric perturbation of constant state equal to 4100. On the left, the asymptotic state for u is displayed for different AHO schemes and on the right, the same for v. We used as space and time steps h = 1/90, k = 5 × 10 -4 and final time T = 20.

Fig. 7 .

 7 Fig. 7. Log-log plot of the error, namely the quantity ||u h -u h/2 || in L 1 norm as a function of the space step for three different schemes (upwind, (AHO) 2 , (AHO) 3 ). The error for the functions u and v are displayed on te left and on the right respectively. The initial condition for u and φ is a non-symmetric perturbation of the constant state equal to 1135. We took different space steps between 5 × 10 -4 and 10 -2 , time steps equal to k = h/2λ and convergent final time equal to T = 50.

  2.

	Conditions (4.6), (4.8) and (4.7a) are clearly satisfied. Conditions (4.7b) lead to
	k ≤	4h h + 4q

Table 1

 1

  Table of numerical orders of Figure

Table 2

 2 Table of numerical orders of Figure

Acknowledgments. This work has been partially supported by the project "Mathematical problems for the biological damage of monuments" in the CNR-CNRS 2008-2009 agreement, the INdAM-GNAMPA project 2008 "Hyperbolic models for


S Lacroix

G Gaval1

M Bouix

J J Pernelle

Differential response to the filamentous induction among various strains of Sphaerotilus natans

Background and aims:

The good efficiency of the activated sludge process is based on the gravity separation of purified water and bacterial biomass in the clarifier. The proliferation of filamentous bacteria preferentially occurring at low F/M ratio can be promoted by inadequate shortage of oxygen, by C / N imbalances or by the nature of the influent. Filamentous bulking caused by an overabundance of filamentous micro-organisms, severely disrupts the ability of sludge to settle and degrade the effluent quality. The most commonly employed curative treatments use oxidizing compounds such as chlorine. They are not selective. A targeted strategy could be considered when molecular mechanisms of the filamentous growth of the major strains involved will be better known.

Among the filamentous bacteria described in activated sludge, Sphaerotilus natans is especially favoured by shortage of oxygen, nitrogen or phosphorus, as well as by high concentrations of readily metabolizable soluble substrates [START_REF] Eikelboom | Process control of activated sludge plants by microscopic investigation[END_REF]. This filament is dominant in 12% of bulking cases in the United States [START_REF] Jenkins | Manual of the Causes and Control of Activated Sludge bulking, Foaming and other solids separation problems[END_REF]. It occurs mainly in waste water treatment plants (WWTP) treating industrial influents as from agro-food industry. This rod shape bacteria can grow as single cells when the media is rich. The filamentous form needs to synthesize a polysaccharide sheath to grow [START_REF] Suzuki | Identification of a gene essential for sheathed structure formation in Sphaerotilus natans, a filamentous sheated bacterium[END_REF]. This bacterium was chosen as test micro-organism because its filamentous form, can be induced reversibly from a culture of dispersed cells.

The emergence of a bulking could depend on the mode of filament elongation from a single cell pre-existing population. This study was first conducted to determine whether the S. natans filamentous development proceeds by linkage of pre-existing single cells or by cell divisions within proto filaments. The different easily degradable carbon substrates in the common culture media used to cultivate S. natans, do not show the same ability to induce filamentous growth. This ability has been characterized by experience plan in the purpose of defining an effective inducer media.

It was finally determined whether different strains of S. natans from collection, exhibited an equivalent capacity to grow as filaments. This last point will allow to select the most adequate strains on which filamentous growth mechanisms could be subsequently studied.

Methods and results

In order to determine the S. natans filament elongation mode, a culture of dispersed cells was transferred in filamentous inducing media [START_REF] Takeda | Isolation and Chemical Composition of the Sheath of Sphaerotilus natans[END_REF], to which was added an inhibitor of cell division: nalidixic acid 20 mg. L -1 . The cells were stained with cFDA / SE. After 48 hours of incubation at 30°C under agitation, the intensity of fluorescence per cell was identical to the original state, indicating that the cells have not entered into division. No filament structure was observed. The cells were locally clustered.

To determine their involvement in the morphological differentiation of S. natans, the main sources of organic carbon present in the culture media, an experiment plan involving 7 factors, based on Hadamard matrix was used in broth and solid agar plates (Larpent-Gourgaud and Sanglier, 1992). The 7 factors combined in this experiment correspond to the following: peptone, casitone, tryptone, glucose, glycerol, yeast extract and mineral salts. They were examined at two different concentrations on 3 different strains of S. natans: ATCC 13338T, 13929 and 15291. They were isolated respectively from river water, activated sludge and paper mill effluents. The 3 strains studied show differences in behaviour.

Nearly all tested agar media tend to promote dispersed growth form. In liquid broth, at low concentration, all carbon sources tested induce a filamentous growth. However, at high concentration, growth as single cells is favoured by the association of peptone and tryptone for strain 15291. This strain present the highest degree of differentiation. 13338T strain grows mainly as single cells and strain 13929 as filaments.

To clarify these behaviour differences, 5 different strains of S. natans have been studied by varying access to oxygen or carbon substrate concentration of nutrient medium. These strains were: 13338T, 13929, 15291, 29329 and 29330. The last two were isolated from pond. It appears that low concentrations of oxygen or substrate promote growth of filamentous form. However all strains show variable degrees of differentiation. Strains 29329 and 29330 grow mainly as single cells form very few filaments at low concentration of oxygen or substrate. Conversely, the strain 13929 grows mainly as filaments and high concentrations of substrate or oxygen are necessary for the emergence of single cells.

Conclusions:

The pre-existence of a large amount of single cells of bacterial species able to form filaments in activated sludge at the onset of inducing factors, represents a potential risk to the WWTP. This biomass is currently very difficult to detect. It is shown here that S. natans dispersed single cells, could cause a filamentous proliferation by cell divisions and not by simple cell linkage. In pure culture, this morphological duality expressed preferentially with peptones and tryptone as carbon substrates. The ATCC 15291 strain presents broader dynamics of morphological differentiation, closely followed by the strain 13338T. The strains 29329 and 29330, growing mainly as single cells and strain 13929 growing mainly as filaments, are less suited for a study on filamentous growing mechanisms. For instance, a comparative study of proteomes requires cultures of the same strain in both morphological states. These studies should be complemented by using strains of filamentous bacteria frequently encountered in bulking. Unfortunately, the latter are generally much more difficult to grow.