 # Discontinuous Galerkin time stepping methods for second order hyperbolic problems

Discontinuous Galerkin methods, based on piecewise polynomials of degree q≥ 0, are investigated for temporal semi-discretization for second order hyperbolic equations. Energy identities and stability estimates of the discrete problem are proved for a slightly more general problem, that are used to prove optimal order a priori error estimates with minimal regularity requirement. Uniform norm in time error estimates are proved for the constant and linear cases. Numerical experiments are performed to verify the theoretical rate of convergence.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

We study discontinuous Galerkin methods of order , dG(q), for temporal semi-discretization of the second order hyperbolic problems

 (1.1) ¨u+Au=f,t∈(0,T), with u(0)=u0, ˙u(0)=v0,

where is a self-adjoint, positive definite, uniformly elliptic second-order operator on a Hilbert space .

We may consider, as a prototype equation for such second order hyperbolic equations, with homogeneous Dirichlet boundary conditions. That is, the classical wave equation,

 (1.2) ¨u(x,t)−Δu(x,t)=f(x,t)inΩ×(0,T),u(x,t)=0onΓ×(0,T),u(x,0)=u0(x),˙u(x,0)=v0(x)inΩ,

where is a bounded and convex polygonal domain in , , with boundary . We denote and

. The present work applies also to wave phenomena with vector valued solution

, such as wave elasticity.

We recall that for the homogeneous parabolic problem

 ˙u+Au=0,t∈(0,T), with u(0)=u0,

the solution is given by , where the analytic semigroup , , is uniformly bounded. Furthermore, we have the smoothing properties,

 (1.3) ∥AγE(t)v∥≤Ct−γ∥v∥,∀v∈H, t>0, 0≤γ≤1,

that has a crucial impact on analysis of parabolic type problems, in particular, for non-smooth data analysis, . However, such smoothing properties does not hold for the wave equation.

The dG methods can be traced back to the ’s where they were proposed as variational methods for numerically solving initial-value problem and transport problems, see, e.g., [6, 12, 17, 21]

. The discontinuous Galerkin type methods for time and/or space discretization have been studied extensively in the literature for ordinary differential equations and parabolic/hyperbolic partial differential equations; see, for example,

[1, 2, 3, 4, 5, 6, 7, 9, 15, 16, 18, 24, 25] and the references therein. In particular, several discontinuous and continuous Galerkin finite element methods, both in time and space variables, for solving second order hyperbolic equations have appeared in the literature, see e.g., [2, 10, 11, 13, 22] and the references therein.

Uniform in time stability analysis, also so-called strong stability or -stability, has been studied for parabolic problems, see, e.g., [8, 16, 24]. An important tool for such analysis is based on the smoothing property (1.3) of the solution operator for parabolic problems, that is due to analytic semigroups.

In , uniform in time stability and error estimates for dG(q), , have been proved using Dunford-Taylor formula based on smoothing properties of the analytic semigroups. For parabolic problems which is perturbed by a memory term, such analysis has been done for dG(0) and dG(1), using the linearity of the basis functions, see . Another way to analyse uniform in time stability is using a lifting operator technique to write the dG(q) formulation in a strong (pointwise) form, see  and the references therein. For the second order hyperbolic problems, say the wave equation, there is no smoothing property for the solution operator, lacking an analytic semigroup.

In the present work, we formulate the discontinuous Galerkin method of arbitrary integer order for (1.1), in particular, for (1.2). We prove energy identity and stability estimates for the discrete problem of a more general form, by considering an extra (artificial) load term in the so called displacement-velocity formulation. This gives the flexibility to obtain optimal order a priori error estimates with minimal regularity requirement on the solution. We prove optimal order a priori error estimates in and norms for the displacement and -norm of the velocity . The stability constants are independent of the length of the time interval, , that means long-time integration without error accumulation is possible. Similar idea has been used for error analysis of the finite element approximation of the second order hyperbolic equations, see, e.g., [14, 23]. Here, we show that how it can be applied for dG(q) time stepping methods, using energy argument, for stability and error analysis. The other contribution of this work is to prove uniform in time optimal order a priori error estimates for dG(0) and dG(1) based on linearity of the approximate solutions, lacking smoothing property for the solution operator. We also think that our techniques for stability and error analysis is straight forward.

The outline of this paper is as follows. We provide some preliminaries and the weak formulation of the model problem, in . In section 3, we formulate the dG(q) method, and we obtain energy identity and stability estimates for the discrete problem of a slightly more general form. Then, in , we prove optimal order a priori error estimates in and norms for the displacement and -norm of the velocity, with minimal regularity requirement on the solution. We also prove uniform in time a priori error estimates for dG(0) and dG(1), based on linearity of the approximate solutions. Finally, numerical experiments are presented in section 5 in order to illustrate the theory.

## 2. Preliminaries

We let with the inner product and the induced norm . Denote with the energy inner product and the induced norm . Let be defined with homogeneous Dirichlet boundary conditions on , and be the eigenpairs of , i.e.,

 (2.1) Aφk=λkφk,k∈N.

It is known that with

and the eigenvectors

form an orthonormal basis for . Then

 (Alu,v)=∞∑m=1λlm(u,φm)(v,φm),

and we introduce the fractional order spaces

 ˙Hα=dom(Aα2),∥v∥2α:=∥Aα2v∥2=∞∑k=1λαk(v,φk)2,α∈R, v∈˙Hα.

We note that and .

Defining the new variables and , we can write the velocity-displacement form of (1.2) as

 −Δ˙u1+Δu2=0inΩ×(0,T),˙u2−Δu1=finΩ×(0,T),u1=u2=0onΓ×(0,T),u1(⋅,0)=u0, u2(⋅,0)=v0inΩ,

for which, the weak form is to find and such that

 (2.2) a(˙u1(t),v1)−a(u2(t),v1)=0,(˙u2(t),v2)+a(u1(t),v2)=(f(t),v2),∀v1,v2∈V,t∈(0,T),u1(0)=u0, u2(0)=v0.

This equation is used for dG(q) formulation.

## 3. The discontinuous Galerkin time discretization

Here, we apply the dG method in time variable using piecewise polynomials of degree at most , and we investigate the stability.

### 3.1. dG(q) formulation

Let be a temporal mesh with time subintervals and steps , and the maximum step-size by . Let and define the finite element space for each space-time ’Slab’ .

We follow the usual convention that a function is left-continuous at each time level and we define , writing

 U−i,n=Ui(t−n),U+i,n=Ui(t+n),[Ui]n=U+i,n−U−i,nfori=1,2.

The dG method determines on for by setting , and then

 (3.1) ∫In(a(˙U1,V1)−a(U2,V1))dt+a(U+1,n−1,V+1,n−1)=a(U−1,n−1,V+1,n−1),∫In((˙U2,V2)+a(U1,V2))dt+(U+2,n−1,V+2,n−1)=(U−2,n−1,V+2,n−1)+∫In(f,V2)dt,∀V=(V1,V2)∈Pq×Pq.

Now, we define the function space consists of functions which are piecewise smooth with respect to the temporal mesh with values in . We note that . Then we define the bilinear form and the linear form by

 (3.2)

Then , the solution of discrete problem (3.1), satisfies

 (3.3) B(U,V)=L(V),∀V=(V1,V2)∈Vq,U−0=(U−1,0,U−2,0)=(u0,v0).

We note that the solution of (2.2) also satisfies

 (3.4) B(u,v)=L(v),∀v=(v1,v2)∈W,(u1(0),u2(0))=(u0,v0).

These imply the Galerkin’s orthogonality for the error ,

 (3.5)

Integration by parts yields an alternative expression for the bilinear form (3.2), as

 (3.6) B∗(u,v)=N∑n=1∫In{−a(u1,˙v1)−a(u2,v1)−(u2,˙v2)+a(u1,v2)}dt−N−1∑n=1{a(u−1,n,[v1]n)+(u−2,n,[v2]n)}+a(u−1,N,v−1,N)+(u−2,N,v−2,N).
###### Remark 3.1.

We note that the framework applies also to spatial finite dimensional function spaces , such as, a continuous Galerkin finite element method of order for discretization in space variable.

### 3.2. Stability

In this section we present a stability (energy) identity and stability estimate, that are used in a priori error analysis. In our error analysis we need a stability identity for a slightly more general problem, that is such that

 (3.7) B(U,V)=^L(V),∀V∈Vq,

where the linear form is defined by

 ^L((v1,v2))=N∑n=1∫In{a(f1,v1)+(f2,v2)}dt+a(u0,v+1,0)+(v0,v+2,0).

That is, instead of (2.2), we study stability of the dG(q) discretization of a more general problem

 a(˙u1(t),v1)−a(u2(t),v1)=a(f1(t),v1),(˙u2(t),v2)+a(u1(t),v2)=(f(t),v2),∀v1,v2∈V,t∈(0,T),u1(0)=u0, u2(0)=v0.

See Remark 4.2.

We define the -projection by

 ∫In(Pk,nu−u)χdt=0,∀χ∈Pq(In),

and denote .

###### Theorem 3.2.

Let be a solution of (3.7). Then for any and , we have the energy identity

 (3.8) ∥U−1,N∥2l+1+∥U−2,N∥2l+N−1∑n=0{∥[U1]n∥2l+1+∥[U2]n∥2l}=∥u0∥2l+1+∥v0∥2l+2∫T0{a(Pkf1,AlU1)+(Pkf2,AlU2)}dt.

Moreover, for some (independent of ), we have the stability estimate

 (3.9) ∥U−1,N∥l+1+∥U−2,N∥l≤C(∥u0∥l+1+∥v0∥l+∫T0{∥f1∥l+1+∥f2∥l}dt).
###### Proof.

We set in (3.7) to obtain

 N∑n=1∫In{a(˙U1,AlU1)−a(U2,AlU1)+(˙U2,AlU2)+a(U1,AlU2)}dt+N−1∑n=1{a([U1]n,AlU+1,n)+([U2]n,AlU+2,n)}+a(U+1,0,AlU+1,0)+(U+2,0,AlU+2,0)=∫T0{a(f1,AlU1)+(f2,AlU2)}dt+a(u0,AlU+1,0)+(v0,AlU+2,0),

that implies

 12N∑n=1∫In∂∂t∥U1∥2l+1dt+12N∑n=1∫In∂∂t∥U2∥2ldt+N−1∑n=1{a([U1]n,AlU+1,n)+([U2]n,AlU+2,n)}+a(U+1,0,AlU+1,0)+(U+2,0,AlU+2,0)=∫T0{a(Pkf1,AlU1)+(Pkf2,AlU2)}dt+a(u0,AlU+1,0)+(v0,AlU+2,0).

Now writing the first two terms at the left side as

we have

 N−1∑n=1{12∥U−1,n∥2l+1−12∥U+1,n∥2l+1+a([U1]n,AlU+1,n)}+12∥U−1,N∥2l+1+12∥U+1,0∥2l+1+N−1∑n=1{12∥U−2,n∥2l−12∥U+2,n∥2l+([U2]n,AlU+2,n)}+12∥U−2,N∥2l+12∥U+2,0∥2l=N∑n=1∫In{a(Pkf1,AlU1)+(Pkf2,AlU2)}dt+a(U−1,0,AlU+1,0)+(U−2,0,AlU+2,0).

Then, using (for )

 12∥U−1,n∥2l+1−12∥U+1,n∥2l+1+a([U1]n,AlU+1,n)=12∥[U1]n∥2l+1,12∥U−2,n∥2l−12∥U+2,n∥2l+a([U2]n,AlU+2,n)=12∥[U2]n∥2l,

we conclude

 12N−1∑n=1∥[U1]n∥2l+1+12∥U1,N∥2l+1+12∥U+1,0∥2l+1−a(U−1,0,AlU+1,0)+12N−1∑n=1∥[U2]n∥2l+12∥U2,N∥2l+12∥U+2,0∥2l−(U−2,0,AlU+2,0)=∫T0{a(Pkf1,AlU1)+(Pkf2,AlU2)}dt.

Hence, using

 12∥U+1,0∥2l+1−a(Al2U−1,0,Al2U+1,0)=12∥[U1]0]∥2l+1−12∥U−1,0∥2l+1,12∥U+2,0∥2l−(Al2U−2,0,Al2U+2,0)=12∥[U2]0]∥2l−12∥U−2,0∥2l.

We have the identity

 12∥U−1,N∥2l+1+12∥U−2,N∥2l+12N−1∑n=0∥[U1]n∥2l+1+12N−1∑n=0∥[U2]n∥2l=12∥u0∥2l+1+12∥v0∥2l+∫T0{a(Pkf1,AlU1)+(Pkf2,AlU2)}dt.

Finally, to prove the stability estimate (3.9), recalling that all terms on the left side of the stability identity (3.8) are non-negative, we have

 ∥U−1,N∥2l+1+∥U−2,N∥2l≤∥u0∥2l+1+∥v0∥2l+2N∑n=1∫In{a(Pk,nf1,AlU1)+(Pk,nf2,AlU2)}dt,

that, using Couchy-Schwarz inequality, in a classical way we conclude the stability estimate (3.9). ∎

## 4. A priori error estimates

For a given function we define the interpolatant by

 (4.1) Πku(t−n)=u(t−n),forn≥0,∫In(Πku(t)−u(t))χdt=0,forχ∈Pq−1,n≥1,

where the latter condition is not used for . By standard arguments we then have

 (4.2) ∫In∥Πku−u∥jdt≤Ckq+1n∫In∥u(q+1)∥jdt,forj=0,1,

where , see .

### 4.1. Estimates at the nodes

###### Theorem 4.1.

Let and be the solutions of and respectively. Then with and we have

 (4.3) ∥e−1,N∥1+∥e−2,N∥≤CN∑n=1kq+1n∫In{∥u(q+1)2∥1+∥u(q+1)1∥2}dt,
 (4.4) ∥e−1,N∥≤CN∑n=1kq+1n∫In{∥u(q+1)2∥+∥u(q+1)1∥1}dt.
###### Proof.

1. We split the error into two terms, recalling the interpolant

in (4.1),

 e=(e1,e2)=(U1,U2)−(u1,u2)=((U1,U2)−(Πku1,Πku2))+((Πku1,Πku2)−(u1,u2))=(θ1,θ2)+(η1,η2)=θ+η.

We can estimate by (4.2), so we need to find estimates for . Recalling Galerkin’s orthogonality (3.5), we have

 B(θ,V)=−B(η,V),∀V∈Vq.

Then, using the alternative expression (3.6), we have

 B(θ,V)=−B(η,V)=−B∗(η,V)=N∑n=1∫In{(η1,˙V1)+a(η2,V1)+(η2,˙V2)−a(η1,V2)}dt+N−1∑n=1{a(η−1,n,[V1]n)+(η−2,n,[V2]n)}−a(η−1,N,V−1,N)−(η−2,N,V−2,N).

Now, by the fact that vanishes at the time nodes and using the definition of , it follows that and are of degree on and hence they are orthogonal to the interpolation error. We conclude that satisfies the equation

 (4.5) B(θ,V)=∫tN0{a(η2,V1)−(Aη1,V2)}dt.

That is, satisfies (3.7) with and .

2. Then applying the stability estimate (3.9) and recalling , we have

 (4.6)

To prove the first a priori error estimate (4.3), we set . In view of and , we have

 ∥e−1,N∥1+∥e−2,N∥≤C∫T0{∥η2∥1+∥Aη1∥}dt.

Now, using (4.2) and by the elliptic regularity , the first a priori error estimate (4.3) is obtained.

For the second error estimate, we choose in (4.6). In veiw of and , we have

 ∥e−1,N∥+∥e−2,N∥−1≤C∫T0{∥η2∥+∥Aη1∥−1}dt.

Now, using (4.2) and by the fact that , implies the second a priori error estimate (4.4). ∎

###### Remark 4.2.

We note that (4.5), means that and in (3.7), which is the reason for considering an extra load term in the first equation of (2.2). This way, we can balance between the right operators and suitable norms to get optimal order of convergence with minimal regularity requirement on the solution. Indeed, in , it has been proved that the minimal regularity that is required for optimal order convergence for finite element discretization of the wave equation is one extra derivative compare to the order of convergence, and it cannot be relaxed. This means that the regularity requirement on the solution in our error estimates are minimal.

### 4.2. Interior estimates

Now, we prove uniform in time a priori error estimates for dG(0) and dG(1), based on the linearity of the basis functions. We define the following norms

 ∥u∥In=supIn∥u(t)∥,∥u∥IN=sup(0,tN)∥u(t)∥and∥u∥s,In=supIn∥u(t)∥s.
###### Theorem 4.3.

Let and and be the solutions of and , respectively. Then there exists a constant such that for ,

 (4.7) ∥e1∥1,IN+∥e2∥IN≤C(kq+1∥u(q+1)1∥1,IN+kq+1∥u(q+1)2∥IN+N∑n=1kq+2n∥u(q+1)2∥1,In+N∑n=1kq+2n∥u(q+1)1∥2,In),
 (4.8) ∥e1∥IN≤C(kq+1∥u(q+1)1∥IN+N∑n=1kq+2n∥u(q+1)2∥In