E Creuse Lamav

M Farhloul
email: farhlom@umoncton.ca

L Paquet
email: luc.paquet@univ-valenciennes.fr

A POSTERIORI ERROR ESTIMATION FOR THE DUAL MIXED FINITE ELEMENT METHOD FOR THE p-LAPLACIAN IN A POLYGONAL DOMAIN

Keywords: p-Laplacian, dual mixed FEM, a posteriori error estimators, Helmholtz decomposition AMS Subject Classification: 65N15, 65N30, 65N50, 35J60

For the discrete solution of the dual mixed formulation for the p-Laplace equation, we define two residues R and r. Then we bound the norm of the errors on the two unknowns in terms of the norms of these two residues. Afterwards, we bound the norms of these two residues by functions of two error estimators whose expressions involve at the very most the datum and the computed quantities. We next explain how the discretized dual mixed formulation is hybridized and solved. We close our paper by numerical tests for p = 1.8 and p = 3 firstly to corroborate the orders of convergence established by M. Farhloul and H. Manouzi [1], and secondly to experimentally verify the reliability of our a posteriori error estimates.

Introduction

Let us fix a bounded plane domain Ω with a polygonal boundary. More precisely, we assume that Ω is a simply connected domain and that its boundary Γ is the union of a finite number of linear segments Γ j , 1 ≤ j ≤ n e (Γ j is assumed to be an open segment). In this domain we consider the p-Laplace equation

   -div (|∇ u| p-2 ∇ u) = f in Ω, u = 0 on Γ, (0.1)
where 1 < p < ∞, ∇ u is the vector gradient of u,

|∇ u| = 2 i=1 ∂u ∂x i 2 1 2
, and f is a given function in L q (Ω) (q conjugate of p). This equation occurs in many mathematical models of physical processes such as nonlinear diffusion and filtration, power-law materials, and viscoelastic materials.

In Farhloul and Manouzi [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF], these authors have introduced a mixed finite element method for problem (0.1) based on the introduction of σ = |∇ u| p-2 ∇ u as a new unknown and its approximation by the lowest degree Raviart-Thomas finite element (cf. Thomas [START_REF] Thomas | Sur l'analyse numérique des méthodes d'éléments finis hybrides et mixtes[END_REF], Raviart and Thomas [START_REF] Raviart | A mixed finite element method for second order elliptic problems, Mathematical aspects of finite element methods[END_REF]). The old unknown u is approximated by a piecewise constant function on the triangulation of Ω. They have proved a priori error estimates without any restriction on p. In this paper, we are concerned by the construction of reliable a posteriori error estimators for σ h the approximation of σ and u h the approximation of u.

In section 1, we firstly introduce two residues R and r and we bound the norms of the errors ||σ -σ h || 0,q,Ω and ||u-u h || 0,p,Ω in terms of these two residues (|| . || 0,p,Ω denotes the L p -norm in Ω and 1 q = 1 -1 p ). In this part of the work, we distinguish two cases : p ≥ 2 and 1 < p < 2. Secondly, using a Helmholtz type decomposition for arbitrary vector field τ ∈ [L q (Ω)] 2 , we establish a bound of the norm of R in terms of a local error estimator

η 1 = K∈T h η 1 (K) p 1 p
, where T h denotes the triangulation of Ω and K an arbitrary triangle of T h . η 1 (K) is function of A (σ h ) = |σ h | q-2 σ h , curl A (σ h ), the tangential jump of A (σ h ) at each interface of two triangles A (σ h ).t and the diameter h K of the triangle K.

The norm of r is easily bounded in terms of the local error estimator

η 2 = K∈T h η 2 (K) q 1 q
, where η 2 (K) = ||f -P 0 h f || 0,q,K . P 0 h f is the function constant on each K ∈ T h equal on K to the mean of f on K.

In section 2 of our paper, we introduce the hybrid formulation of the discrete mixed formulation by expressing that the normal component of σ h is continuous at each interface between triangles of T h . We obtain in that way a nonlinear system in the unknowns σ h , u h and λ h where λ h is in some sense an approximation of u along the edges of the triangulation T h . We next want to solve this nonlinear system by a simple Picard scheme by replacing the non linear term

Ω A (σ h ). τ h dx by Ω |σ m h | q-2 σ m+1 h
. τ h dx. Introducing the "generalized barycenter" of K:

x m K = K |σ m h (x )| q-2 x dx K |σ m h (x )| q-2 dx ∈ K, (0.2)
and writing the lowest degree Raviart-Thomas vectorfield on

K σ m+1 K = σ m+1 h |K , K ∈ T h , in the form σ m+1 K (x ) = a m+1 K + c m+1 K (x -x m K ), ∀x ∈ K, (0.3
σ m+1 K (x ) = ∂K λ m+1 h n K ds K |σ m K | q-2 dx - 1 2|K| K
f dx (x -x m K ) ∀x ∈ K, (0.4)

u m+1 K = 1 4|K| 2 K f dx K |σ m K (x )| q-2 |x -x m K | 2 dx (0.5) + 1 2|K| ∂K λ m+1 h (x -x m K ). n K ds,
n K denoting the unitary normal vectorfield along the boundary of K. Expressing the continuity of the normal component of the vectorfield σ m+1 h at each interface (i.e. at each interior edge e) of the triangulation T h , and applying the Picard scheme, we obtain at step m + 1 a linear system in the unknowns λ m+1 e = λ m+1 h |e , e running over the set of all interior edges of the triangulation T h of the polygonal domain Ω. Once λ m+1 h has been calculated, we then use it in formulas (0.4) and (0.5) in order to obtain σ m+1 h and u m+1 h , and continue the iterations. Finally, in section 3, we present a numerical test in the square Ω =] -1/2; 1/2[ 2 . Denoting by r(x ) = r(x 1 , x 2 ) the distance to the origin, we consider the bounded function u defined in the square by

u(x ) =        1 4 -r 2 (x ) 2 e -r 2 (x )/ε if r(x ) ≤ 1 2 , 0 if r(x ) > 1 2 . (0.6) f = -div (|∇ u| p-2 ∇ u)
is easily computed and shown to be equal to:

f (x ) = -2 p-1 g p-2 (x ) 2 + 1 ε g(x ) p-2 e -(p-1) r 2 (x ) ε r p-2 (x ) × 2(p -1) 2 + 1 ε g(x ) r 2 (x ) + 2(p -1)g(x ) 3 + 1 ε g(x ) r 2 (x ) ε -pg(x ) 2 + 1 ε g(x ) , (0.7)
where

g(x ) = 1 4 -r 2 (x ). (0.8) Ones verifies that f ∈ L q (Ω) if the conjugate p > 1+ √ 5 2
= 1.618... We test this example for p = 1.8 and p = 3. For each of these two cases, four graphics are given: ||σ -σ h || 0,q,Ω , ||u -u h || 0,p,Ω , the ratio between ||σ -σ h || 0,q,Ω and the estimated error on σ h and finally the ratio between ||u -u h || 0,p,Ω and the estimated error on u h , all with respect to the number of degrees of freedom. One verifies the predicted order of convergence given respectively by theorem 3.2 for p = 1.8 and theorem 3.1 for p = 3 of Farhloul and Manouzi [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] and that the ratios between the errors and the estimated errors are asymptotically independent of the number of degrees of freedom.

Definition of residues and bounds of the errors

Let us first recall the weak formulation of the mixed method for the p-Laplacian introduced by Farhloul and Manouzi in [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]. Let M denotes the Banach space L p (Ω) (1 < p < +∞) and X the following:

X = τ ∈ [L q (Ω)] 2 ; div τ ∈ L q (Ω) , (1.9)
where q denotes the harmonic conjugate of p. Endowed with the norm

||τ || X = ||τ || q 0,q,Ω + ||div τ || q 0,q,Ω 1 q , X is also a Banach space. For every τ ∈ [L q (Ω)] 2 , let us set A (τ ) = |τ | q-2 τ . It is straightforward to check that when τ ∈ [L q (Ω)] 2 , A (τ ) ∈ [L p (Ω)] 2 .
The weak formulation of the mixed method related to the Dirichlet problem for the p-Laplacian (0.1) reads as follows [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]: Find (σ , u) ∈ X × M such that the following equations are satisfied:

         Ω A (σ ). τ dx + Ω u div τ dx = 0 ∀τ ∈ X (i) Ω (div σ ) v dx + Ω f v dx = 0 ∀v ∈ M (ii) (1.10)
It is proved in [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] (theorem 2.1) that (1.10) possesses one and only one solution (σ , u) and that u is the weak solution of (0.1) and σ = |∇ u| p-2 ∇ u. Now let us consider a regular family of triangulations (T h ) h>0 on Ω [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. The following finite element spaces are defined in [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]:

X h = τ h ∈ X; τ h|K ∈ RT 0 (K), ∀K ∈ T h , (1.11)
M h = v h ∈ M ; v h|K ∈ P 0 (K), ∀K ∈ T h , (1.12)
where P 0 (K) denotes the space of constant functions on the triangle K and [START_REF] Thomas | Sur l'analyse numérique des méthodes d'éléments finis hybrides et mixtes[END_REF][START_REF] Raviart | A mixed finite element method for second order elliptic problems, Mathematical aspects of finite element methods[END_REF]:

RT 0 (K) = P 0 (K) 2 ⊕ P 0 (K) x , x = (x 1 , x 2 ). (1.13)
Then the finite element approximation to problem (1.10) is given by [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]:

Find (σ h , u h ) ∈ X h × M h satisfying          Ω A (σ h ). τ h dx + Ω u h div τ h dx = 0, ∀τ h ∈ X h Ω (div σ h ) v h dx + Ω f v h dx = 0, ∀v h ∈ M h (1.14)
The existence and uniqueness of the solution of the discrete problem (1.14) follows from the uniform inf-sup condition [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]. A priori error estimates are also proved in Farhloul and Manouzi [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]. Now, we define two residues R on X and r on M by:

< R, τ > = (A (σ h ), τ ) + (div τ , u h ), ∀τ ∈ X, (1.15) < r, v > = (div σ h , v) + (f, v), ∀v ∈ M. (1.16)
R is a continuous linear form on X and r is a continuous linear form on M . When we will speak of ||R|| (respectively ||r||), it will mean the norm of R (respectively r) as a continuous linear form on X (respectively M ). We have the following bound of (A (σ h ) -A (σ ), σ h -σ ) in terms of the norms of these two residues:

Proposition 1.1. (A (σ h ) -A (σ ), σ h -σ ) ||R|| ||σ h -σ || 0,q,Ω + ||r|| sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X + 2||r|| ||R|| (1.17)
Before proceeding to the proof, let us make the following remark: Proof. Using the definition of R (1.15) and the equation (1.10 (i)), we obtain:

(A (σ h ) -A (σ ), τ ) + (div τ , u h -u) =< R, τ >, ∀τ ∈ X.
(1.18)

Substracting (1.10 (ii)) form (1.16), we obtain:

(div (σ h ) -div (σ ), v) =< r, v >, ∀v ∈ M. (1.19) Substracting (1.19) with v = u h -u from (1.18) with τ = σ h -σ , we obtain: (A (σ h ) -A (σ ), σ h -σ ) =< R, σ h -σ > -< r, u h -u > . (1.20)
By the inf-sup condition stated in proposition 2.1 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] and (1.18) it follows:

||u h -u|| M sup τ ∈X (div τ , u h -u) ||τ || X sup τ ∈X < R, τ > ||τ || X + sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X = ||R|| + sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X . (1.21)
From (1.20) and (1.21) it follows that:

(A (σ h ) -A (σ ), σ h -σ ) ||R|| ||σ h -σ || X + ||r|| ||R|| + ||r|| sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X (1.22) But: ||σ h -σ || X ∼ ||σ h -σ || 0,q,Ω + ||div (σ h -σ )|| 0,q,Ω = ||σ h -σ || 0,q,Ω + sup v∈L p (Ω) < div (σ h -σ ), v > ||v|| 0,p,Ω = ||σ h -σ || 0,q,Ω + ||r|| (1.23)
by (1.19). The result follows now from (1.22) and (1.23).

Our purpose now is to bound the norms of the errors ||σ h -σ || 0,q,Ω and ||u h -u|| 0,p,Ω in terms of the norms of the residues ||R|| and ||r||. We will have to distinct two cases: p ≥ 2 and 1 < p < 2. We begin with the case p ≥ 2.

Firstly, we need the following inequality due to Sandri [START_REF] Sandri | Sur l'approximation numérique des écoulements quasi-Newtoniens dont la viscosité suit la loi puissance ou la loi de Carreau[END_REF] (see also [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] p. 74):

Lemma 1.1. We suppose p ≥ 2. Then:

(A (σ h ) -A (σ ), σ h -σ ) ||σ h -σ || 2 0,q,Ω ||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω + Ω |A (σ h ) -A (σ )| |σ h -σ | dx. (1.24)
We will also need the following inequality:

Lemma 1.2. Supposing p ≥ 2, we have ∀τ ∈ X:

(A (σ ) -A (σ h ), τ ) Ω |A (σ h ) -A (σ )| |σ h -σ | dx 1 p ||τ || 0,q,Ω , (1.25)
the multiplicative constant hidden in the sign being independent of τ .

The proof of inequality (1.25) relies on the following technical lemma:

Lemma 1.3. ∀p ≥ 2, ∀x , y ∈ R 2 , |x | q-2 x -|y | q-2 y p-1 |x -y |.
(1.26)

Proof.

( (2) One verifies easily that the left-hand side of inequality (1.26) is a homogeneous function of degree 1, i.e. that

) 1
||tx | q-2 tx -|ty | q-2 ty | p-1 = t ||x | q-2 x -|y | q-2 y |
for every t > 0 which implies using the first step of our proof that inequality (1.26) is still true under the weaker condition that |x | = |y |.

(

||y | q-2 y -|λx | q-2 λx | p-1 |y -λx |. (1.27) But |y -λx | ≤ |y -x | + |x ||1 -λ| = |y -x | + |x | |1 - |y | |x | | ≤ 2 |y -x |.
(1.28)

From (1.27) and (1.28), it follows that

||y | q-2 y -|λx | q-2 λx | p-1 |y -x |. (1.29)
On the other hand:

||λx | q-2 λx -|x | q-2 x | p-1 = |y | |x | q-1 -1 |x | q-2 x p-1 = ||y | q-1 -|x | q-1 | p-1 . (1.30)
But it is elementary to prove that for s ≥ 1 and α ≥ 0, β ≥ 0:

|α 1/s -β 1/s | s ≤ |α -β|. ( 1
||y | q-1 -|x | q-1 | p-1 ≤ ||y | -|x || ≤ |y -x |. (1.32)
From (1.30) and (1.32), it follows that:

||λx | q-2 λx -|x | q-2 x | p-1 ≤ |x -y |. (1.33)
From (1.29) and (1.33) we obtain:

||y | q-2 y -|x | q-2 x | p-1 ≤ 2 p-1 ||y | q-2 y -|λx | q-2 λx | p-1 + ||λx | q-2 λx -|x | q-2 x | p-1 |x -y |,
what was to be proved.

Now we give the proof of lemma 1.2:

Proof.

(A (σ ) -A (σ h ), τ ) Ω |A (σ ) -A (σ h )| p dx 1/p ||τ || 0,q,Ω = Ω ||σ | q-2 σ -|σ h | q-2 σ h | p dx 1/p ||τ || 0,q,Ω = Ω ||σ | q-2 σ -|σ h | q-2 σ h | p-1 |A (σ ) -A (σ h )| dx 1/p ||τ || 0,q,Ω Ω |σ -σ h ||A (σ ) -A (σ h )| dx 1/p ||τ || 0,q,Ω (1.34)
due to lemma 1.3. This proves the result.

Using lemma 1.1, 1.2 and inequality (1.22), we are now in a position to prove the following result: Proposition 1.2. Supposing p ≥ 2, we have the following bounds in terms of the norms of the residues:

||σ h -σ || 2 0,q,Ω ||R|| 2 + ||r|| q + ||r|| ||R|| (1.35)
and

Ω |A (σ h ) -A (σ )| |σ h -σ | dx ||R|| 2 + ||r|| q + ||r|| ||R||. (1.36)
Proof. Reading inequality (1.24) from right to left and using inequality (1.17), we obtain:

||σ h -σ || 2 0,q,Ω ||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω + Ω |A (σ h ) -A (σ )| |σ h -σ | dx ||R|| ||σ h -σ || 0,q,Ω + ||r|| sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X + 2||r|| ||R|| ||R|| ||σ h -σ || 0,q,Ω + ||r|| Ω |A (σ h ) -A (σ )| |σ h -σ | dx 1/p + 2||r|| ||R|| (1.37)
by lemma 1.2. Now using Young's inequality, i.e. ∀a ≥ 0, ∀b ≥ 0, ab ≤ 1 p a p + 1 q b q , we obtain ∀ ε > 0 and ∀ ε > 0 from the previous inequality:

||σ h -σ || 2 0,q,Ω ||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω + Ω |A (σ h ) -A (σ )| |σ h -σ | dx ε -1 (||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω )||R|| 2 + ε ||σ h -σ || 2 0,q,Ω ||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω +(ε) -q ||r|| q + (ε) p Ω |A (σ h ) -A (σ )| |σ h -σ | dx + ||r|| ||R||. (1.38)
Choosing adequately ε > 0 and ε > 0, it follows from the previous inequality the two following ones:

||σ h -σ || 2 0,q,Ω (||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω ) × (||σ h || 2-q 0,q,Ω + ||σ || 2-q 0,q,Ω ) ||R|| 2 + ||r|| q + ||r|| ||R|| (1.39) and Ω |A (σ h )-A (σ )| |σ h -σ | dx (||σ h || 2-q 0,q,Ω +||σ || 2-q 0,q,Ω )||R|| 2 +||r|| q +||r|| ||R||.(1.40)
On the other hand, following the proof of theorem 2.1 of [START_REF] Farhloul | A mixed finite element method for a nonlinear Dirichlet problem[END_REF] (see also theorem 2.1 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]), we have

||σ || X ≤ ||f || 0,q,Ω and ||σ h || X ≤ ||f || 0,q,Ω . (1.41)
Using this fact in the right-hand side of (1.39) (respectively (1.40)), we obtain the inequality (1.35) (respectively (1.36)).

From (1.21), lemma 1.2 and (1.36), it follows that:

Corollary 1.1. Supposing p ≥ 2, we have the following bound for the norm of the error u -u h in terms of the norms of the residues:

||u h -u|| 0,p,Ω ||R|| + (||R|| 2 + ||r|| q + ||r|| ||R||) 1/p . (1.42)
We now turn to the case 1 < p < 2. Firstly, we are going to bound the norm of the error on σ h , i.e. ||σ h -σ || 0,q,Ω .

Proposition 1.3. Supposing 1 < p < 2, we have the following bound on the norm of the error on σ h :

||σ h -σ || 0,q,Ω ||R|| p/q + ||r|| 2/q + (||r|| ||R||) 1/q .
(1.43)

Proof. By inequality (3.14) of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] (see also Sandri [START_REF] Sandri | Sur l'approximation numérique des écoulements quasi-Newtoniens dont la viscosité suit la loi puissance ou la loi de Carreau[END_REF]), we have:

||σ h -σ || q 0,q,Ω + Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx (A (σ h ) -A (σ ), σ h -σ ).(1.44)
Using inequality (1.17) from proposition 1.1, it follows from the previous inequality

||σ h -σ || q 0,q,Ω + Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx ||R|| ||σ h -σ || 0,q,Ω + ||r|| sup τ ∈X (A (σ ) -A (σ h ), τ ) ||τ || X + 2||r|| ||R||. (1.45)
By inequality (3.15) of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF],

||A (σ )-A (σ h )|| 0,p,Ω Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx 1/2 [||σ || 0,q,Ω + ||σ h || 0,q,Ω ] q-2 2
, and Young's inequality, we obtain from inequality (1.45) that, ∀ε > 0 and ∀ε > 0,

||σ h -σ || q 0,q,Ω + Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx ε q ||σ h -σ || q 0,q,Ω + 1 ε p ||R|| p + ε Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx + 1 ε [||σ || 0,q,Ω + ||σ h || 0,q,Ω ] q-2 ||r|| 2 + ||r|| ||R||. (1.46)
Choosing adequately ε > 0 and ε > 0, il follows from the previous inequality the following one:

||σ h -σ || q 0,q,Ω ||R|| p + [||σ || 0,q,Ω + ||σ h || 0,q,Ω ] q-2 ||r|| 2 + ||r|| ||R||. (1.47) Using (1.41), it follows that ||σ h -σ || q 0,q,Ω ||R|| p + ||r|| 2 + ||r|| ||R||,
and taking the two-sides of this inequality to the power 1 q , we obtain the result.

Proposition 1.4. Supposing 1 < p < 2, we have the following bound on the norm of the error on u h :

||u -u h || 0,p,Ω ||R|| + ||r|| + ||R|| p/2 . (1.48)
Proof. By the definition of R (see (1.15)) and equation (1.10 (i)), we have

(A (σ h ) -A (σ ), τ ) + (div τ , u h -u) =< R, τ >, ∀τ ∈ X. (1.49)
By the inf-sup inequality (cf. proposition 2.1 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]) and the previous equation, we have the following bound on ||u h -u|| 0,p,Ω :

||u h -u|| 0,p,Ω sup τ ∈X < R, τ > ||τ || X + sup τ ∈X (A (σ h ) -A (σ ), τ ) ||τ || X ||R|| + ||A (σ h ) -A (σ )|| 0,p,Ω .
(1.50)

By inequality (3.15) of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] and (1.41), it follows that

||A (σ h ) -A (σ )|| 0,p,Ω Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx 1/2
.

(1.51)

On the other hand, inequality (1.46) implies

Ω (|σ | + |σ h |) q-2 |σ h -σ | 2 dx ||R|| p + ||r|| 2 + ||r|| ||R||. (1.52)
From (1.50), (1.51) and (1.52), it follows that

||u -u h || 0,p,Ω ||R|| + (||R|| p + ||r|| 2 + ||r|| ||R||) 1/2 ||R|| + ||R|| p/2 + ||r|| + ||r|| 1/2 ||R|| 1/2 ||r|| + ||R|| + ||R|| p/2 .
Now, in both cases p ≥ 2 and 1 < p < 2, we are reduced to bound the norms of the residues R and r in terms of computed or known quantities. This is immediate for r.

Proposition 1.5. ||r|| = η 2 where η 2 = K∈T h η 2 (K) q 1/q and η 2 (K) = ||f -P 0 h f || 0,q,K , ∀K ∈ T h .
Proof. r has been defined by equation (1.16) as being div σ h +f . But by the second equation of the discretized mixed formulation (1.14), we have div σ h = -P 0 h f , where P 0 h denotes the projection operator defined by:

P 0 h : L q (Ω) -→ M h v -→ P 0 h v, (1.53)
P 0 h v being the function constant on each triangle K equal to the mean of v on K. Thus r = f -P 0 h f and

||r|| = sup v∈M \{0} < r, v > ||v|| M = ||f -P 0 h f || L q (Ω) = K∈T h η 2 (K) q 1/q .
We turn now, to the far more serious problem to bound the norm of R which has been defined by equation (1.15). To do that we need the Helmholtz decomposition of arbitrary vectorfields τ ∈ [L q (Ω)] 2 with estimates.

Theorem 1.1. Let q ∈]1, +∞[. For an arbitrary vectorfield τ ∈ [L q (Ω)] 2 such that div τ ∈ L q (Ω), there exists a function z ∈ W 2,q (Ω) and a function Φ ∈ W 1,q (Ω) such that:

τ = ∇ z + curl Φ (1.54) ||z|| W 2,q (Ω) + ||Φ|| W 1,q (Ω) ||τ || 0,q,Ω + ||div τ || 0,q,Ω . (1.55)
Proof. Let D be a bounded domain of the plane, with a regular boundary and containing Ω. Let us denote by (div τ ) ∼ the extension of div τ to D by 0 out of Ω.

Of course, (div τ ) ∼ ∈ L q (D). In D, we consider the boundary value problem:

   ∆S = (div τ ) ∼ in D S |∂D = 0. (1.56)
This problem possesses one and only one solution in H 1 0 (D). Moreover by theorem 2.4.2.5 of [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF] p 124-125, the variational solution S ∈ W 2,q (D). It follows from the closed graph theorem that also:

||S|| W 2,q (D) ||(div τ ) ∼ || 0,q,D = ||div τ || 0,q,Ω . (1.57)
Let us set z = S |Ω . Consequently, we have:

||z|| W 2,q (
Ω) ||div τ || 0,q,Ω .

(1.58)

By (1.56), we have that div (∇ z) = div τ . Let us set:

v = τ -∇ z. (1.59) v is a divergence-free vectorfield in Ω, v ∈ [L q (Ω)] 2 and
||v || 0,q,Ω ||τ || 0,q,Ω + ||div τ || 0,q,Ω (1.60) by inequality (1.58). We now suppose that q ≥ 2 (the proof in the case 1 < q < 2 will be done later). This hypothesis implies that v ∈ [L 2 (Ω)] 2 . This allows us to apply theorem I.3.1 p 37 of [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and algorithms[END_REF] which tells us that there exists one and only one function Φ ∈ H 1 (Ω) ∩ L 2 0 (Ω) such that the divergence-free vectorfield v may be written in the form v = curl Φ.

(1.61)

Obviously, |Φ| W 1,q (Ω) = ||v || 0,q,Ω as v 1 = ∂Φ ∂x2 and v 2 = -∂Φ ∂x1 .
Moreover by Sobolev's imbedding theorem and Poincaré's inequality for functions in H 1 (Ω) of mean zero [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF]:

||Φ|| 0,q,Ω ||Φ|| H 1 (Ω) |Φ| H 1 (Ω) |Φ| W 1,q (Ω) .
Putting these two facts together, it follows that:

||Φ|| W 1,q (Ω) ||v || 0,q,Ω .
(1.62)

Using inequalities (1.62) and (1.60), we obtain that ||Φ|| W 1,q (Ω) ||τ || 0,q,Ω + ||div τ || 0,q,Ω .

(1.63)

Adding the two inequalities (1.58) and (1.63), we obtain (1.55). Also from the equalities (1.59) and (1.61) follows (1.54).

We now turn to the case 1 < q < 2. The proof of this case requires the following density lemma.

Lemma 1.4. Let us suppose that q ∈]1, 2[. The vectorial subspace {v ∈ [L 2 (Ω)] 2 ; div v = 0} is dense in the vectorial space {v ∈ X; div v = 0} for the norm of X.
Proof. By lemma 1.1.2 p. 16 of Mghazli's thesis [START_REF] Mghazli | Une méthode mixte pour les problèmes d'hydrodynamique[END_REF], we know that [D(Ω)] 2 is dense in X. Let v ∈ X such that div v = 0 be given. There exists a sequence of vectorfields (v * k ) k≥1 ⊂ [D(Ω)] 2 such that (v * k ) k≥1 tends to v in the sense of the norm of X. Proceeding like in (1.56)-(1.58) but with the vectorfield v * k instead of τ in the right-hand side of the first equation of (1.56), we construct a sequence (z k ) k≥1 ⊂ W 2,q (Ω) (theorem 2.4.2.5 of [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF] p 124-125 that we have used is also true in the case 1 < q < 2) such that

||z k || W 2,q (Ω) ||div v * k || 0,q,Ω (1.64)
and div (∇

z k ) = div v * k , ∀k ≥ 1. Let us set v k = v * k -∇ z k .
v k is a square integrable vectorfield and also divergence-free. Thus v k belongs to the vectorial subspace {v ∈ [L 2 (Ω)] 2 ; div v = 0}. It remains to prove that the sequence (v k ) k≥1 tends to v in the sense of the norm of X.

||v k -v || X ≤ ||v * k -v || X + ||∇ z k || X .
(1.65)

v * k tends to v in the sense of the norm of X. On the other hand by inequality (1.64) :

||∇ z k || X ||div v * k || 0,q,Ω (1.66)
and as div v * k tends to div v = 0 in L q (Ω) as k → +∞, it follows that ||∇ z k || X → 0 as k → +∞. From these two facts and inequality (1.65) it follows that ||v k -v || X → 0 as k → +∞.

We now turn back to the proof of theorem 1.1 when 1 < q < 2.

(Proof of theorem 1.1 when 1 < q < 2.) By the preceding lemma, there exists a sequence (v k ) k≥1 of square-integrable divergence-free vectorfields on Ω converging to the vectorfield v defined by (1.59) in the sense of the norm of X. Applying theorem I.3.1 p. 37 of [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and algorithms[END_REF] to the vectorfield v k , there exists one and only one function

φ k ∈ H 1 (Ω) ∩ L 2 0 (Ω) such that curl φ k = v k . As v k → v in [L q (Ω)] 2 , the sequence (∇ φ k ) k≥1 is a Cauchy sequence in [L q (Ω)] 2 .
Moreover, by a compacity argument similar to that used to prove lemma 1.8 in [START_REF] Sossa | Refined mixed finite element method for the Poisson problem on a polygonal domain with a reentrant corner[END_REF], it is easily seen that the mapping:

W 1,q (Ω) → R + : Ψ → ∂Ψ ∂x 1 0,q,Ω + ∂Ψ ∂x 2 0,q,Ω + Ω Ψ dx (1.67)
is an equivalent norm on W 1,q (Ω). In view of (1.67) and the fact that the functions φ k , k ≥ 1, are all of mean 0, it follows that (φ k ) k≥1 is a Cauchy sequence in W 1,q (Ω). Let us set Φ = lim k→+∞ φ k in W 1,q (Ω). Φ is also of mean 0 and curl Φ = v . By these facts and (1.67), it follows:

||Φ|| W 1,q (Ω) |Φ| W 1,q (Ω) = ||v || 0,q,Ω . (1
X: < R, τ > = K∈T h (A (σ h ), ∇ z -π h ∇ z) + K∈T h (curl (A (σ h )), Φ -I cl (Φ)) - E∈E h < A (σ h ). t E , Φ -I cl (Φ) > E (1.69)
where:

• (z, Φ) ∈ W 2,q (Ω) × W 1,q (Ω) denotes the Helmholtz decomposition of the vectorfield τ ∈ X, • I cl (Φ) is the Clément interpolate of Φ (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and algorithms[END_REF] p. 109-111),

• σ h is defined by (1.14)

, • E h denotes the set of all edges of the triangulation T h , • A (σ h )). t E denotes the tangential jump of A (σ h ),

• π h ∇ z is the Raviart-Thomas interpolate of lowest degree of ∇ z.

Proof. By (1.15) for every τ ∈ X, < R, τ >= (A (σ h ), τ ) + (div τ , u h ).

(1.70) By theorem 1.1, for every τ ∈ X, there exists (z, Φ) ∈ W 2,q (Ω) × W 1,q (Ω) such that τ = ∇ z + curl Φ and the estimates (1.55) holds. Using this Helmholtz decomposition, (1.70) may be rewritten:

< R, τ >= (A (σ h ), ∇ z) + (A (σ h ), curl Φ) + (div (∇ z), u h ). (1.71)
Now, ∇ z ∈ [W 1,q (Ω)] 2 and thus its Raviart-Thomas interpolate of the lowest degree is well defined (see (14) of [START_REF] Sossa | Refined mixed finite element method for the Poisson problem on a polygonal domain with a reentrant corner[END_REF]). Let us denote it π h (∇ z). We have:

(div π h (∇ z), v h ) = (div (∇ z), v h ), ∀v h ∈ M h .
Thus (1.71) may be rewritten:

< R, τ >= (A (σ h ), ∇ z) + (A (σ h ), curl Φ) + (div (π h (∇ z)), u h ). (1.72)
Now by the first equation of the discrete formulation (1.14):

(A (σ h ), τ h ) + (u h , div τ h ) = 0, ∀τ h ∈ X h . (1.73)
Applying (1.73) with τ h = π h (∇ z) and τ h = curl (I cl (Φ)) successively, we obtain the following equations:

(div (π h (∇ z)), u h ) = -(A (σ h ), π h (∇ z)) (1.74) and (A (σ h ), curl (I cl (Φ))) = 0. (1.75)
Injecting the relations (1.74) and (1.75) in the right-hand side of (1.72), we obtain :

< R, τ > = (A (σ h ), ∇ z -π h (∇ z)) + (A (σ h ), curl (Φ -I cl (Φ))) = (A (σ h ), ∇ z -π h (∇ z)) + K∈T h (A (σ h ), curl (Φ -I cl (Φ))) = (A (σ h ), ∇ z -π h (∇ z)) + K∈T h {(curl A (σ h ), Φ -I cl (Φ))-< A (σ h ).t , Φ -I cl (Φ) > ∂K }
(by Green's formula by adapting (1.31) p. 19 of [START_REF] Mghazli | Une méthode mixte pour les problèmes d'hydrodynamique[END_REF]) We now define the error estimator η 1 :

= K∈T h (A (σ h ), ∇ z -π h (∇ z)) + K∈T h (curl A (σ h ), Φ -I cl (Φ)) - E∈E h < A (σ h ). t E , Φ -I cl (Φ) > E . ( 1
Definition 1.1. For each K ∈ T h , we define η 1 (K) by:

η 1 (K) p = h p K ||A (σ h )|| p 0,p,K + h p K ||curl (A (σ h ))|| p 0,p,K + E⊂∂K h E || A (σ h ).t E || p 0,p,E
and we define η 1 by

η 1 = K∈T h η 1 (K) p 1/p .
We are now in a position to bound ||R|| by a constant times the error estimator η 1 .

Theorem 1.2.

||R|| η 1 .

(1.77)

Proof. It follows from lemma 1.5 that for every τ ∈ X:

| < R, τ > | K∈T h ||A (σ h )|| 0,p,K ||∇ z -π h ∇ z|| 0,q,K + K∈T h ||curl A (σ h )|| 0,p,K ||Φ -I cl (Φ)|| 0,q,K + E∈E h || A (σ h ). t E || 0,p,E ||Φ -I cl (Φ)|| 0,q,E .
(1.78) By (3.6) p. 72 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] (see also lemma 3.1 p. 125 of [START_REF] Farhloul | A mixed finite element method for a nonlinear Dirichlet problem[END_REF]):

||∇ z -π h ∇ z|| 0,q,K h K |∇ z| 1,q,K .
(1.79) By lemma 3.1 p. 57 of [START_REF] Verfürth | A Review of a posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]:

||Φ -I cl (Φ)|| 0,q,K h K |Φ| 1,q,ωK (1.80)
and

||Φ -I cl (Φ)|| 0,q,E h 1/p E |Φ| 1,q,ωE . (1.81)
Let us recall that ω K denotes the union of K with all the triangles from the triangulation T h adjacent to the triangle K and that ω E denotes the union of the at most two triangles of T h admitting E as an edge. Using inequalities (1.79)-(1.81) to bound the right-hand side of inequality (1.78), we obtain:

| < R, τ > | K∈T h h K ||A (σ h )|| 0,p,K |∇ z| 1,q,K + K∈T h h K ||curl A (σ h )|| 0,p,K |Φ| 1,q,ωK + E∈E h h 1/p E || A (σ h ).t E || 0,p,E |Φ| 1,q,ωE K∈T h h p K ||A (σ h )|| p 0,p,K 1/p |∇ z| 1,q,Ω + K∈T h h p K ||curl A (σ h )|| p 0,p,K 1/p |Φ| 1,q,Ω + E∈E h h E || A (σ h ). t E || p 0,p,E 1/p |Φ| 1,q,Ω
and so:

| < R, τ > | K∈T h h p K ||A (σ h )|| p 0,p,K + h p K ||curl A (σ h )|| p 0,p,K + E⊂∂K h E || A (σ h ). t E || p 0,p,E 1/p |∇ z| q 1,q,Ω + |Φ| q 1,q,Ω 1/q .
Using (1.55) and definition 1.1, we obtain:

| < R, τ > | K∈T h η 1 (K) p 1/p (||τ || 0,q,Ω + ||div τ || 0,q,Ω ) η 1 ||τ || X (1.82)
This proves (1.77).

Remark 1.3. As already said in the introduction:

||r|| = ||f -P 0 h f || 0,q,Ω = K∈T h η 2 (K) q 1/q , where η 2 (K) = ||f -P 0 h f || 0,q,K , ∀K ∈ T h . (1.83)
Thus ||r|| is immediately bounded in terms of the data. We may now turn to the hybridization of the discretized mixed formulation (1.14).

Hybridization of the discrete mixed formulation

The construction of a basis of vectorfields in X h presents the difficulty that for a vectorfield τ h to be in X h , its normal component must be continuous at each interface between two adjacent triangles of the triangulation T h . To overcome this difficulty, we introduce the enlarged space

Xh = τ h ∈ [L q (Ω)] 2 ; τ h|K ∈ RT 0 (K), ∀K ∈ T h (2.84)
and also the space of Lagrange multipliers

Mh = µ h : ∂T h → R; µ h | •
e is constant for every edge e of T h and µ h | • e = 0 if the edge e ⊂ ∂Ω .

(2.85)

In (2.85), ∂T h denotes the "skeleton" of the triangulation T h , i.e. the union of all the edges of the triangulation T h and for an edge e of T h ,

• e denotes its "interior", that is to say the edge e without its extremities. We are now in a position to introduce the "hybrid formulation" of (1.14): find ( σ h , ũh , λ h ) ∈ Xh × M h × Mh such that:

                           Ω A ( σ h ). τ h dx + K∈T h K ũh div τ h dx - K∈T h ∂K λ h τ h . n K ds = 0 ∀ τ h ∈ Xh , K∈T h K (div σ h )v h dx + Ω f v h dx = 0 ∀v h ∈ M h , K∈T h ∂K µ h σ h . n K ds = 0 ∀µ h ∈ Mh .
(2.86)

In (2.86), K denotes the unit exterior normal to the triangle K. Arguing as in [START_REF] Sossa | Refined mixed finite element method for the Poisson problem on a polygonal domain with a reentrant corner[END_REF] p. 527, it is easily seen that if ( σ h , ũh , λ h ) ∈ Xh × M h × Mh is a solution of the hybrid formulation (2.86), then necessarily σ h = σ h and ũh = u h . Now let us prove that (2.86) possesses at most one solution. Thus, let ( σ h , ũh , λ h ) ∈ Xh ×M h × Mh be a solution of (2.86). We know already that this implies that σ h = σ h and ũh = u h . Let us consider for τ h ∈ Xh the vectorfield τ K,e ∈ Xh such that for the edge e of the triangle K ∈ T h , τ K,e . n K equals 1 on e and 0 on the two other sides of the triangle K, and τ K,e is null on all other triangles K ′ = K of the triangulation T h . Putting it into the first equation of (2.86), we obtain:

|e|λ h|e = K A (σ h ). τ K,e dx + K u h div τ K,e dx.
(2.87) Equation (2.87) shows that the Lagrange multiplier λ h is also completely determined. Note that the existence of the above vectorfield τ h is proved in lemma 2.2 of [START_REF] Sossa | Refined mixed finite element method for the Poisson problem on a polygonal domain with a reentrant corner[END_REF]. Thinking a bit to equation (2.87) shows that it proves also the existence of a solution to the nonlinear system (2.86). Taking σ h = σ h and ũh = u h , we can define λ h by equation (2.87) noting that these special vectorfields τ K,e form a basis of the space Xh . Let us observe that if e has two adjacent triangles K 1 and K 2 , we may choose any in the right-hand side of equation (2.87) as it gives the same result. This follows from equation (1.14)(i) with τ h equals to τ K1,e on K 1 , -τ K2,e on K 2 and 0 on the other triangles of T h . Note also that if the edge e is contained in the boundary of Ω, then τ K,e ∈ X h and thus by the first equation of (1.14), the right hand side of (2.87) is zero showing that λ h|e = 0 in this case. We have thus established the following result:

Proposition 2.1. The hybrid formulation (2.86) possesses one and only one solution

( σ h , ũh , λ h ) ∈ Xh × M h × Mh . Moreover, σ h = σ h and ũh = u h .
Numerically, we may try to solve the nonlinear problem (2.86) by a simple Picard scheme also called the method of successive substitutions. Thus we construct a sequence (σ m , u m , λ m ) ∈ Xh ×M h × Mh by the following iterative algorithm: Given

(σ m , u m , λ m ) ∈ Xh × M h × Mh , compute (σ m+1 , u m+1 , λ m+1 ) ∈ Xh × M h × Mh solution of the linear system                            Ω |σ m | q-2 σ m+1 . τ h dx + K∈T h K u m+1 div τ h dx - K∈T h ∂K λ m+1 τ h . n K ds = 0 ∀ τ h ∈ Xh , K∈T h K (div σ m+1 )v h dx + Ω f v h dx = 0 ∀v h ∈ M h , K∈T h ∂K µ h σ m+1 . n K ds = 0 ∀µ h ∈ Mh .
(2.88)

It is not evident that the linear system (2.88) possesses one and only one solution

(σ m+1 , u m+1 , λ m+1 ) ∈ Xh × M h × Mh . Let us observe that if (σ m+1 , u m+1 , λ m+1
) is a solution to (2.88), then a fortiori (σ m+1 , u m+1 ) ∈ X h × M h due to the third equation of (2.88) and is solution of the following linear system:

         Ω |σ m | q-2 σ m+1 . τ h dx + Ω u m+1 div τ h dx = 0 ∀τ h ∈ X h , Ω div σ m+1 v h dx + Ω f v h dx = 0 ∀v h ∈ M h . (2.89)
We first give a sufficient condition for the linear system (2.89) to have one and only one solution.

Proposition 2.2. If for every triangle K of the triangulation T h , σ m |K ≡ 0, then the linear system (2.89) possesses one and only one solution (σ m+1 , u m+1 ) ∈ X h × M h .

Proof. By corollary I.4.1 of [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and algorithms[END_REF] it suffices to prove that the bilinear form a(., .) :

X h × X h → R : (τ ′ h , τ ′′ h ) → Ω |σ m | q-2 τ ′ h . τ ′′ h dx is coercive on V h = {τ h ∈ X h ; Ω div τ h v h dx = 0, ∀v h ∈ M h }. One must in fact check also that the bilinear form b(., .) : X h × M h → R : (τ h , v h ) → Ω div τ h v h dx
satisfies the inf-sup condition but this results from proposition 3.1 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF]. Let us consider τ h ∈ V h . As τ h|K ∈ RT 0 (K), ∀K ∈ T h , it follows immediately from the condition τ h ∈ V h that τ h must be constant on each triangle K of the triangulation T h . Thus for such a vectorfield, one has:

a(τ h , τ h ) = Ω |σ m | q-2 |τ h | 2 dx = K∈T h |τ h|K | 2 K |σ m | q-2 dx > 0.
Due to our hypothesis on σ m , it follows that:

c = inf K∈T h 1 |K| K |σ m | q-2 dx > 0.
Thus:

a(τ h , τ h ) ≥ c   K∈T h |K| |τ h|K | 2 1/2   2 ∼ ||τ h || 2 ,
all norms being equivalent on finite dimensional vector spaces. This proves the coercivity of the bilinear form a(., .) on V h . Now, it remains to prove that under the same condition on σ m , that the linear system (2.88) possesses one and only one solution (σ m+1 , u m+1 , λ m+1 ) ∈ Xh ×M h × Mh . This is the purpose of the following corollary.

Corollary 2.1. If for every triangle K of the triangulation T h , σ m |K ≡ 0, then the linear system (2.88) possesses one and only one solution (σ m+1 , u m+1 , λ m+1 ) ∈ Xh × M h × Mh .

Proof. We know already that the linear system (2.89) possesses one and only one solution (σ m+1 , u m+1 ) ∈ X h × M h . Now, we put that solution in system (2.88) and define λ m+1 by the first equation of (2.88). More precisely, we define λ m+1 |e by:

|e|λ m+1 |e = K |σ m | q-2 σ m+1 . τ K,e dx + K u m+1 div τ K,e dx, (2.90)
where K is an arbitrary adjacent triangle to the edge e and τ K,e is the same vectorfield that appears in formula (2.87). Let us note that if the edge e is contained in the boundary of Ω, then τ K,e ∈ X h and thus the first equation of (2.89) implies that λ m+1 |e defined by (2.90) must be 0. Thus λ m+1 ∈ Mh which proves that the linear system (2.88) possesses one and only one solution.

Definition 2.1. In the sequel, we will denote for a triangle K ∈ T h by x m K the barycenter of the probability measure of density

|σ m |K | q-2 K |σ m |K | q-2 dx with respect to
Lebesgue measure on R 2 restricted to K that is ( [START_REF] Dieudonné | Eléments d'Analyse[END_REF] p. 146):

x m K = 1 K |σ m |K | q-2 dx K x |σ m (x )| q-2 dx.
It is well defined for every K ∈ T h under the hypothesis of proposition 2.2 on σ m . Proposition 2.3. Under the hypothesis of proposition 2.2 on σ m , for every triangle K of the triangulation T h we have:

(σ m+1 |K )(x ) = ∂K λ m+1 n K ds K |σ m | q-2 dx -K f dx 2|K| (x -x m K ) ∀x ∈ K, (2.91)
where |K| denotes the area of the triangle K.

Proof. Let K be an arbitrary triangle of the triangulation T h . Putting in the first equation of (2.88), any vector field τ h constant on K and nul outside K, we obtain :

K |σ m |K | q-2 σ m+1 |K dx = ∂K λ m+1 n K ds. (2.92)
Writing the unknown vectorfield σ m+1 |K in the form:

σ m+1 |K (x ) = a + c(x -x m K ) (2.93)
where a ∈ R 2 , c ∈ R and x m K has been defined in definition (2.1), it follows from equation (2.92):

a K |σ m |K | q-2 dx = ∂K λ m+1 n K ds.
(2.94)

Thus we have:

a = ∂K λ m+1 n K ds K |σ m |K | q-2 dx . ( 2
u m+1 |K = K |σ m |K | q-2 |x -x m K | 2 dx K f dx 4|K| 2 + 1 2|K| ∂K λ m+1 (x -x m K ). n K ds.
(2.98)

Proof. Let us take for the vectorfield τ h in the first equation of (2.88), the vectorfield equals to x -x m K in the triangle K and to 0 outside of K. We obtain using (2.91) and the fact that K |σ m (x )| q-2 (x -x m K )dx = 0:

- K |σ m | q-2 (x )|x -x m K | 2 dx . 1 2|K| K f dx + 2|K|u m+1 |K = ∂K λ m+1 (x -x m K ). n K ds.
From this last equation, (2.98) follows easily.

Remark 2.1. Equations (2.91) and (2.98) show that σ m+1 |K and u m+1 |K are completely defined by the restriction of the Lagrange multiplier λ m+1 to the boundary of K. Thus, we are reduced to the determination of the Lagrange multiplier λ m+1 . To obtain the linear system from which λ m+1 is the solution, it suffices to express that σ m+1 has its normal component continuous across each interior edge e of the triangulation T h .

Proposition 2.5. The Lagrange multiplier λ m+1 is the solution of the linear system obtained by associating to each interior edge e of the triangulation the equation:

1 K1 |σ m K1 | q-2 dx |e|λ m+1 |e - |d 1 | 2 + |e| 2 -|g 1 | 2 2|e| λ m+1 |d1 - |g 1 | 2 + |e| 2 -|d 1 | 2 2|e| λ m+1 |g1 + 1 K2 |σ m K2 | q-2 dx |e|λ m+1 |e - |d 2 | 2 + |e| 2 -|g 2 | 2 2|e| λ m+1 |d2 - |g 2 | 2 + |e| 2 -|d 2 | 2 2|e| λ m+1 |g2 = 1 2|K 1 | K1 f (y)dy (x -x m K1 ). n K1|e + 1 2|K 2 | K2 f (y)dy (x -x m K2 ). n K2|e (2.99)
In equation (2.99), K 1 and K 2 denote the two triangles of T h sharing in common the edge e. d 1 and g 1 (resp. d 2 and g 2 ) denote the two other sides than e of the triangle K 1 (resp. K 2 ) (for an arbitrary edge b of T h , |b| denotes its length). In the right-hand side of equation (2.99), x denotes an arbitrary point of the edge e (usually we will take its midpoint).

Proof. Before beginning really the proof, let us explain why x is an arbitrary point of e in the right-hand side of (2.99). Let us take two arbitrary points x ′ and x ′′ on e. Then (x ′ -x m Ki ). n Ki|e = (x ′′ -x m Ki ). n Ki|e (2.100) for i = 1, 2 as x ′ -x ′′ is orthogonal to n Ki|e . Now let us prove (2.99). To do that, let us express that

σ m+1 |K1 . n K1|e + σ m+1 |K2
. n K2|e = 0 using formula (2.91). To alleviate the notations, we will write in the sequel n e Ki for n Ki|e , (i = 1, 2). We obtain from (2.91):

∂K1 λ m+1 n K1 . n e K1 ds K1 |σ m | q-2 dx + ∂K2 λ m+1 n K2 .n e K2 ds K2 |σ m | q-2 dx = K1 f dx 2|K 1 | (x -x m K1 ). n e K1 + K2 f dx 2|K 2 | (x -x m K2 ). n e K2
(2.101) for every x ∈ e; but due to (2.100) it suffices to express it for one x ∈ e. Decomposing each of the integral above on ∂K i (i = 1, 2) as sums of integrals on each side, we obtain:

1 K1 |σ m | q-2 dx |e|λ m+1 |e -|d 1 | cos α 1 λ m+1 |d1 -|g 1 | cos β 1 λ m+1 |g1 + 1 K2 |σ m | q-2 dx |e|λ m+1 |e -|d 2 | cos α 2 λ m+1 |d2 -|g 2 | cos β 2 λ m+1 |g2 = K1 f dx 2|K 1 | (x -x m K1 ). n e K1 + K2 f dx 2|K 2 | (x -x m K2 ). n e K2 , (2.102)
where α i (resp. β i ) denotes the angle in the triangle K i between the sides d i and e (resp. g i and e), (i = 1, 2). Using the law of cosines we have:

cos α i = |d i | 2 + |e| 2 -|g i | 2 2|d i ||e| and cos β i = |g i | 2 + |e| 2 -|d i | 2 2|g i ||e| ,
and replacing in (2.102) we obtain (2.99).

Numerical Experiments

The following experiments will underline and confirm our theoretial predictions. These numerical tests consists in solving the mixed formulation of the p-Laplace equation (0.1) for the value of f given by (0.7) with ε = 0.05 on the unit square Ω = ]-1/2, 1/2[ 2 , using the Raviart-Thomas finite element of degree 0. The correponding exact solution is given by (0.6) and displayed on Figure 1. The computations are performed for p = 1.8 and p = 3.0, each of these two values ensuring that f ∈ L q (Ω), with 1 p + 1 q = 1. Firstly, with the knowledge of the exact solution u given by (0.6), the order of convergence (function of p) is verified, using theorem 3.2 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] for p = 1.8 and theorem 3.1 of [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF] for p = 3.0 (one verifies that (σ , u) ∈ W 1,q (Ω) 2 × W 1,p (Ω) for p = 1.8 or p = 3.0). Secondly, the ratios between the errors and the estimators are plotted, to verify the reliability of the estimators. The numerical tests are performed using quasi-uniform meshes similar to the one displayed on Figure 2, and characterized by the number of nodes n in each direction of space. If we note h the diameter of any triangle and DoF the number of degree of freedom associated to the mesh (i.e. the number of segments), we have h = √ 2/(n -1) and DoF = (n -1)(3n -1).

Order of convergence of the methods

Let us quantify how the numerical solution converges towards the exact one. To this end, we plot the curves :

• ||σ h -σ || 0,q,Ω as a function of DoF (see Fig. 3),

• ||u -u h || 0,p,Ω as a function of DoF (see Fig. 4).

The predicted orders of convergence in h are given in [START_REF] Farhloul | On a mixed finite element method for the p-Laplacian[END_REF], and summarized in table For p = 1.8, the numerical results always give slightly better orders of convergence asymptotically whatever the error considered. For p = 3.0, the good order of convergence is asymptotically achieved for the error on σ h , and a better order of convergence than the one expected by the theory is obtained for the error on u h .

Reliability of the estimators

Proposition 1.5 proves that ||r|| = η 2 , and theorem 1.2 allows to bound ||R|| by a constant times the error estimator η 1 . As a consequence, following proposition 1.3 (resp. proposition 1.4), we plot on figure 5 the ratio r 1 (resp. r 2 ) defined for p = 1.8 by :

             r 1 =
||σ h -σ || 0,q,Ω η p/q 1 + η 2/q 2 + (η 2 η 1 ) 1/q , r 2 = ||u -u h || 0,p,Ω

η 1 + η 2 + η p/2 1
.

Similarly, following proposition 1.2 (resp. corollary 1.1), we plot on figure 6 the ratio r 3 (resp. r 4 ) defined for p = 3.0 by :

           r 3 = ||σ h -σ || 0,q,Ω η 2 1 + η q 2 + η 2 η 1 , r 4 = ||u -u h || 0,p,Ω η 1 + (η 2
1 + η q 2 + η 2 η 1 ) 1/p . The meshes used are the same than the ones involved in part (3.1). We see that whatever the ratio r i chosen (1 ≤ i ≤ 4), the ratio is asymptotically constant. In other words, it appears that all these ratios are asymptotically independent of the number of degrees of freedom DoF associated to the mesh. As a consequence, it shows that the estimators are actually reliable.

) where a m+1 K∈ R 2 |K

 m+12 and c m+1 K ∈ R, we obtain the following expressions of σ m+1 K in terms of the Lagrange multiplier λ m+1 h |∂K :

Remark 1 . 1 .

 11 For a, b ∈ R + , we use the notation a b to mean that there exists a positive constant C such that a C b. Similarly, we write a ∼ b to mean that there exists two strictly positive constants C 1 and C 2 such that C 1 b ≤ a ≤ C 2 b.

)

  Thus it remains to prove inequality (1.26) when |x | = |y |. x and y playing a symmetric role in inequality (1.26), we may suppose |x | > |y |. Considering y and λx with λ = |y | |x | so that |λx | = |y |, it follows by the previous step of our work that

  .31) Applying (1.31) with α = |y |, β = |x | and s = p -1, we get:

Fig. 1 .

 1 Fig. 1. The exact solution u.

Fig. 2 . 8 Fig. 3 . 4 Fig. 4 .

 28344 Fig. 2. Quasi-uniform mesh with n = 5.

Fig. 5 .

 5 Fig. 5. r 1 (left) and r 2 (right) as a function of DoF for p = 1.8.

Fig. 6 .

 6 Fig. 6. r 3 (left) and r 4 (right) as a function of DoF for p = 3.0.

  Under the hypothesis of proposition 2.2 on σ m , for every triangle K of the triangulation T h we have:

							.95)
	By the second equation of (2.88):
	By (2.93), we have also	div σ m+1	|K = -	1 |K| K	f dx.	(2.96)
		div σ m+1	|K = 2c		(2.97)
	From (2.93), (2.95), (2.96) and (2.97) follows our claim (2.91).
	For u m+1 |K , we have the following formula:
	Proposition 2.4.

Table 1 .

 1 1. The corresponding order of convergence in DoF are mentionned into brackets (DoF = O(h -2 )).||σ h -σ || 0,q,Ω ||u -u h || 0,p,ΩThe predicted orders of convergence in h (and in DoF ).

	p = 1.8	8 9	-	4 9	1	-	1 2
	p = 3	3 4	-	3 8	1 2	-	1 4

Remark 3.1. Our computations show that the effectivity index, which is equal to 1/r i , 1 ≤ i ≤ 4, is very high for some of the examples. This fact is due to the presence of unknown constants that enter into the calculation. Using recovery type estimators, one can adjust these constants properly to match the predictions of the residual and recovery estimators. Remark 3.2. Computations have been performed with the use of the code Simula+ developed by the LAMAV laboratory (University of Valenciennes, France) and the LPMM laboratory (ENSAM of Metz, France).


Emmanuel Creusé
email: emmanuel.creuse@univ-valenciennes.fr

Serge Nicaise
email: serge.nicaise@univ-valenciennes.fr

Discrete compactness for a discontinuous Galerkin approximation of Maxwell's system

Keywords: DG method, Maxwell's system, discrete compactness, eigenvalue approximation AMS Subject Classification: 65N25, 65N30

In this paper we prove the discrete compactness property for a discontinuous Galerkin approximation of Maxwell's system on quite general tetrahedral meshes. As a consequence, a discrete Friedrichs inequality is obtained and the convergence of the discrete eigenvalues to the continuous ones is deduced using the theory of collectively compact operators. Some numerical experiments confirm the theoretical predictions.

Introduction

In this paper we deal with the discrete compactness property for a discontinuous Galerkin approximation of Maxwell's system on quite general tetrahedral meshes. This property has been studied for standard Galerkin approximation for a quite large family of edge elements on two and three dimensional 1 domains [START_REF] Kikuchi | On a discrete compactness property for the Nédélec finite elements[END_REF][START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF][START_REF] Boffi | Discrete compactness for p and hp 2d edge finite elements[END_REF][START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF][START_REF] Demkowicz | Maxwell eigenvalues and discrete compactness in two dimensions[END_REF]. But to our knowledge this property is not yet proved for the discontinuous Galerkin method. We here concentrate on the interior penalty method introduced in [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] (see also [START_REF] Lohrengel | A discontinuous Galerkin method on refined meshes for the 2d time-harmonic Maxwell equations in composite materials[END_REF] for two-dimensional domains).

The success of DG methods is today well recognized and is mainly due to its flexibility in the choice of the approximation space, since it allows the use of meshes with hanging nodes and local spaces of different orders. This renders this method well suited for h -p adaptivity. Furthermore the implementation of DG methods is quite easy by using elementwise shape functions, a great advantage for elements of high order, such a property is not satisfied for standard edge elements.

Our proof of the discrete compactness property is based on the same property for the standard Galerkin approximation proved in [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF] and the use of a decomposition of the discontinuous approximation space into a continuous one and its orthogonal for an appropriate inner product similar to [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF] (but different from the one used in these references). The discrete Friedrichs inequality follows from this discrete compactness property and a contradiction argument. The convergence of the discrete eigenvalues to the continuous ones is deduced using the theory of collectively compact operators [START_REF] Anselone | Collectively compact operator approximation theory[END_REF], which requires pointwise convergence of the sequence of the discrete operators (see [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF] for the use of this approach for the standard Galerkin approximation of Maxwell's system). In our case, the collectively compact property is deduced from the discrete compactness property and the pointwise convergence is obtained by introducing mixed formulations and using a variant of the second Strang lemma.

Let us notice that we restrict ourselves to the h-version of the method, even if we consider a discontinuous Galerkin method for polynomials of order less than k, for any positive integer k, we do not estimate the dependence of the constant with respect to k. The extension of our analysis to the so-called h -p version needs more investigations (see [START_REF] Boffi | Discrete compactness for the hp version of rectangular edge finite elements[END_REF] for rectangular meshes). The convergence analysis of the DG method for the source problem by using the results of this paper also requires some further analysis.

For the numerical experiments, since the null space of the operator is relatively large, we have used a discrete regularization method that allows us to work in the setting of positive definite matrices (see [START_REF] Arbenz | Eigenvalue solvers for electromagnetic fields in cavities[END_REF][START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF] for the standard edge elements).

The paper is organized as follows. Section 2 introduces the problem and some notation. The discretization and the discontinuous Galerkin method are given in Section 3. Section 4 contains the main results of the paper, namely the proof of the discrete compactness property and of the discrete Friedrichs inequality. The discrete eigenvalue problem is presented in section 5. In section 6, we first recall some results about collectively compact operators and then used them to deduce the convergence of the eigenvalues. Finally some numerical experiments are presented in section 7, that confirm the theoretical predictions. The obtained numerical eigenvalues are furthermore comparable with the ones given in the existing literature.

Preliminaries and notation

Let us fix a bounded domain Ω of R 3 with a Lipschitz polyhedral boundary Γ. We denote by n the unit outward normal vector along Γ. For the sake of simplicity we assume that Γ is simply connected.

On this domain we consider the following eigenvalue problem: Find an electric field E = 0 and an electric eigenvalue λ such that

curl curl E = λE in Ω, (1) div E = 0 in Ω, (2)
E × n = 0 on Γ. (3)
This eigenvalue problem has been studied extensively in the literature, see e. g. [START_REF] Leis | Initial boundary value problems in Mathematical Physics[END_REF][START_REF] Monk | Finite element methods for Maxwell's equations[END_REF]. It is well known that under the above assumptions, the electric eigenvalues are real and positive and that the set of eigenvalues is discrete. Moreover if λ is an electric eigenvalue, then the eigenspace is finite dimensional.

For further purposes, we recall the variational formulation of this problem. It requires the use of the following (real) function spaces:

H 0 (curl; Ω) = {u ∈ L 2 (Ω) 3 : curl u ∈ L 2 (Ω) 3 and u × n = 0 on Γ}, X = {u ∈ H 0 (curl; Ω) : div u = 0 in Ω},
both spaces being equipped with the norm

u H 0 (curl;Ω) = ( u 2 + curl u 2 ) 1/2 ,
where here and below • means the L 2 (Ω) 3 (or L 2 (Ω) according to the context) norm. Furthermore for any u, v ∈ L 2 (Ω) 3 , we will denote

(u, v) = Ω u • v,
their standard inner product.

In view of (2), the weak formulation of the electric eigenvalue problem is to find E ∈ X, E = 0 and a real number λ such that

(curl E, curl v) = λ(E, v), ∀v ∈ X. (4)
Let us reformulate this problem in a more operator form. More precisely, introduce the operator

A : L 2 (Ω) 3 → X : u → Au,
where Au is the unique element in X that satisfies

(curl Au, curl v) = (u, v), ∀v ∈ X. ( 5
)
This problem is meaningful since the bilinear form a(u, v) = (curl u, curl v) is coercive on X [START_REF] Krizek | On the validity of Friedrichs' inequalities[END_REF][START_REF] Monk | Finite element methods for Maxwell's equations[END_REF] and therefore problem (5) has a unique solution by the Lax-Milgram Lemma.

Using the operator A, the eigenvalue problem (4) is equivalent to the problem of finding E ∈ X, E = 0 such that AE = µE, [START_REF] Birman | L 2 -theory of the Maxwell operator in arbitrary domains[END_REF] with µ = 1 λ . Let us finish this introduction with some notation used in the whole paper: For shortness the L 2 (D) 3 -norm will be denoted by • D . In the case D = Ω, we will drop the index Ω. The usual norm and seminorm of H t (D) are denoted by • t,D and | • | t,D .

Discontinuous Galerkin discretization

In this section, we introduce the interior penalty DG discretization of problem (4). To this end, let us introduce the next notation.

The domain Ω is discretized by a discrete family of conforming meshes T h , h ∈ Λ, made of tetrahedra. The discrete set Λ is equal to {h n : n ∈ N}, where we suppose that the sequence (h n ) n∈N is strictly decreasing and satisfies h n > 0, for all n and h n → 0 as n → ∞. In other words the meshes are built progressively and become finer and finer.

The family is supposed to be regular in Ciarlet's sense [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], i.e. there exists σ > 0 such that the ratio

h T ρ T ≤ σ, ∀T ∈ T h , h ∈ Λ,
where we recall that h T is the diameter of T and ρ T is the diameter of the largest ball inscribed into T . As usual h = max T ∈T h h T .

As in [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF], we further assume that the meshes are weakly quasi-uniform, i.e., if

h min = min T ∈T h h T ,
then there exists µ ⋆ ∈ (0, 1) such that

hh -µ ⋆ min → 0 as h → 0.
As pointed out in [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF], quasi-uniform meshes are weakly quasi-uniform. Moreover graded meshes of Raugel's type used in the presence of corner and edge singularities (see for instance [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF]) are also weakly quasi-uniform.

Elements will be denoted by T , its faces are denoted by E. The set of all (interior and boundary) faces of the mesh will be denoted by E. The measure of an element or face is denoted by |T | := meas 3 (T ) and |E| := meas 2 (E), respectively. For any face E we will denote by h E its diameter (which is equivalent to the diameter h T of any element T containing E by the regularity of the mesh). For each element T ∈ T h , denote by n T the unit outward normal vector along ∂T .

For our further analysis we need to define some jumps and means through any E ∈ E of the mesh. For E ∈ E such that E ⊂ Ω, denote by T + and T -the two elements of T h containing E. Let v be a vector-valued function defined on T + ∪ T -, and which is in H 1 inside each element T ± . We denote by v ± the traces of v on E taken from T ± , respectively. Then we define the mean of v on E by

{{v}} = v + + v - 2 ,
while its tangential jump on E is defined as follows:

[[v]] T = n T + × v + + n T -× v -.
For a boundary face E, i. e., E ⊂ ∂Ω, there exists a unique element T + ∈ T h such that E ⊂ ∂T + . Therefore the mean and the tangential jumps are defined as before by taking v -= 0.

Following [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF], we consider the following discontinuous Galerkin approximation of the continuous eigenvalue problem: Given a mesh T h and a polynomial degree k ≥ 1, we consider the approximation space

V h = {v ∈ L 2 (Ω) 3 : v |T ∈ P k (T ) 3 , ∀T ∈ T h },
where P k (T ) denotes the space of polynomials of total degree at most k on T .

On this space, we define the discontinuous curl operator curl h :

(curl h u) |T = curl(u |T ), ∀T ∈ T h ,
and the bilinear form a h (., .):

a h (u, v) := (curl h u, curl h v) - E [[u]] T • {{ curl h v}} - E [[v]] T • {{ curl h u}} + α E h -1 [[u]] T • [[v]] T ,
where the positive parameter α is the interior penalty stabilization parameter and that will be chosen large enough to ensure coerciveness of the bilinear form a h (see Lemma 5.1 below). Here we use the notation

E ϕ := E∈E E ϕ(x)ds(x),
and h is the function defined on each face

E ∈ E by h(x) = h E , ∀x ∈ E.
The discontinous Galerkin approximation of the eigenvalue problem ( 4)

reads now: Find E h ∈ V h , E h = 0 and λ h ∈ R such that a h (E h , v h ) = λ h (E h , v h ), ∀v h ∈ V h . (7)

The discrete Friedrichs inequality

Following [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF] we consider the space

V (h) = H 0 (curl; Ω) + V h ,
that we equipped with the norm

u 2 DG,h = u 2 + |u| 2 DG,h ,
where the semi-norm is defined by

|u| 2 DG,h = curl h u 2 + E h -1 |[[u]] T | 2 . Now we consider V c h = V h ∩ H 0 (curl; Ω) and decompose V h into V h = V c h ⊕ V ⊥ h , (8)
where the decomposition is orthogonal with respect to the following inner product

(u, v) DG,h new = h -2 (u, v) + (curl h u, curl h v) + E h -1 [[u]] T • [[v]] T .
In other words, V ⊥ h is the orthogonal complement of V c h into V h for this inner product. Note that our orthogonal decomposition is different from the one from [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF] but is motivated by the following result proved in Proposition 4.5 of [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] (or in the Appendix of [START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF]):

Theorem 4.1 There exists a positive constant C (independent of h) such that for any v ∈ V h , there exists v c ∈ V c h satisfying v -v c DG,h ≤ C E h -1 |[[v]] T | 2 1/2 , v -v c ≤ Ch E h -1 |[[v]] T | 2 1/2 .
The following approximation property and norm equivalence (compare with Proposition 4.6 from [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] and Theorem 2 of [START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF]) will be the basis of our analysis below.

Corollary 4.2 There exists a positive constant C (independent of h) such that for all v ∈ V h , the orthogonal projection

π h v of v on V c h wrt to the inner product (•, •) DG,h new satisfies v -π h v DG,h new ≤ C E h -1 |[[v]] T | 2 ) 1/2 , where • DG,h new = (•, •) 1/2 DG,h new .
Proof: The requested estimate directly follows from Theorem 4.1 because it shows that there exists

v c ∈ V c h such that v -v c DG,h new ≤ C E h -1 |[[v]] T | 2 ) 1/2
,

and clearly v -π h v DG,h new ≤ v -v c DG,h new . Note that this corollary shows in particular that v -π h v ≤ Ch E h -1 |[[v]] T | 2 ) 1/2
, which is not the case if we take the orthogonal projection wrt the natural inner product associated with • DG,h .

Corollary 4.3 There exists a positive constant C (independent of h) such that C v h DG,h ≤ v h V ⊥ h ≤ v h DG,h , ∀v h ∈ V ⊥ h , where the norm • V ⊥ h is defined v h 2 V ⊥ h := E h -1 |[[v h ]] T | 2 , ∀v h ∈ V ⊥ h .
Proof: For v h ∈ V ⊥ h , its projection π h v h = 0 and by Corollary 4.2, we obtain

v h DG,h new ≤ C v h V ⊥ h , for some C > 0. The proof is completed by the simple estimate v h DG,h ≤ max{1, h 0 } v h DG,h new .
Let us finally introduce the spaces:

S h = {p h ∈ H 1 0 (Ω) : p h|T ∈ P k+1 (T ), ∀T ∈ T h }, X h = {v h ∈ V h : (v h , ∇p h ) = 0, ∀p h ∈ S h }.
The space X h corresponds to the subspace of V h of functions which are "discrete divergence free". Note further that

∇S h ⊂ V c h . (9)
Now we are ready to prove a property which is similar to the standard discrete compactness property.

Theorem 4.4 Let (h n ) n∈N ⊂ Λ N satisfy h n → 0 as n goes to infinity. Let (u n ) n∈N be a sequence such that u n belongs to X hn and satisfying

u n DG,hn ≤ C, ∀n ∈ N, (10)
for some C > 0. Then using the orthogonal decomposition

u n = u c n + u ⊥ n (11) from (8), with u c n ∈ V c hn and u ⊥ n ∈ V ⊥ hn . Then u c n ∈ X c hn := {v h ∈ V c hn : (v h , ∇p hn ) = 0, ∀p hn ∈ S hn }, (12)
and there exists a subsequence, still denoted by (u n ) n∈N , and u ∈ X such that

u n → u strongly in L 2 (Ω) 3 as n → ∞, ( 13
)
u c n → u weakly in H 0 (curl; Ω) as n → ∞. ( 14
)
Proof: As u n ∈ X hn , it satisfies (u n , ∇q) = 0, ∀q ∈ S hn , and by [START_REF] Chatelin | Spectral approximation of linear operators[END_REF], we get

(u c n , ∇q) + (u ⊥ n , ∇q) = 0, ∀q ∈ S hn .
For q ∈ S hn , by [START_REF] Boffi | Discrete compactness for p and hp 2d edge finite elements[END_REF] we have curl hn ∇q = curl ∇q = 0 and [[∇q]] T = 0 on each edge E, therefore

h -2 n (u ⊥ n , ∇q) = (u ⊥ n , ∇q) DG,hnnew = 0, ∀q ∈ S hn ,
again thanks to [START_REF] Boffi | Discrete compactness for p and hp 2d edge finite elements[END_REF]. The two above identities imply that (u c n , ∇q) = 0, ∀q ∈ S hn , which shows the inclusion [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Now the decomposition (11) yields

u c n DG,hn ≤ u n DG,hn + u ⊥ n DG,hn , ∀n ∈ N,
and by Corollary 4.3, there exists a positive constant

C (independent of h n ) such that u c n DG,hn ≤ u n DG,hn + C u ⊥ n V ⊥ hn , ∀n ∈ N. As [[u ⊥ n ]] T = [[u n ]] T , we deduce that u c n DG,hn ≤ (1 + C) u n DG,hn , ∀n ∈ N.
By [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] we then have

u c n H 0 (curl;Ω) = u c n DG,hn ≤ C, ∀n ∈ N.
Consequently the sequence (u c n ) n∈N is bounded in H 0 (curl; Ω). Hence Theorem 4.1 of [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF] (also valid for the second family of edge elements on tetrahedra, see the beginning of section 3 of [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF]) yields a subsequence, still denoted by (u c n ) n∈N , and u ∈ X satisfying ( 14) and

u c n → u strongly in L 2 (Ω) 3 as n → ∞. ( 15
)
But according to Corollary 4.2,

u ⊥ n ≤ Ch n u n DG,hn ≤ Ch n , which shows that u ⊥ n → 0 strongly in L 2 (Ω) 3 as n → ∞. (16)
The property [START_REF] Costabel | Singularities of electromagnetic fields in polyhedral domains[END_REF] is proved.

This "discontinuous discrete compactness" result allows us to prove the discrete Friedrichs inequality: Theorem 4.5 There exists a positive constant C > 0 such that for all h ∈ Λ, one has

u h ≤ C|u h | DG,h , ∀u h ∈ X h . ( 17
)
Proof: Assume that (17) does not hold, then there exists a sequence u n ∈ X hn , n ∈ N such that

u n = 1, ∀n ∈ N, ( 18
)
|u n | DG,hn = 1 n , ∀n ∈ N. ( 19
)
By the definition of the semi-norm | • | DG,hn , we then have

E h -1 |[[u n ]] T | 2 ≤ 1 n 2 .
Using the (orthogonal) decomposition [START_REF] Chatelin | Spectral approximation of linear operators[END_REF] of u n , we obtain

E h -1 |[[u ⊥ n ]] T | 2 ≤ 1 n 2 , because [[u c
n ]] T = 0. Now with the help of Corollary 4.3, we deduce that

u ⊥ n DG,hn ≤ C 1 n , (20)
for some positive constant C 1 .

From ( 18) and ( 19), the sequence (u n ) n∈N satisfies

u n DG,hn ≤ C 2 , ∀n ∈ N,
for some positive constant C 2 . Therefore by Theorem 4.4, there exists a subsequence, still denoted by (u n ) n∈N , and u ∈ X satisfying ( 13) and [START_REF] Dauge | Benchmark computations for Maxwell equations for the approximation of highly singular solutions[END_REF]. By ( 13) and ( 18), we deduce that

u = 1. ( 21
)
Now from [START_REF] Dauge | Benchmark computations for Maxwell equations for the approximation of highly singular solutions[END_REF], we may write

curl u ≤ lim n→∞ curl u c n .
But the triangular inequality leads to

curl u c n ≤ curl hn u n + curl hn u ⊥ n ≤ |u n | DG,hn + u ⊥ n DG,hn
and by [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] and [START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF], we conclude that curl u c n → 0 as n → ∞.

Hence curl u = 0 in Ω, and since u belongs to X, u = 0 by the standard Friedrichs inequality. This is a contradiction with [START_REF] Kikuchi | On a discrete compactness property for the Nédélec finite elements[END_REF].

The discrete eigenvalue problem

Now we can come back to the discrete eigenvalue problem [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] and show its well posedness.

Lemma 5.1 There exists a parameter α min > 0, independent of the mesh size h such that for α ≥ α min > 0, and h ∈ Λ, the bilinear form a h is strongly coercive on X h , namely there exists β > 0 independent of the mesh size h such that

a h (u h , u h ) ≥ β u h 2 DG,h , ∀u h ∈ X h .
Proof: Using standard inverse estimates (see Lemma 3.1 of [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF]), for α ≥ α min > 0, there exists β 0 > 0 independent of the mesh size h such that

a h (u h , u h ) ≥ β 0 |u h | 2 DG,h , ∀u h ∈ V h .
The conclusion follows from Theorem 4.5.

In view of this lemma, we reduce problem [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] to the space X h .

Lemma 5.2 If λ h = 0, then (7) is equivalent to find E h ∈ X h solution of a h (E h , v h ) = λ h (E h , v h ), ∀v h ∈ X h . ( 22
)
Proof: ⇒ Let E h ∈ V h be a solution of [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF]. We only need to show that it belongs to X h . To this end, we simply take v h = ∇q h with q h ∈ S h as test function in [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] to obtain

a h (E h , ∇q h ) = λ h (E h , ∇q h ).
Since a h (E h , ∇q h ) = 0 and λ h = 0, we conclude that (E h , ∇q h ) = 0, ∀q h ∈ S h .

⇐ Let E h ∈ X h be a solution of [START_REF] Krizek | On the validity of Friedrichs' inequalities[END_REF]. As X h is a closed subspace of V h , any v h ∈ V h may be split up v h = v 0 h + ∇q h , where v 0 h belongs to X h and q h in S h . Consequently, we obtain

a h (E h , v h ) = a h (E h , v 0 h ) = λ h (E h , v 0 h ) = λ h {(E h , v h ) -(E h , ∇q h )} = λ h (E h , v h ),
using successively the property a h (E h , ∇q h ) = 0, the identity ( 22) and the fact that E h ∈ X h .

From Lemma 5.1, the eigenvalues problem ( 22) has only positive eigenvalues. Hence from Lemma 5.2, problem (7) has only nonegative eigenvalues. Thus from the computational point of view, we can solve problem [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] and keep the positive eigenvalues, since only these ones have a physical meaning. Owing to Lemma 5.2, our theoretical analysis can be based on [START_REF] Krizek | On the validity of Friedrichs' inequalities[END_REF].

As in the continuous case, we rewrite problem [START_REF] Krizek | On the validity of Friedrichs' inequalities[END_REF] as a discrete operator equation. Namely we introduce the discrete operator

A h : L 2 (Ω) 3 → X h : u → A h u, where A h u ∈ X h is the unique solution of a h (A h u, v h ) = (u, v h ), ∀v h ∈ X h . ( 23
)
Thanks to Lemma 5.1 this problem is well-posed and hence A h is well-defined.

Lemma 5.3 E h ∈ X h is solution of (22) for λ h > 0 if and only if E h ∈ L 2 (Ω) 3 satisfies A h E h = µ h E h , (24)
with µ h = 1 λ h . Proof: ⇒ Let E h ∈ X h be a solution of [START_REF] Krizek | On the validity of Friedrichs' inequalities[END_REF]. By the definition of A h , we may write for any

v h ∈ X h a h (A h E h , v h ) = (E h , v h ) = 1 λ h a h (E h , v h ).
This leads to the conclusion since a h (•, •) is an inner product on X h . ⇐ Let E h ∈ L 2 (Ω) 3 be a solution of [START_REF] Lohrengel | A discontinuous Galerkin method on refined meshes for the 2d time-harmonic Maxwell equations in composite materials[END_REF]. Then E h = 1 µ h A h E h , which belongs to X h . The conclusion then follows from this identity and (23).

Collective compactness and pointwise convergence

Following [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF], we prove the convergence of the positive eigenvalues of the discrete problem [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] to the positive eigenvalues of the continuous problem (4) by using the theory of collectively compact and pointwise convergent operators [START_REF] Anselone | Collectively compact operator approximation theory[END_REF][START_REF] Osborn | Spectral approximation for compact operators[END_REF][START_REF] Chatelin | Spectral approximation of linear operators[END_REF]. For completeness, we first recall a basis result that will be used later on. Definition 6.1 Let X, Y be two normed spaces and Λ a subset of R. A set A = {A h : X → Y, h ∈ Λ} of linear operators from X to Y is called collectively compact if for each bounded set U of X, the range

A(U) := {A h u : u ∈ U, h ∈ Λ} is relatively compact in Y .
With this definition, we can recall the next convergence result proved in [START_REF] Osborn | Spectral approximation for compact operators[END_REF]. Theorem 6.2 Let X be a Hilbert space of inner product (•, •) X and norm

• X , and let Λ = {h n : n ∈ N} be a discrete subset of R such that h n → 0 as n → ∞. Assume that the set A = {A h : X → X, h ∈ Λ} of linear selfadjoint operators in X is collectively compact. Assume furthermore that there exists a selfadjoint and compact operator A in X such that A h converges pointwise to A (i.e., for all x ∈ X, A h x → Ax in X as h → 0). Let µ be an eigenvalue of A of multiplicity m and denote by {φ j } m j=1 the set of associated orthogonal eigenvectors. Then for any ε > 0 such that the disk B(µ, ε) of radius ε and center µ contains no other eigenvalues of A, there exists h ε such that for all h < h ε , A h has exactly m eigenvalues (repeated according to their multiplicity) in B(µ, ε). Moreover for h < h ε , if we denote by µ h,j , j = 1, • • • , m, the set of the eigenvalues of A h in B(µ, ε), then for all j = 1, • • • , m, there exists a positive constant C such that

|µ -µ h,j | ≤ C m j,l=1 |((A -A h )φ j , φ l )) X | + ( m j=1 (A -A h )φ j 2 X ) 1/2 . ( 25
)
The remainder of the section consists in checking that the set of operators A h defined by [START_REF] Leis | Initial boundary value problems in Mathematical Physics[END_REF] satisfies the hypotheses of the above Theorem. Theorem 6.3 The set A := {A h : h ∈ Λ} is collectively compact.

Proof: Let U be a bounded set of L 2 (Ω) 3 . By the definition of A h , for any

u ∈ U, A h u ∈ X h satisfies a h (A h u, v h ) = (u, v h ), ∀v h ∈ X h .

By Lemma 5.1, we then have

A h u 2 DG,h ≤ β -1 a h (A h u, A h u) = β -1 (u, A h u).
By Cauchy-Schwarz's inequality and the boundedness assumption of U, we obtain

A h u DG,h ≤ C, ∀h ∈ Λ, ( 26
)
for some C > 0 (independent of h).

This estimate and Theorem 4.4 imply that the

W = {A h u : u ∈ U, h ∈ Λ}, is precompact in L 2 (Ω) 3 . Indeed let us fix a sequence (A hn u n ) n∈N in W .
Then either h n → 0 as n → ∞ and since [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF] guarantees that ( 10) holds, we conclude by Theorem 4.4; or h n does not tend to zero as n → ∞ and in that case the set {h n : n ∈ N} is finite, in this last case we conclude by finite dimensionality.

Let us go on with the pointwise convergence:

Theorem 6.4 Let u ∈ L 2 (Ω) 3 , then A h u → Au strongly in L 2 (Ω) 3 , as h → 0,
where we recall that A was defined by [START_REF] Assous | Characterization of the singular part of the solution of Maxwell's equations in a polyhedral domain[END_REF].

Proof: We follow the proof of Theorem 4.5 of [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF] with the necessary adaptations due to the discontinuous scheme. First following some recent ideas from [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF][START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation[END_REF], we introduce the lifting operator

L : V (h) → V h : v → L(v),
where

L(v) is the unique element in V h such that (L(v), w h ) = E [[v]] T • {{w h }}, ∀w h ∈ V h .
Introduce further the modified bilinear form

ãh (u, v) := (curl h u, curl h v) -(L(u), curl h v) -(L(v), curl h u) + α E h -1 [[u]] T • [[v]] T .
We may remark that

ãh (u h , v h ) = a h (u h , v h ), ∀u h , v h ∈ V h , ãh (u, v) = (curl u, curl v), ∀u, v ∈ H 0 (curl; Ω).
Now we transform the discrete problem ( 23) and the continuous problem ( 5) in mixed problems: If A h u ∈ X h is solution of ( 23), then the pair

(A h u, p h ) ∈ V h × S h is the unique solution of a h (A h u, v h ) + b(v h , p h ) = (u, v h ), ∀v h ∈ V h , b(A h u, q h ) = 0, ∀q h ∈ S h , ( 27
)
where b(v, p) = (∇p, v). Indeed the second identity means that A h u ∈ X h and the first one reduces to [START_REF] Leis | Initial boundary value problems in Mathematical Physics[END_REF], when the test functions are taken in X h because b(v h , p h ) = 0 for p h ∈ S h and v h ∈ X h . This mixed problem has a unique solution because a h is coercive on X h , which is precisely the kernel of b in V h , and b satisfies the discrete inf-sup condition

sup v h ∈V h b(v h , p h ) v h DG,h ≥ |p h | 1,Ω , for all p h ∈ S h , since v h = ∇p h yields b(v h , p h ) v h DG,h = |p h | 1,Ω .
Note that in the above problem we can replace a h by ãh . In a similar manner if Au ∈ X is solution of (5), then the pair (Au, p) ∈ H 0 (curl; Ω) × H 1 0 (Ω) is the unique solution of

(curl Au, curl v) + b(v, p) = (u, v), ∀v ∈ H 0 (curl; Ω), b(Au, q) = 0, ∀q ∈ H 1 0 (Ω). (28)
Since V h is not included into H 0 (curl; Ω), we may look at [START_REF] Osborn | Spectral approximation for compact operators[END_REF] as a nonconforming approximation of (28). Therefore by Proposition II.2.16 of [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF],

we may write

Au -Au h DG,h ≤ C inf v h ∈V h Au -v h DG,h + inf p h ∈S h |p -p h | 1,Ω + sup w h ∈V h |ã h (Au, w h ) + b(w h , p) -(u, w h )| w h DG,h .
The two first terms of this right-hand side are so-called approximation error terms and the third one is the consistency error term. Since V c h is included into V h we can estimate the first approximation error term by

inf v h ∈V h Au -v h DG,h ≤ inf v h ∈V c h Au -v h DG,h = inf v h ∈V c h Au -v h H 0 (curl;Ω) .
Since V c h is dense in H 0 (curl; Ω) (see Theorem 4.5 of [START_REF] Monk | Discrete compactness and the approximation of Maxwell's equations in R 3[END_REF]) the above righthand side tends to zero as h → 0, hence

inf v h ∈V h Au -v h DG,h → 0 as h → 0.
Similarly the density of S h in H 1 0 (Ω) leads to inf

p h ∈S h |p -p h | 1,Ω → 0 as h → 0.
It remains to estimate the consistency term. Let us denote by

r h (w h ) = ãh (Au, w h ) + b(w h , p) -(u, w h ).
By the definition of ãh , it can be written

r h (w h ) = (curl Au, curl h w h ) -(L(w h ), curl Au) + b(w h , p) -(u, w h ).
Introducing Π h the L 2 -orthogonal projection on V h , we may write

r h (w h ) = (curl Au, curl h w h ) -(L(w h ), Π h curl Au) + b(w h , p) -(u, w h ).
Elementwise integration by parts leads to

r h (w h ) = T ∈T h T curl curl Au • w h + E [[w h ]] T • {{ curl Au}} -(L(w h ), Π h curl Au) + b(w h , p) -(u, w h ).
As (28) implies that curl curl Au + ∇p = u, we finally arrive at

r h (w h ) = E [[w h ]] T • {{ curl Au -Π h curl Au}}.
Cauchy-Schwarz's inequality and Lemma 4.3 of [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] yield

|r h (w h )| ≤ E∈E [[w h ]] T E {{ curl Au -Π h curl Au}} E ≤ C T ∈T h h t-1/2 T curl Au t,T E∈E,E⊂T [[w h ]] T E , for any t ∈]1/2, 1]
. Using discrete Cauchy-Schwarz's inequality and the definition of the norm • DG,h , we obtain

|r h (w h )| ≤ Ch t curl Au t,Ω w h DG,h .
The Theorem will be proved if one can show that curl Au belongs to H t (Ω) for some t > 1/2. As Au belongs to X and curl curl Au = u -∇p belongs to L 2 (Ω) 3 , we see that curl Au ∈ X T (Ω), where X T (Ω) is defined by (see [START_REF] Costabel | Singularities of electromagnetic fields in polyhedral domains[END_REF])

X T (Ω) := {v ∈ L 2 (Ω) 3 : curl v ∈ L 2 (Ω) 3 , div v ∈ L 2 (Ω) and v • n = 0 on Γ}.
The conclusion follows from the well known embedding X T (Ω) ֒→ H t (Ω) for some t > 1/2, consequence of Theorem 1.1 from [START_REF] Costabel | Singularities of electromagnetic fields in polyhedral domains[END_REF] for instance (see also [START_REF] Birman | L 2 -theory of the Maxwell operator in arbitrary domains[END_REF][START_REF] Assous | Characterization of the singular part of the solution of Maxwell's equations in a polyhedral domain[END_REF]).

Corollary 6.5 Let λ be a positive eigenvalue of the continuous problem (4) of multiplicity m. Then for h small enough the discrete problem [START_REF] Boffi | Fortin operator and discrete compactness for edge elements[END_REF] has exactly m positive eigenvalues λ h,j , j = 1, • • • , m (repeated according to their multiplicity) such that

λ h,j → λ as h → 0, ∀j = 1, • • • , m.
Proof: We have verified in Theorems 6.3 and 6.4 the assumptions of Theorem 6.2 for the set A introduced in Theorems 6.3. Denote by µ = 1 λ the corresponding eigenvalue of A (see section 2). By Theorem 6.2, for h small enough, A h has exactly m positive eigenvalues µ h,j , j = 1, • • • , m, moreover the estimate [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF] implies

|µ -µ h,j | ≤ C m j=1 (A -A h )φ j 2 X 1/2 .
By the pointwise convergence of A h to A, we deduce that

µ h,j → µ as h → 0, ∀j = 1, • • • , m,
and we conclude by Lemma 5.3.

Some numerical tests

For the sake of simplicity, our theoretical results are tested on two-dimensional domains. In fact using the results from [START_REF] Demkowicz | Maxwell eigenvalues and discrete compactness in two dimensions[END_REF][START_REF] Lohrengel | A discontinuous Galerkin method on refined meshes for the 2d time-harmonic Maxwell equations in composite materials[END_REF], our theoretical results also hold for two-dimensional domains, meshed by a regular family (in Ciarlet's sense) of triangles. For our tests we take k = 1 and compute the first fifth eigenvalues for the square domain (0, 1) 2 and the L-shape domain (-1, 1) 2 \ (0, 1) 2 , as well as the first four eigenvalues for the sector domain given in polar coordinates by {(r, θ) ∈ R

2 ; 0 ≤ r ≤ 1, 0 ≤ θ ≤ 3π 2 }.
For the first domain, analytical solutions are known (see for instance [START_REF] Boffi | Discrete compactness for the hp version of rectangular edge finite elements[END_REF]) and the non zero eigenvalues of problem (1) are given by λ j,l = (j 2 + l 2 )π 2 , j + l ≥ 1, in other words, for j ≥ 1, the number j 2 π 2 is a double eigenvalue, while for j, l ≥ 1, j = l, the number (j 2 + l 2 )π 2 is a simple eigenvalue. For the second domain, no analytical solution exists but we compare our results with the computational eigenvalues presented by the benchmark of M. Dauge [START_REF] Dauge | Benchmark computations for Maxwell equations for the approximation of highly singular solutions[END_REF] obtained by using a Galerkin approximation of the Neumann eigenvalues with a geometrical refined mesh near the corner and polynomials of high degree.

For the third domain, analytical solutions are known (see for instance [START_REF] Demkowicz | Maxwell eigenvalues and discrete compactness in two dimensions[END_REF]).

Since the discrete problem (7) has a quite large kernel and therefore does not correspond to an eigenvalue problem for positive definite matrices we use the discrete regularization technique (compare with [START_REF] Arbenz | Eigenvalue solvers for electromagnetic fields in cavities[END_REF][START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF]): First we remark that the kernel K h of problem ( 7) is explicitly known, namely

K h = ∇S h .
Taking into account Lemma 5.2 we further remark that λ h is a nonzero eigenvalue of (7) if and only if there exists E h ∈ X h solution of

a h (E h , v h ) = λ h (E h , v h ), ∀v h ∈ V h . (29)
Now we fix a basis {λ j } J j=1 of S h and a basis {v i } I i=1 of V h (clearly J < I, but J is quite large). Writing

E h = I i=1 x i v i ,
we see that the eigenvalue problem (29) is equivalent to

AX = λ h BX, (30)
where X = (x i ) I i=1 is the unknown vector, A = (a(v i , v i ′ )) 1≤i,i ′ ≤I is the "rigidity" matrix and B = ((v i , v i ′ )) 1≤i,i ′ ≤I is the mass matrix. Using the basis of S h we can say that E h belongs to X h if and only if

C ⊤ X = 0, (31)
where C = (c ij ) 1≤i≤I,1≤j≤J , with c ij = (v i , ∇λ j ).

In summary the eigenvalue problem is reduced to find X satisfying (31) and solution of (30).

At this stage taking an arbitrary weight s > 0, we see that X satisfying (31) and solution of (30) also satisfies

(A + sCC ⊤ )X = λ h BX. (32)
The main point is that A + sCC ⊤ is a positive definite matrix. Indeed we clearly have

X ⊤ (A + sCC ⊤ )X ≥ 0, ∀X ∈ R I .
On the other hand if X ⊤ (A + sCC ⊤ )X = 0, then X ⊤ AX = 0 and C ⊤ X = 0, and consequently X belongs to the kernel of A and to its orthogonal, so that X = 0.

The second advantage of the eigenvalue problem (32) is that the zero eigenvalues of problem (30) become positive but depend linearly on s (while the non zero eigenvalues of problem (30) are unchanged). Consequently the nonzero eigenvalues of interest can be captured by choosing the parameter s large enough (or by taking different values of s). In that respect, our method is similar to the discrete regularization method for standard edge elements [START_REF] Arbenz | Eigenvalue solvers for electromagnetic fields in cavities[END_REF][START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF] and can be compared with the so-called grad-div regularizated method for the continuous Maxwell's problem where the term s Ω div u div v is added to the standard bilinear form of the curl [START_REF] Hazard | On the solution of time-harmonic scattering problems for Maxwell's equations[END_REF][START_REF] Costabel | Singularities of electromagnetic fields in polyhedral domains[END_REF].

The approximated eigenvalues are then produced by solving the problem (32) using Matlab; the matrices being generated by our finite element code Simula+. For both problems we have fixed the penalty parameter α equal to 100 and the weight s = 3200 in order that the first eigenvalues are the expected ones.

In Figure 1 we have plotted the first fifth eigenvalues for the square in comparison with the exact ones for different meshes, while in Figure 2 the error is given as a function of the DoF . This last figure demonstrates the convergence of the discrete eigenvalues to the exact ones with a convergence rate of 1 in DoF -1 (i.e. of order 2 in h). Figures 3 and4 present the same results for the L-shape domain, by taking as exact values the one given by M. Dauge [START_REF] Dauge | Benchmark computations for Maxwell equations for the approximation of highly singular solutions[END_REF]. Here again, a convergence rate of order 1 in DoF -1 (or of order 2 in h) is obtained. We should only notice the particular behaviour of the first eigenvalue convergence rate : The graph of the error for this smallest eigenvalue exhibits a sudden rise after a dramatic fall, and the convergence rate is no more equal to 1 with respect to DoF -1 . This phenomenom occurs for another domain (the slit domain) in [START_REF] Demkowicz | Maxwell eigenvalues and discrete compactness in two dimensions[END_REF], for which the smallest eigenvalue also corresponds to a singular eigenfield. Moreover, also for the L-shape, numerical results obtained in [START_REF] Hesthaven | High-order nodal discontinuous Galerkin methods for the Maxwell eigenvalue problem[END_REF] confirm that there exists a strong correlation between the convergence rate and the regularity of the corresponding eigenfunction. In Figure 5 we have plotted the first four eigenvalues for the sector domain in comparison with the exact ones for different meshes, while in Figure 6 the error is given as a function of the DoF . For this last example, the convergence properties are the same as the ones obtained for the L-shape domain. Note further that these results are comparable with those from [START_REF] Demkowicz | Maxwell eigenvalues and discrete compactness in two dimensions[END_REF].

The tests presented in this section have been performed with the help of the SIMULA+ finite element code (MACS, University of Valenciennes and LPMM, University and ENSAM of Metz, both in France)

Figure 1 :

 1 Figure 1: Eigenvalues wrt DoF for the unit square domain

Figure 2 :

 2 Figure 2: Convergence rates in DoF for the unit square domain

Figure 3 :

 3 Figure 3: Eigenvalues wrt DoF for the L-shape domain

Figure 4 :

 4 Figure 4: Convergence rates in DoF for the L-shape domain

Figure 5 :

 5 Figure 5: Eigenvalues wrt DoF for the sector domain

Figure 6 :

 6 Figure 6: Convergence rates in DoF for the sector domain