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

We construct and analyze a new family of rectangular (two-dimensional) or cubic (three-dimensional) mixed finite elements for the approximation of the acoustic wave equations. The main advantage of this element is that it permits us to obtain through mass lumping an explicit scheme even in an anisotropic medium. Nonclassical error estimates are given for this new element.


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 [12] 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.[12,18,22]).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.[4,11,13]).This model is very efficient compared with the traditional absorbing boundary conditions and can be naturally extended to any first order system [14].
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 [19]) 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 [2] 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.[23,8,9,10,26,24]).
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.[6,5,1]) (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 [3].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 [21] 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 [15] 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.
2.1.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 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 subject to the initial conditions 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 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: where (2.6) The bilinear form a( They satisfy the following properties (see, for instance, [6]): (i) The continuous inf-sup condition: (ii) The coercivity of the form a(•, •) on V ≡ KerB : (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
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: 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 subject to the initial conditions The usual choice consists in taking for X the lowest order Raviart-Thomas element (see [21]), X = RT [0] = P 10 × P 01 , (2.10) and for M the piecewise constants: Here, P k is the space of polynomials of degree not greater than k and P kl is defined by 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 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 N1 ) and [V ] = (V 1 , . . ., V N2 ) the coordinates of p h and v h in the bases B N1 and B N2 .In these bases, problem (2.9) can be written in the following form: (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 [10,25]).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.

Now consider the integral
with ( M q ) q=1,...,N being the quadrature points and (ω q ) q=1,...,N the associated weights to approximate (2.13), we obtain This would lead to a diagonal matrix if 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 where Q 1 is the space of piecewise bilinear functions.This element was initially introduced by Nédélec in [19].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 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 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 [25,9,8] and to Cohen and Monk [10] and Elmkies and Joly [16] 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.
2.3.Extension to higher order and mass lumping.The natural generalization of the lowest order element, presented in the previous section, consists of taking 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 [10] and Tordjman [25], 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 .qh belongs to Q 2k+2 .However, following [7,25,17] 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 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.
3. Analysis of the new mixed finite element for an elliptic problem.Following [5], 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 We know that (3.1) admits a unique solution u ∈ H 1 0 (Ω) and that there exists c > 0 such that As for the time dependent problem, we set p = A −1 (x)∇u, (3.3) and this gives The mixed formulation associated to (3.3) where a(•, •) and b(•, •) are defined by (2.6) and satisfy properties (2.7).As proven in [6], 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 The elliptic problem (3.5) and its approximation (3.6) have been studied by several authors (see, for example, [21,6]), and we know that it admits a unique solution (p h , u h ) in X h × M h with the convergence property when the following assumptions are satisfied: The uniform discrete inf-sup condition: where 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 [21].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) We can easily see that 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 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 The reader can have in mind that for our application we shall have 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 The kernels of B and B t are defined as follows: We make the assumptions where the norm in the quotient space is defined by We shall identify the quotient space M/Ker B t with the orthogonal complement of We are interested in the numerical approximation of the solution (p, u) to the following problem: with f ∈ M ′ the dual space of M .Under these assumptions, we have the following classical result (see [6]).Theorem 3.1.For all f ∈ Im B, problem (3.11) has a unique solution (p, u) in X × (M/Ker B t ).Moreover, Suppose X h ⊂ X and M h ⊂ M are finite dimensional approximation spaces.We consider then the approximate problem and make the following hypotheses: (H1) Orthogonal decomposition of X h : (H2) "Strong" discrete uniform inf-sup condition: there exists a constant c > 0, independent of h, such that (H3) "Weak" coercivity: there exists a constant α > 0, independent of h, such that Remark 3.2.Note that, since X is dense in H, (H4) implies that We can introduce, as for the continuous problem, an operator B h from X h to M h defined by and the kernels of B h (see (3.13)) and of its transpose It may be more convenient to characterize hypothesis (H0) with one of the following equivalent statements: 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 , Theorem 3.2.Under hypotheses (H0)-(H4), problem (3.12) admits a unique solution such that ), and the following convergence result holds: More precisely, Proof.First note that hypothesis (H0) and nonuniform discrete coercivity on the kernel Ker B h , i.e., 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 h − q 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 Using (H1) in (3.16), we can write, with obvious notation, We have used the fact that b(w h , p r h − q 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 h − q r h ) (which would give p r h − q r h X from the continuity of b).Then, by (3.11), We choose now or, by using the orthogonality of X r h and X s h from (H1), Furthermore, from the inequality and from hypothesis (H1) and (H3), (3.20) leads to We deduce the existence of a constant C depending on f, p, a , b , and α such that This gives, using (3.9), To conclude, let us recall that the inf-sup condition (H2) implies (cf.[6]) that inf Finally, it remains to prove estimates for u − u h M/Ker B t h .Let us subtract the first equation of (3.12) from the first equation of (3.11).We get Using this and the inf-sup condition, we have It follows from the triangle inequality that 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
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 [21]; 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
Let us begin by checking hypothesis (H1).We take as approximation spaces We define X s h as the lowest order Raviart-Thomas element RT [0] defined in (2.10) 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) where It is then easy to prove the following lemma.
h can be generated from the four functions: .

Note now that p
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.[20]) we obtain the following theorem (here we have Ker B t h = {0}).Theorem 3.5.The problem (3.6) admits a unique solution 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 where |.| H 1 (resp., |.| (H 1 ) 2 ) denotes the usual semi-norm in H 1 (Ω) (resp., (H 1 (Ω)) 2 ).Remark 3.7.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 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: where a d (•, •) is given by (2.6) after replacing A(x) by A d (x).Indeed, in this particular case, 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 : 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,
We have seen in section 2.3 that the natural generalization of the lowest order element consists of taking and thus to introduce the spaces (2.16) 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 which obviously satisfies the inclusion The main property of the space Ψ k is the following.Lemma 3.6.For any ψ in Proof.For simplicity we consider the reference element We begin by a characterization of Ψ k .Let σ k be the polynomial in one variable of degree k + 1 such that Equivalently, σ k is defined, up to a multiplicative constant, by Then, we claim that To see that, we first remark that (see [21]) RT [k] is generated by vector fields of the form τ 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 As q x and q y are in P k , by the definition of σ k , (ii) To prove (3.23) we first use Green's formula Let us decompose ∂K as as shown in Figure 3.3.Then, we check that for 1 ≤ j ≤ 4 and ∀q ∈ R k (T j ), 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 To conclude, it remains to notice that, if v belongs to Q k+1 , then v| Tj belongs to R k (T j ).
This result suggests that we define X s h and X r h by 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 [21].(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.
We can now apply Theorem 3.1 to the approximation problem (3.6), and by using the usual interpolation results (cf.[20]) 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 Furthermore, if we assume that the solution (p, u) of where 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 We set For w ∈ D(B t ), we have b(w, q) = (B t w, q) H ∀ q ∈ X.
Let us introduce the notation 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 (ii) For all (p, u) ∈ X × D(B t ), there exists a constant C, independent of h, such that where In particular, |||p − p h ||| and u − u 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 h − q r h ) = −a(p, p r h − q r h ) ≡ −(Ap, p r h − q 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 [15] 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 [7,25,17] 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 Or, equivalently, in an operator form, if we assume enough regularity on the solution in time, 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 From the classical theory of ODEs, we have the following result.
Following [12,15] 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 This primitive is unique up to a constant element of Ker B t .
) and there exists a constant C independent of h such that where E h is defined by In Remark 4.2.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 (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 h − q 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 h − q 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 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 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 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): 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,    and we choose as approximate initial conditions the elliptic projection of the exact initial condition, (4.8), so that at time t = 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 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, (recalling that η 1 X = |η 1 | H ∀η 1 ∈ Ker B h ).We set η h = p h − p h .The term |η 1 | H is 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 to        there exists a constant C > 0, independent of h, such that ∀ q s h ∈ X s h , sup We also know that η 2 X = η s h X/Ker B h .Thus, by taking q s h = η 2 in (4.19), sup Now use the second equation of (4.16) to obtain Till now, we have used only the C 1 regularity of the solution.In order to bound ∂ t ( v h − v h ) M , we need C 2 .Indeed, we want to apply (4.11)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 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 We apply (4.6) for k = 1 and get Again, using estimate (4.11) to bound v h − v h M , we easily get the second inequality of (4.9).
• Estimate in X-norm.We write p − p s h X ≤ p − p s h X + p s h − p s h X .
4.2.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.
particular, |||p − p h ||| and u − u 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) ) 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 Let (p, be the solution of (4.2) and (p h , v h ) the solution of the approximate problem (4.3) with the initial conditions .10)In order to prove Theorem 4.3 we need the following lemma.Lemma 4.4.Let (p, v) be the solution of (4.2) and (p h , v h ) be the solution of the approximate problem (4.3) with the initial conditions (4.8).Let Π h (p, u) = ( p s h + p r h , u h ) be the elliptic projection defined in Lemma 4.2, and let v h be defined by (4.7).Set δ h w h ).Furthermore, by taking q h = p h − p h and w h = v h − v h in(4.15)and by adding the two equations, we geta(∂ 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) p h − p h , p h − p h ) + ( v h − v h , v h − v h )) (t).(t) ≥ C | p h − p h | 2 H (t) + v h − v h ≤ C ′ ∂ t (v − v h ) M (t) ≡ C ′ ∂ t δ h M (t), (4.18)where we set δ h = v − v 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 h with v h replaced by ∂ t v h , v h by ∂ t v h , and so on.More precisely, we have∂ t ( v h − v h ) M (t) ≤ CFinally, combining (4.20) and (4.21), we getη 2 X (t) ≤ C ′ ∂ t δ h M (t) + Proof (ofTheorem 4.3).We combine results given in Lemma 4.2 and in Lemma 4.4.•Estimates in H-and M -norms.We have|p − p s h | H + |p r h | H ≤ |p − p s h | H + | p r h | H + | p s h − p s h | H + | p r h − p r h | H .