ENERGY CONSERVING LOCAL DISCONTINUOUS GALERKIN METHODS FOR WAVE PROPAGATION PROBLEMS

Wave propagation problems arise in a wide range of applications. The energy conserving property is one of the guiding principles for numerical algorithms, in order to minimize the phase or shape errors after long time integration. In this paper, we develop and analyze a local discontinuous Galerkin (LDG) method for solving the wave equation. We prove optimal error estimates, superconvergence toward a particular projection of the exact solution, and the energy conserving property for the semi-discrete formulation. The analysis is extended to the fully discrete LDG scheme, with the centered second-order time discretization (the leap-frog scheme). Our numerical experiments demonstrate optimal rates of convergence and superconvergence. We also show that the shape of the solution, after long time integration, is well preserved due to the energy conserving property.


1.
Introduction. Wave propagation problems arise in science, engineering and industry, and they are significant to geoscience, petroleum engineering, telecommunication, and the defense industry (see [14,21] and the references therein). Among all the equations describing wave propagation, the constant coefficient wave equation is the simplest one and has been extensively studied. One of the most important properties of the wave equation is the conservation of energy. Experiences reveal that energy conserving numerical methods, which conserve the discrete approximation of energy, are favorable because they are able to maintain the phase and shape of the waves accurately. Numerical methods without this property may result in substantial phase and shape errors after long time integration.
A vast amount of literature can be found on the numerical approximation of the wave equation. All types of numerical methods, including finite difference, finite element, finite volume, spectral methods and integral equation based methods, have their proponents. Among the partial differential equation techniques, the finite difference method provides an efficient solver, but is usually limited by the shape of the domain. The finite element method can handle complex geometry, but often dissipates energy. Here, we will confine our attention in finite element methods, in particular, discontinuous Galerkin (DG) methods. The DG methods belong to a class of finite element methods using discontinuous piecewise polynomial spaces for both the numerical solution and the test functions. They were originally devised to solve hyperbolic conservation laws with only first order spatial derivatives, e.g. [8,9,10,12,13]. They allow arbitrarily unstructured meshes, and have compact stencils. Moreover, they easily accommodate arbitrary h-p adaptivity. The DG methods were later generalized to the local discontinuous Galerkin (LDG) methods by Cockburn and Shu to solve convection-diffusion equations [11], motivated by successful numerical experiments from Bassi and Rebay [3] for the compressible Navier-Stokes equations.
Methods for solving the wave equation can be divided into two categories. The first category is to rewrite the wave equation into a first order hyperbolic system. Therefore, various DG methods designed for first order hyperbolic systems can be applied, for example the standard DG method proposed by Cockburn and Shu [10]. High order nodal DG methods for Maxwell's equations of a first order system are proposed in [20]. Space-time DG methods are studied by Falk and Richter in [15], and later by Monk and Richter in [22]. A suboptimal DG method using central fluxes, which is also energy conserving, has been presented in [16] for Maxwell's equations. More recently, Chung and Engquist [5,6] have proposed an optimal, energy conserving DG method for the wave equation using staggered grids.
The second category of numerical methods is to tackle the second derivatives in the wave equation directly, without introducing more unknowns to reduce the order of the system. Two decades ago, Safjan and Oden [24] introduced a family of unconditionally stable high order Taylor-Galerkin schemes for acoustic and elastic wave propagation. A non-symmetric interior penalty discontinuous Galerkin (IPDG) method has been presented by Rivière and Wheeler in [23]. Later Grote et al. [18,19] proposed a symmetric IPDG method for the second order equation, which conserves a specifically defined energy. Adjerid and Temimi [1] presented a Galerkin method which uses finite element spaces of continuous functions in space, and discontinuous functions in time, and their method achieves the optimal convergence rate and superconvergence property. Local discontinuous Galerkin method has been applied to the wave equation by Baccouch [2] to study its superconvergence property at the Radau points.
In this paper, we primarily consider the one-dimensional linear wave equation subject to the initial conditions We will consider both the homogeneous Dirichlet boundary conditions and the periodic boundary condition u(a, t) = u(b, t).
DG method for this problem with the upwind flux has optimal high order accuracy, but is dissipative. On the other hand, DG method with central flux is energy conserving, but is only suboptimal accurate. Usually it is difficult to obtain DG schemes for wave equations which are non-dissipative (energy conserving for the physical energy) and optimal high order accurate, without going to staggered mesh [5,6] or using penalty with a parameter to adjust [18,19]. Here, we consider the local discontinuous Galerkin method defined on a single mesh for this linear wave problem. The proposed semi-discrete scheme is shown to be energy conserving. It has the optimal convergence rates in both the energy and L 2 norms, and the upper bound of the errors grows in time only in a linear fashion. This method can achieve superconvergence, towards a particular projection of the exact solution. The order of superconvergence is proved to be k+3/2 when piecewise P k polynomials are used. Coupled with the centered second-order time discretization (the leap-frog method) for the temporal derivatives, we present the fully discrete LDG method, which is explicit in time, and we also show the energy conservation property. To demonstrate the importance of energy conservation, we test an example with long time integration, and the numerical results reveal that our method stays very accurate after long time integration, in contrast to numerical methods without this property. We remark here that since our scheme is non-dissipative, it is more oscillatory than the commonly used upwind (energy-dissipative) DG method when applied to problems with discontinuities. The advantage of energy-conserving methods is to solve smooth wave problems, with the attempt to resolve all waves for long time periods. This paper is organized as follows. In Section 2, we present the semi-discrete LDG method, and prove its energy conserving property. The optimal error estimates, both in the energy norm and the L 2 norm, are analyzed in Section 3, and therein, the upper bound of errors is proved to grow linearly in time. In Section 4, we prove the superconvergence of the LDG method toward a particular projection of the exact solution. The fully discrete LDG method, with the leap-frog time discretization, and its energy conserving properties are presented in Section 5. Section 6 contains numerical experiments that demonstrate the optimal convergence rates and energy conservation of the proposed LDG method, and a comparison with the IPDG method proposed in [18]. Finally, we give the concluding remarks in Section 7.

2.1.
Notations. We divide the interval I = [a, b] into N subintervals and denote the cells by , and the mesh size is denoted by h j = x j+ 1 2 − x j− 1 2 , with h = max 1≤j≤N h j being the maximal mesh size. We assume that the mesh is regular, namely, the ratio between the maximal and the minimal mesh sizes stays bounded during mesh refinement. The piecewise polynomial space V k h is defined as the space of polynomials of degree up to k in each cell I j , that is, (5) V k h = {v : v| Ij ∈ P k (I j ), j = 1, 2, · · · , N }. Note that functions in V k h are allowed to have discontinuities across element interfaces.
The solution of the numerical scheme is denoted by u h , which belongs to the finite element space V k h . We denote by (u h ) + to respectively represent the jump and the mean of the function u h at the element interfaces. The L 2 norm over the interval I is denoted by · .
2.2. The LDG method. In this subsection, we define the semi-discrete LDG method for the wave equation (1), by discretizing the space with the LDG method and leaving the time dependence continuous. First, we write the wave equation into a first order system by substituting u x with variable q: The LDG method for (6) is then formulated as follows: for all test functions v, w ∈ V k h . The hatted terms,q h andû h , in (7)-(8) are the cell boundary terms obtained from integration by parts, and they are the so-called numerical fluxes. These numerical fluxes are single-valued functions defined on the cell boundaries and should be designed according to guiding principles for different PDEs to ensure numerical stability. Here we use the simple alternating fluxes: h , in which we have omitted the half-integer indices j + 1 2 , as all quantities in (9) are computed at the same points (i.e., the cell interface). In the case of Dirichlet boundary conditions (3), the numerical fluxes (9) at the two boundaries become = 0, at the two end boundaries. In the proofs of the following sections, the flux terms at the boundary cells (j = 1/2 and N + 1/2) will vanish because u vanishes at the boundary, and therefore, we will retain the flux notation of internal cells for simplicity. We remark that the choice of the fluxes (9) is not unique. We can, for example, alternatively choose the numerical fluxes to be (11)

Energy conservation.
As is well-known, the one-dimensional linear wave equation (1) admits an important conserved quantity -the energy, E = I (u 2 t + u 2 x )dx. Experiences show that schemes conserving the discrete analogs of energy often produce approximations that behave better for long time simulation. In this subsection, we will show that the proposed semi-discrete LDG method conserves energy. Proposition 1. The (continuous) energy is conserved by the the semi-discrete LDG method (7)-(8) for all time.
Proof. We first take the time derivative of Eq. (8) and choose the test function w = q h to obtain In Eq. (7), we choose the test function to be (u h ) t : Adding Eq. (13) to Eq. (14), one obtains Using integration by parts on the term Ij q h (u h ) tx dx, we get If the Dirichlet boundary conditions are used, the flux terms of the boundary cells (j = 1/2 and N + 1/2) in the above equations will vanish, and therefore, we keep the notation of the internal cells for simplicity. By summing up Eq. (15) over all cells and using the periodic or Dirichlet boundary conditions, we have Therefore the quantity E h is invariant in time.
3. Error estimate. In this section, we derive the optimal error estimates for the energy conserving LDG method (7)-(8) proposed in Section 2. The error estimate in the energy norm will be carried out first, and then the analysis will be extended to the L 2 norm. We will also show that these error bounds are both linear in time.
We start by introducing the projections and other notations to be used throughout this paper. The standard L 2 projection of a function ω(x) with k +1 continuous derivatives into space V k h is denoted by P h , i.e., for any v ∈ P k−1 on I j , and Similarly, the projection P + h ω is defined as the projection of ω into V k h , such that for any v ∈ P k−1 on I j , and For these projections, it is easy to show (see [7]): h ω, and Γ h denotes the set of boundary points of all cells. The constant C depends on the function ω, but is independent of the mesh size h.
Let us denote the errors by which, from left to right, respectively represent the errors between the exact solution and the numerical solution, the projection errors, and the errors between the numerical solution and the particular projection of the exact solution. Note that the signs of the projection P ± h of u and q in (17) are consistent with the choice of the numerical fluxes in (9). So if the other set of numerical fluxes (11) is chosen, the signs of P ± h in (17) should be changed accordingly. To obtain the error estimates and the superconvergence property of the proposed LDG method, the projections of the initial conditions for the numerical scheme need to be carefully chosen. Note that we have two initial conditions in (2), one for u and the other for u t . We take the initial condition u h (x, 0) as P + h u(x, 0) = P + h u 0 (x), which is consistent with the choice of the numerical fluxes (9). The other initial condition (u h ) t (x, 0) is given by the standard L 2 projection. Thus, we have the following lemma.
there holds the following error estimate and Proof. Here we only give the proof for the error estimate of ē q (0) . All the other results can be obtained similarly.
Subtracting (8) of the LDG method from the weak formulation satisfied by the exact solution q yields for any test function w ∈ V k h . By the properties of the projection P + h , we can derive Sinceē and therefore, This completes the proof.
Based on initial conditions (18), we have the following error estimate in the energy norm.
Proposition 2. Let u and q be the exact solutions of the wave equation (6), and u h , q h be the numerical solutions of the semi-discrete LDG method (7)-(8), with the numerical fluxes defined in (9) and the initial conditions (18). Consider a regular discretization of I and the piecewise polynomial finite element space V k h , there holds the following error estimates: where the constant C depends on u k+3 and u t k+2 .
Proof. By subtracting the LDG method (7)-(8) from the weak formulation satisfied by the exact solutions u and q, we can derive the error equations for all test functions v, w ∈ V k h . Using the properties of the projections P ± h , the error equations are equivalent to Along the same line in the proof of Proposition 1, we first take the time derivative of Eq. (28), choose the test functions w =ē q , v = (ē u ) t , and sum up the resulting two equations to get If the Dirichlet boundary conditions are used, the flux terms of the boundary cells (j = 1/2 and N + 1/2) in the above equations will vanish, and therefore, we keep the notation of the internal cells for simplicity. By summing up the above equation over all cells and using the periodic or Dirichlet boundary conditions, we have which leads to Combining this inequality with the property of the initial condition (19), we conclude that in which the constant C only depends on u k+3 and u t k+2 . Together with the optimal projection error (16), the error estimate (24) follows.
Next, we consider the error estimate with respect to the L 2 norm.
Proposition 3. Let u and q be the exact solutions of the wave equation (6), and u h , q h be the numerical solutions of the semi-discrete LDG method (7)- (8), with the numerical fluxes defined in (9) and the initial conditions (18). Consider a regular discretization of I and the piecewise polynomial finite element space V k h , there holds the following error estimates: where the constant C only depends on the solution u.
Proof. We split (e u ) tt as the summation of (ē u ) tt and (ε u ) tt in Eq. (27), and use chain rule in time derivative to obtain For any time τ ≤ T , we denote the time integral of the error bȳ Taking the integral of Eq. (28) in time, from t to τ , gives If we choose the test functions to be v =Ē u (t) and w =ē q (t) in (30)-(31), and use the fact that v t = −ē u (t), we have If the Dirichlet boundary conditions are used, the flux terms of the boundary cells (j = 1/2 and N + 1/2) in the above equations will vanish, and therefore, we keep the notation of the internal cells for simplicity. Adding up these two equations and summing over all cells, and using the periodic or Dirichlet boundary conditions, one obtains By integrating the inequality from 0 to τ , we get in which the factĒ q (τ ) =Ē u (τ ) = 0 is used. By the properties of the projection (16), we have (ε u ) t = Ch k+1 . Note that and therefore we can conclude that E ε q = Ch k+1 . Combining with the property of the L 2 projection (20), we have Since this is true for any τ < T , we have Hence, where the constant C only depends on the solution u.

4.
Superconvergence. The superconvergence property of the LDG method is studied in this section. We will prove superconvergence towards a particular projection of the exact solution, with the order k + 3/2, in which an increase of 1/2 over the optimal error estimate of Section 3 can be observed.
To obtain the superconvergence property of the method, two functionals related to the L 2 norm of a function f (x) on I j , as defined in [4], are needed: Properties of these functionals in the following lemma are essential to the proof of the superconvergence. The proof can be found in [4] and is therefore omitted.

Lemma 4.1.
For any function f (x) ∈ C 1 on I j , we have

(See [4] for the detailed proof )
We prove the superconvergence of the numerical solution in the following Proposition.
Proposition 4. Let u and q be the exact solutions of the wave equation (6), and u h , q h be the numerical solutions of the semi-discrete LDG method (7)-(8), with the numerical fluxes defined in (9) and the initial conditions (18). Consider a regular discretization of I and the piecewise polynomial finite element space V k h , there holds the following error estimate: where the constant C only depends on the solution u.
Proof. In the proof of Proposition 3, we have derived an upper bound for the error term ē u (τ ) 2 in Eq. (32). Utilizing the property (20) due to the L 2 projection of the initial condition (u h ) t (0), Eq. (32) becomes for any time τ ≤ T . Firstly, we consider the first term on the right-hand side and tackle the termē u . We rewrite the error equation (28) as by performing integration by parts on the term Ijē u w x dx.
Secondly, we repeat similar analysis to the termĒ q in (36). From Eq. (27), we can derive Integrating the above equation in time from t to τ gives Similar to the previous analysis, we defineĒ q = b j + s j (x)(x − x j )/h j on the cell I j , where b j is a constant and s j (x) ∈ P k−1 . Choose the test function v in (39) to be s j (x)(x − x j+ 1 2 )/h j on I j , and use the definition of B + j (f ) and Lemma 4.1 to obtain Hence, Then we define piecewise polynomials s(x) and φ 2 (x), such that s(x) = s j (x), and φ 2 (x) = x − x j+ 1 2 on I j , and finally get By Proposition 2, we conclude that Lastly, we will use the above results to bound the right-hand side of Eq. (36). By the definition of projections, (ε u ) t and E ε q are both orthogonal to piecewise constant function, hence Define the function φ 3 (x) = (x − x j )/h j on I j , then φ 3 ∞ = 1 2 . Therefore, Combining this inequality with the initial condition (19), and replacing τ by t, we have and in particular ē u (t) ≤ C(t + 1)h k+ 3 2 , where the constant C only depends on the solution u.
Although the proofs for both the error estimate in Section 3 and the superconvergence in Section 4 are shown with the numerical fluxes (9), the same results hold for the other choice of fluxes (11).
5. Time discretization. In order to extend the energy conservation property of the the semi-discrete method to the fully discrete method, it is natural to employ time stepping methods conserving discrete energy. In this paper, we consider the second order leap-frog method for time discretization, which is well-known to be energy conserving.
Let 0 ≤ t 0 < t 1 < · · · < t N = T be a partition of the interval [0, T ] with time step ∆t n = t n+1 − t n . Here uniform time step ∆t is used. The fully discrete approximations u n h to u(·, t n ) are constructed as follows: for n = 1, . . . , N − 1, for all test functions v and w in V k h , and the numerical fluxesû h ,q h are defined as in (9), (10) or (11).
In Proposition 1, we showed the semi-discrete LDG method conserves the continuous energy E h (t). Along the same line of analysis, we can exhibit the following conserved quantity for the fully discrete method. Proof. Suppose the numerical fluxes (9) are used, as the following proof can be carried out similarly for fluxes (11). Choose the test function v = Consider Eq. (42) at time levels t n−1 and t n+1 , and let the test function w be 1 2∆t q n h . Subtracting these two equations yields Adding Eq. (44) to (45) and summing over all cells gives Therefore, by the definition of E n h in (43), we have E n+1 h = E n h for all n, which completes the proof. Remark 1. The discrete energy E n h can be rewritten as: which is a consistent approximation of the continuous energy (12).
6. Numerical experiments. Numerical experiments designed to gauge the performance of our conservative schemes are reported in this section. Attention is given particularly to two issues: 1. Verification of the theoretical results, including a study of the convergence and superconvergence rates. 2. Investigating the long time behavior of our energy conserving scheme, which is also compared with the performance of the energy conserving IPDG method proposed in [19]. This includes a comparison of the errors as a function of time, and the solutions after long time integration.
6.1. Convergence rates. In this subsection, we test the order of convergence and superconvergence of the proposed LDG scheme. The results here present the case of uniform meshes, in which the domain is uniformly divided into N cells. Since the second order leap-frog time discretization is employed and our interest is in the effect of the spatial discretization, we determine the time-step by the relation ∆t = Ch 2 . This relation guarantees that the error will be dominated by the spatial discretization.
We implemented the LDG method with the alternating fluxes (9) and take ∆t = 0.01h 2 . Since leap-frog method requires initial conditions for two time steps, we consider Taylor expansion of u at t = 0: and convert the higher derivatives of t to derivatives of x by repeatedly using the wave equation, while u and u t are given by initial conditions. To obtain the desired order of convergence of u, following the initial conditions (18), we take the projection P + h of u(x, 0), and L 2 projection P h of u t (x, 0); that is, with u(x, 0) denoted by u 0 and u t (x, 0) by v 0 , to use the initial conditions: Tables 1 -3 list the numerical errors and the orders of convergence for P k spaces, k = 1, 2, 3. In each table, the L 2 -norm of the errors e u ,ē u , and e q at final time T = 1 are presented. The (k + 1)th order for e u can clearly be observed. However, the order of convergence ofē u and e q may not be fully revealed from the tables. In order to closely examine the order of convergence, we plotted the numerical errors against mesh sizes in Figs. 1-3. The data sets are fitted by straight line in least-squares sense, and the slopes of fitting lines are calculated, as shown atop the figures. Altogether, the slopes of fitting lines forē u are approximately (k + 2), and the slopes for e q are close to (k + 1). We would like to mention that, although the superconvergence rate is proved to be of order k + 3/2 in Section 4, we can observe the superconvergence of order k + 2 in our numerical experiments. The same phenomenon has been observed in [4].
6.2. Time history of the L 2 error. In this subsection, we investigate the long time evolution of the L 2 error of the LDG method. An energy conserving IPDG method was proposed in [18], which conserves a specifically defined energy. Here we will compare the performance of these two methods. We consider again the wave equation   with a periodic boundary condition u(0, t) = u(2π, t) for all t ≥ 0, and initial conditions u(x, 0) = e sin x , u t (x, 0) = −e sin x cos x. This problem has the exact solution u = e sin (x−t) .
The LDG and IPDG methods are implemented with a uniform mesh with N cells, and the leap-frog time discretization, with ∆t = 0.6h 2 . For the IPDG method, the interior penalty stabilization parameter is taken to be 40000/π. In order to examine the potential shape difference resulted from long time integration, both methods are run until T = 1000, with finite element spaces P 2 and P 3 , and N = 40, 80, respectively.
In Fig. 4, the time evolution of L 2 errors of both schemes are shown, with red color representing LDG method and blue color for IPDG method, and the errors  are recorded every 10,000 iterations. The L 2 errors of both schemes grow in a linear fashion, but the slope for IPDG method is much larger than that for LDG method, which almost stays as constant and is close to zero. In Fig. 4, the errors are plotted in log scale just for easy visualization. From the figure, one can see that for LDG method, the level of the errors are reduced by refining the mesh from N = 40 to N = 80, but the mesh refinement does not substantially reduce the errors of IPDG method due to the rapid growth.
It can be observed from Fig. 4 that, up to T = 1000, the L 2 error of IPDG method is greater than 10 −1 , and this large error can easily be visualized by directly comparing the solutions of both methods. Fig. 5 displays the exact solution (red), the solution of LDG method (green) and the solution of IPDG method (blue) at T = 1000, for spaces P 2 and P 3 with N = 40. It can be seen that solution of LDG method overlaps with the exact solution, while the solution of IPDG method preserves the shape but has a phase shift. A finer mesh with N = 80 for IPDG was also tested, and the results look almost identical to Fig. 5, therefore not displayed here. This suggests that a much finer mesh is required to obtain a more accurate numerical solution. We would like to remark that tuning the interior penalty stabilization parameter in IPDG method does result in solutions with different levels of phase shift. However, evident phase error similar to Fig. 5 can still be observed within a large range of the parameter.    7. Concluding remarks. In this paper, we developed and analyzed the LDG method for solving the wave equation. We have proved the optimal error estimates, superconvergence towards a particular projection of the exact solution, and its energy conserving property for the semi-discrete formulation. The leap-frog time discretization was used to obtain a fully discrete method, and the fully discrete method is shown to conserve discrete energy. Usually it is difficult to obtain DG schemes for wave equations which are non-dissipative and optimal high order accurate, without going to staggered mesh [5,6] or using penalty with a parameter to adjust [18,19]. Our scheme does not have parameters to tune and is defined on a single mesh. Numerical tests have demonstrated the optimal convergence rates and superconvergence which were proved in this paper. Numerical simulations with long time integration, in comparison with results from IPDG methods, confirmed that energy conservation property is important in order to preserve the phase and shape.
Same techniques in the proofs of this paper can be generalized to variable coefficient problem and also to multi-dimensions, which constitute our current work in progress. Future work also includes designing a higher order time stepping method for the proposed LDG scheme, which conserves energy as the leap-frog discretization. Another direction to pursue is to consider the space-time DG method for the wave propagation problem.