Éliane Bécache

AND Patrick Joly

Chrysoula Tsogka

E Becache

An Analysis of New Mixed Finite Elements for the Approximation of Wave Propagation Problems

Keywords: mixed finite elements, mass lumping, anisotropic waves

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction.

In this paper we develop and analyze a new family of mixed finite elements for the scalar anisotropic wave equation. It is essential, for the applications we have in mind, to use the mixed formulation involving both velocity and pressure and not the traditional displacement formulation. Let us first explain our motivations for introducing such mixed formulation.

This work falls within the more general framework of developing efficient numerical methods for approximating wave propagation in complex media such as anisotropic heterogeneous media with cracks or obstacles of arbitrary shapes. As a first step toward the more complicated case of elastic waves with a free boundary condition (on free surface or crack), we consider the scalar anisotropic wave equation, with a Neumann boundary condition on the boundary Γ of an obstacle (or crack).

The characteristics of wave propagation problems (large-scale computations, complex geometries for the cracks or the obstacles, unbounded domains) together with the need to construct an efficient method lead us to adopt the following criteria.

To facilitate the implementation and reduce the calculation time, we use regular meshes, squares in two dimensions, and cubes in three dimensions. The question with such a choice is how to take into account in an accurate way complex geometries. A staircase approximation of the geometry introduces artificial singularities on the boundary which produce spurious diffractions (see [START_REF] Collino | Fictitious Domain Method for Unsteady Problems: Application to Electromagnetic Scattering[END_REF] for a numerical illustration of this phenomenon). To reduce these spurious diffractions a very fine mesh in space is required. Moreover, because of the stability condition, the time discretization step is then very small for an explicit scheme.

To avoid this, we intend to use the fictitious domain method (cf. [START_REF] Collino | Fictitious Domain Method for Unsteady Problems: Application to Electromagnetic Scattering[END_REF][START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF][START_REF] Glowinski | A fictitious domain method for Dirichlet problem and applications[END_REF]). The basic principle of this method is to extend the solution to a domain of simple geometry that ignores the obstacle or the crack (typically a rectangle in two dimensions), the boundary condition being seen as an essential condition, i.e., an equality constraint in a given functional space. To apply this method to the wave equation with a Neumann boundary condition, we need to rewrite the Neumann condition as an essential one. This can be achieved by introducing the unknown p = ∇u, where u is the solution of the scalar wave equation. The Neumann condition on u can then be written as an essential condition on p. This justifies our choice to consider the wave equation as a first order system in (v = ∂ t u, p). The boundary condition is then taken into account via the introduction of a Lagrange multiplier, which can be interpreted as the jump of v through the boundary. We have in that case two computational domains, a domain of simple geometry in which a uniform mesh can be used for the unknowns (v, p) and the boundary in which a conforming mesh is considered for the Lagrange multiplier. This is our main reason to use a mixed variational formulation. However, in this paper, we will not consider the presence of the crack, and we refer for that to a future paper.

A second advantage of the first order system is that it permits the use of a new absorbing layer model, inspired by the perfectly matched layer (PML) introduced by Bérenger for the two-dimensional (2D) Maxwell problem (cf. [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF][START_REF] Collino | Conditions et couches absorbantes pour les équations de Maxwell[END_REF]). This model is very efficient compared with the traditional absorbing boundary conditions and can be naturally extended to any first order system [START_REF] Collino | Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heteregeneous media[END_REF].

Finally, in order to get an efficient method we wish to obtain an explicit time discretization scheme through mass lumping. This is why we consider in this paper nonstandard mixed finite elements (inspired by Nédélec's second family [START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF]) which we shall show to be compatible with mass lumping. These elements are in fact the scalar version of a new family of mixed finite elements we have introduced in [START_REF] Bécache | Eléments finis mixtes et condensation de masse en élastodynamique linéaire. (I) Construction[END_REF] for the velocity-stress formulation of elastodynamics.

We construct arbitrary higher order elements which are of great importance in the case of large-scale problems. In this sense, this work can be considered as the continuation of many articles about higher order methods for wave propagation problems (cf. [START_REF] Alford | Accuracy of finite difference modeling of the acoustic wave equation[END_REF][START_REF] Cohen | Construction analysis of fourth-order finite difference schemes for the acoustic wave equation in nonhomogeneous media[END_REF][START_REF] Cohen | Construction and analysis of higher order finite elements with mass lumping for the wave equation[END_REF][START_REF] Cohen | Gauss point mass lumping schemes for Maxwell's equations[END_REF][START_REF] Mulder | A comparison between higher-order finite elements and finite differences for solving the wave equation[END_REF][START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF]).

The main objective of the present paper is the analysis of these new elements. We shall show that they do not satisfy the classical Ladyzhenskaya-Babuska-Brezzi approximation condition (cf. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[END_REF][START_REF] Babuska | The finite element method with Lagrangian multipliers[END_REF]) (defect of coerciveness). Similar difficulties have been raised by Nédélec in [19, p. 80]. Our main contribution is to present a method for overcoming this difficulty. The abstract results concerning the elliptic problem have been announced in a short note [START_REF] Bécache | Etude d'un nouvel élément fini mixte permettant la condensation de masse[END_REF]. Let us mention that the generalization of this finite element analysis to the case of elastic waves presents a second major difficulty linked to the symmetry condition of the stress tensor. We intend to treat this difficulty in a forthcoming paper.

The present paper is organized as follows. In section 2 we recall why the classical RT [k] mixed finite elements introduced by Raviart and Thomas in [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF] do not lead to an explicit scheme in the case of an anisotropic medium and we propose instead the use of new mixed finite elements. The lowest order element is presented in section 2.2, the higher orders in section 2.3. Section 3 is the main part of the paper and is concerned with the mixed approximation of the elliptic problem which is nothing but the stationary problem associated to the evolution problem presented in section 2. This analysis will be used to study an elliptic projection operator that will be useful for the analysis of the approximation of the evolution problem. We shall show in section 3.1 why the analysis of the new element does not fit the classical theory. That is why we develop in section 3.2 a novel abstract theory leading to nonclassical error estimates. In section 3.3, we show that the mixed finite elements that we have constructed enter this framework and error estimates are given. Section 4 is devoted to relating the error estimates on the time domain solution to the error estimates obtained in the previous section on the elliptic problem. This part of the analysis is inspired by [START_REF] Dupont | L 2 -estimates for Galerkin methods for second order hyperbolic equations[END_REF] and essentially relies on energy estimates.

Although in this paper we restrict ourselves to the 2D case, the extension of this work to three dimensions is straightforward.

Presentation of the new mixed finite elements.

The model problem:

The anisotropic wave equation. Let Ω be a bounded domain of R 2 and let A(x) be a positive definite symmetric matrix satisfying

A(x) ξ • ξ ≥ α |ξ| 2 , α > 0, ∀ξ ∈ R 2 , a.e. x ∈ Ω. (2.1)
In practice, Ω will be the simple geometrical domain (let us say a rectangle) we referred to in the introduction. We consider the scalar evolution problem

⎧ ⎪ ⎨ ⎪ ⎩ Find u : [0, T ] -→ H 1 0 (Ω) such that ∂ 2 u ∂t 2 -div (A -1 (x)∇u) = f, f ∈ C 0 (0, T ; L 2 (Ω)), (2.2)
subject to the initial conditions

u(t = 0) = u 0 ∈ H 1 (Ω) ; ∂u ∂t (t = 0) = u 1 ∈ L 2 (Ω).
Remark 2.1. In this section we are not interested in the regularity of u with respect to time variable t. This is why we simply write u : [0, T ] -→ H 1 0 (Ω). Remark 2.2. The homogeneous Dirichlet condition on ∂Ω has been considered for simplicity only. This has nothing to do with a Neumann condition on the boundary of an obstacle interior to Ω that will be taken into account with the ficitious domain method as explained in the introduction.

Let

p = A -1 (x)∇u and v = ∂u ∂t . Substituting into (2.2) yields ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂v ∂t -div p = f, A ∂p ∂t -∇v = 0 (2.3) with initial conditions p(0) = p 0 = A -1 (x)∇u 0 ; v(0) = v 0 = u 1 . (2.4)
Remark 2.3. Throughout this paper, the generic point of R 2 will be indifferently denoted by the letter x or the couple (x, y) (depending on the context).

A mixed formulation associated to (2.3) is given by the following problem:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p, v) : [0, T ] -→ X × M ≡ H(div ; Ω) × L 2 (Ω) such that d dt a(p, q) + b(v, q) = 0 ∀q ∈ X, d dt (v, w) -b(w, p) = (f, w) ∀w ∈ M, (2.5)
where

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a(p, q) = Ω A(x) p • q dx ∀(p, q) ∈ X × X, b(w, q) = Ω w div q dx ∀(w, q) ∈ M × X, (f, w) = Ω f w dx ∀w ∈ M.
(2.6)

The bilinear form a(

•, •) (resp., b(•, •)) is continuous on H × H (H = (L 2 (Ω)) 2 ) (resp., on X × M ). The bilinear form a(•, •) (resp., b(•, •))
thus defines a linear continuous operator A : H → H by Ap, q H ×H = a(p, q) (resp., B : X → M by Bq, w M ×M = b(w, q)). They satisfy the following properties (see, for instance, [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]):

(i)

The continuous inf-sup condition:

∃c > 0 / ∀ w ∈ M, ∃ q ∈ X / b(w, q) ≥ c w M q X .
(ii) The coercivity of the form a(•, •) on V ≡ KerB :

∃α > 0 / ∀ p ∈ V, a(p, p) ≥ α p 2 X .
(2.7)

In the following, we consider only the semidiscretization in space of this problem, but with the specific objective of obtaining a discretization for which mass lumping is possible.

Presentation of the Q div

1

-Q 0 finite element in the lowest order. We suppose now that Ω is a union of rectangles in such a way that we can consider a regular mesh (T h ) with square elements K of edge h > 0. We introduce the following approximation spaces:

⎧ ⎪ ⎨ ⎪ ⎩ X h = q h ∈ X / ∀K ∈ T h , q h | K ∈ X , M h = w h ∈ M / ∀K ∈ T h , w h | K ∈ M , (2.8)
where X (resp., M ) denotes a finite dimensional space of vector (resp., scalar) functions. The discrete problem associated to (2.5) and (2.4) is

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p h , v h ) : [0, T ] → X h × M h such that d dt a(p h , q h ) + b(v h , q h ) = 0 ∀q h ∈ X h , d dt (v h , w h ) -b(w h , p h ) = (f, w h ) ∀w h ∈ M h , (2.9)
subject to the initial conditions

p h (0) = p 0,h ; v h (0) = v 1,h .
The usual choice consists in taking for X the lowest order Raviart-Thomas element (see [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]), X = RT [0] = P 10 × P 01 , (2.10) and for M the piecewise constants:

M = Q 0 . (2.11)
Here, P k is the space of polynomials of degree not greater than k and P kl is defined by

P kl = ⎧ ⎨ ⎩ p(x, y) | p(x, y) = i≤k,j≤l a ij x i y j ⎫ ⎬ ⎭ .
This choice, however, does not lead to an explicit scheme when one considers the evolution problem of anisotropic waves corresponding to (2.9). We introduce here

B N 1 = {τ i } N 1 i=1 , B N 2 = {φ i } N 2
i=1 the bases of X h and M h , resp., where N 1 = dimX h and N 2 = dimM h . We denote then by [P ] = (P 1 , . . . , P N 1 ) and [V ] = (V 1 , . . . , V N 2 ) the coordinates of p h and v h in the bases B N 1 and B N 2 . In these bases, problem (2.9) can be written in the following form:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (P, V ) : [0, T ] -→ R N 1 × R N 2 such that M p dP dt + C T V = 0, M v dV dt -CP = F, + initial conditions with (i) (M p ) i,j = (Aτ i , τ j ) (L 2 (Ω)) 2 , 1 ≤ i, j ≤ N 1 , (ii) (M v ) i,j = (φ i , φ j ) L 2 (Ω) , 1 ≤ i, j ≤ N 2 , (iii) (C) i,j = (φ i , div τ j ) L 2 (Ω) , 1 ≤ i ≤ N 2 , 1 ≤ j ≤ N 1 , (iv) (F ) j = (f, φ j ) L 2 (Ω) , 1 ≤ j ≤ N 2 ;
(2.12)

C T denotes the transpose of C. If we use a centered finite difference approximation for the time discretization, the solution at each time step is obtained by inverting the mass matrices M v and M p . Although they are symmetric and sparse, this inversion can become costly (for large systems). To avoid that, we want to reduce them to diagonal (or block diagonal) matrices by using a mass lumping technique (see [START_REF] Cohen | Gauss point mass lumping schemes for Maxwell's equations[END_REF][START_REF] Tordjman | Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes[END_REF]). This consists of using adequate quadrature formulas to approximate the integrals in (2.12 (i) and (ii)). One can remark that, M h being chosen discontinuous, the matrix M v does not really need to be mass lumped (it is already diagonal here); so, we focus our attention on the mass lumping of M p . Let now (τ I ) be the RT [0] global basis functions. In Figure 2.1, we represent the four local basis vector functions in RT [0] ( K), K being the reference element. Each global basis function τ I is associated to an edge of the mesh and is supported by the two elements adjacent to this edge; see Figure 2.2. The restriction of τ I on one element K of its support corresponds to a local basis function, denoted τ i , on the reference element K. In fact τ i corresponds to τ x i (resp., τ y i ) when τ I is associated to a vertical (resp., horizontal) edge of the mesh.

τ y 1 = 0 1 -y , τ y 2 = 0 y τ x 1 = 1 -x 0 , τ x 2 = x 0 K τ y 2 τ x 1 τ y 1 τ x 2 Fig. 2.1. Local vector basis functions in RT [0] ( K).

Now consider the integral

a(τ I , τ J ) = K∈T h K Aτ I | K • τ J | K dx = K∈T h mes(K) K A τ i • τ j d x. (2.13) If we use the quadrature formula K f ( x)d x = N q=1 ω q f ( M q ),
with ( M q ) q=1,...,N being the quadrature points and (ω q ) q=1,...,N the associated weights to approximate (2.13), we obtain

K A τ i • τ j d x ≈ q ω q A τ i ( M q ) • τ j ( M q ).
This would lead to a diagonal matrix if

A τ i ( M q ) • τ j ( M q ) = 0 ∀ i = j. (2.14)
For the isotropic wave equation (A = Id), it is easy to find a quadrature formula satisfying (2.14); take, for instance, the lowest order Gauss-Lobatto quadrature formula using the vertices of the element K as quadrature points. But for the anisotropic wave equation, there is no usual quadrature formula that satisfies (2.14). The difficulty is that for two orthogonal edges of the element K, the functions A τ y 1 and τ x 1 are no longer orthogonal.

The alternative solution that we propose consists of changing the approximation space X. Let

X = Q 1 × Q 1 , (2.15)
where Q 1 is the space of piecewise bilinear functions. This element was initially introduced by Nédélec in [START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF]. In Figure 2.3 we represent the eight local basis vector functions on the reference element K. We will call this choice (2.15) combined with (2.11)

the Q div 1 -Q 0 element. K τ y 4 τ y 3 τ y 1 τ x 4 = ((1 -x)y, 0) t , τ x 3 = (xy, 0) t τ y 2 τ x 1 τ x 4 τ x 2 τ x 3 M 1 M 2 M 3 M 4 τ y 1 = (0, (1 -x)(1 -y)) t , τ y 2 = (0, x(1 -y)) t τ y 4 = (0, (1 -x)y) t , τ y 3 = (0, xy) t τ x 1 = ((1 -x)(1 -y), 0) t , τ x 2 = (x(1 -y), 0) t Fig. 2.3. The local basis functions in X = Q 1 × Q 1 .
The lowest order Gauss-Lobatto quadrature formula now satisfies (2.14) with the new basis functions; the key point is that the four quadrature points M j coincide with the locations in space of the degrees of freedom and the new basis functions satisfy

τ x i (M j ) = δ ij (1, 0) t , τ y i (M j ) = δ ij (0, 1) t .
Note that, over each element, there are four quadrature points but eight degrees of freedom. For more details on quadrature formulae and mass lumping techniques we refer to the work of Cohen, Joly, and Tordjman for the acoustic wave equation [START_REF] Tordjman | Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes[END_REF][START_REF] Cohen | Construction and analysis of higher order finite elements with mass lumping for the wave equation[END_REF][START_REF] Cohen | Construction analysis of fourth-order finite difference schemes for the acoustic wave equation in nonhomogeneous media[END_REF] and to Cohen and Monk [START_REF] Cohen | Gauss point mass lumping schemes for Maxwell's equations[END_REF] and Elmkies and Joly [START_REF] Elmkies | Éléments finis et condensation de masse pour les équations de Maxwell: le cas 2D[END_REF] for Maxwell's equations. Remark 2.4. With Gauss-Lobatto's rule and the Q 1 × Q 1 element presented above, the mass matrix a(p h , q h ) is block diagonal, the size of each block being equal to the number of degrees of freedom at each vertex, which is four, see Figure 2.4; thus after an inversion of local 4 × 4 matrices, we obtain an explicit scheme. In the particular case of the isotropic wave equation, A = Id, the mass matrix becomes 2 × 2 block diagonal, because the basis functions τ x i and τ y j are orthogonals. Remark 2.5. From the physical point of view, having four degrees of freedom for p at each vertex is coherent with the actual discontinuities of p x or p y that occur when A is dicontinuous across the edges of the mesh.

Before analyzing this element, let us show how it can easily be extended to higher order.

1 0 0 0 1 1 1 0000 1111 i + 1 (τ y ij ) + j j + 1 i i -1 i -1 i i + 1 j + 1 j -1 j j -1 (τ x ij ) + (τ x ij ) - (τ y ij ) - Fig. 2.4. The Q 1 × Q 1 element.

Extension to higher order and mass lumping.

The natural generalization of the lowest order element, presented in the previous section, consists of taking

⎧ ⎨ ⎩ X h = {q h ∈ X / ∀K ∈ T h , q h | K ∈ Q k+1 × Q k+1 } , M h = {w h ∈ M / ∀K ∈ T h , w h | K ∈ Q k } , ( 2.16)
and we will call it the Q div k+1 -Q k element. This choice still satisfies our requirement with respect to mass lumping. We represent in Figures 2.5 and 2.6 the degrees of freedom for k = 1 and k = 2; the exact locations of the quadrature points and the associated weights are given in Appendix A. We also indicate in these figures the number of degrees of freedom per node. The locations of the degrees of freedom correspond to tensor products of 1D quadrature points associated to Gauss-Lobatto (for p h ) or Gauss-Legendre (for v h ) quadrature formulas. Indeed, following the approach of Cohen and Monk [START_REF] Cohen | Gauss point mass lumping schemes for Maxwell's equations[END_REF] and Tordjman [START_REF] Tordjman | Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes[END_REF], we obtain mass lumping by approximating the integrals in (2.12(i), (ii)) with adequate quadrature formulas, for which the position of the quadrature points coincide with the position of the degrees of freedom. For our elements they consist of • using the Gauss-Lobatto quadrature formulas for the approximation of the M p matrix. The resulting matrix is now block diagonal. Each block is associated to one quadrature point and its dimension is equal to the number of degrees of freedom at this point. (The "worst" case concerns the vertices of the elements K, where the dimension of the local block is 4 × 4.) • using the Gauss-Legendre quadrature formula for the approximation of the M v matrix, so that the resulting matrix is diagonal. Remark 2.6. The Gauss-Legendre quadrature formula is exact for the integration of the product v h w h which is in Q 2k . On the other hand, considering for simplicity the case where A is constant, the product of Ap h .q h belongs to Q 2k+2 . However, following [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Tordjman | Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes[END_REF][START_REF] Baker | The effect of quadrature errors on finite element approximations for the second-order hyperbolic equations[END_REF] one can prove that one does not lose any order of accuracy provided that the quadrature formula is exact in Q k which is the case with the Gauss-Lobatto formula we use.

Remark 2.7. The Raviart-Thomas kth order approximation consists of the choice

XRT = RT [k] = P k+1,k × P k,k+1 and M = Q k .
It is clear that the new approximate space X h contains the space X RT h . Thus, we have enriched the p-approximation, while keeping the same approximation space for v.

Analysis of the new mixed finite element for an elliptic problem.

Following [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[END_REF], we will study in this section the mixed approximation of the elliptic problem which is in fact the stationary problem associated to the evolution problem (2.5). Actually, we give in section 3.2 an abstract result for a class of elliptic problems posed in a more general framework and we show in section 3.3 that the model problem enters this framework. This analysis will be used in section 3.4 to study an elliptic projection operator that will be useful for the analysis of the approximation of the evolution problem done in section 4.

The elliptic problem.

The elliptic problem we consider here is

⎧ ⎨ ⎩ Find u ∈ H 1 0 (Ω) such that -div (A -1 (x)∇u) = f, f ∈ L 2 (Ω). (3.1)
We know that (3.1) admits a unique solution u ∈ H 1 0 (Ω) and that there exists c > 0 such that

u H 1 (Ω) ≤ c f L 2 (Ω) . (3.2)
As for the time dependent problem, we set p = A -1 (x)∇u, (3.3) and this gives

-div p = f. (3.4)
The mixed formulation associated to (3.3)

and (3.4) is ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p, u) ∈ X × M = H(div ; Ω) × L 2 (Ω) such that a(p, q)+ b(u, q) = 0 ∀q ∈ X, b(w, p) = -(f, w) ∀w ∈ M, (3.5)
where a(•, •) and b(•, •) are defined by (2.6) and satisfy properties (2.7). As proven in [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF], there exists a unique solution (p, u) in X × M to problem (3.5), where u is also the solution of the initial problem (3.1). (In fact the abstract result yields the uniqueness of u only in M/ Ker B t , but it is easy to check that for the divergence operator Ker B t = {0}.) For the approximation of this problem, we again consider the finite dimensional spaces X h and M h defined by (2.8). The discrete problem associated to

(3.5) is ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p h , u h ) ∈ X h × M h such that a(p h , q h )+ b(u h , q h ) = 0 ∀q h ∈ X h , b(w h , p h ) = -(f, w h ) ∀w h ∈ M h . ( 3.6)
The elliptic problem (3.5) and its approximation (3.6) have been studied by several authors (see, for example, [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF][START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]), and we know that it admits a unique solution

(p h , u h ) in X h × M h with the convergence property (p h , u h ) -→ (p, u) ∈ X × M (3.7)
when the following assumptions are satisfied:

(i)
The uniform discrete inf-sup condition:

∃c > 0 independent of h such that ∀ w h ∈ M h , ∃ q h ∈ X h / b(w h , q h ) ≥ c w h M q h X ; (ii) The coercivity of the form a(•, •) on V h : ∃α > 0 independent of h such that ∀ p h ∈ V h , a(p h , p h ) ≥ α p h 2 X ,
where

V h = KerB h = {q h ∈ X h / b(w h , q h ) = 0 ∀ w h ∈ M h } . (3.8) K 1 K 2 (0, 0) (-h, 0) (h, 0) (h, h) (0, h) (-h, h) Fig. 3.1. The function f h .
These assumptions are satisfied by the couple of spaces (X h , M h ) defined in (2.8), X being the lowest order Raviart-Thomas element (2.10) and M the piecewise constants (2.11); see [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]. With our new choice for X it is easy to verify that property (3.8(i)) is still true but we no longer satisfy (3.8(ii)). Indeed, consider the function f h given by (see Figure 3.1)

f h | K 1 = ⎛ ⎝ 1 - x h 1 -2 y h 0 ⎞ ⎠ and f h | K 2 = ⎛ ⎝ 1 + x h 1 -2 y h 0 ⎞ ⎠ .
We can easily see that

f h ∈ V h (if M = Q 0 ) and a(f h , f h ) -→ h→0 0, f h X -→ h→0 2 3 ,
so that we cannot expect to verify (3.8(ii)). However, we will prove (see section 3) that this choice gives a good approximate solution and we will show a new convergence result. Remark 3.1. In order to preserve (3.8(ii)) one could change the approximation space M . The natural choice is

M = P 1 ,
the key point being that the divergence operator sends Q 1 × Q 1 into P 1 . 1 We have eliminated this choice because it is rather expensive in terms of computational time and memory requirements.

An abstract result.

The first point in our new theory is that we need to introduce a third Hilbert space. More precisely, let M, X, H be three Hilbert spaces with

X ⊂ H, X dense in H and |•| H ≤ • X . (3.9)
The reader can have in mind that for our application we shall have

H = (L 2 (Ω)) 2 , X = H(div; Ω) and M = L 2 (Ω).
Let a(•, •) and b(•, •) be two continuous bilinear forms in H × H and M × X. In the same way as in section 2.1, the bilinear form a(•, •) defines an operator A in L(H), such that a(p, q) = (Ap, q) H ∀ (p, q) ∈ H × H and the bilinear form b(•, •) defines an operator B : X → M (and its transpose B t : M → X ) such that

Bp, w M ×M = p, B t w X×X = b(w, p) ∀(p, w) ∈ X × M.
The kernels of B and B t are defined as follows:

⎧ ⎨ ⎩ V ≡ Ker B = {p ∈ X / b(w, p) = 0 ∀w ∈ M } , Ker B t = {w ∈ M / b(w, p) = 0 ∀p ∈ X} .
We make the assumptions

(i) ∃ c > 0 / ∀ w ∈ M, ∃ q ∈ X / b(w, q) ≥ c w M/Ker B t q X , (ii) ∃ α > 0 / ∀ p ∈ V, a(p, p) ≥ α p 2 X , (3.10)
where the norm in the quotient space is defined by

w M/Ker B t = inf w 0 ∈Ker B t w + w 0 M .
We shall identify the quotient space M/Ker B t with the orthogonal complement of

Ker B t M/Ker B t ≡ (Ker B t ) ⊥ ≡ w ∈ M / (v, w) M = 0 ∀ v ∈ Ker B t .
We are interested in the numerical approximation of the solution (p, u) to the following problem:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p, u) ∈ X × M such that a(p, q)+ b(u, q) = 0 ∀q ∈ X, b(w, p) = -f, w ∀w ∈ M, (3.11)
with f ∈ M the dual space of M . Under these assumptions, we have the following classical result (see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]). Theorem 3.1. For all f ∈ Im B, problem (3.11) has a unique solution (p, u) in X × (M/Ker B t ). Moreover,

u M/Ker B t + p X ≤ C f M .
Suppose X h ⊂ X and M h ⊂ M are finite dimensional approximation spaces. We consider then the approximate problem

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p h , u h ) ∈ X h × M h such that a(p h , q h )+ b(u h , q h ) = 0 ∀q h ∈ X h , b(w h , p h ) = -f, w h ∀w h ∈ M h . (3.12) We set ⎧ ⎨ ⎩ V h (f ) = {q h ∈ X h / b(w h , q h ) = -f, w h ∀w h ∈ M h } = B -1 h (f ), V h = V h (0) = Ker B h (see below for the meaning of B h ) (3.13)
and make the following hypotheses:

(H0) ∀f ∈ Im B, V h (f ) = ∅.
(H1) Orthogonal decomposition of X h :

X h = X s h ⊕ X r h (p h = p s h + p r h ) , X r h ⊂ V h ; ∀(p s h , p r h ) ∈ X s h × X r h (p s h , p r h ) H = 0.
(H2) "Strong" discrete uniform inf-sup condition:

there exists a constant c > 0, independent of h, such that

∀ w h ∈ M h , ∃ q s h ∈ X s h / b(w h , q s h ) ≥ c w h M/Ker B t h q s h X .
(H3) "Weak" coercivity:

there exists a constant α > 0, independent of h, such that

∀ p h ∈ V h , a(p h , p h ) ≥ α p s h 2 X + |p r h | 2 H . (H4) Approximation properties ⎧ ⎪ ⎨ ⎪ ⎩ lim h→0 inf q s h ∈X s h p -q s h X = 0 ∀p ∈ X, lim h→0 inf w h ∈M h u -w h M = 0 ∀u ∈ M. Remark 3.2. Note that, since X is dense in H, (H4) implies that lim h→0 inf q s h ∈X s h |p -q s h | H = 0 ∀p ∈ H.
We can introduce, as for the continuous problem, an operator B h from X h to M h defined by

(B h p h , w h ) M ×M = p h , B t h w h X×X = b(p h , w h ) ∀p h ∈ X h ∀w h ∈ M h
and the kernels of B h (see (3.13)) and of its transpose

KerB t h = {w h ∈ M h / b(p h , w h ) = 0 ∀p h ∈ X h } . Remark 3.3.
It may be more convenient to characterize hypothesis (H0) with one of the following equivalent statements:

(H0ii) ∀p ∈ X, ∃p h ∈ X h such that b(p -p h , w h ) = 0 ∀w h ∈ M h ; (H0iii) Ker B t h = Ker B t ∩ M h ⊂ Ker B t .
Remark 3.4. We call hypothesis (H2) "strong" the existence of a q s h in X s h instead of in X h , which would be the classical assumption; more precisely X s h ⊂ = X h . On the contrary, hypothesis (H3) is "weaker" in the sense that

∀ p h ∈ X s h or p h ∈ X r h , p s h 2 X + |p r h | 2 H ≤ p h 2 X .
Theorem 3.2. Under hypotheses (H0)-(H4), problem (3.12) admits a unique solution such that

(p h = p s h + p r h , u h ) ∈ X h × (M h /Ker B t h
), and the following convergence result holds:

• (p s h , u h ) → (p, u) in X × M, • p r h → 0 in H.

More precisely,

|p r h | H + p -p s h X + u -u h M/Ker B t h ≤ C inf q s h ∈X s h p -q s h X + inf w h ∈M h u -w h M + inf z s h ∈X s h |Ap -z s h | H .
Proof. First note that hypothesis (H0) and nonuniform discrete coercivity on the kernel Ker B h , i.e.,

∃α h > 0 ∀p h ∈ V h a(p h , p h ) ≥ α h p h 2 X , (3.14)
ensure the existence and uniqueness of the solution (p h , u h ) of the discrete problem (3.12) in X h × (M h /Ker B t h ). Since in finite dimensions all norms are equivalent, the nonuniform discrete coercivity (3.14) on the kernel Ker B h is a consequence of (H3).

It remains to prove the error estimate. The convergence result will follow from assumption (H4) and Remark 3.2.

The second equation of (3.12) shows that p h ∈ V h (f ) (cf. (3.13)). If we also take q h ∈ V h (f ), the difference is in the kernel, i.e., p hq h ∈ V h . We can write

a(p h -q h , p h -q h ) = a(p -q h , p h -q h ) + a(p h -p, p h -q h ). (3.15)
Taking the difference between the first equation of the continuous (3.11) and the discrete (3.12) problem we have

a(p h -p, p h -q h ) = b(u -w h , p h -q h ) ∀w h ∈ M h . (3.16)
Using (H1) in (3.16), we can write, with obvious notation,

a(p h -p, p h -q h ) = b(u -w h , p s h -q s h ) + b(u, p r h -q r h ). (3.17)
We have used the fact that b(w h , p r hq r h ) = 0, since X r h ⊂ V h . Now, since we only have a "weak" coercivity, we want to keep only terms leading to H-norm on the X r h part; thus we want to eliminate b(u, p r hq r h ) (which would give p r hq r h X from the continuity of b). Then, by (3.11),

a(p h -p, p h -q h ) = b(u -w h , p s h -q s h ) -a(p, p r h -q r h ) = b(u -w h , p s h -q s h ) -(Ap, p r h -q r h ) H .
(3.18)

We choose now

q h = q s h ∈ V s h (f ) = V h (f ) ∩ X s h , q r h = 0. ( 3
a(p h -q s h , p h -q s h ) = a(p -q s h , p h -q s h ) + b(u -w h , p s h -q s h ) -(Ap, p r h ),
or, by using the orthogonality of X r h and X s h from (H1),

a(p h -q s h , p h -q s h ) = a(p -q s h , p h -q s h ) + b(u -w h , p s h -q s h ) -(Ap -z s h , p r h ) H ∀ z s h ∈ X s h . (3.20)
Furthermore, from the inequality

|p h -q s h | H ≤ p s h -q s h X + |p r h | H
and from hypothesis (H1) and (H3), (3.20) leads to

α( p s h -q s h 2 X + |p r h | 2 H ) ≤ a |p -q s h | H ( p s h -q s h X + |p r h | H ) + b u -w h M p s h -q s h X + |Ap -z s h | H |p r h | H .
We deduce the existence of a constant C depending on f, p, a , b , and α such that

∀(q s h , w h , z s h ) ∈ V s h (f ) × M h × X s h , p s h -q s h X + |p r h | H ≤ C (|p -q s h | H + u -w h M + |Ap -z s h | H ) .
This gives, using (3.9),

p -p s h X ≤ p -q s h X + p s h -q s h X ∀q s h ∈ V s h (f ) ≤ (1 + C) inf q s h ∈V s h (f ) p -q s h X + inf w h ∈M h u -w h M + inf z s h ∈X s h |Ap -z s h | H , |p r h | H ≤ C inf q s h ∈V s h (f ) p -q s h X + inf w h ∈M h u -w h M + inf z s h ∈X s h |Ap -z s h | H .
To conclude, let us recall that the inf-sup condition (H2) implies (cf. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]) that inf

q s h ∈V s h (f ) p -q s h X ≤ c 1 inf q s h ∈X s h p -q s h X .
Finally, it remains to prove estimates for uu h M/Ker B t h . Let us subtract the first equation of (3.12) from the first equation of (3.11). We get

a(p -p h , q h ) + b(u -u h , q h ) = 0 ∀q h ∈ X h , so that, for any w h in M h , b(u h -w h , q h ) = a(p -p h , q h ) + b(u -w h , q h ) ∀q h ∈ X h . Choosing now q h = q s h implies that b(u h -w h , q s h ) = a(p -p h , q s h ) + b(u -w h , q s h ) ∀q s h ∈ X s h .
Using this and the inf-sup condition, we have

u h -w h M/Ker B t h ≤ 1 c sup q s h ∈X s h b(u h -w h , q s h ) q s h X ≤ 1 c sup q s h ∈X s h a(p -p h , q s h ) + b(u -w h , q s h ) q s h X ≤ 1 c ( a |p -p h | H + b u -w h M ) .
It follows from the triangle inequality that

u -u h M/Ker B t h ≤ C inf w h ∈M h u -w h M + |p r h | H + |p -p s h | H .
Remark 3.5. One could easily remark that the space X s h satisfies all assumptions (H0) to (H4). Why not use X s h as approximate space instead of X h ? In fact, we should imagine the case in which either we cannot characterize the space X s h or we prefer X h for other reasons. This occurs in particular for the evolution problem, where we prefer X h in order to achieve mass lumping.

Application to the model problem. For problem (3.5) we have

H = (L 2 (Ω)) 2 ; X = H(div ; Ω); M = L 2 (Ω).
The operator B is the divergence operator, which is surjective from X onto M , so that Ker B t = {0}. From Remark 3.3, it follows that hypothesis (H0) is equivalent to the surjectivity of B h from X h onto M h . This is true for the classical Raviart-Thomas approximations [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]; i.e., B h is surjective from X RT h onto M h . Since the new space X h contains the space X RT h , as noted in Remark 2.7, it is straightforward to check (H0); therefore, Ker B t h = {0}.

Approximation with the lowest order finite element

Q div 1 -Q 0 .
Let us begin by checking hypothesis (H1). We take as approximation spaces

⎧ ⎨ ⎩ X h = {q h ∈ X / ∀K ∈ T h , q h | K ∈ Q 1 × Q 1 } , M h = {w h ∈ M / ∀K ∈ T h , w h | K ∈ Q 0 } .
We define X s h as the lowest order Raviart-Thomas element RT [0] defined in (2.10)

X s h = q h ∈ X / ∀K ∈ T h , q h | K ∈ RT [0] = P 10 × P 01 ⊂ X h , dim RT [0] = 4.
In order to describe its orthogonal complement X r h , denote by -(i, j), the node of T h with coordinates (ih, jh), -(i + 1 2 , j), the horizontal side joining the nodes (i, j) and (i + 1, j), -(i, j + 1 2 ), the vertical side joining the nodes (i, j) and (i, j + 1). Over the square K, the four basis functions of RT [0] can be written as (see Figure 3.2 for an illustration)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ϕ i,j+ 1 2 = ⎛ ⎝ ϕ i,j+ 1 2 (x) 0 ⎞ ⎠ , ϕ i+ 1 2 ,j = ⎛ ⎝ 0 ϕ i+ 1 2 ,j (y) ⎞ ⎠ , ϕ i+1,j+ 1 2 = ⎛ ⎝ ϕ i+1,j+ 1 2 (x) 0 ⎞ ⎠ , ϕ i+ 1 2 ,j+1 = ⎛ ⎝ 0 ϕ i+ 1 2 ,j+1 (y) ⎞ ⎠ , K (i + 1, j + 1) ϕ i+1,j+1/2 (i, j + 1) ϕ i,j+1/2 ϕ i+1/2,j+1 ϕ i+1/2,j (i, j) (i + 1, j) Fig. 3.2. The RT [0] element.
where

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ϕ i,j+ 1 2 (x) = (x i+1 -x) h ∈ P 10 , ϕ i+ 1 2 ,j (y) = (y j+1 -y) h ∈ P 01 , ϕ i+1,j+ 1 2 (x) = (x -x i ) h ∈ P 10 , ϕ i+ 1 2 ,j+1 (y) = (y -y j ) h ∈ P 01 .
It is then easy to prove the following lemma.

Lemma 3.3. If X r h = (X s h ) ⊥ , the space X r h | K of the restrictions to K of elements of X r
h can be generated from the four functions: .

ψ i+ 1 2 ,j = ⎛ ⎜ ⎝ 0 x -x i+ 1 2 ϕ i+ 1 2 ,j (y) ⎞ ⎟ ⎠ , ψ i,j+ 1 2 = ⎛ ⎝ y -y j+ 1 2 ϕ i,j+ 1 2 (x) 0 ⎞ ⎠ , ψ i+ 1 2 ,j+1 = ⎛ ⎜ ⎝ 0 x -x i+ 1 2 ϕ i+ 1 2 ,j+1 (y) ⎞ ⎟ ⎠ , ψ i+1,j+ 1 2 = ⎛ ⎝ y -y j+ 1 2 ϕ i+1,j+ 1 2 (x) 0 ⎞ ⎠ , where x i+ 1 2 = (i + 1 2 )h and y j+ 1 2 = (j + 1 2 )h. Moreover, X r h ⊂ V h . Remark 3.6. Note that the dimension of Q 1 × Q 1 (=8) is the sum of the dimen- sions of X s h | K (= dim RT [0] = 4) and X r h | K (=4

Note now that p

h ∈ V h and X r h ⊂ V h imply p s h ∈ V h ∩ X s h . Hence it follows (cf. [6]) that ⎧ ⎨ ⎩ (i) divp s h ∈ M h , (ii) (divp s h , w h ) = 0 ∀w h ∈ M h .

Recalling (i), we can take

w h = div p s h in (ii), which gives div p s h = 0. The claim is thus proved.

We can now apply Theorem 3.2 to the approximation problem (3.6), and by using the usual interpolation results (cf. [START_REF] Roberts | Mixed and hybrid methods[END_REF]) we obtain the following theorem (here we have Ker B t h = {0}). Theorem 3.5. The problem (3.6) admits a unique solution

(p h , u h ) ∈ X h × M h , which satisfies • (p s h , u h ) -→ (p, u) in H(div ; Ω) × L 2 (Ω), • p r h -→ 0 in (L 2 (Ω)) 2 .
Furthermore, if the solution is assumed to be sufficiently regular, i.e., (p, u)

∈ (H 1 (Ω)) 2 × H 1 (Ω), div p ∈ L 2 (Ω), and Ap ∈ (H 1 (Ω)) 2 , then p r h (L 2 (Ω)) 2 + p -p s h H(div;Ω) + u -u h L 2 (Ω) ≤ Ch |u| H 1 + |p| (H 1 ) 2 + |div p| H 1 + |Ap| (H 1 ) 2 ,
where

|.| H 1 (resp., |.| (H 1 ) 2 ) denotes the usual semi-norm in H 1 (Ω) (resp., (H 1 (Ω)) 2 ). Remark 3.

Let us now consider the case of the isotropic wave equation, which corresponds in taking A(x) in (3.1) as a diagonal matrix (A d (x)). For the discretization we assume that A d (x) is approximated by piecewise constant functions per element. It is then easy to prove that the approximate problem

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p h , u h ) ∈ X h × M h such that a d (p h , q h )+ b(u h , q h ) = 0 ∀q h ∈ X h , b(w h , p h ) = -f, w h ∀w h ∈ M h (3.21)
admits a unique solution (p h = p s h +p r h , u h ) with p r h ≡ 0 and (p s h , u h ) being the solution of the following problem:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p s h , u h ) ∈ X s h × M h such that a d (p s h , q s h )+ b(u h , q s h ) = 0 ∀q s h ∈ X s h , b(w h , p s h ) = -f, w h ∀w h ∈ M h
where a d (•, •) is given by (2.6) after replacing A(x) by A d (x). Indeed, in this particular case,

⎧ ⎨ ⎩ a d (p r h , q s h ) = 0 ∀q s h ∈ X s h , a d (p s h , q r h ) = 0 ∀q r h ∈ X r h .
By using the fact that X r h ⊂ V h , we can decompose (3.21) into two independent problems in X s h and X r h :

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p s h , u h ) ∈ X s h × M h such that a d (p s h , q s h )+ b(u h , q s h ) = 0 ∀q s h ∈ X s h , b(w h , p s h ) = -f, w h ∀w h ∈ M h , and ⎧ ⎨ ⎩ Find p r h ∈ X r h such that a d (p r h , q r h ) = 0 ∀q r h ∈ X r h . (3.22)
It is obvious then from (3.22) that p r h ≡ 0. This remark is no longer true when A(x) is not diagonal.

Approximation with higher order finite elements,

Q div k+1 -Q k .
We have seen in section 2.3 that the natural generalization of the lowest order element consists of taking

X = X k = Q k+1 × Q k+1 and M = M k = Q k ;
and thus to introduce the spaces (2.16)

⎧ ⎨ ⎩ X h = {q h ∈ X / ∀K ∈ T h , q h | K ∈ X k } , M h = {w h ∈ M / ∀K ∈ T h , w h | K ∈ M k } .
We are going to show that we can apply the abstract theory of section 3.2 to these spaces, with the main difficulty lying in the construction of an orthogonal decomposition of X h such that assumptions (H0) to (H4) are satisfied. In fact we shall deduce such a global decomposition from a local decomposition of the space X, considered as a subspace of L 2 (K), where K is a single element. First, we recall that

RT [k] = P k+1,k × P k,k+1 ,
which obviously satisfies the inclusion

RT [k] ⊂ X k .
Let us define Ψ k as the orthogonal complement in X k of RT [k] (for the inner product in L 2 (K)):

Ψ k = ψ ∈ X / K ψ φdx = 0 ∀φ ∈ RT [k] .
Note that

dimΨ k = 2(k + 2).
The main property of the space Ψ k is the following. Lemma 3.6. For any

ψ in Ψ k and v in M k , K div ψ v dx = 0. (3.23)
Proof. For simplicity we consider the reference element K = [0, 1] × [0, 1]. (i) We begin by a characterization of Ψ k . Let σ k be the polynomial in one variable of degree k + 1 such that

P k+1 = P k ⊕ [σ k ] .
Equivalently, σ k is defined, up to a multiplicative constant, by

⎧ ⎪ ⎨ ⎪ ⎩ 1 0 σ k (x)p(x)dx = 0 ∀p ∈ P k , σ k ∈ P k+1 , σ k / ≡ 0.
Then, we claim that

Ψ k = Ψ k ≡ ⎧ ⎨ ⎩ ⎛ ⎝ p x (x)σ k (y) p y (y)σ k (x) ⎞ ⎠ , (p x , p y ) ∈ P k+1 × P k+1 ⎫ ⎬ ⎭ . Indeed, as dim Ψ k = 2(k + 2), it is sufficient to check that Ψ k is orthogonal to RT [k] .
To see that, we first remark that (see [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]) RT [k] is generated by vector fields of the form As q x and q y are in P k , by the definition of σ k ,

φ(x, y) = ⎡ ⎣ τ x (x)q x (y) τ y (y)q y (x) ⎤ ⎦ with (τ x , τ y ) ∈ P k+1 × P k+1 , (q x , q y ) ∈ P k × P k . Let us consider ψ(x, y) = ⎡ ⎣ p x (x)σ k (y) p y (y)σ k (x) ⎤ ⎦ with (p x , p y ) ∈ P k+1 × P k+1 . (3.24) We have (φ, ψ) L 2 (K) = 1 0 τ x (x)p x (x)dx 1 0 σ k (y)q x (y)dy + 1 0 σ k (x)q y (x)dx
1 0 σ k (y)q x (y)dy = 1 0 σ k (x)q y (x)dx = 0 ⇒ (φ, ψ) L 2 (K) = 0. K T 1 T 3 T 4 T 2 Fig. 3.3. The edges of K.
(ii) To prove (3.23) we first use Green's formula

K div ψ vdx = - K ψ • ∇vdx + ∂K (ψ • n) v dγ. If v belongs to Q k+1 , ∇v belongs to RT [k] . Therefore, since ψ belongs to Ψ k , K div ψ v dx = ∂K (ψ • n) v dγ.
Let us decompose ∂K as

∂K = T 1 ∪ T 2 ∪ T 3 ∪ T 4 ,
as shown in Figure 3.3. Then, we check that for 1 ≤ j ≤ 4 and ∀q ∈ R k (T j ),

T j (ψ • n)qdγ = 0,
where R k (T j ) is the set of polynomials of degree k with respect to the abscissa along T j . Let us consider, for instance, j = 1 and assume that ψ is given by (3.24). Then

T 1 (ψ • n)qdγ = - 1 0 σ k (x)q(x) p y (0) = 0.
To conclude, it remains to notice that, if v belongs to Q k+1 , then v| T j belongs to R k (T j ).

This result suggests that we define X s h and X r h by

⎧ ⎨ ⎩ X s h = p h ∈ H(div; Ω) / ∀K ∈ T h , p h | K ∈ RT [k] , X r h = {p h ∈ H(div; Ω) / ∀K ∈ T h , p h | K ∈ Ψ k } . ( 3.25)
Let us now state the main result of this section, Theorem 3.7. The spaces (X h , M h ), defined in (2.16), satisfy assumptions (H0) to (H4). The spaces X s h and X r h that achieve the orthogonal decomposition (H1) of X h are the ones defined in (3.25).

Proof. (H0), (H3), and (H4) are classical properties of the Raviart-Thomas approximation spaces [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]. (H1) and (H2) are straightforward consequences (decomposing the integrals over Ω as the sum of the integrals over elements K) of the definition of Ψ k and Lemma 3.6. Remark 3.8. For k = 0 we obtain the orthogonal decomposition of X h described in section 3.3.1 that we summarize in Figure 3.4. In the same way we can summarize the orthogonal decomposition of X h for any k in Figure 3.5.

Q = 1 RT 0 Q x Q 1 P 0,1 1,1 P σ 0
= 1 Fig. 3.4. Orthogonal decomposition of Q 1 × Q 1 . k = Q k+1 1 2 k k+1 1
We can now apply Theorem 3.1 to the approximation problem (3.6), and by using the usual interpolation results (cf. [START_REF] Roberts | Mixed and hybrid methods[END_REF]) in X s h and M h defined by (2.16), (3.25) we obtain the following result.

Theorem 3.8. Problem (3.6) admits a unique solution

(p h , u h ) ∈ X h × M h , which satisfies • (p s h , u h ) → (p, u) in H(div ; Ω) × L 2 (Ω), • p r h → 0 in (L 2 (Ω)) 2 .
Furthermore, if we assume that the solution (p, u) of

(3.5) is such that (p, u) ∈ (H m (Ω)) 2 × (H m (Ω)), div p ∈ H m (Ω), and Ap ∈ (H m (Ω)) 2 for m = k + 1, then ⎧ ⎪ ⎨ ⎪ ⎩ p r h (L 2 ) 2 + p -p s h H(div;Ω) + u -u h L 2 ≤ C h m |u| H m + |p| (H m ) 2 + |div p| H m + |Ap| (H m ) 2 , (3.26)
where

|.| H m (resp., |.| (H m ) 2 ) denotes the usual seminorm in H m (Ω) (resp., (H m (Ω)) 2 ).
3.4. Application to the elliptic projection operator. We come back to the abstract framework described in section 3.2. We consider the problem of finding

( p h , u h ) ≡ Π h (p, u) such that ⎧ ⎨ ⎩ a(p -p h , q h )+ b(u -u h , q h ) = 0 ∀q h ∈ X h , b(w h , p -p h ) = 0 ∀w h ∈ M h . (3.27)
We set

D(B t ) = {w ∈ M / b(w, q) ≤ C(w) |q| H ∀q ∈ X} = w ∈ M / B t w ∈ H .
For w ∈ D(B t ), we have b(w, q) = (B t w, q) H ∀ q ∈ X.

Let us introduce the notation

⎧ ⎨ ⎩ |||p -p h ||| = p -p s h X + | p r h | H , |||(p, u) -Π h (p, u)||| C = |||p -p h ||| + u -u h M .
(3.28)

It is straightforward, by application of the abstract result, to get interpolation estimates. Theorem 3.9. Assume (H0)-(H4).

(i) For all (p, u) ∈ X × M , problem (3.27) admits a unique solution

Π h (p, u) = ( p s h + p r h , u h ) ∈ X h × (M h /Ker B t h ).
(ii) For all (p, u) ∈ X × D(B t ), there exists a constant C, independent of h, such that

|||(p, u) -Π h (p, u)||| C ≤ C E h (p, u),
where

E h (p, u) = inf q s h ∈X s h p -q s h X + inf w h ∈M h u -w h M + inf z s h ∈X s h |B t u -z s h | H .
In particular, |||pp h ||| and uu h M tend to zero.

Proof. (i) The existence and uniqueness still come from hypotheses (H0) and (H3) (with f = Bp which is, of course, in Im B).

(ii) The proof of this point is a minor modification of the proof of Theorem 3.2. The only difference concerns the treatment of the term b(u, p r h -q r h ) in (3.17). In order to get the error estimate, it is necessary to relate this term to the H-scalar product. In the proof of Theorem 3.2, this was done using the first continuous equation which implies that b(u, p r hq r h ) = -a(p, p r hq r h ) ≡ -(Ap, p r hq r h ) H . Now we assume that u ∈ D(B t ), so that we have

b(u, p r h -q r h ) = (B t u, p r h -q r h ) H .
The result is thus obtained by replacing -Ap with B t u.

4.

Error estimates for the evolution problem. Let us now come back to the initial evolution problem (2.5), (2.4) and see how we can relate the error estimates to the one obtained for the elliptic problem (3.5). This part of the analysis is inspired by [START_REF] Dupont | L 2 -estimates for Galerkin methods for second order hyperbolic equations[END_REF] developed for the approximation of second order hyperbolic problems with conforming finite elements. Although we have constructed this element in order to be able to do mass lumping, we analyze the error for the discrete problem without mass lumping. Of course, if mass lumping is applied, one should add to this error the quadrature error due to the numerical integration (see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Tordjman | Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes[END_REF][START_REF] Baker | The effect of quadrature errors on finite element approximations for the second-order hyperbolic equations[END_REF] and Remark 2.6).

4.1.

From error estimates for the elliptic problem to error estimates for the evolution problem: The abstract case. In this part, we use the same notation and hypotheses as in section 3.2, and we consider the evolution problem

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p, v) : [0, T ] → X × M such that d dt a(p, q) + b(v, q) =0 ∀q ∈ X, d dt (v, w) -b(w, p) =( f, w) ∀w ∈ M, p(0) = p 0 ; v(0) = v 0 . (4.1)
Or, equivalently, in an operator form, if we assume enough regularity on the solution in time,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p, v) ∈ (C 1 (0, T ; H) ∩ C 0 (0, T ; X)) × C 1 (0, T ; M ) such that A dp dt + B t v = 0 in X , dv dt -Bp = f in M , p(0) = p 0 ; v(0) = v 0 . (4.2)
In the following, we use the notation C m,r = C m (0, T ; H) ∩ C r (0, T ; X). X h and M h are, as in section 3.2, finite dimensional spaces satisfying hypotheses (H0) to (H4). We consider the approximate problem

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (p h , v h ) : [0, T ] → X h × M h such that d dt a(p h , q h ) + b(v h , q h ) =0 ∀q h ∈ X h , d dt (v h , w h ) -b(w h , p h ) =( f, w h ) ∀w h ∈ M h , p h (0) = p 0,h ; v h (0) = v 0,h . (4.3)
From the classical theory of ODEs, we have the following result.

(p h , v h ) ∈ C 1 (0, T ; X h ) × C 1 (0, T ; M h ).
Following [START_REF] Collino | Fictitious Domain Method for Unsteady Problems: Application to Electromagnetic Scattering[END_REF][START_REF] Dupont | L 2 -estimates for Galerkin methods for second order hyperbolic equations[END_REF] we introduce the elliptic operator defined in (3.27). By application of Theorem 3.9, we get the following interpolation results. Lemma 4.2. Let (p, v) be the solution of (4.2) and assume that (p, v) ∈ C 1,0 × C 1 (0, T ; M ). Then we have the following:

(i) There exists a primitive of v, u ∈ C 1 (0, T ; M ), satisfying

⎧ ⎪ ⎨ ⎪ ⎩ du dt = v, Ap 0 + B t u(0) = 0. (4.4)
This primitive is unique up to a constant element of Ker B t .

(ii) ∀t ∈ [0, T ], problem (3.27) admits a unique solution

Π h (p, u)(t) = ( p h , u h )(t) ∈ X h × (M h /Ker B t h
) and there exists a constant C independent of h such that

|||(p, u) -Π h (p, u)||| C (t) ≤ C E h (p, u)(t), (4.5)
where E h is defined by

E h (p, u)(t) = inf q s h ∈X s h p(t) -q s h X + inf w h ∈M h u(t) -w h M + inf z s h ∈X s h |Ap(t) -z s h | H .
In

∈ C k (0, T ; X) × C k (0, T ; M ), k ≥ 1, there exists a constant C independent of h such that ∂ k t p, ∂ k t u -Π h ∂ k t p, ∂ k t u (t) ≤ C E h ∂ k t p, ∂ k t u (t), (4.6
p ∈ C k while v ∈ C k-1
). This is a specificity of the new element and is due to the fact that in order to eliminate the terms involving an X-norm of elements in X r h , we used the first equation satisfied by the solution, which relates the embarrassing term b(v, p r hq r h ) to a term which is like an H-scalar product, or more precisely, for the evolution problem, which is the time derivative of something like an H-scalar product, say d dt a(p, p r hq r h ) (see proof of Theorem 3.2).

Proof (of Lemma 4.2). (i) We set f 0 = -Bp 0 ∈ Im B. From hypothesis (3.10), we know that there is a unique (p 0 , u 0

) ∈ X × (M/Ker B t ) such that ⎧ ⎨ ⎩ a(p 0 , q)+ b(u 0 , q) = 0 ∀q ∈ X, b(w, p 0 ) = -(f 0 , w) ∀w ∈ M,
which means that, p 0 being fixed, there is a unique u 0 ∈ (M/Ker B t ) such that Ap 0 + B t u 0 = 0. Now we define u as

u(t) = u 0 + t 0 v(s)ds.
It is clear that u ∈ C 1 (0, T ; M ) and is the unique solution of (4.4).

(ii) Let u ∈ C 1 (0, T ; M ) be the primitive of v; substituting this into the first equation of (4.2) gives

d dt (Ap + B t u) = 0 =⇒ (Ap + B t u)(t) = Ap 0 + B t u(0) = 0; thus, (p, u) ∈ C 1 (0, T ; X) × C 1 (0, T ; M ) satisfies a(p, q) + b(u, q) = 0 ∀q ∈ X,
so that u ∈ D(B t ) and B t u = -Ap. Applying Theorem 3.9, we get the existence and uniqueness of the elliptic projection, for t fixed, and also the error estimate (4.5).

(iii) If (p, u) is sufficiently regular in time, we can differentiate with respect to t and get

A∂ k t p + B t ∂ k t u = 0.
We again apply Theorem 3.9 to get (4.6).

We now give the main result.

(p 0,h , v 0,h ) = Π h (p 0 , v 0 ). (4.8) • If (p, v) ∈ C 2 (0, T ; X) × C 1 (0, T ; M ), then, ∀t ∈ [0, T ], |p -p s h | H (t) → 0 ; |p r h | H (t) → 0 ; v -v h M (t) → 0. More precisely, ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ |p -p s h | H (t) + |p r h | H (t) ≤ C E h (p, u)(t) + t 0 E h ∂ 2 t p, ∂ t v (s)ds , v -v h M (t) ≤ C E h (∂ t p, v) (t) + t 0 E h ∂ 2 t p, ∂ t v (s)ds . (4.9) • If, in addition, (p, v) ∈ C 3 (0, T ; X) × C 2 (0, T ; M ) and (p h , v h ) ∈ C 2 (0, T ; X h ) × C 2 (0, T ; M h ), ∀t ∈ [0, T ], p -p s h X (t) → 0. More precisely, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ p -p s h X (t) ≤ C(E h (p, u)(t) + E h ∂ 2 t p, ∂ t v (t) + t 0 E h ∂ 2 t p, ∂ t v (s) + E h ∂ 3 t p, ∂ 2 t v (s) ds). ( 4
= v -v h . (i) If (p, v) ∈ C 1,0 × C 1 (0, T ; M ), there exists a constant C 1 , independent of h, such that, ∀t ∈ [0, T ], | p s h -p s h | H (t) + | p r h -p r h | H (t) + v h -v h M (t) ≤ C 1 t 0 ∂ t δ h M (s) ds. (4.11) (ii) Moreover, if (p, v) ∈ C 2,1 ×C 2 (0, T ; M ) and (p h , v h ) ∈ C 2 (0, T ; X h )×C 2 (0, T ; M h ), there exists a constant C 2 , independent of h such that, ∀t ∈ [0, T ], p s h -p s h X (t) ≤ C 2 ∂ t δ h M (t) + t 0 ∂ t δ h M (s) + ∂ 2 t δ h M (s) ds . (4.12)
Proof (of Lemma 4.4). (i) To prove estimation (4.11), we begin by rewriting (4.1) with test functions q = q h ∈ X h ⊂ X and w = w h ∈ M h ⊂ M and by subtracting it from (4.3):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d dt a(p -p h , q h ) + b(v -v h , q h ) = 0 ∀q h ∈ X h , d dt (v -v h , w h ) -b(w h , p -p h ) = 0 ∀w h ∈ M h , (p -p h )(0) = p 0 -p 0,h ; (v -v h )(0) = v 0 -v 0,h . (4.13)
Introducing the elliptic projection Π h (p, u) = ( p h , u h ), we split the error between the approximate solution and the exact solution into two parts,

⎧ ⎨ ⎩ (p -p h )(t) = (p -p h )(t) + ( p h -p h )(t), (v -v h )(t) = (v -v h )(t) + ( v h -v h )(t), (4.14)
and we choose as approximate initial conditions the elliptic projection of the exact initial condition, (4.8), so that at time t = 0

( p h -p h )(0) = 0 ; v h -v h (0) = 0.
By (4.13) and (4.14), for any (q h , w By differentiating the first equation of (3.27) (written for (p, u)) with respect to t, we see that

h ) ∈ X h × M h , ⎧ ⎨ ⎩ a(∂ t ( p h -p h ), q h ) + b( v h -v h , q h ) = -a(∂ t (p -p h ), q h ) -b(v -v h , q h ), (∂ t ( v h -v h ), w h ) -b(w h , p h -p h ) = -(∂ t (v -v h ), w h ) + b(w h , p -p h ).
⎧ ⎨ ⎩ a(∂ t (p -p h ), q h ) + b(v -v h , q h ) = 0 ∀q h ∈ X h , b(w h , p -p h ) = 0 ∀w h ∈ M h .
Substituting into (4.15) gives, for any (q h , w h )

∈ X h × M h , ⎧ ⎨ ⎩ a(∂ t ( p h -p h ), q h ) + b( v h -v h , q h ) = 0, (∂ t ( v h -v h ), w h ) -b(w h , p h -p h ) = -(∂ t (v -v h ), w h ).
(4.16)

Furthermore, by taking q h = p hp h and w h = v hv h in (4.15) and by adding the two equations, we get

a(∂ t ( p h -p h ), p h -p h ) + (∂ t ( v h -v h ), v h -v h ) = -(∂ t (v -v h ), v h -v h ). (4.17) Next, set E h (t) = 1 2 (a( p h -p h , p h -p h ) + ( v h -v h , v h -v h )) (t).
Since for some constant C > 0, The second term is estimated using (4.11) of Lemma 4.4. This requires ∂ t δ h M to be estimated. For this, we use (4.6), for k = 2, which requires (p, v) ∈ C 2 (0, T ; X) × C 1 (0, T ; M ). We get

E 1/2 h (t) ≥ C | p h -p h | 2 H (t) + v h -v h 2 M (t)
∂ t (v -v h ) M (s) ≤ C E h ∂ 2
t p, ∂ t v (s); (4.24) then (4.22), (4.23), and (4.24) lead to the first inequality of (4.9). Now, for v, we write

v -v h M ≤ v -v h M + v h -v h M .
We apply (4.6) for k = 1 and get

v -v h M ≤ C E h (∂ t p, v) .
Again, using estimate (4.11) to bound v hv h M , we easily get the second inequality of (4.9).

• Estimate in X-norm. We write pp s h X ≤ pp s h X + p s hp s h X .

The first term is again bounded by |pp h | and can be estimated from (4.5). For the second term, we now use (4.12) of Lemma 4.4, which requires (p, v) ∈ C 2,1 × C 2 (0, T ; M ) and (p h , v h ) ∈ C 2 (0, T ; X h ) × C 2 (0, T ; M h ). In the right-hand side of (4.12) appears the second derivative of vv h . We thus use estimate (4.6) for k = 3, which requires (p, v) ∈ C 3 (0, T ; X) × C 2 (0, T ; M ), and we get (4.10).

Application to the approximation of the anisotropic wave equation with the new finite element, Q div

k+1 -Q k . We come back to the original problem described in section 2. We consider the approximate spaces given in (2.16). From section 3.3 we know that the Q div k+1 -Q k element is covered by the abstract framework. It is then straightforward to apply Theorem 4.3 and the following theorem results.

Theorem 4.5. Let (p, v) be the solution of (2.5), (2.4), and let u ∈ C 2 (0, T ; M ) be the primitive of v, defined by du dt = v ; B t u(0) = -Ap 0 .

Fig. 2 . 2 .

 22 Fig. 2.2. Support of the basis functions in RT[0] .

Fig. 2 . 5 . 2 -

 252 Fig. 2.5. Degrees of freedom in the Q div 2 -Q 1 element.

Fig. 2 . 6 . 3 -

 263 Fig. 2.6. Degrees of freedom in the Q div 3 -Q 2 element.

1 0τ

 1 y (y)p y (y)dy .

Fig. 3 . 5 .

 35 Fig. 3.5. Orthogonal decomposition of Q k+1 .

Theorem 4 . 1 .

 41 If f ∈ C 0 (0, T ; M h ), then problem (4.3) has a unique solution

  particular, |||pp h ||| and uu h M tend to zero uniformly in time (t ∈ [0, T ]) (|||.||| C and |||.||| are defined in (3.28)). (iii) In the same way, if (p, u)

  ) It is quite astonishing that the error estimate does not require the same regularity in time on v and on p, as one could expect (

	where ∂ k t g = Remark 4.1. Operators Π h and ∂ k d k g dt k . t commute, and we set
	(4.7)	v h	.
	Remark 4.2.

h = ∂ t ( u h ) = ∂ t u

  Theorem 4.3. Assume (H0)-(H4). Let (p, v) be the solution of (4.2) and (p h , v h ) the solution of the approximate problem (4.3) with the initial conditions

  ≤ C ∂ t (vv h ) M (t) ≡ C ∂ t δ h M (t),(4.18) where we set δ h = vv h . From the choice (4.8) for the initial conditions, we deduceE h (0) = 0.It is then easy to see that (4.18) leads to (4.11) (from the orthogonalityX s h ⊥ X r h , |q h | 2 H = |q s h | 2 H + |q r h | 2 H ∀q h ∈ X h ). (ii)To get the estimate (4.12) in the X-norm, we recall that, ∀ η s h ∈ X s h , we haveη s h = η 1 + η 2 with η 1 ∈ Ker B h and η 2 ∈ (Ker B h ) ⊥ , so that (recalling that η 1 X = |η 1 | H ∀η 1 ∈ Ker B h ). We set η h = p hp h . The term |η 1 | His already estimated from the inequality (4.11); therefore, in order to get the second inequality (4.12), we only need to estimate η 2 X . To do so, we start by recalling that the inf-sup condition (hypothesis (H2)) is equivalent toWe also know that η 2 X = η s h X/Ker B h . Thus, by taking q s h = η 2 in (4.[START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF]), Now use the second equation of (4.16) to obtainη 2 X ≤ C { ∂ t (vv h ) M + ∂ t ( v hv h ) M } . (4.20)Till now, we have used only the C 1 regularity of the solution. In order to bound ∂ t ( v hv h ) M , we need C 2 . Indeed, we want to apply (4.11) with v h replaced by ∂ t v h , v h by ∂ t v h , and so on. More precisely, we have∂ t ( v hv h ) M (t) ≤ CFinally, combining (4.20) and (4.21), we getη 2 X (t) ≤ C ∂ t δ h M (t) +and the proof is achieved. Proof (of Theorem 4.3). We combine results given in Lemma 4.2 and in Lemma 4.4.• Estimates in H-and M -norms. We have|pp s h | H + |p r h | H ≤ |pp s h | H + | p r h | H + | p s hp s h | H + | p r hp r h | H . (4.22) The first term on the right is bounded by |||pp h ||| and thus by |||(p, u) -Π h (p, u)||| C . Assuming that (p, v) ∈ C 1,0 × C 1 (0, T ; M ), by (4.5) it can be estimated that |||pp h ||| (t) ≤ C E h (p, u)(t). (4.23)

						t
						0	∂ 2 t δ h M (s) ds ,
						1/2 ,
	it follows that
			1/2 h dt dE	(t) η s h	2 X = |η 1 | 2 H + η 2	2 X
		⎧ ⎪ ⎪ ⎨	there exists a constant C > 0, independent of h, such that
	(4.19)	⎪ ⎪ ⎩	∀ q s h ∈ X s h , sup w h ∈M h	b(w h , q s w h M h )	≥ C q s h X/KerB h .
				sup w h ∈M h	b(w h , η 2 ) w h M	≥ C η 2 X .
						t
	(4.21)					0	∂ 2 t δ h M (s) ds.

This implies V h = {q h ∈ X h / div q h = 0} ⊂ Ker B, so that (3.8(ii)) is a consequence of (2.7(ii)).

Let (p h , u h ) be the solution of the approximate problem (2.9) with initial conditions (p 0,h , v 0,h ) = Π h (p 0 , v 0 ).

Proof. We apply Theorem 4.3 which relates the errors to quantities such as E h (∂ m t p, ∂ l t u), i.e., to the error due to the approximation of H(div; Ω) with the Raviart-Thomas RT [k] space and to the error due to the approximation of L 2 (Ω) with the Q k discontinuous elements.

Appendix. Gauss-Lobatto and Gauss-Legendre quadrature formulae in two dimensions. We give the locations of the quadrature points and the associated weights used respectively for the Q div 2 -Q 1 element (k = 1, Figure 2.5) and the

6).

For k = 1 (Figure 2.5):

• The Gauss-Lobatto points and weights:

;

9 .

• The Gauss-Legendre points and weights:

For k = 2 (Figure 2.6):

• The Gauss-Lobatto points and weights: • The Gauss-Legendre points and weights:

;


Keywords: integer programming, separation, cut generation, disjunctive cut, convex analysis, reverse polar programmation entière, séparation, génération de coupes, coupes disjonctives, analyse convexe, polaire inverse

The problem of separation is to find an affine hyperplane, or "cut", that lies between the origin O and a given closed convex set Q in a Euclidean space. We focus on cuts which are deep for the Euclidean distance, and facet-defining. The existence of a unique deepest cut is shown and cases when it is decomposable as a combination of facet-defining cuts are characterized using the reverse polar set. When Q is a split polyhedron, a new description of the reverse polar is given. A theoretical successive projections algorithm is proposed that could be used to compute deep facet-defining split cuts.

Introduction

The problem of separation is essential in combinatorial optimization, occuring for instance in the following context: a function must be optimized on a "complicated" set in R n , typically a set of integer points; this problem being too hard, an easier relaxation is solved; the problem is then to separate the solution of the relaxation from the complicated set, so as to tighten the relaxation.

Separating hyperplanes, called cuts in the community of combinatorial optimization, are used in many practical methods aimed at solving mixed integer programs, often embedded in a branch-and-bound framework like the branch-and-cut algorithm. An overview of the existing general techniques to compute such cuts can be found in [START_REF] Cornuéjols | Valid inequalities for mixed integer linear programs[END_REF]. We will focus on particular cuts, called split cuts which belong to the family of disjunctive cuts [START_REF] Balas | Disjunctive programming[END_REF]. In [START_REF] Balas | A lift-and-project cutting plane algorithm for mixed 0-1 programs[END_REF], lift-and-project cuts are defined and shown to be a particular case of split cuts for the mixed 0-1 case. In [START_REF] Perregaard | Generating cuts from multiple-term disjunctions[END_REF], a procedure is devised to generate one facet-defining lift-and-project cut. In this paper, we address the problem of generating facet-defining cuts which are deep; besides, we consider the general case of split cuts for the mixed integer case.

The key idea to address the general problem of separation between O and some convex set Q is to use a precise correspondence between facets of Q and extreme points of another convex set in the dual space: the reverse polar of Q. Moreover, a facet is deep in some sense when the corresponding extreme point of the reverse polar is close to the origin in the dual space. These two facts lead us to a theoretical algorithm aimed at computing deep facetdefining cuts by solving an optimization problem, more precisely a quadratic programming problem, on the reverse polar set. We then particularize the theory to the special case, of great practical importance, where the point to be separated is the solution of the linear relaxation of a mixed integer program, and the convex set is a disjunctive (split) polyhedron.

The article is organized as follows: section 1 introduces our main object, the reverse polar of a convex set. Section 2 advocates deep, facet-defining cuts as a good choice for the applications. A characterization of such cuts is given which uses the projection of O onto the reverse polar Q -and its decomposition as a convex combination of extreme points of Q -. Section 3 is devoted to the problem of computing the projection of a point onto a closed convex set, and the above-mentionned decomposition. Section 4 gives a characterization of Q -when Q is an explicitly described polyhedron, which shows that computing the projection of O onto Q -is tractable in this case. Section 5 introduces disjunctive programming and section 6 generalizes the method of section 4 to the case where Q is a split polyhedron. We finally propose a theoretical algorithm to compute deep, facet-defining cuts split cuts for general mixed integer programs.

RR n

¡ 6177 1 Reverse polar and normal cone of a convex set Let us first recall a few definitions and properties to introduce our key object: the reverse polar of a convex set. Let Q R n be a nonempty closed convex set, such that O / ∈ Q. We aim at separating O from Q. The elements of R n are identified with column vectors, and the standard Euclidean scalar product of vectors x and y is denoted by x • y. We will use the so-called support function σ S : R → R ∪ +∞ of a convex set S ⊂ R n :

σ S (d) := sup x∈S d • x. (1)
For this and other elementary concepts of convex analysis such as recession cones and exposed faces, we refer to [START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF].

Definition (Reverse polar) Let Q R n be a nonempty closed convex set with O / ∈ Q. The set Q -:= {d ∈ R n ; σ Q (d) -1} = {d ∈ R n ; ∀x ∈ Q : d • x -1} (2)
is called the reverse polar of Q.

This set was introduced by Balas [START_REF] Balas | Disjunctive programming[END_REF] and is very convenient for cutting: indeed, d separates O from Q if and only if d ∈ Q -(up to multiplication by a positive constant), via the hyperplane of equation d

• x = σ Q (d). Definition (Normal cone) Let Q R n be a nonempty closed convex set with O / ∈ Q. The set N Q := {d ∈ R n ; σ Q (d) 0} = {d ∈ R n ; ∀x ∈ Q : d • x 0} (3)
is called the normal cone of Q.

It is a (small) generalization of the usual notion of normal cone, since usually O ∈ Q is required in this definition. Note that this definition is very similar to the one of Q -: the -1 in the right hand side was simply replaced with 0. The normal cone will be useful in the description of Q -via the following theorem.

Theorem 1.1 (Recession cone) Let Q R n be a nonempty closed convex set with O / ∈ Q. The recession cone of Q -is exactly N Q : (Q -) ∞ = N Q . ( 4
)
Proof By definition of Q -and

N Q we have Q -⊂ N Q so that the recession cone of Q -is a subset of (N Q ) ∞ = N Q . Conversely, ∀d 1 ∈ Q -and ∀d 2 ∈ N Q , we have : ∀t 0, ∀x ∈ Q, (d 1 + td 2 ) • (x -q) -1 so d 1 + td 2 ∈ Q -, which implies that d 2 ∈ (Q -) ∞ .
The case when Q is (the closed convex hull of) the union of two convex sets will appear in the sequel, therefore we now recall the following result, taken from [START_REF] Cornuéjols | A convex-analysis perspective on disjunctive cuts[END_REF].

Theorem 1.2 (Reverse polar of a union) Let Q 0 and Q 1 be two nonempty closed con- vex sets in R n , with O / ∈ Q := conv(Q 0 ∪ Q 1 )
. There holds

Q -= Q - 0 ∩ Q - 1 . ( 5
) INRIA O Q Q - N Q Figure 1:
Normal cone and reverse polar of a polyhedron

"Good" cuts

There are infinitely many possible cuts, and this section is devoted to choosing good ones.

A first natural requirement for a cut is to touch Q. This section introduces two more requirements: to be deep and facet-defining.

Facet-defining cuts

The particular case where Q is a polyhedron deserves special attention. Recall that by definition, a facet of a polyhedron Q is a face (that is, the intersection of Q with some affine hyperplane supporting it) of maximal dimension, distinct from Q. When Q is a polyhedron, the inequalities that define a facet of Q play a particular role, as stated by the following theorem taken from chapter 8 in [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF].

Theorem 2.1 Let Q ⊂ R n be a full-dimensional polyhedron described by irredundant inequalities Ax b. Then Ax b is the unique minimal representation of Q, up to multiplication of the inequalities by positive scalars. Moreover, there exists a one-to-one correspondence between the facets of Q and the inequalities in Ax b.

As a consequence, when Q is a polyhedron, we are especially interested in cuts which expose a facet of Q. This motivates the next result (Theorem 6.2 of [START_REF] Cornuéjols | A convex-analysis perspective on disjunctive cuts[END_REF]).

Theorem 2.2 Let Q R n be a polyhedron such that lin(Q) = R n and O / ∈ Q. Then Q - is a polyhedron, and the face of Q exposed by d ∈ Q -is a facet of Q if and only if d * := - d σ Q (d) (6) is an extreme point of Q -. RR n ¡ 6177 F. Cadoux Q O Q - Figure 2: When lin(Q) R n . . .
Remark Theorem 2.2 does not hold without the assumption that lin(Q) = R n . For instance on Figure 2, in two dimensions, the linear hull lin(Q) is the x axis which is not full-dimensional. The orthogonal direction Q ⊥ is the y axis and Q -is invariant along the y axis, therefore having no extreme point. Theoretically, this problem can be avoided by working in the subspace lin(Q) instead of R n , but in practice we may not know anything about lin(Q). This problem will be dealt with in section 6.

Deep cuts

An idea supported in [START_REF] Cornuéjols | A convex-analysis perspective on disjunctive cuts[END_REF] is to consider additionally the depth of the cut, defined as the Euclidean distance from O to the cutting hyperplane. The following theorem states that there exists a deepest cut. It touches Q at the projection q * of O onto the convex Q.

Theorem 2.3 Let Q R n be a nonempty convex set with O / ∈ Q, q * the projection of O onto the convex set Q and H = {x ∈ R n : d • x = β} an hyperplane that separates O from Q. Let H be the hyperplane orthogonal to q * that contains q * . There holds

d(O, H) d(O, H ) and d(O, H) = d(O, H ) ⇐⇒ H = H . (7)
Proof Let h be the projection of O onto H and q the point on the line segment [O, q * ] which is contained in H (its existence is ensured by the Theorem of intermediate values, and uniqueness holds since O / ∈ H). There holds

d(O, H) = d(O, h) d(O, q ) d(O, q * ) = d(O, H )
and the inequalities are equalities if and only if h = q (by uniqueness of the projection on H) and q = q * (O, q and q * are aligned), that is to say, if and only if H = H . This result reveals the importance of the projection onto a closed convex set C (C = Q in Theorem 2.3). Recall the following characterization of a projection.

Deep facet-defining cuts

The two properties introduced above, to be deep and facet-defining, are antagonistic. For example, Figure 3 clearly shows that the deepest cut (in full line) does not expose a facet; and the two facet-defining cuts (in dashed line) are moderately deep. On the other hand, they imply the deepest cut. At this point, it is quite natural to look for such cuts, which expose facets of Q, and which imply the deepest cut as a consequence. This question can be more conveniently formulated in the dual space, thanks to the next result. Theorem 2.5 Let q * be the projection of O onto Q and d * its projection onto Q -. There holds

d * = - q * |q * | 2 . ( 8
)
Proof By definition, for all

d ∈ Q -, -1 σ Q (d) d • q * -d q * ;
so we have proved

d 1 q * , for all d ∈ Q -. ( 9
) Now set d * := -q * |q * | 2 : for all x ∈ Q, d * • x = - q * • (x -q * + q * ) q * 2 = - q * • (x -q * ) q * 2 -1; RR n ¡ 6177 F. Cadoux but q * • (x -q * ) 0 (Theorem 2.4
). So we have proved

d * • x -1 , for all x ∈ Q hence d * ∈ Q -. Since obviously |d * | = 1 |q * | , (9) establishes that d * has shortest norm in Q -.
This result is important with relation to our previous discussion: computing facet-defining cuts implying the deepest cut amounts to computing extreme points of Q -having the projection d * as a convex combination. Indeed, consider a set of hyperplanes defined by the equation

d k • x = σ Q (d k ); each of them separates 0 from Q if σ Q (d k ) < 0, and defines a facet of Q if d k is an extreme point of Q -(Theorem 2.2; then σ Q (d k ) = -1)
; and they altogether imply the deepest cut if q * =k α k d k , for some coefficients α k ≥ 0 (Theorem 2.3). The above result says that this latter property means that k α k d k is collinear to d * ; we may assume that the α k sum up to 1, so that

d * = k α k d k .
3 Projection algorithm

The previous section showed the importance of minimizing the quadratic function • 2 (the square of the Euclidean distance) in R n to compute deep cuts. It also showed the importance of decomposing the projection in extreme points. In this section, we present an algorithm to solve this problem where the extreme points are generated by a sequence of linear programs.

In section 6 we will see a situation where these linear programs are tractable.

The algorithm

Let C be a closed bounded convex set. The following algorithm computes the projection of the origin onto C by a column-generation mechanism, illustrated by Fig. 4. Note that, being bounded, C has an extreme point so that the initialization makes sense. Step 1 (Master problem). Compute the projection d * k of the origin onto the convex hull of d 1 , . . . , d k .

Step 2 (Oracle call). Compute an extreme point d k+1 ∈ C by minimizing over C the linear

function d → d * k • d. Step 3 (Stopping test). If d * k • d k+1 < d * k 2
, increase k by 1 and loop to Step 1. Otherwise stop.

Step 1 consists in solving for α

∈ R k min 1 2 k i=1 α i d i 2 , k i=1 α i = 1 , α i 0, i = 1, . . . , k . ( 10
) INRIA d1 • d = d1 2 0 d2 d 2 d 1 = d1 C Figure 4:
The first two iterations of algorithm 1 Any standard QP algorithm can be used to solve this standard QP problem; Ph. Wolfe described one in [START_REF] Ph | Finding the nearest point in a polytope[END_REF], adapted to the particular form (10). The output of any such algorithm is a set of positive α i 's; the corresponding d i 's are extreme points of C (by construction in Step 2). Algorithm 1 will become completely defined and implementable if, in addition to a QP solver, an oracle is available to minimize a linear function over C. This oracle must always return extreme points.

Remark In its original formulation [START_REF] Ph | Finding the nearest point in a polytope[END_REF], Wolfe's algorithm explicitly uses the list of points (here d 1 , . . . , d k ) onto wich the projection is computed. This list is scanned whenever the algorithm needs the smallest of k numbers of the form

X • d 1 , . . . , X • d k .
Instead of repeated calls to a QP solver, a variant of Algorihm 1 could call Wolfe's algorithm just once, to project directly the origin onto C. Instead of the (unknown) list of extreme points of C, Wolfe's algorithm would be input with the above-mentioned oracle, to minimize X • d over C.

If the algorithm stops at some iteration K, then Steps 2 and 3 imply that

d * K ∈ C satisfies d * K • d d * K
2 for all d ∈ C; in view of Theorem 2.4, this exactly means that d * K is the projection of the origin onto C. Now suppose Algorithm 1 were applied to C = Q -and did stop at some iteration K. In the language of the present paper, d * K would be the deepest cut (Theorem 2.5). More importantly, a number of cuts d i with α i > 0 would be produced, which would be facet-defining (Theorem 2.2); because d * K = i α i d i , they would imply the deepest cut as a consequence. In a word, Algorithm 1 would do just the job mentioned in section 2. However, the identification C = Q -violates the inherently unbounded character of Q -; subsection 3.3 will address this difficulty.

Convergence

The d k+1 produced by Step 2 of Algorithm 1 cannot be any of the previous d i 's; the algorithm is therefore finite if C has finitely many extreme points. However, a more intrinsic conver-

d(t) 2 d * k 2 -2t d * -d * k 2 + t 2 d k+1 -d * k 2 .
Proof In the development

d(t) 2 = d * k 2 + 2td * k • (d k+1 -d * k ) + t 2 d k+1 -d * k 2 ,
we work out the coefficient of t 0: from the definition of

d k+1 in Step 2 (d * ∈ C!): d * k • (d k+1 -d * k ) d * k • (d * -d * k ) = d * • (d * -d * k ) -d * k -d * 2 -d * k -d * 2 ,
where the last inequality comes from Theorem 2.4 (d * k ∈ C!).

d * k+1 2 d(t) 2 d * k 2 -2tδ + t 2 M for all t ∈ [0, δ/M ]. Setting in particular t = δ/M : d * k+1 2 d * k 2 - δ 2 M , for all k ∈ K δ .
The whole sequence { d * k 2 } is decreasing, so the above inequality can hold for finitely many k only: K δ is a finite set.

In other words: for δ arbitrarily small,

d * -d * k 2 δ if k is large enough; this means that d * k -d * 2 → 0.

The boundedness assumption

When C is unbounded, the proof of Theorem 3.2 suggests that algorithm 1 may not converge. Moreover, if C is unbounded, linear functions may be unbounded from below: Step 2 may produce no d k+1 . Indeed, the projection d * need not be a convex combination of extreme points: extreme rays may be necessary to describe it. In our present framework, this corresponds to a deepest cut impossible to describe by facet-defining cuts; see Fig. 5. In the unbounded case, Algorithm 1 may therefore be trying to compute non-existing objects.

INRIA q * Q - d * Q Figure 5:
The deepest cut need not be implied by facet cuts Admittedly, these arguments disappear when C is polyhedral: it can then be written

C = p i=1 α i d i + r j=1 λ j r j , (11)
where α varies in the unit-simplex and λ 0. The projection d * can thus be described as a combination of extreme points d i and extreme rays r j . With appropriate modifications, Algorithm 1 can still make sense in the unbounded but polyhedral case. We note, however, that extreme rays of C = Q -play a little role in our framework (Fig. 5 reveals that they do not define facet-cuts of Q). Besides, numerical difficulties still exist, as illustrated by Fig. 6: -At the first iteration, the oracle answers the horizontal ray r The above reasons explain that, in our application where C = Q -is definitely unbounded, we will artificially bound Q -by adding a normalization constraint.

RR n

¡ 6177

Case of a polyhedron defined by inequalities

We assume in this section that Q is a nonempty polyhedron explicitly described by inequalities. Then we characterize the reverse polar Q -as the envelope of its extreme points and extreme rays.

Characterization of the normal cone

The following lemma ([4], Theorem 6.7) characterizes the normal cone of Q.

N Q = {d = A u ; u 0 ; b • u 0}. ( 12
)
Said otherwise, N Q is the image by A of the cone:

K = {u 0 ; b • u 0}. ( 13
)
Remark Theorem 4.1 does not hold without the assumption that Q is nonempty. In practice however, we may not know whether Q is empty or not. This problem will be dealt with in section 6 and a counterexample to Theorem 4.1 will be shown when

Q = ∅. Let I -= {i ∈ 1...m : b i < 0} I 0 = {i ∈ 1...m : b i = 0} I + = {i ∈ 1...m : b i > 0}. (14)
The following result gives the extreme rays of K and generators of N Q .

Theorem 4.2 Let I -, I 0 , I + as in (14).

1. Let (e 1 , . . . , e m ) be the canonical basis of R m . The extreme rays of the cone K are exactly (up to multiplication by a positive scalar) the elements of

E = {e i ; i ∈ I -∪ I 0 } ∪ {b j e i -b i e j ; (i, j) ∈ I -× I + }. (15)
2. The image by A of E generates the normal cone N Q :

cone(A E) = N Q . ( 16
)
Proof Inspection of (13) shows that the extreme rays of K ⊂ R m are obtained as follows:

-Extract from {e 1 , . . . , e m , b} all possible sets of m -1 linearly independent vectors.

-Solve the corresponding systems of equations, thus producing a line in R m . Note that this line cannot be contained in K because K ⊂ R m + .

INRIA

-If this line has only 0 as feasible point in (13), the extreme point 0 ∈ K is produced; but no extreme ray. -Otherwise, the feasible half of that line is an extreme ray of K.

To extract an (m -1)-uple, i.e. to delete a couple from {e 1 , . . . , e m , b}, there are two possibilities.

1. Either one deletes e i (i ∈ [1, m]) and b, then the m -1 remaining vectors e k with k = i are automatically independent, and the intersection of K with the half-line R + e i is not the singleton {0} if and only if e i • b 0.

2. Or one deletes e i and e j (i = j). In this case the family is linearly independent if and only if (e i • b, e j • b) = (0, 0) (we can even assume that e i • b = 0 and e j • b = 0 otherwise we are back to the previous case where the line is directed by one of the basis vectors), and then the line defined by these m -1 vectors is

{u = αe i + βe j ; u • b = 0}
A directing vector of this line is (b • e j )e i -(b • e i )e j . The intersection of K with the line is nontrivial if and only if b • e i and b • e j have opposite sign (we assumed they were both nonzero).

Equation ( 16) comes from

N Q = A K = A cone(E) = cone(A E).
Remark Since the vectors in E contain only one or two nonzeros, the generators A E of N Q computed with this method are (the transpose of) either rows of A or linear combinations of two rows of A: if A is sparse, so will be the computed generators.

Remark Note that we considered the Cartesian product of I -and I + and E has a priori O(m 2 ) elements. More precisely, the number of computed generators is:

|E| = |I -| + |I 0 | + (|I -|)(|I + |). ( 17
)
This technique thus provides us with a set of vectors containing the extreme rays of N Q (q), but possibly also some "useless" vectors which are not extreme.

Example Consider for Q the square with vertices (-3,-1) , (-3,1) , (-1,-1) , and (-1,1). Let

A =     0 1 1 0 0 -1 -1 0     and b =     1 -1 1 3    
and since q = (0, 0) there holds w := (Aqb) = -b = (-1, 1, -1, -3). Thus I + = {1, 3, 4} and I -= {2}. The extreme rays of K are the vectors e 2 = (0, 1, 0, 0), a 2,1 = (1, 1, 0, 0), a 2,3 = (0, 1, 1, 0) and a 2,4 = (0, 3, 0, 1). Their image by A are respectively (1, 0), (1, 1), (1, -1) and (2, 0). Observe that the extreme rays of N Q are given by the second and the third vectors ((1, 1) and (1, -1)), whereas the two others are not extreme.

RR n

¡ 6177 F. Cadoux (1) (2) (3) (4) d 2 d 2,3 d 2,1 d 2,4

Characterization of the reverse polar

When the polyhedron Q is full dimensional and a description by inequalities of Q is available, it is easy to find the extreme points of Q -.

Theorem 4.3 Let Q R n be a full dimensional polyhedron explicitly described by a set of inequalities Ax b, with O / ∈ Q. Denote by a i the rows of A, so that

Ax b ⇐⇒ a i • x b i , i = 1, . . . , m.
1. If the description Ax b is minimal, then the extreme points of Q -are exactly the

P i := - a i b i with i ∈ I -. (18)
2. If the description Ax b is not minimal, the P i defined in (18) contains all the extreme points of Q -, plus some points which lie in Q -but are not extreme.

Proof When the description Ax b of Q is minimal, we know the facets of Q. A direct application of Theorem 2.2 yields the result. When it is not minimal, a minimal subset of inequalities in Ax b yields the extreme points of Q -.

Remark In practice, Q may not be full dimensional. We will see in section 6 how to get rid of this problem.

We are now able to state an important theorem which characterizes Q -when Q is a polyhedron defined by inequalities. Theorem 4.4 Let Q R n be a full dimensional polyhedron explicitly described by a set of inequalities Ax b, with O / ∈ Q. Let I -, I 0 , I + be defined by (14), the set E defined in Theorem 4.2 and the P i for i ∈ I -defined in Theorem 4.3. There holds:

Q -= conv(P i ; i ∈ I -) + cone(A E). (19)

INRIA

Proof Being a polyhedron, Q -is the sum of the convex hull of its extreme points and of its recession cone (Q -) ∞ . Theorem 4.3 gives the extreme points; Theorem 1.1 states that (Q -) ∞ = N Q and Theorem 4.2 that N Q = cone(A E).

Application to disjunctive programming

The results given so far apply to the problem of disjunctive programming, defined in [START_REF] Balas | Disjunctive programming[END_REF], with applications in integer programming. A first result is that the closed convex hull of the union of two (or more genrally, finitely many) polyhedra is a closed polyhedron. Indeed, it is the image by a linear transform (a projection) of a polyhedron in a higher dimensional space (lemma 6.4 of [START_REF] Cornuéjols | A convex-analysis perspective on disjunctive cuts[END_REF]).

Q Q 0 Q 1 O Figure 8: Disjunctive programming Definition (Disjunctive programming) Let Q 0 and Q 1 be two nonempty closed poly- hedra in R n , such that O / ∈ conv(Q 0 ∪ Q 1 ). The polyhedron Q = conv(Q 0 ∪ Q 1 ) is called a disjunctive polyhedron and a cut that separates O from Q is called a disjunctive cut.
The next subsection illustrates the idea of disjunctive programming. It is devoted to a particular and common type of disjunctions, the split disjunctions.

Split disjunctions

Consider a general mixed integer linear program of the form

       min c • x (x 1 , . . . , x n ) ∈ Z n (x n+1 , . . . , x n+p ) ∈ R p Ax b (20) RR n ¡ 6177 F. Cadoux with n, p, m ∈ N, A ∈ M m,n+p (R), c ∈ R n and b ∈ R m .
Practical algorithms used to solve such optimization problems often use the following relaxation of (20) :

   min c • x (x 1 , . . . , x n+p ) ∈ R n+p Ax b. (21)
which is much easier to solve. The problem of cut generation is to separate an optimal solution x of (21) from the feasible set in (20). Let R be the relaxed polyhedron:

R = {x ∈ R n+p ; Ax b} and let i ∈ 1 . . . n be such that x i / ∈ Z (if no such i exists, (20) is solved). Let π = e i and π 0 = x i . (22)
The feasible set of ( 20) is entirely contained in the union of the two following polyhedra

Q 0 = R ∩ {x ∈ R n ; π • x π 0 } Q 1 = R ∩ {x ∈ R n ; π • x π 0 + 1} (23) whereas x is not. Admitting that x is an extreme point of R, it cannot lie in conv(Q 0 ∪ Q 1 ).
To separate it from the feasible set in (20), it suffices to separate it from conv(Q 0 ∪ Q 1 ).

x Theorem 5.1 Let R R n be a full dimensional polyhedron explicitly described by m inequalities,

Q 0 Q 1
x ∈ R, (π, π 0 ) ∈ R n × R, Q 0 and Q 1 defined by (23) and Q = conv(Q 0 ∪ Q 1 ). Assume that x / ∈ Q and Q 0 = ∅, Q 1 = ∅. Then, for i = 0, 1, Q - i is a translated polyhedral cone: Q - i = P i + K i ,
where

P 0 = - π π 0 , P 1 = - π π 0 + 1 , (24)
and each K i is generated by m + 1 vectors a j i given by Theorem 4.4. As a result,

Q -= i=0,1 (P i + cone(a j i , j = 1 . . . m + 1)) = i=0,1 ({ m+1 j=1 α j i a j i , α j i 0}) (25)
Proof Note that x satisfies all the constraints defining each Q i except the splitting inequality: I -of ( 14) is a singleton, each Q - i has only one extreme point described by Theorem 4.3. The number of generators given by Theorem 4.2 is

|E| = |I -| + |I 0 | + |I -| |I + | = 1 + |I 0 | + |I + | = 1 + m.

Thus the description of

Q -= Q - 0 ∩ Q - 1 when Q 0 and Q 1 come
from a split cut is fairly simple, since I -is a singleton.

A particular cut

The description given in section 5.2 allows us to find a particular point of Q -, that is, a cut. Let {a i , i = 1 . . . n} be rows of A corresponding to n independent constraints active at x. Note that the a i are among the generators of both cones Q - 0 and Q - 1 computed in Theorem 5.1. Let P i be the only extreme point of Q - i (i = 1, 2). Consider the following linear system:

n i=1 α i a i = P 1 -P 0 = ( 1 π 0 - 1 π 0 + 1 )π (26)
which has a unique solution (α i ) i=1...n . Then d defined by

P 0 + i : αi 0 α i a i = P 1 - i : αi<0 α i a i =: d lies in Q - 0 ∩ Q - 1 ; therefore it is in Q -,
and it defines a cut. To compute it, it suffices to solve a linear system in dimension n. We will use this particular element of Q -in section 6.2.

Remark

This cut is the one used in [START_REF] Karamanov | Branching on general disjunctions (submitted)[END_REF] to estimate the quality of a branching direction.

Algorithm

The previous sections naturally suggest to compute split cuts by projecting the origin onto Q -, decomposing the projection as a convex combination of extreme points, and using these extreme points of Q -as cutting directions. But before using this approach, we need to deal with several issues. Theorem 2.2 fails when lin(Q) = R n .

The characterization of the normal cone in Theorem 4.2 fails in general for an empty polyhedron.

Algorithm 1 can only be used on a bounded polyhedron (a polytope); otherwise, the linear program in step 2 may be unbounded.

According to (25), Q -is described via variables α ∈ R 2m+2 . If the LP in step 2 of algorithm 1 has several solutions, an extremal such solution in this extended space need not correspond to an extremal solution in R n .

When lin

(Q) = R n
When lin(Q) = R n , Theorem 2.2 fails: facets of Q no longer correspond to extreme points of Q -. Theoretically, this can be avoided by working in the Euclidean space lin(Q) instead of R n . In the applications, we will ignore this problem and act as if lin(Q) was equal to the whole space1 R n . For this reason, Q -may have no extreme point at all. However, thanks to the normalization constraint presented in the next subsection, we will never run into this case: Q -will be replaced by a bounded set Q-.

Normalization constraint

It is clear from its definition that Q -is unbounded. Unfortunately, the black box used by the successive projection algorithm needs the poyhedron to be bounded, so we must somehow bound Q -artificially to use our method. The description (25) suggests to bound Q -by replacing it with Q-, defined as

Q-=    d = m+1 j=1 α j 0 a j 0 = m+1 j=1 α j 1 a j 1 , α j i 0, i,j α j i N    ( 27
)
where N is a sufficiently large normalization constant. Note that it suffices to compute the particular cut described in section 5.3 to get a value of N (say, N = i |α i | where α is the solution of the linear system (26

)) ensuring that Q-= ∅. O a 1 a 2 a 3 Q A B F Q - Q- Figure 10: A degenerated example where lin(Q) = R n
Remark Replacing Q -with Qmay introduce new extreme points which do not correspond to facets of Q. The corresponding cuts will therefore not be facet-defining. This seem to be unavoidable, since we have seen in Figure 5 that the deepest cut need not be implied by facet-defining cuts.

Figure 10 illustrates the problem raised by lin(Q) = R n and the role of the normalization constraint. Q is defined by the constraints a i • x b i (i = 1, . . . , 3) and lin(Q) is the horizontal axis. The projection of O onto Qis decomposed as a convex combination of points A and B, and two cuts are generated which expose the same facet F of Q.

Normal cone of an empty polyhedron

In this subsection, we show that Lemma 4.1 does not hold whithout the assumption that the polyhedron is nonempty. Let n = 2 and

A = 1 0 -1 0 and b = -1 -1 . Then Q := {(x = (x 1 , x 2 ) ∈ R 2 ; Ax b} = ∅. The set of u ∈ R 2 + such that b • u 0 is simply R 2
+ and the set A u for such u is R × {0}. If Lemma 4.1 held, the normal cone of Q would be R × {0}, but since Q is empty, its normal cone is the whole space R 2 .

In our situation, Q = conv(Q 0 ∪ Q 1 ) may be assumed nonempty (otherwise the whole problem makes little sense). However, either Q 0 or Q 1 could well be empty, an event which might be computationally expensive to check. Fortunately, this difficulty is minor. In fact, assume for example

Q 0 = ∅, so that Q -= Q - 1 .
No matter which wrong set (call it Q-0 ) is obtained via Lemma 4.1 and Theorem 4.9 for the reverse polar of Q 0 , elements of Q-

0 ∩ Q - 1
will lie in Q - 1 = Q -: they will provide valid cuts for Q. In other words, in our application to integer programming, there is no risk to cut a feasible solution.

The black box may return a non-extreme point

In algorithm 1, the oracle called in step 2 was supposed to return an extreme point of Q -. However, in our application to disjunctive programming, this cannot be guaranteed. The situation of figure 11, where Q is the intersection of two cones as in (25), may happen: the whole line segment AB is optimal for the minimization of x → c • x over Q -, and the points A, B and C (note that C is not extreme in Q -) can all be returned by the solver which is provided the description (25) of Q -. This problem is not too serious however: it occurs in degenerate cases only (c has to be very particular) and if it does, the possible existence of a non-extreme point of Q -in the decomposition of the projection of O onto Q -corresponds to generating a cut which is not facet-defining, but at least it is valid.

Algorithm

Summing up, we suggest the following method to separate the origin from the closed convex hull of the union of two disjunctive polyhedra.

-Consider a mixed integer program (20) and obtain an extreme point, optimal solution x to the linear relaxation (21). Translate the origin O to x. -Choose a splitting direction i such that x i is not integer although x i is subject to integrality in (20). Split the polyhedron of the linear constraints in (20) into Q 0 and Q 1 . -Use Theorem 5.1 to build Q 0 -and Q 1 -.

-Intersect Q 0 -and Q 1 -to get Q -.

-Normalize Q -into Q-, for example with the help of section 5.

INRIA

This cut generation algorithm was implemented and tested on small random instances of integer programs. The results are encouraging, but more numerical experiments are still needed to evaluate its practical interest. In particular, it would be interesting to embed such a cut generator in an existing branch-and-cut algorithm to see how it performs on classical problem libraries. This extensive numerical work is beyond the scope of this paper and is a natural direction for future work.

Conclusion

We aim at finding "good" hyperplanes separating O from a convex set Q. The reverse polar Q -of Q determines cutting directions, and moreover, when Q is a polyhedron, its extreme points (up to intersection with a convenient subspace) determine cuts that expose a facet of Q. Besides, the problem of projecting O onto Q to find the deepest cut amounts to projecting it onto Q -. As a consequence, we suggest to compute the projection of O onto Q -along with its decomposition in extreme points of Q -. We then use these points to generate several cuts that separate O from Q. This technique applies in particular to disjunctive programming, when the polyhedra are explicitly described by inequalities : it is then possible to do all the computations explicitly. This results in an implementable algorithm, the practical interest of which still needs to be studied.

INRIA Theorem 2 . 4 QFigure 3 :

 243 Figure 3: Deepest cut and facets of Q

Algorithm 1 (

 1 Projection by column generation) Choose an extreme point d 1 of C; set k = 1.

  gence proof points out the role of boundedness of C and lays the ground for a complexity study.Lemma 3.1 Setting d(t) := d * k + t(d k+1d * k ), there holds for all t 0

Theorem 3 . 2 2 M

 322 The sequence {d * k } converges to d * . Proof Because d k+1 and d * k vary in the bounded set C, d k+1d * k for some number M . Now take δ > 0 and call K δ the set of k such that d *d * k 2 > δ. Diminishing δ if necessary, we may assume δ/M 1; then d(δ/M ) lies in the segment [d * k , d k+1 ] and, from Step 1:

Figure 6 :

 6 Figure 6: Numerical instablity in the unbounded case

Lemma 4 . 1

 41 Let Q = {x ∈ R n ; Ax b} where m ∈ N, A ∈ M m,n (R) and b ∈ R m . Assume Q nonempty and O / ∈ Q.Then the normal cone of the polyhedron Q is

Figure 7 :

 7 Figure 7: Generators of a normal cone

Figure 9 : 5 . 2

 952 Figure 9: Particular case of a split cut

Figure 11 :

 11 Figure 11: The point returned by the black box may not be extreme

3 ,

 3 and design a black box which minimizes over Q-. -Use the successive projection algorithm 1 to compute the projection of the origin O onto Q -and its decomposition in extreme points d 1 , . . . , d k of Q -. -Use the constraints d i • x -1, i = 1, . . . , k as cuts.

Or lin(Q) = {Aex = be} if the problem has equality constraints Aex = be.INRIA

Acknowledgement. The author wishes to thank his advisor Claude Lemaréchal for his help, comments and suggestions all along this research.