SIAM J. NUMER. ANAL.
c
2006 Society for Industrial and Applied Mathematics
Vol. 44, No. 5, pp. 2131–2158
OPTIMAL DISCONTINUOUS GALERKIN METHODS
FOR WAVE PROPAGATION
∗
ERIC T. CHUNG
†
AND BJ
¨
ORN ENGQUIST
‡
Abstract. We have developed and analyzed a new class of discontinuous Galerkin methods (DG)
which can be seen as a compromise between standard DG and the finite element (FE) method in the
way that it is explicit like standard DG and energy conserving like FE. In the literature there are
many methods that achieve some of the goals of explicit time marching, unstructured grid, energy
conservation, and optimal higher order accuracy, but as far as we know only our new algorithms
satisfy all the conditions. We propose a new stability requirement for our DG. The stability analysis
is based on the careful selection of the two FE spaces which verify the new stability condition. The
convergence rate is optimal with respect to the order of the polynomials in the FE spaces. Moreover,
the convergence is described by a series of numerical experiments.
Key words. discontinuous Galerkin, wave propagation, optimal rate of convergence
AMS subject classifications. 65M12, 65M15, 65M60, 78M10
DOI. 10.1137/050641193
1. Introduction. Many applications involve the solution of wave equations. Ex-
amples are electromagnetic waves for radar and communication as well as acoustic
and seismic wave propagation. Let Ω ⊂ R
2
be a two dimensional polygonal domain
with outward normal vector n and let T>0 be a fixed time. Given two positive
constants a
1
> 0, a
2
> 0 and two given functions F
1
(x, t), F
2
(x, t), we will consider,
for (x, t) ∈ Ω ×(0,T), the following wave propagation problem: find a function u(x, t)
and a vector field v(x, t) ∈ R
2
such that
a
1
∂u
∂t
+ Bv = F
1
,(1.1)
a
2
∂v
∂t
− B
∗
u = F
2
,(1.2)
where the two operators B and B
∗
satisfy
Ω
0
(B
∗
φ)ψdx−
Ω
0
(Bψ)φdx=
∂Ω
0
(Lψ)φdσ(1.3)
for all subset Ω
0
⊂ Ω. Here L is some operator depending on the two operators B,
B
∗
and the subdomain Ω
0
. We also denote by L
⊥
the operator such that |ψ|
2
=
|Lψ|
2
+ |L
⊥
ψ|
2
. Assume that there is an operator B
⊥
such that BB
⊥
p = 0 for all p.
Acting (B
⊥
)
∗
to (1.2) and using (B
⊥
)
∗
B
∗
= 0, we have the following:
∂
∂t
(a
2
(B
⊥
)
∗
v)=(B
⊥
)
∗
F
2
.(1.4)
∗
Received by the editors September 26, 2005; accepted for publication (in revised form) April 6,
2006; published electronically November 3, 2006.
http://www.siam.org/journals/sinum/44-5/64119.html
†
Applied and Computational Mathematics, California Institute of Technology, Pasadena, CA
91125 (tschung@acm.caltech.edu).
‡
Department of Mathematics, University of Texas at Austin, Austin, TX 78712 (engquist@math.
utexas.edu).
2131
2132 ERIC T. CHUNG AND BJ
¨
ORN ENGQUIST
The condition (1.4) usually has important physical significance. For example, in the
case of electromagnetic wave propagation (see (E) in the following), v represents the
electric field, (B
⊥
)
∗
is the divergence operator, and (1.4) is just the continuity equation
expressing the conservation of charges. We supplement the system (1.1)–(1.2) with
boundary condition
Lv =0 ∀x ∈ ∂Ω(1.5)
and initial conditions
u(x, 0) = u
0
(x) and v(x, 0) = v
0
(x) ∀x ∈ Ω,(1.6)
where u
0
(x) and v
0
(x) are given. In particular, we are interested in the acoustic and
the electromagnetic wave equations, which correspond to the following choice of B
and B
∗
:
(A) Acoustic: Bv = −∇ · v, B
∗
u = ∇u, and Lv = v · n.
(E) Electromagnetic: Bv = ∇×v, B
∗
u = ∇×u, and Lv = v × n.
Notice that, in (1.1)–(1.2), u(x, t) is a function while v(x, t)=(v
1
(x, t),v
2
(x, t)) is a
vector having two components. So the operator ∇× is defined as ∇×v = ∂
1
v
2
− ∂
2
v
1
for any vector field v and ∇×u =(∂
2
u, −∂
1
u) for any function u. Here v × n =
v
1
n
2
−v
2
n
1
with n =(n
1
,n
2
). In (A) and (E), we use n to denote generically the unit
normal vector of the corresponding subdomain which defines L and L
⊥
. Moreover,
we have L
⊥
ψ = ψ × n for (A) while L
⊥
ψ = ψ · n for (E). Furthermore, we have
B
⊥
p = ∇×p and (B
⊥
)
∗
p = ∇×p for (A) while B
⊥
p = ∇p and (B
⊥
)
∗
p = −∇ · p for
(E). Wave propagation problems can be solved by partial differential equation (PDE)
techniques, integral equation techniques, and asymptotic techniques. Among PDE
techniques, finite difference (FD) method, finite volume (FV) method, finite element
(FE) method, and discontinuous Galerkin (DG) method are the most popular choices.
The FD method provides a simple way to solve wave propagation problems, but it
is typically low order and applies only to structured grids. The FV method can be
seen as a generalization of the FD method to unstructured grids, but it is still low
order. The FE and DG methods provide high order solvers for the time dependent
wave equations on unstructured grids.
N´ed´elec [15] introduces a curl-conforming FE method for solving Maxwell’s equa-
tions. Geveci [8] proposed a mixed FE method for the scalar wave equation. The
inversion of the mass matrix at each time step causes some possible drawback in the
efficiency of those methods. Mass lumping techniques can be used to avoid solving
linear systems. In Cohen and Monk [6], a mass lumping method for rectangular grids
is developed. In B´ecache, Joly, and Tsogka [1], a new class of mixed FE method, which
is suitable for mass lumping, is developed for the scalar wave equation. Cohen, Joly,
Torjman, and Roberts [5] design a mass lumping technique for triangular grids for
polynomial order up to five. Discontinuous Galerkin methods provide explicit schemes
in the sense that only block diagonal mass matrices have to be inverted. Hesthaven and
Warburton [10] proposed a DG method based on upwind flux and Cockburn, Li, and
Shu [4] proposed a DG method based on locally divergence free basis and upwind flux.
While the schemes are successful, energy is not conserved due to the upwinding. Fe-
zoui, Lanteri, Lohrengel, and Piperno [7] proposed a DG method based on central flux.
This method preserves energy, but the convergence rate of the scheme is suboptimal.
Recently, a new DG method has been developed for the wave equation in second order
form; see Grote, Schneebeli, and Sch¨otzau [9]. The method is also energy conserving
OPTIMAL DISCONTINUOUS GALERKIN METHODS 2133
DG
both u and v are
discontinuous
Our DG
u and v are continuous
at different points
v
u
FE
both u and v are
continuous
Fig. 1.1. Comparison among standard DG, our new DG, and FE methods.
in the sense of a newly defined energy. A space-time DG method has also been
developed in Monk and Richter [14].
In this paper, we will develop and analyze a new class of DG methods which can be
seen as a compromize between FE and DG methods. Our new DG method combines
the advantages of FE and DG methods in the sense that it is both energy conserving
and explicit. The idea is to use discontinuous functions with extra continuity. In
the velocity-potential formulation of the scalar wave equation (A), we will add extra
continuity to the velocity where the potential is discontinuous and add extra continuity
to the potential where the velocity is discontinuous. For Maxwell’s equations (E), a
similar idea can be applied to the electric and magnetic fields. As a result, the flux
integrals are evaluated exactly, which is the basis of energy conservation. However, the
addition of the extra continuity cannot be done arbitrarily due to stability concerns.
It has to be done in such a way that some inf-sup conditions are satisfied. In Figure
1.1, we illustrate this idea in one space dimension. For standard DG, both unknown
functions u and v, which are velocity and potential for scalar wave equation and
are electric and magnetic fields for Maxwell’s equations, are discontinuous at cell
boundaries. For FE methods, both u and v are continuous. For our new DG, the two
functions are continuous at different points.
Yee’s scheme [16] has been a very popular numerical method for computational
electromagnetics. It is a second order central FD method on structured grids. The
success of the scheme is due to the use of a staggered grid. Our new DG method is
a FE method on staggered grids and can be seen as a higher generalization of Yee’s
scheme on unstructured grids. In particular, in one space dimension, our new DG
method with piecewise constant approximation is the same as Yee’s scheme. In two
space dimensions, our new DG method in the lowest order is some averaged version
of Yee’s scheme.
The rest of the paper is organized as follows. In section 2, we will introduce
the new FE spaces and prove the corresponding unisolvence and interpolation error
estimates. The new DG is then derived in section 3. In section 4, under the assumption
of some inf-sup conditions, the stability and convergence of the method are proved.
The inf-sup conditions are then verified in section 5. Furthermore, some numerical
experiments are presented in section 6. The paper ends with a conclusion.
Remark. We consider only two space dimensions in this paper. For three space
dimensions, a careful choice of the two FE spaces U
h
and V
h
that verify (3.1) and
2134 ERIC T. CHUNG AND BJ
¨
ORN ENGQUIST
(3.2) as well as the two inf-sup conditions (4.1)–(4.2) are required. This work will be
developed in a forthcoming paper.
2. FE spaces. Assume the domain Ω is triangulated by a family of triangles T
so that Ω = ∪{τ | τ ∈T}. Let τ ∈T. We define h
τ
as the diameter of τ and ρ
τ
as the
supremum of the diameters of the circles inscribed in τ . The mesh size h is defined
as h = max
τ∈T
h
τ
. We will assume the set of triangles T forms a regular family of
triangulation of Ω so that there exist a uniform constant K independent of the mesh
size such that [3]
h
τ
≤ Kρ
τ
∀τ ∈T.
In addition, we will assume the triangulation satisfies the inverse assumption [3].
Let E be the set of all edges and let E
0
⊂E be the set of all interior edges of the
triangles in T . The length of σ ∈E will be denoted by h
σ
. We also denote by N the
set of all interior nodes of the triangles in T . Here, by interior edge and interior node,
we mean any edge and node that does not lie on the boundary ∂Ω. Let ν ∈N.We
define
S(ν)=∪{τ ∈T |ν ∈ τ}.(2.1)
That is, S(ν) is the union of all triangles having vertex ν. We will assume the
triangulation of Ω satisfies the following condition.
Assumption on triangulation: There exists a subset N
1
⊂N such that
(A1) Ω = ∪{S(ν) | ν ∈N
1
}.
(A2) S(ν
i
) ∩S(ν
j
) ∈E
0
for all distinct ν
i
,ν
j
∈N
1
.
Let ν ∈N
1
. We define
E
u
(ν)={σ ∈E|ν ∈ σ}.(2.2)
That is, E
u
(ν) is the set of all edges that have ν as one of their endpoints. We further
define
E
u
= ∪{E
u
(ν) | ν ∈N
1
} and E
v
= E\E
u
.(2.3)
Notice that E
u
contains only interior edges since one of the endpoints of edges in E
u
has a vertex from N
1
. On the other hand, E
v
has both interior and boundary edges.
So, we also define E
0
v
= E
v
∩E
0
which contains elements from E
v
that are interior
edges. Notice that we have E
v
\E
0
v
= E∩∂Ω. Furthermore, for σ ∈E
0
v
, we will let
R(σ) be the union of the two triangles sharing the same edge σ.Forσ ∈E
v
\E
0
v
,we
will let R(σ) be the only triangle having the edge σ.
In practice, triangulations that satisfy assumptions (A1)–(A2) are not difficult
to construct. In Figure 2.1, we illustrate how this kind of triangulation is generated.
First, the domain Ω is triangulated by a family of triangles, called
˜
T . Each triangle in
this family is then subdivided into three subtriangles by connecting a point inside the
triangle with its three vertices. Then we define the union of all these subtriangles to
be our triangulation T . Each triangle in
˜
T corresponds to an S(ν) for some ν inside
the triangle. In Figure 2.1, we show two of the triangles, enclosed by solid lines, in
this family
˜
T . This corresponds to 6 triangles in the triangulation T . The dotted
lines represent edges in the set E
u
while solid lines represent edges in the set E
v
.
Lemma 2.1. Each τ ∈T has exactly two edges that belong to E
u
.
OPTIMAL DISCONTINUOUS GALERKIN METHODS 2135
•
•
S(ν
1
)
S(ν
2
)
ν
1
ν
2
Fig. 2.1. Triangulation.
Proof. First of all, τ has at least one interior vertex. We will show that there is
exactly one vertex of τ that belongs to N
1
. If none of the three vertices of τ belong
to N
1
, then τ
0
∩S(ν) is an empty set for all ν ∈N
1
, where τ
0
is the interior of τ .
Then, ∪{S(ν) | ν ∈N
1
}∩τ
0
is an empty set. So, ∪{S(ν) | ν ∈N
1
} = Ω, which
violates assumption (A1). If τ has two vertices, ν
i
and ν
j
, that belong to N
1
, then
S(ν
i
) ∩S(ν
j
) contains τ . So, it violates assumption (A2). The case that τ has all
vertices belonging to N
1
can be discussed in the same way. In conclusion, τ has
exactly one vertex which belongs to N
1
. So, by the definition of E
u
, the two edges
having the vertex in N
1
belong to E
u
.
Given τ ∈T, we will denote by ν(τ)
1
, ν(τ )
2
, and ν(τ )
3
the three vertices of
τ. Moreover, ν(τ)
1
is the vertex that is one of the endpoints of the two edges of τ
that belong to E
u
. Then ν(τ)
2
and ν(τ)
3
are named in a counterclockwise direction.
In addition, λ
τ,1
(x), λ
τ,2
(x), and λ
τ,3
(x) are the barycentric coordinates on τ with
respect to the three vertices ν(τ )
1
, ν(τ )
2
, and ν(τ )
3
.
Now, we will discuss the FE spaces. Let k ≥ 0 be a nonnegative integer. Let
τ ∈T. We define P
k
(τ) as the space of polynomials of degree less than or equal to k
on τ. We also define
R
k
(τ)=P
k
(τ) ⊕
˜
P
k+1
(τ),(2.4)
where
˜
P
k+1
(τ) is the space of homogeneous polynomials of degree k+1 on τ in the two
variables λ
τ,2
and λ
τ,3
such that the sum of the coefficients of λ
k+1
τ,2
and λ
k+1
τ,3
is equal to
zero. That is, any function in
˜
P
k+1
(τ) can be written as
i+j=k+1,i≥0,j≥0
a
i,j
λ
i
τ,2
λ
j
τ,3
such that a
k+1,0
+ a
0,k+1
= 0. Now, we define
U
h
= {φ | φ|
τ
∈ R
k
(τ); φ is continuous at the k + 1 Gaussian points of σ ∀σ ∈E
u
}.
For any edge σ, we use P
k
(σ) to represent the space of one dimensional polynomials
of degree less than or equal to k on σ. We define the following degrees of freedom:
(UD1) For each edge σ ∈E
u
,wehave
σ
φp
k
dσ
for all p
k
∈ P
k
(σ).
(UD2) For each triangle τ ∈T,wehave
τ
φp
k−1
dx
for all p
k−1
∈ P
k−1
(τ) (for k ≥ 1).
