A Priori Error Estimates for Mixed Finite Element $\theta$-Schemes for the Wave Equation

A family of implicit-in-time mixed finite element schemes is presented for the numerical approximation of the acoustic wave equation. The mixed space discretization is based on the displacement form of the wave equation and the time-stepping method employs a three-level one-parameter scheme. A rigorous stability analysis is presented based on energy estimation and sharp stability results are obtained. A convergence analysis is carried out and optimal a priori $L^\infty(L^2)$ error estimates for both displacement and pressure are derived.


Introduction
The acoustic wave equation is used to model the effects of wave propagation in heterogeneous media. Solving this equation efficiently is of fundamental importance in many real-life problems. In geophysics, it helps for instance in the interpretation of the seismograph field data and to predict damage patterns due to earthquakes. Using finite element methods for its approximation is attractive because of the ability to handle complex discretizations and design adaptive grid refinement strategies based on error indicators.
Previous attempts on wave simulation by finite elements have used continuous Galerkin methods [1,2,5,8,18,22], mixed finite element methods [6,7,9,13,20,21,26], and discontinuous Galerkin methods [10,14,24,25]. In a mixed finite element formulation both displacements and stresses are approximated simultaneously. This approach provides higher-order approximations to the stresses. This property is important in many problems, in particular in modeling boundary controlability of the wave equation [11]. One of the main difficulties of the mixed finite element techniques is the requirement of compatibility of the approximating spaces for convergence and stability.
The function f f f represents a general source term and u u u 0 and v v v 0 are initial conditions on displacements and velocities.
We assume that f f f , u u u 0 and v v v 0 are smooth enough so that there is a unique solution u u u ∈ C 2 ((0, T ) × Ω) to (1.1)-(1.5), see [17]. The limiting case of (1.1) with µ = 0 is referred to as the acoustic wave equation, which is It is assumed that ρ and λ are bounded below and above by the positive constants ρ 0 , ρ 1 , λ 0 , and λ 1 , respectively. This vector equation is equivalent to the scalar wave equation after making the substitution p = λ∇ · u. The mixed method is established by using this relationship, leading to the coupled system with the appropriate boundary and initial conditions. A priori error estimates for solving (1.7)-(1.8) were obtained in [6,7,9,13]. In [9], Geveci derived L ∞ -in-time, L 2 -in-space error bounds for the continuous-in-time mixed finite element approximations of velocity and stress. In [6,7], a priori error estimates were obtained for the mixed finite element approximation of displacement which requires less regularity than was needed in [9]. Stability for a family of discrete-in-time schemes was also demontratred. In [13], an alternative mixed finite element displacement formulation was proposed reducing requirement on the regularity on the displacement variable. For the explicit discrete-in-time problem, stability results were established and error estimates were obtained. The effectiveness of the method analyzed in [13] was demonstrated in [12] by providing simulations using both lowest-order and next-to-lowest-order Raviart-Thomas elements on rectangles [23].
The purpose of this paper is to analyze an implicit time-stepping method combined with the mixed finite element discretization proposed in [13]. We prove the stability of the proposed method by using energy estimation, and show in particular that it conserves certain energy. We also invertigate the convergence of the method and prove optimal a priori L ∞ (L 2 ) error estimates for both displacement and pressure. The rest of the paper is organized as follows. In sections 2 and 3, we introduce notations and describe the weak formulation of the problem. The fully discrete mixed finite element method is presented in section 4. Stability results are established in section 4 and optimal a priori error estimates are obtained in section 5. Conclusions are given in the last section.

Notation
We shall use the following inner products and norms in this paper. The L 2 -inner product over Ω is defined by inducing the L 2 -norm over Ω, ||v|| L 2 (Ω) = (v, v) 1/2 . The inner product over the boundary ∂Ω is denoted by for u, v ∈ H 1 2 +ε (Ω) with ε > 0. We introduce the time-space norm: The time-space norm || · || L ∞ (L 2 ) is similarly defined. In addition to the L 2 spaces, we use the standard Sobolev space for mixed methods: For the time discretization, we adopt the following notation. Let N be a positive integer, ∆t = T /N , and t n = n∆t. For any function v of time, let v n denote v(t n ). We shall use this notation for functions defined for all times as well as those defined only at discrete times. Set where 0 ≤ θ ≤ 1. We also define the discrete l ∞ -norm for time-discrete functions by ||v|| l ∞ ∆t (0,T ;L 2 (Ω)) = ||v|| l ∞ (L 2 ) = max 0≤n≤N ||v n || L 2 (Ω) .

Weak Formulation
The finite element approximation of the wave problem is based on its weak formulation which is derived in the usual manner. Integrating by parts and using the data on the boundary of Ω, we obtain the weak formulation [13]: For any t ≥ 0, find (u u u(t), where V V V and W are given by The present formulation requires less regularity on the displacement than standard approaches. For instance in [6,7] it is necessary that ∇p ∈ H H H(Ω, div) so that ∇ · u u u ∈ H 2 (Ω). Here, it is only required that ∇ · u u u ∈ H 1 2 , and it can be verified that the solution u u u of problem (1.1)-(1.5) with p = λ∇ · u u u is a solution to (3.4)-(3.5), see [13].
Differentiate (3.5) with respect to time to obtain We next assume f f f = 0 and choose v v v = u u u t and w = p in (3.4) and (3.6), respectively, so that By adding the two equations, we find that Thus, in the absence of forcing, the (continuous) energy is conserved for all time. It will be shown that a similar form of energy is conserved by the numerical solution of the wave problem.

Finite Element Approximation
For the finite element approximation, we let {E h } h>0 be a quasi-uniform family of finite element partitions of Ω, where h is the maximum element diameter. Let V V V h × W h be any of the usual mixed finite element approximating subspaces of V V V ×W , that is, the Raviart-Thomas-Nedelec spaces [19,23], Brezzi-Douglas-Marini spaces [4], or Brezzi-Douglas-Fortin-Marini spaces [3]. For each of these mixed spaces there is a projection Π h : where k is associated with the degree of polynomial and || · || s is the standard Sobolev norm on (H s (Ω)) m .
Here and in what follows, C is a generic positive constant which is independent of h and ∆t.
For φ ∈ W , we denote by P h φ the L 2 -projection of φ onto W h defined by requiring that The semidiscrete mixed finite element approximation to (u u u(t), Existence and uniqueness of a solution (U U U (t), P (t)) to the variational problem (4.5)-(4.9) is shown in [13]. The fully discrete mixed finite element θ-scheme is then defined by finding a sequence of pairs (4.14) Equation (4.12) is derived from the following expansion: The present θ-scheme is explicit in time if θ = 0 and implicit otherwise. The existence and uniqueness of a solution to the resulting linear system for a nonzero value of θ follows from the unisolvancy of the mixed formulation of the following elliptic problem: The explicit case has been considered in [13]. As expected from an explicit scheme, the method is conditionally stable. As a stability constraint, it requires to choose ∆t = O(h).
In the next sections, stability and convergence properties of the proposed θ-scheme are analyzed.

Stability Analysis
We derive sharp stability bounds based on the energy technique and show that the proposed scheme conserves certain energy. We consider (4.13) and (4.14) for the homogeneous case We will make use of the inverse assumption, which states that there exists a constant C 0 independent of h, such that for all φ ∈ W h . The following stability result holds.

4)
and conserves the discrete energy The scheme is unconditionally stable if θ ≥ 1/4.
The case with θ = 1/4 is interesting because the form of the discrete energy in this case is similar to that of the continuous problem. In addition, one can verify that the time truncation error is minimized over the set of all θ ≥ 1/4 when θ = 1/4.

Convergence Analysis
In this section, we prove optimal convergence of the fully discrete finite element solution in the L ∞ (L 2 ) norm. Some of the techniques used in the proofs can be found in previous works [15,16]. In order to estimate the errors in the finite element approximation, we define the auxiliary functions where Π h and P h are defined in Section 4. From (3.4)-(3.5) and (4.13)-(4.14), and the properties of the projections Π h and P h , we arrive at where r r r n = ρ(u u u n;θ tt −∂ tt u u u n ). Another equation has to be derived for the initial errors χ χ χ 1 and ξ 1 . Consider (3.4) at n = 0 and n = 1, respectivey, and subtract the resulting equations so that A use of Taylor's formula with integral remainder yields Using (6.3) and (6.4), we readily obtain Subtracting (6.5) from (4.12) and taking into account (4.11) and (3.4) to arrive that Note that ξ 0 = 0 and χ χ χ 0 = 0. We now state and prove our convergence result.

8)
where r is associated with the degree of the finite element polynomial.
We now distinguish the cases where θ ≥ 1 4 and θ < 1 4 . In the first case, we sum on (6.14) over time levels and multiply through by 2∆t. This results in .

(6.16)
Applying the algebraic inequality: ab ≤ ǫ 2 a 2 + 1 2ǫ b 2 to the right-hand side of (6.16) shows that If we take the supremum over n on the left-hand side and use the fact that N ∆t = T , we conclude that For the case θ < 1 4 , we can follow the analysis presented in [15,16] to derive error estimates similar to (6.18) under condition (6.7).
To complete the proof, we need to bound each term on the right-hand side of (6.18). The first term can be bounded using the approximation properties. Similarly, we have .
For the last term on the right-hand side of (6.18), we have .
Finally, using the approximation property (4.2) and combining all the bounds, we arrive at which completes the proof of the desired estimate.
Remarks. It is worthwhile to mention that the time discretization method is fourth-order accurate when θ = 1/12. To preserve the temporal accuracy of the finite element scheme one has to modify (4.12) carefully to obtain an appropriate initial value U U U 1 . The analysis presented in [15] can be used to derive optimal a priori error estimates in this case.