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

Éliane Bécache, Patrick Joly, Chrysoula Tsogka

<figure><img src="image_1.png" /></figure>

## To cite this version:

Éliane Bécache, Patrick Joly, Chrysoula Tsogka. An Analysis of New Mixed Finite Elements for the Approximation of Wave Propagation Problems. SIAM Journal on Numerical Analysis, 2000, 37 (4), pp.1053 - 1084. 10.1137/S0036142998345499. hal-01443182

HAL Id: hal-01443182

https://hal.science/hal-01443182

Submitted on 22 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
AN ANALYSIS OF NEW MIXED FINITE ELEMENTS FOR THE APPROXIMATION OF WAVE PROPAGATION PROBLEMS*

E. BECACHE, P. JOLY, AND C. TSOGKA

**Abstract.** 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.

**Key words.** mixed finite elements, mass lumping, anisotropic waves

**1. 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 $\Gamma$ 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 = \nabla 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 = \partial_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.

## 2. Presentation of the new mixed finite elements.

**2.1. The model problem: The anisotropic wave equation.** Let $\Omega$ be a bounded domain of $\mathbb{R}^2$ and let $A(x)$ be a positive definite symmetric matrix satisfying

$$ (2.1) \qquad A(x) \xi \cdot \xi \geq \alpha |\xi|^2, \quad \alpha > 0, \quad \forall \xi \in \mathbb{R}^2, \quad \text{a.e. } x \in \Omega. $$

In practice, $\Omega$ will be the simple geometrical domain (let us say a rectangle) we referred to in the introduction. We consider the scalar evolution problem

$$ (2.2) \qquad \left\{ \begin{array}{ll} \text{Find } u : [0, T] \mapsto H_0^1(\Omega) \text{ such that} \\ \frac{\partial^2 u}{\partial t^2} - \text{div} (A^{-1}(x)\nabla u) = f, & f \in C^0(0, T; L^2(\Omega)), \end{array} \right. $$

subject to the initial conditions

$$ u(t=0) = u_0 \in H^1(\Omega); \quad \frac{\partial u}{\partial t}(t=0) = u_1 \in L^2(\Omega). $$

**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] \mapsto H_0^1(\Omega)$.

**REMARK 2.2.** The homogeneous Dirichlet condition on $\partial\Omega$ has been considered for simplicity only. This has nothing to do with a Neumann condition on the boundary of an obstacle interior to $\Omega$ that will be taken into account with the fictitious domain method as explained in the introduction.

Let

$$ p = A^{-1}(x)\nabla u \quad \text{and} \quad v = \frac{\partial u}{\partial t}. $$

Substituting into (2.2) yields

$$ (2.3) \qquad \left\{ \begin{array}{ll} \frac{\partial v}{\partial t} - \text{div } p = f, \\ A \frac{\partial p}{\partial t} - \nabla v = 0 \end{array} \right. $$

with initial conditions

$$ (2.4) \qquad p(0) = p_0 = A^{-1}(x)\nabla u_0; \quad v(0) = v_0 = u_1. $$

**REMARK 2.3.** Throughout this paper, the generic point of $\mathbb{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:

$$ (2.5) \qquad \left\{ \begin{array}{ll} \text{Find } (p, v) : [0, T] \mapsto X \times M \equiv H(\text{div} ; \Omega) \times L^2(\Omega) \text{ such that} \\ \frac{d}{dt}a(p, q) + b(v, q) = 0 & \forall q \in X, \\ \frac{d}{dt}(v, w) - b(w, p) = (f, w) & \forall w \in M, \end{array} \right. $$
where

$$
(2.6) \qquad \left\{
\begin{array}{ll}
a(p, q) = \int_{\Omega} A(x) \, p \cdot q \, dx & \forall (p, q) \in X \times X, \\
b(w, q) = \int_{\Omega} w \operatorname{div} q \, dx & \forall (w, q) \in M \times X, \\
(f, w) = \int_{\Omega} f w \, dx & \forall w \in M.
\end{array}
\right.
$$

The bilinear form $a(\cdot, \cdot)$ (resp., $b(\cdot, \cdot)$) is continuous on $H \times H$ ($H = (L^2(\Omega))^2$) (resp., on $X \times M$). The bilinear form $a(\cdot, \cdot)$ (resp., $b(\cdot, \cdot)$) thus defines a linear continuous operator $\mathcal{A}: H \to H'$ by $\langle \mathcal{A}p, q \rangle_{H' \times H} = a(p, q)$ (resp., $B: X \to M'$ by $\langle Bq, w \rangle_{M' \times M} = b(w, q)$). They satisfy the following properties (see, for instance, [6]):

$$
(2.7) \qquad \left|
\begin{array}{ll}
\text{(i)} & \text{The continuous inf-sup condition:} \\
& \exists c > 0 / \forall w \in M, \exists q \in X / b(w, q) \ge c \|w\|_M \|q\|_X. \\[1em]
\text{(ii)} & \text{The coercivity of the form } a(\cdot, \cdot) \text{ on } V \equiv \operatorname{Ker} B: \\
& \exists \alpha > 0 / \forall p \in V, a(p, p) \ge \alpha \|p\|_X^2.
\end{array}
\right.
$$

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.

**2.2. Presentation of the $Q_1^{\text{div}} - Q_0$ finite element in the lowest order.**

We suppose now that $\Omega$ is a union of rectangles in such a way that we can consider a
regular mesh ($\mathcal{T}_h$) with square elements $K$ of edge $h > 0$. We introduce the following
approximation spaces:

$$
(2.8) \qquad \left\{
\begin{array}{@{}l@{}}
X_h = \{ q_h \in X / \forall K \in \mathcal{T}_h, q_h|_K \in \hat{X} \}, \\
M_h = \{ w_h \in M / \forall K \in \mathcal{T}_h, w_h|_K \in \hat{M} \},
\end{array}
\right.
$$

where $\hat{X}$ (resp., $\hat{M}$) denotes a finite dimensional space of vector (resp., scalar) func-
tions. The discrete problem associated to (2.5) and (2.4) is

$$
(2.9) \qquad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\text{Find } (p_h, v_h) : [0, T] \to X_h \times M_h & \text{such that} \\
\displaystyle\frac{d}{dt}a(p_h, q_h) + b(v_h, q_h) = 0 & \forall q_h \in X_h, \\
\displaystyle\frac{d}{dt}(v_h, w_h) - b(w_h, p_h) = (f, w_h) & \forall w_h \in M_h,
\end{array}
\right.
$$

subject to the initial conditions

$p_h(0) = p_{0,h} ; v_h(0) = v_{1,h}.$

The usual choice consists in taking for $\hat{X}$ the lowest order Raviart–Thomas element
(see [21]),

$$
(2.10) \qquad \hat{X} = RT_{[0]} = P_{10} \times P_{01},
$$
and for $\hat{M}$ the piecewise constants:

$$ (2.11) \qquad \hat{M} = Q_0. $$

Here, $P_k$ is the space of polynomials of degree not greater than $k$ and $P_{kl}$ is defined by

$$ P_{kl} = \left\{ p(x, y) \mid p(x, y) = \sum_{i \le k, j \le l} a_{ij} x^i y^j \right\}. $$

This choice, however, does not lead to an explicit scheme when one considers the evolution problem of anisotropic waves corresponding to (2.9). We introduce here $B_{N_1} = \{\tau_i\}_{i=1}^{N_1}$, $B_{N_2} = \{\phi_i\}_{i=1}^{N_2}$ the bases of $X_h$ and $M_h$, resp., where $N_1 = \dim X_h$ and $N_2 = \dim M_h$. We denote then by $[P] = (P_1, \dots, P_{N_1})$ and $[V] = (V_1, \dots, V_{N_2})$ the coordinates of $p_h$ and $v_h$ in the bases $B_{N_1}$ and $B_{N_2}$. In these bases, problem (2.9) can be written in the following form:

$$ \left\{ \begin{array}{l} \text{Find } (P, V) : [0, T] \longrightarrow \mathbb{R}^{N_1} \times \mathbb{R}^{N_2} \text{ such that} \\ M_p \frac{dP}{dt} + C^T V = 0, \\ M_v \frac{dV}{dt} - CP = F, \\ + \text{initial conditions} \end{array} \right. $$

with

$$ \begin{array}{ll} (i) & (M_p)_{i,j} = (A\tau_i, \tau_j)_{(L^2(\Omega))^2}, \quad 1 \le i, j \le N_1, \\ (ii) & (M_v)_{i,j} = (\phi_i, \phi_j)_{L^2(\Omega)}, \quad 1 \le i, j \le N_2, \\ (iii) & (C)_{i,j} = (\phi_i, \operatorname{div} \tau_j)_{L^2(\Omega)}, \quad 1 \le i \le N_2, 1 \le j \le N_1, \\ (iv) & (F)_j = (f, \phi_j)_{L^2(\Omega)}, \quad 1 \le j \le N_2; \end{array} \qquad (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 $(\tau_I)$ be the $RT_{[0]}$ global basis functions. In Figure 2.1, we represent the four local basis vector functions in $RT_{[0]}(\hat{K})$, $\hat{K}$ being the reference element. Each global basis function $\tau_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 $\tau_I$ on one element $K$ of its support corresponds to a local basis function, denoted $\hat{\tau}_i$, on the reference element $\hat{K}$. In fact $\hat{\tau}_i$ corresponds to $\hat{\tau}_i^x$ (resp., $\hat{\tau}_i^y$) when $\tau_I$ is associated to a vertical (resp., horizontal) edge of the mesh.
<figure><img src="image_3.png" /><figcaption>FIG. 2.1. Local vector basis functions in RT$_{[0]}$($\hat{K}$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2.2. Support of the basis functions in RT$_{[0]}$.</figcaption></figure>

Now consider the integral

$$ (2.13) \quad a(\tau_I, \tau_J) = \sum_{K \in \mathcal{T}_h} \int_K A \tau_I|_K \cdot \tau_J|_K dx = \sum_{K \in \mathcal{T}_h} \mathrm{mes}(K) \int_{\hat{K}} A \hat{\tau}_i \cdot \hat{\tau}_j d\hat{x}. $$

If we use the quadrature formula

$$ \int_{\hat{K}} \hat{f}(\hat{x}) d\hat{x} = \sum_{q=1}^{N} \omega_q \hat{f}(\hat{M}_q), $$

with $(\hat{M}_q)_{q=1,...,N}$ being the quadrature points and $(\omega_q)_{q=1,...,N}$ the associated weights to approximate (2.13), we obtain

$$ \int_{\hat{K}} A \hat{\tau}_i \cdot \hat{\tau}_j d\hat{x} \approx \sum_q \omega_q A \hat{\tau}_i(\hat{M}_q) \cdot \hat{\tau}_j(\hat{M}_q). $$

This would lead to a diagonal matrix if

$$ (2.14) \qquad A \hat{\tau}_i(\hat{M}_q) \cdot \hat{\tau}_j(\hat{M}_q) = 0 \quad \forall i \neq j. $$
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 $\hat{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 $\hat{K}$, the functions $A\hat{\tau}_1^y$ and $\hat{\tau}_1^x$ are no longer orthogonal.

The alternative solution that we propose consists of changing the approximation space $\hat{X}$. Let

$$ (2.15) \qquad \hat{X} = Q_1 \times Q_1, $$

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 $\hat{K}$. We will call this choice (2.15) combined with (2.11) the $Q_1^{div} - Q_0$ element.

<figure><img src="image_5.png" /><figcaption>FIG. 2.3. The local basis functions in $\hat{X} = Q_1 \times Q_1$.</figcaption></figure>

$$ \hat{\tau}_4^y = (0, (1-x)y)^t, \quad \hat{\tau}_3^y = (0, xy)^t $$

$$ \hat{\tau}_4^x = ((1-x)y, 0)^t, \quad \hat{\tau}_3^x = (xy, 0)^t $$

$$ \hat{\tau}_1^x = ((1-x)(1-y), 0)^t, \quad \hat{\tau}_2^x = (x(1-y), 0)^t $$

$$ \hat{\tau}_1^y = (0, (1-x)(1-y))^t, \quad \hat{\tau}_2^y = (0, x(1-y))^t $$

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

$$ \tau_i^x(M_j) = \delta_{ij}(1, 0)^t, \quad \tau_i^y(M_j) = \delta_{ij}(0, 1)^t. $$

Note that, over each element, there are four quadrature points but eight degrees of freedom. For more details on quadrature formulae and mass lumping techniques we refer to the work of Cohen, Joly, and Tordjman for the acoustic wave equation [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 \times 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 $\tau_i^x$ and $\tau_j^y$ 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.
<figure><img src="image_6.png" /><figcaption>FIG. 2.4. The $Q_1 \times Q_1$ element.</figcaption></figure>

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

$$ (2.16) \qquad \left\{ \begin{array}{l} X_h = \{q_h \in X / \forall K \in \mathcal{T}_h, q_h|_K \in Q_{k+1} \times Q_{k+1}\}, \\ M_h = \{w_h \in M / \forall K \in \mathcal{T}_h, w_h|_K \in Q_k\}, \end{array} \right. $$

and we will call it the $Q_{k+1}^{\text{div}} - 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.

<figure><img src="image_7.png" /><figcaption>FIG. 2.5. Degrees of freedom in the $Q_2^{\text{div}} - Q_1$ element.</figcaption></figure>

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
<figure><img src="image_8.png" /><figcaption>FIG. 2.6. Degrees of freedom in the $Q_3^{div} - Q_2$ element.</figcaption></figure>

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 \times 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 $A_{p_h} \cdot q_h$ 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

$$ \hat{X}^{RT} = RT_{[k]} = P_{k+1,k} \times P_{k,k+1} \quad \text{and} \quad \hat{M} = Q_k. $$

It is clear that the new approximate space $X_h$ contains the space $X_h^{RT}$. 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.

**3.1. The elliptic problem.** The elliptic problem we consider here is

$$ (3.1) \qquad \left\{ \begin{array}{l} \text{Find } u \in H_0^1(\Omega) \text{ such that} \\ -\text{div} (A^{-1}(x)\nabla u) = f, \quad f \in L^2(\Omega). \end{array} \right. $$

We know that (3.1) admits a unique solution $u \in H_0^1(\Omega)$ and that there exists $c > 0$ such that

$$ (3.2) \qquad \|u\|_{H^1(\Omega)} \le c \|f\|_{L^2(\Omega)}. $$
As for the time dependent problem, we set

$$ (3.3) \qquad p = A^{-1}(x)\nabla u, $$

and this gives

$$ (3.4) \qquad -\operatorname{div} p = f. $$

The mixed formulation associated to (3.3) and (3.4) is

$$ (3.5) \qquad \left\{ \begin{array}{ll} \text{Find } (p, u) \in X \times M = H(\operatorname{div}; \Omega) \times L^2(\Omega) \text{ such that} \\ a(p, q) + b(u, q) = 0 & \forall q \in X, \\ b(w, p) = -(f, w) & \forall w \in M, \end{array} \right. $$

where $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ are defined by (2.6) and satisfy properties (2.7). As proven in [6], there exists a unique solution $(p, u)$ in $X \times 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/\operatorname{Ker} B^t$, but it is easy to check that for the divergence operator $\operatorname{Ker} B^t = \{0\}$.) For the approximation of this problem, we again consider the finite dimensional spaces $X_h$ and $M_h$ defined by (2.8). The discrete problem associated to (3.5) is

$$ (3.6) \qquad \left\{ \begin{array}{ll} \text{Find } (p_h, u_h) \in X_h \times M_h \text{ such that} \\ a(p_h, q_h) + b(u_h, q_h) = 0 & \forall q_h \in X_h, \\ b(w_h, p_h) = -(f_h, w_h) & \forall w_h \in M_h. \end{array} \right. $$

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 \times M_h$ with the convergence property

$$ (3.7) \qquad (p_h, u_h) \longrightarrow (p, u) \quad \in X \times M $$

when the following assumptions are satisfied:

(i) The uniform discrete inf-sup condition:

$\exists c > 0$ independent of $h$ such that

$\forall w_h \in M_h, \exists q_h \in X_h / b(w_h, q_h) \ge c \|w_h\|_M \|q_h\|_X$;

(3.8) (ii) The coercivity of the form $a(\cdot, \cdot)$ on $V_h$:

$\exists \alpha > 0$ independent of $h$ such that

$\forall p_h \in V_h, a(p_h, p_h) \ge \alpha \|p_h\|_X^2,$

where $V_h = \operatorname{Ker} B_h = \{q_h \in X_h / b(w_h, q_h) = 0 \forall w_h \in M_h\}$.
<figure><img src="image_9.png" /><figcaption>FIG. 3.1. The function $\vec{f}_h$.</figcaption></figure>

These assumptions are satisfied by the couple of spaces $(X_h, M_h)$ defined in (2.8), $\hat{X}$ being the lowest order Raviart-Thomas element (2.10) and $\hat{M}$ the piecewise constants (2.11); see [21]. With our new choice for $\hat{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 $\vec{f}_h$ given by (see Figure 3.1)

$$ \vec{f}_h|_{K_1} = \begin{pmatrix} \left(1 - \frac{x}{h}\right) \left(1 - 2\frac{y}{h}\right) \\ 0 \end{pmatrix} \quad \text{and} \quad \vec{f}_h|_{K_2} = \begin{pmatrix} \left(1 + \frac{x}{h}\right) \left(1 - 2\frac{y}{h}\right) \\ 0 \end{pmatrix}. $$

We can easily see that $\vec{f}_h \in V_h$ (if $\hat{M} = Q_0$) and

$$ a(f_h, f_h) \underset{h \to 0}{\longrightarrow} 0, $$

$$ \|f_h\|_X \underset{h \to 0}{\longrightarrow} \frac{2}{3}, $$

so that we cannot expect to verify (3.8(ii)). However, we will prove (see section 3) that this choice gives a good approximate solution and we will show a new convergence result.

**REMARK 3.1.** In order to preserve (3.8(ii)) one could change the approximation space $\hat{M}$. The natural choice is

$$ \hat{M} = P_1, $$

the key point being that the divergence operator sends $Q_1 \times Q_1$ into $P_1$.<fnref n="1" /> We have eliminated this choice because it is rather expensive in terms of computational time and memory requirements.

**3.2. 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

$$ (3.9) \qquad X \subset H, \quad X \text{ dense in } H \quad \text{and} \quad |\cdot|_H \le \|\cdot\|_X. $$

The reader can have in mind that for our application we shall have

$$ H = (L^2(\Omega))^2, \quad X = H(\operatorname{div}; \Omega) \text{ and } M = L^2(\Omega). $$

Let $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ be two continuous bilinear forms in $H \times H$ and $M \times X$. In the same way as in section 2.1, the bilinear form $a(\cdot, \cdot)$ defines an operator $\mathcal{A}$ in $\mathcal{L}(H)$,

<footnote n="1">¹This implies $V_h = \{q_h \in X_h / \operatorname{div} q_h = 0\} \subset \operatorname{Ker} B$, so that (3.8(ii)) is a consequence of (2.7(ii)).</footnote>
such that $a(p, q) = (Ap, q)_H \forall (p, q) \in H \times H$ and the bilinear form $b(\cdot, \cdot)$ defines an operator $B: X \to M'$ (and its transpose $B^t: M \to X'$) such that

$$ \langle Bp, w \rangle_{M' \times M} = \langle p, B^t w \rangle_{X \times X'} = b(w, p) \quad \forall (p, w) \in X \times M. $$

The kernels of $B$ and $B^t$ are defined as follows:

$$ \left\{ \begin{array}{l} V \equiv \text{Ker } B = \{p \in X / b(w, p) = 0 \ \forall w \in M\}, \\ \text{Ker } B^t = \{w \in M / b(w, p) = 0 \ \forall p \in X\}. \end{array} \right. $$

We make the assumptions

$$ (3.10) \quad \left| \begin{array}{ll} \text{(i)} & \exists c > 0 / \ \forall w \in M, \ \exists q \in X / b(w, q) \ge c \|w\|_{M/\text{Ker } B^t} \|q\|_X, \\ \text{(ii)} & \exists \alpha > 0 / \ \forall p \in V, \quad a(p, p) \ge \alpha \|p\|_X^2, \end{array} \right. $$

where the norm in the quotient space is defined by

$$ \|w\|_{M/\text{Ker } B^t} = \inf_{w_0 \in \text{Ker } B^t} \|w + w_0\|_M. $$

We shall identify the quotient space $M/\text{Ker } B^t$ with the orthogonal complement of $\text{Ker } B^t$

$$ M/\text{Ker } B^t \equiv (\text{Ker } B^t)^\perp \equiv \{w \in M / (v, w)_M = 0 \ \forall v \in \text{Ker } B^t\}. $$

We are interested in the numerical approximation of the solution $(p, u)$ to the following problem:

$$ (3.11) \quad \left\{ \begin{array}{ll} \text{Find } (p, u) \in X \times M \text{ such that} \\ & a(p, q) + b(u, q) = 0 \quad \forall q \in X, \\ & b(w, p) = -\langle f, w \rangle \quad \forall w \in M, \end{array} \right. $$

with $f \in M'$ the dual space of $M$. Under these assumptions, we have the following classical result (see [6]).

**THEOREM 3.1.** For all $f \in \text{Im } B$, problem (3.11) has a unique solution $(p, u)$ in $X \times (M/\text{Ker } B^t)$. Moreover,

$$ \|u\|_{M/\text{Ker } B^t} + \|p\|_X \le C \|f\|_{M'}. $$

Suppose $X_h \subset X$ and $M_h \subset M$ are finite dimensional approximation spaces. We consider then the approximate problem

$$ (3.12) \quad \left\{ \begin{array}{ll} \text{Find } (p_h, u_h) \in X_h \times M_h \text{ such that} \\ & a(p_h, q_h) + b(u_h, q_h) = 0 \quad \forall q_h \in X_h, \\ & b(w_h, p_h) = -\langle f, w_h \rangle \quad \forall w_h \in M_h. \end{array} \right. $$

We set

$$ (3.13) \quad \left\{ \begin{array}{l} V_h(f) = \{q_h \in X_h / b(w_h, q_h) = -\langle f, w_h \rangle \ \forall w_h \in M_h\} = B_h^{-1}(f), \\ V_h = V_h(0) = \text{Ker } B_h \text{ (see below for the meaning of } B_h) \end{array} \right. $$

and make the following hypotheses:
(H0) $\forall f \in \operatorname{Im} B, V_h(f) \neq \emptyset$.

(H1) Orthogonal decomposition of $X_h$:

$$ 
\left\| 
\begin{array}{l}
X_h = X_h^s \oplus X_h^r \ (p_h = p_h^s + p_h^r) \ , \ X_h^r \subset V_h; \\
\forall (p_h^s, p_h^r) \in X_h^s \times X_h^r \ (p_h^s, p_h^r)_H = 0.
\end{array}
\right. 
$$

(H2) “Strong” discrete uniform inf-sup condition:

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

$$ \forall w_h \in M_h, \exists q_h^s \in X_h^s \mid b(w_h, q_h^s) \geq c \|w_h\|_{M/\ker B_h^t} \|q_h^s\|_X. $$

(H3) “Weak” coercivity:

there exists a constant $\alpha > 0$, independent of $h$, such that

$$ \forall p_h \in V_h, a(p_h, p_h) \geq \alpha (\|p_h^s\|_X^2 + |p_h^r|_H^2). $$

(H4) Approximation properties

$$ 
\left\{ 
\begin{array}{ll}
\lim_{h \to 0} \inf_{q_h^s \in X_h^s} \|p - q_h^s\|_X = 0 & \forall p \in X, \\
\lim_{h \to 0} \inf_{w_h \in M_h} \|u - w_h\|_M = 0 & \forall u \in M.
\end{array}
\right. 
$$

REMARK 3.2. Note that, since $X$ is dense in $H$, (H4) implies that

$$ \lim_{h \to 0} \inf_{q_h^s \in X_h^s} |p - q_h^s|_H = 0 \quad \forall p \in H. $$

We can introduce, as for the continuous problem, an operator $B_h$ from $X_h$ to $M_h$ defined by

$$ (B_h p_h, w_h)_{M' \times M} = (p_h, B_h^t w_h)_{X \times X'} = b(p_h, w_h) \quad \forall p_h \in X_h \forall w_h \in M_h $$

and the kernels of $B_h$ (see (3.13)) and of its transpose

$$ \operatorname{Ker} B_h^t = \{ w_h \in M_h \mid b(p_h, w_h) = 0 \ \forall p_h \in X_h \}. $$

REMARK 3.3. It may be more convenient to characterize hypothesis (H0) with one of the following equivalent statements:

$$ (H0ii) \quad \forall p \in X, \exists p_h \in X_h \text{ such that } b(p - p_h, w_h) = 0 \ \forall w_h \in M_h; $$

$$ (H0iii) \quad \operatorname{Ker} B_h^t = \operatorname{Ker} B^t \cap M_h \subset \operatorname{Ker} B^t. $$

REMARK 3.4. We call hypothesis (H2) “strong” the existence of a $q_h^s$ in $X_h^s$ instead of in $X_h$, which would be the classical assumption; more precisely $X_h^s \subsetneq X_h$. On the contrary, hypothesis (H3) is “weaker” in the sense that $\forall p_h \in X_h^s$ or $p_h \in X_h^r$,

$$ \|p_h^s\|_X^2 + |p_h^r|_H^2 \leq \|p_h\|_X^2. $$
THEOREM 3.2. Under hypotheses (H0)-(H4), problem (3.12) admits a unique solution such that

$$ (p_h = p_h^s + p_h^r, u_h) \in X_h \times (M_h/\ker B_h^t), $$

and the following convergence result holds:

$$ \bullet \quad (p_h^s, u_h) \to (p, u) \quad \text{in } X \times M, $$

$$ \bullet \quad p_h^r \to 0 \quad \text{in } H. $$

More precisely,

$$ |p_h^r|_H + \|p - p_h^s\|_X + \|u - u_h\|_{M/\ker B_h^t} \le C \left( \inf_{q_h^s \in X_h^s} \|p - q_h^s\|_X + \inf_{w_h \in M_h} \|u - w_h\|_M + \inf_{z_h^s \in X_h^s} |Ap - z_h^s|_H \right). $$

Proof. First note that hypothesis (H0) and nonuniform discrete coercivity on the kernel $\ker B_h$, i.e.,

$$ (3.14) \qquad \exists \alpha_h > 0 \quad \forall p_h \in V_h \ a(p_h, p_h) \ge \alpha_h \|p_h\|_X^2, $$

ensure the existence and uniqueness of the solution $(p_h, u_h)$ of the discrete problem (3.12) in $X_h \times (M_h/\ker B_h^t)$. 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 \in V_h(f)$ (cf. (3.13)). If we also take $q_h \in V_h(f)$, the difference is in the kernel, i.e., $p_h - q_h \in V_h$. We can write

$$ (3.15) \qquad 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). $$

Taking the difference between the first equation of the continuous (3.11) and the discrete (3.12) problem we have

$$ (3.16) \qquad a(p_h - p, p_h - q_h) = b(u - w_h, p_h - q_h) \quad \forall w_h \in M_h. $$

Using (H1) in (3.16), we can write, with obvious notation,

$$ (3.17) \qquad a(p_h - p, p_h - q_h) = b(u - w_h, p_h^s - q_h^s) + b(u, p_h^r - q_h^r). $$

We have used the fact that $b(w_h, p_h^r - q_h^r) = 0$, since $X_h^r \subset V_h$. Now, since we only have a “weak” coercivity, we want to keep only terms leading to $H$-norm on the $X_h^r$ part; thus we want to eliminate $b(u, p_h^r - q_h^r)$ (which would give $\|p_h^r - q_h^r\|_X$ from the continuity of $b$). Then, by (3.11),

$$ (3.18) \qquad \begin{aligned} a(p_h - p, p_h - q_h) &= b(u - w_h, p_h^s - q_h^s) - a(p, p_h^r - q_h^r) \\ &= b(u - w_h, p_h^s - q_h^s) - (Ap, p_h^r - q_h^r)_H. \end{aligned} $$

We choose now

$$ (3.19) \qquad \left| \begin{array}{l} q_h = q_h^s \in V_h^s(f) = V_h(f) \cap X_h^s, \\ q_h^r = 0. \end{array} \right. $$
From (3.15), (3.18), and (3.19), we obtain

$$a(p_h - q_h^s, p_h - q_h^s) = a(p - q_h^s, p_h - q_h^s) + b(u - w_h, p_h^s - q_h^s) - (\mathcal{A}p, p_h^r),$$

or, by using the orthogonality of $X_h^r$ and $X_h^s$ from (H1),

$$(3.20) \qquad \left| \begin{array}{l} a(p_h - q_h^s, p_h - q_h^s) = a(p - q_h^s, p_h - q_h^s) + b(u - w_h, p_h^s - q_h^s) \\ \qquad \qquad - (\mathcal{A}p - z_h^s, p_h^r)_H \quad \forall z_h^s \in X_h^s. \end{array} \right.$$

Furthermore, from the inequality

$$|p_h - q_h^s|_H \le \|p_h^s - q_h^s\|_X + |p_h^r|_H$$

and from hypothesis (H1) and (H3), (3.20) leads to

$$\left| \begin{array}{l} \alpha(\|p_h^s - q_h^s\|_X^2 + |p_h^r|_H^2) \le \|a\| |p - q_h^s|_H (\|p_h^s - q_h^s\|_X + |p_h^r|_H) \\ \qquad + \|b\| \|u - w_h\|_M \|p_h^s - q_h^s\|_X + |\mathcal{A}p - z_h^s|_H |p_h^r|_H. \end{array} \right.$$

We deduce the existence of a constant $C$ depending on $f$, $p$, $\|a\|$, $\|b\|$, and $\alpha$ such that $\forall (q_h^s, w_h, z_h^s) \in V_h^s(f) \times M_h \times X_h^s,$

$$\|p_h^s - q_h^s\|_X + |p_h^r|_H \le C (|p - q_h^s|_H + \|u - w_h\|_M + |\mathcal{A}p - z_h^s|_H).$$

This gives, using (3.9),

$$\left| \begin{array}{l} \|p - p_h^s\|_X \le \|p - q_h^s\|_X + \|p_h^s - q_h^s\|_X \quad \forall q_h^s \in V_h^s(f) \\ \le (1 + C) \left( \inf_{q_h^s \in V_h^s(f)} \|p - q_h^s\|_X + \inf_{w_h \in M_h} \|u - w_h\|_M + \inf_{z_h^s \in X_h^s} |\mathcal{A}p - z_h^s|_H \right), \\ |p_h^r|_H \le C \left( \inf_{q_h^s \in V_h^s(f)} \|p - q_h^s\|_X + \inf_{w_h \in M_h} \|u - w_h\|_M + \inf_{z_h^s \in X_h^s} |\mathcal{A}p - z_h^s|_H \right). \end{array} \right.$$

To conclude, let us recall that the inf-sup condition (H2) implies (cf. [6]) that

$$\inf_{q_h^s \in V_h^s(f)} \|p - q_h^s\|_X \le c_1 \inf_{q_h^s \in X_h^s} \|p - q_h^s\|_X.$$

Finally, it remains to prove estimates for $\|u - w_h\|_{M/\mathrm{Ker} B_h^t}$. Let us subtract the first equation of (3.12) from the first equation of (3.11). We get

$$a(p - p_h, q_h) + b(u - u_h, q_h) = 0 \quad \forall q_h \in X_h,$$

so that, for any $w_h$ in $M_h$,

$$b(u_h - w_h, q_h) = a(p - p_h, q_h) + b(u - w_h, q_h) \quad \forall q_h \in X_h.$$

Choosing now $q_h = q_h^s$ implies that

$$b(u_h - w_h, q_h^s) = a(p - p_h, q_h^s) + b(u - w_h, q_h^s) \quad \forall q_h^s \in X_h^s.$$
Using this and the inf-sup condition, we have

$$
\begin{align*}
\|u_h - w_h\|_{M/\mathrm{Ker} B_h^t} &\le \frac{1}{c} \sup_{q_h^s \in X_h^s} \frac{b(u_h - w_h, q_h^s)}{\|q_h^s\|_X} \\
&\le \frac{1}{c} \sup_{q_h^s \in X_h^s} \frac{a(p - p_h, q_h^s) + b(u - w_h, q_h^s)}{\|q_h^s\|_X} \\
&\le \frac{1}{c} (\|a\| |p - p_h|_H + \|b\| \|u - w_h\|_M).
\end{align*}
$$

It follows from the triangle inequality that

$$
\|u - u_h\|_{M/\mathrm{Ker} B_h^t} \le C' \left\{ \inf_{w_h \in M_h} \|u - w_h\|_M + |p_h^r|_H + |p - p_h^s|_H \right\}. \quad \Box
$$

REMARK 3.5. One could easily remark that the space $X_h^s$ satisfies all assumptions (H0) to (H4). Why not use $X_h^s$ as approximate space instead of $X_h$ ? In fact, we should imagine the case in which either we cannot characterize the space $X_h^s$ 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.

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

$$
H = (L^2(\Omega))^2; \quad X = H(\mathrm{div}; \Omega); \quad M = L^2(\Omega).
$$

The operator $B$ is the divergence operator, which is surjective from $X$ onto $M$, so that $\mathrm{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_h^{RT}$ onto $M_h$. Since the new space $X_h$ contains the space $X_h^{RT}$, as noted in Remark 2.7, it is straightforward to check (H0); therefore, $\mathrm{Ker} \, B_h^t = \{0\}$.

3.3.1. Approximation with the lowest order finite element $Q_1^{div} - Q_0$.

Let us begin by checking hypothesis (H1). We take as approximation spaces

$$
\left\{
\begin{array}{l}
X_h = \{q_h \in X / \forall K \in \mathcal{T}_h, q_h|_K \in Q_1 \times Q_1\}, \\
M_h = \{w_h \in M / \forall K \in \mathcal{T}_h, w_h|_K \in Q_0\}.
\end{array}
\right.
$$

We define $X_h^s$ as the lowest order Raviart–Thomas element $RT_{[0]}$ defined in (2.10)

$$
X_h^s = \{ q_h \in X / \forall K \in \mathcal{T}_h, q_h|_K \in \mathrm{RT}_{[0]} = P_{10} \times P_{01} \} \subset X_h, \dim \mathrm{RT}_{[0]} = 4.
$$

In order to describe its orthogonal complement $X_h^r$, denote by

- (i, j), the node of $\mathcal{T}_h$ with coordinates $(ih, jh)$,

- $(i + \frac{1}{2}, j)$, the horizontal side joining the nodes $(i, j)$ and $(i + 1, j)$,

- $(i, j + \frac{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)

$$
\left\{
\begin{array}{ll}
\vec{\varphi}_{i,j+\frac{1}{2}} = \begin{pmatrix} \varphi_{i,j+\frac{1}{2}}(x) \\ 0 \end{pmatrix}, & \vec{\varphi}_{i+\frac{1}{2},j} = \begin{pmatrix} 0 \\ \varphi_{i+\frac{1}{2},j}(y) \end{pmatrix}, \\[1em]
\vec{\varphi}_{i+1,j+\frac{1}{2}} = \begin{pmatrix} \varphi_{i+1,j+\frac{1}{2}}(x) \\ 0 \end{pmatrix}, & \vec{\varphi}_{i+\frac{1}{2},j+1} = \begin{pmatrix} 0 \\ \varphi_{i+\frac{1}{2},j+1}(y) \end{pmatrix},
\end{array}
\right.
$$
<figure><img src="image_10.png" /><figcaption>FIG. 3.2. The RT$_{[0]}$ element.</figcaption></figure>

where

$$
\left\{
\begin{array}{l}
\varphi_{i,j+\frac{1}{2}}(x) = \frac{(x_{i+1}-x)}{h} \in P_{10}, \quad \varphi_{i+\frac{1}{2},j}(y) = \frac{(y_{j+1}-y)}{h} \in P_{01}, \\
\\
\varphi_{i+1,j+\frac{1}{2}}(x) = \frac{(x-x_i)}{h} \in P_{10}, \quad \varphi_{i+\frac{1}{2},j+1}(y) = \frac{(y-y_j)}{h} \in P_{01}.
\end{array}
\right.
$$

It is then easy to prove the following lemma.

LEMMA 3.3. If $X_h^r = (X_h^s)^\perp$, the space $X_h^r|_K$ of the restrictions to $K$ of elements of $X_h^r$ can be generated from the four functions:

$$
\vec{\psi}_{i+\frac{1}{2},j} = \begin{pmatrix} 0 \\ (x - x_{i+\frac{1}{2}}) \varphi_{i+\frac{1}{2},j}(y) \end{pmatrix}, \quad \vec{\psi}_{i,j+\frac{1}{2}} = \begin{pmatrix} (y - y_{j+\frac{1}{2}}) \varphi_{i,j+\frac{1}{2}}(x) \\ 0 \end{pmatrix}, \\
\vec{\psi}_{i+\frac{1}{2},j+1} = \begin{pmatrix} 0 \\ (x - x_{i+\frac{1}{2}}) \varphi_{i+\frac{1}{2},j+1}(y) \end{pmatrix}, \quad \vec{\psi}_{i+1,j+\frac{1}{2}} = \begin{pmatrix} (y - y_{j+\frac{1}{2}}) \varphi_{i+1,j+\frac{1}{2}}(x) \\ 0 \end{pmatrix},
$$

where $x_{i+\frac{1}{2}} = (i + \frac{1}{2})h$ and $y_{j+\frac{1}{2}} = (j + \frac{1}{2})h$. Moreover, $X_h^r \subset V_h$.

REMARK 3.6. Note that the dimension of $Q_1 \times Q_1$ (=8) is the sum of the dimensions of $X_h^s|_K$ (= dim $RT_{[0]}$ = 4) and $X_h^r|_K$ (=4).

Using now the well-known properties of $X_h^s$ and $M_h$ [21, 20, 6], hypotheses (H2) and (H4) are trivial; only (H3) needs to be checked.

LEMMA 3.4. For all $p_h$ in $V_h$,

$$
a(p_h, p_h) \geq \alpha \left( \|p_h^s\|_X^2 + \|p_h^r\|_H^2 \right),
$$

$\alpha$ being defined by (2.1).

Proof. From (2.1),

$$
a(p_h, p_h) \geq \alpha \|p_h\|_0^2 = \alpha \left( \|p_h^s\|_0^2 + \|p_h^r\|_0^2 \right).
$$
Note now that $p_h \in V_h$ and $X_h^r \subset V_h$ imply $p_h^s \in V_h \cap X_h^s$. Hence it follows (cf. [6]) that

$$ egin{cases} 
(i) & \mathrm{div} p_h^s \in M_h, \\
(ii) & (\mathrm{div} p_h^s, w_h) = 0 \quad \forall w_h \in M_h.
\end{cases} $$

Recalling (i), we can take

$$ w_h = \mathrm{div} \, p_h^s $$

in (ii), which gives

$$ \mathrm{div} \, p_h^s = 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_h^t = \{0\}$).

**THEOREM 3.5.** *The problem (3.6) admits a unique solution*

$$ (p_h, u_h) \in X_h \times M_h, $$

which satisfies

$$ \bullet \quad (p_h^s, u_h) \quad \longrightarrow \quad (p, u) \quad \text{in } H(\mathrm{div} ; \Omega) \times L^2(\Omega), $$

$$ \bullet \quad p_h^r \quad \longrightarrow \quad 0 \quad \text{in } (L^2(\Omega))^2. $$

Furthermore, if the solution is assumed to be sufficiently regular, i.e., $(p, u) \in (H^1(\Omega))^2 \times H^1(\Omega)$, $\mathrm{div} \, p \in L^2(\Omega)$, and $\mathcal{A}p \in (H^1(\Omega))^2$, then

$$ \|p_h^r\|_{(L^2(\Omega))^2} + \|p - p_h^s\|_{H(\mathrm{div};\Omega)} + \|u - u_h\|_{L^2(\Omega)} $$

$$ \leq C h \left( |u|_{H^1} + |p|_{(H^1)^2} + |\mathrm{div} \, p|_{H^1} + |\mathcal{A}p|_{(H^1)^2} \right), $$

where $|\cdot|_{H^1}$ (resp., $|\cdot|_{(H^1)^2}$) denotes the usual semi-norm in $H^1(\Omega)$ (resp., $(H^1(\Omega))^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

$$ (3.21) \quad \left\{ \begin{array}{ll} \text{Find } (p_h, u_h) \in X_h \times M_h \text{ such that} \\ a_d(p_h, q_h) + b(u_h, q_h) = 0 & \forall q_h \in X_h, \\ b(w_h, p_h) = - \langle f, w_h \rangle & \forall w_h \in M_h \end{array} \right. $$

admits a unique solution $(p_h = p_h^s + p_h^r, u_h)$ with $p_h^r = 0$ and $(p_h^s, u_h)$ being the solution of the following problem:

$$ \left\{ \begin{array}{ll} \text{Find } (p_h^s, u_h) \in X_h^s \times M_h \text{ such that} \\ a_d(p_h^s, q_h^s) + b(u_h, q_h^s) = 0 & \forall q_h^s \in X_h^s, \\ b(w_h, p_h^s) = - \langle f, w_h \rangle & \forall w_h \in M_h \end{array} \right. $$
where $a_d(\cdot, \cdot)$ is given by (2.6) after replacing $A(x)$ by $A_d(x)$. Indeed, in this particular case,

$$ \left\{ \begin{array}{ll} a_d(p_h^r, q_h^s) = 0 & \forall q_h^s \in X_h^s, \\ a_d(p_h^s, q_h^r) = 0 & \forall q_h^r \in X_h^r. \end{array} \right. $$

By using the fact that $X_h^r \subset V_h$, we can decompose (3.21) into two independent problems in $X_h^s$ and $X_h^r$:

$$ \left\{ \begin{array}{l} \text{Find } (p_h^s, u_h) \in X_h^s \times M_h \text{ such that} \\ \qquad a_d(p_h^s, q_h^s) + b(u_h, q_h^s) = 0 \quad \forall q_h^s \in X_h^s, \\ \qquad b(w_h, p_h^s) = -\langle f, w_h \rangle \quad \forall w_h \in M_h, \end{array} \right. $$

and

$$ (3.22) \qquad \left\{ \begin{array}{l} \text{Find } p_h^r \in X_h^r \text{ such that} \\ \qquad a_d(p_h^r, q_h^r) = 0 \quad \forall q_h^r \in X_h^r. \end{array} \right. $$

It is obvious then from (3.22) that $p_h^r \equiv 0$. This remark is no longer true when $A(x)$ is not diagonal.

### 3.3.2. Approximation with higher order finite elements, $Q_{k+1}^{\text{div}} - Q_k$.
We have seen in section 2.3 that the natural generalization of the lowest order element consists of taking

$$ \hat{X} = X_k = Q_{k+1} \times Q_{k+1} \quad \text{and} \quad \hat{M} = M_k = Q_k; $$

and thus to introduce the spaces

$$ (2.16) \qquad \left\{ \begin{array}{l} X_h = \{q_h \in X / \forall K \in \mathcal{T}_h, q_h|_K \in X_k\}, \\ M_h = \{w_h \in M / \forall K \in \mathcal{T}_h, w_h|_K \in M_k\}. \end{array} \right. $$

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 $\hat{X}$, considered as a subspace of $L^2(K)$, where $K$ is a single element. First, we recall that

$$ RT_{[k]} = P_{k+1,k} \times P_{k,k+1}, $$

which obviously satisfies the inclusion

$$ RT_{[k]} \subset X_k. $$

Let us define $\Psi_k$ as the orthogonal complement in $X_k$ of $RT_{[k]}$ (for the inner product in $L^2(K)$):

$$ \Psi_k = \left\{ \psi \in \hat{X} / \int_K \psi \, \phi \, dx = 0 \quad \forall \phi \in RT_{[k]} \right\}. $$
Note that

$$ \dim \Psi_k = 2(k + 2). $$

The main property of the space $\Psi_k$ is the following.

LEMMA 3.6. For any $\psi$ in $\Psi_k$ and $v$ in $M_k$,

$$ (3.23) \qquad \int_K \operatorname{div} \psi v \,dx = 0. $$

*Proof.* For simplicity we consider the reference element $K = [0, 1] \times [0, 1]$. (i) We begin by a characterization of $\Psi_k$. Let $\sigma_k$ be the polynomial in one variable of degree $k+1$ such that

$$ P_{k+1} = P_k \oplus [\sigma_k]. $$

Equivalently, $\sigma_k$ is defined, up to a multiplicative constant, by

$$ \left\{ \begin{array}{l} \displaystyle \int_0^1 \sigma_k(x)p(x)dx = 0 \quad \forall p \in P_k, \\ \sigma_k \in P_{k+1}, \sigma_k \neq 0. \end{array} \right. $$

Then, we claim that

$$ \Psi_k = \tilde{\Psi}_k \equiv \left\{ \begin{pmatrix} p_x(x)\sigma_k(y) \\ p_y(y)\sigma_k(x) \end{pmatrix}, (p_x, p_y) \in P_{k+1} \times P_{k+1} \right\}. $$

Indeed, as $\dim \tilde{\Psi}_k = 2(k+2)$, it is sufficient to check that $\tilde{\Psi}_k$ is orthogonal to $RT_{[k]}$. To see that, we first remark that (see [21]) $RT_{[k]}$ is generated by vector fields of the form

$$ \phi(x, y) = \begin{bmatrix} \tau_x(x)q_x(y) \\ \tau_y(y)q_y(x) \end{bmatrix} \quad \text{with} \quad (\tau_x, \tau_y) \in P_{k+1} \times P_{k+1}, (q_x, q_y) \in P_k \times P_k. $$

Let us consider

$$ (3.24) \qquad \psi(x, y) = \begin{bmatrix} p_x(x)\sigma_k(y) \\ p_y(y)\sigma_k(x) \end{bmatrix} \quad \text{with} \quad (p_x, p_y) \in P_{k+1} \times P_{k+1}. $$

We have

$$ (\phi, \psi)_{L^2(K)} = \left( \int_0^1 \tau_x(x) p_x(x) dx \right) \left( \int_0^1 \sigma_k(y) q_x(y) dy \right) \\ + \left( \int_0^1 \sigma_k(x) q_y(x) dx \right) \left( \int_0^1 \tau_y(y) p_y(y) dy \right). $$

As $q_x$ and $q_y$ are in $P_k$, by the definition of $\sigma_k$,

$$ \int_0^1 \sigma_k(y) q_x(y) dy = \int_0^1 \sigma_k(x) q_y(x) dx = 0 \Rightarrow (\phi, \psi)_{L^2(K)} = 0. $$
<figure><img src="image_11.png" /><figcaption>FIG. 3.3. The edges of K.</figcaption></figure>

(ii) To prove (3.23) we first use Green's formula

$$ \int_K \operatorname{div} \psi v dx = - \int_K \psi \cdot \nabla v dx + \int_{\partial K} (\psi \cdot n) v d\gamma. $$

If $v$ belongs to $Q_{k+1}$, $\nabla v$ belongs to $RT_{[k]}$. Therefore, since $\psi$ belongs to $\Psi_k$,

$$ \int_K \operatorname{div} \psi v \, dx = \int_{\partial K} (\psi \cdot n) v \, d\gamma. $$

Let us decompose $\partial K$ as

$$ \partial K = T_1 \cup T_2 \cup T_3 \cup T_4, $$

as shown in Figure 3.3. Then, we check that

$$ \text{for } 1 \le j \le 4 \text{ and } \forall q \in R_k(T_j), \quad \int_{T_j} (\psi \cdot n) q d\gamma = 0, $$

where $R_k(T_j)$ is the set of polynomials of degree $k$ with respect to the abscissa along $T_j$. Let us consider, for instance, $j = 1$ and assume that $\psi$ is given by (3.24). Then

$$ \int_{T_1} (\psi \cdot n) q d\gamma = - \left( \int_0^1 \sigma_k(x) q(x) \right) p_y(0) = 0. $$

To conclude, it remains to notice that, if $v$ belongs to $Q_{k+1}$, then $v|_{T_j}$ belongs to $R_k(T_j)$. $\square$

This result suggests that we define $X_h^s$ and $X_h^r$ by

$$ (3.25) \qquad \begin{cases} X_h^s = \{p_h \in H(\operatorname{div}; \Omega) / \forall K \in \mathcal{T}_h, p_h|_K \in RT_{[k]}\}, \\ X_h^r = \{p_h \in H(\operatorname{div}; \Omega) / \forall K \in \mathcal{T}_h, p_h|_K \in \Psi_k\}. \end{cases} $$

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_h^s$ and $X_h^r$ 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 $\Omega$ as the sum of the integrals over elements K) of the definition of $\Psi_k$ and Lemma 3.6. $\square$
<figure><img src="image_12.png" /><figcaption>FIG. 3.4. Orthogonal decomposition of $Q_1 \times Q_1$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 3.5. Orthogonal decomposition of $Q_{k+1}$.</figcaption></figure>

### 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_h^s$ and $M_h$ defined by (2.16), (3.25) we obtain the following result.

### THEOREM 3.8.
Problem (3.6) admits a unique solution $(p_h, u_h) \in X_h \times M_h$, which satisfies

$$ \bullet \quad (p_h^s, u_h) \to (p, u) \quad \text{in } H(\mathrm{div}; \Omega) \times L^2(\Omega), $$

$$ \bullet \quad p_h^r \to 0 \quad \text{in } (L^2(\Omega))^2. $$

Furthermore, if we assume that the solution $(p, u)$ of (3.5) is such that $(p, u) \in (H^m(\Omega))^2 \times (H^m(\Omega))$, $\mathrm{div}\,p \in H^m(\Omega)$, and $A p \in (H^m(\Omega))^2$ for $m = k + 1$, then

$$ (3.26) \quad \left\{ \begin{array}{l} \|p_h^r\|_{(L^2)^2} + \|p - p_h^s\|_{H(\mathrm{div};\Omega)} + \|u - u_h\|_{L^2} \\ \le C \, h^m \left( |u|_{H^m} + |p|_{(H^m)^2} + |\mathrm{div}\,p|_{H^m} + |Ap|_{(H^m)^2} \right), \end{array} \right. $$
where $|\cdot|_{H^m}$ (resp., $|\cdot|_{(H^m)^2}$) denotes the usual seminorm in $H^m(\Omega)$ (resp., $(H^m(\Omega))^2$).

**3.4. Application to the elliptic projection operator.** We come back to the abstract framework described in section 3.2. We consider the problem of finding $(\hat{p}_h, \hat{u}_h) \equiv \Pi_h(p, u)$ such that

$$ (3.27) \qquad \left\{ \begin{array}{ll} a(p - \hat{p}_h, q_h) + b(u - \hat{u}_h, q_h) & = 0 \quad \forall q_h \in X_h, \\ b(w_h, p - \hat{p}_h) & = 0 \quad \forall w_h \in M_h. \end{array} \right. $$

We set

$$ D(B^t) = \{w \in M / b(w, q) \le C(w) |q|_H \quad \forall q \in X\} = \{w \in M / B^t w \in H\}. $$

For $w \in D(B^t)$, we have

$$ b(w, q) = (B^t w, q)_H \quad \forall q \in X. $$

Let us introduce the notation

$$ (3.28) \qquad \left\{ \begin{array}{l} |||p - \hat{p}_h||| = \|p - \hat{p}_h^s\|_X + |\hat{p}_h^r|_H, \\ |||(p, u) - \Pi_h(p, u)|||_C = |||p - \hat{p}_h||| + \|u - \hat{u}_h\|_M. \end{array} \right. $$

It is straightforward, by application of the abstract result, to get interpolation estimates.

**THEOREM 3.9.** Assume (H0)-(H4).

(i) For all $(p, u) \in X \times M$, problem (3.27) admits a unique solution

$$ \Pi_h(p, u) = (\hat{p}_h^s + \hat{p}_h^r, \hat{u}_h) \in X_h \times (M_h / \text{Ker } B_h^t). $$

(ii) For all $(p, u) \in X \times D(B^t)$, there exists a constant $C$, independent of $h$, such that

$$ |||(p, u) - \Pi_h(p, u)|||_C \le C \mathcal{E}_h(p, u), $$

where

$$ \mathcal{E}_h(p, u) = \inf_{q_h^s \in X_h^s} \|p - q_h^s\|_X + \inf_{w_h \in M_h} \|u - w_h\|_M + \inf_{z_h^s \in X_h^s} |B^t u - z_h^s|_H. $$

In particular, $|||p - \hat{p}_h|||$ and $\|u - \hat{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 $\text{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_h^r - q_h^r)$ 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_h^r - q_h^r) = -a(p, p_h^r - q_h^r) \equiv -(Ap, p_h^r - q_h^r)_H$. Now we assume that $u \in D(B^t)$, so that we have

$$ b(u, p_h^r - q_h^r) = (B^t u, p_h^r - q_h^r)_H. $$

The result is thus obtained by replacing $-Ap$ with $B^t u$. $\square$

<figure><img src="image_14.png" /></figure>
**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

$$ (4.1) \qquad \left\{ \begin{array}{ll} \text{Find } (p, v) : [0, T] \to X \times M \text{ such that} \\ \frac{d}{dt}a(p, q) + b(v, q) & = 0 \quad \forall q \in X, \\ \frac{d}{dt}(v, w) - b(w, p) & = (f, w) \quad \forall w \in M, \\ p(0) = p_0; \quad v(0) = v_0. \end{array} \right. $$

Or, equivalently, in an operator form, if we assume enough regularity on the solution in time,

$$ (4.2) \qquad \left\{ \begin{array}{ll} \text{Find } (p, v) \in (C^1(0, T; H) \cap C^0(0, T; X)) \times C^1(0, T; M) \text{ such that} \\ \mathcal{A} \frac{dp}{dt} + B^t v & = 0 \quad \text{in } X', \\ \frac{dv}{dt} - Bp & = f \quad \text{in } M', \\ p(0) = p_0; \quad v(0) = v_0. \end{array} \right. $$

In the following, we use the notation $C^{m,r} = C^m(0, T; H) \cap 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

$$ (4.3) \qquad \left\{ \begin{array}{ll} \text{Find } (p_h, v_h) : [0, T] \to X_h \times M_h \text{ such that} \\ \frac{d}{dt}a(p_h, q_h) + b(v_h, q_h) & = 0 \quad \forall q_h \in X_h, \\ \frac{d}{dt}(v_h, w_h) - b(w_h, p_h) & = (f_h, w_h) \quad \forall w_h \in M_h, \\ p_h(0) = p_{0,h}; \quad v_h(0) = v_{0,h}. \end{array} \right. $$

From the classical theory of ODEs, we have the following result.

**THEOREM 4.1.** If $f \in C^0(0, T; M_h)$, then problem (4.3) has a unique solution $(p_h, v_h) \in C^1(0, T; X_h) \times C^1(0, T; M_h)$.

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) \in C^{1,0} \times C^1(0, T; M)$. Then we have the following:

(i) There exists a primitive of $v$, $u \in C^1(0, T; M)$, satisfying

$$ (4.4) \qquad \left\{ \begin{array}{l} \displaystyle \frac{du}{dt} = v, \\ \displaystyle A p_0 + B^t u(0) = 0. \end{array} \right. $$

This primitive is unique up to a constant element of Ker $B^t$.

(ii) $\forall t \in [0, T]$, problem (3.27) admits a unique solution $\Pi_h(p, u)(t) = (\hat{p}_h, \hat{u}_h)(t) \in X_h \times (M_h/\text{Ker } B_h^t)$ and there exists a constant $C$ independent of $h$ such that

$$ (4.5) \qquad |||(p, u) - \Pi_h(p, u)|||_C(t) \le C \mathcal{E}_h(p, u)(t), $$

where $\mathcal{E}_h$ is defined by

$$ \mathcal{E}_h(p, u)(t) = \inf_{q_h^s \in X_h^s} \|p(t) - q_h^s\|_X + \inf_{w_h \in M_h} \|u(t) - w_h\|_M + \inf_{z_h^s \in X_h^s} |Ap(t) - z_h^s|_H. $$

In particular, $|||p - \hat{p}_h|||$ and $\|u - \hat{u}_h\|_M$ tend to zero uniformly in time ($t \in [0, T]$) ($|||\cdot|||_C$ and $|||\cdot|||$ are defined in (3.28)).

(iii) In the same way, if $(p, u) \in C^k(0, T; X) \times C^k(0, T; M)$, $k \ge 1$, there exists a constant $C$ independent of $h$ such that

$$ (4.6) \qquad |||(\partial_t^k p, \partial_t^k u) - \Pi_h (\partial_t^k p, \partial_t^k u)||| (t) \le C \mathcal{E}_h (\partial_t^k p, \partial_t^k u) (t), $$

where $\partial_t^k g = \frac{d^k g}{dt^k}$.

REMARK 4.1. Operators $\Pi_h$ and $\partial_t^k$ commute, and we set

$$ (4.7) \qquad \hat{v}_h = \partial_t (\hat{u}_h) = (\widehat{\partial_t u})_h. $$

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 \in C^k$ while $v \in 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_h^r$, we used the first equation satisfied by the solution, which relates the embarrassing term $b(v, p_h^r - q_h^r)$ 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 $\frac{d}{dt}a(p, p_h^r - q_h^r)$ (see proof of Theorem 3.2).

Proof (of Lemma 4.2). (i) We set $f_0 = -Bp_0 \in \text{Im } B$. From hypothesis (3.10), we know that there is a unique $(p_0, u_0) \in X \times (M/\text{Ker } B^t)$ such that

$$ \left\{ \begin{array}{ll} a(p_0, q) + b(u_0, q) & = 0 \quad \forall q \in X, \\ b(w, p_0) & = -(f_0, w) \quad \forall w \in M, \end{array} \right. $$

which means that, $p_0$ being fixed, there is a unique $u_0 \in (M/\text{Ker } B^t)$ such that $Ap_0 + B^t u_0 = 0$. Now we define $u$ as

$$ u(t) = u_0 + \int_0^t v(s)ds. $$
It is clear that $u \in C^1(0, T; M)$ and is the unique solution of (4.4).

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

$$ \frac{d}{dt}(Ap + B^t u) = 0 \implies (Ap + B^t u)(t) = Ap_0 + B^t u(0) = 0; $$

thus, $(p, u) \in C^1(0, T; X) \times C^1(0, T; M)$ satisfies

$$ a(p, q) + b(u, q) = 0 \quad \forall q \in X, $$

so that $u \in D(B^t)$ and $B^t u = -Ap$. Applying Theorem 3.9, we get the existence and uniqueness of the elliptic projection, for $t$ fixed, and also the error estimate (4.5).

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

$$ \mathcal{A}\partial_t^k p + B^t \partial_t^k u = 0. $$

We again apply Theorem 3.9 to get (4.6). □

We now give the main result.

**THEOREM 4.3.** Assume (H0)–(H4). Let $(p, v)$ be the solution of (4.2) and $(p_h, v_h)$ the solution of the approximate problem (4.3) with the initial conditions

$$ (4.8) \qquad (p_{0,h}, v_{0,h}) = \Pi_h(p_0, v_0). $$

• If $(p, v) \in C^2(0, T; X) \times C^1(0, T; M)$, then, $\forall t \in [0, T]$,

$$ |p - p_h^s|_H(t) \to 0; \quad |p_h^r|_H(t) \to 0; \quad \|v - v_h\|_M(t) \to 0. $$

More precisely,

$$ (4.9) \qquad \begin{cases} |p - p_h^s|_H(t) + |p_h^r|_H(t) \le C \left( \mathcal{E}_h(p, u)(t) + \int_0^t \mathcal{E}_h(\partial_t^2 p, \partial_t v)(s) ds \right), \\ \|v - v_h\|_M(t) \le C \left( \mathcal{E}_h(\partial_t p, v)(t) + \int_0^t \mathcal{E}_h(\partial_t^2 p, \partial_t v)(s) ds \right). \end{cases} $$

• If, in addition, $(p, v) \in C^3(0, T; X) \times C^2(0, T; M)$ and $(p_h, v_h) \in C^2(0, T; X_h) \times C^2(0, T; M_h)$,

$$ \forall t \in [0, T], \quad \|p - p_h^s\|_X(t) \to 0. $$

More precisely,

$$ (4.10) \qquad \left\{ \begin{array}{l} \|p - p_h^s\|_X(t) \le C(\mathcal{E}_h(p, u)(t) + \mathcal{E}_h(\partial_t^2 p, \partial_t v)(t)) \\ \qquad + \int_0^t (\mathcal{E}_h(\partial_t^2 p, \partial_t v)(s) + \mathcal{E}_h(\partial_t^3 p, \partial_t^2 v)(s)) ds). \end{array} \right. $$

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 $\Pi_h(p, u) = (\hat{p}_h^s + \hat{p}_h^r, \hat{u}_h)$ be the elliptic projection defined in Lemma 4.2, and let $\hat{v}_h$ be defined by (4.7). Set $\delta_h = v - \hat{v}_h$.
(i) If $(p, v) \in C^{1,0} \times C^1(0, T; M)$, there exists a constant $C_1$, independent of $h$, such that, $\forall t \in [0, T]$,

$$ (4.11) \quad |\hat{p}_h^s - p_h^s|_H(t) + |\hat{p}_h^r - p_h^r|_H(t) + \|\hat{v}_h - v_h\|_M(t) \le C_1 \int_0^t \|\partial_t \delta_h\|_M(s) \, ds. $$

(ii) Moreover, if $(p, v) \in C^{2,1} \times C^2(0, T; M)$ and $(p_h, v_h) \in C^2(0, T; X_h) \times C^2(0, T; M_h)$, there exists a constant $C_2$, independent of $h$ such that, $\forall t \in [0, T]$,

$$ (4.12) \quad \|\hat{p}_h^s - p_h^s\|_X(t) \le C_2 \left\{ \|\partial_t \delta_h\|_M(t) + \int_0^t \left( \|\partial_t \delta_h\|_M(s) + \|\partial_t^2 \delta_h\|_M(s) \right) ds \right\}. $$

*Proof (of Lemma 4.4).* (i) To prove estimation (4.11), we begin by rewriting (4.1) with test functions $q = q_h \in X_h \subset X$ and $w = w_h \in M_h \subset M$ and by subtracting it from (4.3):

$$ (4.13) \quad \left\{ \begin{array}{ll} \displaystyle \frac{d}{dt}a(p - p_h, q_h) + b(v - v_h, q_h) = 0 & \forall q_h \in X_h, \\ \\ \displaystyle \frac{d}{dt}(v - v_h, w_h) - b(w_h, p - p_h) = 0 & \forall w_h \in M_h, \\ \\ (p - p_h)(0) = p_0 - p_{0,h}; & (v - v_h)(0) = v_0 - v_{0,h}. \end{array} \right. $$

Introducing the elliptic projection $\Pi_h(p, u) = (\hat{p}_h, \hat{u}_h)$, we split the error between the approximate solution and the exact solution into two parts,

$$ (4.14) \quad \left\{ \begin{array}{l} (p - p_h)(t) = (p - \hat{p}_h)(t) + (\hat{p}_h - p_h)(t), \\ (v - v_h)(t) = (v - \hat{v}_h)(t) + (\hat{v}_h - v_h)(t), \end{array} \right. $$

and we choose as approximate initial conditions the elliptic projection of the exact initial condition, (4.8), so that at time $t = 0$

$$ (\hat{p}_h - p_h)(0) = 0 \quad ; \quad \hat{v}_h - v_h(0) = 0. $$

By (4.13) and (4.14), for any $(q_h, w_h) \in X_h \times M_h$,

$$ (4.15) \quad \left\{ \begin{array}{l} a(\partial_t(\hat{p}_h - p_h), q_h) + b(\hat{v}_h - v_h, q_h) = -a(\partial_t(p - \hat{p}_h), q_h) - b(v - \hat{v}_h, q_h), \\ (\partial_t(\hat{v}_h - v_h), w_h) - b(w_h, \hat{p}_h - p_h) = -(\partial_t(v - \hat{v}_h), w_h) + b(w_h, p - \hat{p}_h). \end{array} \right. $$

By differentiating the first equation of (3.27) (written for $(p, u)$) with respect to $t$, we see that

$$ \left\{ \begin{array}{l} a(\partial_t(p - \hat{p}_h), q_h) + b(v - \hat{v}_h, q_h) = 0 \quad \forall q_h \in X_h, \\ \\ b(w_h, p - \hat{p}_h) = 0 \quad \forall w_h \in M_h. \end{array} \right. $$

Substituting into (4.15) gives, for any $(q_h, w_h) \in X_h \times M_h$,

$$ (4.16) \quad \left\{ \begin{array}{l} a(\partial_t(\hat{p}_h - p_h), q_h) + b(\hat{v}_h - v_h, q_h) = 0, \\ (\partial_t(\hat{v}_h - v_h), w_h) - b(w_h, \hat{p}_h - p_h) = -(\partial_t(v - \hat{v}_h), w_h). \end{array} \right. $$
Furthermore, by taking $q_h = \hat{p}_h - p_h$ and $w_h = \hat{v}_h - v_h$ in (4.15) and by adding the two equations, we get

$$ (4.17) \quad a(\partial_t(\hat{p}_h - p_h), \hat{p}_h - p_h) + (\partial_t(\hat{v}_h - v_h), \hat{v}_h - v_h) = -(\partial_t(v - \hat{v}_h), \hat{v}_h - v_h). $$

Next, set

$$ E_h(t) = \frac{1}{2} (a(\hat{p}_h - p_h, \hat{p}_h - p_h) + (\hat{v}_h - v_h, \hat{v}_h - v_h)) (t). $$

Since for some constant $C > 0$,

$$ E_h^{1/2}(t) \geq C \left( |\hat{p}_h - p_h|^2_H(t) + \|\hat{v}_h - v_h\|_M^2(t) \right)^{1/2}, $$

it follows that

$$ (4.18) \quad \frac{dE_h^{1/2}}{dt}(t) \leq C' \|\partial_t(v - \hat{v}_h)\|_M(t) \equiv C' \|\partial_t \delta_h\|_M(t), $$

where we set $\delta_h = v - \hat{v}_h$. From the choice (4.8) for the initial conditions, we deduce

$$ E_h(0) = 0. $$

It is then easy to see that (4.18) leads to (4.11) (from the orthogonality $X_h^s \perp X_h^r$, $|q_h|^2_H = |q_h^s|^2_H + |q_h^r|^2_H \forall q_h \in X_h$).

(ii) To get the estimate (4.12) in the $X$-norm, we recall that, $\forall \eta_h^s \in X_h^s$, we have $\eta_h^s = \eta_1 + \eta_2$ with $\eta_1 \in \text{Ker } B_h$ and $\eta_2 \in (\text{Ker } B_h)^\perp$, so that

$$ \|\eta_h^s\|_X^2 = |\eta_1|^2_H + \|\eta_2\|_X^2 $$

(recalling that $\|\eta_1\|_X = |\eta_1|_H \quad \forall \eta_1 \in \text{Ker } B_h$). We set $\eta_h = \hat{p}_h - p_h$. The term $|\eta_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 $\|\eta_2\|_X$. To do so, we start by recalling that the inf-sup condition (hypothesis (H2)) is equivalent to

$$ (4.19) \quad \left\{ \begin{array}{ll} \text{there exists a constant } C > 0, \text{ independent of } h, \text{ such that} \\ \forall q_h^s \in X_h^s, \quad \sup_{w_h \in M_h} \frac{b(w_h, q_h^s)}{\|w_h\|_M} \geq C \|q_h^s\|_{X/\text{Ker } B_h}. \end{array} \right. $$

We also know that $\|\eta_2\|_X = \|\eta_h^s\|_{X/\text{Ker } B_h}$. Thus, by taking $q_h^s = \eta_2$ in (4.19),

$$ \sup_{w_h \in M_h} \frac{b(w_h, \eta_2)}{\|w_h\|_M} \geq C \|\eta_2\|_X. $$

Now use the second equation of (4.16) to obtain

$$ (4.20) \quad \|\eta_2\|_X \leq C' \{\|\partial_t(v - \hat{v}_h)\|_M + \|\partial_t(\hat{v}_h - v_h)\|_M\}. $$

Till now, we have used only the $C^1$ regularity of the solution. In order to bound $\|\partial_t(\hat{v}_h - v_h)\|_M$, we need $C^2$. Indeed, we want to apply (4.11) with $v_h$ replaced by $\partial_t v_h$, $\hat{v}_h$ by $\partial_t \hat{v}_h$, and so on. More precisely, we have

$$ (4.21) \quad \|\partial_t(\hat{v}_h - v_h)\|_M(t) \leq C \int_0^t \|\partial_t^2 \delta_h\|_M(s) ds. $$
Finally, combining (4.20) and (4.21), we get

$$ \|\eta_2\|_{X}(t) \leq C' \left\{ \|\partial_t \delta_h\|_{M}(t) + \int_{0}^{t} \|\partial_t^2 \delta_h\|_{M}(s) \, ds \right\}, $$

and the proof is achieved. □

*Proof (of Theorem 4.3).* We combine results given in Lemma 4.2 and in Lemma 4.4.

• Estimates in H- and M-norms. We have

$$ (4.22) \qquad |p - p_h^s|_H + |p_h^r|_H \leq |p - \hat{p}_h^s|_H + |\hat{p}_h^r|_H + |\hat{p}_h^s - p_h^s|_H + |\hat{p}_h^r - p_h^r|_H. $$

The first term on the right is bounded by $|||p - \hat{p}_h|||$ and thus by $|||(p, u) - \Pi_h(p, u)|||_C$. Assuming that $(p, v) \in C^{1,0} \times C^1(0, T; M)$, by (4.5) it can be estimated that

$$ (4.23) \qquad |||p - \hat{p}_h||| (t) \leq C \mathcal{E}_h(p, u)(t). $$

The second term is estimated using (4.11) of Lemma 4.4. This requires $\|\partial_t \delta_h\|_M$ to be estimated. For this, we use (4.6), for $k=2$, which requires $(p, v) \in C^2(0, T; X) \times C^1(0, T; M)$. We get

$$ (4.24) \qquad \|\partial_t(v - \hat{v}_h)\|_M(s) \leq C \mathcal{E}_h(\partial_t^2 p, \partial_t v)(s); $$

then (4.22), (4.23), and (4.24) lead to the first inequality of (4.9). Now, for $v$, we write

$$ \|v - v_h\|_M \leq \|v - \hat{v}_h\|_M + \|\hat{v}_h - v_h\|_M. $$

We apply (4.6) for $k=1$ and get

$$ \|v - \hat{v}_h\|_M \leq C \mathcal{E}_h(\partial_t p, v). $$

Again, using estimate (4.11) to bound $\|\hat{v}_h - v_h\|_M$, we easily get the second inequality of (4.9).

• Estimate in X-norm. We write

$$ \|p - p_h^s\|_X \leq \|p - \hat{p}_h^s\|_X + \|\hat{p}_h^s - p_h^s\|_X. $$

The first term is again bounded by $|||p - \hat{p}_h|||$ and can be estimated from (4.5). For the second term, we now use (4.12) of Lemma 4.4, which requires $(p, v) \in C^{2,1} \times C^2(0, T; M)$ and $(p_h, v_h) \in C^2(0, T; X_h) \times C^2(0, T; M_h)$. In the right-hand side of (4.12) appears the second derivative of $v - \hat{v}_h$. We thus use estimate (4.6) for $k=3$, which requires $(p, v) \in C^3(0, T; X) \times C^2(0, T; M)$, and we get (4.10). □

**4.2. Application to the approximation of the anisotropic wave equation with the new finite element, $Q_{k+1}^{div} - 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_{k+1}^{div} - Q_k$ element is covered by the abstract framework. It is then straightforward to apply Theorem 4.3 and the following theorem results.

**THEOREM 4.5.** Let $(p, v)$ be the solution of (2.5), (2.4), and let $u \in C^2(0, T; M)$ be the primitive of $v$, defined by

$$ \frac{du}{dt} = v; \quad B^t u(0) = -A p_0. $$
Let $(p_h, u_h)$ be the solution of the approximate problem (2.9) with initial conditions

$$ (p_{0,h}, v_{0,h}) = \Pi_h(p_0, v_0). $$

• Convergence in the $L^2$-norm:

(i) If $(p, v) \in C^2(0, T; X) \times C^1(0, T; M)$, $\forall t \in [0, T]$,

$$ |p - p_h^s|_H(t) \to 0; \quad |p_h^r|_H(t) \to 0; \quad \|v - v_h\|_M(t) \to 0. $$

(ii) Furthermore, if $(p, u) \in C^2(0, T; H^{k+1}(\text{div}; \Omega)) \times C^2(0, T; H^{k+1}(\Omega))$ and $\mathcal{A}p \in C^2(0, t; (H^{k+1}(\Omega))^2)$, then

$$ |p - p_h^s|_H(t) + |p_h^r|_H(t) + \|v - v_h\|_M(t) \le C_1(t)h^k, $$

with $C_1(t) = O(\|p\|_{C^2(0,t;H^{k+1}(\text{div};\Omega))} + \|\mathcal{A}p\|_{C^2(0,t;(H^{k+1}(\Omega))^2)} + \|u\|_{C^2(0,t;H^{k+1}(\Omega))})$.

• Convergence in the $H(\text{div})$-norm:

(iii) If $(p, v) \in C^3(0, T; X) \times C^2(0, T; M)$, $\forall t \in [0, T]$,

$$ \|p - p_h^s\|_X(t) \to 0. $$

(iv) Furthermore, if $(p, u) \in C^3(0, T; H^{k+1}(\text{div}; \Omega)) \times C^3(0, T; H^{k+1}(\Omega))$ and $\mathcal{A}p \in C^3(0, t; (H^{k+1}(\Omega))^2)$, then

$$ \|p - p_h^s\|_X(t) \le C_2(t)h^k, $$

with $C_2(t) = O(\|p\|_{C^3(0,t;H^{k+1}(\text{div};\Omega))} + \|\mathcal{A}p\|_{C^3(0,t;(H^{k+1}(\Omega))^2)} + \|u\|_{C^3(0,t;H^{k+1}(\Omega))})$.

*Proof.* We apply Theorem 4.3 which relates the errors to quantities such as $\mathcal{E}_h(\partial_t^m p, \partial_t^l u)$, i.e., to the error due to the approximation of $H(\text{div}; \Omega)$ with the Raviart-Thomas $RT_{[k]}$ space and to the error due to the approximation of $L^2(\Omega)$ with the $Q_k$ discontinuous elements. □

**Appendix. Gauss–Lobatto and Gauss–Legendre quadrature formulae in two dimensions.** We give the locations of the quadrature points and the associated weights used respectively for the $Q_2^{\text{div}} - Q_1$ element ($k=1$, Figure 2.5) and the $Q_3^{\text{div}} - Q_2$ element ($k=2$, Figure 2.6).

For $k=1$ (Figure 2.5):

• The Gauss–Lobatto points and weights:

$$ \left\{ \begin{array}{l} x_1^p = y_1^p = 0, \quad x_2^p = y_2^p = \frac{1}{2}, \quad x_3^p = y_3^p = 1; \\ \\ \omega_{11}^p = \omega_{31}^p = \omega_{33}^p = \omega_{13}^p = \frac{1}{36}; \\ \\ \omega_{21}^p = \omega_{32}^p = \omega_{23}^p = \omega_{12}^p = \frac{1}{9}; \\ \\ \omega_{22}^p = \frac{4}{9}. \end{array} \right. $$

• The Gauss–Legendre points and weights:

$$ \left\{ \begin{array}{l} x_1^v = y_1^v = \frac{3 - \sqrt{3}}{6}, \quad x_2^v = y_2^v = \frac{3 + \sqrt{3}}{6}, \\ \\ \omega_{11}^v = \omega_{21}^v = \omega_{22}^v = \omega_{12}^v = \frac{1}{4}. \end{array} \right. $$
For $k = 2$ (Figure 2.6):

* The Gauss-Lobatto points and weights:

$$
\left\{
\begin{array}{l}
x_1^p = y_1^p = 0, \quad x_2^p = y_2^p = \frac{5-\sqrt{5}}{10}, \quad x_3^p = y_3^p = \frac{5+\sqrt{5}}{10}, \quad x_4^p = y_4^p = 1; \\
\\
\omega_{11}^p = \omega_{41}^p = \omega_{44}^p = \omega_{14}^p = \frac{1}{144}; \\
\\
\omega_{21}^p = \omega_{31}^p = \omega_{42}^p = \omega_{43}^p = \omega_{34}^p = \omega_{24p} = \omega_{13}^p = \omega_{12}^p = \frac{5}{144}; \\
\\
\omega_{22}^p = \omega_{32}^p = \omega_{33}^p = \omega_{23}^p = \frac{25}{144}.
\end{array}
\right.
$$

* The Gauss-Legendre points and weights:

$$
\left\{
\begin{array}{l}
x_1^v = y_1^v = \frac{5 - \sqrt{15}}{10}, \quad x_2^v = y_2^v = \frac{1}{2}, \quad x_3^v = y_3^v = \frac{5 + \sqrt{15}}{10}; \\
\\
\omega_{11}^v = \omega_{31}^v = \omega_{33}^v = \omega_{13}^v = \frac{25}{324}; \\
\\
\omega_{21}^v = \omega_{32}^v = \omega_{23}^v = \omega_{12}^v = \frac{10}{81}; \\
\\
\omega_{22}^v = \frac{16}{81}.
\end{array}
\right.
$$

REFERENCES

[1] I. BABUSKA, *The finite element method with Lagrangian multipliers*, Numer. Math., 20 (1973), pp. 179–192.

[2] E. BÉCACHE, P. JOLY, AND C. TSOGKA, *Eléments finis mixtes et condensation de masse en élastodynamique linéaire*. (I) *Construction*, C. R. Acad. Sci. Paris Sér. I Math., 325 (1997), pp. 545–550.

[3] E. BÉCACHE, P. JOLY, AND C. TSOGKA, *Etude d'un nouvel élément fini mixte permettant la condensation de masse*, C. R. Acad. Sci. Paris Sér. I Math., 324 (1997), pp. 1281–1286.

[4] J. BÉRENGER, *A perfectly matched layer for the absorption of electromagnetic waves*, J. Comput. Phys., 114 (1994), pp. 185–200.

[5] F. BREZZI, *On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers*, Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge, 8 (1974), pp. 129–151.

[6] F. BREZZI AND M. FORTIN, *Mixed and Hybrid Finite Element Methods*, Springer-Verlag, New York, 1991.

[7] P. G. CIARLET, *The Finite Element Method for Elliptic Problems*, North-Holland, Amsterdam, The Netherlands, 1978.

[8] G. COHEN AND P. JOLY, *Construction analysis of fourth-order finite difference schemes for the acoustic wave equation in nonhomogeneous media*, SIAM J. Numer. Anal., 33 (1996), pp. 1266–1302.

[9] G. COHEN, P. JOLY, AND N. TORDJMAN, *Construction and analysis of higher order finite elements with mass lumping for the wave equation*, in Proceedings of the Second International Conference on Mathematical and Numerical Aspects of Wave Propagation, Newark, DE, SIAM, Philadelphia, 1993, pp. 152–160.

[10] G. COHEN AND P. MONK, *Gauss point mass lumping schemes for Maxwell's equations*, Numer. Methods Partial Differential Equations, 14 (1998), pp. 63–88.

[11] F. COLLINO, *Perfectly matched absorbing layers for the paraxial equations*, J. Comput. Phys., 131 (1996), pp. 164–180.

[12] F. COLLINO, P. JOLY, AND F. MILLOT, *Fictitious Domain Method for Unsteady Problems: Application to Electromagnetic Scattering*, Tech. Report 2963, INRIA, Rocquencourt, France, 1996.
[13] F. COLLINO AND P. MONK, *Conditions et couches absorbantes pour les équations de Maxwell*, in Aspects récents en méthodes numériques pour les équations de Maxwell, G. Cohen and P. Joly, eds., École des Ondes, INRIA, Rocquencourt, France, 1998.

[14] F. COLLINO AND C. TSOGKA, *Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media*, Geophysics, to appear.

[15] T. DUPONT, *$L^2$-estimates for Galerkin methods for second order hyperbolic equations*, SIAM J. Numer. Anal., 10 (1973), pp. 880–889.

[16] A. ELMKIES AND P. JOLY, *Éléments finis et condensation de masse pour les équations de Maxwell: le cas 2D*, Tech. Report 3035, INRIA, Rocquencourt, France, 1996.

[17] G. A. BAKER AND V. A. DOUGALIS, *The effect of quadrature errors on finite element approximations for the second-order hyperbolic equations*, SIAM J. Numer. Anal., 13 (1976), pp. 577–598.

[18] V. GIRAULT AND R. GLOWINSKI, *Error analysis of a fictitious domain method applied to a Dirichlet problem*, Japan J. Indust. Appl. Math., 12 (1995), pp. 487–514.

[19] J. C. NÉDÉLEC, *A new family of mixed finite elements in $\mathbb{R}^3$*, Numer. Math., 50 (1986), pp. 57–81.

[20] J. E. ROBERTS AND J. M. THOMAS, *Mixed and hybrid methods*, in Handbook of Numerical Analysis, vol. 2, P. G. Ciarlet and J. L. Lions, eds., North-Holland, Amsterdam, The Netherlands, 1989.

[21] P. A. RAVIART AND J. M. THOMAS, *A mixed finite element method for 2nd order elliptic problems*, in Proceedings of Mathematical Aspects on the Finite Element Method, Lecture Notes in Math 606, Springer, Berlin, Heidelberg, New York, 1977, pp. 292–315.

[22] R. GLOWINSKI AND T. W. PAN, AND J. PERIAUX, *A fictitious domain method for Dirichlet problem and applications*, Comput. Methods Appl. Mech. Engrg., 111 (1994), pp. 283–304.

[23] R. M. ALFORD AND K. R. KELLY, AND D. M. BOORE, *Accuracy of finite difference modeling of the acoustic wave equation*, Geophysics, 39 (1974), pp. 834–842.

[24] G. SERIANI AND E. PRIOLO, *Spectral element method for acoustic wave simulation in heterogeneous media*, in Proceedings of ICOSAHOM '92, Montpellier, France, 1992.

[25] N. TORDJMAN, *Éléments finis d'ordre élevé avec condensation de masse pour l'équation des ondes*, Ph.D. thesis, University de Paris IX, Paris, 1995.

[26] W. A. MULDER, *A comparison between higher-order finite elements and finite differences for solving the wave equation*, Numer. Methods Engrg., (1996), pp. 344–350.