Abstract
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 -norm of with polynomial of degree and grid size . Numerical experiments are given to verify the theoretical results.
1. Introduction
The 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. There 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 one-dimensional second-order wave equation. The 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 second-order scalar wave equation in one space dimension are presented, and the results show convergence rate in 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. This 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. The 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. They obtained the optimal error estimate based on the construction and analysis of the newly designed generalized Gauss–Radau projections. The 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 -norm of with polynomial of degree 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. The 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 norm in one space dimension. The 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. The 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. Then, 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 -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.
2. The LDG Method
In this paper, we consider the following one-dimensional second-order wave equations:with the initial conditionand the boundary conditionswhere the constant is the wave speed. Here, “” represents the displacement of an elastic bar from its equilibrium position. The functions “” and “h” in the initial condition are initial displacement and velocity. In our analysis, the functions . The initial condition and .
We first introduce two auxiliary variables: and , 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. The two forms are equivalent through and .
The computational domain is divided into cells , for , in whichand we set and set and as the length of the largest and the smallest subintervals, respectively. The associated finite element space is defined as piecewise polynomial space:where denotes the set of all polynomials of degree up to on . It is worth noting that polynomials in the space are allowed to have discontinuities across element boundaries.
Under the above conditions, the numerical solution in the element can be written as follows:where denotes the Legendre basis functions and denote the degrees of freedom (DOF) in cell . The global degrees of freedom in are the set of as .
and are values of at from the right cell and the left cell , respectively:and at each element boundary point , the jump term is denoted byand the weighted average is defined bywhere . In particular, when , and are Gauss–Radau projection and , respectively.
The local discontinuous Galerkin method (LDG method) for (4) and (5) is then formulated as follows: find so that for all test functions and all are satisfied:where the numerical fluxes are defined as follows:for .
2.1. 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 that is smooth enough, the generalized Gauss–Radau projection of , , is defined as the unique element in satisfying
Similarly, is defined as
Under the definition of the above projections, we can obtain the following lemma.
Lemma 1. Let projection be either or , and assume that the function is smooth enough. For , the projection error satisfieswhere is denoted a generic constant independent of and .
The 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].
2.2. Stability Analysis
Now, we state the stability analysis of the above LDG method based on the generalized alternating numerical fluxes.
Theorem 1. The numerical solutions obtained by using the LDG method to solve wave equations (12) and (13) satisfy
If , is a constant for all time.
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 havewhereSuppose , ; then, taking them into (23), we haveFor the term , by using partial integration, we can obtainFor , we have two ways of dealing with it as follows:Substituting (25)–(27) into (24), we can get . The following equation can be set up by (22):It is obvious thatIntegrate both ends of the equation simultaneously:By the Gronwall inequality, we can obtainIt is worth noting that can be established when .
2.3. Error Estimate
In this section, we will discuss the error of the LDG method ((4) and (5)) with the periodic boundary.
Theorem 2. For the numerical solution of the wave equation obtained by (4) and (5) and the exact solution , the following error estimates satisfyif .
Proof. First, we will introduce new alternative notations to represent the error terms as follows:By Lemma 1, we can see thatSo, in the following, we just need to prove that .
For the exact solutions, it is true thatSubtract (22) from (35) to getTaking into (36),According to the result of the stability analysis, the left hand side of (37) can be rewritten intoFor the term in the right hand side of (37), it can be derived thatNoticing the definition of Gauss–Radau projections (15)–(18), the term in (39) since and . Similarly, , , and . In this case, we can get the conclusion that .
So the right hand side of (37) can be transformed intoNow equation (37) is changed intoConsidering the basic inequality, we can obtain thatAccording to Lemma 1, we can getCombining equation (41), we can getAssume that ; when the left and right ends are simultaneously multiplied by the integral factor , (44) can be converted intoFor the integration of (45), we can getwhere .
Finally, we can obtain that
3. 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 and are solution vectors containing the degrees of freedom of and . Because the mass matrices and are block diagonal, the mass matrices can be inverted easily to get
Set and , . Then, (49) can be rewritten into
Assume the final time is and let time step , . We multiply (50) by the integration factor and integrate over one time step from to to obtain
Using the trapezoidal integration to calculate the integral in (52), the second-order exponential integration factor formula is obtained:
It is noticed that if , the exponential integration factor is exact. To construct such a scheme, we need to formulate an efficient algorithm for computing the products of matrices and vector . In this paper, we accomplish the product by using Krylov subspace projections. The matrix and vector is projected on a Krylov subspace:
The orthonormal basis in the subspace 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 of the matrix is constructed in the process of orthogonalization. Then, the obtained results are used to evaluate the approximation to in . In this section, we take .
4. Numerical Experiments
In this section, the solutions of the wave equation will be investigated by using the proposed LDG method. The LDG method is carried out using the spaces with 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:The exact solution of this problem can be given by .
We solve this problem using the LDG method on uniform meshes having , and 256. The time step is chosen as . Tables 1 and 2 list the errors and corresponding convergence orders for and . In each table, the parameter is taken as , and 1.0. The tables show different convergence orders with different . For , the optimal accuracy of -th order for is obtained. For , the suboptimal accuracy of -th order for and optimal -th order with for can be observed. This coincides with Theorem 2 and shows the sharpness of theoretical results.
We present the true and LDG solution at time in Figure 1 with finite element spaces and and . 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 at time and . These computations are also performed with finite element spaces and and . The obtained results indicate that the errors do not grow in time and match well with the numerical results reported in [18]. The example shows that the proposed LDG scheme has very good dissipation and dispersion properties.
(a)
(b)
(a)
(b)
(a)
(b)
Example 2. We then consider the nonhomogeneous wave problem:The exact solution of this problem is given by . The initial condition and are determined by the exact solution.
We start by solving (55) on uniform meshes having , and 128. We choose the same parameters of as Example 1. The norm errors and corresponding convergence orders at time are shown in Table 3 for . The error for 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. Then, we show that the shape of the solution by LDG method will preserve after long time integration is well preserved. The shape of LDG solutions at with finite element spaces and is shown in Figure 4. The evolution of error curve is presented in Figure 5 with and elements.
(a)
(b)
(a)
(b)
5. 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. The second-order wave equation is transformed into the first-order wave equation by introducing a variable. The first-order wave equation is then discretized by the LDG method which results in a stiff system of first-order ordinary differential equations. This system of ordinary differential equations is solved by the exponential matrix method analytically. Numerical results are compared with exact solutions at different times. The obtained results confirm that our LDG method is a powerful and reliable method for capturing the shock phenomenon. The methods extended to higher dimensional case will be our future work.
Data Availability
The numerical data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was partially supported by the Foundation of LCP (no. 6142A0502020717) and Natural Science Foundation of Liaoning Province (grant no. 20180550996).