Discontinuous Galerkin Method

Praveen Chandrashekar

ATM Workshop on PDE and Mechanics

Kerala School of Mathematics, Kozhikode
4-6 February, 2016
[Link]
Contents

Chapter One: Scalar conservation law

1 Scalar conservation law . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Linear convection equation . . . . . . . . . . . . . . . . . . . . . 1
1.2 Burger’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Weak solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Rankine-Hugoniot condition . . . . . . . . . . . . . . . . . . . . 3
3 Kruzkov’s result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

Chapter Two: DG scheme in 1-D

1 Mesh and approximation space . . . . . . . . . . . . . . . . . . . . . 5
2 Semi-discrete DG scheme . . . . . . . . . . . . . . . . . . . . . . . . . 6
3 Relation to finite volume scheme . . . . . . . . . . . . . . . . . . . . 7
4 Conservation property . . . . . . . . . . . . . . . . . . . . . . . . . . 7

Chapter Three: Numerical flux

1 Properties of numerical flux . . . . . . . . . . . . . . . . . . . . . . . 9
2 Riemann solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3 Numerical flux: Linear convection equation . . . . . . . . . . . . . . 10
4 Numerical flux: Non-linear conservation law . . . . . . . . . . . . . . 10
4.1 Godunov flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.2 Lax-Friedrich flux . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.3 Rusanov flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.4 Roe flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Chapter Four: Stability

1 Linear convection equation . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1 Periodic boundary conditions . . . . . . . . . . . . . . . . . . . . 14
1.2 Dirichlet condition . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Non-linear conservation law . . . . . . . . . . . . . . . . . . . . . . . 15
2.1 Entropy condition . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Chapter Five: Error estimates

1 Error estimate: ut + cux = 0, semi-discrete scheme . . . . . . . . . . 19
Chapter Six: Basis functions
1 Nodal basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2 Taylor basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Orthogonal polynomials (Modal approach) . . . . . . . . . . . . . . . 23

Chapter Seven: Implementation in 1-D

1 ODE system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2 Quadrature rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3 Setting initial condition . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5 Strong stability preserving RK schemes . . . . . . . . . . . . . . . . . 28
6 CFL condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
7 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
8 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Chapter Eight: Limiters and TVD property

1 Limiter for DG scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2 Limiters: Implementation . . . . . . . . . . . . . . . . . . . . . . . . 38
3 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4 TVB Limiter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Chapter A: Quadrature rules

1 Quadrature in 1-D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

∗
Starred subsections are optional.
 Chapter One
Scalar conservation law

1 Scalar conservation law

A scalar conservation law is a partial differential equation of the form
ut + f(u)x = 0
where
• u is called the conserved variable

• f(u) is the flux of u

Such an equation is a conservation law because, if we integrate the equation on any
interval [a, b], we obtain
Z
d b
u(x, t)dx + f(u(b, t)) − f(u(a, t)) = 0
dt a
This says that
Rate of change of total u in [a, b] = Net flux of u into [a, b]
If the net flux is zero, then the total quantity is conserved. Moreover, these are
hyperbolic conservation laws since their solutions have wave-like properties.

1.1 Linear convection equation

The flux is a linear function of the conserved variable, f(u) = cu, where c = constant
ut + cux = 0, u(x, 0) = u0 (x)
The exact solution is given by
u(x, t) = u0 (x − ct)

u0 (x)
u(x, t)

ct
The initial condition is transported with velocity c without change of form.
2 Chapter One. Scalar conservation law

1.2 Burger’s equation

If the flux is of the form f(u) = 12 u2 , then the conservation law is

1
u t + ( u 2 )x = 0 or ut + uux = 0 if u is smooth
2
As long as the solution is smooth, it is given by

u(x, t) = u0 (x − u(x, t)t)

which is an implicit solution. In many cases, even if the initial condition is infinitely
smooth, the solution at future times cannot remain smooth and discontinuities can
develop in finite time.
Solution at times t=0, t=tc and t>tc
2.0 2.0

1.5 1.5

1.0 1.0

0.5 0.5

0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5

2 Weak solution

Take a smooth test function with compact support, say φ ∈ C10 (R × R+ )

Z∞ Z  
∂u ∂f
+ φdxdt = 0
0 R ∂t ∂x

t Supp(φ)

a O b x

Perform an integration by parts on both terms in the equation so that the derivatives
are transfered to the test function. This motivates the following definition.
Section 3. Kruzkov’s result 3

1.1 Definition (Weak solution) A function u : R × R+ → R is a weak solution of the

IVP
ut + f(u)x = 0, (x, t) ∈ R × R+ , u(x, 0) = u0 (x)
together with locally integrable initial data u0 , if u is locally integrable and satisfies
Z ∞
∞ Z Z
∞

(uφt + f(u)φx )dxdt + u0 (x)φ(x, 0)dx = 0, ∀φ ∈ C10 (R × R+ )

0 −∞ −∞

1.2 Lemma (Classical solution) Let u ∈ C1 (R × R+ ) be a weak solution. Then it is a

classical solution.

2.1 Rankine-Hugoniot condition

Suppose there is a discontinuity at x(t) and let

d
s= x(t) = shock speed
dt
Then the two values u(x− (t), t) and u(x+ (t), t) satisfy the RH condition

f(u(x+ (t), t)) − f(u(x− (t), t)) = s(u(x+ (t), t) − u(x− (t), t))

1.1 Definition (Weak solution) A weak solution u is a piecewise smooth solution which
satisfies the RH condition at the points of discontinuity of u.

3 Kruzkov’s result

The scalar Cauchy problem for

ut + f(u)x = 0, f ∈ C1 (R)

with initial condition

u(0, x) = u0 (x), u0 ∈ L∞ (R)

has a unique entropy solution

u ∈ L∞ (R+ × R)

which fulfills (important for numerics)

(1) Stability: ||u(t, ·)||L∞ 6 ||u0 ||L∞ , a.e. in t ∈ R+

(2) Monotone: if u0 > v0 a.e. in R, then

u(t, ·) > v(t, ·) a.e. in R, a.e. in t ∈ R+

4 Chapter One. Scalar conservation law

(3) TV-diminishing: if u0 ∈ BV(R) then

u(t, ·) ∈ BV(R) and T V(u(t, ·)) 6 T V(u0 )

(4) Conservation: if u0 ∈ L1 (R) then

Z Z
u(t, x)dx = u0 (x)dx, a.e. in t ∈ R+
R R

(5) Finite domain of dependence: if u, v are two entropy solutions corresponding

to u0 , v0 ∈ L∞ and

M = max{|f0 (φ)| : |φ| 6 max(||u0 ||L∞ , ||v0 ||L∞ )}

φ

then Z Z
|u(t, x) − v(t, x)|dx 6 |u0 (x) − v0 (x)|dx
|x|6R |x|6R+Mt
 Chapter Two
DG scheme in 1-D

In finite volume methods, there is one solution variable per cell, the cell average
value, and the solution is assumed to be piecewise constant. To obtain high order
accuracy, we reconstruct the solution by a polynomial inside each cell by making
use of the cell averages in a small stencil around the current cell. IN DG methods
we start with a polynomials solution in each cell which evolved forward in time
by the scheme. DG methods are finite element methods, where we approximate
some function in terms of certain basis functions with compact support. The basis
functions are taken to be polynomials and are allowed to be discontinuous. The
compact support property ensures that the stencil of the scheme is small which
leads to efficient methods for solution. Let us consider a general conservation law of
the form
∂u ∂f
+ =0
∂t ∂x
and try to construct the DG scheme for this problem.

1 Mesh and approximation space

Divide domain Ω = [0, 1] into disjoint cells Ii = [xi− 21 , xi+ 21 ] using the parition
0 = x 12 < x 32 < . . . < xN+ 21 = 1
as shown in figure (Two.1). Define the cell center and cell size as
1
xi = (xi− 21 + xi+ 12 ), ∆xi = xi+ 12 − xi− 21 , h = max ∆xi
2 i

Define the space of broken polynomials

Vhk = {v ∈ L2 (Ω) : v|Ii ∈ Pk (Ii ), 1 6 i 6 N}

x=0 Ii x=1

1 i−1 i i+1 N
1 1 1 1
2
i− 2
i+ 2
N+ 2

Figure Two.1: Mesh for DG scheme

6 Chapter Two. DG scheme in 1-D

Ii−1 Ii Ii−1
1 1
i− 2
i+ 2

Figure Two.2: Discontinuous solution

Note that these functions can be discontinuous on the boundary of the elements as
shown in the figure (Two.2). Define the left and right limits

vh (x− ) = lim vh (x − ), vh (x+ ) = lim vh (x + )

&0 &0

2 Semi-discrete DG scheme
Multiply conservation law by a smooth test function v
Z  
∂u ∂f
+ vdx = 0
Ii ∂t ∂x
and integrate by parts on flux derivative term
Z Z
∂u ∂v
vdx− f(u) dx
Ii ∂t Ii ∂x
+ f(xi+ 12 , t)v(x−
i+ 1
) − f(xi− 21 , t)v(x+
i− 1
)=0
2 2

We now want to replace u with uh and v with vh . At x = xi+ 12 , uh may be

discontinuous, i.e., uh (x− i+ 12
, t) 6= uh (x+
i+ 21
, t). In this case, how to compute the flux
f(xi+ 21 , t) ? Following the finite volume method, we will approximate this flux by a
numerical flux function denoted by f̂i+ 12 (t) = f̂(uh (x− i+ 1
, t), uh (x+
i+ 1
, t)) leading to
2 2

2.1 Definition (Semi-discrete DG scheme) Find uh (·, t) ∈ Vhk such that for all vh ∈ Vhk
Z Z
∂uh ∂vh
vh dx− f(uh ) dx
Ii ∂t Ii ∂x (Two.1)
− +
+ f̂i+ 21 (t)vh (xi+ 1 ) − f̂i− 21 (t)vh (xi− 1 ) = 0
2 2

The initial condition is obtained by an L2 projection onto the finite element

solution space Vhk , i.e.,
Z Z
uh (x, 0)vh (x)dx = u(x, 0)vh (x)dx, ∀ vh ∈ Vhk
Ii Ii
Section 3. Relation to finite volume scheme 7

Note that the numerical flux couples the solution in Ii to those in the neighbouring
elements.

3 Relation to finite volume scheme

If the degree is k = 0, then we have only the constant basis function vh = 1, and
the solution is piecewise constant

uh (x) = ui , x ∈ (xi− 12 , xi+ 12 )

The DG scheme takes the form

dui
∆xi + f̂(ui , ui+1 ) − f̂(ui−1 , ui ) = 0
dt
which is the finite volume scheme. Thus the DG scheme can be considered as a
generalization of finite volume scheme to higher order of accuracy. In FVM only
the cell average values are updated by the scheme and to achieve higher order
accuracy, the solution has to be reconstructed by piecewise polynomials inside each
cell using the cell average values. In DG schemes, we start with a polynomial
solution representation inside each cell and evolve the entire polynomial forward in
time by using the DG scheme.

4 Conservation property
Finite volume schemes are conservative which is necessary to compute correct weak
solutions. We will now show that DG schemes are also conservative. Let degree
k > 0. Take the test function vh ∈ Vhk of the form

1 x ∈ Ii
vh =
0 otherwise

The DG scheme (Two.1) becomes

Z
d
uh dx + f̂i+ 21 (t) − f̂i− 21 (t) = 0
dt Ii

which is a statement of conservation of u in the element Ii .

 Chapter Three
Numerical flux

The numerical flux is an important ingredient of the DG scheme. We can make

use of the developments in finite volume methods to compute the numerical flux
and there are many possible schemes for this. In any case, we desire some common
properties of the numerical flux as follows.

1 Properties of numerical flux

(1) Consistency
f̂(u, u) = f(u)

(2) Locally Lipschitz continuous

|f̂(a2 , b2 ) − f̂(a1 , b1 )| 6 L1 |a2 − a1 | + L2 |b2 − b1 |

(3) monotone

f̂(a, b), increasing in a and decreasing in b, f(↑, ↓)

In terms of the derivatives, we can state this property as

∂ ∂
f̂(a, b) > 0, f̂(a, b) 6 0
∂a ∂b

2 Riemann solver
One way to obtain a numerical flux is to solve a Riemann problem. To compute the
flux f̂i+ 21 , we solve the Riemann problem

∂w ∂f(w) u(x− 1 , t) x < xi+ 12
i+
+ = 0, w(x, τ = t) = 2
, τ>t
∂τ ∂x u(x+ 1 , t) x > xi+ 12
i+ 2

This has a self-similar solution

w(x, τ) = wR ((x − xi+ 12 )/(t − τ); u(x−

i+ 1
, t), u(x+
i+ 1
, t))
2 2
10 Chapter Three. Numerical flux

and the flux across x = xi+ 21 is given by

f̂i+ 12 = f(wR (0; u(x−

i+ 1
, t), u(x+
i+ 1
, t)) (Godunov flux)
2 2

In practice, the Riemann problem is solved approximately.

3 Numerical flux: Linear convection equation

The flux is linear f(u) = au where a = constant and the equation is of the form

∂u ∂u
+a =0
∂t ∂x

The upwind numerical flux is given by


au− if a > 0
f̂(u− , u+ ) =
au+ if a < 0

This can also be written as

1 1
f̂(u− , u+ ) = [f(u− ) + f(u+ )] − |a|(u+ − u− )
|2 {z } |2 {z }
centered flux dissipative flux

The second part of the flux contributes to the numerical dissipation that gives rise
to stability.

4 Numerical flux: Non-linear conservation law

Here we list some standard numerical fluxes and the reader must consult the literature
on finite volume schemes for more details.

4.1 Godunov flux

This is obtained by exactly solving the Riemann problem. The general formula is


 min f(u), if u− 6 u+
−
u∈[u ,u+ ]
f̂(u− , u+ ) =

 max f(u), if u − > u+
+ −
u∈[u ,u ]

For convex flux, we can give simpler expression as follows. Let u∗ be the only sonic
point, i.e., f 0 (u∗ ) = 0. Then the Godunov flux is given by

f̂ = max{f(max{u∗ , u− }), f(min{u∗ , u+ })}

Section 4. Numerical flux: Non-linear conservation law 11

4.2 Lax-Friedrich flux

This may be considered as a generalization of the upwind flux formula by using the
wave speed to be ∆x/∆t

1 1 ∆x +
f̂(u− , u+ ) = [f(u− ) + f(u+ )] − (u − u− )
2 2 ∆t

4.3 Rusanov flux

This flux is also refered to as a local Lax-Friedrich flux and it makes use of a local
wave speed estimate in the dissipative flux
1 1
f̂ = [f(u− ) + f(u+ )] − λ(u+ − u− )
2 2
where
λ= max |f 0 (ξ)|
ξ∈(u− ,u+ )

A simple choice, which is exact for convex fluxes, is

λ = max{|f 0 (u− )|, |f 0 (u+ )|}

4.4 Roe flux

This flux is based on a approximate Riemann solver. The non-linear Riemann

problem is replaced by the linear problem
∂u ∂u
+λ =0
∂t ∂x
and the exact solution of this problem gives the flux
1 1
f̂ = [f(u− ) + f(u+ )] − λ(u+ − u− )
2 2
where
f(u+ ) − f(u− )
λ=
u+ − u−
We have to modify the Roe scheme to satisfy entropy condition.
3.1 Remark Many of the numerical fluxes have the upwind property in the sense
that if f 0 > 0 in the Riemann problem, then f̂(u− , u+ ) = f(u− ). This is not the case
for Lax-Friedrich and Rusanov fluxes, which are hence called central fluxes.
 Chapter Four
Stability

1 Linear convection equation

Let us consider the linear convection equation
∂u ∂u
+a =0 x ∈ (0, 1)
∂t ∂x
The semi-discrete DG scheme is given by
Z Z
∂uh ∂vh
vh dx − auh dx + f̂i+ 21 (t)vh (x−
i+ 12
) − f̂i− 21 (t)vh (x+
i− 12
)=0
Ii ∂t Ii ∂x
Taking vh = uh
Z h i
1 d a
u2h dx− 2 − 2 +
uh (xi+ 1 , t) − uh (xi− 1 , t)
2 dt Ii 2 2 2

+f̂i+ 12 (t)uh (x−

i+ 12
, t) − f̂i− 12 (t)uh (x+
i− 1
, t) = 0
2

Summing over all the cells, we obtain the energy equation

1 d a X ha
N−1 i
kuh k2 + u2h (x 12 , t) − f̂ 21 (t)uh (x 12 , t) + Ju2h Ki+ 21 − f̂i+ 21 (t)Juh Ki+ 12
2 dt 2 i=1
2
a
− u2h (xN+ 12 , t) + f̂N+ 21 (t)uh (xN+ 12 , t) = 0
2
where the jump term is defined as

JvKi+ 21 = v(x+
i+ 1
) − v(x−
i+ 1
)
2 2

For the upwind numerical flux, which is given by

a |a|
f̂i+ 12 = [uh (x−
i+ 1 , t) + uh (x
+
i+ 1 , t)] − Juh Ki+ 12
2 2 2 2
we can write
a 2 |a|
Juh Ki+ 12 − f̂i+ 21 (t)Juh Ki+ 12 = Juh K2i+ 1
2 2 2

so that the energy equation is

1 d 2 a 2
X |a|
N−1
kuh k + uh (x 21 , t) − f̂ 21 (t)uh (x 12 , t) + Juh K2i+ 1
2 dt 2 i=1
2 2
(Four.1)
a
− u2h (xN+ 21 , t) + f̂N+ 12 (t)uh (xN+ 21 , t) = 0
2
14 Chapter Four. Stability

1.1 Periodic boundary conditions

In this case
f̂ 21 = f̂N+ 12 = f̂(uh (xN+ 21 , t), uh (x 12 , t))
and (Four.1) simpifies to

|a| X
N
1 d
kuh k2 + Juh K2i+ 1 = 0
2 dt 2 i=1 2

where
Juh KN+ 12 = uh (xN+ 21 ) − uh (x 12 )
Integrating over time
N ZT
X
2
kuh (T )k + |a| Juh K2i+ 1 dt = kuh (0)k2
2
i=1 0

This immediately implies that

kuh (t)k 6 kuh (0)k

The use of upwind scheme together with discontinuous basis functions leads to L2
stability. If we use the central flux
a
f̂i+ 21 = [uh (x−
i+ 12
, t) + uh (x+
i+ 12
, t)]
2
then we conclude that
1 d
kuh k2 = 0
2 dt
and in this case, the energy is conserved, kuh (T )k = kuh (0)k. The scheme does not
have any dissipation and is neutrally stable. However, the central flux should not
be used in practice since it gives less accurate solutions than the upwind flux and
will be unstable for non-linear problems.

1.2 Dirichlet condition

Assume that a > 0. In this case, we can specify boundary conditions only at
x = x 21 = 0, i.e.,
u(0, t) = g(t)
We take the numerical flux at the boundaries as

f̂ 21 (t) = ag(t), f̂N+ 21 (t) = auh (xN+ 12 , t)

which is an upwind flux. This corresponds to a weak implementation of the boundary

condition since the solution does not satisfy uh (x 12 , t) = g(t). The energy equation
is

|a| X
N−1
1 d a a
kuh k2 + u2h (x 21 , t) − ag(t)uh (x 12 , t) + Juh K2i+ 1 + u2h (xN+ 12 , t) = 0
2 dt 2 2 i=1 2 2
Section 2. Non-linear conservation law 15

The boundary terms can be rearranged as

|a| X
N−1
1 d 2 a 2 a 2 a
kuh k + [g(t) − uh (x 21 , t)] − g (t) + Juh K2i+ 1 + u2h (xN+ 12 , t) = 0
2 dt 2 2 2 i=1 2 2

Integrating in time, we obtain the energy equation

ZT  X
N−1

a |a| a
kuh (T )k2 + [g(t) − uh (x 21 , t)]2 + Juh K2i+ 1 + u2h (xN+ 12 , t) dt
0 2 2 i=1
2 2
Z
2 a T 2
= kuh (T )k + g (t)dt
2 0
The jump terms at the boundary and interior faces lead to dissipation of energy
giving the inequality
Z Z
2 2 a T 2 2 a T 2
kuh (T )k 6 kuh (0)k + g (t)dt 6 ku(0)k + g (t)dt
2 0 2 0
The energy in the numerical solution is bounded by the energy in the exact solution.
If g(t) ≡ 0 then the energy will decrease with time and we have L2 stability.

2 Non-linear conservation law

Weak solutions of non-linear conservation laws can be non-unique. To obtain a
unique weak solution, we need to impose an entropy condition.

2.1 Entropy condition

Let U(u) be a convex entropy function and let F(u) be an associated entropy flux
such that
F 0 (u) = U 0 (u)f 0 (u) (Four.2)
If u is a smooth solution, then it satisfies the equation
∂u ∂u
+ f 0 (u) =0
∂t ∂x
Multiplying throughout by U 0 (u)
∂u ∂u ∂U ∂F
U 0 (u) + U 0 (u) f 0 (u) =0 =⇒ + =0
∂t | {z ∂x} ∂t ∂x
∂f
∂x

we see that the smooth solution satisfies an additional conservation law. For a
discontinuous solution, we will demand that it satisfy the entropy inequality
∂U ∂F
+ 60
∂t ∂x
in the sense of distributions with equality in smooth regions. Then Kruzkov theory
shows that the weak solution is unique.
16 Chapter Four. Stability

4.1 Theorem (Cell entropy inequality) The solution uh of the semi-discrete DG scheme
satisfies Z
d
U(uh )dx + F̂i+ 21 (t) − F̂i− 12 (t) 6 0
dt Ii
for the square entropy U(u) = 21 u2 with some consistent numerical entropy flux
F̂i+ 12 (t) = F̂(uh (x−
i+ 1
, t), uh (x+
i+ 1
, t)).
2 2

Proof: Take vh = uh in the DG scheme

Z Z
∂uh ∂uh
uh dx− f(uh ) dx
Ii ∂t Ii ∂x
+ f̂i+ 12 (t)uh (x−
i+ 1
, t) − f̂i− 12 (t)uh (x+
i− 1
, t) = 0
2 2

Define Zu
F̃(u) = f(s)ds =⇒ F̃ 0 (u) = f(u)
0

Integrating the compatibility condition (Four.2), we get

Zu
F(u) − F(0) = uf(u) − f(s)ds = uf(u) − F̃(u)
0

Ignoring the constant term F(0) we have F(u) = uf(u) − F̃(u). Then
Z Z
∂uh ∂uh
− f(uh ) dx = − F̃ 0 (uh ) dx = −F̃(uh (x−
i+ 12
, t)) + F̃(uh (x+
i− 12
, t))
Ii ∂x Ii ∂x

so that the entropy equation becomes

Z
∂U(uh )
−F̃(uh (x−
i+ 21
, t)) + F̃(uh (x+
i− 12
, t))
Ii ∂t
+ f̂i+ 12 (t)uh (x−
i+ 1
, t) − f̂i− 21 (t)uh (x+
i− 1
, t) = 0
2 2

This can be re-written as

Z
d
U(uh )dx + F̂i+ 21 − F̂i− 12 + Θi− 21 = 0
dt Ii

with consistent numerical entropy flux

F̂i+ 12 = −F̃(uh (x−

i+ 1
)) + f̂i+ 12 uh (x−
i+ 1
)
2 2

and

Θi− 12 = −F̃(uh (x−

i− 1
)) + f̂i− 12 uh (x−
i− 1
) + F̃(uh (x+
i− 1
)) − f̂i− 21 uh (x+
i− 1
)
2 2 2 2

Ignoring the suffixes, the quantity Θ can be written as

Θ = F̃(u+ ) − F̃(u− ) + (u− − u+ )f̂

= (u+ − u− )(F̃ 0 (ξ) − f̂), min{u− , u+ } 6 ξ 6 max{u− , u+ }
= (u+ − u− )(f(ξ) − f̂)
Section 2. Non-linear conservation law 17

Now assume that u+ > ξ > u− . Then since f̂ is a monotone flux

f(ξ) = f̂(ξ, ξ) > f̂(u− , ξ) > f̂(u− , u+ )

and hence Θ > 0. In the case u+ 6 ξ 6 u− we can again show that Θ > 0. Thus
the semi-discrete DG scheme satisfies the entropy condition for any order of
the basis functions k. QED
4.2 Remark To obtain entropy inequality, we can also use the E-flux condition

(u+ − u− )(f(ξ) − f̂(u− , u+ )) > 0, ∀ξ between u− , u+

This condition can be extended to system of conservation laws (Barth).

4.3 Corollary (L2 stability) For periodic or compactly supported boundary conditions,
the semi-discrete DG scheme satisfies
Z
d
u2 dx 6 0
dt Ω h
which implies that
kuh (t)k 6 kuh (0)k 6 ku(0)k

Proof: Adding the cell entropy inequality from all the cells

XN Z XN
d
U(uh )dx + [F̂i+ 12 (t) − F̂i− 12 (t)] 6 0
i=1
dt I i i=1

the internal fluxes cancel one another, leading to

Z
1 d
u2h dx + F̂N+ 12 − F̂ 12 6 0
2 dt Ω

For periodic case, F̂ 12 = F̂N+ 12 , while for compactly supported case F̂ 12 = F̂N+ 12 = 0,
we obtain desired result. QED
 Chapter Five
Error estimates

1 Error estimate: ut + cux = 0, semi-discrete scheme

(B. Cockburn, Lecture notes on Discontinuous Galerkin methods for convection

dominated problem)
 Chapter Six
Basis functions

We have to construct basis functions for Vhk for which there are two approaches:
nodal and modal. The DG solution has the form
X
N
x ∈ Ii : uh (x, t) = uij (t)φij (x), N = dim(Vhk ) = k + 1
j=1

and we refer to the set of values {uij : 1 6 j 6 N} as the degrees of freedom or dof
associated with the i’th cell.

1 Nodal basis functions

A degree k polynomial is determined by N = k + 1 values. In the nodal approach,
we choose N distinct nodes in each cell

xij ∈ Ii , j = 1, 2, . . . , N

These nodes can be used to define the Lagrange polynomials of degree k which have
the interpolation property

φij (xil ) = δjl , 1 6 j, l 6 N

This property implies that

uh (xij ) = uij , 16j6N

so that the dofs in this case as the solution values at the nodes as shown in
figure (Six.1). The location of the N nodal points xij can be

• uniformly distributed inside Ii

– We may encounter Runge phenomenon for high degree polynomials

– The mass matrix is full and ill-conditioned for large degree k

• based on Gauss-Legendre or Gauss-Lobatto integration points

Once the nodal points are chosen, the basis functions can be obtained from Lagrange
interpolation.
22 Chapter Six. Basis functions

Ii−1 Ii Ii−1

Figure Six.1: Nodal solution representation

• k = 0: One dof per element


1 x ∈ Ii
φi1 (x) =
0 otherwise

• k = 1: Two dof per element; xi1 = xi− 12 , xi2 = xi+ 12

 x −x  x−x
 i+ 12 x ∈ Ii  i− 1
2
x ∈ Ii
∆xi ∆xi
φi1 (x) = , φi2 (x) =
0 otherwise 0 otherwise

• k = 2: Three dof per element; xi1 = xi− 21 , xi2 = xi , xi3 = xi+ 12

• In general: We choose N distinct nodes {xi1 , xi2 , . . . , xi,N } ⊂ Ii

(x − xi,1 ) . . . (x − xi,j−1 )(x − xi,j+1 ) . . . (x − xi,N )

φij (x) =
(xi,j − xi,1 ) . . . (xi,j − xi,j−1 )(xi,j − xi,j+1 ) . . . (xi,j − xi,N )

6.1 Remark It is efficient to compute the shape functions on a reference cell. Let us
map cell Ii to [−1, +1] by

x − xi 1−ξ 1+ξ
ξ= 1
, x= xi− 12 + xi+ 12
2
∆xi 2 2

Choose a set of distinct nodes ξ0 , ξ1 , . . . , ξk ∈ [−1, +1]. The j’th basis function is
given by

(ξ − ξ0 ) . . . (ξ − ξj−1 )(ξ − ξj+1 ) . . . (ξ − ξk )

φij (x) = ϕ̂j (ξ) =
(ξj − ξ0 ) . . . (ξj − ξj−1 )(ξj − ξj+1 ) . . . (ξj − ξk )

6.2 Remark If nodes are located at the element boundaries, then we have multiple
dofs at the boundary since the solution is in general discontinuous.
Section 2. Taylor basis functions 23

2 Taylor basis functions

Here we do not use nodal basis functions, but use Taylor series to generate the basis
functions. The cell average value ui is one of the degrees freedom. The other dof
are the gradient, hessian, etc.
Define the moments
Z  s
1 x − xi
mis = dx, s = 1, 2, . . .
s!∆xi Ii ∆xi
Then the solution for different degree are taken as follows.
k = 1: Dof are (ui , si )
x − xi
uh (x, t) = ui (t) + si (t)
∆xi
 
x−xi
1 x ∈ Ii ∆xi
x ∈ Ii
φi1 (x) = , φi2 (x) =
0 otherwise 0 otherwise
k = 2: Dof are (ui , si , qi )
"  2 #
x − xi 1 x − xi
uh (x, t) = ui (t) + si (t) + − mi2 qi (t)
∆xi 2 ∆xi
  
 1 x−xi 2 − m x ∈ Ii
i2
φi3 (x) = 2 ∆xi
0 otherwise
In general, the basis functions are chosen such that
Z Z
φi1 = ∆xi , φij = 0, j = 2, 3, . . .
Ii Ii

• Can be extended to higher degrees using Taylor series

• Heirarchical representation

• Diagonal mass matrix for k 6 2

• Extension to multi-dimensions on arbitrary polygonal elements; orthogonalize

using Gram-Schmidt process

3 Orthogonal polynomials (Modal approach)

The Legendre polynomials are the solution of Legendre’s differential equation
 
d 2 d
(1 − ξ ) Pn (ξ) + n(n + 1)Pn (ξ) = 0, n = 0, 1, 2, . . .
dξ dξ
The first few polynomials are given by
P0 (ξ) = 1 P1 (ξ) = ξ
1 2
P2 (ξ) = 2 (3ξ − 1) P3 (ξ) = 21 (5ξ3 − 3ξ)
P4 (ξ) = 18 (35ξ4 − 30ξ2 + 3) P5 (ξ) = 81 (63ξ5 − 70ξ3 + 15ξ)
24 Chapter Six. Basis functions

These polynomials satisfy the following recurrence relation which is useful for
numerical computation
   
2n − 1 n−1
Pn (ξ) = ξPn−1 (ξ) − Pn−2 (ξ), n = 2, 3, . . .
n n

A very useful property of these polynomials is that they are mutually orthogonal,
i.e., 
Z +1
0 j 6= k
Pj (ξ)Pk (ξ)dξ =
2
−1
2j+1
j=k
Using these polynomials, we can define our basis functions as: j = 0, 1, 2, . . .
  Z 
p x − xi 0 j 6= k
φij (x) = 2j + 1Pj , φij φik dx =
∆xi /2 Ii ∆xi j = k

In [Link], the FE_DGP space makes use of these basis functions. It is convenient
to define the functions on the reference cell [−1, +1] as
p x − xi
φij (x) = ϕ̂j (ξ) = 2j + 1Pj (ξ), ξ= 1
2
∆xi

Many authors use the following definition for the basis functions
  Z 
x − xi 0 j 6= k
φij (x) = Pj , φij φik dx = ∆x
∆xi /2 Ii i
j=k
2j+1
 Chapter Seven
Implementation in 1-D

The solution inside cell Ii is a polynomial of degree k and is of the form

X
N
uh (x, t) = uij (t)φij (x), N=k+1
j=1

The semi-discrete DG scheme is given by

Z Z
∂uh ∂φij
φij dx− f(uh ) dx
Ii ∂t Ii ∂x
+ f̂i+ 12 (t)φij (x−
i+ 1
) − f̂i− 12 (t)φij (x+
i− 1
) = 0, 16j6N
2 2

1 ODE system
Plugging the solution uh in the DG scheme, the first term is of the form
Z XN Z XN
∂uh duik (i) duik
φij dx = φij φik dx = Mjk
Ii ∂t k=1
dt Ii k=1
dt
(i)
The quantities Mjk form the elements of the mass matrix M(i) ∈ RN×N . This is
evaluated using a q-point quadrature rule as
Z X
(i)
Mjk = φij φik = ωiq φij (xiq )φik (xiq ), 1 6 j, k 6 N
Ii q

where {xiq } are the quadrature nodes in element Ii and ωiq are the corresponding
weights. The remaining terms which we will put on the right hand side are computed
as
Z
(i) 0
L (U(t))j = f(uh )φij − f̂i+ 21 (t)φij (x−
i+ 12
) + f̂i− 12 (t)φij (x+
i− 12
)
Ii
X
0
≈ ωiq f(uh (xiq , t))φij (xiq ) − f̂i+ 12 (t)φij (x−i+ 1
)
2
q

+f̂i− 12 (t)φij (x+

i− 1
)
2

Let us put the degrees of freedom associated with Ii in a vector

U(i) = [ui1 , . . . , uiN ]> ∈ RN

26 Chapter Seven. Implementation in 1-D

Then we obtain the following system of ODE for element Ii

dU(i) (t)
M(i) = L(i) (U(t))
dt
The mass matrix is obviously symmetric, and it is also positive definite. For any
vector U = [u1 , . . . , uN ]> ∈ RN
Z X
n
>
U M U= (i)
| uj φij |2 dx > 0, U 6= 0
Ii j=1

Hence we can invert the mass matrix and write the ODE as
dU(i) (t)
= [M(i) ]−1 L(i) (U(t)) = R(i) (U(t))
dt

2 Quadrature rules
The quadrature rule must be chosen such that mass matrix is evaluated exactly.
Moreover, in order to achieve the optimal convergence rate, the integral inside the
cell must be computed with a quadrature rule which is exact for polynomials of
degree 2k. These requirements are satisfied by using

• (k + 1)-point Gauss-Legendre, which is exact for P2k+1

• (k + 2)-point Gauss-Lobatto-Legendre, which is exact for P2k+1

The quadrature rules are explained in section (1). The elements of the mass matrix
can be computed as
Z X
(i)
Mjk = φij φik = ωiq φij (xiq )φik (xiq ), 1 6 j, k 6 N
Ii q

where
1 1 − ξq 1 + ξq
ωiq = ωq ∆xi , xiq = xi− 12 + xi+ 12
2 2 2
and {ξq } are the quadrature points in the reference element [−1, +1]. If we use the
orthogonal Legendre basis, then the mass matrix can be computed directly as
Z
(i)
Mjk = φij φik dx = ∆xi δjk
Ii

In this case, the mass matrix is diagonal, i.e., M(i) = ∆xi IN×N . If we use nodal
Lagrange basis at Gauss-Legendre points and if the quadrature is exact, then the
mass matrix is again diagonal. Using (k + 1)-point Gauss-Legendre quadrature,
which should give the exact mass matrix, we get
Z X
(i)
Mjk = φij φik = ωiq φij (xiq )φik (xiq ), 1 6 j, k 6 N
Ii q
X
= ωiq δjq δkq = ωij δjk
q
Section 3. Setting initial condition 27

The right hand side is also computed using a quadrature rule

Z
(i) 0
L (U(t))j = f(uh )φij dx − f̂i+ 21 (t)φij (x−i+ 21
) + f̂i− 21 (t)φij (x+i− 21
)
Ii
X
0
≈ ωiq f(uh (xiq , t))φij (xiq ) − f̂i+ 21 (t)φij (x− i+ 1
)
2
q

+f̂i− 12 (t)φij (x+

i− 1
)
2

and if f(u) is non-linear, then the quadrature is not exact. The derivative of the
test function is given by
∂ dξ d 1
φij (xiq ) = ϕ̂j (ξq ) = ∆xi ϕ̂j0 (ξq )
∂x dx dξ 2
and hence the quadrature is given by
Z X
0 0 1X
f(uh )φij dx ≈ ωiq f(uh (xiq , t))φij (xiq ) = ωq f(uh (ξq ))ϕ̂j0 (ξq )
Ii q
2 q

The solution at the quadrature points are given by

X
N
uh (ξq ) = uij ϕ̂j (ξq )
j=1

Thus, all of computations involving basis functions can be done on the reference
cell which is common to all the cells.

3 Setting initial condition

Suppose the initial condition is
u(x, 0) = u0 (x)
• Nodal basis: We can interpolate initial condition
uij (0) = uh (xij , 0) = u0 (xij )
Alternately, we can also set the initial condition by performing an L2 projection.

• Modal basis: Do an L2 projection, i.e.,

Z
min (uh (x, 0) − u0 (x))2 dx wrt ui1 , . . . , uiN
Ii
The first order optimality condition is
Z
d
(uh (x, 0) − u0 (x))2 dx = 0, k = 1, 2, . . . , N
duik Ii
Z Z
=⇒ uh φik dx = u0 φik dx, k = 1, 2, . . . , N
Ii Ii
Using the orthogonality of the basis functions, and a quadrature rule for the
right hand side, we get
X
uik ∆xi = u0 (xiq )φik (xiq )ωiq
q

which determines the dofs of the initial solution.

28 Chapter Seven. Implementation in 1-D

4 Boundary condition
We can specify Dirichlet boundary condition at some boundary point if the charac-
teristics are entering the domain (inflow boundary) at that point. For example, at
x = 0 if f 0 > 0, then we can specify the boundary condition on u. In general let us
take the boundary conditions

u(0, t) = a(t), u(1, t) = b(t)

The boundary conditions are incorporated in the DG scheme via the boundary fluxes

f̂ 12 (t) = f̂(a(t), uh (0+ , t)), f̂N+ 21 (t) = f̂(uh (1− , t), b(t))

In case of periodic boundaries, the fluxes at the boundary are computed as

f̂ 12 (t) = f̂(uh (x−

N+ 1
, t), uh (x+1 , t)), f̂N+ 12 (t) = f̂(uh (x−
N+ 1
, t), uh (x+1 , t))
2 2 2 2

5 Strong stability preserving RK schemes

Consider an ODE of the form
dU
= R(U, t)
dt
A standard method to solve ODEs is the Runge-Kutta method. We will consider
the strong stability preserving RK schemes, examples of which are given below.

• 1-stage, first order scheme

Un+1 = Un + ∆t R(Un , tn )

• 2-stage, second order RK

U(1) = Un + ∆t R(Un , tn )
1 n 1 (1)
Un+1 = U + [U + ∆t R(U(1) , tn + ∆t)]
2 2
• 3-stage, third order RK

U(1) = Un + ∆t R(Un , tn )
3 n 1 (1)
U(2) = U + [U + ∆t R(U(1) , tn + ∆t)]
4 4
n+1 1 n 2 (2)
U = U + [U + ∆t R(U(2) , tn + 12 ∆t)]
3 3

7.1 Lemma (SSP Runge-Kutta) If the forward Euler scheme is stable in some norm
and under some time step condition, i.e.,

∆t 6 ∆tc (U) =⇒ kU + ∆t R(U)k 6 kUk

then the SSPRK scheme is stable under a CFL condition ∆t 6 α∆tc . For the second
and third order schemes, α = 1.
Section 6. CFL condition 29

Proof: Let us show this for the 2-stage scheme. Since the forward Euler scheme is
stable if ∆t 6 ∆tc , and as first stage resembles a forward Euler scheme, we get

kU(1) k = kUn + ∆t R(Un )k 6 kUn k

For the second stage, we apply triangle inequality

1 1 1 1
kUn+1 k 6 kUn k + kU(1) + ∆t R(U(1) )k 6 kUn k + kU(1) k 6 kUn k
2 2 2 2
since the second term has the form of a forward Euler scheme. QED
7.2 Remark The allowed time step ∆tc depends on the current solution. The above
proof is correct provided

∆t 6 min{∆tc (Un ), ∆tc (U(1) )}

But we do not know U(1) at the beginning of the time step. In practice, we can put
a margin of safety, e.g., ∆t = α∆tc (Un ) where α = 0.9 for example, which usually
works in practice. If it turns out after the first stage, that the ∆t we have chosen is
larger than ∆tc (U(1) ), then we can reduce it further and restart the time step from
the first stage.
7.3 Remark The above SSPRK scheme require storage for three steps of vectors,
• Un

• current stage solution U(s)

• Residual R
7.4 Remark There exist higher order SSPRK schemes but they need more steps.
There is a 5-stage, 4-th order SSPRK scheme [2].

6 CFL condition
Explicit time integration schemes are stable only under a restriction on the time step
∆t. For k = 0, the scheme is L2 stable under CFL number of one. For k > 1, the
forward Euler scheme (RK1) is known to be unconditionally unstable in L2 (Chavent
and Cockburn, 1989) if the CFL number is of order unity. It is L2 stable for finite



time intervals if the CFL number is O h1/2 which means that ∆t = O h3/2 which
is very restrictive. For DG space discretizations using polynomials of degree k, and
a (k + 1)-stage RK method of order k + 1, a Von-Neumann stability analysis for the
one-dimensional linear case
f(u) = cu
with upwind flux gives the CFL condition
∆t 1
|c| 6
∆x 2k + 1
Theoretical proof of this is available only for k = 0, 1, 2. For k > 3 the above
condition is close to the numerically determined values of CFL numbers.
stable and high-order accurate in time method. We must say also that there
are anecdotal reports of other time discretizations that seem to work just
fine. However, the fact remains that only with this special class of TVD–
30 Chapter Seven. Implementation in 1-D
RK methods can the non-linear stability of the method actually be proven.
Let us The
7.1• Remark alsoCFLstress the fact
number that the
for SSPRK generalized
scheme to be TVDM slope limiterthan
is higher is likewise
that
an required
essentialforingredient
L [Link] However,
2
the method without
to control which
round-off its stability
errors, the smaller CFL be
cannot
2
guaranteed. Although
condition from our has
L stability numerical
to be usedexperience
in practical indicates that second-order
computations.

Table 2.3. CFL Numbers for RKDG Methods of Order k+1

k 0 1 2

CFL TV 1 1/2 1/2

CFL L 2 1 1/3 1/5

7 Algorithm

Let us now summarize the main steps in the DG scheme.

• Compute and store the mass matrix if needed.

• Determine u0h from initial condition u0 by an L2 -projection or interpolation

• Set time counter t = 0

• For n = 0, 1, . . .

– Compute time step ∆t from CFL condition

– Set un,0
h = un
h

– RK stages: For r = 0, 1, . . . , Nrk − 1

∗ Compute right hand side Lh (un,r
h )

∗ Update solution to next RK stage

un,r n,r+1
h → uh

– Increment time counter t = t + ∆t

8 Numerical example

Let us apply the DG scheme to linear and non-linear problems with smooth and
discontinuous initial conditions.

• Linear convection equation: smooth initial condition with periodic boundary

conditions
u(x, 0) = sin(πx), x ∈ [−1, +1]
Section 8. Numerical example 31

• Linear convection equation: continuous initial condition with periodic bound-

ary conditions



1
1 + 2x, − 2 6 x 6 0
u(x, 0) = 1 , x ∈ [−1, +1]
1 − 2x, 06x6


2
0, otherwise

• Linear convection equation: discontinuous initial condition with periodic

boundary conditions

1, |x| < 14
u(x, 0) = , x ∈ [−1, +1]
0, otherwise

• Burgers equation: smooth initial condition

u(x, 0) = sin(2πx), x ∈ [0, 1]

 Chapter Eight
Limiters and TVD property

When the solution is discontinuous or has large gradients, the higher order DG
scheme produces oscillatory solution. This situation is similar to high order finite
volume schemes and is related to loss of TVD property. The oscillatory numerical
solution has more total variation than the initial condition. In the case of finite
volume schemes, this problem is resolved by reducing the slope of the reconstructed
solution by appropriate limiter functions so that the scheme becomes TVD. In DG
schemes we do not have to perform any reconstruction since we have a polynomial
inside each cell. But we can borrow the limiter idea and reduce the slope of the
solution in each cell to achieve TVD property. The approach we will take is to
construct a limiter so that the DG scheme with forward Euler discretization is TVD.
The use of an SSPRK scheme then automatically gives TVD property for higher
order versions of the DG scheme.

1 Limiter for DG scheme

Consider the forward Euler time discretization, i.e., find wn+1
h such that
Z Z
wn+1 − un
h ∂vh
h
vh dx − f(un
h) dx + f̂n v (x− ) − f̂n
i+ 12 h i+ 12
v (x+ ) = 0
i− 21 h i− 12
Ii ∆t Ii ∂x

The solution wn+1

h may be oscillatory. We will treat this as a provisional solution
and limit it to obtain the solution at the next time step

un+1
h = ΛΠh (wn+1
h )

We require that the limiter ΛΠh () satisfy some basic properties as follows.
Properties of ΛΠh (·)

(1) It should not change the cell average value.

(2) It should not affect the accuracy in smooth regions.

Define the cell average value Z

1
ūi = uh dx
∆xi Ii
 34 Chapter Eight. Limiters and TVD property

the forward and backward differences

ûi = uh (x−
i+ 1
) − ūi , ǔi = ūi − uh (x+
i− 1
)
2 2

∆+ ūi = ūi+1 − ūi , ∆− ūi = ūi − ūi−1

For k = 1 we have ûi = ǔi . We cannot modify the cell average value but we can
modify the slopes
(m) (m)
ûi = m(ûi , ∆+ ūi , ∆− ūi ), ǔi = m(ǔi , ∆+ ūi , ∆− ūi )

where m is the minmod function


s min(|a1 |, . . . , |al |) s = sign(a1 ) = . . . = sign(al )
m(a1 , . . . , al ) =
0 otherwise
If all arguments have same sign, the minmod function returns the one with smallest
magnitude, otherwise it returns zero. The trace values are recomputed using the
limited slopes
(m) (m) (m) (m)
uh (x−
i+ 1
) = ūi + ûi , uh (x+
i− 1
) = ūi − ǔi
2 2

(m) (m)
ath. Mech. 83, No. 11For k=
(2003) 1, we have ûi = ǔi
/ [Link] . The effect of the limiter is to reduce
749 the slope of

0.75

0.5

0.25

-0.25
0.1 0.125 0.15 0.075 0.1 0.125 0.15

Example of slope limiters: The MUSCL limiter (left) and the less restrictive ΛΠ1h limiter (right). Displayed are the
ans of uh (thick line),Figure
the linearEight.1:
function uh Effect
in the element of the middle
of applying TVD before limitingis
limiter (dotted line), andthe
to reduce the resulting
slope of the solution
after limiting (solid line).
so that it is less oscillatory

olutions that are not necessarily piecewise-constant, the above result still holds provided that the following
atisfied: the solution in case it is larger than the finite difference slope of the cell average
"values as# shown in $ figure (Eight.1). If the cell averages are monotone around cell i,
u+ +
j+1/2 − uj−1/2 = sign uj+1 − uj ,
"
then the limited linear polynomial satisfies
# $
u− −
j+1/2 − uj−1/2 = sign uj − uj−1 . (m)
min{ūi−1 , ūi+1 } 6 uh (x) 6 max{ūi−1 , ūi+1 }, x ∈ Ii
ns are not necessarily satisfied it is necessary to enforce them by means of what will be called a generalized
For k = 0, 1, 2, this procedure uniquely determines a new polynomial of degree
e what happened in k.
the previous
For k =section,
1, letthe
usintroduction
denote the of the numerical
limited traces is not
function by enough tohguarantee
ΛΠ1h (u ). For k > 3 there is
DG method. For non-linear hyperbolic problems, the use of a generalized slope limiter is indispensable,
more freedom since the cell average and the two trace values do not completely
own for the so-called high-resolution methods. See also [27] for a motivation of the introduction of this
determine the polynomial. One approach is to determine the remaining dofs by an
L2 projection. A more simple approach that works well in practice is the following.
he generalized slope limiter
functions

v j + (x − xj ) vx,j ,
Section 1. Limiter for DG scheme 35

(1) If the limiter does not modify the trace values, i.e.,
(m) (m)
uh (x+
i− 1
) = uh (x+
i− 1
) and uh (x−
i+ 1
) = uh (x−
i+ 1
)
2 2 2 2

(m)
then take uh = uh for x ∈ Ii .
(m)
(2) Otherwise, let u1h ∈ P1 (Ii ) be the L2 projection of uh |Ii . Take uh |Ii =
ΛΠ1h (u1h ).
If the DG solution has been processed by the above limiter, then we can prove a
TVD property. Let us first recall Harten’s sufficient conditions for a finite volume
scheme to be total variation diminishing.

8.1 Lemma (Harten) If a scheme can be written as

ūn+1
i = ūn n n n n
i + Ci+ 1 ∆+ ūi − Di− 1 ∆− ūi
2 2

Assume that boundary conditions are periodic or compactly supported. If

Ci+ 12 > 0, Di+ 21 > 0, Ci+ 21 + Di+ 21 6 1

then the scheme is TVD

TV(un+1
h ) 6 TV(un
h)

where the total variation is defined as

X
TV(uh ) = |∆+ ūi |
i

Proof: Take the update equation at i and i + 1

ūn+1
i+1 = ūn
i+1 + C
n n n n
3 ∆+ ūi+1 − Di+ 1 ∆− ūi+1
i+ 2
2
ūn+1
i = ūn
i + Cn ∆
i+ 1 + i
ūn − Dn ∆ ūn
i− 1 − i
2 2

and subtracting them yields

∆+ ūn+1
i = ∆+ ūn
i +C
n n
3 ∆+ ūi+1 − Dn ∆ ūn − Cn
i+ 1 − i+1
∆ ūn + Dn
i+ 1 + i
∆ ūn
i− 1 − i
i+ 2 2 2
2
= ∆+ ūn
i + C 3 ∆+ ūn
n
i+1 − Dn ∆ ūn
i+ 21 + i
− Cn ∆ ūn
i+ 21 + i
+ Dn ∆ ūn
i− 12 + i−1
i+
2
= (1 − Ci+ 12 − Di+ 12 )∆+ ūn
i +C
n n
3 ∆+ ūi+1 + Dn ∆ ūn
i− 1 + i−1
i+ 2
2

Applying triangle inequality and noting that all the coefficients are positive according
to our assumption

|∆+ ūn+1
i | 6 (1 − Ci+ 21 − Di+ 21 )|∆+ ūn
i |+C
n
3 |∆+ ūi+1 |
n
+ Dn
i− 1
|∆+ ūn
i−1 |
i+ 2
2

Summing over all cells, all the terms on the right cancel except the first one, yielding
X X
|∆+ ūn+1
i |6 |∆+ ūn
i |
i i

which proves the lemma. QED

36 Chapter Eight. Limiters and TVD property

8.2 Remark Consider the first order upwind finite volume scheme for ut + aux = 0
which is identical to the DG scheme for degree k = 0. This scheme can be written
in the incremental form with
∆t |a| − a ∆t |a| + a
Ci+ 12 = , Di+ 21 =
∆x 2 ∆x 2
The coefficients and the condition
|a|∆t
Ci+ 12 + Di+ 21 = 61
∆x
leads to the CFL condition.

8.3 Definition (TVM) For DG solution uh , let us define the total variation of the means
by X
TVM(uh ) = |∆+ ūi |
i

We will show that the TVM does not increase with time if the limiter is applied.
For this we make use of the monotone flux property. Since f̂(u, v) is increasing
function of u and decreasing function of v, we have

f̂(a, b) − f̂(c, b) f̂(a, b) − f̂(a, c)

> 0, 60
a−c b−c
We can bound the above ratios as follows. Since
f̂(a, b) − f̂(c, b) ∂
= f̂(ξ, b), ξ ∈ I(a, c)
a−c ∂u
we have
f̂(a, b) − f̂(c, b) ∂
6 max f̂(ξ, b) =: L1 (b)
a−c ξ∈I(b,c) ∂u

Similarly
f̂(a, b) − f̂(a, c) ∂
− = − f̂(a, η), η ∈ I(b, c)
b−c ∂v
so that
f̂(a, b) − f̂(a, c) ∂
− 6 max f̂(a, η) =: L2 (a)
b−c η∈I(b,c) ∂v

8.4 Theorem For periodic or compactly supported boundary conditions, the DG

scheme with the minmod limiter is TVD in the means, i.e.,

TVM(un+1
h ) 6 TVM(un
h)

Proof: Taking vh = 1 for x ∈ Ii

∆t
wi = ūi − λi [f̂i+ 12 − f̂i− 12 ], λi =
∆xi
= ūi − λi [f̂(ūi + ûi , ūi+1 − ǔi+1 ) − f̂(ūi−1 + ûi−1 , ūi − ǔi )]
Section 1. Limiter for DG scheme 37

We can write this in incremental form with

f̂(ūi + ûi , ūi+1 − ǔi+1 ) − f̂(ūi + ûi , ūi − ǔi )
Ci+ 12 = −λi
∆+ ūi

f̂(ūi + ûi , ūi − ǔi ) − f̂(ūi−1 + ûi−1 , ūi − ǔi )

Di− 12 = λi
∆− ūi
Rewrite the coefficient
 
ǔi+1 ǔi
Ci+ 12 = −λi f̂2 1 − +
| {z } ∆+ ūi ∆+ ūi
>0

where
f̂(ūi + ûi , ūi+1 − ǔi+1 ) − f̂(ūi + ûi , ūi − ǔi )
0 6 −f̂2 = − 6 L2
(ūi+1 − ǔi+1 ) − (ūi − ǔi )
Since un
h has been pre-processed by the minmod limiter, we have

ǔi+1 ǔi ǔi+1 ǔi

06 6 1, 06 61 =⇒ 1− + 62
∆+ ūi ∆+ ūi ∆+ ūi ∆+ ūi
and hence
0 6 Ci+ 12 6 2λi L2
and similarly
0 6 Di+ 12 6 2λi+1 L1
If the time step satisfies the condition
1
Ci+ 21 + Di+ 12 = 2(λi L2 + λi+1 L1 ) 6 1 or ∆t 6
2( ∆xLi+1
1
+ L2
∆xi
)

then from Hartens’s Lemma, we have

TVM(wn+1
h ) 6 TVM(un
h)

We now apply the limiter to obtain un+1

h = ΛΠh (wn+1
h ), and we know that un+1
h
and wn+1
h have the same cell average values. Hence

TVM(un+1
h ) = TVM(wn+1
h ) 6 TVM(un
h)

which proves the theorem. QED

8.5 Remark In the case of linear convection equation and upwind flux, we have
L1 = 12 (|a| + a), L2 = 12 (|a| − a) so that the time step on a uniform mesh should
satisfy the condition |a|∆t
∆x
6 12 .
8.6 Remark On non-uniform grids, we can define the limiter as
!
(m) ûi ∆+ ūi ∆− ūi
ûi = ∆xi m ,1 ,1
∆xi 2 (∆xi + ∆xi+1 ) 2 (∆xi + ∆xi−1 )
!
(m) ǔi ∆+ ūi ∆− ūi
ǔi = ∆xi m ,1 ,1
∆xi 2 (∆xi + ∆xi+1 ) 2 (∆xi + ∆xi−1 )
38 Chapter Eight. Limiters and TVD property

The factor in Ci+ 12 is

ǔi+1 ǔi ∆xi + ∆xi+1
1− + 61+
∆+ ūi ∆+ ūi 2∆xi
and in Di+ 21 is
ûi ûi+1 ∆xi + ∆xi+1
1− + 61+
∆− ūi+1 ∆− ūi+1 2∆xi+1
The CFL condition is
∆xi + ∆xi+1 ∆xi + ∆xi+1
Ci+ 12 + Di+ 12 6 (1 + )λi L1 + (1 + )λi+1 L2 6 1
2∆xi 2∆xi+1

2 Limiters: Implementation
Degree k = 1: Let us write the solution in terms of Taylor or Legendre basis
!
x − xi
uh = ūi + 1 si
2
∆xi
Note that uh (x−
i+ 21
) = ūi + si and uh (x+
i− 12
) = ūi − si . We limit the slope with the
minmod function
(m)
si = minmod (si , ūi − ūi−1 , ūi+1 − ūi )
The limited solution is
!
x − xi (m)
ΛΠ1h (uh ) = ūi + 1
si
2
∆xi
Degree k > 1: Let us write the solution in terms of Taylor or Legendre basis
!
x − xi
uh = ūi + 1 si + HOT
2
∆xi
We obtain a limited slope
(m)
si = minmod (si , ūi − ūi−1 , ūi+1 − ūi )
(m)
If si = si then
ΛΠh (uh ) = uh ∈ Pk
else !
x − xi (m)
ΛΠh (uh ) = ūi + 1
si ∈ P1
2
∆xi

8.1 Remark If we use the Legendre polynomials to form the basis functions as
√
ϕ̂j (ξ) = 2j + 1Pj (ξ), then the linear solution is of the form
!
√ x − xi
uh = ūi + 3 1 si
2
∆x i

The limited slope is given by

1 √ 
(m)
si = √ minmod 3si , ūi − ūi−1 , ūi+1 − ūi
3
Section 3. Numerical example 39

3 Numerical example
Discontinuous solution and smooth solution
Numerical example
If the solution in cell Ii has an extremum, see figure (xxx), then the slope in the
cell has a different sign compared to ∆− ūi , ∆+ ūi and the limiter returns zero slope.
The limited solution becomes constant and equal to the cell average value. Thus
there is loss of accuracy at smooth extrema due to TVD limiter. We have to modify
the limiter so that it does not change the solution at smooth extrema.

4 TVB Limiter
In smooth regions of the solution, the differences inside the cell are
1
1

ûi = ux (xi )∆xi + O h2 , ǔi = ux (xi )∆xi + O h2
2 2
while the differences of cell averages are
1

∆+ ūi = ux (xi )(∆xi + ∆xi+1 ) + O h2
2
1

∆− ūi = ux (xi )(∆xi + ∆xi−1 ) + O h2
2

If the solution is smooth and monotone around Ii , and since the O h2 term is
small, all the above quantities have the same sign, so that the limiter yields
(m) (m)
ûi = ûi , ǔi = ǔi
Since the solution is not modified by the limiter, we obtain the full accuracy of the
scheme. However, if there is a smooth extremum in cell Ii , then ux (xi ) ≈ 0 and the
sign is determined by the second derivative. If M is the magnitude of the second


derivative at smooth extrema, then ûi = O Mh2 . This motivates the definition of
the TVB limiter function as

a1 if |a1 | 6 Mh2
m̃(a1 , a2 , . . . , al ) =
m(a1 , a2 , . . . , al ) otherwise
If we are near a smooth extremum, the TVB limiter returns the original slope and
the solution is not modified in that cell. With the above limiter, the scheme is no
longer TVDM and the TVM can increase. However this violation is small and of
the other of the mesh size.

8.1 Lemma (TVB property) With the TVB limiter, if the CFL condition
1
∆t 6 , ∀i
2( ∆xLi+1
1
+ L2
∆xi
)
is satisfied, then
TVM(un+1
h ) 6 TVM(un
h ) + CMh
40 Chapter Eight. Limiters and TVD property

Proof: See [1]

8.2 Remark If we are interested in the solution in a finite time interval [0, T ], then
the TVB limiter ensures that

TVM(un
h ) 6 TV(u0 ) + CMT, ∀ n∆t 6 T

so that the numerical solutions have bounded variation. This is sufficient to prove
convergence to a weak solution.
8.3 Remark The quantity M is an estimate of the second derivative of the solution
at smooth extrema. This can be based on the initial condition, e.g.,

M = max{|u000 (x)| : u00 (x) = 0}

x

Ideally M should be estimated from the numerical solution but there is no reliable
way to do this. The solution may have several extrema with different magnitude of
the second derivatives. In practice people choose the value of M by some trial and
error. But this is still a weak point of the TVB limiter.

Numerical example

5 Algorithm
The DG scheme together with the TVD/TVB limiter is as follows.

• Compute and store the mass matrix

• Determine w0h from initial condition u0 by an L2 -projection

• Find u0h from w0h by applying the limiter, u0h = ΛΠh (w0h )

• For n = 0, 1, . . .

– Compute time step from CFL condition

– Set un,0
h = un
h

– RK stages: For r = 0, 1, . . . , Nrk − 1

∗ Compute right hand side Lh (un,r
h )

∗ Update solution to next RK stage

un,r n,r+1
h → wh

∗ Apply limiter
un,r+1
h = ΛΠh (wn,r+1
h )
 Appendix A
Quadrature rules

1 Quadrature in 1-D
The DG scheme involves integrals which must be approximated by quadrature. Let

f : [−1, +1] → R

Choose n quadrature nodes {ξ1 , ξ2 , . . . , ξn } ⊂ [−1, +1].

Z +1 X
n
f(ξ)dξ ≈ ωq f(ξq )
−1 q=1

• Gauss-Legendre quadrature

– Nodes {ξq } are roots of Legendre polynomial Pn (ξ)

– n-point rule is exact for any f ∈ P2n−1

• Gauss-Lobatto-Legendre quadrature

– Nodes include {−1, +1} and the roots of Pn−1

0
(ξ)
– n-point rule is exact for any f ∈ P2n−3

For function on general interval f : [a, b] → R, do change of variable

1−ξ 1+ξ
x(ξ) = a+ b, ξ ∈ [−1, +1]
2 2
and
Zb Z +1
1
f(x)dx = (b − a) f(x(ξ))dξ
a 2 −1

1 X
n
≈ (b − a) ωq f(x(ξq ))
2 q=1
X
n
1
= ω̃q f(x(ξq )), ω̃q = (b − a)ωq
q=1
2
Bibliography

[1] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discon-
tinuous Galerkin finite element method for conservation laws II: General
framework, Mathematics of Computation, 52 (1989), pp. pp. 411–435.

[2] R. Spiteri and S. Ruuth, A new class of optimal high-order strong-stability-

preserving time discretization methods, SIAM Journal on Numerical Analysis,
40 (2002), pp. 469–491.