Discontinuous Galerkin Finite Element Method for the Wave Equation

. The symmetric interior penalty discontinuous Galerkin ﬁnite element method is presented for the numerical discretization of the second-order wave equation. The resulting stiﬀness matrix is symmetric positive deﬁnite, 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,(cid:2) } ) with respect to the mesh size h , the polynomial degree (cid:2) , 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 (cid:2) +1 ). Numerical results conﬁrm the expected convergence rates and illustrate the versatility of the method.


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 inhomogenous 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.† 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).
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 in Ω, (2.4) where J = (0, T ) is a finite time interval and Ω is a bounded domain in R d , d = 2, 3.For simplicity, we assume that Ω is a polygon (d = 2) or a polyhedron (d = 3).The (known) source term f lies in L 2 (J; L 2 (Ω)), while u 0 ∈ H 1 0 (Ω) and v 0 ∈ L 2 (Ω) are prescribed initial conditions.We assume that the speed of propagation, c(x), is piecewise smooth and satisfies the bounds The standard variational form of (2.1)- (2.4) Here, the time derivatives are understood in a distributional sense, •, • denotes the duality pairing between H −1 (Ω) and H 1 0 (Ω), (•, •) is the inner product in L 2 (Ω), and a(•, •) is the elliptic bilinear form given by 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, 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.

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 T h that partition the domain Ω into disjoint elements {K} such that Ω = ∪ K∈T h 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∈T h h K .We assume that the partition is aligned with the discontinuities of the wave speed √ 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 κ, depending only on the shape-regularity of the mesh, such that for all neighboring elements K and K .
An interior face of T h is the (nonempty) interior of ∂K + ∩∂K − , where K + and K − are two adjacent elements of T h .Similarly, a boundary face of T h is the (nonempty) interior of ∂K ∩ ∂Ω, which consists of entire faces of ∂K.We denote by F I h the set of all interior faces of T h and by F B h the set of all boundary faces and set Here we generically refer to any element of 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 ∈ F I h be an interior face shared by two neighboring elements K + and K − and let x ∈ F ; we write n ± to denote the unit outward normal vectors on the boundaries ∂K ± .Denoting by v ± the trace of v taken from within K ± , we define the jump and average of v at x ∈ F by respectively.On every boundary face F ∈ F B h , we set [[v]] := vn and { {v} } := v. Here, n is the unit outward normal vector on ∂Ω.
For a piecewise smooth vector-valued function q, we analogously define the average across interior faces by { {q} } := (q + + q − )/2, and on boundary faces we set { {q} } := q.The jump of a vector-valued function will not be used.For a vector-valued function q with continuous normal components across a face f , the trace identity immediately follows from the definitions.

Discretization in space.
For a given partition T h of Ω and an approximation order ≥ 1, we wish to approximate the solution u(t, •) of (2.1)-(2.4) in the finite element space where S (K) is the space P (K) of polynomials of total degree at most on K if K is a triangle or a tetrahedra, or the space Q (K) of polynomials of degree at most 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 Here, Π h denotes the L 2 -projection onto V h , and the discrete bilinear form a h on V h × V h is given by (3.6) The last three terms in (3.6) correspond to jump and flux terms at element boundaries; they vanish when u, v ∈ H 1 0 (Ω) ∩ H 1+σ (Ω) for σ > 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 T h .It is referred to as the interior penalty stabilization function and is defined as follows.We first introduce the function h by For x ∈ F , we further define c by Then, on each F ∈ F h , we set where α 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 α min > 0 which depends only on the shape-regularity of the mesh, the approximation order , the dimension d, and the bounds in (2.5), such that for α ≥ α min where the constant C 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 α 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 α ≥ α min , so that by Lemma 3.1 the semidiscrete problem (3.3)-(3.5)has a unique solution.
We remark that the condition α ≥ α 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].

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.

Main results.
To state our a priori error bounds, we define the space On V (h), we define the DG energy norm Furthermore, for 1 ≤ p ≤ ∞ we will make use of the Bochner space L p (J; V (h)), endowed with the norm Our first main result establishes an optimal error estimate of the energy norm • h of the error.It also gives a bound in the L 2 (Ω)-norm on the error in the first time derivative.Theorem 4.1.Let the analytical solution u of (2.1)-(2.4)satisfy for a regularity exponent σ > 1  2 , and let u h be the semidiscrete discontinuous Galerkin approximation obtained by (3.3)-(3.5),with α ≥ α min .Then, the error e = u − u h satisfies the estimate with a constant C that is independent of T and h.
We remark that the fact that u t ∈ L ∞ (J; H 1+σ (Ω)) implies that u is continuous in time on J with values in H 1+σ (Ω).Similarly, u tt ∈ L 1 (J; H σ (Ω)) implies the continuity of u t on J with values in H σ (Ω).In Theorem 4.1 we thus implicitly assume that the initial conditions satisfy u 0 ∈ H 1+σ (Ω) and v 0 ∈ H σ (Ω).Hence, standard approximation properties imply that see also Lemma 4.6 below.As a consequence, Theorem 4.1 yields optimal convergence in the (DG) energy norm 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 λ ∈ L 2 (Ω) the solution of the problem belongs to H 2 (Ω) and satisfies the stability bound 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 for a regularity exponent σ > 1 2 .Let u h be the semidiscrete DG approximation obtained by (3.3)-(3.5)with α ≥ α min .Then, the error e = u − u h satisfies the estimate 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: 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.

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) × V (h).Indeed the average { {c∇v} } on a face F ∈ F h is not well defined in general for v ∈ H 1 (Ω).To circumvent this difficulty, we shall extend the form a h in a nonstandard and nonconsistent way to the space V (h) × V (h) by using the lifting operators from [2] and the approach in [18].Thus, for v ∈ V (h) we define the lifted function, where c is the material coefficient from (2.1).We shall now show that the lifting operator 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 inv which depends only on the shaperegularity of the mesh, the approximation order , and the dimension d such that Moreover, if the speed of propagation c 1 2 is piecewise constant, with discontinuities aligned with the finite element mesh T h , then Proof.We have Here, we have used the Cauchy-Schwarz inequality, the definition of a in (3.7), and the upper bound for c in (2.5).We recall the inverse inequality with a constant C inv that depends only on the shape-regularity of the mesh, the approximation order , and the dimension d.Using this bound, we obtain ≤ C inv w 0 , which shows the first statement.
With c which completes the proof.Next, we introduce the auxiliary bilinear form The following result establishes that a h is continuous and coercive on the entire space V (h) × V (h); hence it is well defined.Furthermore, since the form a h can be viewed as an extension of the two forms a h and a to the space V (h) × V (h).
Lemma 4.4.Let the interior penalty parameter a be defined as in (3.7), and set for a general piecewise smooth c, and for a piecewise constant c, with discontinuities aligned with the finite element mesh T h .C inv is the constant from Lemma 4.3.
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 For α ≥ α min , the continuity of a h immediately follows.The case of piecewise constant c follows analogously.
To show the coercivity of the form a h , we note that By using the weighted Cauchy-Schwarz inequality, the geometric-arithmetic inequality ab ≤ εa 2 2 + b 2 2ε , valid for any ε > 0, the bounds in (2.5), and the stability bound for the lifting operator in Lemma 4.3, we obtain for general c for a parameter ε > 0 still at our disposal.We conclude that For ε = 1 2 and α ≥ α min , we obtain the desired coercivity bound.For a piecewise constant c we use the bound for c − 1 2 L c (u) 2 0 from Lemma 4.4 and proceed analogously.
Error equation.Because a h coincides with a h on V h × V h , the semidiscrete scheme in (3.3)-(3.5) is equivalent to the following: Find 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 ∈ H 1+σ (Ω) with σ > 1/2 Here Π h denotes the L 2 -projection onto (V h ) d .The assumption u ∈ H 1+σ (Ω) 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 ∈ H 1 0 (Ω).Lemma 4.5.Let the analytical solution u of (2.1)-(2.4)satisfy Let u h be the semidiscrete DG approximation obtained by (4.6).Then, the error e = u − u h satisfies ) almost everywhere in J. Hence, using the discrete formulation in (3.3)-(3.5),we obtain that a.e. in J.
Now, by definition of a h , the fact that L c (u) = 0 and that [[u]] = 0 on all faces, the defining properties of the L 2 -projection Π h , and the definition of the lifted element L c (v), we obtain almost everywhere in J, which implies that c ∇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 From the definition of r h (u, v) in (4.7), we therefore conclude that and obtain where we have used the differential equation (2.1).

Approximation properties.
Let Π h and Π 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 ∈ T h .Then the following hold: 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 , the dimension d, and the regularity exponent t.
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 , the dimension d, and the regularity exponent σ.As a consequence of the approximation properties in Lemma 4.6, we obtain the following results.
Lemma 4.7.Let u ∈ H 1+σ (Ω), σ > 1 2 .Then the following hold: with a constant C A that is independent of the mesh size and depends only on α, the constant κ in (3.1), the bounds in (2.5), and the constants in Lemma 4.6.(ii) For v ∈ V (h), the form r h (u; v) in (4.7) can be bounded by with a constant C R independent of h, which depends only on α, 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 a, and the bounded variation property (3.1).To show the bound in (ii), we apply the Cauchy-Schwarz inequality and obtain .
Applying the approximation properties in Lemma 4.6 completes the proof.

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 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 Lemma 4.7 then implies the two estimates and which concludes the proof of the lemma.
To complete the proof of Theorem 4.1, we now set e = u − u h and recall that Π h is the L 2 -projection onto V h .Because of (2.8), we have Next, we use the symmetry of a h and the error equation in Lemma 4. 5 Integration by parts of the third term on the right-hand side yields From the stability properties of a h in Lemma 4.4 and standard Hölder's inequalities, we conclude that Since this inequality holds for any s ∈ J, it also holds for the maximum over J, that is with Using the geometric-arithmetic mean inequality |ab| ≤ 1 2ε a 2 + ε 2 b 2 , valid for any ε > 0, and the approximation results in Lemma 4.6, we conclude that with a constant C that depends only on the constants in Lemma 4.6.Similarly, where the constant C depends on C coer , C 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 The triangle inequality, the geometric-arithmetic mean, and the approximation properties of Π h in Lemma 4.7 then yield with a constant C that depends only on C coer , C R , and C A .Combining the above estimates for T 1 , T 2 , and T 3 then shows that with a constant that is independent of T and the mesh size.This concludes the proof of Theorem 4.1.

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 ∈ H 1+σ (Ω) with σ > 1 2 , let w h ∈ V h be the solution of Then, we have 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, From the approximation properties of Π h in Lemma 4.7, we immediately infer that

It remains to bound Π h u − w h
h .From the coercivity and continuity of a h in Lemma 4.4, the definition of w h , and the bound in Lemma 4.7, we conclude that Thus, which proves the bound for u − w h h .We shall now prove the L 2 -bound.To do so, let z ∈ H 1 0 (Ω) be the solution of Then, the elliptic regularity assumption in (4.1) and (4.2) implies that Next, we multiply (4.10) by u − w h and integrate the resulting expression by parts.Since c ∇z has continuous normal components across all interior faces, we have with n K denoting the unit outward normal on ∂K.By definition of a h and r h , we immediately find that From the symmetry of a h , the definition of w h , and the fact that [[z]] = 0 on all faces, we conclude that 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 a h , the approximation result in Lemma 4.7 with σ = 1, and the bound in (4.11).Thus, By using Lemma 4.7 and the stability bound in (4.11), we can estimate T 2 by Similarly, The use of these bounds for T 1 , T 2 , and T 3 in (4.12) then leads to which completes the proof of the lemma, since u − w h h ≤ Ch min{σ, } u 1+σ .Now, let u be defined by the exact solution of (2.1)-(2.4).We may define w h (t, •) ∈ V h almost everywhere in J by as well as 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 Moreover, if elliptic regularity as defined in (4.1) and (4.2) holds, we have the L 2bounds The constants C E and C L are as in Lemma 4.9.
To complete the proof of Theorem 4.2, let w h ∈ L ∞ (J; V (h)) be defined by (4.13) and consider 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 ∈ L ∞ (J; V h ) and assume that v t ∈ L ∞ (J; V h ).From the definition of w h in (4.13) and the error equation in Lemma 4.5, we have We rewrite this identity as Let τ ∈ (0, T ] be fixed, and consider the function Since the DG form a h (•, •) is symmetric, we obtain 1 2 Integration over (0, τ ) and using that v(τ, •) = 0 then yield Since u t (0) = v 0 , u h t (0) = Π h v 0 , and v(0) belongs to V h , we conclude that Hence, the first term on the right-hand side of (4.16) vanishes.Moreover, the coercivity of the form a h in Lemma 4.4 ensures that a h ( v(0), v(0)) ≥ 0. This leads to the inequality By using the Cauchy-Schwarz inequality and the geometric-arithmetic mean inequality, we obtain Because this upper bound is independent of τ , it also holds for the supremum over τ ∈ J, which yields the estimate Next, we use this estimate in (4.14) to obtain From the L 2 -approximation properties in Lemma 4.6, Lemma 4.9, and Lemma 4.10, we finally conclude that Here, C is the constant from Lemma 4.6.This completes the proof of Theorem 4.2.

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, 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 α = 20.

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 Here, M denotes the mass matrix and 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 • k.Then the Newmark method consists in finding approximations {u h n } n to u h (t n ) such that 1 See www.dealii.org.
for n = 1, . . ., N − 1.Here, f n := f (t n ), while β ≥ 0 and γ ≥ 1/2 are free parameters that still can be chosen.We recall that for γ = 1/2 the Newmark scheme is secondorder accurate in time, whereas it is only first-order accurate for γ > 1/2.For β = 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 γ = 1/2, the explicit Newmark method corresponds to the standard leap-frog scheme.
For β > 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 γ = 1/2 is unconditionally stable for β ≥ 1/4, whereas for 1/4 > β ≥ 0 the time step k has to be restricted by a CFL condition.In the case β = 0, that condition is k 2 λ max (A) ≤ 4(1 − ε), ε ∈ (0, 1), where λ max (A) is the maximal eigenvalue of the DG stiffness matrix A (which is of the order O(h −2 ) and also depends on α).

Example 1: Smooth solution.
First, we consider the two-dimensional wave equation (2.1)-(2.4) in J × Ω = (0, 1) × (0, 1) 2 , with c ≡ 1 and data f, u 0 , and v 0 chosen such that the analytical solution is given by u(x 1 , x 2 , t) = t 2 sin(πx 1 ) sin(πx 2 ).(5.5)This solution is arbitrarily smooth so that all our theoretical regularity assumptions are satisfied.We discretize this problem using the polynomial spaces Q (K), = 1, 2, 3, on a sequence {T h } i≥1 of square meshes of size h i = 2 −i .With increasing polynomial degree 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 ) for the energy norm and of O(h +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 u(x 1 , x 2 , t) = sin(t 2 ) sin(πx 1 ) sin(πx 2 ).(5.6)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, = 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.everywhere and choose the data f, u 0 , and v 0 such that the analytical solution u is given by u(r, φ, t) = t 2 r 2/3 sin(2/3 φ) (5.7) in polar coordinates (r, φ).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 ∈ C ∞ (J; H 5/3 (Ω)).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 Ω, 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, •) : Here, g is the boundary data and n is the outward unit normal vector on ∂Ω.
We discretize (5.8) by using bilinear polynomials ( = 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 with σ = 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+σ 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 λ ∈ L 2 (Ω) we assume that the solution of the problem belongs to H 1+s (Ω) for a parameter s ∈ (1/2, 1] and satisfies the stability bound z 1+s ≤ C S λ 0 (5.10) 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 ∞ (J;L 2 (Ω)) ≤ Ch min{σ, }+s u 0 1+σ For the L-shaped domain Ω and c ≡ 1, the (weaker) regularity assumption in (5.9)-(5.10)holds with s = 2/3, which underpins the rate σ +s = 4/3 observed in Table 5.1.The wave is locally excited until t = 0.2 by the source term f (x 1 , x 2 , t) = 1, 0.2 < x 1 < 0.4 and t < 0.2, 0, else.
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 ≈ 0.016.The hanging nodes are naturally incorporated in the DG method without any difficulty.Here, the time step k = 0.002,  that is, k ≈ h min /8, proved to be sufficiently small to ensure the stability of the explicit Newmark method (β = 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 ≤ 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.

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: 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.

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