# Discontinuous Galerkin Finite Element Method for the Wave Equation

Marcus Grote, Anna Schneebeli, Dominik Schötzau

► To cite this version:

Marcus Grote, Anna Schneebeli, Dominik Schötzau. Discontinuous Galerkin Finite Element Method for the Wave Equation. SIAM Journal on Numerical Analysis, 2006, 44, pp.2408 - 2431. 10.1137/05063194X. hal-01443184

HAL Id: hal-01443184

https://hal.science/hal-01443184

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DISCONTINUOUS GALERKIN FINITE ELEMENT METHOD FOR THE WAVE EQUATION

MARCUS J. GROTE$^{†}$, ANNA SCHNEEBELI$^{†}$, AND DOMINIK SCHÖTZAU$^{‡}$

$^{†}$Department of Mathematics, University of Basel, Rheinsprung 21, 4051 Basel, Switzerland (Marcus.Grote@unibas.ch, anna.schneebeli@unibas.ch). The second author was supported by the Swiss National Science Foundation.

$^{‡}$Mathematics Department, University of British Columbia, 121-1984 Mathematics Road, Vancouver V6T 1Z2, BC, Canada (schoetzau@math.ubc.ca). This author was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC).

**Abstract.** The symmetric interior penalty discontinuous Galerkin finite element method is presented for the numerical discretization of the second-order wave equation. The resulting stiffness matrix is symmetric positive definite, and the mass matrix is essentially diagonal; hence, the method is inherently parallel and leads to fully explicit time integration when coupled with an explicit time-stepping scheme. Optimal a priori error bounds are derived in the energy norm and the $L^2$-norm for the semidiscrete formulation. In particular, the error in the energy norm is shown to converge with the optimal order $O(h^{\min\{s,l\}})$ with respect to the mesh size $h$, the polynomial degree $l$, and the regularity exponent $s$ of the continuous solution. Under additional regularity assumptions, the $L^2$-error is shown to converge with the optimal order $O(h^{l+1})$. Numerical results confirm the expected convergence rates and illustrate the versatility of the method.

**Key words.** discontinuous Galerkin finite element methods, wave equation, acoustic waves, second-order hyperbolic problems, a priori error analysis, explicit time integration

**1. Introduction.** The numerical solution of the wave equation is of fundamental importance to the simulation of time dependent acoustic, electromagnetic, or elastic waves. For such wave phenomena the scalar second-order wave equation often serves as a model problem. Finite element methods (FEMs) can easily handle inhomogeneous media or complex geometry. However, if explicit time-stepping is subsequently employed, the mass matrix arising from the spatial discretization by standard continuous finite elements must be inverted at each time step: a major drawback in terms of efficiency. For low-order Lagrange ($P^1$) elements, so-called mass lumping overcomes this problem [6, 15], but for higher-order elements this procedure can lead to unstable schemes unless particular finite elements and quadrature rules are used [11]. In addition, continuous Galerkin methods impose significant restrictions on the underlying mesh and discretization; in particular, they do not easily accommodate hanging nodes.

To avoid these difficulties, we consider instead discontinuous Galerkin (DG) methods. Based on discontinuous finite element spaces, these methods easily handle elements of various types and shapes, irregular nonmatching grids, and even locally varying polynomial order; thus, they are ideally suited for *hp*-adaptivity. Here continuity is weakly enforced across mesh interfaces by adding suitable bilinear forms, so-called numerical fluxes, to standard variational formulations. These fluxes are easily included within an existing conforming finite element code.
Because individual elements decouple, DGFEMs are also inherently parallel; see [8, 9, 10, 7] for further details and recent reviews. Moreover, the mass matrix arising from the spatial DG discretization is block-diagonal, with block size equal to the number of degrees of freedom per element; it can therefore be inverted at very low computational cost. In fact, for a judicious choice of (locally orthogonal) shape functions, the mass matrix is diagonal. When combined with explicit time integration, the resulting time marching scheme will be fully explicit.

The origins of DG methods can be traced back to the 1970s, when they were proposed for the numerical solution of hyperbolic neutron transport equations, as well as for the weak enforcement of continuity in Galerkin methods for elliptic and parabolic problems; see Cockburn, Karniadakis, and Shu [8] for a review of the development of DG methods. When applied to second-order hyperbolic problems, most DG methods first require the problem to be reformulated as a first-order hyperbolic system, for which various DG methods are available. In [9], for instance, Cockburn and Shu used a DGFEM in space combined with a Runge-Kutta scheme in time to discretize hyperbolic conservation laws. Hesthaven and Warburton [13] used the same approach to implement high-order methods for Maxwell's equations in first-order hyperbolic form. Space-time DG methods for linear symmetric first-order hyperbolic systems were presented by Falk and Richter in [12] and later generalized by Monk and Richter in [17] and by Houston, Jensen, and Süli in [14]. A first DG method for the acoustic wave equation in its original second-order formulation was recently proposed by Rivière and Wheeler [21]; it is based on a nonsymmetric interior penalty formulation and requires additional stabilization terms for optimal convergence in the $L^2$-norm [20].

Here we propose and analyze the symmetric interior penalty DG method for the spatial discretization of the second-order scalar wave equation. In particular, we shall derive optimal a priori error bounds in the energy norm and the $L^2$-norm for the semidiscrete formulation. Besides the well-known advantages of DG methods mentioned above, a symmetric discretization of the wave equation in its second-order form offers an additional advantage, which also pertains to the classical continuous Galerkin formulation: since the stiffness matrix is positive definite, the semidiscrete formulation conserves (a discrete version of) the energy for all time; thus, it is free of any (unnecessary) damping. The dispersive properties of the symmetric interior penalty DG method were recently analyzed by Ainsworth, Monk, and Muniz [1].

The outline of our paper is as follows. In section 2 we describe the setting of our model problem. Next, we present in section 3 the symmetric interior penalty DG method for the wave equation. Our two main results, optimal error bounds in the energy norm and the $L^2$-norm for the semidiscrete scheme, are stated at the beginning of section 4 and proved subsequently. The analysis relies on an idea suggested by Arnold et al. [2] together with the approach presented by Perugia and Schötzau in [18] to extend the DG bilinear forms by suitable lifting operators. In section 5, we demonstrate the sharpness of our theoretical error estimates by a series of numerical experiments. By combining our DG method with the second-order Newmark scheme, we obtain a fully discrete method. To illustrate the versatility of our method, we also propagate a wave across an inhomogenous medium with discontinuity, where the underlying finite element mesh contains hanging nodes. Finally, we conclude with some remarks on possible extensions of our DG method to electromagnetic and elastic waves.

## 2. Model problem. We consider the (second-order) scalar wave equation
$$ (2.1) \qquad u_{tt} - \nabla \cdot (c \nabla u) = f \qquad \text{in } J \times \Omega, $$

$$ (2.2) \qquad u = 0 \qquad \text{on } J \times \partial\Omega, $$

$$ (2.3) \qquad u|_{t=0} = u_0 \qquad \text{in } \Omega, $$

$$ (2.4) \qquad u_t|_{t=0} = v_0 \qquad \text{in } \Omega, $$

where $J = (0, T)$ is a finite time interval and $\Omega$ is a bounded domain in $\mathbb{R}^d$, $d = 2, 3$. For simplicity, we assume that $\Omega$ is a polygon ($d = 2$) or a polyhedron ($d = 3$). The (known) source term $f$ lies in $L^2(J; L^2(\Omega))$, while $u_0 \in H_0^1(\Omega)$ and $v_0 \in L^2(\Omega)$ are prescribed initial conditions. We assume that the speed of propagation, $\sqrt{c(x)}$, is piecewise smooth and satisfies the bounds

$$ (2.5) \qquad 0 < c_* \le c(x) \le c^* < \infty, \qquad x \in \bar{\Omega}. $$

The standard variational form of (2.1)–(2.4) is to find $u \in L^2(J; H_0^1(\Omega))$, with $u_t \in L^2(J; L^2(\Omega))$ and $u_{tt} \in L^2(J; H^{-1}(\Omega))$, such that $u|_{t=0} = u_0$, $u_t|_{t=0} = v_0$, and

$$ (2.6) \qquad \langle u_{tt}, v \rangle + a(u, v) = (f, v) \qquad \forall v \in H_0^1(\Omega) \quad \text{a.e. in } J. $$

Here, the time derivatives are understood in a distributional sense, $\langle \cdot, \cdot \rangle$ denotes the duality pairing between $H^{-1}(\Omega)$ and $H_0^1(\Omega)$, $(\cdot, \cdot)$ is the inner product in $L^2(\Omega)$, and $a(\cdot, \cdot)$ is the elliptic bilinear form given by

$$ (2.7) \qquad a(u, v) = (c \nabla u, \nabla v). $$

It is well known that problem (2.6) is well posed [16]. Moreover, the weak solution $u$ can be shown to be continuous in time; that is,

$$ (2.8) \qquad u \in C^0(\bar{J}; H_0^1(\Omega)), \qquad u_t \in C^0(\bar{J}; L^2(\Omega)); $$

see [16, Chapter III, Theorems 8.1 and 8.2] for details. In particular, this result implies that the initial conditions in (2.3) and (2.4) are well defined.

**3. Discontinuous Galerkin discretization.** We shall now discretize the wave equation (2.1)–(2.4) by using the interior penalty discontinuous Galerkin finite element method in space, while leaving the time dependence continuous.

**3.1. Preliminaries.** We consider shape-regular meshes $\mathcal{T}_h$ that partition the domain $\Omega$ into disjoint elements $\{K\}$ such that $\bar{\Omega} = \cup_{K \in \mathcal{T}_h} \bar{K}$. For simplicity, we assume that the elements are triangles or parallelograms in two space dimensions, and tetrahedra or parallelepipeds in three dimensions, respectively. The diameter of element $K$ is denoted by $h_K$, and the mesh size $h$ is given by $h = \max_{K \in \mathcal{T}_h} h_K$. We assume that the partition is aligned with the discontinuities of the wave speed $\sqrt{c}$. Generally, we allow for irregular meshes with hanging nodes. However, we assume that the local mesh sizes are of bounded variation; that is, there is a positive constant $\kappa$, depending only on the shape-regularity of the mesh, such that

$$ (3.1) \qquad \kappa h_K \le h_{K'} \le \kappa^{-1} h_K $$

for all neighboring elements $K$ and $K'$.

An interior face of $\mathcal{T}_h$ is the (nonempty) interior of $\partial K^+ \cap \partial K^-$, where $K^+$ and $K^-$ are two adjacent elements of $\mathcal{T}_h$. Similarly, a boundary face of $\mathcal{T}_h$ is the (nonempty)
interior of $\partial K \cap \partial \Omega$, which consists of entire faces of $\partial K$. We denote by $\mathcal{F}_h^I$ the set of all interior faces of $\mathcal{T}_h$ and by $\mathcal{F}_h^B$ the set of all boundary faces and set $\mathcal{F}_h = \mathcal{F}_h^I \cup \mathcal{F}_h^B$. Here we generically refer to any element of $\mathcal{F}_h$ as a “face,” in both two and three dimensions.

For any piecewise smooth function $v$ we now introduce the following trace operators. Let $F \in \mathcal{F}_h^I$ be an interior face shared by two neighboring elements $K^+$ and $K^-$ and let $x \in F$; we write $\mathbf{n}^{\pm}$ to denote the unit outward normal vectors on the boundaries $\partial K^{\pm}$. Denoting by $v^{\pm}$ the trace of $v$ taken from within $K^{\pm}$, we define the jump and average of $v$ at $x \in F$ by

$$ \llbracket v \rrbracket := v^+ \mathbf{n}^+ + v^- \mathbf{n}^-, \quad \llbracket v \rrbracket := (v^+ + v^-)/2, $$

respectively. On every boundary face $F \in \mathcal{F}_h^B$, we set $\llbracket v \rrbracket := v\mathbf{n}$ and $\llbracket v \rrbracket := v$. Here, $\mathbf{n}$ is the unit outward normal vector on $\partial \Omega$.

For a piecewise smooth vector-valued function $\mathbf{q}$, we analogously define the average across interior faces by $\llbracket \mathbf{q} \rrbracket := (\mathbf{q}^+ + \mathbf{q}^-)/2$, and on boundary faces we set $\llbracket \mathbf{q} \rrbracket := \mathbf{q}$. The jump of a vector-valued function will not be used. For a vector-valued function $\mathbf{q}$ with continuous normal components across a face $f$, the trace identity

$$ v^+(\mathbf{n}^+ \cdot \mathbf{q}^+) + v^-(\mathbf{n}^- \cdot \mathbf{q}^-) = \llbracket v \rrbracket \cdot \llbracket \mathbf{q} \rrbracket \quad \text{on } f $$

immediately follows from the definitions.

**3.2. Discretization in space.** For a given partition $\mathcal{T}_h$ of $\Omega$ and an approximation order $\ell \ge 1$, we wish to approximate the solution $u(t, \cdot)$ of (2.1)–(2.4) in the finite element space

$$ (3.2) \qquad V^h := \{v \in L^2(\Omega) : v|_{K} \in S^\ell(K) \ \forall K \in \mathcal{T}_h\}, $$

where $S^\ell(K)$ is the space $P^\ell(K)$ of polynomials of total degree at most $\ell$ on $K$ if $K$ is a triangle or a tetrahedra, or the space $Q^\ell(K)$ of polynomials of degree at most $\ell$ in each variable on $K$ if $K$ is a parallelogram or a parallelepiped.

Then, we consider the following (semidiscrete) DG approximation of (2.1)–(2.4): find $u^h: \bar{J} \times V^h \to \mathbb{R}$ such that

$$ (3.3) \qquad (u_{tt}^h, v) + a_h(u^h, v) = (f, v) \quad \forall v \in V^h, \quad t \in J, $$

$$ (3.4) \qquad u^h|_{t=0} = \Pi_h u_0, $$

$$ (3.5) \qquad u_t^h|_{t=0} = \Pi_h v_0. $$

Here, $\Pi_h$ denotes the $L^2$-projection onto $V^h$, and the discrete bilinear form $a_h$ on $V^h \times V^h$ is given by

$$ (3.6) \qquad a_h(u, v) := \sum_{K \in \mathcal{T}_h} \int_K c \nabla u \cdot \nabla v \, dx - \sum_{F \in \mathcal{F}_h} \int_F \llbracket u \rrbracket \cdot \llbracket c \nabla v \rrbracket \, dA \\ \qquad \qquad - \sum_{F \in \mathcal{F}_h} \int_F \llbracket v \rrbracket \cdot \llbracket c \nabla u \rrbracket \, dA + \sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[\llbracket u \rrbracket] \cdot \llbracket v \rrbracket \, dA. $$

The last three terms in (3.6) correspond to jump and flux terms at element boundaries; they vanish when $u, v \in H_0^1(\Omega) \cap H^{1+\sigma}(\Omega)$ for $\sigma > \frac{1}{2}$. Hence the above semidiscrete DG formulation (3.3) is consistent with the original continuous problem (2.6).
In (3.6) the function a penalizes the jumps of u and v over the faces of $\mathcal{T}_h$. It is referred to as the interior penalty stabilization function and is defined as follows. We first introduce the function h by

$$h|_F = \begin{cases} \min\{h_K, h_{K'}\}, & F \in \mathcal{F}_h^I, F = \partial K \cap \partial K', \\ h_K, & F \in \mathcal{F}_h^B, F = \partial K \cap \partial \Omega. \end{cases}$$

For $x \in F$, we further define c by

$$c|_F(x) = \begin{cases} \max\{c|_K(x), c|_{K'}(x)\}, & F \in \mathcal{F}_h^I, F = \partial K \cap \partial K', \\ c|_K(x), & F \in \mathcal{F}_h^B, F = \partial K \cap \partial \Omega. \end{cases}$$

Then, on each $F \in \mathcal{F}_h$, we set

$$(3.7) \qquad a|_F := \alpha c h^{-1},$$

where $\alpha$ is a positive parameter independent of the local mesh sizes and the coefficient c.

To conclude this section we recall the following stability result for the DG form $a_h$.

**LEMMA 3.1.** There exists a threshold value $\alpha_{\min} > 0$ which depends only on the shape-regularity of the mesh, the approximation order $\ell$, the dimension $d$, and the bounds in (2.5), such that for $\alpha \ge \alpha_{\min}$

$$a_h(v, v) \ge C_{\text{coer}} \left( \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla v\|_{0, K}^2 + \sum_{F \in \mathcal{F}_h} \|a^{\frac{1}{2}}[v]\|_{0, F}^2 \right), \quad v \in V^h,$$

where the constant $C_{\text{coer}}$ is independent of c and h.

The proof of this lemma follows readily from the arguments in [2]. However, to make explicit the dependence of $\alpha_{\min}$ on the bounds in (2.5), we present the proof of a slightly more general stability result in Lemma 4.4 below. Throughout the rest of the paper we shall assume that $\alpha \ge \alpha_{\min}$, so that by Lemma 3.1 the semidiscrete problem (3.3)–(3.5) has a unique solution.

We remark that the condition $\alpha \ge \alpha_{\min}$ can be omitted by using other symmetric DG discretizations of the div-grad operator, such as the local discontinuous Galerkin (LDG) method; see, e.g., [2] for details. It can also be avoided by using the nonsymmetric interior penalty method proposed in [20]. However, since the symmetry of $a_h$ is crucial in the analysis below, our error estimates (section 4) do not hold for the nonsymmetric DG method in [20].

**Remark 3.2.** Because the bilinear form $a_h$ is symmetric and coercive, for $\alpha \ge \alpha_{\min}$, the semidiscrete DG formulation (3.3)–(3.5) with $f = 0$ conserves the (discrete) energy

$$E_h(t) := \frac{1}{2} \|u_t^h(t)\|_0^2 + \frac{1}{2} a_h(u^h(t), u^h(t)).$$

**4. A priori error estimates.** We shall now derive optimal a priori error bounds for the DG method (3.3)–(3.5), first with respect to the DG energy norm and then with respect to the $L^2$-norm. These two key results are stated immediately below, while their proofs are postponed to subsequent sections.
**4.1. Main results.** To state our a priori error bounds, we define the space

$$V(h) = H_0^1(\Omega) + V^h.$$

On $V(h)$, we define the DG energy norm

$$\|v\|_h^2 := \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla v\|_{0,K}^2 + \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2.$$

Furthermore, for $1 \le p \le \infty$ we will make use of the Bochner space $L^p(J; V(h))$, endowed with the norm

$$\|v\|_{L^p(J;V(h))} = \begin{cases} (\int_J \|v\|_h^p dt)^{1/p}, & 1 \le p < \infty, \\ \ess \sup_{t \in J} \|v\|_h, & p = \infty. \end{cases}$$

Our first main result establishes an optimal error estimate of the energy norm $\| \cdot \|_h$ of the error. It also gives a bound in the $L^2(\Omega)$-norm on the error in the first time derivative.

**THEOREM 4.1.** Let the analytical solution $u$ of (2.1)–(2.4) satisfy

$$u \in L^{\infty}(J; H^{1+\sigma}(\Omega)), \quad u_t \in L^{\infty}(J; H^{1+\sigma}(\Omega)), \quad u_{tt} \in L^1(J; H^{\sigma}(\Omega))$$

for a regularity exponent $\sigma > \frac{1}{2}$, and let $u^h$ be the semidiscrete discontinuous Galerkin approximation obtained by (3.3)–(3.5), with $\alpha \ge \alpha_{\min}$. Then, the error $e = u - u^h$ satisfies the estimate

$$\|e_t\|_{L^{\infty}(J;L^2(\Omega))} + \|e\|_{L^{\infty}(J;V(h))} \le C \left[ \|e_t(0)\|_0 + \|e(0)\|_h \right] \\ + C h^{\min\{\sigma,\ell\}} \left[ \|u\|_{L^{\infty}(J;H^{1+\sigma}(\Omega))} + T \|u_t\|_{L^{\infty}(J;H^{1+\sigma}(\Omega))} + \|u_{tt}\|_{L^1(J;H^{\sigma}(\Omega))} \right],$$

with a constant $C$ that is independent of $T$ and $h$.

We remark that the fact that $u_t \in L^\infty(J; H^{1+\sigma}(\Omega))$ implies that $u$ is continuous in time on $\bar{J}$ with values in $H^{1+\sigma}(\Omega)$. Similarly, $u_{tt} \in L^1(J; H^\sigma(\Omega))$ implies the continuity of $u_t$ on $\bar{J}$ with values in $H^\sigma(\Omega)$. In Theorem 4.1 we thus implicitly assume that the initial conditions satisfy $u_0 \in H^{1+\sigma}(\Omega)$ and $v_0 \in H^\sigma(\Omega)$. Hence, standard approximation properties imply that

$$\|e_t(0)\|_0 = \|v_0 - \Pi_h v_0\|_0 \le C h^{\min\{\sigma,\ell+1\}} \|v_0\|_{\sigma},$$

$$\|e(0)\|_h = \|u_0 - \Pi_h u_0\|_h \le C h^{\min\{\sigma,\ell\}} \|u_0\|_{1+\sigma};$$

see also Lemma 4.6 below. As a consequence, Theorem 4.1 yields optimal convergence in the (DG) energy norm

$$\|e_t\|_{L^{\infty}(J;L^2(\Omega))} + \|e\|_{L^{\infty}(J;V(h))} \le C h^{\min\{\sigma,\ell\}},$$

with a constant $C = C(T)$ that is independent of $h$.

Next, we state an optimal error estimate with respect to the $L^2$-norm (in space). To do so, we need to assume elliptic regularity; that is, we assume that there is a stability constant $C_S$ such that for any $\lambda \in L^2(\Omega)$ the solution of the problem

$$(4.1) \qquad -\nabla \cdot (c\nabla z) = \lambda \quad \text{in } \Omega, \quad z = 0 \quad \text{on } \Gamma,$$
belongs to $H^2(\Omega)$ and satisfies the stability bound

$$ (4.2) \qquad \|z\|_2 \le C_S \|\lambda\|_0. $$

This condition is certainly satisfied for convex domains and smooth coefficients. Then, the following $L^2$-error bound holds.

**THEOREM 4.2.** Assume elliptic regularity as in (4.1)–(4.2), and let the analytical solution $u$ of (2.1)–(2.4) satisfy

$$ u \in L^\infty(J; H^{1+\sigma}(\Omega)), \quad u_t \in L^\infty(J; H^{1+\sigma}(\Omega)), \quad u_{tt} \in L^1(J; H^\sigma(\Omega)) $$

for a regularity exponent $\sigma > \frac{1}{2}$. Let $u^h$ be the semidiscrete DG approximation obtained by (3.3)–(3.5) with $\alpha \ge \alpha_{\min}$. Then, the error $e = u - u^h$ satisfies the estimate

$$ \|e\|_{L^\infty(J;L^2(\Omega))} \le C h^{\min\{\sigma,\ell\}+1} \left[ \|u_0\|_{1+\sigma} + \|u\|_{L^\infty(J;H^{1+\sigma}(\Omega))} + T \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))} \right], $$
with a constant $C$ that is independent of $T$ and the mesh size.

For smooth solutions, Theorem 4.2 thus yields optimal convergence rates in the $L^2$-norm:

$$ \|e\|_{L^\infty(J;L^2(\Omega))} \le C h^{\ell+1}, $$

with a constant $C$ that is independent of $h$.

The rest of this section is devoted to the proofs of Theorems 4.1 and 4.2. We shall first collect preliminary results in section 4.2. In section 4.3, we present the proof of Theorem 4.1. Following an argument by Baker [3] for conforming finite element approximations, we shall then derive the estimate of Theorem 4.2 in section 4.4.

## 4.2. Preliminaries.

**Extension of the DG form $a_h$.** The DG form $a_h$ in (3.6) does not extend in a standard way to a continuous form on the (larger) space $V(h) \times V(h)$. Indeed the average $\|c\nabla v\|$ on a face $F \in \mathcal{F}_h$ is not well defined in general for $v \in H^1(\Omega)$. To circumvent this difficulty, we shall extend the form $a_h$ in a nonstandard and nonconsistent way to the space $V(h) \times V(h)$ by using the lifting operators from [2] and the approach in [18]. Thus, for $v \in V(h)$ we define the lifted function, $\mathcal{L}_c(v) \in (V^h)^d$, $d = 2, 3$, by requiring that

$$ \int_{\Omega} \mathcal{L}_c(v) \cdot \mathbf{w} \, dx = \sum_{F \in \mathcal{F}_h} \int_F \|v\| \cdot \|c\mathbf{w}\| \, dA, \quad \mathbf{w} \in (V^h)^d, $$

where $c$ is the material coefficient from (2.1). We shall now show that the lifting operator $\mathcal{L}_c$ is stable in the DG norm; see [18] for a similar result for the LDG method.

**LEMMA 4.3.** There exists a constant $C_{\text{inv}}$ which depends only on the shape-regularity of the mesh, the approximation order $\ell$, and the dimension $d$ such that

$$ \|\mathcal{L}_c(v)\|_0^2 \le \alpha^{-1} c^* C_{\text{inv}}^2 \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2 $$

for any $v \in V(h)$.
Moreover, if the speed of propagation $c^{\frac{1}{2}}$ is piecewise constant, with discontinuities aligned with the finite element mesh $\mathcal{T}_h$, then

$$ \|c^{-\frac{1}{2}} \mathcal{L}_c(v)\|_0^2 \leq \alpha^{-1} C_{\text{inv}}^2 \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2. $$

*Proof.* We have

$$ \begin{aligned} \|\mathcal{L}_c(v)\|_0 &= \max_{\mathbf{w} \in (V^h)^d} \frac{\sum_{F \in \mathcal{F}_h} \int_F [v] \cdot \{\!\{c\mathbf{w}\}\!\} dA}{\|\mathbf{w}\|_0} \\ &\leq \max_{\mathbf{w} \in (V^h)^d} \frac{\left(\sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[v]^2 dA\right)^{\frac{1}{2}} \left(\sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}^{-1}|\{\!\{c\mathbf{w}\}\!\}|^2 dA\right)^{\frac{1}{2}}}{\|\mathbf{w}\|_0} \\ &\leq \alpha^{-\frac{1}{2}} \max_{\mathbf{w} \in (V^h)^d} \frac{\left(\sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[v]^2 dA\right)^{\frac{1}{2}} \left(\sum_{F \in \mathcal{F}_h} \int_F \mathbf{h}\mathbf{c}^{-1}|\{\!\{c\mathbf{w}\}\!\}|^2 dA\right)^{\frac{1}{2}}}{\|\mathbf{w}\|_0} \\ &\leq \alpha^{-\frac{1}{2}} (c^*)^{\frac{1}{2}} \max_{\mathbf{w} \in (V^h)^d} \frac{\left(\sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[v]^2 dA\right)^{\frac{1}{2}} \left(\sum_{K \in \mathcal{T}_h} h_K \int_{\partial K} |\mathbf{w}|^2 dA\right)^{\frac{1}{2}}}{\|\mathbf{w}\|_0}. \end{aligned} $$

Here, we have used the Cauchy–Schwarz inequality, the definition of $\mathbf{a}$ in (3.7), and the upper bound for $c$ in (2.5). We recall the inverse inequality

$$ (4.3) \qquad \|\mathbf{w}\|_{0,\partial K}^2 \leq C_{\text{inv}}^2 h_K^{-1} \|\mathbf{w}\|_{0,K}^2, \quad \mathbf{w} \in (S^\ell(K))^d, $$

with a constant $C_{\text{inv}}$ that depends only on the shape-regularity of the mesh, the approximation order $\ell$, and the dimension $d$. Using this bound, we obtain

$$ \left( \sum_{K \in \mathcal{T}_h} h_K \int_{\partial K} |\mathbf{w}|^2 dA \right)^{\frac{1}{2}} \leq C_{\text{inv}} \|\mathbf{w}\|_0, $$

which shows the first statement.

With $c^{\frac{1}{2}}$ piecewise constant, we have $c^{-\frac{1}{2}}z \in (V^h)^d$ for all $z \in (V^h)^d$. Hence, we obtain as before

$$ \begin{aligned} \|c^{-\frac{1}{2}} \mathcal{L}_c(v)\|_0 &= \max_{\mathbf{w} \in (V^h)^d} \frac{\sum_{F \in \mathcal{F}_h} \int_F [v] \cdot \{\!\{c^{\frac{1}{2}}\mathbf{w}\}\!\} dA}{\|\mathbf{w}\|_0} \\ &\leq \alpha^{-1} C_{\text{inv}}^2 \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2, \end{aligned} $$

which completes the proof. □

Next, we introduce the auxiliary bilinear form

$$ (4.4) \qquad \begin{aligned} \tilde{a}_h(u, v) := & \sum_{K \in \mathcal{T}_h} \int_K c \nabla u \cdot \nabla v \,dx - \sum_{K \in \mathcal{T}_h} \int_K \mathcal{L}_c(u) \cdot \nabla v \,dx \\ & - \sum_{K \in \mathcal{T}_h} \int_K \mathcal{L}_c(v) \cdot \nabla u \,dx + \sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[u] \cdot [v] \,dA. \end{aligned} $$

The following result establishes that $\tilde{a}_h$ is continuous and coercive on the entire space $V(h) \times V(h)$; hence it is well defined. Furthermore, since

$$ (4.5) \qquad \tilde{a}_h = a_h \quad \text{on } V^h \times V^h, \quad \tilde{a}_h = a \quad \text{on } H_0^1(\Omega) \times H_0^1(\Omega), $$
the form $\tilde{a}_h$ can be viewed as an extension of the two forms $a_h$ and $a$ to the space $V(h) \times V(h)$.

LEMMA 4.4. Let the interior penalty parameter $a$ be defined as in (3.7), and set

$$ \alpha_{\min} = 4 c_{\star}^{-1} c^{\star} C_{\text{inv}}^{2} $$

for a general piecewise smooth $c$, and

$$ \alpha_{\min} = 4 C_{\text{inv}}^{2} $$

for a piecewise constant $c$, with discontinuities aligned with the finite element mesh $T_h$. $C_{\text{inv}}$ is the constant from Lemma 4.3.

Setting $C_{\text{cont}} = 2$ and $C_{\text{coer}} = 1/2$, we have for $\alpha \ge \alpha_{\min}$

$$ \begin{aligned} |\tilde{a}_h(u, v)| &\le C_{\text{cont}} \|u\|_h \|v\|_h, & u, v \in V(h), \\ \tilde{a}_h(u, u) &\ge C_{\text{coer}} \|u\|_h^2, & u \in V(h). \end{aligned} $$

In particular, the coercivity bound implies the result in Lemma 3.1.

*Proof.* By taking into account the bounds in (2.5) and Lemma 4.3, application of the Cauchy–Schwarz inequality readily gives in the general case

$$ |\tilde{a}_h(u, v)| \le \max\{2, \alpha^{-1}c_{\star}^{-1}c^{\star}C_{\text{inv}}^2 + 1\} \|u\|_h \|v\|_h. $$

For $\alpha \ge \alpha_{\min}$, the continuity of $\tilde{a}_h$ immediately follows. The case of piecewise constant $c$ follows analogously.

To show the coercivity of the form $\tilde{a}_h$, we note that

$$ \tilde{a}_h(u, u) = \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla u\|_{0,K}^2 - 2 \sum_{K \in \mathcal{T}_h} \int_K \mathcal{L}_c(u) \cdot \nabla u \, dx + \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[\llbracket u \rrbracket]\|_{0,F}^2. $$

By using the weighted Cauchy–Schwarz inequality, the geometric-arithmetic inequality $ab \le \frac{\varepsilon a^2}{2} + \frac{b^2}{2\varepsilon}$, valid for any $\varepsilon > 0$, the bounds in (2.5), and the stability bound for the lifting operator in Lemma 4.3, we obtain for general $c$

$$ \begin{aligned} 2 \sum_{K \in \mathcal{T}_h} \int_K \mathcal{L}_c(u) \cdot \nabla u \, dx &= 2 \sum_{K \in \mathcal{T}_h} \int_K c^{-\frac{1}{2}} \mathcal{L}_c(u) \cdot c^{\frac{1}{2}} \nabla u \, dx \\ &\le 2 \sum_{K \in \mathcal{T}_h} \|c^{-\frac{1}{2}} \mathcal{L}_c(u)\|_{0,K} \|c^{\frac{1}{2}} \nabla u\|_{0,K} \\ &\le \varepsilon \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla u\|_{0,K}^2 + \varepsilon^{-1} c_{\star}^{-1} \sum_{K \in \mathcal{T}_h} \|\mathcal{L}_c(u)\|_{0,K}^2 \\ &\le \varepsilon \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla u\|_{0,K}^2 + \varepsilon^{-1} \alpha^{-1} c_{\star}^{-1} c^{\star} C_{\text{inv}}^2 \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[\llbracket u \rrbracket]\|_{0,F}^2 \end{aligned} $$

for a parameter $\varepsilon > 0$ still at our disposal. We conclude that

$$ \tilde{a}_h(u, u) \ge (1 - \varepsilon) \sum_{K \in \mathcal{T}_h} \|c^{\frac{1}{2}} \nabla u\|_{0,K}^2 + (1 - \varepsilon^{-1} \alpha^{-1} c_{\star}^{-1} c^{\star} C_{\text{inv}}^2) \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[\llbracket u \rrbracket]\|_{0,F}^2. $$

For $\varepsilon = \frac{1}{2}$ and $\alpha \ge \alpha_{\min}$, we obtain the desired coercivity bound.

For a piecewise constant $c$ we use the bound for $\|c^{-\frac{1}{2}} \mathcal{L}_c(u)\|_0^2$ from Lemma 4.4 and proceed analogously. $\square$
**Error equation.** Because $\tilde{a}_h$ coincides with $a_h$ on $V^h \times V^h$, the semidiscrete scheme in (3.3)-(3.5) is equivalent to the following:

Find $u^h: \bar{J} \times V^h \to \mathbb{R}$ such that $u^h|_{t=0} = \Pi_h u_0$, $u_t^h|_{t=0} = \Pi_h v_0$, and

$$ (4.6) \qquad (u_{tt}^h, v) + \tilde{a}_h(u^h, v) = (f, v) \quad \forall v \in V^h. $$

We shall use the formulation in (4.6) as the basis of our error analysis.

To derive an error equation, we first define for $u \in H^{1+\sigma}(\Omega)$ with $\sigma > 1/2$

$$ (4.7) \qquad r_h(u; v) = \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \nabla u - c \Pi_h(\nabla u)\}\} dA, \quad v \in V(h). $$

Here $\Pi_h$ denotes the $L^2$-projection onto $(V^h)^d$. The assumption $u \in H^{1+\sigma}(\Omega)$ ensures that $r_h(u; v)$ is well defined. From the definition in (4.7) it is immediate that $r_h(u; v) = 0$ when $v \in H_0^1(\Omega)$.

**LEMMA 4.5.** Let the analytical solution $u$ of (2.1)-(2.4) satisfy

$$ u \in L^\infty(J; H^{1+\sigma}(\Omega)), \quad u_{tt} \in L^1(J; L^2(\Omega)). $$

Let $u^h$ be the semidiscrete DG approximation obtained by (4.6). Then, the error $e = u - u^h$ satisfies

$$ (e_{tt}, v) + \tilde{a}_h(e, v) = r_h(u; v) \quad \forall v \in V^h \text{ a.e. in } J, $$

with $r_h(u; v)$ given in (4.7).

**Proof.** Let $v \in V^h$. Since $u_{tt} \in L^1(J; L^2(\Omega))$, we have $\langle u_{tt}, v \rangle = (u_{tt}, v)$ almost everywhere in $J$. Hence, using the discrete formulation in (3.3)-(3.5), we obtain that

$$ (e_{tt}, v) + \tilde{a}_h(e, v) = (u_{tt}, v) + \tilde{a}_h(u, v) - (f, v) \quad \text{a.e. in } J. $$

Now, by definition of $\tilde{a}_h$, the fact that $\mathcal{L}_c(u) = 0$ and that $[[u]] = 0$ on all faces, the defining properties of the $L^2$-projection $\Pi_h$, and the definition of the lifted element $\mathcal{L}_c(v)$, we obtain

$$ \tilde{a}_h(u, v) = \sum_{K \in \mathcal{T}_h} \int_K c \nabla u \cdot \nabla v \, dx - \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \Pi_h(\nabla u)\}\} dA. $$

Since $u_{tt} \in L^1(J; L^2(\Omega))$ and $f \in L^2(J; L^2(\Omega))$, we have that $\nabla \cdot (c \nabla u) \in L^2(\Omega)$ almost everywhere in $J$, which implies that $c \nabla u$ has continuous normal components across all interior faces. Therefore, elementwise integration by parts combined with the trace operators defined in section 3.1 yields

$$ \begin{aligned} \tilde{a}_h(u, v) = & - \sum_{K \in \mathcal{T}_h} \int_K \nabla \cdot (c \nabla u) v \, dx + \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \nabla u\}\} dA \\ & - \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \Pi_h(\nabla u)\}\} dA. \end{aligned} $$

From the definition of $r_h(u, v)$ in (4.7), we therefore conclude that

$$ (u_{tt}, v) + \tilde{a}_h(u, v) = (u_{tt} - \nabla \cdot (c \nabla u), v) + r_h(u; v) $$

and obtain

$$ (e_{tt}, v) + \tilde{a}_h(e, v) = (u_{tt} - \nabla \cdot (c \nabla u) - f, v) + r_h(u; v) = r_h(u; v), $$

where we have used the differential equation (2.1). $\square$
**Approximation properties.** Let $\Pi_h$ and $\Pi_h$ denote the $L^2$-projections onto $V^h$ and $(V^h)^d$, respectively. We recall the following approximation properties; see [6].

LEMMA 4.6. Let $K \in \mathcal{T}_h$. Then the following hold:

(i) For $v \in H^t(K)$, $t \ge 0$, we have

$$ \|v - \Pi_h v\|_{0,K} \le C h_K^{\min\{t, \ell+1\}} \|v\|_{t,K}, $$

with a constant $C$ that is independent of the local mesh size $h_K$ and depends only on the shape-regularity of the mesh, the approximation order $\ell$, the dimension $d$, and the regularity exponent $t$.

(ii) For $v \in H^{1+\sigma}(K)$, $\sigma > \frac{1}{2}$, we have

$$ \|\nabla v - \nabla(\Pi_h v)\|_{0,K} \le C h_K^{\min\{\sigma, \ell\}} \|v\|_{1+\sigma,K}, $$

$$ \|v - \Pi_h v\|_{0,\partial K} \le C h_K^{\min\{\sigma, \ell\} + \frac{1}{2}} \|v\|_{1+\sigma,K}, $$

$$ \|\nabla v - \Pi_h(\nabla v)\|_{0,\partial K} \le C h_K^{\min\{\sigma, \ell+1\} - \frac{1}{2}} \|v\|_{1+\sigma,K}, $$

with a constant $C$ that is independent of the local mesh size $h_K$ and depends only on the shape-regularity of the mesh, the approximation order $\ell$, the dimension $d$, and the regularity exponent $\sigma$.

As a consequence of the approximation properties in Lemma 4.6, we obtain the following results.

LEMMA 4.7. Let $u \in H^{1+\sigma}(\Omega)$, $\sigma > \frac{1}{2}$. Then the following hold:

(i) We have

$$ \|u - \Pi_h u\|_h \le C_A h^{\min\{\sigma, \ell\}} \|u\|_{1+\sigma}, $$

with a constant $C_A$ that is independent of the mesh size and depends only on $\alpha$, the constant $\kappa$ in (3.1), the bounds in (2.5), and the constants in Lemma 4.6.

(ii) For $v \in V(h)$, the form $r_h(u; v)$ in (4.7) can be bounded by

$$ |r_h(u; v)| \le C_R h^{\min\{\sigma, \ell\}} \left( \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2 \right)^{\frac{1}{2}} \|u\|_{1+\sigma}, $$

with a constant $C_R$ independent of $h$, which depends only on $\alpha$, the bounds in (2.5), and the constants in Lemma 4.6.

*Proof.* The estimate in (i) is an immediate consequence of Lemma 4.6, the definition of $\mathbf{a}$, and the bounded variation property (3.1). To show the bound in (ii), we apply the Cauchy–Schwarz inequality and obtain

$$ \begin{aligned} |r_h(u; v)| & \le \left( \sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}[v]^2 ds \right)^{\frac{1}{2}} \left( \sum_{F \in \mathcal{F}_h} \int_F \mathbf{a}^{-1} \left| \{c \nabla u - c \Pi_h(\nabla u)\} \right|^2 ds \right)^{\frac{1}{2}} \\ & \le \alpha^{-\frac{1}{2}} c_{\star}^{-\frac{1}{2}} c^{\star} \left( \sum_{F \in \mathcal{F}_h} \|\mathbf{a}^{\frac{1}{2}}[v]\|_{0,F}^2 \right)^{\frac{1}{2}} \left( \sum_{K \in \mathcal{T}_h} h_K \|\nabla u - \Pi_h(\nabla u)\|_{0,\partial K}^2 \right)^{\frac{1}{2}}. \end{aligned} $$

Applying the approximation properties in Lemma 4.6 completes the proof. $\square$
**4.3. Proof of Theorem 4.1.** We are now ready to complete the proof of Theorem 4.1. We begin by proving the following auxiliary result.

LEMMA 4.8. Let the analytical solution $u$ of (2.1)–(2.4) satisfy

$$u \in L^{\infty}(J; H^{1+\sigma}(\Omega)), \quad u_t \in L^{\infty}(J; H^{1+\sigma}(\Omega))$$

for $\sigma > \frac{1}{2}$. Let $v \in C^0(\bar{J}; V(h))$ and $v_t \in L^1(J; V(h))$. Then we have

$$\int_J |r_h(u; v_t)| dt \le C_R h^{\min\{\sigma, \ell\}} \|v\|_{L^\infty(J; V(h))} \\ \cdot \left[ 2 \|u\|_{L^\infty(J; H^{1+\sigma}(\Omega))} + T \|u_t\|_{L^\infty(J; H^{1+\sigma}(\Omega))} \right],$$

where $C_R$ is the constant from the bound (ii) in Lemma 4.7.

*Proof.* From the definition of $r_h$ in (4.7) and integration by parts, we obtain

$$\begin{align*}
\int_J r_h(u; v_t) dt &= \int_J \sum_{F \in \mathcal{F}_h} \int_F [[v_t]] \cdot \{\{c \nabla u - c \Pi_h(\nabla u)\}\} dA dt \\
&= - \int_J \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \nabla u_t - c \Pi_h(\nabla u_t)\}\} dA dt \\
&\quad + \left[ \sum_{F \in \mathcal{F}_h} \int_F [[v]] \cdot \{\{c \nabla u - c \Pi_h(\nabla u)\}\} dA \right]_{t=0}^{t=T} \\
&= - \int_J r_h(u_t; v) dt + \left[ r_h(u; v) \right]_{t=0}^{t=T}.
\end{align*}$$

Lemma 4.7 then implies the two estimates

$$\left| \int_J r_h(u_t; v) dt \right| \le C_R h^{\min\{\sigma, \ell\}} T \|v\|_{L^\infty(J; V(h))} \|u_t\|_{L^\infty(J; H^{1+\sigma}(\Omega))}$$

and

$$\left| \left[ r_h(u; v) \right]_{t=0}^{t=T} \right| \le 2C_R h^{\min\{\sigma, \ell\}} \|v\|_{L^\infty(J; V(h))} \|u\|_{L^\infty(J; H^{1+\sigma}(\Omega))},$$

which concludes the proof of the lemma. □

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

To complete the proof of Theorem 4.1, we now set $e = u - u^h$ and recall that $\Pi_h$ is the $L^2$-projection onto $V^h$. Because of (2.8), we have

$$e \in C^0(\bar{J}; V(h)) \cap C^1(\bar{J}; L^2(\Omega)).$$

Next, we use the symmetry of $\tilde{a}_h$ and the error equation in Lemma 4.5 to obtain

$$\begin{align}
\frac{1}{2} \frac{d}{dt} [\|e_t\|_0^2 + \tilde{a}_h(e, e)] &= (e_{tt}, e_t) + \tilde{a}_h(e, e_t) \nonumber \\
&= (e_{tt}, (u - \Pi_h u)_t) + \tilde{a}_h(e, (u - \Pi_h u)_t) \nonumber \\
&\quad + r_h(u; (\Pi_h u - u^h)_t). \tag{4.8}
\end{align}$$

We fix $s \in J$ and integrate (4.8) over the time interval $(0, s)$. This yields

$$\begin{align*}
\frac{1}{2} \|e_t(s)\|_0^2 + \frac{1}{2} \tilde{a}_h(e(s), e(s)) &= \frac{1}{2} \|e_t(0)\|_0^2 + \frac{1}{2} \tilde{a}_h(e(0), e(0)) \\
&\quad + \int_0^s (e_{tt}, (u - \Pi_h u)_t) dt + \int_0^s \tilde{a}_h(e, (u - \Pi_h u)_t) dt \\
&\quad + \int_0^s r_h(u; (\Pi_h u - u^h)_t) dt.
\end{align*}$$
Integration by parts of the third term on the right-hand side yields

$$ \int_{0}^{s} (e_{tt}, (u - \Pi_{h}u)_{t}) dt = - \int_{0}^{s} (e_{t}, (u - \Pi_{h}u)_{tt}) dt + \left[ (e_{t}, (u - \Pi_{h}u)_{t}) \right]_{t=0}^{t=s} $$

From the stability properties of $\tilde{a}_h$ in Lemma 4.4 and standard Hölder's inequalities, we conclude that

$$ \begin{aligned} \frac{1}{2} \|e_t(s)\|_0^2 + \frac{1}{2} C_{\text{coer}} \|e(s)\|_h^2 &\le \frac{1}{2} \|e_t(0)\|_0^2 + \frac{1}{2} C_{\text{cont}} \|e(0)\|_h^2 \\ &\quad + \|e_t\|_{L^\infty(J;L^2(\Omega))} \left( \|(u - \Pi_h u)_{tt}\|_{L^1(J;L^2(\Omega))} + 2\|(u - \Pi_h u)_t\|_{L^\infty(J;L^2(\Omega))} \right) \\ &\quad + C_{\text{cont}} T \|e\|_{L^\infty(J;V(h))} \|(u - \Pi_h u)_t\|_{L^\infty(J;V(h))} \\ &\quad + \left| \int_J r_h(u; (\Pi_h u - u^h)_t) dt \right|. \end{aligned} $$

Since this inequality holds for any $s \in J$, it also holds for the maximum over $J$, that is

$$ \|e_t\|_{L^\infty(J;L^2(\Omega))}^2 + C_{\text{coer}} \|e\|_{L^\infty(J;V(h))}^2 \le \|e_t(0)\|_0^2 + C_{\text{cont}} \|e(0)\|_h^2 + T_1 + T_2 + T_3, $$

with

$$ T_1 = 2 \|e_t\|_{L^\infty(J;L^2(\Omega))} \left( \|(u - \Pi_h u)_{tt}\|_{L^1(J;L^2(\Omega))} + 2 \|(u - \Pi_h u)_t\|_{L^\infty(J;L^2(\Omega))} \right), $$

$$ T_2 = 2C_{\text{cont}} T \|e\|_{L^\infty(J;V(h))} \|(u - \Pi_h u)_t\|_{L^\infty(J;V(h))}, $$

$$ T_3 = 2 \left| \int_J r_h(u; (\Pi_h u - u^h)_t) dt \right|. $$

Using the geometric-arithmetic mean inequality $|ab| \le \frac{1}{2\epsilon}a^2 + \frac{\epsilon}{2}b^2$, valid for any $\epsilon > 0$, and the approximation results in Lemma 4.6, we conclude that

$$ \begin{aligned} T_1 &\le \frac{1}{2} \|e_t\|_{L^\infty(J;L^2(\Omega))}^2 + 2 \left( \|(u - \Pi_h u)_{tt}\|_{L^1(J;L^2(\Omega))} + 2 \|(u - \Pi_h u)_t\|_{L^\infty(J;L^2(\Omega))} \right)^2 \\ &\le \frac{1}{2} \|e_t\|_{L^\infty(J;L^2(\Omega))}^2 + 4 \|(u - \Pi_h u)_{tt}\|_{L^1(J;L^2(\Omega))}^2 + 16 \|(u - \Pi_h u)_t\|_{L^\infty(J;L^2(\Omega))}^2, \\ &\le \frac{1}{2} \|e_t\|_{L^\infty(J;L^2(\Omega))}^2 + C h^{2 \min\{\sigma, \ell\}} \left( \|u_{tt}\|_{L^1(J;H^\sigma(\Omega))}^2 + h^2 \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))}^2 \right), \end{aligned} $$

with a constant $C$ that depends only on the constants in Lemma 4.6. Similarly,

$$ \begin{aligned} T_2 &\le \frac{1}{4} C_{\text{coer}} \|e\|_{L^\infty(J;V(h))}^2 + 4 \frac{C_{\text{cont}}^2}{C_{\text{coer}}} T^2 \|(u - \Pi_h u)_t\|_{L^\infty(J;V(h))}^2 \\ &\le \frac{1}{4} C_{\text{coer}} \|e\|_{L^\infty(J;V(h))}^2 + T^2 C h^{2 \min\{\sigma, \ell\}} \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))}^2, \end{aligned} $$

where the constant $C$ depends on $C_{\text{coer}}$, $C_{\text{cont}}$, and the constant $C_A$ in Lemma 4.7. It remains to bound the term $T_3$. To do so, we use Lemma 4.8 to obtain

$$ T_3 \le 2C_R \mathcal{R} h^{\min\{\sigma, \ell\}} \|\Pi_h u - u^h\|_{L^\infty(J;V(h))}, $$

with

$$ \mathcal{R} := \left[ 2 \|u\|_{L^\infty(J;H^{1+\sigma}(\Omega))} + T \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))} \right]. $$
The triangle inequality, the geometric-arithmetic mean, and the approximation properties of $\Pi_h$ in Lemma 4.7 then yield

$$ 
\begin{aligned}
T_3 &\le 2C_R \mathcal{R} h^{\min\{\sigma, \ell\}} \left[ \|e\|_{L^\infty(J;V(h))} + \|u - \Pi_h u\|_{L^\infty(J;V(h))} \right] \\
&\le \frac{1}{4} C_{\text{coer}} \|e\|_{L^\infty(J;V(h))}^2 + C h^{2\min\{\sigma, \ell\}} \left[ \|u\|_{L^\infty(J;H^{1+\sigma}(\Omega))}^2 + \mathcal{R}^2 \right],
\end{aligned}
$$

with a constant $C$ that depends only on $C_{\text{coer}}$, $C_R$, and $C_A$. Combining the above estimates for $T_1$, $T_2$, and $T_3$ then shows that

$$ 
\begin{aligned}
\frac{1}{2} \|e_t\|_{L^\infty(J;L^2(\Omega))}^2 + \frac{1}{2} C_{\text{coer}} \|e\|_{L^\infty(J;V(h))}^2 &\le \|e_t(0)\|_0^2 + C_{\text{cont}} \|e(0)\|_h^2 \\
&\quad + C h^{2\min\{\sigma, \ell\}} \left[ \|u_{tt}\|_{L^1(J;H^\sigma(\Omega))}^2 + T^2 \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))}^2 + \|u\|_{L^\infty(J;H^{1+\sigma}(\Omega))}^2 \right],
\end{aligned}
$$

with a constant that is independent of $T$ and the mesh size. This concludes the proof of Theorem 4.1.

**4.4. Proof of Theorem 4.2.** To prove the error estimate in Theorem 4.2, we first establish the following variant of [3, Lemma 2.1].

LEMMA 4.9. For $u \in H^{1+\sigma}(\Omega)$ with $\sigma > \frac{1}{2}$, let $w^h \in V^h$ be the solution of

$$ \tilde{a}_h(w^h, v) = \tilde{a}_h(u, v) - r_h(u; v) \quad \forall v \in V^h. $$

Then, we have

$$ \|u - w^h\|_h \le C_E h^{\min\{\sigma, \ell\}} \|u\|_{1+\sigma}, $$

with a constant $C_E$ that is independent of $h$ and depends only on $C_{\text{coer}}$, $C_{\text{cont}}$ in Lemma 4.4 and $C_A$, $C_R$ in Lemma 4.7.

Moreover, if the elliptic regularity defined in (4.1) and (4.2) holds, we have the $L^2$-bound

$$ \|u - w^h\|_0 \le C_L h^{\min\{\sigma, \ell\}+1} \|u\|_{1+\sigma}, $$

with a constant $C_L$ that is independent of $h$ and depends only on the stability constant $C_S$ in (4.2); $C_{\text{coer}}$, $C_{\text{cont}}$ in Lemma 4.4; and $C_A$, $C_R$ in Lemma 4.7.

*Proof.* We first remark that the approximation $w^h$ is well defined because of the stability properties in Lemma 4.4 and the estimates in Lemma 4.7. To prove the estimate for $\|u - w^h\|_h$, we first use the triangle inequality,

$$ (4.9) \qquad \|u - w^h\|_h \le \|u - \Pi_h u\|_h + \|\Pi_h u - w^h\|_h. $$

From the approximation properties of $\Pi_h$ in Lemma 4.7, we immediately infer that

$$ \|u - \Pi_h u\|_h \le C_A h^{\min\{\sigma, \ell\}} \|u\|_{1+\sigma}. $$

It remains to bound $\|\Pi_h u - w^h\|_h$. From the coercivity and continuity of $\tilde{a}_h$ in Lemma 4.4, the definition of $w^h$, and the bound in Lemma 4.7, we conclude that

$$ 
\begin{aligned}
C_{\text{coer}} \|\Pi_h u - w^h\|_h^2 &\le \tilde{a}_h (\Pi_h u - w^h, \Pi_h u - w^h) \\
&= \tilde{a}_h (\Pi_h u - u, \Pi_h u - w^h) + \tilde{a}_h (u - w^h, \Pi_h u - w^h) \\
&= \tilde{a}_h (\Pi_h u - u, \Pi_h u - w^h) + r_h (u; \Pi_h u - w^h) \\
&\le C_{\text{cont}} \|\Pi_h u - u\|_h \|\Pi_h u - w^h\|_h + C_R h^{\min\{\sigma, \ell\}} \|u\|_{1+\sigma} \|\Pi_h u - w^h\|_h.
\end{aligned}
$$
Thus,

$$ \|\Pi_h u - w^h\|_h \le \left( \frac{C_{\text{cont}} C_A + C_R}{C_{\text{coer}}} \right) h^{\min\{\sigma, \ell\}} \|u\|_{1+\sigma}, $$

which proves the bound for $\|u - w^h\|_h$.

We shall now prove the $L^2$-bound. To do so, let $z \in H_0^1(\Omega)$ be the solution of
$$ (4.10) \qquad -\nabla \cdot (c \nabla z) = u - w^h \quad \text{in } \Omega, \quad z = 0 \quad \text{on } \Gamma. $$

Then, the elliptic regularity assumption in (4.1) and (4.2) implies that

$$ (4.11) \qquad z \in H^2(\Omega), \qquad \|z\|_2 \le C_S \|u - w^h\|_0. $$

Next, we multiply (4.10) by $u - w^h$ and integrate the resulting expression by parts. Since $c \nabla z$ has continuous normal components across all interior faces, we have

$$ \begin{aligned} \|u - w^h\|_0^2 &= \sum_{K \in \mathcal{T}_h} \left[ \int_K c \nabla z \cdot \nabla (u - w^h) \,dx - \int_{\partial K} c \nabla z \cdot \mathbf{n}_K (u - w^h) \,dA \right] \\ &= \sum_{K \in \mathcal{T}_h} \int_K c \nabla z \cdot \nabla (u - w^h) \,dx - \sum_{F \in \mathcal{F}_h} \int_F \{c \nabla z\} \cdot [u - w^h] \,dA, \end{aligned} $$

with $\mathbf{n}_K$ denoting the unit outward normal on $\partial K$. By definition of $\tilde{a}_h$ and $r_h$, we immediately find that

$$ \|u - w^h\|_0^2 = \tilde{a}_h(z, u - w^h) - r_h(z; u - w^h). $$

From the symmetry of $\tilde{a}_h$, the definition of $w^h$, and the fact that $[z] = 0$ on all faces, we conclude that

$$ (4.12) \qquad \begin{aligned} \|u - w^h\|_0^2 &= \tilde{a}_h(u - w^h, z - \Pi_h z) - r_h(u; z - \Pi_h z) - r_h(z; u - w^h) \\ &=: T_1 + T_2 + T_3. \end{aligned} $$

We shall now derive upper bounds for each individual term $T_1$, $T_2$, and $T_3$ in (4.12).

To estimate the term $T_1$, we use the continuity of $\tilde{a}_h$, the approximation result in Lemma 4.7 with $\sigma = 1$, and the bound in (4.11). Thus,

$$ \begin{aligned} T_1 &\le C_{\text{cont}} \|u - w^h\|_h \|z - \Pi_h z\|_h \\ &\le C_{\text{cont}} C_A h \|u - w^h\|_h \|z\|_2 \\ &\le C_{\text{cont}} C_A C_S h \|u - w^h\|_h \|u - w^h\|_0. \end{aligned} $$

By using Lemma 4.7 and the stability bound in (4.11), we can estimate $T_2$ by

$$ \begin{aligned} T_2 &\le C_R h^{\min\{\sigma, \ell\}} \|z - \Pi_h z\|_h \|u\|_{1+\sigma} \\ &\le C_R C_A h^{\min\{\sigma, \ell\}+1} \|z\|_2 \|u\|_{1+\sigma} \\ &\le C_R C_A C_S h^{\min\{\sigma, \ell\}+1} \|u - w^h\|_0 \|u\|_{1+\sigma}. \end{aligned} $$

Similarly,

$$ T_3 \le C_R h \|z\|_2 \|u - w^h\|_h \le C_R C_S h \|u - w^h\|_0 \|u - w^h\|_h. $$
The use of these bounds for $T_1$, $T_2$, and $T_3$ in (4.12) then leads to

$$ \|u - w^h\|_0 \le Ch\|u - w^h\|_h + Ch^{\min\{\sigma, \ell\}+1}\|u\|_{1+\sigma}, $$

which completes the proof of the lemma, since $\|u - w^h\|_h \le Ch^{\min\{\sigma, \ell\}}\|u\|_{1+\sigma}$. $\square$

Now, let $u$ be defined by the exact solution of (2.1)–(2.4). We may define $w^h(t, \cdot) \in V^h$ almost everywhere in $J$ by

$$ (4.13) \qquad \tilde{a}_h(w^h(t, \cdot), v) = \tilde{a}_h(u(t, \cdot), v) - r_h(u(t, \cdot); v) \quad \forall v \in V^h. $$

If $u \in L^\infty(J; H^{1+\sigma}(\Omega))$, it can be readily seen that $w^h \in L^\infty(J; V(h))$. Moreover, if we also have $u_t \in L^\infty(J; H^{1+\sigma}(\Omega))$, then $w_t^h \in L^\infty(J; V(h))$ and

$$ \tilde{a}_h(w_t^h, v) = \tilde{a}_h(u_t, v) - r_h(u_t; v), \quad v \in V^h \text{ a.e. in } J, $$

as well as

$$ \tilde{a}_h(w^h(0), v) = \tilde{a}_h(u_0, v) - r_h(u_0; v), \quad v \in V^h. $$

Therefore Lemma 4.9 immediately implies the following estimates.

**LEMMA 4.10.** Let $w^h$ be defined by (4.13). Under the regularity assumptions of Theorem 4.2, we have

$$ \|(u - w^h)_t\|_{L^\infty(J; V(h))} \le C_E h^{\min\{\sigma, \ell\}} \|u_t\|_{L^\infty(J; H^{1+\sigma}(\Omega))}, $$

$$ \|(u - w^h)(0)\|_h \le C_E h^{\min\{\sigma, \ell\}} \|u_0\|_{1+\sigma}. $$

Moreover, if elliptic regularity as defined in (4.1) and (4.2) holds, we have the $L^2$-bounds

$$ \|(u - w^h)_t\|_{L^\infty(J; L^2(\Omega))} \le C_L h^{\min\{\sigma, \ell\}+1} \|u_t\|_{L^\infty(J; H^{1+\sigma}(\Omega))}, $$

$$ \|(u - w^h)(0)\|_0 \le C_L h^{\min\{\sigma, \ell\}+1} \|u_0\|_{1+\sigma}. $$

The constants $C_E$ and $C_L$ are as in Lemma 4.9.

To complete the proof of Theorem 4.2, let $w^h \in L^\infty(J; V(h))$ be defined by (4.13) and consider

$$ (4.14) \qquad \|e\|_{L^\infty(J; L^2(\Omega))}^2 \le 2\|u - w^h\|_{L^\infty(J; L^2(\Omega))}^2 + 2\|w^h - u^h\|_{L^\infty(J; L^2(\Omega))}^2. $$

The first term can be estimated from the $L^2$-bounds in Lemma 4.9. We shall now derive an estimate for the second term. First, we fix $v \in L^\infty(J; V^h)$ and assume that $v_t \in L^\infty(J; V^h)$. From the definition of $w^h$ in (4.13) and the error equation in Lemma 4.5, we have

$$ \begin{aligned} \left\| (u^h - w^h)_{tt}, v \right\| + \tilde{a}_h(u^h - w^h, v) &= (u_{tt}^h, v) + \tilde{a}_h(u^h, v) - \tilde{a}_h(w^h, v) - (w_{tt}^h, v) \\ &= (u_{tt}^h, v) + \tilde{a}_h(u^h - u, v) + r_h(u; v) - (w_{tt}^h, v) \\ &= (u_{tt}^h, v) - (w_{tt}^h, v). \end{aligned} $$

We rewrite this identity as

$$ \frac{d}{dt}((u^h - w^h)_t, v) - ((u^h - w^h)_t, v_t) + \tilde{a}_h(u^h - w^h, v) = \frac{d}{dt}((u - w^h)_t, v) - ((u - w^h)_t, v_t), $$
which yields

$$ (4.15) \quad -((u^h - w^h)_t, v_t) + \tilde{a}_h(u^h - w^h, v) = \frac{d}{dt}((u - u^h)_t, v) - ((u - w^h)_t, v_t). $$

Let $\tau \in (0, T]$ be fixed, and consider the function

$$ \hat{v}(t, \cdot) = \int_t^\tau (u^h - w^h)(s, \cdot) ds, \quad t \in \bar{J}. $$

Note that

$$ \hat{v}(\tau, \cdot) = 0, \quad \hat{v}_t(t, \cdot) = -(u^h - w^h)(t, \cdot) \quad \text{a.e. } t \in \bar{J}. $$

Next, choose $v = \hat{v}$ in (4.15) which yields

$$ ((u^h - w^h)_t, u^h - w^h) - \tilde{a}_h(\hat{v}_t, \hat{v}) = \frac{d}{dt}((u - u^h)_t, \hat{v}) + ((u - w^h)_t, u^h - w^h). $$

Since the DG form $\tilde{a}_h(\cdot, \cdot)$ is symmetric, we obtain

$$ \frac{1}{2} \frac{d}{dt} \|u^h - w^h\|_0^2 - \frac{1}{2} \frac{d}{dt} \tilde{a}_h(\hat{v}, \hat{v}) = \frac{d}{dt}((u - u^h)_t, \hat{v}) + ((u - w^h)_t, u^h - w^h). $$

Integration over $(0, \tau)$ and using that $\hat{v}(\tau, \cdot) = 0$ then yield

$$ \begin{aligned} & \|(u^h - w^h)(\tau)\|_0^2 - \|(u^h - w^h)(0)\|_0^2 + \tilde{a}_h(\hat{v}(0), \hat{v}(0)) \\ & \quad = -2((u - u^h)_t(0), \hat{v}(0)) + 2 \int_0^\tau ((u - w^h)_t, u^h - w^h) dt. \end{aligned} \quad (4.16) $$

Since $u_t(0) = v_0$, $u_t^h(0) = \Pi_h v_0$, and $\hat{v}(0)$ belongs to $V^h$, we conclude that

$$ ((u - u^h)_t(0), \hat{v}(0)) = (v_0 - \Pi_h v_0, \hat{v}(0)) = 0. $$

Hence, the first term on the right-hand side of (4.16) vanishes. Moreover, the coercivity of the form $\tilde{a}_h$ in Lemma 4.4 ensures that $\tilde{a}_h(\hat{v}(0), \hat{v}(0)) \ge 0$. This leads to the inequality

$$ (4.17) \quad \|(u^h - w^h)(\tau)\|_0^2 \le \|(u^h - w^h)(0)\|_0^2 + 2 \int_0^\tau \|(u - w^h)_t\|_0 \|u^h - w^h\|_0 dt. $$

By using the Cauchy–Schwarz inequality and the geometric-arithmetic mean inequality, we obtain

$$ \begin{aligned} & 2 \int_0^\tau \|(u - w^h)_t\|_0 \|u^h - w^h\|_0 dt \le 2T \|(u - w^h)_t\|_{L^\infty(J;L^2(\Omega))} \|u^h - w^h\|_{L^\infty(J;L^2(\Omega))} \\ & \quad \le \frac{1}{2} \|u^h - w^h\|_{L^\infty(J;L^2(\Omega))}^2 + 2T^2 \|(u - w^h)_t\|_{L^\infty(J;L^2(\Omega))}^2. \end{aligned} $$

Because this upper bound is independent of $\tau$, it also holds for the supremum over $\tau \in J$, which yields the estimate

$$ \begin{aligned} & \frac{1}{2} \|u^h - w^h\|_{L^\infty(J;L^2(\Omega))}^2 \le \|(u^h - w^h)(0)\|_0^2 + 2T^2 \|(u - w^h)_t\|_{L^\infty(J;L^2(\Omega))}^2 \\ & \quad \le 2 \|(u^h - u)(0)\|_0^2 + 2 \|(u - w^h)(0)\|_0^2 + 2T^2 \|(u - w^h)_t\|_{L^\infty(J;L^2(\Omega))}^2. \end{aligned} $$
Next, we use this estimate in (4.14) to obtain

$$ \|e\|_{L^{\infty}(J;L^2(\Omega))}^2 \le 2\|u - w^h\|_{L^{\infty}(J;L^2(\Omega))}^2 \\ + 8 \|u_0 - \Pi_h u_0\|_0^2 + 8\|u_0 - w^h(0)\|_0^2 + 8T^2 \|(u - w^h)_t\|_{L^{\infty}(J;L^2(\Omega))}^2. $$

From the $L^2$-approximation properties in Lemma 4.6, Lemma 4.9, and Lemma 4.10, we finally conclude that

$$ \|e\|_{L^{\infty}(J;L^2(\Omega))}^2 \le h^{2 \min\{\sigma, \ell\} + 2} \left[ \max\{8C, 8C_L^2\} \|u_0\|_{1+\sigma}^2 \right. \\ \left. \qquad + 2C_L^2 \|u\|_{L^{\infty}(J;H^{1+\sigma}(\Omega))}^2 + 8C_L T^2 \|u_t\|_{L^{\infty}(J;H^{1+\sigma}(\Omega))}^2 \right]. $$

Here, $C$ is the constant from Lemma 4.6. This completes the proof of Theorem 4.2.

**5. Numerical results.** We shall now present a series of numerical experiments which verify the sharpness of the theoretical error bounds stated in Theorems 4.1 and 4.2. Furthermore, we shall demonstrate the robustness and flexibility of our DG method by propagating a pulse through an inhomogeneous medium with discontinuity on a finite element mesh with hanging nodes.

To obtain a fully discrete discretization of the wave equation, we choose to augment our DG spatial discretization with the second-order Newmark scheme in time; see, e.g., [19, sections 8.5–8.7]. The resulting scheme has been implemented using the general purpose finite element library deal.II,<fnref n="1" /> which provides powerful $C^{++}$ classes to handle the finite element mesh and the degrees of freedom, and to solve the linear system of equations; see [5, 4]. In all our examples, the DG stabilization parameter is set to $\alpha = 20$.

**5.1. Time discretization.** The discretization of (2.1)–(2.4) in space by the DG method (3.3)–(3.5) leads to the linear second-order system of ordinary differential equations

$$ (5.1) \qquad \mathbf{M} \ddot{\mathbf{u}}^h(t) + \mathbf{A} \mathbf{u}^h(t) = \mathbf{f}^h(t), \quad t \in J, $$

with initial conditions

$$ (5.2) \qquad \mathbf{M} \mathbf{u}^h(0) = \mathbf{u}_0^h, \quad \mathbf{M} \dot{\mathbf{u}}^h(0) = \mathbf{v}_0^h. $$

Here, $\mathbf{M}$ denotes the mass matrix and $\mathbf{A}$ the stiffness matrix. To discretize (5.1) in time, we employ the Newmark time-stepping scheme; see, e.g., [19]. We let $k$ denote the time step and set $t_n = n \cdot k$. Then the Newmark method consists in finding approximations $\{\mathbf{u}_n^h\}_n$ to $\mathbf{u}^h(t_n)$ such that

$$ (5.3) \qquad (\mathbf{M} + k^2 \beta \mathbf{A}) \mathbf{u}_1^h = \left[ \mathbf{M} - k^2 \left( \frac{1}{2} - \beta \right) \mathbf{A} \right] \mathbf{u}_0^h + k \mathbf{M} \mathbf{v}_0^h + k^2 \left[ \beta \mathbf{f}_1^h + \left( \frac{1}{2} - \beta \right) \mathbf{f}_0^h \right] $$

and

$$ (\mathbf{M} + k^2 \beta \mathbf{A}) \mathbf{u}_{n+1}^h = \left[ 2\mathbf{M} - k^2 \left( \frac{1}{2} - 2\beta + \gamma \right) \mathbf{A} \right] \mathbf{u}_n^h - \left[ \mathbf{M} + k^2 \left( \frac{1}{2} + \beta - \gamma \right) \mathbf{A} \right] \mathbf{u}_{n-1}^h \\ (5.4) \qquad + k^2 \left[ \beta \mathbf{f}_{n+1}^h + \left( \frac{1}{2} - 2\beta + \gamma \right) \mathbf{f}_n + \left( \frac{1}{2} + \beta - \gamma \right) \mathbf{f}_{n-1} \right] $$

<footnote n="1">¹See www.dealii.org.</footnote>
for $n = 1, \dots, N - 1$. Here, $f_n := f(t_n)$, while $\beta \ge 0$ and $\gamma \ge 1/2$ are free parameters that still can be chosen. We recall that for $\gamma = 1/2$ the Newmark scheme is second-order accurate in time, whereas it is only first-order accurate for $\gamma > 1/2$. For $\beta = 0$, the Newmark scheme (5.3)–(5.4) requires at each time step the solution of a linear system with the mass matrix $M$. However, because individual elements decouple, $M$ is block-diagonal with a block size equal to the number of degrees of freedom per element. It can be inverted at very low computational cost, and the scheme is essentially fully explicit. In fact, if the basis functions are chosen mutually orthogonal, $M$ reduces to the identity; see [8] and the references therein. Then, with $\gamma = 1/2$, the explicit Newmark method corresponds to the standard leap-frog scheme.

For $\beta > 0$, the resulting scheme is implicit and involves the solution of a linear system with the symmetric positive definite stiffness matrix $A$ at each time step. We finally note that the second-order Newmark scheme with $\gamma = 1/2$ is unconditionally stable for $\beta \ge 1/4$, whereas for $1/4 > \beta \ge 0$ the time step $k$ has to be restricted by a CFL condition. In the case $\beta = 0$, that condition is $k^2\lambda_{\max}(A) \le 4(1 - \varepsilon)$, $\varepsilon \in (0, 1)$, where $\lambda_{\max}(A)$ is the maximal eigenvalue of the DG stiffness matrix $A$ (which is of the order $O(h^{-2})$ and also depends on $\alpha$).

In all our tests, we will employ the explicit second-order Newmark scheme, setting $\gamma = 1/2$ and $\beta = 0$ in (5.3)–(5.4).

**5.2. Example 1: Smooth solution.** First, we consider the two-dimensional wave equation (2.1)–(2.4) in $J \times \Omega = (0, 1) \times (0, 1)^2$, with $c \equiv 1$ and data $f, u_0$, and $v_0$ chosen such that the analytical solution is given by

$$ (5.5) \qquad u(x_1, x_2, t) = t^2 \sin(\pi x_1) \sin(\pi x_2). $$

This solution is arbitrarily smooth so that all our theoretical regularity assumptions are satisfied. We discretize this problem using the polynomial spaces $Q^\ell(K)$, $\ell = 1, 2, 3$, on a sequence $\{T_h\}_{i \ge 1}$ of square meshes of size $h_i = 2^{-i}$. With increasing polynomial degree $\ell$ and decreasing mesh size $h_i$, smaller time steps $k_i$ are necessary to ensure stability. We found that the choice $k_i = h_i/20$ provides a stable time discretization on every mesh. Because our numerical scheme is second-order accurate in time, the time integration of (5.5) is exact so that the spatial error is the only error component in the discrete solution.

In Figure 5.1 we show the relative errors at time $T = 1$ in the energy norm and in the $L^2$-norm, as we decrease the mesh size $h_i$. The numerical results corroborate with the expected theoretical rates of $O(h^\ell)$ for the energy norm and of $O(h^{\ell+1})$ for the $L^2$-norm; see Theorems 4.1 and 4.2.

Next, we modify the data so that the analytical solution $u$ is given by

$$ (5.6) \qquad u(x_1, x_2, t) = \sin(t^2) \sin(\pi x_1) \sin(\pi x_2). $$

Although $u$ remains arbitrarily smooth, it is no longer integrated exactly in time by (5.3)–(5.4). Since the Newmark scheme is only second-order accurate, we repeat the above experiment only for the lowest order spatial discretization, $\ell = 1$. Again, we set $k_i = h_i/20$. In Figure 5.2, the relative errors for the fully discrete approximation of (5.6) show convergence rates of order $h$ in the energy norm and order $h^2$ in the $L^2$-norm, thereby confirming the theoretical estimates of Theorems 4.1 and 4.2.

**5.3. Example 2: Singular solution.** Here, we consider the two-dimensional wave equation (2.1)–(2.4) on the L-shaped domain $\Omega = (-1, 1)^2 \setminus [0, 1)^2$. We set $c = 1$
<figure><img src="image_3.png" /><figcaption>FIG. 5.1. Example 1: Convergence of the error at time T = 1 in the energy norm and the L²-norm for ℓ = 1, 2, 3.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 5.2. Example 1: Convergence of the error at time T = 1 in the energy norm and the L²-norm for ℓ = 1.</figcaption></figure>

everywhere and choose the data $f, u_0$, and $v_0$ such that the analytical solution $u$ is given by

$$ (5.7) \qquad u(r, \phi, t) = t^2 r^{2/3} \sin(2/3 \phi) $$

in polar coordinates $(r, \phi)$. Although $u$ is smooth in time (and can even be integrated exactly in time), it has a spatial singularity at the origin, such that $u \in C^\infty(\bar{J}; H^{5/3}(\Omega))$. Hence, this example is well suited to establishing the sharpness of the regularity assumptions in our theoretical results. Since $u$ is inhomogeneous at the boundary of $\Omega$, we need to impose inhomogeneous Dirichlet conditions within our DG discretization. We do so in straightforward fashion by modifying the semidiscrete formulation as follows: find $u^h(t, \cdot) : J \to V^h$ such that

$$ (5.8) \qquad (u_{tt}^h, v) + a_h(u^h, v) = (f, v) + \sum_{F \in \mathcal{F}_h^B} \int_F g (av - c \nabla v \cdot \mathbf{n}) dA. $$

Here, $g$ is the boundary data and $\mathbf{n}$ is the outward unit normal vector on $\partial\Omega$.

We discretize (5.8) by using bilinear polynomials ($\ell = 1$) on the same sequence of meshes as before. Again, we set $k_i = h_i/20$ and integrate the problem up to $T = 1$. For the analytical solution $u$ in (5.7), the regularity assumptions in Theorem 4.1 hold
<caption>TABLE 5.1</caption>

*Example 2: Relative errors at time T = 1 in the energy norm and L²-norm, and corresponding numerical convergence rates.*

<table><thead><tr><th rowspan="2">i</th><th rowspan="2">Cells</th><th colspan="2">Energy-error</th><th colspan="2">L²-error</th></tr><tr><th>1.11e-01</th><th>-</th><th>1.61e-02</th><th>-</th></tr></thead><tbody><tr><td>1</td><td>12</td><td>-</td><td>-</td><td>1.61e-02</td><td>-</td></tr><tr><td>2</td><td>48</td><td>7.18e-02</td><td>0.62</td><td>5.96e-03</td><td>1.43</td></tr><tr><td>3</td><td>192</td><td>4.61e-02</td><td>0.64</td><td>2.27e-03</td><td>1.40</td></tr><tr><td>4</td><td>768</td><td>2.94e-02</td><td>0.65</td><td>8.72e-04</td><td>1.38</td></tr><tr><td>5</td><td>3072</td><td>1.87e-02</td><td>0.66</td><td>3.38e-04</td><td>1.37</td></tr><tr><td>6</td><td>12288</td><td>1.18e-02</td><td>0.66</td><td>1.32e-04</td><td>1.36</td></tr></tbody></table>

with $\sigma = 2/3$. Thus, Theorem 4.1 predicts numerical convergence rates of 2/3 in the energy norm, as confirmed by our numerical results in Table 5.1.

As the elliptic regularity assumptions (4.1)-(4.2) from Theorem 4.2 are violated, we do not expect $L^2$-error rates of the order $1+\sigma$ for this problem. Indeed, in Table 5.1 we observe convergence rates close to 4/3. To explain this behavior, let us consider the following weaker elliptic regularity assumption: for any $\lambda \in L^2(\Omega)$ we assume that the solution of the problem

$$ (5.9) \qquad -\nabla \cdot (c\nabla z) = \lambda \quad \text{in } \Omega, \quad z = 0 \quad \text{on } \partial\Omega, $$

belongs to $H^{1+s}(\Omega)$ for a parameter $s \in (1/2, 1]$ and satisfies the stability bound

$$ (5.10) \qquad \|z\|_{1+s} \le C_S \|\lambda\|_0 $$

for a stability constant $C_S$. The results from Lemmas 4.9 and 4.10 can be easily adapted to this case. As a consequence, the $L^2$-bound for $e = u - u^h$ from Theorem 4.2 can then be generalized to this weaker setting as

$$ \|e\|_{L^\infty(J;L^2(\Omega))} \le C h^{\min\{\sigma,\ell\}+s} \left[ \|u_0\|_{1+\sigma} + \|u\|_{L^\infty(J;H^{1+\sigma}(\Omega))} + T \|u_t\|_{L^\infty(J;H^{1+\sigma}(\Omega))} \right]. $$

For the L-shaped domain $\Omega$ and $c \equiv 1$, the (weaker) regularity assumption in (5.9)-(5.10) holds with $s = 2/3$, which underpins the rate $\sigma+s = 4/3$ observed in Table 5.1.

**5.4. Example 3: Inhomogeneous medium.** Finally, we consider (2.1)-(2.4) on the rectangular domain $\Omega = (-1, 2) \times (-1, 1)$, with homogeneous initial and boundary conditions and the piecewise constant material coefficient

$$ c(x_1, x_2) = \begin{cases} 0.1, & x_1 \le 0, \\ 1, & \text{else.} \end{cases} $$

The wave is locally excited until $t = 0.2$ by the source term

$$ f(x_1, x_2, t) = \begin{cases} 1, & 0.2 < x_1 < 0.4 \text{ and } t < 0.2, \\ 0, & \text{else.} \end{cases} $$

We discretize the problem by the DG method (3.3)-(3.5) on a fixed mesh $T_h$ that consists of *nonmatching components*, which are adapted to the discontinuity $c$; see Figure 5.3. The mesh $T_h$ is composed of 9312 nonuniform squares, where the smallest local mesh size is given by $h_{min} \approx 0.016$. The hanging nodes are naturally incorporated in the DG method without any difficulty. Here, the time step $k = 0.002$,
<figure><img src="image_5.png" /><figcaption>FIG. 5.3. Example 3: Domain $\Omega$ with a finite element mesh $T_h$ adapted to the values of the piecewise constant wave speed $\sqrt{c} = 0.1$ (left), $\sqrt{c} = 1$ (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5.4. Example 3: The approximate DG solutions $u_h^n$ shown at times $t_n = 0.2, 0.6, 1.8, 4$ display the behavior of a wave propagating through an inhomogeneous medium with homogeneous Dirichlet boundary.</figcaption></figure>

that is, $k \approx h_{\min}/8$, proved to be sufficiently small to ensure the stability of the explicit Newmark method ($\beta = 0$).

In Figure 5.4, the numerical solution is shown after $n = 100, 300, 900$, and $2000$ time steps. The initial pulse splits into two planar wave fronts propagating in opposite directions to either side of the domain. After $n = 300$ time steps, the left-moving wave hits the much slower medium in the region $x_1 \le 0$, resulting in a much steeper and narrower wave front. Meanwhile, the right-moving wave rapidly reaches the boundary at $x_1 = 2$, where it is reflected and eventually enters the slow medium, too. The discontinuous interface at $x_1 = 0$ generates multiple reflections, which interact with each other at later times.
**6. Conclusion.** We have presented and analyzed the symmetric interior penalty discontinuous Galerkin finite element method (DGFEM) for the numerical solution of the (second-order) scalar wave equation. Taking advantage of the symmetry of the method, we have carried out an a priori error analysis of the semidiscrete method and derived optimal error bounds in the energy norm and, under additional regularity assumptions, optimal error bounds in the $L^2$-norm. Our numerical results confirm the expected convergence rates and demonstrate the versatility of the method. The error analysis of the fully discrete scheme is the subject of ongoing work.

Based on discontinuous finite element spaces, the proposed DG method easily handles elements of various types and shapes, irregular nonmatching grids, and even locally varying polynomial order. As continuity is only weakly enforced across mesh interfaces, domain decomposition techniques immediately apply. Since the resulting mass matrix is essentially diagonal, the method is inherently parallel and leads to fully explicit time integration schemes. Moreover, as the stiffness matrix is symmetric positive definite, the DG method shares the following two important properties with the classical continuous Galerkin approach. First, the semidiscrete formulation conserves (a discrete version of) the energy for all time and therefore is nondissipative. Second, if implicit time integration is used to overcome CFL constraints, the resulting linear system to be solved at each time step will also be symmetric positive definite.

The symmetric interior penalty DGFEM, applied here to the scalar wave equation, can also be utilized for other second-order hyperbolic equations, such as in electromagnetics or elasticity. In fact, our error analysis for the semidiscrete (scalar) case readily extends to the second-order (vector) wave equation for time dependent elastic waves. The same DG approach also extends to Maxwell's equations in second-order form:

$$ \varepsilon \mathbf{u}_{tt} + \sigma \mathbf{u}_t + \nabla \times (\mu^{-1} \nabla \times \mathbf{u}) = \mathbf{f}, \quad \sigma \geq 0. $$

Here, DG discretizations with standard discontinuous piecewise polynomials indeed offer an alternative to edge elements, as typically used in conforming discretizations of Maxwell's equations. The corresponding error analysis, however, is more involved and will be reported elsewhere in the near future.

## REFERENCES

[1] M. AINSWORTH, P. MONK, AND W. MUNIZ, *Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second order wave equation*, J. Sci. Comput., 27 (2006), pp. 5–40.

[2] D. N. ARNOLD, F. BREZZI, B. COCKBURN, AND L. D. MARINI, *Unified analysis of discontinuous Galerkin methods for elliptic problems*, SIAM J. Numer. Anal., 39 (2002), pp. 1749–1779.

[3] G. A. BAKER, *Error estimates for finite element methods for second order hyperbolic equations*, SIAM J. Numer. Anal., 13 (1976), pp. 564–576.

[4] W. BANGERTH, R. HARTMANN, AND G. KANSCHAT, *deal.II Differential Equations Analysis Library, Technical Reference*, IWR, Universität Heidelberg, Heidelberg, Germany, http://www.dealii.org.

[5] W. BANGERTH AND G. KANSCHAT, *Concepts for Object-Oriented Finite Element Software—The deal.II Library*, Report 99-43, Sonderforschungsbereich 3-59, IWR, Universität Heidelberg, Heidelberg, Germany, 1999.

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

[7] B. COCKBURN, *Discontinuous Galerkin methods for convection-dominated problems*, in High-Order Methods for Computational Physics, Lect. Notes Comput. Sci. Eng. 9, T. Barth and H. Deconink, eds., Springer-Verlag, Berlin, 1999, pp. 69–224.
[8] B. COCKBURN, G. E. KARNIADAKIS, AND C.-W. SHU, *The development of discontinuous Galerkin methods*, in *Discontinuous Galerkin Methods: Theory, Computation and Applications*, Lect. Notes Comput. Sci. Eng. 11, B. Cockburn, G. E. Karniadakis, and C.-W. Shu, eds., Springer-Verlag, Berlin, 2000, pp. 3–50.

[9] B. COCKBURN AND C.-W. SHU, *TVB Runge-Kutta local projection discontinuous Galerkin method for conservation laws*. II: *General framework*, Math. Comp., 52 (1989), pp. 411–435.

[10] B. COCKBURN AND C.-W. SHU, *Runge-Kutta discontinuous Galerkin methods for convection-dominated problems*, J. Sci. Comput., 16 (2001), pp. 173–261.

[11] G. COHEN, P. JOLY, J. E. ROBERTS, AND N. TORDJMAN, *Higher order triangular finite elements with mass lumping for the wave equation*, SIAM J. Numer. Anal., 38 (2001), pp. 2047–2078.

[12] R. S. FALK AND G. R. RICHTER, *Explicit finite element methods for symmetric hyperbolic equations*, SIAM J. Numer. Anal., 36 (1999), pp. 935–952.

[13] J. S. HESTHAVEN AND T. WARBURTON, *Nodal high-order methods on unstructured grids*. I. *Time-domain solution of Maxwell's equations*, J. Comput. Phys., 181 (2002), pp. 186–221.

[14] P. HOUSTON, M. JENSEN, AND E. SÜLI, *hp-discontinuous Galerkin finite element methods with least-squares stabilization*, J. Sci. Comput., 17 (2002), pp. 3–25.

[15] T. HUGHES, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Prentice-Hall, Englewood Cliffs, NJ, 1987.

[16] J.-L. LIONS AND E. MAGENES, *Non-Homogeneous Boundary Value Problems and Applications*, Vol. I, Springer-Verlag, New York, 1972.

[17] P. MONK AND G. R. RICHTER, *A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media*, J. Sci. Comput., 22–23 (2005), pp. 443–477.

[18] I. PERUGIA AND D. SCHÖTZAU, *An hp-analysis of the local discontinuous Galerkin method for diffusion problems*, J. Sci. Comput., 17 (2002), pp. 561–571.

[19] P. A. RAVIART AND J.-M. THOMAS, *Introduction à l'Analyse Numérique des Équations aux Dérivées Partielles*, Masson, Paris, 1983.

[20] B. RIVIÈRE AND M. F. WHEELER, *Discontinuous Finite Element Methods for Acoustic and Elastic Wave Problems*, I: *Semidiscrete Error Estimates*, Technical report 01-02, TICAM, University of Texas, Austin, TX, 2001.

[21] B. RIVIÈRE AND M. F. WHEELER, *Discontinuous finite element methods for acoustic and elastic wave problems*, in ICM2002-Beijing Satellite Conference on Scientific Computing, Contemp. Math. 329, AMS, Providence, RI, 2003, pp. 271–282.