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

Authors: Rongpei Zhang; Jia Liu; Shaohua Jiang; Di Wang
Journal: Complexity
 (2020)

---

## 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 inL2-norm of Ohk+1 with polynomial of degree k and grid size h. Numerical experiments are given to verify the theoretical results.


---

## Body

## 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 Ohp+1 convergence rate in L2 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 L2-norm of Ohk+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. 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 theL2 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 L2-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:(1)utt+c2uxx=fx,t,−1≤x≤1,0≤t≤T,with the initial condition(2)ux,0=gx,−1≤x≤1utx,0=hx,and the boundary conditions(3)u−1,t=u1,t,ux−1,t=ux1,t,0≤t≤T,where the constant c is the wave speed. Here, “ux,t” represents the displacement of an elastic bar from its equilibrium position. The functions “g” and “h” in the initial condition are initial displacement and velocity. In our analysis, the functions f∈L20,T;−1,1. The initial condition g∈H01−1,1 and h∈L2−1,1.We first introduce two auxiliary variables:p=cux and w=ut, and then we can rewrite (1) as a first-order system in space:(4)wt=cpx+f,(5)pt=cwx.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 p=cux and w=ut.The computational domainΩ=−1,1 is divided into J cells Ij=xj−1/2,xj+1/2, for 1≤j≤J, in which(6)−1=x1/2<x3/2<⋯<xJ+1/2=1,and we set hj=xj+1/2−xj−1/2 and set hmax=maxjhj and hmin=minjhj as the length of the largest and the smallest subintervals, respectively. The associated finite element space is defined as piecewise polynomial space:(7)Vhk=v:vIj∈PkIj,j=1,…,J,where PkIj denotes the set of all polynomials of degree up to k on Ij. It is worth noting that polynomials in the space Vhk are allowed to have discontinuities across element boundaries.Under the above conditions, the numerical solution in the elementIj can be written as follows:(8)whx,t=∑m=0kwj,mtvjmx=∑m=0kwj,mtPmξ=VTxWjt,where V=P0ξ,P1ξ,…,PkξT denotes the Legendre basis functions and Wjt=wj,0,wj,1,…,wj,kT denote the degrees of freedom (DOF) in cell Ij. The global degrees of freedom in Ω are the set of Wjt as W=W1;W2;…,WJ.wj+1/2+ and wj+1/2− are values of w at xj+1/2 from the right cell Ij+1 and the left cell Ij, respectively:(9)wj+1/2+=limϵ⟶0+wxj+1/2+ϵ,wj+1/2−=limϵ⟶0+wxj+1/2−ϵ,and at each element boundary point xj+1/2, the jump term is denoted by(10)vj+1/2=vj+1/2+−vj+1/2−,and the weighted average is defined by(11)vj+1/2θ=θvj+1/2−+θ˜vj+1/2+,j=1,2,…,J−1,where θ˜=1−θ. In particular, when θ=1, P and Q are Gauss–Radau projection P− and Q−, respectively.The local discontinuous Galerkin method (LDG method) for (4) and (5) is then formulated as follows: find wh,ph∈Vhk so that for all test functions vh,qh∈Vhk and all 1≤j≤J are satisfied:(12)∫Ijwtvdx=−c∫Ijpvxdx+cp^v−j+1/2−cp^v+j−1/2+∫Ijfvdx,(13)∫Ijptqdx=−c∫Ijwqxdx+cw^q−j+1/2−cw^q+j−1/2,where the numerical fluxes w^handp^h are defined as follows:(14)w^hj+1/2=whθ˜,p^hj+1/2=phθ,for j=0,1,2,…,J.


### 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 functionz that is smooth enough, the generalized Gauss–Radau projection of z, Pθz, is defined as the unique element in PkIj satisfying(15)∫IjPθzvdx=∫Ijzvdx,∀v∈Pk−1Ij,j=1,2,…,J,(16)Pθzθxj+1/2=zθxj+1/2,j=1,2,…,J−1.Similarly,Qθ˜z is defined as(17)∫IjQθ˜zqdx=∫Ijzqdx,∀q∈Pk−1Ij,j=1,2,…,J,(18)Qθ˜zθ˜xj+1/2=zθ˜xj+1/2,j=1,2,…,J−1.Under the definition of the above projections, we can obtain the following lemma.Lemma 1.
Let projectionZh be either Pθ or Qθ˜, and assume that the function z is smooth enough. For θ≠1/2, the projection error η=Zhz−z satisfies(19)ηL2Ωh+h1/2ηL2Γh≤Chk+1zW∞k+1Ωh,where C>0 is denoted a generic constant independent of h and z.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 solutionswh,ph obtained by using the LDG method to solve wave equations (12) and (13) satisfy(20)wh2+ph2≤ect∫0tf2dx.Iff=0, wh2+ph2=0 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:(21)∫Ijwhtvhdx=−c∫Ijphvhxdx+cph^vhj−1/2j+1/2+∫Ijfvhdx,∫Ijphtqhdx=−c∫Ijwhqhxdx+cwh^qhj−1/2j+1/2.
Add the above two equations on the whole area, and for the numerical solutions, we have(22)∫Ωphtqhdx+∫Ωwhtvhdx+cBwh,vh,ph,qh=∫Ωfwhdx,where(23)Bwh,vh,ph,qh=∑j=1J∫Ijwhqhxdx+∑j=1J∫Ijphvhxdx+∑j=1J−1wh^qhj+1/2+∑j=1J−1ph^vhj+1/2.
Supposevh=wh, qh=ph; then, taking them into (23), we have(24)Bwh,wh,ph,ph=∑j=1J∫Ijwhphxdx+∑j=1J∫Ijphwhxdx+∑j=1J−1wh^phj+1/2+∑j=1J−1ph^whj+1/2.
For the term∑j=1J∫Ijwhphxdx, by using partial integration, we can obtain(25)∑j=1J∫Ijwhphxdx=−∑j=1J∫Ijphwhxdx−∑j=1Jph+wh+−ph−wh−j+1/2.
For∑j=1Jph+wh+−ph−wh−j+1/2, we have two ways of dealing with it as follows:(26)∑j=1Jph+wh+−ph−wh−j+1/2=∑j=1Jph+wh+−ph+wh−+ph+wh−−ph−wh−j+1/2=∑j=1Jph+wh+wh−phj+1/2,(27)∑j=1Jph+wh+−ph−wh−j+1/2=∑j=1Jph+wh+−wh+ph−+wh+ph−−ph−wh−j+1/2=∑j=1Jwh+ph+ph−whj+1/2.
Substituting (25)–(27) into (24), we can get Bh=0. The following equation can be set up by (22):(28)∫Ωwhtwhdx+∫Ωphtphdx=∫Ωfwhdx.
It is obvious that(29)12ddtwh2+ph2≤fwh.
Integrate both ends of the equation simultaneously:(30)wh2+ph2≤∫0tf2dx+∫0twh2dx.
By the Gronwall inequality, we can obtain(31)wh2+ph2≤ect∫0tf2dx.
It is worth noting thatwh2+ph2=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.Theorem 2.
For the numerical solutionwh of the wave equation obtained by (4) and (5) and the exact solution w, the following error estimates satisfy(32)wx,t−whx,t2≤Ch2k+2,if θ≠1/2.Proof.
First, we will introduce new alternative notations to represent the error terms as follows:(33)ew=w−wh=w−Pθw+Pθw−wh=ηw−εw,ep=p−ph=p−Qθp+Qθp−ph=ηp−εp.
By Lemma1, we can see that(34)ηwL2Ωh+h1/2ηwL2Γh≤Chk+1wk+1,Ωh.
So, in the following, we just need to prove thatεwL2Ωh≤Ch2k+1wk+1,Ωh.
For the exact solutions, it is true that(35)∫Ωptqhdx+∫Ωwtvhdx+cBw,vh,p,qh=∫Ωfvhdx.
Subtract (22) from (35) to get(36)∫Ωηp−εptqhdx+∫Ωηw−εwtvhdx+cBηw−εw,vh,ηp−εp,qh=0.
Takingvh=εw,qh=εp into (36),(37)∫Ωεptεpdx+∫Ωεwtεwdx+cBεw,εw,εp,εp=∫Ωηptεpdx+∫Ωηwtεwdx+cBηw,εw,ηp,εp.
According to the result of the stability analysis, the left hand side of (37) can be rewritten into(38)LHS=12ddtεw2+εp2.
For the termBηw,εw,ηp,εp in the right hand side of (37), it can be derived that(39)Bηw,εw,ηp,εp=∫Ωηpεwxdx+∫Ωηwεpxdx−ηp^εwj−1/2j+1/2−ηw^εpj−1/2j+1/2=∫Ωηpεwxdx+∫Ωηwεpxdx−εp−ηwθ˜j+1/2+εp+ηwθ˜j−1/2−εw−ηpθj+1/2+εw+ηpθj−1/2.
Noticing the definition of Gauss–Radau projections (15)–(18), the term ∫Ωηwεpxdx=0 in (39) since εp∈PkIj and εpx∈Pk−1Ij. Similarly, ∫Ωηpεwxdx=0, εp−ηwθ˜j+1/2=εp+ηwθ˜j−1/2, and εw−ηpθ˜j+1/2=εw+ηpθ˜j−1/2. 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(40)RHS=∫Ωηptεpdx+∫Ωηwtεwdx.
Now equation (37) is changed into(41)12ddtεw2+εp2=∫Ωηptεpdx+∫Ωηwtεwdx.
Considering the basic inequality, we can obtain that(42)∫Ωηptεpdx+∫Ωηwtεwdx≤ηwt2+εw2+ηpt2+εp2.
According to Lemma1, we can get(43)ηwt2≤Ch2k+2,ηpt2≤Ch2k+2.
Combining equation (41), we can get(44)12ddtεw2+εp2≤Ch2k+2+εw2+εp2.
Assume thatϕ=εw2+εp2; when the left and right ends are simultaneously multiplied by the integral factor e−t, (44) can be converted into(45)e−tdϕdt≤e−tCh2k+2+e−tϕ.
For the integration of (45), we can get(46)e−tϕ≤Ch2k+2∫0te−τdτ,where ϕ≤Ch2k+2eT.
Finally, we can obtain that(47)w−wh2≤ηw2+εw2≤Ch2k+2.


## 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 functionz that is smooth enough, the generalized Gauss–Radau projection of z, Pθz, is defined as the unique element in PkIj satisfying(15)∫IjPθzvdx=∫Ijzvdx,∀v∈Pk−1Ij,j=1,2,…,J,(16)Pθzθxj+1/2=zθxj+1/2,j=1,2,…,J−1.Similarly,Qθ˜z is defined as(17)∫IjQθ˜zqdx=∫Ijzqdx,∀q∈Pk−1Ij,j=1,2,…,J,(18)Qθ˜zθ˜xj+1/2=zθ˜xj+1/2,j=1,2,…,J−1.Under the definition of the above projections, we can obtain the following lemma.Lemma 1.
Let projectionZh be either Pθ or Qθ˜, and assume that the function z is smooth enough. For θ≠1/2, the projection error η=Zhz−z satisfies(19)ηL2Ωh+h1/2ηL2Γh≤Chk+1zW∞k+1Ωh,where C>0 is denoted a generic constant independent of h and z.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 solutionswh,ph obtained by using the LDG method to solve wave equations (12) and (13) satisfy(20)wh2+ph2≤ect∫0tf2dx.Iff=0, wh2+ph2=0 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:(21)∫Ijwhtvhdx=−c∫Ijphvhxdx+cph^vhj−1/2j+1/2+∫Ijfvhdx,∫Ijphtqhdx=−c∫Ijwhqhxdx+cwh^qhj−1/2j+1/2.
Add the above two equations on the whole area, and for the numerical solutions, we have(22)∫Ωphtqhdx+∫Ωwhtvhdx+cBwh,vh,ph,qh=∫Ωfwhdx,where(23)Bwh,vh,ph,qh=∑j=1J∫Ijwhqhxdx+∑j=1J∫Ijphvhxdx+∑j=1J−1wh^qhj+1/2+∑j=1J−1ph^vhj+1/2.
Supposevh=wh, qh=ph; then, taking them into (23), we have(24)Bwh,wh,ph,ph=∑j=1J∫Ijwhphxdx+∑j=1J∫Ijphwhxdx+∑j=1J−1wh^phj+1/2+∑j=1J−1ph^whj+1/2.
For the term∑j=1J∫Ijwhphxdx, by using partial integration, we can obtain(25)∑j=1J∫Ijwhphxdx=−∑j=1J∫Ijphwhxdx−∑j=1Jph+wh+−ph−wh−j+1/2.
For∑j=1Jph+wh+−ph−wh−j+1/2, we have two ways of dealing with it as follows:(26)∑j=1Jph+wh+−ph−wh−j+1/2=∑j=1Jph+wh+−ph+wh−+ph+wh−−ph−wh−j+1/2=∑j=1Jph+wh+wh−phj+1/2,(27)∑j=1Jph+wh+−ph−wh−j+1/2=∑j=1Jph+wh+−wh+ph−+wh+ph−−ph−wh−j+1/2=∑j=1Jwh+ph+ph−whj+1/2.
Substituting (25)–(27) into (24), we can get Bh=0. The following equation can be set up by (22):(28)∫Ωwhtwhdx+∫Ωphtphdx=∫Ωfwhdx.
It is obvious that(29)12ddtwh2+ph2≤fwh.
Integrate both ends of the equation simultaneously:(30)wh2+ph2≤∫0tf2dx+∫0twh2dx.
By the Gronwall inequality, we can obtain(31)wh2+ph2≤ect∫0tf2dx.
It is worth noting thatwh2+ph2=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.Theorem 2.
For the numerical solutionwh of the wave equation obtained by (4) and (5) and the exact solution w, the following error estimates satisfy(32)wx,t−whx,t2≤Ch2k+2,if θ≠1/2.Proof.
First, we will introduce new alternative notations to represent the error terms as follows:(33)ew=w−wh=w−Pθw+Pθw−wh=ηw−εw,ep=p−ph=p−Qθp+Qθp−ph=ηp−εp.
By Lemma1, we can see that(34)ηwL2Ωh+h1/2ηwL2Γh≤Chk+1wk+1,Ωh.
So, in the following, we just need to prove thatεwL2Ωh≤Ch2k+1wk+1,Ωh.
For the exact solutions, it is true that(35)∫Ωptqhdx+∫Ωwtvhdx+cBw,vh,p,qh=∫Ωfvhdx.
Subtract (22) from (35) to get(36)∫Ωηp−εptqhdx+∫Ωηw−εwtvhdx+cBηw−εw,vh,ηp−εp,qh=0.
Takingvh=εw,qh=εp into (36),(37)∫Ωεptεpdx+∫Ωεwtεwdx+cBεw,εw,εp,εp=∫Ωηptεpdx+∫Ωηwtεwdx+cBηw,εw,ηp,εp.
According to the result of the stability analysis, the left hand side of (37) can be rewritten into(38)LHS=12ddtεw2+εp2.
For the termBηw,εw,ηp,εp in the right hand side of (37), it can be derived that(39)Bηw,εw,ηp,εp=∫Ωηpεwxdx+∫Ωηwεpxdx−ηp^εwj−1/2j+1/2−ηw^εpj−1/2j+1/2=∫Ωηpεwxdx+∫Ωηwεpxdx−εp−ηwθ˜j+1/2+εp+ηwθ˜j−1/2−εw−ηpθj+1/2+εw+ηpθj−1/2.
Noticing the definition of Gauss–Radau projections (15)–(18), the term ∫Ωηwεpxdx=0 in (39) since εp∈PkIj and εpx∈Pk−1Ij. Similarly, ∫Ωηpεwxdx=0, εp−ηwθ˜j+1/2=εp+ηwθ˜j−1/2, and εw−ηpθ˜j+1/2=εw+ηpθ˜j−1/2. 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(40)RHS=∫Ωηptεpdx+∫Ωηwtεwdx.
Now equation (37) is changed into(41)12ddtεw2+εp2=∫Ωηptεpdx+∫Ωηwtεwdx.
Considering the basic inequality, we can obtain that(42)∫Ωηptεpdx+∫Ωηwtεwdx≤ηwt2+εw2+ηpt2+εp2.
According to Lemma1, we can get(43)ηwt2≤Ch2k+2,ηpt2≤Ch2k+2.
Combining equation (41), we can get(44)12ddtεw2+εp2≤Ch2k+2+εw2+εp2.
Assume thatϕ=εw2+εp2; when the left and right ends are simultaneously multiplied by the integral factor e−t, (44) can be converted into(45)e−tdϕdt≤e−tCh2k+2+e−tϕ.
For the integration of (45), we can get(46)e−tϕ≤Ch2k+2∫0te−τdτ,where ϕ≤Ch2k+2eT.
Finally, we can obtain that(47)w−wh2≤ηw2+εw2≤Ch2k+2.


## 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:(48)M1dwdt=A1P+F,M2dPdt=A2W,where U and P are solution vectors containing the degrees of freedom of u and p. Because the mass matrices M1 and M2 are block diagonal, the mass matrices can be inverted easily to get(49)ddtWP=0M1−1A1M2−1A20WP+M−1F10.SetΦ=WP and F=M−1F10, A=0M1−1A1M2−1A20. Then, (49) can be rewritten into(50)ddtΦ=AΦ+F.Assume the final time ist=T and let time step Δt=T/N, tn=nΔt,0≤n≤N. We multiply (50) by the integration factor e−At and integrate over one time step from tn to tn+1 to obtain(51)Φn+1=eAΔtΦn+eAΔt∫0Δte−AτFτdτ.Using the trapezoidal integration to calculate the integral in (52), the second-order exponential integration factor formula is obtained:(52)Φn+1=eAΔtΦn+Δt2Fn+Δt2Fn+1.It is noticed that iff=0, the exponential integration factor is exact. To construct such a scheme, we need to formulate an efficient algorithm for computing the products of matrices eAΔt and vector Φn. In this paper, we accomplish the product by using Krylov subspace projections. The matrix A and vector Φn is projected on a Krylov subspace:(53)KmA,Φn=spanΦn,AΦn,…,Am−1Φn.The orthonormal basisv1,v2,…,vm in the subspace KmA,Φ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 Hm of the matrix A is constructed in the process of orthogonalization. Then, the obtained results are used to evaluate the approximation to eAΔt in KmA,Φn. In this section, we take m=25.


## 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 spacesPk 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:(54)utt−uxx=0,x∈−1,1,t∈0,T,ux,0=sin10πx,utx,0=10πcos10πx,x∈−1,1,u−1,t=u1,t,ux−1,t=ux1,t.
The exact solution of this problem can be given byux,t=sin10πx+t.
We solve this problem using the LDG method on uniform meshes havingJ=32,64,128, and 256. The time step is chosen as Δt=10−3. Tables 1 and 2 list the L2 errors and corresponding convergence orders for wh and ph. In each table, the parameter θ is taken as θ=0.2,0.5, and 1.0. The 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. This coincides with Theorem 2 and shows the sharpness of theoretical results.
We present the true and LDG solution at timeT=100 in Figure 1 with finite element spaces P2 and P3 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 wh at time T=1 and T=100. These computations are also performed with finite element spaces P2 and P3 and J=100. The obtained results indicate that the L2 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.Table 1
TheL2 spatial errors and order of convergence of wh for Example 1.
mθ=0.2θ=0.5θ=1.0ErrorOrderErrorOrderErrorOrderP262.15E − 1—2.08E − 1—3.08E − 1—71.80E − 23.581.72E − 23.591.88E − 24.0383.99E − 32.172.07E − 33.052.22E − 33.0593.74E − 43.412.59E − 42.992.66E − 43.09P362.23E − 2—2.44E − 2—1.45E − 2—71.00E − 34.473.87E − 32.968.94E − 44.0288.70E − 53.524.96E − 42.968.45E − 53.4092.92E − 64.895.78E − 53.104.84E − 64.13Table 2
TheL2 spatial errors and order of convergence of ph for Example 1.
mθ=0.2θ=0.5θ=1.0ErrorOrderErrorOrderErrorOrderP261.74E − 1—2.08E − 1—3.49E − 1—71.98E − 23.131.72E − 23.594.06E − 23.1083.01E − 32.722.07E − 33.054.00E − 33.3493.73E − 43.012.59E − 42.993.59E − 43.48P362.25E − 2—2.44E − 2—1.25E − 2—71.87E − 33.583.87E − 32.967.27E − 44.1181.33E − 43.814.96E − 42.966.10E − 53.5895.06E − 64.725.78E − 53.104.36E − 64.16Figure 1
The exact solutionw (solid line) and LDG solution wh (dashed red dots) at T=100 by P2 (a) and P3 (b) for Example 1.
(a)(b)Figure 2
The errorwx,t−whx,t versus x at T=5 by P2 (a) and P3 (b) for Example 1.
(a)(b)Figure 3
The errorwx,t−whx,t versus x at T=100 by P2 (a) and P3 (b) for Example 1.
(a)(b)Example 2.
We then consider the nonhomogeneous wave problem:(55)utt−uxx=fx,t,x∈−1,1,t∈0,T.
The exact solution of this problem is given byux,t=sinπxe−0.1t. The initial condition and fx,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. The L2 norm errors and corresponding convergence orders at time T=1 are shown in Table 3 for ph. The error for wh 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 T=10 with finite element spaces P2 and P3 is shown in Figure 4. The evolution of error curve is presented in Figure 5 with P2 and P3 elements.Table 3
TheL2 spatial errors and order of convergence of ph for Example 2.
mθ=0.2θ=0.5θ=1.0ErrorOrderErrorOrderErrorOrderP261.19E − 3—7.80E − 4—1.67E − 3—71.41E − 43.079.59E − 53.022.20E − 42.9281.71E − 53.051.20E − 53.001.95E − 53.4992.30E − 62.891.52E − 62.971.59E − 63.60P344.54E − 5—1.78E − 4—1.99E − 5—53.37E − 63.752.82E − 52.661.17E − 64.0963.02E − 73.483.60E − 62.978.13E − 83.8572.79E − 83.445.29E − 72.774.69E − 94.11Figure 4
The exact solutionw (solid line) and LDG solution wh (dashed red dots) at T=10 by P2 (a) and P3 (b) for Example 2.
(a)(b)Figure 5
The errorwx,t−whx,t versus x at T=10 by P2 (a) and P3 (b) for Example 2.
(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.



---

*Source: 8464153-2020-08-19.xml*