The Local Discontinuous Galerkin Method with Generalized Alternating Flux Applied to the Second-Order Wave Equations

In this paper, we propose the local discontinuous Galerkin method based on the generalized alternating numerical flux for solving the one-dimensional second-order wave equation with the periodic boundary conditions. Introducing two auxiliary variables, the second-order equation is rewritten into the first-order equation systems. We prove the stability and energy conservation of this method. By virtue of the generalized Gauss–Radau projection, we can obtain the optimal convergence order in L2-norm of O(hk+1) with polynomial of degree k and grid size h. Numerical experiments are given to verify the theoretical results.


Introduction
e wave equation is an important second-order linear partial differential equation for the description of many traveling wave phenomena, such as the sound waves and seismic waves. Accurate and efficient numerical methods to solve the wave equation are of fundamental importance to the simulation of wave propagation. ere are many numerical methods to solve the wave equation, such as the finite element methods (FEMs) [1,2], spectral methods [3,4], and finite volume method [5].
In this paper, we consider the local discontinuous Galerkin (LDG) methods for numerical solution of onedimensional second-order wave equation. e LDG method was first introduced by Cockburn and Shu [6] which was motivated by Bassi and Rebay's work on Navier-Stokes problems [7]. With discontinuous finite element basis function, the LDG method possesses certain flexibility and advantage, such as local (elementwise) conservation, capability of capturing discontinuity, and high accuracy. It can also be easily designed for any order of accuracy since the order of accuracy can be locally determined in each cell.
Many DG methods have been developed for the wave equation in both first-order and second-order forms [2,8,9], and some of them are also energy conserving. For instance, in [10], the authors solved the time-dependent acoustic and elastic wave equations by using the discontinuous Galerkin method (DG) for spatial discretization and Crank-Nicolson methods for time integration. And a new class of discontinuous Galerkin (DG) methods for the acoustic wave equation in mixed form is developed and analyzed in [11]. In [9], the new superconvergence results for the local discontinuous Galerkin (LDG) method applied to the secondorder scalar wave equation in one space dimension are presented, and the results show O(h p + 1) convergence rate in L 2 norm. Chou et al. solved the linear second-order wave equation in heterogeneous media by energy conserving the LDG method [12].
In this paper, we study a new local discontinuous Galerkin (LDG) method for solving the second-order wave equation. We first introduce two auxiliary variables to rewrite the second-order equation as the velocity-stress form of the wave equation. is system is then discretized by the discontinuous Galerkin (DG) method to get the solution of the introduced variable. In the design of the LDG method, the numerical flux should be chosen to ensure the stability and high order accuracy. e LDG methods in the previous papers [9,12] applied the purely alternating numerical flux. In this paper, we consider the generalized alternating flux which is motivated by work of Meng et al. [13]. Recently, Cheng et al. [14,15] developed the LDG method with generalized numerical flux for convection-diffusion equations. ey obtained the optimal error estimate based on the construction and analysis of the newly designed generalized Gauss-Radau projections. e local discontinuous Galerkin method based on the generalized alternating numerical flux is also used for solving the one-dimensional nonlinear Burger's equation [16]. By using the generalized Gauss-Radau projection, we can prove that the proposed method can get optimal rate of convergence in L 2 -norm of O(h k+1 ) with polynomial of degree k for wave equation.
After the LDG discretization of the wave equation, we get a system of ordinary differential equations (ODEs) which will be solved by the exponential integration factor method. e DG discretization leads to a relatively large number of degrees of freedom in comparison to other discretization methods. To efficiently evaluate the product of the matrix exponential and a vector, we perform Krylov subspace approximations.
Compared with the existing numerical method for wave equation, we believe our method has a number of other attractive properties. First, the generalized flux in our paper admits a wide variety of numerical flux choices which can ensure the energy conservation, including central flux and alternating flux. Second, we prove the optimal convergence in the L 2 norm in one space dimension. e optimal convergence for a variety of flux choices is observed experimentally. Finally, We use the Krylov exponential integration method which has a much better performance in computer time. e outline of the paper is as follows. In Section 2, we will present the LDG method with generalized alternating numerical flux for the first-order system transformed from the second-order wave equation. en, the construction and analysis of the generalized Gauss-Radau projections are stated in Section 2.1. After that, we will give the stability and a priori error estimates in the L 2 -norm in Sections 2.2 and 2.3. In Section 3, we then move to the time discretization of semidiscrete form by the matrix exponential method. In Section 4, we present some numerical experiments that confirm our theoretical analysis.

The LDG Method
In this paper, we consider the following one-dimensional second-order wave equations: with the initial condition and the boundary conditions where the constant c is the wave speed. Here, "u(x, t)" represents the displacement of an elastic bar from its equilibrium position. e functions "g" and "h" in the initial condition are initial displacement and velocity. In our analysis, the functions f ∈ L 2 (0, T; [− 1, 1]). e initial condition g ∈ H 1 0 ([− 1, 1]) and h ∈ L 2 ([− 1, 1]). We first introduce two auxiliary variables: p � cu x and w � u t , and then we can rewrite (1) as a first-order system in space: Acoustic wave equations (4) and (5) are known as velocity-stress form, while equation (1) is known as displacement form.
e two forms are equivalent through and we set h j � x j+(1/2) − x j− (1/2) and set h max � max j h j and h min � min j h j as the length of the largest and the smallest subintervals, respectively. e associated finite element space is defined as piecewise polynomial space: where P k (I j ) denotes the set of all polynomials of degree up to k on I j . It is worth noting that polynomials in the space V k h are allowed to have discontinuities across element boundaries. Under the above conditions, the numerical solution in the element I j can be written as follows: where V � (P 0 (ξ), P 1 (ξ), . . . , P k (ξ)) T denotes the Legendre basis functions and W j (t) � (w j,0 , w j,1 , . . . , w j,k ) T denote the degrees of freedom (DOF) in cell I j . e global degrees of freedom in Ω are the set of W j (t) as W � [W 1 ; W 2 ; . . . , W J ]. w + j+(1/2) and w − j+(1/2) are values of w at x j+(1/2) from the right cell I j+1 and the left cell I j , respectively: and at each element boundary point x j+ (1/2) , the jump term is denoted by 2 Complexity and the weighted average is defined by where θ � 1 − θ. In particular, when θ � 1, P and Q are Gauss-Radau projection P − and Q − , respectively. e local discontinuous Galerkin method (LDG method) for (4) and (5) is then formulated as follows: find where the numerical fluxes w h and p h are defined as follows: for j � 0, 1, 2, . . . , J.

Generalized Gauss-Radau Projection.
In order to obtain the optimal error estimates, we first introduce two forms of generalized Gauss-Radau projections: for any function z that is smooth enough, the generalized Gauss-Radau projection of z, P θ z, is defined as the unique element in P k (I j ) satisfying Similarly, Q θ z is defined as Under the definition of the above projections, we can obtain the following lemma. Lemma 1. Let projection Z h be either P θ or Q θ , and assume that the function z is smooth enough. For θ ≠ (1/2), the projection error η � Z h z − z satisfies where C > 0 is denoted a generic constant independent of h and z.
e proof of this lemma is not repeated in this paper. For the concrete proof process, the reader is referred to the analysis of Lemma 3.2 in [14].

Stability
Analysis. Now, we state the stability analysis of the above LDG method based on the generalized alternating numerical fluxes.

Theorem 1.
e numerical solutions w h , p h obtained by using the LDG method to solve wave equations (12) and (13) satisfy Proof. First, considering the definition of Gauss-Radau projections, (12) and (13) can be transformed into the following form: Add the above two equations on the whole area, and for the numerical solutions, we have where we have two ways of dealing with it as follows: Substituting (25)-(27) into (24), we can get B h � 0. e following equation can be set up by (22): It is obvious that Integrate both ends of the equation simultaneously: By the Gronwall inequality, we can obtain It is worth noting that ‖w h ‖ 2 + ‖p h ‖ 2 � 0 can be established when f � 0. □ 2.3. Error Estimate. In this section, we will discuss the error of the LDG method ((4) and (5)) with the periodic boundary. (4) and (5) and the exact solution w, the following error estimates satisfy

Theorem 2. For the numerical solution w h of the wave equation obtained by
Proof. First, we will introduce new alternative notations to represent the error terms as follows: By Lemma 1, we can see that So, in the following, we just need to prove that For the exact solutions, it is true that

Subtract (22) from (35) to get
Taking According to the result of the stability analysis, the left hand side of (37) can be rewritten into For the term B(η w , ε w , η p , ε p ) in the right hand side of (37), it can be derived that Noticing the definition of Gauss-Radau projections In this case, we can get the conclusion that B(η w , ε w , η p , ε p ) � 0.
So the right hand side of (37) can be transformed into Now equation (37) is changed into Considering the basic inequality, we can obtain that (42) According to Lemma 1, we can get Combining equation (41), we can get Assume that ϕ � ‖ε w ‖ 2 + ‖ε p ‖ 2 ; when the left and right ends are simultaneously multiplied by the integral factor e − t , (44) can be converted into For the integration of (45), we can get where ϕ ≤ Ch 2k+2 e T . Finally, we can obtain that

Time Discretization: Krylov Implicit Integration Factor Method
After the LDG discretization of (4) and (5) in space, a set of ordinary differential equations is obtained: where U and P are solution vectors containing the degrees of freedom of u and p. Because the mass matrices M 1 and M 2 are block diagonal, the mass matrices can be inverted easily to get d dt en, (49) can be rewritten into Assume the final time is t � T and let time step Δt � T/N, t n � nΔt, 0 ≤ n ≤ N. We multiply (50) by the integration factor e − At and integrate over one time step from t n to t n+1 to obtain Δt 0 e − Aτ F(τ)dτ. (51) Using the trapezoidal integration to calculate the integral in (52), the second-order exponential integration factor formula is obtained: It is noticed that if f � 0, the exponential integration factor is exact. To construct such a scheme, we need to formulate an efficient algorithm for computing the products     6 Complexity of matrices e AΔt and vector Φ n . In this paper, we accomplish the product by using Krylov subspace projections. e matrix A and vector Φ n is projected on a Krylov subspace: e orthonormal basis v 1 , v 2 , . . . , v m in the subspace K m (A, Φ n ) is constructed using an algorithm based on the Gram-Schmidt orthogonalization. A numerically preferable version of this algorithm is the Arnoldi modified Gram-Schmidt procedure. We refer the reader to the paper [17] for the detailed description of the Arnoldi algorithm. A Krylov subspace representation H m of the matrix A is constructed in the process of orthogonalization. en, the obtained results are used to evaluate the approximation to e AΔt in K m (A, Φ n ). In this section, we take m � 25.

Numerical Experiments
In this section, the solutions of the wave equation will be investigated by using the proposed LDG method. e LDG method is carried out using the spaces P k with k � 2, 3 in (7). All experiments were carried out via Matlab on a PC with Intel Core i7 1.8 GHz processor and 16 GB RAM. Example 1. We first consider the following homogeneous wave problem with the periodic boundary conditions: (1, t), e exact solution of this problem can be given by u(x, t) � sin(10π(x + t)).
We solve this problem using the LDG method on uniform meshes having J � 32, 64, 128, and 256. e time step is chosen as Δt � 10 − 3 . Tables 1 and 2 list the L 2 errors and corresponding convergence orders for w h and p h . In each table, the parameter θ is taken as θ � 0.2, 0.5, and 1.0. e tables show different convergence orders with different θ. For θ ≠ (1/2), the optimal accuracy of (k + 1)-th order for k � 2, 3 is obtained. For θ � 1/2, the suboptimal accuracy of k-th order for k � 3 and optimal (k + 1)-th order with for k � 2 can be observed. is coincides with eorem 2 and shows the sharpness of theoretical results.
We present the true and LDG solution at time T � 100 in Figure 1 with finite element spaces P 2 and P 3 and J � 100. Figure 1 shows that the shape of the solution after long time integration is well preserved by using the LDG method. In Figures 2 and 3, we show the time evolution of the error curves of w h at time T � 1 and T � 100. ese computations are also performed with finite element spaces P 2 and P 3 and J � 100. e obtained results indicate that the L 2 errors do not grow in time and match well with the numerical results reported in [18]. e example shows that the proposed LDG scheme has very good dissipation and dispersion properties.
Example 2. We then consider the nonhomogeneous wave problem: e exact solution of this problem is given by u(x, t) � sin(πx)e − 0.1t . e initial condition and f(x, t) are determined by the exact solution.
We start by solving (55) on uniform meshes having J � 16, 32, 64, and 128. We choose the same parameters of θ as Example 1.
e L 2 norm errors and corresponding      Table 3 for p h . e error for w h is similar and we omit it. We can get the similar convergence order as in Example 1, which verifies the theoretical results in Section 2. en, we show that the shape of the solution by LDG method will preserve after long time integration is well preserved. e shape of LDG solutions at T � 10 with finite element spaces P 2 and P 3 is shown in Figure 4. e evolution of error curve is presented in Figure 5 with P 2 and P 3 elements.

Summary
In this paper, we have applied the LDG method with generalized alternating flux for solving the one-dimensional second-order wave equation with the periodic boundary conditions. e second-order wave equation is transformed into the first-order wave equation by introducing a variable. e first-order wave equation is then discretized by the LDG method which results in a stiff system of first-order ordinary differential equations. is system of ordinary differential equations is solved by the exponential matrix method analytically. Numerical results are compared with exact solutions at different times. e obtained results confirm that our LDG method is a powerful and reliable method for capturing the shock phenomenon. e methods extended to higher dimensional case will be our future work.

Data Availability
e numerical data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.