Pierre Degond

Fabrice Deluzet

Claudia Negulescu

An asymptotic preserving scheme for strongly anisotropic elliptic problems

come

An asymptotic preserving scheme for strongly anisotropic elliptic problems

Introduction

The objective of this paper is to introduce an efficient and accurate numerical scheme to solve a strongly anisotropic elliptic problem of the form

-∇ • (A∇φ) = f , in Ω φ = 0 on ∂Ω D , ∂ z φ = 0 on ∂Ω z , (1)
where Ω ⊂ R 2 or Ω ⊂ R 3 is a domain, with boundary ∂Ω = ∂Ω D ∪ ∂Ω z and the diffusion matrix A is given by

A = A ⊥ 0 0 1 ε A z .
The terms A ⊥ and A z are of the same order of magnitude, whereas the parameter 0 < ε < 1 can be very small, provoking thus the high anisotropy of the problem. In the present paper the considered anisotropy direction is fixed and is aligned with the z-axis of a Cartesian coordinate system. The method presented here is extended in some forthcoming works to more general anisotropies [START_REF] Brull | An asymptotic preserving scheme in the drift limit for the Euler-Lorentz system for a variable magnetic field[END_REF]. Anisotropic problems are common in mathematical modeling and numerical simulation. Indeed they occur in several fields of applications such as flows in porous media [START_REF] Ashby | A numerical solution of groundwater flow and contaminant transport on the CRAY T3D and C90 supercomputers[END_REF][START_REF] Hou | A Multiscale Finite Element Method for Elliptic Problems in Composite Materials and Porous Media[END_REF], semiconductor modelling [START_REF] Manku | Electrical properties of silicon under nonuniform stress[END_REF], quasi-neutral plasma simulations [START_REF] Degond | Asymptotic-Preserving Particle-In-Cell method for the Vlasov-Poisson system near quasineutrality[END_REF], image processing [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF][START_REF] Wang | Anisotropic diffusion with nonlinear structure tensor[END_REF], atmospheric or oceanic flows [START_REF] Tréguier | Modélisation numérique pour l'océanographie physique[END_REF], and so on, the list being not exhaustive. More specifically high anisotropy aligned with one direction may occur in shell problems or simulation in stretched media. The initial motivation for the present work is closely related to magnetized plasma simulations such as atmospheric [START_REF] Kelley | Mid-Latitude ionospheric fluctuation spectra due to secondary E ×B instabilities[END_REF][START_REF] Keskinen | Three-dimensional nonlinear evolution of equatorial ionospheric spread-F bubbles[END_REF] or inertial fusion plasmas [START_REF] Beer | Field-aligned coordinates for nonlinear simulations of tokamak turbulence[END_REF][START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF] or plasma thrusters [START_REF] Adam | Physics, simulation, and diagnostics of Hall effect thrusters[END_REF]. In this context, the medium is structured by the magnetic field. Indeed, the motion of charged particles in planes perpendicular to the magnetic field is governed by a fast gyration around the magnetic field lines. This explains the large number of collisions the particles encounter in the perpendicular plane, whereas the dynamic in the parallel direction is rather undisturbed. As a consequence the particle mobilities in the perpendicular and parallel directions differ by many orders of magnitude. In the context of ionospheric plasma modelling [START_REF] Besse | Numerical simulations of the ionospheric striation model in a non-uniform magnetic field[END_REF][START_REF] Hysell | An overview and synthesis of plasma irregularities in equatorial spread F[END_REF], the ratio of the aligned and transverse mobilities (denoted in this paper by ε -1 ) can be as huge as ten to the power ten. The relevant boundary conditions in many fields of application are periodic (for instance in simulations of tokamak plasmas on a torus) or Neumann boundary conditions (see for instance [START_REF] Besse | A model Hierarchy for Ionospheric Plasma modeling[END_REF] for atmospheric plasmas). The system [START_REF] Adam | Physics, simulation, and diagnostics of Hall effect thrusters[END_REF] is thus a good model to elaborate a robust numerical method.

The main difficulties with the resolution of problem (1) are of numerical nature, as solving this singular perturbation problem for small 0 < ε ≪ 1 is rather delicate. Indeed, replacing in the anisotropic elliptic equation ε by zero, yields an ill-posed problem, which has an infinite number of solutions (namely all functions which are constant in the z-direction). This feature is translated in the discrete case (after the discretization of the problem) into a linear system which is very ill-conditioned for ε ≪ 1, due to the different order of magnitudes of the various terms. As a consequence standard numerical methods for the resolution of linear systems lead to important numerical costs and unacceptable numerical errors. More generally, this numerical difficulty arises when the boundary conditions supplied to the dominant O(1/ε) operator lead to an ill-posed problem with a multiplicity of solutions. This is the case for Neumann boundary conditions, but also of periodic boundary conditions. If instead, the boundary conditions are such that the dominant operator gives a well-posed problem with a unique solution, this difficulty vanishes as the leading operator alone will suffice to completely determine the limit solution. In this case, one can resort to standard methods. This is the case of Dirichlet or Robin boundary conditions. In spite of the fact that the problem addressed in the present paper arises only with specific boundary conditions, it has a considerable impact in many physics problem, such as plasmas, geophysical flows, plate and shells, etc. In this paper, we will focus on Neumann boundary conditions because they represent a larger range of physical applications, but we could address periodic boundary conditions in a similar way.

Numerical methods for anisotropic elliptic problems have been extensively investigated in the literature. Depending on the underlying physics, distinct numerical methods are developed. For example domain decomposition (Schur complement) and multigrid techniques, using multiple coarse grid corrections are adapted to anisotropic equations in [START_REF] Giraud | Schur complement preconditioners for anisotropic problems[END_REF][START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF] and [START_REF] Gee | A new smoothed aggregation multigrid method for anisotropic problems[END_REF][START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF]. For anisotropy aligned with one (or two directions), point (or plane) smoothers are shown to be very efficient [START_REF] Llorente | Robust multigrid smoothers for three dimensional elliptic equations with strong anisotropies[END_REF]. A problem very similar to (1) is addressed in [START_REF] Guillaume | Numerical Convergence of a Parameterisation Method for the Solution of a Highly Anisotropic Two-Dimensional Elliptic Problem[END_REF], treated via a parametrisation technique, and seems to give good results for rather large anisotropy ratios. However, these techniques are only developed in the context of an elliptic operator with a dominant part supplemented with Dirichlet boundary conditions. An alternative approach for dealing with highly anisotropic problems is based on a mathematical reformulation of the continuous problem, in order to obtain a more harmless problem, which can be solved numerically in an uncomplicated manner. In this category can be situated for example asymptotic models, describing for small values of the asymptotic parameter ε the evolution of an approximation φ of the solution of (1) [START_REF] Besse | A model Hierarchy for Ionospheric Plasma modeling[END_REF][START_REF] Keskinen | Nonlinear theory of the E × B instability with an inhomogeneous electric field[END_REF]. However, these asymptotic models are precise only for ε ≪ 1, and cannot be used on the whole range of values covered by the physical parameter ε. Thus model coupling methods have to be employed. In sub-domains where the limit model is no longer valid, the original model has to be used, which means that a model coupling strategy has to be developed. However the coupling strategy requires the existence of an area where both models are valid and still demands an accurate numerical method for the resolution of the original model (i.e. the anisotropic elliptic problem) with large anisotropies. This can be rather undesirable.

In this paper, we present an original numerical algorithm belonging to the second approach. A reformulation of the continuous problem (1) will permit us to solve this problem in an inexpensive way and accurately enough, independently of the parameter ε. This scheme is related to the Asymptotic Preserving numerical method introduced in [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF]. These techniques are designed to provide computations in various regimes without any restriction on the discretization meshes and with the additional property to converge towards the solution of the limit problem when the asymptotic parameter goes to zero. The derivation of such Asymptotic Preserving methods requires first the identification of the limit model. For singular perturbation problems, a reformulation of the problem is required in order to derive a set of equations containing both the initial and the limit model with a continuous transition from one regime to another, according to the values of the parameter ε. This reformulated system of equations sets the foundation of the AP-scheme. Other singular perturbations have already been explored in previous studies, for instance quasi-neutral or gyro-fluid limits [START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF]. These techniques have been first introduced for non-stationary systems of equations, for which the time discretization must be studied with care in order to guarantee the asymptotic preserving property. For the anisotropic elliptic equation investigated in this article, we only need to precise the reformulated system and provide a discretization of this one.

The outline of this paper is the following. Section 2 of this article presents first the initial anisotropic elliptic model. In the remainder of this paper, it will be referred to as the Singular-Perturbation model (P-model). The reformulated system (referred to as the Asymptotic Preserving formulation or AP-formulation) is then derived. It relates on a decomposition of the solution φ(x, z) according to its mean part φ(x) along the z coordinate and a fluctuation φ ′ (x, z) consisting of a correction to the mean part needed to recover the full solution. The mean part φ(x) is solution of an ε-independent elliptic problem, and the fluctuation φ ′ (x, z) = φ(x, z) -φ(x) is given by a well-posed ε-dependent elliptic problem. The advantage is that the ε-dependent problem for the fluctuation is well-posed and solvable in an inexpensive way, and this uniformly in ε. In the limit ε → 0 the AP-formulation reduces to the so called Limit model (L-model), whose solution is an acceptable approximation of the P-model solution for ε ≪ 1. The present derivation is carried out in the framework of an anisotropy aligned along one axis of a Cartesian coordinate system. In the context of magnetized plasma simulations, this initial work is extended in a forthcoming work for the three dimensional case in curvilinear coordinates, designed to fit a more complex magnetic field topology (i.e. anisotropy direction) [START_REF] Besse | Numerical simulations of the ionospheric striation model in a non-uniform magnetic field[END_REF]. The main constraints of this method reside in the construction of the mean part which necessitates the integration of the solution along the anisotropy direction. This operation is easily carried out in the context of coordinates adapted with the anisotropy direction. However, an extension of the techniques presented here is currently developed for non-adapted coordinates [START_REF] Brull | An asymptotic preserving scheme in the drift limit for the Euler-Lorentz system for a variable magnetic field[END_REF]. Section 3 is devoted to the numerical implementation of the AP-formulation. Numerical results are then presented for a test case, and the three approaches (AP-formulation, straight discretization and resolution of the P-model and L-model) are compared according to the precision of the approximation for different values of ε. In section 4 we shall rigorously analyse the convergence of the AP-scheme. Error estimates will be established which underline the advantages of the AP-scheme as compared to the initial Singular Perturbation model and the Limit model.

Current research directions are concerned with the adaptation of the present technique to the case of arbitrary spatially varying anisotropies, without adaptation of the coordinate system to the direction of the anisotropy. These developments will allow the treatment on nonlinear problems, when the diffusion tensor (and its principal directions) depend on the solution itself. This treatment will involve iterative methods which, at each iterate, will reduce the problem to the solution of a linear anisotropic diffusion problem.

The asymptotic preserving formulation

For simplicity we shall consider in this paper the two-dimensional problem, posed on a rectangular domain Ω = Ω x × Ω z , where Ω x ⊂ R and Ω z ⊂ R are intervals. The ideas exposed here can be extended without any problems to the more physical three-dimensional domain, with two transverse directions (x, y) and an anisotropy direction aligned with the z-direction. In this section we introduce the Singular Perturbation Model, the Limit Model and the Asymptotic Preserving formulation.

The Singular Perturbation Model (P-model)

The main concern of this paper is the numerical resolution of the following anisotropic, elliptic problem, called in the sequel Singular Perturbation Model (P )

     -∇ • (A∇φ) = f , in Ω , ∂φ ∂z = 0 on Ω x × ∂Ω z , φ = 0 on ∂Ω x × Ω z . (2)
The anisotropy of the media is modeled via the definition of the diffusion matrix

A A = A ⊥ 0 0 1 ε A z , (3)
where A ⊥ (x, z) and A z (x, z) are given functions with comparable order of magnitudes. The source term f (x, z) is given and the parameter ε is small compared to both A ⊥ as well as A z . The medium becomes more anisotropic as the value of ε goes to zero.

The limit regime (L-model)

In this section we establish that in the limit ε → 0 the solution of the perturbation model converges towards φ, solution of the L-model defined by

(L)    - ∂ ∂x Ā⊥ ∂ φ ∂x = f (x) , in Ω x , φ = 0 on ∂Ω x , (4)
where overlined quantities designate averages over the z-coordinate :

f (x) = 1 |Ω z | Ωz f (x, z) dz.
First we can rewrite the P-model as

(P )        - ∂ ∂x A ⊥ ∂φ ∂x - 1 ε ∂ ∂z A z ∂φ ∂z = f , in Ω , ∂φ ∂z = 0 on Ω x × ∂Ω z , φ = 0 on ∂Ω x × Ω z , (5)
and integrating along the z-coordinate gives

∂ ∂x A ⊥ ∂φ ∂x = f (x) . (6)
This equation holds for any ε > 0. Now, letting formally ε tend to zero in (5) yields the reduced model (R-model) (R)

       - ∂ ∂z A z ∂φ ∂z = 0 , in Ω , ∂φ ∂z = 0 on Ω x × ∂Ω z , φ = 0 on ∂Ω x × Ω z . (7)
The functions verifying this ill-posed R-model are constant along the z-coordinate. Thus including this asymptotic limit property into equation [START_REF] Besse | Numerical simulations of the ionospheric striation model in a non-uniform magnetic field[END_REF] gives rise to the L-model ( 4), verified by the solution of the Singular Perturbation model in the limit ε → 0.

Remark 2.1 The L-model is the singular limit of the original P-model [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF]. It provides an accurate approximation of the P-solution only for small values of ε. The P-model is valid for all 0 < ε < 1, but numerically impracticable for ε ≪ 1. Indeed working with a finite precision, the asymptotic model degenerates into the R-model defined by (7) as ε vanishes. This R-model is ill-posed since it exhibits an infinite amount of solutions φ = φ(x), depending only on the variable x. This implies that the discretization matrix derived from the P-model is very ill-conditioned for small 0 < ε ≪ 1. This point is addressed by the numerical experiments of section 3.2. Consequently, in a domain where ε varies significantly, a model coupling method has to be developed in order to exploit the validity of each model, the P-and L-model. This can be rather undesirable. In the next section we shall present an alternative approach, which is based on a reformulation of the Singular-Perturbation model providing a means of computing an accurate numerical approximation of the solution for all values 0 < ε < 1.

Remark 2.2

The asymptotics is totally different in the case of Dirichlet boundary conditions. In this case, the R-model is well posed, with a unique solution, and there is no difficulty anymore. Any standard numerical solution of the P-model will converge to that of the R-model. In other words, with Dirichlet boundary conditions, the perturbation becomes regular and the limit solution is fully determined by the formal limit system. The situation and the difficulty addressed in the present paper require that the R-model be ill-posed. This is the case with Neumann boundary conditions (which is the framework chosen here) but also with periodic boundary conditions, or any other boundary condition which would result in an ill-posed R-model.

The Asymptotic Preserving reformulation (AP-formulation)

In order to circumvent the just described numerical difficulties in handling the Singular Perturbation model, we introduce a reformulation, which permits a transition from the initial P -model to its singular limit (L-model), as ε → 0.

For this, we shall decompose each quantity f (x, z) into its mean value f (x) along the z coordinate and a fluctuation part f ′ (x, z). For simplicity reasons let in the following Ω x := (0, L x ) and Ω z := (0, L z ). Then

f (x, z) = f (x) + f ′ (x, z) , (8)
with

f (x) := 1 L z Lz 0 f (x, z)dz , f ′ (x, z) := f (x, z) -f (x) . (9)
Note that we have the following properties

f ′ = 0 , ( ∂f /∂x) = ∂f /∂x , f g = f ḡ + f ′ g ′ , ( 10
) ∂f /∂z = ∂f ′ /∂z , ( ∂f /∂x) ′ = ∂ f ′ /∂x , (f g) ′ = f ′ g ′ -f ′ g ′ + f g ′ + f ′ ḡ . ( 11
)
Taking now the mean of the elliptic equation ( 5) along the z-coordinate, we get thanks to ( 10) and ( 11), an equation for the evolution of the mean part φ(x)

(AP 1)      - ∂ ∂x Ā⊥ ∂ φ ∂x = f + ∂ ∂x A ′ ⊥ ∂φ ′ ∂x , in Ω x , φ = 0 on ∂Ω x . (12)
Substracting from [START_REF] Besse | A model Hierarchy for Ionospheric Plasma modeling[END_REF] this mean equation ( 12), gives rise to the evolution equation for the fluctuation part φ ′ (x, z)

(AP 2)                      - ∂ ∂z A z ∂φ ′ ∂z -ε ∂ ∂x A ⊥ ∂φ ′ ∂x + ε ∂ ∂x A ′ ⊥ ∂φ ′ ∂x = εf ′ + ε ∂ ∂x A ′ ⊥ ∂φ ∂x , in Ω , ∂φ ′ ∂z = 0 on Ω x × ∂Ω z , φ ′ = 0 on ∂Ω x × Ω z , φ ′ = 0 , in Ω x . (13)
Thus we have replaced the resolution of the initial Singular Perturbation model ( 5) by the resolution of the system ( 12)-( 13), which will be done iteratively. Starting from a guess function φ ′ , equation [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF] gives the mean value φ(x), which inserted in ( 13) shall give the fluctuation part φ ′ (x, z) and so on.

The constraint φ ′ = 0 in (13) (which is automatic for ε > 0, as explained in Remark 2.3) has the essential consequence that the conditioning of the discretized system becomes ε-independent, because the problem (13) reduces in the limit ε → 0 to the system

             - ∂ ∂z A z ∂φ ′ ∂z = 0 , in Ω , ∂φ ′ ∂z = 0 on Ω x × ∂Ω z , φ ′ = 0 on ∂Ω x × Ω z , φ′ = 0 in Ω x , (14)
which is uniquely solvable, with the solution φ ′ ≡ 0. Inserting this solution in [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF], we conclude that the solution of the AP formulation converges for ε → 0 towards the mean value part φ(x), computed thanks to the Limit model

(L)    - ∂ ∂x Ā⊥ ∂ φ ∂x = f (x) , in Ω x , φ = 0 on ∂Ω x . (15)
The AP reformulation ( 12)-( 13) is equivalent to the Singular Perturbation problem (5) and is therefore valid for all 0 < ε < 1. This new formulation guarantees that, working with a finite precision arithmetic, the computed solution converges in the limit ε → 0 towards the solution of the limit model ( 4). This is a huge difference with the original Singular Perturbation model which degenerates into an ill-posed problem. Thus, by using the APformulation, we expect the computation of the numerical solution to be accurate, uniformly in ε.

For the detailed mathematical proofs, we refer to the next section.

Remark 2.3

The condition φ ′ = 0 in (13) holds automatically for ε > 0, since the righthand side has zero average along the z-coordinate. Indeed, let ψ be the solution of

       - ∂ ∂z A z ∂ψ ∂z -ε ∂ ∂x A ⊥ ∂ψ ∂x + ε ∂ ∂x A ′ ⊥ ∂ψ ∂x = εg ′ , in Ω , ∂ψ ∂z = 0 on Ω x × ∂Ω z , ψ = 0 on ∂Ω x × Ω z , (16)
with g ′ = 0. Taking the average along z, we get

   - ∂ ∂x Ā⊥ ∂ ψ ∂x = 0 , in Ω x , ψ = 0 on ∂Ω x ,
and thus ψ ≡ 0, which is nothing but the constraint added in [START_REF] Gee | A new smoothed aggregation multigrid method for anisotropic problems[END_REF].

The computations of the fluctuating part φ ′ via the equation (13) requires the discretization of an integro-differential operator. This means that the discretization matrix will contain dense blocks. However, using (12) the system (AP2) can be rewritten as

(AP 2 ′ )                      - ∂ ∂z A z ∂φ ′ ∂z -ε ∂ ∂x A ⊥ ∂φ ′ ∂x = εf + ε ∂ ∂x A ⊥ ∂φ ∂x , in Ω , ∂φ ′ ∂z = 0 on Ω x × ∂Ω z , φ ′ = 0 on ∂Ω x × Ω z , φ ′ = 0 , in Ω x . (17)
In this expression the right-hand side has no longer zero mean value along the z-coordinate, but the integro-differential operator has disappeared. The associated discretization matrix is thus sparser than that obtained from the system [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF]. Systems ( 12)-( 13) and ( 12)-( 17) are equivalent.

Mathematical study of the AP-formulation

We establish in this section the mathematical framework of the AP-formulation ( 12)-( 13) and study its mathematical properties. Let us thus introduce the two Hilbert-spaces

V := {ψ(•, •) ∈ H 1 (Ω) / ψ = 0 on ∂Ω x × Ω z } , W := {ψ(•) ∈ H 1 (Ω x ) / ψ = 0 on ∂Ω x } ,
with the corresponding scalar-products

(φ, ψ) V := ε(∂ x φ, ∂ x ψ) L 2 + (∂ z φ, ∂ z ψ) L 2 , (φ, ψ) W := (∂ x φ, ∂ x ψ) L 2 , ( 18
)
and the induced norms || • || V , respectively || • || W . For simplicity reasons, we denote in the sequel the L 2 scalar-product simply by the bracket (•, •). Defining the following bilinear forms

a 0 (φ ′ , ψ ′ ) := Lz 0 Lx 0 A z (x, z) ∂φ ′ ∂z (x, z) ∂ψ ′ ∂z (x, z)dxdz , a 1 (φ ′ , ψ ′ ) := Lz 0 Lx 0 A ⊥ (x, z) ∂φ ′ ∂x (x, z) ∂ψ ′ ∂x (x, z)dxdz , a 2 φ, ψ := Lx 0 Ā⊥ (x) ∂ φ ∂x (x) ∂ ψ ∂x (x)dx , c φ ′ , ψ := Lz 0 Lx 0 A ′ ⊥ (x, z) ∂φ ′ ∂x (x, z) ∂ ψ ∂x (x)dxdz , d(φ ′ , ψ ′ ) := 1 L z Lx 0 Lz 0 Lz 0 A ′ ⊥ (x, z) ∂φ ′ ∂x (x, z) ∂ψ ′ ∂x (x, ζ) dzdζdx , b( P , ψ ′ ) := Lx 0 P (x) Lz 0 ψ ′ (x, z)dzdx , a(φ ′ , ψ ′ ) := a 0 (φ ′ , ψ ′ ) + εa 1 (φ ′ , ψ ′ ) -εd(φ ′ , ψ ′ ) , (19)
permits to rewrite the AP system ( 12)-( 13) under the weak form

(AP )          a 2 φ, ψ = ( f , ψ) - 1 L z c φ ′ , ψ , ∀ ψ ∈ W , a (φ ′ , ψ ′ ) + b( P , ψ ′ ) = ε(f ′ , ψ ′ ) -εc ψ ′ , φ , ∀ψ ′ ∈ V , b( Q, φ ′ ) = 0 , ∀ Q ∈ W , (20)
where φ ′ (x, z) ∈ V, φ(x) ∈ W as well as P (x) ∈ W are the unknowns and ψ ′ ∈ V, ψ ∈ W and Q ∈ W the test functions. It can be observed that the constraint φ′ = 0 was introduced via the Lagrange multiplier P . We will see in the next theorem that the weak formulation (20) is equivalent for ε > 0 to the system

     a 2 φ, ψ = ( f , ψ) - 1 L z c φ ′ , ψ , ∀ ψ ∈ W , (21)
a (φ ′ , ψ ′ ) = ε(f ′ , ψ ′ ) -εc ψ ′ , φ , ∀ψ ′ ∈ V , (22)
where the explicit constraint φ′ = 0 does not appear. Let us assume in the sequel

Hypothesis A Let the diffusion functions A ⊥ ∈ L ∞ (Ω) and A z ∈ L ∞ (Ω) satisfy 0 < c ⊥ ≤ A ⊥ (x, z) ≤ M ⊥ , 0 < c z ≤ A z (x, z) ≤ M z , f.a.a. (x, z) ∈ Ω , with some positive constants c ⊥ , c z , M ⊥ , M z . Let moreover f ∈ L 2 (Ω).
The next theorem analyzes the well-posedness of the AP-formulation.

Theorem 2.4 For every ε > 0 the problem ( 21)-( 22) admits under Hypothesis A a unique solution (φ ′ ε , φ ε ) ∈ V × W, where φ ε := φ ′ ε + φ ε is the unique solution of the Singular Perturbation model [START_REF] Besse | A model Hierarchy for Ionospheric Plasma modeling[END_REF]. The function φ ′ ε has zero mean value along the z-coordinate, i.e. φ ′ ε = 0 for every ε > 0. Consequently, (φ ′ ε , φ ε ) ∈ V×W is the unique solution of ( 21)-( 22) if and only if [START_REF] Keskinen | Nonlinear theory of the E × B instability with an inhomogeneous electric field[END_REF]. In this last case, we have P ε = 0. Finally, these solutions satisfy the bounds

(φ ′ ε , φ ε , P ε ) ∈ V × W × W is a solution of the AP-formulation
||φ ε || H 1 (Ω) ≤ C||f || L 2 (Ω) , ||φ ′ ε || H 1 (Ω) ≤ C||f || L 2 (Ω) , ||φ ε || H 1 (Ωx) ≤ C||f || L 2 (Ω) ,
with an ε-independent constant C > 0. In the limit ε → 0 there exist some functions (φ ′ 0 , φ 0 ) ∈ V × W, such that we have the following weak convergences in

H 1 φ ′ ε ⇀ ε→0 φ ′ 0 in H 1 (Ω) , φ ε ⇀ ε→0 φ 0 in H 1 (Ω x ) ,
and the strong L 2 convergences

φ ′ ε → ε→0 φ ′ 0 in L 2 (Ω) , ∂ z φ ′ ε → ε→0 ∂ z φ ′ 0 in L 2 (Ω) , φ ε → ε→0 φ 0 in L 2 (Ω x ) ,
where φ ′ 0 ≡ 0 and φ 0 is the unique solution of the Limit model ( 4).

Proof: The Singular Perturbation model ( 5) and the Limit model ( 4) are standard elliptic problems and posses under Hypothesis A (and for every ε > 0) unique solutions φ ε ∈ V, respectively φ ∈ W. It is then a simple consequence of the decomposition [START_REF] Brull | An asymptotic preserving scheme in the drift limit for the Euler-Lorentz system for a variable magnetic field[END_REF], that the problem ( 21)-( 22) admits a unique solution (φ ′ ε , φ ε ) ∈ V × W, where φ ε (x) := 1 Lz Lz 0 φ ε (x, z)dz is the mean and φ ′ ε := φ εφ ε the fluctuation part. Thus we have also φ ′ ε = 0. This property can also be understood from the fact that the right-hand side of [START_REF] Gee | A new smoothed aggregation multigrid method for anisotropic problems[END_REF], denoted in the sequel by g g(x, z)

:= f ′ (x, z) + ∂ ∂x A ′ ⊥ (x, z) ∂φ ∂x (x) ,
has zero mean value along the z-coordinate. Indeed, taking in [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF] test functions ψ ′ (x) ∈ V depending only on x, yields immediately that φ ′ ε = 0 for all ε > 0. Standard stability results for elliptic problems yield now the ε-independent estimate for the solution of the Singular Perturbation model ( 5)

||φ ε || 2 H 1 (Ω) ≤ ||∂ x φ ε || 2 L 2 (Ω) + 1 ε ||∂ z φ ε || 2 L 2 (Ω) ≤ C||f || 2 L 2 (Ω) , implying that ||φ ε || 2 H 1 (Ωx) ≤ C||f || 2 L 2 (Ω) and ||φ ′ ε || 2 H 1 (Ω) ≤ C||f || 2 L 2 (Ω)
, with a constant C > 0 independent of ε > 0. Thus there exist some functions (φ ′ 0 , φ 0 ) ∈ V × W, such that, up to a subsequence φ ′ ε ⇀ ε→0 φ ′ 0 in H 1 (Ω) and φ ε ⇀ ε→0 φ 0 in H 1 (Ω x ). Hence we have

Lx 0 Lz 0 φ ′ ε (x, z)ψ(x, z)dx dz → ε→0 Lx 0 Lz 0 φ ′ 0 (x, z)ψ(x, z)dx dz , ∀ψ ∈ V .
Taking here ψ(x) ∈ V depending only on the x-coordinate, we observe that the feature φ ′ ε ≡ 0 yields the crucial property of the limit solution φ ′ 0 ≡ 0. Passing now to the limit ε → 0 in [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF], we get that φ ′ 0 is solution of

a 0 (φ ′ 0 , ψ ′ ) = 0 , ∀ψ ′ ∈ V , with φ ′ 0 = 0 in Ω x ,
which is the weak form of ( 14) and implies φ ′ 0 ≡ 0. Finally, passing to the limit in [START_REF] Keskinen | Three-dimensional nonlinear evolution of equatorial ionospheric spread-F bubbles[END_REF], yields that φ 0 is the unique solution of the Limit model [START_REF] Anderson | LAPACK Users' Guide, Third Edition[END_REF]. Because of the uniqueness of the limit (φ ′ 0 , φ 0 ), we deduce that the whole sequence (φ ′ ε , φ ε ) converges weakly towards this limit. To conclude the first part of the proof, we shall show the strong L 2 convergences. For this, taking in (22) φ ′ ε as test function and passing to the limit ε → 0, yields

∂ z φ ′ ε → 0 in L 2 (Ω). As φ ′ ε ∈ V and φ′ ε = 0, the Poincaré inequality ||φ ′ ε || L 2 ≤ C||∂ z φ ′ ε || L 2 ,
is valid and implies that φ ′ ε → 0 in L 2 (Ω). The convergence φε → φ0 in L 2 (Ω x ) is immediate by compacity. It remains finally to prove the equivalence between ( 20) and ( 21)- [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF]. This is immediate. Indeed, if (φ ′ ε , φ ε ) ∈ V × W is solution of ( 21)- [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF], then (φ ′ ε , φ ε , 0) is solution of [START_REF] Keskinen | Nonlinear theory of the E × B instability with an inhomogeneous electric field[END_REF]. And if (φ ′ ε , φ ε , P ε ) ∈ V × W × W satisfies (20), then P ε ≡ 0 (obvious by taking as test function in (20) ψ ′ (x) ∈ V depending only on x) and (φ ′ ε , φ ε ) solves hence ( 21)- [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF].

The subject of the next section will be the numerical resolution of the AP-formulation ( 12)-( 13) (or [START_REF] Keskinen | Nonlinear theory of the E × B instability with an inhomogeneous electric field[END_REF]) and this shall be done iteratively via a fixed-point application. Let us thus introduce here the fixed-point map, construct an iterative sequence and analyze its convergence. In the rest of this section, the parameter ε > 0 shall be considered as fixed.

Due to the fact that the two systems ( 20) and ( 21)-( 22) are equivalent, we shall concentrate on the simpler one, i.e. ( 21)- [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF]. Let us define the Hilbert space

U := {ψ(•, •) ∈ V / ψ = 0} ,
associated with the scalar product

(φ, ψ) * := Lx 0 Lz 0 A z ∂ z φ ∂ z ψdzdx + ε Lx 0 Lz 0 A ⊥ ∂ x φ ∂ x ψdzdx ,
which is equivalent to the scalar product (•, •) V on V, defined by [START_REF] Kelley | Mid-Latitude ionospheric fluctuation spectra due to secondary E ×B instabilities[END_REF]. The fixed-point map T : U → U is defined as follows: With φ ′ ∈ U we associate φ ∈ W, solution of [START_REF] Keskinen | Three-dimensional nonlinear evolution of equatorial ionospheric spread-F bubbles[END_REF]. Then constructing the right-hand side of ( 22) via this φ ∈ W, we define T (φ ′ ) as the corresponding solution of [START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF]. Denoting by (φ ′ * , φ * ) ∈ V × W the unique solution of ( 21)-( 22), we remark by Theorem 2.4 that φ ′ * ∈ U and that it is the unique fixed-point of the map T .

Theorem 2.5 Let ε > 0 be fixed and let φ ′ * ∈ U be the unique fixed-point of the application T : U → U constructed as follows

φ ′ ∈ U (21) --→ φ ∈ W (22) --→ T (φ ′ ) ∈ U .
Then for every starting point φ ′ 0 ∈ U, the sequence Then the set of fixed-points of T, denoted by F P (T ), is identical with the set of accumulation points of the sequences {T k (φ)} k∈N , with φ ∈ U, set which is denoted by AP (T ). Moreover, these two spaces contain at most one element.

φ ′ k := T (φ ′ k-1 ) = T k (φ ′ 0 ) converges in (U, || • || * ),
Proof of theorem 2.5 :

The linear application T is well-defined. The first step φ ′ ∈ U

--→ φ ∈ W is immediate by the Lax-Milgram theorem. For the second step, we remark that for given φ ∈ W the equation

a(θ, ψ ′ ) = ε(f ′ , ψ ′ ) -εc(ψ ′ , φ) , ∀ψ ′ ∈ V , (23)
has a unique solution θ ∈ U. Indeed, we notice first (by taking test functions only depending on the x-coordinate) that θ = 0. This enables us to consider instead of ( 23), the variational formulation

m(θ, ψ ′ ) = ε(f ′ , ψ ′ ) -εc(ψ ′ , φ) , ∀ψ ′ ∈ V , (24)
where the bilinear form a(•, •), which is not coercive, was replaced by the coercive bilinear form m(•, •), given by

m(θ, ψ ′ ) := a(θ, ψ ′ ) + εM ⊥ L z Lx 0 Lz 0 ∂ x θ(x, z)dz Lz 0 ∂ x ψ ′ (x, z)dz dx . (25)
Indeed, due to the property θ = 0, the two equations ( 23) and ( 24) are equivalent and this time m(•, •) is a continuous, coercive bilinear form, as for all ψ ′ ∈ V we have

m(ψ ′ , ψ ′ ) ≥ Lx 0 Lz 0 A z |∂ z ψ ′ | 2 dzdx + ε Lx 0 Lz 0 A ⊥ |∂ x ψ ′ | 2 dzdx ≥ C||ψ ′ || 2 V .
Thus the Lax-Milgram theorem implies the existence and uniqueness of a solution θ ∈ U of the continuous problem [START_REF] Manku | Electrical properties of silicon under nonuniform stress[END_REF] and hence also of problem [START_REF] Llorente | Robust multigrid smoothers for three dimensional elliptic equations with strong anisotropies[END_REF]. We have shown by this that T is a well-defined mapping. Furthermore we know that T admits, for fixed ε > 0, a unique fixed-point, denoted by φ ′ * ∈ U. Let us now suppose that we have shown that T is contractive. Then lemma 2.6 implies that F P (T ) = AP (T ) = {φ ′ * }. Thus choosing an arbitrary starting point φ ′ 0 ∈ U, and constructing the sequence

φ ′ k := T (φ ′ k-1 ) = T k (φ ′ 0 ),
:= φ ′ 1 -φ ′ 2 , φ := φ 1 -φ 2 (
where φ i ∈ W are the corresponding solutions of ( 21)) and θ ′ := T (φ ′ 1 ) -T (φ ′ 2 ), we have to show that ||θ ′ || * < ||φ ′ || * . First we observe that φ solves

a 2 ( φ, ψ) = - 1 L z c(φ ′ , ψ) , ∀ ψ ∈ W , (26)
and θ ′ is solution of

a(θ ′ , ψ ′ ) = -εc ψ ′ , φ , ∀ψ ′ ∈ V . (27)
Taking in (26) φ as test function, gives rise to

Lx 0 A ⊥ |∂ x φ(x)| 2 dx = - Lx 0 1 L z Lz 0 A ′ ⊥ ∂ x φ ′ (x, z)dz ∂ x φ(x) dx = - Lx 0 1 L z Lz 0 A ⊥ ∂ x φ ′ (x, z)dz ∂ x φ(x) dx ≤ 1 √ L z Lx 0 Lz 0 A ⊥ |∂ x φ ′ | 2 dzdx 1/2 Lx 0 A ⊥ |∂ x φ| 2 dx 1/2 . Thus Lx 0 A ⊥ |∂ x φ(x)| 2 dx 1/2 ≤ 1 √ L z Lx 0 Lz 0 A ⊥ |∂ x φ ′ | 2 dzdx 1/2 .
Equally, taking in (27) θ ′ as test function gives rise to

Lx 0 Lz 0 A z |∂ z θ ′ | 2 dzdx + ε Lx 0 Lz 0 A ⊥ |∂ x θ ′ | 2 dz ≤ -ε Lx 0 Lz 0 A ⊥ ∂ x φ ∂ x θ ′ dzdx ≤ ε Lx 0 Lz 0 A ⊥ |∂ x φ| 2 dzdx 1/2 Lx 0 Lz 0 A ⊥ |∂ x θ ′ | 2 dzdx 1/2 ≤ ε L z Lx 0 A ⊥ |∂ x φ| 2 dx 1/2 Lx 0 Lz 0 A ⊥ |∂ x θ ′ | 2 dzdx 1/2 . ( 28
)
This last inequality yields

Lx 0 Lz 0 A z |∂ z θ ′ | 2 dzdx + ε Lx 0 Lz 0 A ⊥ |∂ x θ ′ | 2 dzdx ≤ εL z Lx 0 A ⊥ |∂ x φ| 2 dx ≤ ε Lx 0 Lz 0 A ⊥ |∂ x φ ′ | 2 dzdx < Lx 0 Lz 0 A z |∂ z φ ′ | 2 dzdx + ε Lx 0 Lz 0 A ⊥ |∂ x φ ′ | 2 dzdx .
In this last step we would have the "equality" if and only if

Lx 0 Lz 0 A z |∂ z φ ′ | 2 dzdx =
0. This is however only possible for functions depending exclusively on the x-coordinate, φ ′ (x), which is in contradiction with the fact that φ ′ = 0 and φ ′ = 0. Thus we have shown that ||T (φ ′ )|| * < ||φ ′ || * for φ ′ = 0, φ ′ ∈ U, which means that T is a contractive application on (U, || • || * ).

Numerical discretization and simulation results

This part of the paper is concerned with the numerical discretization of the AP-scheme ( 12)-( 13) and the comparison of the simulation results with those obtained via the Singular Perturbation model ( 5) and the Limit model (4).

Discretization

The numerical resolution of the Asymptotic Preserving system (12)-( 13) is done by means of the standard finite element method. Let us recall the variational formulation of the AP-formulation

         a 2 φ, ψ = ( f , ψ) - 1 L z c φ ′ , ψ , ∀ ψ ∈ W , a (φ ′ , ψ ′ ) + b( P , ψ ′ ) = ε(f ′ , ψ ′ ) -εc ψ ′ , φ , ∀ψ ′ ∈ V , b( Q, φ ′ ) = 0 , ∀ Q ∈ W , (29)
with the notation of section 2. Here φ ′ (x, z) ∈ V, φ(x) ∈ W as well as P (x) ∈ W are the unknowns and ψ ′ ∈ V, ψ ∈ W and Q ∈ W the test functions. The introduction of the Lagrange multiplier P (x) was explained in a simplistic manner in the preceding sections and will be analyzed in more details in section 4. Due to the equivalence of ( 29) and ( 21)-( 22), one can comment that the introduction of P (x) is superfluous, but this is not the case for the discretized equations. The property φ ′ = 0 is indeed automatically fulfilled since the right-hand side of equation ( 22) has a zero mean value along the z-coordinate. However the discrete implementation of this quantity introduces round-off errors which probably will destroy the zero mean value property and justify the introduction of the Lagrange multiplier.

For simplicity reasons we omitted here the ε-index of the solution (φ ′ ε , φ ε ), the parameter ε > 0 being considered as fixed.

To discretize now the system (29) we introduce the grid

0 = x 0 ≤ • • • ≤ x n ≤ • • • ≤ x Nx+1 = L x , 0 = z 1 ≤ • • • ≤ z k ≤ • • • ≤ z Nz = L z
χ n (x) :=          x -x n-1 x n -x n-1 , x ∈ I n-1 , x n+1 -x x n+1 -x n , x ∈ I n , 0 , else , κ k (x) :=          z -z k-1 z k -z k-1 , z ∈ J k-1 , z k+1 -z z k+1 -z k , z ∈ J k , 0 , else .
Thus we are searching for approximations φ ′ h ∈ V h , φh ∈ W h and Ph ∈ W h , which can be written under the form

φ ′ h (x, z) = Nx n=1 Nz k=1 α nk χ n (x)κ k (z) , φh (x) = Nx n=1 β n χ n (x) , Ph (x) = Nx n=1 γ n χ n (x) .
Inserting these decompositions in the variational formulation [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] and taking as test functions the hat-functions χ n and κ k gives rise to the following linear system to be solved in order to get the unknown coefficients α nk , β n and γ n

A 2 β = w , (30)
A 0 + ε (A 1 -D) B B t 0 α γ = ε v 0 , (31)
where the matrices A 2 ∈ R Nx×Nx , A 0 , A 1 , D ∈ R NxNz×NxNz and B ∈ R NxNz×Nx correspond to the bilinear forms [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF] and the right-hand sides are defined by

w n := ( f , χ n ) - 1 L z c(φ ′ h , χ n ) , v nk := (f ′ , χ n κ k ) -c(χ n κ k , φh ) = (g, χ n κ k ) , for all n = 1, • • • , N x ; k = 1, • • • N z and g(x, z) := f ′ (x, z) + ∂ ∂x A ′ ⊥ (x, z) ∂φ ∂x (x) . ( 32
)
Solving iteratively the linear systems (30)-(31) permits finally to get the unknown function φ h (x, z) = φh (x) + φ ′ h (x, z). The convergence of the iterations was proved for the continuous case in theorem 2.5 and can be identically adapted for the discrete case.

Numerical results

In this section we shall compare the numerical results obtained by the discretization of the Singular Perturbation model, the Limit model and the just presented Asymptotic Preserving reformulation. With this aim, we consider a test case where the exact solution is known. Let thus

φ e (x, z) := sin 2π L x x + ε cos 2π L z z sin 2π L x x , (33)
be the exact solution of problem ( 5), where we choose A ⊥ (x, y) = c 1 + xz 2 and A z (x, z) = c 2 + xz, with two constants c 1 > 0, c 2 > 0. The numerical experiments are performed with L x = L z = 10 and c 1 = c 2 = L z . The exact right-hand side f is computed by inserting (33) in [START_REF] Besse | A model Hierarchy for Ionospheric Plasma modeling[END_REF]. We denote by φ P , φ L and φ A , respectively, the numerical solutions of the Singular Perturbation model ( 5), the Limit model ( 4) and the Asymptotic Preserving formulation ( 12)-( 13). The comparison will be done in the l 2 -norm, that means

||φ e -φ num || 2 = 1 √ N i∈G |φ e (X i ) -φ num,i | 2 1/2 , ( 34
)
where φ num stands for one of the numerical solutions and φ e (X i ) is the exact solution evaluated in the grid point X i . The index i covers all possible grid indices, reassembled in the set G, and N is the total number of grid points. The linear systems obtained after the discretization of either the P-model, the L-model or the AP-formulation are solved thanks to the same numerical algorithm (MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF]). The purpose here is not to design a specific preconditioner for the resolution of these linear systems, but to point out the efficiency of the presently introduced AP-method to deal with a large range of anisotropy ratios. As can be seen from Table 1 and Figure 1, the finite element resolution of the Singular Perturbation model is precise only for large 0 < ε < 1, whereas the Limit model is accurate for small ε ≪ 1. The range of ε-values in which both the Singular Perturbation and the Limit models provide an accurate approximation of the solution shrinks as the mesh size is refined. For a coarse grid (with 50 × 50 points see figure 1(a)) this domain ranges from 10 -12 to 10 -3 while it is reduced to 10 -9 -10 -5 for the refined 500 × 500 grid (figure 1(b)). This question is determinant for the development of a model coupling strategy. Indeed it requires an intermediate area where both discretized models furnish an accurate approximation and we observe that for refined meshes this area may not exist. This reduction of the validity domain can be explained for both the L-model and P-model but for quite different reasons.

ε 10 1 10 -1 10 -4 10 -14 10 -16 AP-scheme 3.4 • 10 -2 7.8 • 10 -3 3.8 • 10 -3 2.7 • 10 -3 2.7 • 10 -3 2.7 • 10 -3 S-scheme 2.8 • 10 -2 4.5 • 10 -3 2.8 • 10 -3 2.7 • 10 -3 6.6 • 10 -2 1.2 L-model 9.9 1.0 • 10 1 1.0 • 10 -1 2.8 • 10 -3 2.7 • 10 -3 2.7 • 10 -3
Table 1: Absolute error in the l ∞ -norm for the approximation computed thanks to the AP-scheme, discretized Singular Perturbation and Limit models (S-scheme and L-model) as compared to the exact solution.

The numerical approximation computed via the Limit model is altered by both the discretization error of the numerical scheme and the approximation error introduced by the reduction of the initial Singular Perturbation problem to the Limit problem. For coarse grids, the global error is rapidly dominated by the scheme discretization error, but as the mesh is refined, the approximation error becomes preponderant, as the Limit model is precise only for small ε-values. The schemes implemented here are of second order, thus when the mesh size is divided by ten, the discretization error is reduced by one hundred. The global error for the L-model displayed in figure 1(a) does not depend on ε as soon as ε < 10 -3 . Below this limit the L-model is able to furnish a better approximation of the solution with vanishing ε, however the numerical scheme is not precise enough and consequently the global error does not decrease. For the refined mesh, this discretization error is lowered by two order of magnitudes and the global error is a function of ε as long as its value is greater than 10 -5 (Fig. 1(b)).

The analysis for the Singular Perturbation model is quite complementary. The accuracy of the approximation provided by the P-model is good for large ε-values and deteriorates rapidly for small ones. This can be explained by the conditioning of the linear system obtained by the P-model discretization. An estimate of the condition number for the matrix is displayed in figure 2 for two different grid sizes. This conditioning deteriorates with vanishing ε-parameter, which is coherent with the fact that, working with a finite-precision arithmetic, the Singular Perturbation model degenerates into an ill-posed problem. This also explains the blow up of the error displayed in figure 1 as soon as the conditioning of the matrix approaches the critical value of the double precision (materialized by the level 10 15 in Fig. 2). This limit is reached on more refined meshes for larger ε-values (ε ≈ 10 -12 on a 50 × 50 grid and ε ≈ 10 -10 on a 200 × 200 grid). As expected, the P-model, though valid for all ε-values, cannot be exploited numerically for small ε. The ε-region where both the P-model and the L-model are accurate all-together, shrinks dramatically with the size of the mesh, fact which motivates the development of the AP-method. The condition number estimate of the linear system providing the approximation of the solution for the AP-scheme is also plotted in Figure 2. The conditioning of the system is rather ε independent and this is due to the introduction of the Lagrange multiplier, which forces the system in the limit to remain well-posed. The accuracy of the AP-scheme is totally comparable to the P-model for the large values of ε and to the L-model for the smallest ones. The AP-formulation is a good tool for computing an approximation for the solution which is accurate uniformly in 0 < ε < 1 and is therefore of great practical interest. Note that this approximation is obtained thanks to an iterative sequence {φ ′ k } k∈N , constructed with the fixed-point mapping T defined in theorem 2.5. The convergence of this iterative process is analysed in figure 3 on a 200 × 200 grid for a large value of ε. The l 2 -absolute error between the mean respectively the fluctuating parts of the exact solution and the approximation provided by the AP-scheme are plotted as a function of the iteration number. The sequence is initiated with the zero function. With the iterative process, both components converge towards the solution until the precision of the schemes is reached. At this point, after roughly 27 iterations, the approximation can not be improved and a plateau is observed. The simulation domain is discretized by a 500 × 500 mesh. For the smallest value of a the support of the function is reduced to 5 cells in each direction. The source term δ h a is normalized, such that the maximal value of δ h a grows with vanishing a-parameter. In table 2 the maxima and minima of the numerical approximations computed by the APscheme (φ A ) and the discretized Singular Perturbation model (φ P ) are gathered for the two source functions δ h a . Only large ε-values are considered to verify the positivity of the numerical approximations. Indeed for very small ε the solution is reduced to its mean part which is the solution of a classical elliptic problem preserving the maximum principle. This means that the relevant question is related to configurations where the fluctuating part φ ′ has a significant contribution to the elliptic problem solution. In this range of large and intermediate ε values, both approximations are comparable. Only slight differences can be observed on the maxima for the smallest ε-parameters. The results of table 2

Numerical analysis of the AP-scheme

In this last part of the paper we shall concentrate on the numerical analysis of the Q 1 finite element scheme introduced in section 3.1 for solving

       - ∂ ∂z A z ∂φ ∂z -ε ∂ ∂x A ⊥ ∂φ ∂x + ε ∂ ∂x A ′ ⊥ ∂φ ∂x = εg , in Ω , ∂φ ∂z = 0 on Ω x × ∂Ω z , φ = 0 on ∂Ω x × Ω z , (35)
where g ∈ L 2 (Ω) is a given function, with mean value along the z-coordinate equal to zero, g = 0. Moreover we shall explain why we have to introduce the Lagrange multiplier in order to solve numerically this equation. We remark that in contrast to section 3 we omitted for simplicity reasons the primes for φ, which indicated the fluctuation functions with zero mean value. The weak form of (35

) is a(φ, ψ) = ε(g, ψ) , ∀ψ ∈ V , (36)
or equivalently m(φ, ψ) = ε(g, ψ) , ∀ψ ∈ V ,

where m(•, •) is the coercive bilinear form defined in [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF]. Let us now consider the corresponding discrete problem

a(φ h , ψ h ) = ε(g, ψ h ) , ∀ψ h ∈ V h , (38)
where the finite dimensional space V h ⊂ V was introduced in section 3.1. It can be seen that the property g = 0 induces also in the discrete case that φ h = 0. Thus, following the same arguments as for the continuous case, we can show that equation ( 38) is equivalent to

m(φ h , ψ h ) = ε(g, ψ h ) , ∀ψ h ∈ V h . (39)
The Lax-Milgram theorem implies then the existence and uniqueness of a discrete solution φ h ∈ V h . The next theorem gives an estimate of the discretization error ||φ -φ h || V . We shall suppose in the sequel, that the diffusion matrices A ⊥ , A z and the function f are regular enough, to be able to use standard regularity/interpolation results. Then we have the following discretization error estimate

||φ -φ h || 2 V = ||∂ z φ -∂ z φ h || 2 L 2 + ε||∂ x φ -∂ x φ h || 2 L 2 ≤ Ch 2 , ( 40
)
with a constant C > 0 independent of ε > 0. Moreover, as φ, φ h ∈ U, we have

||φ -φ h || 2 U ≤ Ch 2 .
Proof: The fact that both solutions φ and φ h belong to the space U is an immediate consequence of the fact that the right-hand side of the equation (36) (resp. (38)) satisfies g = 0. The discretization error estimate is rather standard. Denoting by φ I the interpolant of φ in the finite dimensional space V h , i.e.

φ I (x, z) := Nx n=1 Nz k=1 φ(x n , z k )χ n (x)κ k (z) ,
we have due to the coercivity of the bilinear form m(•, •)

c||φ -φ h || 2 V ≤ m(φ -φ h , φ -φ h ) = m(φ -φ h , φ -φ I ) ≤ c||φ -φ h || V ||φ -φ I || V . Thus ||φ -φ h || V ≤ c||φ -φ I || V .
Standard Q 1 finite element interpolation results [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] yield for the interpolation error

||∂ x φ -∂ x φ I || 2 L 2 + ||∂ z φ -∂ z φ I || 2 L 2 ≤ ch 2 (||∂ xx φ|| 2 L 2 + ||∂ zz φ|| 2 L 2 ) ,
and regularity results for the solution φ of (36), imply

ε 2 ||∂ xx φ|| 2 L 2 + ||∂ zz φ|| 2 L 2 ≤ cε 2 .
This last estimate can be found by applying standard H 2 regularity results on the solution φ ε of the initial Singular Perturbation problem (5) (after a change of variable ξ := √ εx) and then exploiting the decomposition φ ε = φ ′ ε + φε . Thus, we have altogether with a constant c > 0 independent of ε > 0

ε||∂ x φ -∂ x φ h || 2 L 2 + ||∂ z φ -∂ z φ h || 2 L 2 ≤ ch 2 .
What is important to observe from the error estimate (40) is that for ε → 0 the error ||φ -φ h || H 1 in the standard ε-independent H 1 -norm blows up. This is one argument why the Singular Perturbation model is inaccurate for ε ≪ 1. However, in the case where φ and φ h are elements of the space U, we have ||φ -φ h || U ≤ Ch 2 independently of ε, which means that we have convergence of the scheme in (U, || • || U ), uniformly in ε > 0. The Poincaré inequality implies then the uniform convergence in the || • || L 2 norm. The AP-scheme is thus equally accurate for every value of 0 < ε < 1.

The discretization error φ -φ h is not the only error we are introducing when solving numerically (38) instead of (36). Indeed, (38) is nothing but a linear system

M α = v , (41)
to be solved to get the unknowns α nk := φ h (x n , z k ), where v nk := ε(g, χ n κ k ) and the discrete solution of (38) is then reconstructed as

φ h (x, z) = Nx n=1 Nz k=1 α nk χ n (x)κ k (z) .
Unfortunately the implementation of the system (41) introduces round-off as well as approximation errors due for example to the numerical computation of a(χ n κ k , χ r κ p ). Thus the numerical resolution of (41) does not yield the exact solution, but an approximation (α nk ) nk , solution of the slightly perturbed system M α = ṽ .

We are now interested in the error estimate ||φ h -φh || V , as a function of the perturbation ||v -ṽ|| 2 , where || • || 2 denotes the Euclidean norm in R NxNz .

Theorem 4.2 Let α be the exact solution of (41) and α the exact solution of the perturbed system (42). Let φ h ∈ V h and φh ∈ V h denote the corresponding functions Then we have

ε||∂ x φ h -∂ x φh || 2 L 2 + ||∂ z φ h -∂ z φh || 2 L 2 ≤ c ε ||v -ṽ|| 2 2 , (43)
with a constant c > 0 independent of ε > 0 and h > 0. However, if both functions φ h and φh belong to U, then we have the ε-independent estimate Thus we get with a constant c > 0 independent of ε

||e h || V ≤ c √ ε ||M E|| 2 = c √ ε ||v -ṽ|| 2 .
In the case the two functions φ h and φh belong to

) = ε(g, ψ h ) , ∀ψ h ∈ V h , b(Q h , φ h ) = 0 , ∀Q h ∈ W h , (44)
where W h ⊂ W was constructed in section 3.1. As mentioned in the continuous case this problem is equivalent for ε > 0 to the discrete problem (38). If φ h ∈ V h is the unique solution of (38), then (φ h , 0) ∈ V h × W h is a solution of (44). And if (φ h , P h ) ∈ V h × W h solves (44), then P h ≡ 0 and φ h ∈ V h is the unique solution of (38). This last statement is immediately proved by taking in the variational formulation (44) only x-dependent test functions ψ h (x) ∈ V h . By doing this, we can be sure that the numerical solution φh of (44) satisfies φh = 0, such that the error ||φ h -φh || U is uniformly bounded. This proves that the introduction of the constraint φ h = 0 in the AP-formulation is crucial and avoids the numerical difficulties associated with the original P-model.

Conclusion

In this paper we have introduced an Asymptotic Preserving formulation for the resolution of a highly anisotropic elliptic equation. We have shown the advantages of the AP-formulation as compared to the initial Singular Perturbation model and to its limit model, when the asymptotic parameter goes to zero. It came out that the AP-scheme is a powerful tool for the resolution of elliptic problems presenting huge anisotropies along one coordinate, and gives access to the simulation in a very easy and precise manner. The Asymptotic-Preserving method developed here relies on the decomposition of the solution in its mean part along the anisotropy direction, and a fluctuation part. This integration along the anisotropy direction is easily performed in the context of Cartesian coordinate systems with one coordinate aligned with the direction of the anisotropy. In a forthcoming work [START_REF] Brull | An asymptotic preserving scheme in the drift limit for the Euler-Lorentz system for a variable magnetic field[END_REF] this procedure is extended to more general anisotropies.

  and consequently also in (U, || • || V ), towards the fixed-point φ ′ * ∈ U of T . The proof of this theorem is based on the following Lemma 2.6 [8] Let (U, || • || * ) be a normed space and T : U → U a contractive application, i.e. ||T (φ) -T (ψ)|| * < ||φ -ψ|| * , ∀φ, ψ ∈ U with φ = ψ .

  and denote the cells byI n := [x n , x n+1 ] and J k := [z k , z k+1 ].The finite dimensional spaces V h ⊂ V and W h ⊂ W are constructed as usual, by means of the hat functions (Q 1 finite elements)

  Grid with 50 × 50 points.

  Grid with 500 × 500 points.

Figure 1 :

 1 Figure 1: Absolute error in the l 2 -norm between the computed solutions φ P , φ L , φ A and the exact solution φ e , as a function of ε and on different grids. Dashed lines : (S) Standard scheme : discretization of the P-model; Stars : (AP) AP-scheme; Circles : (L) discretization of the L-model.

Figure 2 :Figure 3 :

 23 Figure 2: Condition number estimate for the discretization matrices of the Standard (S) and AP schemes (computed by LAPACK [4]) as a function of ε. Different grids of 50×50 and 200×200 points and different ε-values are used. Dashed/Plain lines : 200 × 200 / 50 × 50 grid ; Stars : AP-scheme.

Theorem 4 . 1

 41 Let φ ∈ V be the unique solution of the continuous problem (36) and φ h ∈ V h the unique solution of the discrete problem (38). Both solutions are elements of the normed space (U, || • || U ), whereU := {ψ(•, •) ∈ V / ψ = 0} with ||ψ|| U := ||∂ z ψ|| L 2 (Ω) .

α

  nk χ n (x)κ k (z) , φh (x, z) = Nx n=1 Nz k=1 αnk χ n (x)κ k (z) .

  ||φ h -φh || U ≤ c||v -ṽ|| 2 .Proof: Let us denote within this proofE nk := α nkαnk for n = 1, • • • , N x , k = 1, • • • , N z and e h (x, z) := φ h (x, z) -φh (x, z), such that e h (x, z) = Nx n=1 Nz k=1 E nk χ n (x)κ k (z) .Moreover let N := N x N z and Y ∈ R N be an arbitrary vector associated with the functiony h (x, z) = Nx n=1 Nz k=1 Y nk χ n (x)κ k (z).Then we have with (•, •) 2 the euclidean scalar product in R N and M the discretization matrix of (41)||M E|| 2 = sup Y ∈R N ,Y =0 (Y, M E) 2 ||Y || 2 = sup

  we deduce that this sequence has a unique accumulation point φ ′

* in U. This means that the sequence {φ ′ k } k∈N converges in (U, || • || * ) towards φ ′ * . Due to the fact that || • || * and || • || V are equivalent norms, we have also the convergence in (U, || • || V ).

It remains to show that T is contractive. For this let φ ′ 1 , φ ′ 2 ∈ U be two given, distinct functions. Denoting by φ ′

  demonstrate the positivity of the approximations computed by either the AP-scheme or the Singular Perturbation model.

		ε	10 2	10	1	10 -1	10 -2	10 -3
	a = 10 -1	max(φ P ) max(φ A ) min(φ P ) min(φ A )	77.58 77.58 1.9 10 -7 2.5 10 -7 2.4 10 -2 2.4 10 -2 2.8 10 -2 2.8 10 -2 3.82 1.63 8.93 7.22 6.93 3.82 1.63 8.93 6.89 6.89 1.9 10 -7 2.5 10 -7 2.4 10 -2 2.4 10 -2 2.8 10 -2 2.8 10 -2
	a = 10 -2	max(φ P ) max(φ A ) min(φ P ) min(φ A )	1.8 10 2 1.8 10 2 1.6 10 -7 2.5 10 -3 2.4 10 -2 2.8 10 -2 2.8 10 -2 2.8 10 -2 7.1 10 1 2.6 10 1 1.2 10 1 8.29 7.34 7.1 10 1 2.6 10 1 1.2 10 1 7.14 7.11 1.6 10 -7 2.5 10 -3 2.4 10 -2 2.8 10 -2 2.8 10 -2 2.8 10 -2

Table 2 :

 2 Maxima and minima of the numerical solutions computed thanks to the AP-scheme (φ A ) and the Singular Perturbation model (φ P ). The elliptic problem is solved with the Dirac δ h a function as a source term on a 500 × 500 mesh.

  Y ∈R N ,Y =0 m(y h , e h ) ||Y || 2 .

	Due to the fact that	||Y || 2 ≤ c||y h || L 2 ≤	c √ ε	||y h || V ,
	we have
	||M E|| 2 =	sup Y ∈R N ,Y =0	m(y h , e h ) ||Y || 2	≥ c	√ ε	sup y h ∈V h ,y h =0	m(y h , e h ) ||y h || V	≥ c	√	ε||e h || V .

  U, i.e. e h ∈ U, we can exploit the fact that in U the Poincaré inequality gives rise to||Y || 2 ≤ c||y h || L 2 ≤ c||y h || U . This yields, as m(•, •) is also coercive on U, that ||M E|| 2 = sup Y ∈R N ,Y =0 m(y h , e h ) ||Y || 2 ≥ c sup y h ∈U ,y h =0 m(y h , e h ) ||y h || U ≥ c||e h || U .and thus the ε-independent estimate is proved.Similarly as for the discretization error, we can deduce from the round-off error estimate (43) that, for ε → 0, the standard H 1 -norm ||φ h -φh || H 1 explodes. However if we impose that both solutions φ h and φh are elements of the space U, space of functions with mean value along the z-coordinate equal to zero, then we have the uniform estimate ||φ h -φh || U ≤ c||v -ṽ|| 2 , and by the Poincaré inequality ||φ h -φh || L 2 ≤ c||v -ṽ|| 2 . Unfortunately even if we know that φ h ∈ U, this is not necessarily true for φh , if we discretize (35). But it can be achieved by forcing the numerical solution φh to satisfy φh = 0. Indeed, this can be done by introducing explicitly in the discrete problem (38) the constraint φ h = 0, such that it is much more ingenious to solve instead a(φ h , ψ h ) + b(P h , ψ h

Acknowledgments

This work has been partially supported by the Marie Curie Actions of the European Commission in the frame of the DEASE project (MEST-CT-2005-021122) and by the CEA-Cesta in the framework of the contracts 'Dynamo-3D' # 4600108543 and 'Magnefig' # 06.31.044.


QUARTZ Roberta C Affonso
email: roberta.costa@supmeca.fr

Florent Couffin
email: florent.couffin@supmeca.fr

Aggregated Indicator for Assessing Station Criticality of Bike Sharing Systems

Keywords: bike sharing system, aggregated performance indicator, users' satisfaction, rebalancing operations

Bike Sharing Systems (BSS) are implemented in numerous cities around the world. The users' satisfaction is directly affected by the availability of bicycles and docking points at the stations. The assessment of users' satisfaction is an important issue to enable the system manager to understand, manage and improve the service quality of the BSS. This paper proposes to define an aggregated performance indicator for BSS stations. This indicator allows to measure the station criticality by integrating the unavailability of bikes and dockings of this station, and also of the neighboring stations. The integration of neighboring stations will be defined according to the distance that users is willing to walk to reach them. This indicator is applied to a part of Vélib' system and will be analyzed in relation to a non-aggregated criticality indicator. The proposed aggregated criticality indicator can be implemented, for example, to define priority stations to be rebalanced by rebalancing operations. Open data, over a six weeks period, from JCDecaux (Vélib' operator) was used to implement this study.

I. INTRODUCTION

The Bike Sharing Systems (BSS) are widely implemented in most of medium and large cities around the world. These systems must meet the needs of users traveling between two stations: a user must find an available bicycle at his departure station and an available docking at his arrival station. Some stations have, for specific time slots, a tendency to completely empty out and others to completely fill up. Bicycle and docking unavailability decreases the service quality of these systems. This is the reason why some operators deploy fleets of rebalancing vehicles which take bicycles from full stations to bring them to empty stations [START_REF] De Chardon | Bike-share rerebalancing strategies, patterns, and purpose[END_REF]. These operations have an impact on the operating costs of these systems and cannot be extensively deployed [START_REF]The bikeshare planning guide[END_REF]. These operations correspond around 3% of users' flow. Thus, from the point of view of the system manager, the network rebalancing must co nsider an assessment of the stations to define which ones must be replenished and emptied as a priority.

In the literature, studies do not take into account the bicycle and docking availability of neighboring stations, while the user can possibly move a short distance to look for or deposit a bike, if the desired station cannot meet the demand. The neighboring stations availability must therefore be integrated in the assessment of the stations.

In this context, we propose an aggregated performance indicator that measures the criticality of a station considering its availability in terms of bicycles and dockings, and also the availability of neighboring stations. Indeed, a station can be considered as critical if for some time, it is out of bikes or dockings. But this criticality is even more important if the stations around it are also out of bikes or dockings respectively. The proposed aggregated indicator and a non-aggregated one will be applied to a set of stations. The objective is to analyze the criticality of a station considering the volume of available bicycles and dockings at a given area around this station, and the criticality of a station considering its own status. The proposed aggregated indicator should identify the stations which must be replenished or emptied as priority.

The paper is organized as follows: in Section II we present a state of the art about the control operation, in particular the rebalancing of the BSS, and the performance indicators considered in these studies. Section III introduces the study methodology. The proposed method to define aggregated criticality indicator of a station integrating neighboring stations performance indicator is defined in section IV. The results and analyzes are presented in section V. Then conclusions and future works are discussed in section VI.

II. SERVICE QUALITY AND BSS CONTROL

The main challenge of BSS is to ensure a high availability of bicycles and dockings in order to satisfy the users. Short rental times, a very strong spatial and temporal irregularity of user demand, lead to imbalances in the distribution of bicycles [START_REF] Vogel | Strategic and operational planning of bike-sharing systems by data mining--a case study[END_REF]. BSS are considered complex stochastic dynamic systems requiring control operations [START_REF] Gavalas | Design and management of vehicle-sharing systems: a survey of algorithmic approaches[END_REF]. The control operations most applied to these systems and most studied in the literature is the rebalancing one. Vehicles, dedicated to rebalancing operations, take bicycles from the full stations to drop them in empty stations. Some studies consider one or more stock location in which bicycles can be taken or stored. When control operations are performed at night, when the system is not in use, it is named static control; when control occurs while the system is in use, it is named dynamic control. Control decisions aim to meet users' demand by providing them with a bike at their departure station and a docking at their arrival station.

Many scientific works aim to solve the problem of rebalancing by vehicle route formalized by a linear program. The objective function consists in minimizing the total unsatisfied demand while seeking to minimize the rebalancing costs (number of vehicles assigned, number or distances of the routes...) [START_REF] Ho | Solving a static repositioning problem in bike-sharing systems using iterated tabu search[END_REF]. For example, in the context of static regulation, [START_REF] Forma | A 3-step math heuristic for the static repositioning problem in bike-sharing systems[END_REF] seeks to minimize a penalty function that measures out of bicycles and dockings. He uses a stock location in which it can take bicycles. [START_REF] Kadri | A branch-and-bound algorithm for solving the static rerebalancing problem in bicycle-sharing systems[END_REF] seeks to minimize the time during which the stations are unbalanced compared to a target stock level in order to meet bicycle and docking demand. [START_REF] Alvarez-Valdes | Optimizing the level of service quality of a bike-sharing system[END_REF], which proposes a rebalancing algorithm, considers "the measure of service quality of the bike-sharing system considered in this study is a measure of cost that combines the two possible situations of unsatisfied demand, namely: (a) A user wishing to withdraw a bicycle finds the station empty. (b) A user returning a bicycle finds the station full". [START_REF] Contardo | Rebalancing a dynamic public bike-sharing system[END_REF] deals with the problem of dynamic rebalancing, considering a succession of periods; he seeks to minimize the unsatisfied demands at the end of each period. [START_REF] Gaspero | Rebalancing bike sharing systems with constraint programming[END_REF] uses constraint programming to solve the routing problem of static rebalancing vehicles by minimizing the gap between user demand and the number of bikes to be brought into each station.

To solve the rebalancing problem, non-deterministic methods are also used. [START_REF] Raviv | Optimal inventory management of a bikesharing station[END_REF] implements a Monte Carlo simulation method to determine the optimal number of bicycles in the stations. Author believes that the shortage of bicycles and dockings are the two main raisons of user dissatisfaction. He defines a penalty function integrating the number of users lost when there is no available bicycle or docking in a station. [START_REF] Wang | Analyzing bike repositioning strategies based on simulations for public bike sharing systems: simulating bike repositioning strategies for bike sharing systems[END_REF] uses simulation to compare several rebalancing policies. To take into account the dissatisfaction of users, he takes into account the out of bicycles or dockings, but also the time customers wait for an available bicycle or docking. [START_REF] Li | A unified framework for analyzing closed queueing networks in bike sharing systems[END_REF] analyzes the BSS performance by applying/implementing a closed queuing network model. He evaluates the steady-state network by computing the probabilities of stations to be out of bicycles or dockings.

All of these works aim to evaluate or improve the availability of bicycles and dockings for effective control. However, the operator's decision to add or remove bicycles from certain stations depends solely on the status of these stations regardless of the status of the surrounding stations. But a user may agree to walk to another station when he cannot find a bicycle at his departure station. Indeed, [START_REF] Ai | A solution to measure travelers transfer tolerance for walking mode and dockless bike-sharing mode[END_REF] studies the behavior of the users, for example, the time that a user is ready to wait for a bicycle or a docking to be released at a station, but also the distance that he agrees to walk to find an available bicycle. The author specifies that "an acceptable walking distance is considered as between 200 and 400 m or between 5 and 10 min". The [START_REF] Watson | Walking for transportation: what do US adults think is a reasonable distance and time?[END_REF] study, which analyzes more generally the distance a user is willing to walk to reach a transport, concludes that "7% adults thought it was not reasonable to walk at all, 57% thought walking up to 1 mile was reasonable, and 43% thought walking a mile or more was reasonable" and it estimates that the average acceptance distance is around 400 m. A study presented in [START_REF] Millward | Active-transport walking behavior: destinations, durations, distances[END_REF] shows that a user is willing to walk an average of 0.51 km to get to a bus stop or ferry terminal. [START_REF] Zhang | Exploring the impact of built environment factors on the use of public bikes at bike stations: case study in Zhongshan, China[END_REF], which studies spatial correlations between nearby BSS stations, verifies that station demand is positively correlated with demand from neighboring stations (within 300 m). [START_REF] Kabra | Bike-share systems: Accessibility and availability[END_REF] estimates the probability of using a station is directly related to the distance a user must walk to reach it (80% of the users of a station come from a point located less than 300 m around the station) and to the availability of bicycles (10% more availability leads to 12% more use).

In this sense, the bicycle or docking availability of a station should not be the only factor to take into account for its rebalancing. Actually, users are ready to walk to a nearby station, especially since operators provide applications indicating the station status [START_REF] Midgley | Bicycle-sharing schemes: enhancing sustainable mobility in urban areas[END_REF]. Thus, to rebalancing a station, it is necessary to take into account its status, but also that of neighboring stations; their distance from the studied station must be taken into account according to the distance that the users agree to walk.

III. STUDY METHODOLOGY

This study is based on the data provided by the operator of the Vélib' network, JCDecaux company. These data concern the flow of bicycles entering and leaving each station of the network, and the status of the stations (number of bicycles and docking available in each station). In order to obtain data corresponding a homogeneous system behavior, we have only used the data of worked days (Saturday and Sunday excluded) for the period September 04, 2017 to October 13, 2017 (30 days). This period was chosen because it does not include school holidays or holidays. In fact, the behavior of the users is very different between the working days and the week-end days and holidays. Also, we were interested in the time slot from 8 am to 10 am, which is the morning peak hours of the system usage [START_REF] Chabchoub | Classification of the vélib stations using Kmeans, Dynamic Time Wraping and DBA averaging method[END_REF].

Our study method is structured in 3 main steps:

-The average time of disruption (shortage)

The average time of bicycle and docking disruption was calculated. This time measures the bicycle and docking unavailability of a station during the studied time slot (from 8 am to 10 am). The 20 stations representing the most important disruption on bicycles and the 20 stations representing the most important disruption on docking were selected.

-The aggregated average time of disruption (shortage) Considering the disruption average time of the previously selected stations, the aggregated associated indicator was defined, integrating the neighboring stations and the distance users accept to walk to reach a bicycle or a docking. The method to select neighboring stations and to calculate the aggregated indicator is presented in the next section.

-The analysis of station criticality

The deviation between both critical indicators' values (non-aggregated and aggregated one) was identified and analyzed. Also, the station criticality ranking in function of both indicators was examined in a rebalancing operation perspective.

IV. PROPOSED AGGREGATED CRITICALITY INDICATOR

As previously mentioned (introduction section), the criticality of a station must take into account the availability (bicycles or dockings) of the neighboring stations and the distance users are willing to walk. For a given availability, a nearby station will have more dominant influence than a station further away. We propose an indicator that integrates this dimension.

To statistically assess the break-even point of distance that users are willing to walk, we carried out a survey of 136 Vélib system users. The selected users are occasional users, who take a bicycle at least once a month on a working day, and regular users, who take a bicycle every working days. The question asked was: "What is the maximum distance you can walk to get a bicycle from a Vélib station?". The results are summarized in Table I. For example, 86 in 136 users (63%) agree to walk till 300 m to get a bicycle. Therefore, for a station, we decompose its neighborhood in concentric circles. A zone may have stations inside the area between two consecutive circles. We define the acceptance rate (arj) of a zonej as the average of the users' proportion (arrj) who agree to walk at a distance between the radius of the inner circle (aupj-1) and the radius of the outer circle (aupj), or [START_REF] De Chardon | Bike-share rerebalancing strategies, patterns, and purpose[END_REF] To avoid associating a station with a large area, the entire neighboring area is limited to 500 meters around a studied station. Also, only 8% of users agree to walk beyond this distance. The Fig. 1 represents the neighboring zones around a studied station. For each zone, the figure presents its radius, the proportion of users agreeing to walk to the outer limit of this zone and its acceptance rate. For example, the zone3 is located between two circles, which respective radii are 200 and 300 m; 63% of users agree to pick up a bicycle to the outer limit of this zone, and the acceptance rate for this zone is equal to 74%. In order to verify that the studied BSS configuration allows us to take advantage of the users' ability to fetch a bike from the nearby stations of an empty station, we have calculated and exploited the distance between the stations through their GPS location. On average over the entire network (1224 stations), we can find 4.62 stations located within 500 m around a given station. This value varies according to the geographical situation of the studied station (155 stations in the Paris center have more than 15 stations within 500 meters (13%), whereas 24 peripheral stations have none (1.96%)).

Table II summarizes the average number of neighboring stations per 100 meters around a given station. The first column characterizes an area between two circles (minimum and maximum radius) concentric to the central station (the same division as for the acceptance areas) and the second one gives the average number of stations in this zone. As a result, users can find in most of cases (98%) another station within the acceptance zone (less than 500 meters). The performance indicator is measured for all stations (studied station and neighboring stations). This measure is associated with the zone acceptance rate (ar) of each station. Then, the aggregated performance indicator of a station i is determined by computing the weighted sum of the performance indicators of all stations (station i and neighboring stations). The aggregated performance indicator must be standardized related to the weighted sum. Indeed, to compare the stations' criticality we must standardize this indicator value because the number of neighboring stations of each zone are not the same.

We want to calculate the aggregated average disruption time of the station i (adtsi). Around the station i five zones can exist, and are associated to the index j. Inside each zone j, there are nj stations and k is the index associated to the k-th station. The weighted average disruption time of the k-th station of the j-th zone is defined in (2).

(2)

Where dtsi,j,k is the average disruption time of the k-th station of the j-th zone.

The aggregated average disruption time of the station i (adtsi) is defined as the weighted sum of dtsi,j,k standardized related to the weighted sum. This indicator is defined in (3).

(3) Table III presents the computation of the aggregated average disruption time of the station 32604, which is one of the 20 stations with most important time disruption of bicycles. The standardized value of this indicator is 4136.5/3.56 which is equal to 1161.9 sec. The computation of this aggregated criticality indicator was performed for the 20 stations representing the most important disruption on bicycles and the 20 stations representing the most important disruption on dockings. The results and analysis are presented in the next section.

V. RESULTS AND ANALYSIS In order to analyze the usefulness of the proposed aggregated criticality indicator, we have compared it with the average time disruption of each station in a rebalancing operation perspective. Indeed, criticality indicators can guide mangers to define the priority stations to be replenished or emptied. -Aggregated indicator | non aggregated one: in this case, we can identify two different situations. In the first one, the disruption time of neighboring station is very similar to the studied station one (Fig. 2, stations 32605, 32603). The second situation concerns isolated stations (with no neighboring station). Some stations near from a big public park or in the suburb of Paris, the end-limit of the network, are in this situation. In Fig. 3, we can see the station 12123. We can notice in Fig. 2 and 3 that the relative deviation between both indicators in the case of dockings (33.6%) is greater than for the bicycles (12.8%). On the entire network, only 166 stations (the average non-aggregated disruption time is 260.8 sec.) have disruptions on dockings while 568 have disruptions on bicycles (the average non-aggregated disruption time is 653.4s). In addition, the curve representing the aggregated average disruption time of dockings decreases much faster than the curve of the same indicator for bicycles. Indeed, over the studied period (8:00 am to 10:00 am), many bicycles are borrowed; the zones where we can find a disruption concentration is more important in the case of bicycles (users not finding a bicycle at a station go to a neighboring station which will be in turn out of bicycle) than dockings' one, for which this effect is time-delayed.

Also, we have ranked the stations in function of their average disruption time and of their aggregated disruption time. Fig. 4 presents this rank for bicycle case. The average rank deviation is 3.9 (it is 8.03 for docking case). In view of this deviation, we have identified for the first 20 stations, with the most important average disruption time, their place in the aggregated indicator ranking, considering all the stations of the system. These 20 stations are ranked between 67 th and 870 th place (their average rank is 496). Thus, we verify that the two indicators (aggregated and non-aggregated) results are very different. This means that the priority stations for rebalancing operations are not the same in function of the applied indicator. Rebalancing operations represent just 3% of the usage flow, and static rebalancing operations (00:00 till 07:59 am) concerns about 100 stations on average. Thus the criticality stations' ranking is very important to rebalance the system as effectiveness as possible. We think so that the proposed aggregated criticality indicator is a good reference for managers to decide the rebalancing operations.

An aspect that should be taken into account by the manager is that the overlapping of neighboring zones. Indeed, some stations integrate the neighboring zones of more than one critical station. Because of this when the manager decides to rebalance a station, the aggregated indicator of all others stations must be recomputed. For a deeply evaluation of the proposed aggregated criticality indicator, an application of this rebalancing policy must be done, for example through a simulator. Some research studies ( [START_REF] Soriguera | A simulation model for public bike-sharing systems[END_REF], [START_REF] Caggiani | A dynamic simulation based model for optimal fleet repositioning in bike-sharing systems[END_REF]) propose simulation models for BSS, and can be adapted to evaluate rebalancing operations policies, taking into account the availability of neighboring stations and the demand transfer from a station out of bicycles (or dockings) toward the neighboring station in function of the acceptance rate. We have applied this indicator to 20 critical real stations (Vélib' system) on bicycles and dockings, and we have compared the results with the average time disruption of each station. We have noted some deviations that show the influence of the neighboring stations' availability (value and rank) in the criticality of the studied station. In order to improve rebalancing operations, we think that the proposed indicator is more relevant to aid managers.

For future works, we will develop a decision support tool to define the priority stations to be rebalanced, considering the proposed indicator, and assess the rebalancing operation. As a first step, this tool will be developed to help managers defining the static rebalancing of the bike sharing system and better prepare it for high-demand schedules. In a second step, we will integrate dynamic rebalancing, which will consider rebalancing vehicle routing, the average travel time on bicycle and the time of walking by users to other stations.

Fig. 1 .

 1 Fig. 1. Neighboring areas of a station.

Fig. 2

 2 Fig. 2 and 3 compare both mentioned indicators for bicycle and docking criticality respectively. For more details, Tables IV and V present the value of these indicators just for the first 10 stations. Through the analysis of these data we have identified 3 different situations: -Aggregated indicator < non aggregated one: the neighboring stations mitigate the criticality level of the studied station. This mitigation can be more (or less) important in function of the number and the status of the stations around the studied station. This is the case of most stations presented in the Fig. 2 (stations 19105, 12126, 16130) and 3 (stations 12006, 12014, 4016).-Aggregated indicator > non aggregated one: the neighboring stations increase the criticality level of the studied station, because the neighboring ones present a longest disruption time. The stations in this situation do not appear in the figures presented here because they are in the middle rank of stations.

Fig. 3 .

 3 Fig. 3. Docking disruption time.

Fig. 4 .

 4 Fig. 4. Criticality ranks of stations

TABLE I .

 I SURVEY RESULTS CONCERNING THE DISTANCE USERS ARE

		WILLING TO WALK
	Distance (m)	Number of Users	User Proportions (aup)
	0	136	1.00
	100	127	0.93
	200	115	0.85
	300	86	0.63
	400	60	0.44
	500	35	0.26
	More than 500	11	0.08

TABLE II

 II

	.	AVERAGE NUMBER OF NEIGHBORING STATIONS
		Zone	Average number of stations
	0 < distance ≤ 100	0.085
	100 < distance ≤ 200	0.69
	200 < distance ≤ 300	1.63
	300 < distance ≤ 400	2.61
	400 < distance ≤ 500	2.98

TABLE III .

 III

		COMPUTATION OF WEIGHTED AVERAGE DISRUPTION TIME
			OF THE STATION
						Weighted
					Average	Average
	Stations'	Stations	Zones	Acceptance	disruption	disruption
	Number	(k)	(j)	Rate (arj)	time of the	time of
					station (sec.)	the station
						(sec.)
	32604		0	1	1164.7	1164.7
	32606	1	3	0.74	1177.5	871.4
	32605	2	3	0.74	1202.1	889.5
	32603	1	4	0.54	1178.6	636.4
	32602	2	4	0.54	1063.9	574.5
			Sum	3.56		4136.5
	Standardized aggregated average disruption time of the station 32604	1161.9