L
2
Stable Discontinuous Galerkin Methods for One-dimensional
Two-way Wave Equations
Yingda Cheng
∗
, Ching-Shan Chou
†
, Fengyan Li
‡
, and Yulong Xing
§
July 15, 2015
Abstract
Simulating wave propagation is one of the fundamental problems in scientific computing.
In this paper , we consider one-dimensional two-way wave equations, and investigate a family
of L
2
stable high order discontinuous Galerkin methods defined through a general form of
numerical fluxes. For these L
2
stable methods, we systema tically establish stability (hence
energy conser va tion), error e stimates (in both L
2
and negative-order norms), and dispersion
analysis. One novelty of this work is to identify a sub-family of the numerical fluxes, termed
αβ-flux e s. Discontinuous Galerkin methods with αβ-fluxes are proven to have optimal L
2
error
estimates and superconvergence properties. Moreover, both the upwind and alternating fluxes
belong to this sub-family. Dispersion a nalysis, which examines both the physical and spur ious
modes, provides insights into the sub-optimal accuracy of the methods using the central flux
and the odd degree polynomials, and demonstrates the importance of numerical initialization
for the proposed non-dissipative schemes. Numerical examples are presented to illustrate the
accuracy and the long-term behavior of the methods under consideration.
1 Introduction
Wave propagation is a fundamental form of energy transmission, which arises in many fields of
science, engineering and industry, and it is significant to geoscience, petroleum engineering and
electromagnetics. A vast amount of research has been done for wave simulations, and the com-
monly used numerical methods range from finite d ifference, finite volume to spectral element and
finite element methods ([18, 17, 24, 15] and references therein, [35, 4, 30]). Among various numer-
ical methods, each with their own advantages, here we will confine our attention to discontinuous
Galerkin (DG) methods. DG method s belong to a class of finite element methods using piecewise
polynomial spaces for both the numerical solution and the test function. They were originally
devised to solve hyperbolic conservation laws, e.g. [26, 12, 14]. The methods can be easily de-
signed to have arbitrary order of accuracy. They are flexible with unstructured meshes, and are
∗
Department of Mathematics, Michigan State Un iversity, East Lansing, MI 48824. Email: ycheng@math.msu.edu.
Research is supported by NSF grants DMS-1217563 and DMS-1318186.
†
Department of Mathematics, The Ohio State University, Columbus, OH 43210. Email: chou@math.osu.edu.
Research is supported by NSF grants DMS-1020625 and DMS-1253481.
‡
Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180. E-mail: lif@rpi.edu.
Supported in part by NSF grants DMS-0847241 and DMS-1318409.
§
Computer Science and Mathematics Division, Oak Ridge Nationalist Laboratory, Oak R idge, TN 37831 and
Department of Mathematics, University of Tennessee, Knoxville, TN 37996. E-mail: xingy@math.utk.edu. Research
is sponsored by NSF grant DMS-1216454, ORNL and the U. S. Department of Energy, Office of Advanced Scientific
Computing Research. The work was partially performed at ORNL, which is managed by UT-Battelle, LLC, under
Contract No. DE-AC05-00OR22725.
1
natural candidates for h-p adap tivity. These methods are known to be highly efficient in parallel
computation, due to the compact stencils. Importantly, DG methods perform well in long-term
wave s imulations [21, 1, 8], given their excellent dispersive and dissipative properties.
In this paper, we consider the one-dimensional linear two-way wave problem,
E
t
= B
x
− S
1
, (1.1a)
B
t
= E
x
− S
2
, (1.1b)
where E = E(t, x) and B = B(t, x) are unknown functions, and S
1
= S
1
(t, x) and S
2
= S
2
(t, x)
are both source terms. Note that the system (1.1) is equivalent to the second-order wave equation.
Moreover, Maxwell’s equations can be viewed as a special case when S
2
= 0. As a hyperbolic
system, (1.1) can be compu ted by DG methods directly, yet the key properties of the methods
lie in the choices of th e numerical fluxes. For example, it is kn own that the use of either an
alternating flux or central fl ux will give an energy preserving scheme, while the use of th e upwind
flux will result in a decreasing discrete energy. Moreover, the L
2
error bounds with alternating
and upw ind fluxes are optimal, while with central flux, the accuracy may be sub-optimal. This
leads to the question that, whether we can find a general principle of selecting num erical fluxes,
and furthermore, understand the accuracy, stability, energy conservation, and other properties of
the resulted numerical methods.
The goal of this paper is to perform a detailed and systematic investigation for a general family
of L
2
stable DG methods. Note that th e L
2
norm of (E, B) for system (1.1) is equivalent to the
total energy, which many DG methods attempt to capture [16, 20, 10, 32, 9]. Th e proposed methods
are defined through a group of numerical fluxes which are certain linear combinations of jumps and
averages of the numerical solutions at the cell interfaces with three parameters α, β
1
, β
2
involved.
In a previous work by Ainsworth and colleagues [3], the same family of fluxes were considered with
emphasis on dispersion analysis. For these L
2
stable methods, we will establish stability and hence
the energy conservation, error estimates in both L
2
and negative-order norms, superconvergence,
and dispersion analysis.
One novelty of this work is that, in search of DG methods with optimal L
2
error estimate, we
identify a sub-family of the numerical fluxes, termed αβ-fluxes, that are d etermined by a specific
relation among α, β
1
, β
2
. This relation, which is satisfied by the widely known upwind and alter-
nating fluxes, was introduced in [3] to characterize different modes in dispersion analysis for DG
methods solving two-dimensional second-order wave equations. The relation is used in this work
to design a new projection op erator, a key component in our proof for the optimal error estimates .
Besides the optimal L
2
accuracy, we further prove superconvergence properties for the DG methods
with αβ-fluxes by following the analysis in [34] for the linear advection equation. Such supercon-
vergence seems to be uniquely enjoyed by DG methods associated with the αβ-fluxes as suggested
numerically.
For the proposed general L
2
stable DG methods, we also systematically perform the dispersion
analysis, and present the negative-order norm error estimates as well as the related post-processing
techniques similar to those in [13, 23]. For long time wave simulations, it is important to understand
the dispersive and dissipative properties of the numerical methods. Our dispersion analysis, which
takes a different viewpoint from [3], examines b oth the physical and spurious modes. In particular,
it gives insights into the sub-optimal accuracy of DG methods with the central flux and the odd
degree polynomials, and demonstrates the importance of numerical initialization for the p roposed
non-dissipative schemes. Related work on the dispersion analysis of semi-discrete or fully discrete
DG methods in literature can be found in [22, 28, 21, 1, 3, 27, 33, 2].
The remaining of this paper is organized as follows. In Section 2, we introduce a general
family of L
2
stable DG methods for one-dimensional two-way wave equations, and also define a
2
sub-family of the methods associated with αβ-fluxes. We provide in Section 2.1 the analysis of L
2
stability and energy conserving property, and this is followed by L
2
error estimates in Section 2.2
for th e general L
2
stable DG methods. In S ection 2.3, we present the superconvergence results and
postprocessing techniques. The dispersion analysis is performed in Section 2.4. Section 3 contains
numerical examples to illustrate the performance of the proposed m ethods, and we conclude with
a few remarks in Section 4.
2 L
2
stable DG methods
In this section, we will formulate a general family of semi-discrete DG methods which is L
2
stable for
the one-dimensional two-way wave equations (1.1). Here we consider periodic boundary conditions
for simplicity. We start with a mesh of the computational domain Ω = [a, b], a = x
1
2
< x
3
2
< ··· <
x
N+
1
2
= b. Each cell is denoted as I
j
= [x
j−
1
2
, x
j+
1
2
], with its center x
j
=
1
2
(x
j−
1
2
+ x
j+
1
2
) and
the length h
j
= x
j+
1
2
− x
j−
1
2
. Let h = max
1≤j≤N
h
j
. The mesh is assumed to be quasi-uniform,
namely, there exists a positive constant σ, such that
h
min
j
h
j
< σ, as the mesh is refined. We now
define a finite dimensional discrete space,
V
r
h
= {v : v|
I
j
∈ P
r
(I
j
), j = 1, 2, ··· , N}, (2.1)
which consists of piecewise polynomials of degree up to r with respect to the mesh. Note that
functions in V
r
h
are allowed to have discontinuities across element interfaces. For any v ∈ V
r
h
, we
denote by v
+
j+
1
2
and v
−
j+
1
2
the limit values of v at x
j+
1
2
from the right cell I
j+1
and fr om the left
cell I
j
, resp ectively. We use th e usual notation [v]
j+
1
2
= v
+
j+
1
2
−v
−
j+
1
2
and {v}
j+
1
2
=
1
2
(v
+
j+
1
2
+ v
−
j+
1
2
)
to represent the jump and the average of th e function v at x
j+
1
2
for any j.
The semi-discrete DG method for the sys tem (1.1) is formulated as follows: find E
h
(t, ·),
B
h
(t, ·) ∈ V
r
h
, su ch that
Z
I
j
(E
h
)
t
φdx +
Z
I
j
B
h
φ
x
dx − (F
B
(B
h
, E
h
)φ
−
)
j+
1
2
+
F
B
(B
h
, E
h
)φ
+
j−
1
2
=
Z
I
j
S
1
φdx, (2.2a)
Z
I
j
(B
h
)
t
ψdx +
Z
I
j
E
h
ψ
x
dx − (F
E
(E
h
, B
h
)ψ
−
)
j+
1
2
+ (F
E
(E
h
, B
h
)ψ
+
)
j−
1
2
=
Z
I
j
S
2
ψdx, (2.2b)
for all test functions φ, ψ ∈ V
r
h
, and for all j. By summing up the two equations in (2.2) over
all mesh elements, we can write the DG method in a more compact form. We look for E
h
(t, ·),
B
h
(t, ·) ∈ V
r
h
, su ch that
a
h
(E
h
, B
h
; φ, ψ) = S(φ, ψ), ∀φ, ψ ∈ V
r
h
, (2.3)
where
a
h
(E
h
, B
h
; φ, ψ) =
Z
Ω
(E
h
)
t
φdx +
X
j
Z
I
j
B
h
φ
x
dx + (F
B
(B
h
, E
h
)[φ])
j−
1
2
!
+
Z
Ω
(B
h
)
t
ψdx +
X
j
Z
I
j
E
h
ψ
x
dx + (F
E
(E
h
, B
h
)[ψ])
j−
1
2
!
, (2.4)
and S(φ, ψ) =
R
Ω
S
1
φ + S
2
ψdx.
3
Both the terms F
B
and F
E
in (2.2) are numerical fluxes, and they are single-valued fu nctions
defined on the cell interfaces and should be designed to ensure numerical s tability and accuracy. In
the present work, we consider the following numerical fluxes,
F
B
(B
h
, E
h
) = {B
h
} + α[B
h
] + β
1
[E
h
], (2.5a)
F
E
(E
h
, B
h
) = {E
h
} − α[E
h
] + β
2
[B
h
]. (2.5b)
Here α, β
1
, β
2
are constants that are taken to be O(1), with β
1
and β
2
being non-negative. These
numerical fluxes were considered in [3], and they are consistent, that is,
F
B
(B, E) = B, F
E
(E, B) = E. (2.6)
The DG methods with fluxes (2.5) define a very general family of L
2
stable DG methods for the
system (1.1). Note that the numerical flu x es (2.5) include several commonly used ones in literature.
For examples, when α = 0, β
1
= β
2
=
1
2
, we have the upwind flux; when α = β
1
= β
2
= 0, we have
the central flux, and the alternating flux is obtained when α = ±
1
2
, β
1
= β
2
= 0.
One novelty of this work is that we further identify a su b-family of the numerical fluxes (2.5),
named αβ-fluxes, and the corresponding DG methods have some important provable resu lts in
terms of L
2
error estimates and superconvergence, which will be carried out in Sections 2.2 and
2.3.
Definition 2.1. An αβ-flux is a numerical flux (2.5) when α and β
i
≥ 0 (i = 1, 2) satisfy
α
2
+ β
1
β
2
=
1
4
. (2.7)
When β
1
= β
2
= β, an αβ-flux can be determined by a single parameter α. (Here β is non-
negative and β =
q
1
4
− α
2
.) It is easy to see that both the upwind and alternating fluxes are
special cases of this αβ-flux family.
Remark 2.1. One can further generalize the numerical flux in (2.5) as follows,
F
B
(B
h
, E
h
) = {B
h
} + α
1
[B
h
] + β
1
[E
h
], (2.8a)
F
E
(E
h
, B
h
) = {E
h
} + α
2
[E
h
] + β
2
[B
h
], (2.8b)
which involves four parameters. In this work, we will only present analysis for the DG methods with
(2.5), and will summarize the L
2
stability and error estimates for the more g eneral DG methods
with (2.8) in Remark 2.3 and Remark 2.8. The parameters β
1
and β
2
are chosen to be non-negative
to ensure that the proposed DG methods are stable (see Theorem 2.2). In addition, there is no
benefit in terms of accuracy or stability if we allow α, β
1
and β
2
to be more general than O(1).
2.1 L
2
stability and energy conservation
In this section, we will establish the L
2
stability, which also informs about the energy conservation
property, for th e semi-discrete DG method with the general numerical flux (2.5). It suffices to
consider S
1
= S
2
= 0.
Theorem 2.2. With S
1
= S
2
= 0, the semi-discrete DG scheme (2.2) (or (2.3) ) with the numerical
flux (2.5) and β
i
≥ 0, i = 1, 2, satisfies
d
dt
E
h
(t) = −
X
j
β
1
[E
h
]
2
+ β
2
[B
h
]
2
j−
1
2
≤ 0, (2.9)
4
where
E
h
(t) =
1
2
Z
Ω
(E
h
(t, x))
2
+ (B
h
(t, x))
2
dx
is the energy of the system (1.1) at time t.
Proof. We fi rst introduce
H
(1)
j
(B
h
, E
h
; φ) = −
Z
I
j
B
h
φ
x
dx + (F
B
(B
h
, E
h
)φ
−
)
j+
1
2
− (F
B
(B
h
, E
h
)φ
+
)
j−
1
2
,
H
(2)
j
(E
h
, B
h
; ψ) = −
Z
I
j
E
h
ψ
x
dx + (F
E
(E
h
, B
h
)ψ
−
)
j+
1
2
− (F
E
(E
h
, B
h
)ψ
+
)
j−
1
2
.
With periodic boundary conditions and the specific definition of the general numerical flux in (2.5),
as well as the identity [φψ] = {ψ}[φ] + {φ}[ψ], the following holds for any φ, ψ ∈ V
r
h
,
X
j
H
(1)
j
(ψ, φ; φ) + H
(2)
j
(φ, ψ; ψ) =
X
j
[φψ] − F
E
(φ, ψ)[ψ] −F
B
(ψ, φ)[φ]
j−
1
2
= −
X
j
β
1
[φ]
2
+ β
2
[ψ]
2
j−
1
2
. (2.10)
Using the definition of a
h
in (2.4), one further has
a
h
(φ, ψ; φ, ψ) =
Z
Ω
(φ
t
φ + ψ
t
ψ)dx −
X
j
H
(1)
j
(ψ, φ; φ) + H
(2)
j
(φ, ψ; ψ)
=
1
2
d
dt
Z
Ω
(φ
2
+ ψ
2
)dx +
X
j
β
1
[φ]
2
+ β
2
[ψ]
2
j−
1
2
. (2.11)
Now in the semi-discrete DG method with S
1
= S
2
= 0, we take φ = E
h
in (2.2a) and ψ = B
h
in
(2.2b), and get a
h
(E
h
, B
h
; E
h
, B
h
) = 0. This, combined with the general result in (2.11), gives the
L
2
stability in (2.9).
Note that all flux choices with β
i
≥ 0, i = 1, 2, produce L
2
stable numerical solutions. I n
particular, the semi-discrete DG method with either the central or alternating, or the more general
flux (2.5) with β
1
= β
2
= 0, preserves the energy of the system. On the other hand, with the
commonly used u pwind flux (α = 0, β
1
= β
2
=
1
2
), the L
2
energy decays with time, as exp ected.
Remark 2.3. For the source free problem, it can be shown that the semi-discrete DG scheme (2.2)
(or (2.3)) with the more general numerical flux (2.8) and β
i
≥ 0, i = 1, 2, (α
1
+ α
2
)
2
≤ 4β
1
β
2
,
satisfies
d
dt
E
h
(t) = −
X
j
β
1
[E
h
]
2
+ β
2
[B
h
]
2
+ (α
1
+ α
2
)[E
h
][B
h
]
j−
1
2
≤ 0.
2.2 L
2
error estimates
In this section, we will establish error estimates in the L
2
-norm for the semi-discrete DG schemes
up to a given time T < ∞ with various choices of numerical fluxes. The following projections,
defined from H
r+1
(Ω) onto V
r
h
, will be used in the analysis.
5