Numerical Hyperbolic Systems

Download as pdf or txt
Download as pdf or txt
You are on page 1of 42
At a glance
Powered by AI
The document discusses several numerical methods for solving hyperbolic systems of conservation laws arising in physics, including finite difference methods, finite volume methods, discontinuous Galerkin methods, and the Godunov method.

Finite difference methods, finite volume methods, discontinuous Galerkin methods, and the method of lines are discussed as approaches for solving hyperbolic systems numerically.

Hyperbolic systems have real eigenvalues and are diagonalizable. Characteristics and weak solutions such as entropy solutions are important concepts. The Riemann problem is also discussed.

Numerical methods for hyperbolic systems

Eric Sonnendr
ucker
Max-Planck-Institut f
ur Plasmaphysik
und
Zentrum Mathematik, TU M
unchen

Lecture notes
Sommersemester 2013
July 16, 2013

Contents
1 Introduction

2 1D linear advection
2.1 Finite Difference schemes for the advection equation
2.1.1 Obtaining a Finite Difference scheme . . . . .
2.1.2 The first order explicit upwind scheme . . . .
2.1.3 The first order upwind implicit scheme . . . .
2.1.4 The method of lines . . . . . . . . . . . . . .
2.1.5 Convergence of finite difference schemes . . .
2.1.6 High-order time schemes . . . . . . . . . . . .
2.1.7 The Cauchy-Kowalevsky procedure . . . . . .
2.2 The Finite Volume method . . . . . . . . . . . . . .
2.2.1 The first order Finite Volume schemes . . . .
2.2.2 Higher order schemes . . . . . . . . . . . . .
2.3 The Finite Element method . . . . . . . . . . . . . .
2.4 The Discontinuous Galerkin (DG) method . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

4
4
4
4
5
6
6
12
15
15
15
17
18
22

3 Linear systems
25
3.1 The Riemann problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 The discontinuous Galerkin method . . . . . . . . . . . . . . . . . . . . . 27
4 Non linear conservation laws
4.1 Characteristics . . . . . . . . . . . . . .
4.2 Weak solutions . . . . . . . . . . . . . .
4.2.1 Definition . . . . . . . . . . . . .
4.2.2 The Rankine-Hugoniot condition
4.2.3 Entropy solution . . . . . . . . .
4.3 The Riemann problem . . . . . . . . . .
4.4 Numerical methods . . . . . . . . . . . .
4.4.1 The Godunov method . . . . . .
4.4.2 Approximate Riemann solvers . .
4.4.3 Higher order methods . . . . . .
4.4.4 Strong stability preserving (SSP)
4.5 Nonlinear systems of conservation laws .
4.5.1 The Rusanov flux . . . . . . . . .
4.5.2 The Roe flux . . . . . . . . . . .

. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
Runge-Kutta methods. .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

29
29
30
30
30
33
34
35
35
37
38
38
39
39
39

Chapter 1

Introduction
Hyperbolic systems arise naturally from the conservation laws of physics. Writing down
the conservation of mass, momentum and energy yields a system of equations that needs
to be solved in order to describe the evolution of the system. In this lecture we will
introduce the classical methods for numerically solving such systems. Up to a few years
ago these were essentially finite difference methods and finite volume methods. But
in the last decades a new class of very efficient and flexible method has emerged, the
Discontinuous Galerkin method, which shares some features both with Finite Volumes
and Finite Elements.
In 1D the systems of conservation laws we consider have the form
u F(u)
+
= 0,
t
x
where u is a vector of unknown values. This can be written also
u
u
+ A(u)
= 0,
t
x
Fi
))i,j . The system is called
where A(u) is the Jacobian matrix with components (( u
j
hyperbolic if for all u the matrix A has only real eigenvalues and is diagonalisable. It is
called strictly hyperbolic if all eigenvalues are distinct.
Examples:

1D Maxwells equation
1D Euler equations

u

2
u +
u + p = 0,
t
x
E
Eu + pu
where , u and E are the density, velocity and energy density of the gas and p is
the pressure which is a known function of .
The idea behind all numerical methods for hyperbolic systems is to use the fact that
the system is locally diagonalisable and thus can be reduced to a set of scalar equations.
For this reason, before going to systems it will be useful to first understand the scalar
case and then see how it can be extended to systems by local diagonalization.

The first part of the lecture will be devoted to the linear case, starting with the
scalar case which boils down to linear advection for which the core methods will be first
introduced. This is fairly straightforward.
Many additional problems arise in the nonlinear case. Indeed in this case even starting from a smooth initial condition discontinuities can appear in the solution. In this
case the concept of weak solutions need to be introduced and there can be several solutions only one of which is physical. We thus need a criterion to find the physical solution
and numerical schemes that capture the physical solution. The notion of conservativity
plays an essential role there. These will be addressed in the second part.

Chapter 2

1D linear advection
2.1

Finite Difference schemes for the advection equation

We consider first the linear 1D advection equation


u
u
+a
=0
t
x

pour x [a, b], t 0.

(2.1)

Let us assume for simplicity that the boundary conditions are periodic. This means that
u and all its derivatives are periodic of period b a. We have in particular u(a) = u(b).
The constant a is given. As the problem is time dependent we also need an initial
condition u(x, 0) = u0 (x).

2.1.1

Obtaining a Finite Difference scheme

We first consider a uniform mesh of the 1D computational domain, i.e. of the interval
[a, b] where we want to compute the solution, see Figure 2.1. The cell size or space step

x0

x1

x2

xN 1 xN

Figure 2.1: Uniform mesh of [a, b]


is defined by x = ba
N where N is the number of cells in the mesh. The coordinates of
the grid points are then defined by xi = a + ix. We then need a time step t and we
will compute approximations of the solution at discrete times tn = nt, n N. As we
assume the solution to be periodic of period b a it will be defined by its values at xi
for 0 i N 1 and we shall have u(xN , tn ) = u(x0 , tn ).
We shall denote by unj = u(xj , tn ).

2.1.2

The first order explicit upwind scheme

A Finite Difference scheme is classically obtained by approximating the derivatives appearing in the partial differential equation by a Taylor expansion up to some given order
which will give the order of the scheme. As we know only the values of the unknown
4

function at the grid points, we use Taylor expansion at different grid points and linearly
combine them so as to eliminate all derivatives up to the needed order.
The same can be done for the time discretisation. For an approximation of order 1
in space and time, we can simply write
u(xj , tn+1 ) u(xj , tn )
u
(xj , tn ) =
+ O(t),
t
t

(2.2)

u(xj , tn ) u(xj1 , tn )
u
(xj , tn ) =
+ O(x).
(2.3)
x
x
Denoting by unj , the approximation of the solution at point xj and time tn and using
the above formulas for the approximation of the partial derivatives we get the following
approximation (2.1) at point xj and time tn :
un+1
unj
j
t

+a

unj unj1
= 0.
x

(2.4)

We thus obtain the following explicit formula which enables to compute un+1
in function
j
of the values of u at time tn and points xj1 , xj and xj1 :
un+1
= unj a
j

t n
(u unj1 ).
x j

(2.5)

Denote by U n the vector of RN whose components are un1 , . . . , unN 1 and

A=

at
x )
at
x

(1

0
..
.
..
.

at
x

..

..

at
x

0
at
(1 x )

The terms at the end of the first line comes from the periodic boundary conditions. We
use that un1 = unN 1 and unN = un0 . Except on the two diagonals all the terms vanish
In this case the scheme (2.5) can be written in matrix form
U n+1 = AU n .

2.1.3

The first order upwind implicit scheme

When using an uncentered difference scheme in the other direction for the time derivative,
we get
u(xj , tn ) u(xj , tn1 )
u
(xj , tn ) =
+ O(t),
(2.6)
t
t
We use the same finite difference approximation for the space derivative. We then get
the following formula
t n
(u unj1 ) = ujn1 .
(2.7)
unj + a
x j
In this case the unj are defined implicitly from the ujn1 as solutions of a linear system.
This is why this scheme is called implicit.

Denote by B the matrix of the linear system:

(1 + at
0
x )

.
..
..
at
.
x

B=
.
.
..
..

at
0
x

at
x

0
(1 + at
)
x

The term at the end of the first line comes from the periodic boundary conditions. We
use that un1 = unN 1 and unN = un0 . The terms not on the two diagonals vanish.
Going now from time step n to n + 1 the implicit scheme in matrix form becomes
BU n+1 = U n .

2.1.4

The method of lines

As we saw, the time discretisation can be performed by finite differences as for the space
discretisation. However it is generally more convenient to separate the space and time
discretisation for a better understanding. The method of lines consists in applying only a
discretisation scheme in space first (this can be Finite Differences or any other scheme).
Then one obtains a system of ordinary differential equations of the form
dU
= AU,
dt
where U (t) is the vector whose components are ui (t) the unknown values at the grid
point at any time t. Then one can use any Ordinary Differential Equation (ODE) solver
for the time discretisation. For example using an explicit Euler method with the upwind
method in space yields the previous explicit upwind scheme and when we use an implicit
Euler method we get the implicit upwind scheme.

2.1.5

Convergence of finite difference schemes

So as to include explicit and implicit schemes, we consider a linear scheme in the following
generic matrix form
LU n+1 = M U n ,
(2.8)
where the matrices L and M , are normalised such that the coefficient of unj , is 1.
Let us denote by

u(x0 , tn )

..
V (tn ) =

.
u(xN 1 , tn ),
where u is the exact solution of the Partial Differential Equation (PDE) that is approximated by the scheme (2.8).
Definition 1 The scheme (2.8) is called consistent of order (q, p) if
kLV (tn+1 ) M V (tn )k = t O(tq + xp ).
Proposition 1 If the exact solution is of class C 3 , the schemes (2.5) and (2.7) are
consistent of order 1 in t and in x.
6

Proof. Lets start with (2.5). Using the Taylor formulas (2.2) and (2.3), we obtain for
the ith line of V (tn+1 ) AV (tn ),
at
(V (tn+1 ) AV (tn ))i = u(xi , tn+1 ) u(xi , tn ) +
(u(xi , tn ) u(xi1 , tn ))
x


u
u
(xi , tn ) + O(t) + a (xi , tn ) + O(x) .
= t
t
x
The result follows as u is a solution of our equation. We thus get the (1,1) consistency
of this scheme.
Then for (2.7), we use (2.6) and (2.3). The ith line of BV (tn ) V (tn1 ),
at
(BV (tn ) V (tn1 ))i = u(xi , tn ) u(xi , tn1 ) +
(u(xi , tn ) u(xi1 , tn ))
x


u
u
= t
(xi , tn ) + O(t) + a (xi , tn ) + O(x) .
t
x
The result follows as u is a solution of the equation. We thus get the (1,1) consistency
of this scheme.
Definition 2 The scheme (2.8) is called stable for some given norm k.k if there exist
constants K and independent of t such that
kU n k KkU 0 k

t such that 0 < t < .

Proposition 2 The scheme (2.5) is stable for the L norm provided


at
1.
x
This condition is called the Courant-Friedrichs-Lewy or CFL condition.
Proof.
Consider the scheme (2.5). Grouping the terms in unj1 , unj and unj+1 , this
scheme (2.5) becomes
at n
at n
un+1
=
uj1 + (1
)u .
j
x
x j
n
n
As, a, t and x are positive, if at
x 1 the two factors of uj1 , uj are positive and
equal to their absolute value. Hence

|un+1
|
j

at
at n
at n
at 
|uj1 | + (1
)|uj |
+ (1
) max |unj |,
j
x
x
x
x

and so
max |un+1
| max |unj |,
j
j

from which it follows that

maxj |unj |

maxj |u0j |

for all n, which yields the L stability.

von Neumann stability analysis. Due to the fact that the discrete Fourier transform conserves the L2 norm because of the discrete Plancherel inequality and that it
diagonalises the Finite Difference operators (provided the original PDE has constant coefficients), it is particularly well adapted for studying the L2 stability. The von Neumann
analysis consists in applying the discrete Fourier transform to the discretised equation.
To make this more precise we first recall the definition and main properties of the discrete
Fourier transform.
Let P be the symmetric matrix formed with the powers of the nth roots of unity
2ijk
2i
the coefficients of which are given by Pjk = 1n e n . Denoting by n = e n , we have
Pjk = 1n njk .
Notice that the columns of P , denoted by Pi , 0 i n 1 are the vectors Xi
normalised so that Pi Pj = i,j . On the other hand the vector Xk corresponds to a
discretisation of the function x 7 e2ikx at the grid points xj = j/n of the interval
[0, 1]. So the expression of a periodic function in the base of the vectors Xk is thus
naturally associated to the Fourier series of a periodic function.
Definition 3 Discrete Fourier Transform.
The dicrete Fourier transform of a vector x Cn is the vector y = P x.
The inverse discrete Fourier transform of a vector y Cn is the vector
x = P 1 x = P x.
Lemma 1 The matrix P is unitary and symmetric, i.e. P 1 = P = P .
Proof. We clearly have P T = P , so P = P . There remains to prove that P P = I.
But we have
n1

n1

l=0

l=0

2i

1 X jl lk
1 X 2i l(jk)
1 1 e n n(jk)
(P P )jk =

=
en
=
,
(jk)
n
n
n 1 e 2i
n
and so (P P )jk = 0 if j 6= k and (P P )jk = 1 if j = k.
= P G, their discrete
Corollary 1 Let F, G Cn and denote by F = P F and G
Fourier transforms. Then we have
the discrete Parseval identity:

= F T G

(F, G) = F T G
= (F , G),

(2.9)

The discrete Plancherel identity:


kF k = kF k,

(2.10)

where (., .) and k.k denote the usual euclidian dot product and norm in Cn .
Proof. The dot product in Cn of F = (f1 , . . . , gn )T and G = (g1 , . . . , gn )T is defined by
(F, G) =

N
X

fi gi = F T G.

i=1

Then using the definition of the inverse discrete Fourier transform, we have F = P F ,
we get
G = P G,

= (P F )T P G
= F T P T P G
FTG
= F T G,
as P T = P and P = P 1 . The Plancherel identity follows from the Parseval identity by
taking G = F .
Remark 1 The discrete Fourier transform is defined as a matrix-vector multiplication.
Its computation hence requires a priori n2 multiplications and additions. But because of
the specific structure of the matrix there exists a very fast algorithm, called Fast Fourier
Transform (FFT) for performing it in O(n log2 n) operations. This makes it particularly
interesting for many applications, and many fast PDE solvers make use of it.
Let us now consider the generic matrix form of the Finite Difference scheme introduced above:
LU n+1 = M U n .
Note that on a uniform grid if the PDE coefficients do not explicitly depend on x the
scheme is identical at all the grid points. This implies that L and M have the same
coefficients on any diagonal including the periodicity. Such matrices, which are of the
form

c0
c1 c2 . . . cn1
cn1 c0 c1
cn2

cn2 cn1 c0
cn3
C=

..
..
..
.
.
.
c1

c2

c3 . . .

c0

with c0 , c1 , . . . , cn1 R are called circulant.


Proposition 3 The eigenvalues of the circulant matrix M are given by
k =

n1
X

cj jk ,

(2.11)

j=0

where = e2i/n .
Proof. Let J be the circulant matrix obtained from M by taking c1 = 1 and cj = 0 for
j 6= 1. We notice that M can be written as a polynomial in J
M=

n1
X

cj J j .

j=0

As J n = I, the eigenvalues of J are the n-th roots of unity that are given by k = e2ik/n .
Looking for Xk such that JXk = k Xk we find that an eigenvector associated to the
eigenvalue k is

1
k

2k
Xk =
.
..
.
(n1)k
9

We then have that

n1
X

M Xk =

cj J Xk =

j=0

n1
X

cj jk Xk ,

j=0

and so the eigenvalues of M associated to the eigenvectors Xk are


k =

n1
X

cj jk .

j=0

Proposition 4 Any circulant matrix C can be written in the form C = P P where


P is the matrix of the discrete Fourier transform and is the diagonal matrix of the
eigenvalues of C. In particular all circulant matrices have the same eigenvectors (which
are the columns of P ), and any matrix of the form P P is circulant.
Corollary 2 We have the following properties:
The product of two circulant matrix is circulant matrix.
A circulant matrix whose eigenvalues are all non vanishing is invertible and its
inverse is circulant.
Proof. The key point is that all circulant matrices can be diagonalized in the same
basis of eigenvectors. If C1 and C2 are two circulant matrices, we have C1 = P 1 P and
C2 = P 2 P so C1 C2 = P 1 2 P .
If all eigenvalues of C = P P are non vanishing, 1 is well defined and
P P P 1 P = I.
So the inverse of C is the circulant matrix P 1 P .
Now applying the discrete Fourier transform to our generic scheme yields:
P LU n+1 = P M U n

P LP P U n+1 = P M P P U n

which is equivalent to
n+1 = M U
n,
L U
where L and M are the diagonal matrices containing the eigenvalues of the circulant matrices M and L which are given explicitly from the matrix coefficients. It fol k = kU k for any vector U that the scheme is L2 stable if and only if
lows because kU
|M,i |
maxi |L,i | 1, where M,i and L,i are the eigenvalues of M and L.
Let us now apply this technique to the scheme (2.7):
Proposition 5 The scheme (2.7) is stable for the L2 norm for all strictly positive values
of x and t.
t
Proof. Let us denote by = a x
. We notice that matrix B is circulant with c0 = 1+,
cn1 = , the other ci being 0.

10

The eigenvalues of B thus are k = c0 + cn1 2ik/n . Which implies that


2k
) 1,
n
as 0. It follows that all eigenvalues of B have a modulus larger or equal to 1, which
implies that B is invertible Moreover the eigenvalues of its inverse all have modulus less
or equal to 1, which implies the L2 stability of the scheme.
<k = 1 + (1 cos

Proposition 6 The centred scheme second order in space and first order in time:
un+1
= unj
j

at n
(u
unj1 ).
2x j+1

is unstable in L2
t
Proof. Let us denote by = a x
. The first order in time centred scheme becomes
n+1
n
in matrix form U
= AU where A is the circulant matrix with three non vanishing
2ik
diagonals corresponding to c0 = 1, c1 = cn1 = 2 e n . Hence its eigenvalues are
2ik

2ik

2k
k = 1 2 (e n e n ) = 1 i sin 2k
n so that |k | > 1 for all k such that sin n 6= 0.
Hence the scheme is unstable.

Theorem 1 (Lax) The scheme (2.8) is convergent if it is stable and consistent.


Proof. Let V (tn ) be the vector whose components are the exact solution at the grid
points at time tn . The, as the scheme is consistent, we have
LV (tn+1 ) = M V (tn ) + tO(tq + xp ).
Note E n = U n V (tn ) the vector containing the errors at each point at time tn , then
as LU n+1 = M U n , we have LE n+1 = M E n + tO(tq + xp ). Hence
E n+1 = L1 M E n + tK1 (tq + xp ),
= L1 M (E n1 + tK1 (tq + xp )) + tK1 (tq + xp ),
= (L1 M )n+1 E 0 + (1 + + (L1 M )n )tK1 (tq + xp ).
Hence
kE n+1 k k(L1 M )n+1 E 0 k + k(1 + + (L1 M )n )tK1 (tq + xp )k.

(2.12)

The stability implies that for any initial condition U 0 , as U n = (L1 M )n U 0 , we have
k(L1 M )n U 0 k KkU 0 k,
which means that k(L1 M )n k K for all n. Plugging this into (2.12), we obtain:
kE n+1 k KkE 0 k + nKK1 t(tq + xp ).
Then as on the one hand E 0 taking as an initial in the scheme U 0 = V (0), and on the
other hand nt T the maximal considered time, we have
kE n+1 k KK1 T (tq + xp ),
whence convergence.
Corollary 3 If the exact solution is of class C 3 , the schemes (2.5) and (2.7) converge.
Proof. This follows immediately from the Lax theorem by applying propositions 1 and
2.
11

2.1.6

High-order time schemes

When working with linear homogeneous equations with no source term, the simplest
way to derive high order time schemes is to use a Taylor expansion in time and plug in
the expression of the successive time derivatives obtained from the differential system
resulting from the semi-discretization in space. Consider for example that after semidiscretization in space using Finite Differences (or any other space discretisation method)
we obtain the differential systems

u0 (t)
dU

= AU, with U = ... ,


dt
un1 (t)
and A the appropriate matrix coming from the semi-discretization in space. Then a
Taylor expansion in time up to order p yields
U (tn+1 ) = U (tn ) + t

dU
tp dp U
(tn ) + +
(tn ) + O(tp+1 ).
dt
p! dtp

Now if A does not depend on time and

dU
dt

= AU , we get that

dp U
= Ap U, for any integer p.
dtp
Hence, denoting U n an approximation of U (tn ), we get a time scheme of order p using
the formula
U n+1 = U n + tAU n + +

tp p n
tp p n
A U = (I + tA + +
A )U .
p!
p!

(2.13)

For p = 1 this boils down to the standard explicit Euler scheme.


Writing U n the solution in vector form at time tn , we define the propagation matrix
A such that
U n+1 = AU n .
Definition 4 The numerical scheme defined by the propagation matrix A is stable if
there exists > 0 such that U n is bounded for all n N and t < .
Proposition 7 The numerical scheme defined by the propagation matrix A is stable if
there exists > 0 such that for all t < all eigenvalues of A are of modulus less or
equal to 1.
Stability of Taylor schemes. For a Taylor scheme of order p applied to dU
dt = AU ,
p . Then denoting by an eigenvalue of A, the
we have A = I + tA + + t
A
p!
p
corresponding eigenvalue of A is = 1 + t + + p t
p! . And one can plot the region
of the complex plane in which |(t)| 1 using for example ezplot in Matlab, which
are the stability regions. This means that the time scheme associate to the semi-discrete
form dU
dt = AU is stable provided all the eigenvalues of A are such that t is in the
stability region.

12

Examples.
1. The Upwind scheme:

dui (t)
dt

i1 (t)
= a ui (t)u
corresponds to the circulant matrix
x
2ijk

a
n ).
= c1 . So its eigenvalues verify k t = at
A with c0 = x
x (1 e
Obviously, for any integer value of k, k t is on a circle in the complex plane of
at
radius at
x centred at ( x , 0). The stability region of the explicit Euler method
is the circle of radius 1 centred at (1, 0), so that in this case we see again that
the scheme is stable provided at
x 1. For the higher order schemes the radius
corresponding to the maximal stability can be found by computing the second real
root (in addition to 0) of the equation |(t)| = 1, see Fig. 2.2. We find that
for the order 2 scheme = 2, so that the stability condition is the same as for
the order 1 scheme. For the order 3 scheme we find that = 2.5127 and for
the order 4 scheme we find that = 2.7853. The value of corresponds to the
diameter of the largest circle of eigenvalues that is still completely enclosed in the
||
stability region. This yields the stability condition at
x 2 . We notice that the
maximal stable time step is larger for the schemes of order 3 and 4.

Maximal stability, order 2

Maximal stability, order 3

3
4

3
4

1
x

1
x

Maximal stability, order 4

3
4

1
x

Figure 2.2: Location of eigenvalues (blue circle) corresponding to maximal stability zone
for explicit time schemes of order 2, 3, 4 (left to right).
i1 (t)
i (t)
2. The centred scheme: dudt
= a ui+1 (t)u
corresponds to the circulant matrix
2x
a
A with c1 = 2x
= cn1 . The corresponding eigenvalues are such that

k t =

2ijk
at 2ijk
iat
2jk
(e n e n ) =
sin
.
2x
x
n

Hence the eigenvalues are all purely imaginary and the modulus of the largest one
is at
x . The stability zones for the schemes of order 1 to 6 are represented in Fig.
2.3. Note that for the order 1 and 2 scheme the intersection of the stability zone
with the imaginary axis is reduced to the point 0. So that when all eigenvalues
are purely imaginary as is the case here, these schemes are not stable for any
positive t. On the other hand the schemes of order 3 and 4 have a non vanishing
stability zone on the imaginary axis, larger for the order 4 scheme. By computing
the intersection of the curve |(t)| = 1 with the
imaginary axis we find the
stability conditions for the order 3 scheme: at

3 and for the order 4 scheme


x

at
x 2 2.
13

Remark 2 The order 5 and 6 schemes are more problematic for eigenvalues of A
on the imaginary axis as the zooms of Figure 2.4 tell us. Even though there is a part
of the imaginary axis in the stability zone, there is also a part in the neighborhood of
0 which is not. Therefore small eigenvalues of A will lead to instability on longer
time scales. This is problematic, as unlike usual Courant condition instability
problems which reveal themselves very fast, this leads to a small growth in time.

1.0

1.5

1.0
0.5

1
0.5

0.0
2.0

1.5

1.0

0.5

0.0

0.0

2.0

1.0

1.5

2.5

0.0

0.5

2.0

1.5

1.0

0.5

0.0

0.5

0.5y
1.0

1.5

1.0

4.0

3
3

3.2
2

2.4

1.6
1

0.8

2.5

2.0

1.5

0.5

1.0

0.0

3.2

2.8

2.4

2.0

1.6

1.2

0.8

0.4

0.0

0.0

0.4

0
0.8
1.6

y
2

2.4
3.2

3
3

4.0

Figure 2.3: Stability zone for Taylor schemes. From top to bottom and left to right
order 1, 2, 3, 4, 5, 6.

2
0.5
1

0
10

0.0
0

5
1

10

103

15

20

1*104

5*105

0*100

x
0.5y

x
y
2

Figure 2.4: Stability zone for Taylor schemes. Zoom around imaginary axis. Left order
5, right order 6.

14

2.1.7

The Cauchy-Kowalevsky procedure

The Cauchy-Kowalevsky procedure consists in writing a Taylor expansion in time to


update the scheme from time tn to tn+1 and then to use the equation to replace the time
derivatives in the Taylor expansion by spatial derivatives. This is particularly simple for
the linear advection, but can be used in principle (with more complicated computations)
for any hyperbolic equation or system. Let us now apply this procedure: (we denote as
usual by uni = u(tn , xi ))
un+1
= uni + t
i
Then using the equation

u
t

t2 2 u
u
(tn , xi ) +
(tn , xi ) + O(t3 ).
t
2 t2

(2.14)

+ a u
x = 0 we get

u
u
(tn , xi ) = a (tn , xi ),
t
x

2
2u

u
u
2 u
(t
,
x
)
=
(a
)

a
=
a
(tn , xi ),
n
i
t2
t
x
x t
x2

and now we use centred schemes for the first and second order spatial derivatives of u:
un uni1
u
(tn , xi ) = i+1
+ O(x2 ),
x
2x

uni+1 2uni + uni1


2u
(t
,
x
)
=
+ O(x2 ).
n i
x2
x2

Plugging this into the Taylor expansion (2.14), we get:


un+1
= uni
i

at n
a2 t2 n
(ui+1 uni1 ) +
(u
2uni + uni1 ) + tO(t2 + x2 ).
2x
2x2 i+1

This scheme, which is second order in space and time, is known as the Lax-Wendroff
scheme.

2.2
2.2.1

The Finite Volume method


The first order Finite Volume schemes

Let us introduce the Finite Volume method on the generic scalar conservation law of the
form
u f (u)
+
= 0.
(2.15)
t
x
In the case of our linear advection equation, we have f (u) = au.
In the Finite Volume method, the computational domain is divided into cells (intervals in 1D) and the unknown quantity that is numerically computed is the cell average
of u on each cell. Recall that for Finite Differences the unknowns were the point values
of u at the grid points. We need to number the cells. In 1D a convenient way to do it
in order to avoid confusion with the grid points, is to assign half integers. Let us denote
by
Z xi+1
1
ui+ 1 (t) =
u(t, x) dx.
2
xi+1 xi xi
The Finite Volume numerical scheme is then obtained by integrating the original
equation on each cell of the domain. As for the time scheme there are at least two
classical ways, the first is to also integrate in time between tn and tn+1 and then use a
quadrature formula to compute the integral, for example, the left hand rectangle rule
15

yields a first order explicit formula. The second is to use the method of lines and separate
space discretisation from time discretisation. Then standard ODE discretisation schemes
can be used. This is what we shall do mostly in this lecture.
So integrating (2.15) on the cell [xi , xi+1 ] and dividing by xi+ 1 = xi+1 xi yields
2

dui+ 1 (t)
2

dt

1
(f (u(t, xi+1 )) f (u(t, xi ))) = 0.
xi+ 1
2

Here we see that a second ingredient is needed in order to define the algorithm. We
only know the cell averages of u, how do we define the value at u at the cell interfaces.
The simplest scheme, which is first order accurate in space consists in assuming that
u is constant on each cell and thus equal to its cell average. But it is not defined at
the cell interface. In order to complete the Finite Volume scheme we need to define a
so called numerical flux at each cell interface denoted by gi that needs to be consistent
with f (xi ), i.e. gi = f (u(xi )) + O(xp ) for some positive p. A numerical flux of order
2 is the centred flux gi = 12 (f (ui 1 ) + f (ui+ 1 )). This yields the following scheme for a
2
2
uniform x:
dui+ 1 (t) (f (ui+ 3 ) f (ui 1 ))
2
2
2
+
= 0.
dt
2x
Coupling it with an explicit Euler scheme in time this becomes, and applying it to the
linear advection (f (u) = au) we get
= uni+ 1
un+1
i+ 1
2

at
(u 3 ui 1 ).
2
2x i+ 2

(2.16)

We recognise here the centred Finite Difference scheme shifted to the cell centres. Remember that this scheme is unstable, so that it cannot be used in practice. In order
to get stable scheme, we need to introduce the notion of unwinding like for Finite Differences. This can be done very easily in the definition of the numerical flux by simply
choosing the value of u in the upwind cell only to define the numerical flux. We have
f (u)
u
0
x = f (u) x . This means that locally at each cell interface the direction of the transport is defined by the sign of f 0 (u) (in the case of the linear advection f 0 (u) = a and
the upwind direction is determined by the sign of a). So the upwind numerical flux is
defined by

ui 1 +ui+ 1

f (ui 1 ) if f 0 ( 2 2 2 ) 0
2
gi =
u 1 +u 1
f (ui+ 1 ) if f 0 ( i 2 2 i+ 2 ) < 0
2

Again, combining the Finite Volume scheme with an upwind flux and an explicit Euler
time discretisation yields for the linear advection with a > 0
un+1
= uni+ 1
i+ 1
2

at
(u 1 ui 1 ).
2
x i+ 2

(2.17)

We also recognise here the first order in time and space upwind scheme shifted to the
cell centres.
Remark 3 Using the midpoint rule
Z xi+1
1
ui+ 1 =
u(x) dx = u(xi+ 1 ) + O(x2 ).
2
2
x xi
16

The we can reinterpret the Finite Volume as a Finite Difference scheme at the cell
centres, which explains that we get the same formulas. However this is not true for
higher orders, for which Finite Volume and Finite Difference schemes are genuinely
different.

2.2.2

Higher order schemes

In order to get high order Finite Volume schemes, the idea is to reconstruct polynomials
of some given degree from the cell averages that are obtained with the Finite Volume
procedure. The main idea for doing this is to contract an interpolation polynomial for
the primitive of the polynomial we are looking for.
At time step tn we know unj+ 1 known average value of un on cell [xj , xj+1 ] of length
2

xj+ 1 = xj+1 xj . We want to construct a polynomial pm (x) of degree m such that


2

1
xj+ 1
2

xj+1

xj

pm (x) dx = unj+ 1 .
2

d
To this aim we look for pm (x) such that dx
pm (x) = pm (x). Then
Z xj+1
xj+ 1 unj+ 1 =
pm (x) dx = pm (xj+1 ) pm (xj ).
2

xj

Rx

u
n (x) dx a primitive of the piecewise constant function u
n with value
Pj
on [xj , xj+1 ]. Then W (xj+1 ) = k=1 hk+ 1 unk+ 1 and

Let W (x) =
unj+ 1

x0

W (xj+1 ) W (xj ) = xj+ 1 unj+ 1 = pm (xj+1 ) pm (xj ).


2

Then we take for pm an interpolating polynomial at points xj of W so that


Z xj+1
1
1
pm (x) dx =
(
pm (xj+1 ) pm (xj ))
xj+ 1 xj
xj+ 1
2

1
=
(W (xj+1 ) W (xj )) = unj+ 1 .
2
xj+ 1
2

There are many ways to choose an interpolating polynomial, one could use spline interpolation or Hermite interpolation, but the simplest and most used choice is to use a
Lagrange interpolation polynomial. This being said, a Lagrange interpolating polynomial of degree k is defined with k + 1 interpolation points. So we need to use as many
values in neighbouring cells as needed.
In order to reconstruct a polynomial of a given degree in a given cell there are many
possible stencils, i.e. ensembles of cells, that can be used. For the reconstruction of a
polynomial of degree k exactly k average values corresponding to k neighbouring cells
are needed. The only constraint is that the value on the cell where the polynomial being
reconstructed is used. High-order methods are prone to oscillations especially around
discontinuities. So one good idea is to use the stencil which minimises the oscillations.
This can be easily done by choosing automatically the stencil based on the Newton
divided differences which can be used to construct the interpolating polynomial. This
method is called ENO (Essentially Non Oscillatory). See for example [4] for a detailed
description.
17

The ENO method can be still improved by taking all possible stencils but putting a
weight on each of the polynomials obtained. This is called the WENO method (Weighted
Essentially Non Oscillatory) A good review of this technique is given in [7].

2.3

The Finite Element method

The Finite Element method is mathematically more involved that the previous methods.
The discretisation idea is based on the Galerkin principle of approximating the function
space in which the solution is defined by a finite dimensional subspace. It uses also a
weak (or variational form of the equation). In order to obtain the weak form, the idea
is to multiply by a smooth test function and integrate over the whole domain, with a
possible integration by parts to eliminate the highest derivatives.
Let us describe it on the advection-diffusion problem (assuming periodic boundary
conditions on the domain [0, L]):
u
u
2u
+a
2 = 0.
t
x
x
Multiplying by a test function v which does not depend on t and integrating, with an
integration by parts in the last integral and periodic boundary conditions, yields
Z
Z L
Z L
d L
u
u v
uv dx + a
v dx +
dx = 0.
dt 0
0 x
0 x x
The natural space in which this formulation is defined is H 1 ([0, L]) = {u L2 ([0, L]) | u
x
L2 ([0, L])}. The variational problem thus reads Find u H 1 such that
d
dt

Z
uv dx + a

u
v dx +
x

u v
dx = 0 v H 1 ([0, L]).
x x

Now in order to define a Finite Element approximation, we need to construct a sequence


of subspaces Vh of H 1 ([0, L]) which is dense in H 1 ([0, L]) (in order to get convergence).
One subspace is of course enough to compute a numerical approximation. There are
many ways to construct such subspaces. The classical Lagrange Finite Element method
consists in building a mesh of the computational domain and to assume that the approximating function is polynomial say of degree k in each cell. For the piecewise polynomial
function space to be a included in H 1 ([0, L]) the only additional requirement is that
the functions are continuous at the cell interfaces. So the subspace Vh on the mesh
x0 < x1 < < xn1 is defined by


Vh = vh C 0 ([a, b]) | vh|[xi ,xi+1 ] Pk ([xi , xi+1 ]) .
In order to express a function vh Vh we need a basis of Vh . This can be constructed
easily combining bases of Pk ([xi , xi+1 ]). In order to impose the continuity requirement at
the cell interface, the simplest is to use a Lagrange basis Pk ([xi , xi+1 ]) with interpolation
points on the edge of the intervals. Given k + 1 interpolation points xi = y0 < y1 <
< yk = xi+1 the Lagrange basis functions of degree k denoted by lk,i , 0 i k, are
the unique polynomials of degree k verifying lk,i (yj ) = i,j . Because of this property,
P
any polynomial p(x) Pk ([xi , xi+1 ]) can be expressed as p(x) = kj=0 p(yj )lj,k (x) and
conversely any polynomial p(x) Pk ([xi , xi+1 ]) is uniquely determined by its values
18

at the interpolation points yj , 0 j k. Hence in order to ensure the continuity


of the piecewise polynomial at the cell interface xi it is enough that the values of the
polynomials on both sides of xi have the same value at xi . This constraint removes one
degree of freedom in each cell, so that the total dimension of Vh is nk and the functions
of Vh are uniquely defined in each cell by their value at the degrees of freedom (which
are the interpolation points) in all the cells. The basis functions denoted of Vh denoted
by (i )0jnk are such that their restriction on each cell is a Lagrange basis function.
Note that for k = 1, corresponding to P1 finite elements, the degrees of freedom are
just the grid points. For higher order finite elements internal degrees of freedom are
needed. For stability and conveniency issues this are most commonly taken to be the
Gauss-Lobatto points on each cell.
In order to obtain a discrete problem that can be solved on a computer, the Galerkin
procedure consist in replacing H 1 by Vh in the variational formulation. The discrete
variational problem thus reads: Find uh Vh such that
d
dt

Z
uv dx + a
0

u
v dx +
x

Z
0

u v
dx = 0 vh Vh .
x x

P
Now expressing uh (and vh ) in the basis of Vh as uh (t, x) = nk
j=1 uj (t)j (x), vh (x) =
Pnk
j=1 vj j (x) and plugging these expression in the variational formulation, denoting by
U = (u0 , u1 , , unk )T and similarly for V yields: Find U Rnk such that
d X
uj vi
dt
i,j

i (x)j (x) dx + a
0

j (x)
i (x) dx
x
0
i,j
Z L
X
i (x) j (x)
+
dx = 0 V Rnk ,
uj vi
x
x
0
X

uj vi

i,j

which can be expressed in matrix form


V T (M

dU
+ DU + AU ) = 0 V Rnk ,
dt

which is equivalent to
dU
+ DU + AU = 0
dt
where the square nk nk matrices are defined by
M

Z
M =(
0

j (x)i (x) dx)i,j ,

Z
D=(
0

j (x)
i (x) dx)i,j ,
x

Z
A=(
0

i (x) j (x)
dx)i,j .
x
x

Note that these matrices can be computed exactly as they involve integration of
polynomials on each cell. Moreover because the Gauss-Lobatto quadrature rule is exact
for polynomials of degree up to 2k 1, both A and D can be computed exactly with the
Gauss-Lobatto quadrature rule. Moreover, approximating the mass matrix M with the
Gauss-Lobatto rule introduces an error which does not decrease the order of accuracy
of the scheme [2] and has the big advantage of yielding a diagonal matrix. This is what
is mostly done in practice.

19

Matrix Assembly. Usually for Finite Elements the matrices M , D and A are computed from the corresponding elementary matrices which are obtained by change of
variables onto the reference element [1, 1] for each cell. So
L

i (x)j (x) dx =
0

=0

and doing the change of variable x =


Z

x+1

n1
X Z x+1

i (x)j (x) dx,

x+1 x
x

x+1 +x
,
2

x+1 x
i (x)j (x) dx =
2

we get

(
x) (
x) d
x,
1

where (
x) = i ( x+12x x
+ x+12+x ). The local indices on the reference element go
from 0 to k and the global numbers of the basis functions not vanishing on element
are j = k + . The are the Lagrange polynomials at the Gauss-Lobatto points in
the interval [1, 1].
The mass matrix in Vh can be approximated with no loss of order of the finite
element approximation using the Gauss-Lobatto quadrature rule. Then because the
products (
x) (
x) vanish for 6= at the Gauss-Lobatto points by definition of the
which are the Lagrange basis functions at these points, the elementary matrix M is
diagonal and we have
Z

(
x)2 d
x

k
X

x )2 = wGL
wGL (

=0

using the quadrature rule, where wGL is the Gauss-Lobatto weight at Gauss-Lobatto
= diag(wGL , . . . wGL ) is the matrix with k + 1
point (
x ) [1, 1]. So that finally M
0
k
lines and columns with the Gauss-Lobatto weights on the diagonal.
Let us now compute the elements of D. As previously we go back to the interval
[1, 1] with the change of variables x = x+12x x
+ x+12+x and we define (
x) =
x+1 x
x+1 +x
i (
x
+
). Note that a global basis function i associated to a grid point
2
2
has a support which overlaps two cells and is associated to two local basis functions.
Thus one needs to be careful to add the two contributions as needed in the final matrix.
We get 0 (
x) = x+12x 0i ( x+12x (
x + 1) + x ). It follows that
Z

x+1

0j (x)i (x) dx

=
1

2
x+1 x
0 (
x) (
x)
d
x=
x+1 x
2

0 (
x) (
x) d
x.

The polynomial (
x) is of degree k so that 0 (
x) is of degree k 1 so that the GaussLobatto quadrature rule with k + 1 points is exact for the product which is of order
2k 1. Using this rule
Z

0 (
x) (
x) d
x=

k
X

GL 0
wm
(
xm ) (
xm ) = wGL 0 (
x ),

m=0

As before, because are the Lagrange polynomials at the Gauss-Lobatto points, only
the value at x in the sum is one and the others are 0. On the other hand evaluating

20

the derivatives of the Lagrange polynomial at the Gauss-Lobatto points at these GaussLobatto points can be done using the formula
0 (
x ) =

X
p /p
for 6= and 0 (
x ) =
0 (
x ),
x
x

6=

where p = 6= (
x x
). This formula is obtained straightforwardly by taking the
derivative of the explicit formula for the Lagrange polynomial
Y
(
xx
)
6=

(
x) = Y

(
x x
)

6=

and using this expression at the Gauss-Lobatto point x


6= x
. We refer to [1] for a
detailed description.
We can now conclude with the computation of the stiffness matrix A. Having already
computed the expression of the change of variable of the derivatives we can quickly go
to the result. We have in each element
Z

x+1

0j (x)0i (x) dx

=
1

2
x+1 x

2
x+1 x

2

0 (
x)0 (
x)

x)0 (
x) d
x=
0 (

x+1 x
d
x
2

k
X
2
GL 0
xm )0 (
xm ).
wm
(
x+1 x
m=0

As the polynomial being integrated is of degree 2(k 1) = 2k 2 the Gauss-Lobatto


quadrature is exact. Here no simplifications occurs in the sum, which has to be computed.
Still the expressions of the derivatives at the Gauss-Lobatto points computed above can
then be plugged in.
Time advance and stability. At the end we get as for the finite difference and finite
volume methods a system of differential equations that can be solved with any ODE
solver. Let us make a few remarks concerning the stability of the scheme (once discretised
in time). As we saw previously this depends on whether the eigenvalues of the matrix
A is included in the stability zone of the ODE solver. Here A = M 1 (D + A). Note
that the matrices M and A are obviously symmetric and thus have only real eigenvalues.
On the other hand, for periodic boundary conditions and integration by parts, yields
R j (x)
R
i (x)
that
x i (x) dx = j (x) x dx. Hence D is skew symmetric and has only
imaginary eigenvalues. Remember that the stability zones of our explicit ODE solvers lie
mostly on the left-hand side of the imaginary axis in the complex plane, and that only
the third and fourth order schemes have a stability zone on the imaginary axis. Hence
in the pure advection case = 0 A has purely imaginary eigenvalues and the order one
and two time schemes are always unstable. In order to stabilise the method a procedure
that is often used with a finite element discretisation of a pure advection problem is to
add a diffusive term that goes to 0 with the cell size, i.e take = x, in this case a
small negative real part is added to the eigenvalues which are thus pushed into the left
half of the complex plane and the stability zone is enhanced.

21

2.4

The Discontinuous Galerkin (DG) method

The DG method represents the unknowns like the Finite Element method by piecewise
polynomial functions, but unlike Finite Element methods the polynomials are discontinuous at the cell interfaces and a numerical flux is defined at the cell interface in the
same way as for Finite Volume methods.
So on each cell the discrete unknown uh is represented as a linear combination of well
chosen basis functions of the space of polynomials of degree k Pk . The dimension of this
space is k +1. As no continuity is enforced at the element interface, there is no constraint
on the basis functions and generally two kinds of basis functions are used: either the
Legendre polynomials which form an orthogonal basis, we then speak of modal DG or
one can use a Lagrange basis defined on a set of interpolation points within the cell we
then speak of nodal DG. The interpolation points are generally chosen to be either the
Gauss points or the Gauss-Lobatto points which are both very convenient as they allow
to express the integrals appearing in the formulation exactly (for Gauss) or almost (for
Gauss-Lobatto) using the associated quadrature rule.
For the derivation of a DG scheme, the equation is multiplied by a polynomial test
function on each cell and and integration by parts is used so that a boundary term
appears which will allow the coupling between two neighbouring cells. Let us apply it
here to the conservation law
u f (u)
+
= 0.
t
x
We then get
Z

x+1

x+1

f (u)
v dx
x
x
x
Z
Z x+1
d x+1
v
=
uv dx
f (u)
dx + (f (u(x+1 ))v(x+1 ) f (u(x ))v(x )) = 0.
dt x
x
x
(2.18)

d
dt

uv dx +

The DG method has then two key ingredients.


1. Choose on each cell a finite dimensional representation, usually by a polynomial
as said previously.
2. Define a unique numerical flux denoted by g = g(uL (x ), uR (x ))) at the cell
interface which is constructed from the two values coming from the cells sharing
the interface on the left and on the right. This is defined using similar recipes as
for the finite volume method.
Choosing as local representation for u and the test function v the Lagrange polynomials at the Gauss-Lobatto points simplifies the computation of the fluxes, as in this case
only the Lagrange polynomial associated to the edge node does not vanish at the edge.
This situation is different when using Legendre polynomials or Lagrange polynomials at
only interior nodes (like the Gauss points). Note however that Legendre polynomials
have the advantage of having exactly a diagonal mass matrix. This is obtained also
with Lagrange polynomials at the Gauss-Lobatto points but in this case at the price of
a small quadrature error.
As opposite to the Finite Element method, only local matrices on each element,
in practice only the elementary matrices on the [1, 1] reference interval need to be
22

on cell on the reference interval has the


assembled. The elementary mass matrix M
components
Z 1

(
x) (
x) dx, 0 , k.
M, =
1

When the basis functions are the Legendre polynomials which form an orthonormal
basis.
as in the previous section using the Gauss-Lobatto quadrature
= diag(wGL , . . . wGL ).
M
0
k
From this matrix, the local mass matrix M+ 1 on cell [x , x+1 ] can be expressed as
2

x+1 x
M.
2

M+ 1 =
2

Let us denote by K+ 1 the local matrix containing the derivative of v. In order to get
2
an expression for the components of K+ 1 we introduce the local basis functions and
2
compute using again the affine change of variable to the reference interval [1, 1]:
Z x+1
Z 1
Z 1
i
2
x+1 x
f (u)
dx =
f (u(x))
d
x=
f (u(x))0 (
x) d
x.
0 (
x)
x
x

x
2
+1

x
1
1
Then using the Gauss-Lobatto quadrature rule this becomes
Z

f (u(
x))0 (
x) d
x

k
X

x ))0 (
x ) =
wGL f (u(

k
X

x ),
wGL f (u )0 (

=0

=0

where u = u(
x ) is the th component of u on the Lagrange basis. Denoting U+ 1 =
2
(u0 , u1 , . . . , uk ) on the cell [x , x+1 ], thus defining the component of matrix K+ 1 at
2

x ), we get that
line and column as being wGL 0 (
Z

f (u(
x))0 (
x) d
x

k
X

x ) = (K+ 1 f (U+ 1 )) ,
wGL f (u )0 (
2

=0

where f (U+ 1 ) = (f (u0 ), f (u1 ), . . . , f (uk )).


2

Remark 4 Because the Gauss-Lobatto quadrature is exact in this case, we notice that
the matrix K+ 1 is exactly the matrix associated to the compotents
2

x+1

j (x)
x

i
dx =
x

(
x))0 (
x) d
x = wGL 0 (
x ).

We also notice that this matrix does not depend on the specific interval and is equal to
= K 1 for all .
the matrix on the reference element K
+
2

Now plugging all this into the formula (2.18) we get on each cell
T
V+
1 M+ 1
2

dU+ 1
2

= V+
1 Kf (U+ 1 ) (g+1 vk g v0 ).
2

23

Then introducing the vector G+ 1 Rk whose only non zero components are the first
2
which is g and the last which is g+1 , we get the following system of ODE
x+1 x dU+ 21
(U 1 ) + G 1 .
M
= Kf
+ 2
+ 2
2
t
The numerical flux g depends on values of u coming from the neighbouring cell, this is
being diagonal there is
where the coupling between the cells takes place. The matrix M
no linear system to solve. Simple examples of fluxes in the linear case f (u) = au are the
same as for the finite volume method with the centred or upwind fluxes, the two values
being used here are the values of u on the interface coming from the two cells sharing
the interface, this will be the local value of uk from the left cell and the local value of
u0 from the right cell.

24

Chapter 3

Linear systems
Let us know consider linear systems of conservation laws in 1D. This can be written in
the general form
U
U
+A
= 0,
t
x
where U (t, x) is a vector in Rn and A a given matrix with constant coefficients. We will
focus in the following on the Discontinuous Galerkin method, which includes the first
order Finite Volume method when polynomials of degree 0, i.e. constants are taken on
each cell.

3.1

The Riemann problem

The main numerical issue when constructing a Finite Volume or Discontinuous Galerkin
scheme is to find a good numerical flux that is consistent (i.e. converges towards the
exact flux when the cell size goes to 0) and stable. As we saw previously in the linear
scalar case enhanced stability is given by upwinding. We now need to generalise the idea
of upwind to the case of systems.
The construction of a numerical flux is a local procedure at the interface between
two cells, where a different value is given on the left side and on the right side from
the polynomial approximation (or reconstruction for Finite Volumes). In order to get
information from the equation itself the idea is to solve it locally using an initial condition which is a step function. The Riemann problem is the corresponding initial value
problem:
U
U
+A
= 0,
t
x

U if x < 0,
U (0, x) = L
UR if x 0,
where UL and UR are two given constant vectors.
The system being hyperbolic implies that A has real eigenvalues and can be diagonalised. Hence A = P P 1 , where is the diagonal matrix containing the eigenvalues.
Then introducing the so-called characteristic variables V = P 1 U , and multiplying the
system by P 1 on the left we get
P 1

U
U
V
V
+ P 1 AP P 1
=
+
= 0.
t
x
t
x
25

So in the variable V the system is diagonal and reduces to the set of linear advection
equations
vi
vi
+ i
= 0, 1 i n
t
x
where the vi are the components of V and the i the eigenvalues of A. The exact
solution of these equations is given by vi (t, x) = vi (0, x i t), where the vi (0, x) are the
components of the initial vector which take the constant values VL = P 1 UL if x < 0
and VR = P 1 UR if x 0. In other terms

v
if x < i t,
vi (t, x) = i,L
vi,R if x i t.
In practice we want to use the Riemann problem to determine the value of V (and U ) at
the cell interface, corresponding to x = 0, the discontinuity point at any strictly positive
time. And we deduce from the previous solution that

v
if 0 < i ,
vi (t, 0) = i,L
vi,R if 0 i .
In order to get a vector expression, we introduce the diagonal matrices + where the
negative eigenvalues are replaced by 0 and where the positive eigenvalues are replaced
by 0. Obviously = + + . Then for t > 0 we have
V (t, 0) = + V (t, 0) + V (t, 0) = + VL + VR ,
as for all positive eigenvalues the corresponding component of V (t, 0) is vi,L and for
all negative eigenvalues the corresponding component of V (t, 0) is vi,R . Note that as
V (t, 0) is multiplied by the components of V (t, 0) corresponding to 0 eigenvalues do
not need to be considered as they are multiplied by 0 anyway. So the side where the
strict inequality is used for the initial condition of the Riemann problem plays no role.
Denoting by A+ = P + P 1 and A = P P 1 the flux AU (t, 0) associated to the
solution of the Riemann problem at the cell interface can also be expressed conveniently
directly in terms of U
AU (t, 0) = P + V (t, 0) + P V (t, 0) = P + VL + P VR = A+ UL + A UR .
This expression AU (t, 0) = A+ UL + A UR can be used to define the numerical flux
at the cell interface, using the value UL coming from the left-hand side of the interface
and UR coming from the right-hand side of the interface. For actual computations,
the matrices A+ and A need to be computed explicitely from the eigenvalues and
eigenvectors of the matrix A. Notice that in the case of a scalar equation the matrix
A is reduced to the scalar a which is then obviously the only eigenvalue of the 1 1
matrix and if a > 0 we have A+ = a and A = 0, so that the numerical flux becomes
au(t, 0) = auL and the same way if a < 0 au(t, 0) = auR , so that the numerical flux
obtained from the solution of the Riemann problem reduces to the upwind flux.
Example.
less units:

We consider the 1D Maxwell equations which can be written in dimensionE B


+
= 0,
t
x
B E
+
= 0.
t
x
26

This can be written in the form of a linear system


 
U
U
E
+A
= 0, with U =
,
B
t
x



0 1
A=
.
1 0

The eigenvalues of A are the solutions of det(AI) = 0, i.e. 2 = 1. So the eigenvalues


are 1 = 1 and 2 = 1. They are real and distinct so that the system is strictly
hyperbolic. Let Vi be a normalised eigenvector associated to the eigenvalue i , i = 1, 2.
We have AV1 = V1 so that V1 = 12 (1, 1)T and AV2 = V2 so that V1 = 12 (1, 1)T . We
define P the matrix whose columns are V1 and V2 . P is obviously orthonormal, so that
its inverse is its transpose. Then we have P A = P . So that we can define:






1 1 1
1 1 1
0 0
1 1
T
A+ = P + P =
=
,
0 1
1 1
2 1 1
2 1 1






1 1 1
1 1 1
1 0
1 1
T
=
A = P P =
0 0
1 1
2 1 1
2 1 1
Hence, the upwind flux is given by
1
AU (t, 0) = A+ UL + A UR =
2


UL,1 + UL,2 + (UR,1 + UR,2 )
.
UL,1 + UL,2 + (UR,1 UR,2 )

Remark 5 As for the scalar problem, the numerical flux can be taken as a linear combination of the centred flux and the upwind flux (solution the Riemann problem):
1
Gj = A(UL + UR ) + (1 )(A+ UL + A UR ),
2

3.2

0 1.

The discontinuous Galerkin method

The discontinuous Galerkin method can be generalised to a system in a straightforward


manner. Each component of the approximate solution vector that we shall denote by Uh
is defined locally in each cell as a polynomial of degree k, i.e. an element of Pk [x , x+1 ].
Denoting by 0 , . . . , k a basis of Pk [x , x+1 ], the restriction of Uh to the cell [x , x+1 ],
+ 21

that we shall denote by Uh


+ 21

Uh

(t, x) =

k
X

can be expressed as

+ 12

j (x)Uj

+ 1
u 2 (t)
j,1 .

+ 21
.
..
with Uj (t) =

(t),

+ 1

j=0

uj,d 2 (t)
+ 1

The unknown function Uh (t, x) is completely determined by its components Uj 2 (t).


In order to determine those we proceed like in the scalar case and plug this expression in
the equation, multiply by each of the basis functions and integrate the derivative term
by part to let the flux through the interface appear:
+ 21

dUj

dt

(t)

x+1

x+1

i (x)j (x) dx
x

+ 21

AUj

(t)0i (x)j (x) dx


+ G+1 i (x+1 ) G i (x ) = 0.

27

The numerical flux G is a consistent approximation of the real flux AU(t, x ) at the
cell interface x which is identical for both cells sharing x .
Choice of the numerical flux: As in the scalar case, there are many possible
choices of the numerical flux.
1. Centred
flux: Use an average from the fluxes coming from both cells G =

1
2 A U 1 (t, x ) + U+ 1 (t, x ) .
2

2. Solution of the Riemann problem corresponding to UL = U 1 (t, x ) and UR =


2
U+ 1 (t, x ). This is the generalisation of the upwind flux to systems. Better
2
stability properties but more diffusive.
3. Linear combination of the above to combine good stability properties of the latter
and lower diffusion properties of former.
4. Generalised Lax-Friedrichs flux (also called Rusanov flux).
G =



i
1h 
A U 1 (t, x ) + U+ 1 (t, x ) U+ 1 (t, x ) U 1 (t, x ) ,
2
2
2
2
2

with = max |k |, where the k are the eigenvalues of A. This solver has good
stability properties without needing to solve the Riemann problem. Especially
interesting in the non linear case when the Riemann problem is hard to solve.

28

Chapter 4

Non linear conservation laws


4.1

Characteristics

We consider a generic scalar conservation law of the form


u f (u)
+
= 0,
t
x

with the initial condition u(0, x) = u0 (x).

Assuming f C 1 (R) and denoting by a(u) = f 0 (u), u also verifies


u
u
+ a(u)
= 0.
t
x
We now define the characteristic curves (or characteristics) associated to the conservation law the solution of the differential equation
dX
= a(u(t, X)).
dt
Then using the chain rule we have
d
u
dX u
u
u
(u(t, X(t)) =
(t, X(t))+
(t, X(t)) =
(t, X(t))+a(u(t, X(t)) (t, X(t)) = 0,
dt
t
dt x
t
x
if u is a solution of the equation. From this we deduce that u(t, X(t)) = u(0, X(0)) =
u0 (X(0)) which is independent of t. It follows that the characteristic curves are the
straight lines in the (t, x) plane of equation:
X(t) = X(0) + a(u0 (X(0)))t.
And it follows that the solutions of the conservation law satisfy u(t, X(t)) = u0 (X(0)).
This allows us to get the solution at a given point x and time t if the characteristic curve
can be traced back in the (t, x) plane to the line t = 0. This is not always the case when
f is non linear.
Remark 6 In the linear case we have f (u) = au, then the characteristics are the parallel
lines of slope a. They obey the equation X(t) = X(0) + at. So they never cross and so
taking X(t) = x, we have X(0) = x at, and we recover the classical solution of the
linear advection equation
u(t, x) = u0 (x at).
29

Example. Burgers equation. It corresponds to f (u) = 21 u2 , so that a(u) = f 0 (u) = u.


Then the characteristic X(t; x) such that X(0) = x satisfies X(t; x) = x + tu0 (x). Hence
if we consider x2 > x1 , we have
X(t; x2 ) X(t; x1 ) = x2 x1 + t(u0 (x2 ) u0 (x1 )).
Then if u0 is non decreasing we have X(t; x2 ) > X(t; x1 ) for all positive times, but if u0 is
x1
strictly decreasing then X(t; x2 ) = X(t; x1 ) for t = u0 (xx12)u
. So the characteristics
0 (x2 )
can cross and in this case the method of characteristics cannot be used to compute the
solution which is then no more C 1 .

4.2

Weak solutions

4.2.1

Definition

As we have seen in the case of Burgers equation, the characteristics associated to a non
linear conservation law can cross, in which case the method of characteristics can non
longer be used to find a solution. Actually when this happens the solution is no longer of
class C 1 and can even become discontinuous. We thus need another form of the equation
f (u)
u
u + x = 0 associated to the initial condition u(0, x) = u0 (x) to make it well defined
for such absolution.
Let us first recall the definition of a classical solution:
Definition 5 Let u0 C 1 (R), then u is called a classical solution if u C 1 (R+ R)
and u satisfies
u f (u)
+
= 0,
and U (0, x) = u0 (x)
u
x
in the classical sense.
Now we can define the notion of weak solution:
Definition 6 Let u0 L (R), then u L (R+ R) is called a weak solution of our
scalar conservation law if u satisfies
Z
0

u
+ f (u)
dt dx +
t
x

u0 (x)(0, x) dx = 0

Cc1 ([0, T [R.

(u)
u
Multiplying the equation u
+ fx
= 0 by and integrating by parts, it is easy to
1
verify that a classical C solution is also a weak solution. So the notion of weak solutions
generalises the notion of classical solutions.

4.2.2

The Rankine-Hugoniot condition

In most physical cases the solution is actually piecewise C 1 and there is only one (or a
few) lines of discontinuities in the space-time plane. The Rankine-Hugoniot condition
gives a constraint on the discontinuity along a line for the piecewise smooth solution to
be a weak solution of the equation.
Theorem 2 Assume the half space R+ R is split into two parts M1 and M2 by a smooth
curve S parametrized by (t, (t)) with C 1 (R+ ). We also assume that u L (R+ R)
30

1 ) and u2 = u|M C 1 (M
2 ) with u1 and u2 two classical
and that u1 = u|M1 C 1 (M
2
solutions of our equation in respectively M1 and M2 .
Then u is a weak solution if and only if
[u1 (t, (t)) u2 (t, (t))] 0 (t) = f (u1 (t, (t))) f (u2 (t, (t)))

t R+ .

(4.1)

Relation (4.1) is called the Rankine-Hugoniot condition. It is often written in the


simplified manner
(u1 u2 )s = f (u1 ) f (u2 ),
where s = 0 (t) is the propagation speed of the discontinuity.
Proof. Assume u is a weak solution of our equation. Then by definition
Z

u
+ f (u)
dt dx +
t
x

u0 (x)(0, x) dx = 0 Cc1 ([0, T [R.

Then we can split the first double integral into an integral on M1 and an integral on M2
and integrate these integrals by parts as u is C 1 on each of these domains. We then get
using Greens divergence formula, with denoting the outward unit normal vector,
Z
Z
Z
v dx = v dx
v ds

that
Z

u
+ f (u)
dt dx =
x
M1 t

M1

u f (u)
+
t
x

Z
dt dx +
M1


u
1 ds
f (u)

Because u = u1 is a classical solution in M1 the first term on the right hand side vanishes.
The boundary term is composed of a part on the t = 0 line which cancels with the part of
the integral on the initial condition which is in M1 and a part that can be parametrized
using . A tangent vector to the parametrized curve S is given by (1, 0 (t)). Hence the
outward unit normal is given by
 0 
1
(t)
1 = p
.
0
2
1
1 + (t)
Recall that the integral over a curve parametrized by : [a, b] Rn is defined by
Z

Z
F ds =

F ((t))k(t))k

2 dt.
a

Thus the part of the boundary integral corresponding to the integral on S which is
parametrized by (t) = (t, (t)) becomes

Z T
 0

1
u1
p
1 ds =
(t)u1 (t, (t)) f (u1 (t, (t))) (t, (t))
f (u1 )
1 + 0 (t)2
0
p
1 + 0 (t)2 dt,
Z T
 0

=
(t)u1 (t, (t)) f (u1 (t, (t))) (t, (t)) dt.

Z 
S

31

In the same way for the part on M2 for which the outward normal is 2 = 1 we get

Z T
Z 
 0

u1
(t)u2 (t, (t)) + f (u2 (t, (t))) (t, (t)) dt.
1 ds =
0
S f (u1 )
Adding the two pieces we find
Z T
 0

(t)(u1 (t, (t)) u2 (t, (t))) + f (u2 (t, (t))) f (u1 (t, (t))) (t, (t)) dt = 0.
0

This is true for all C 1 functions , hence its factor in the integral vanishes. This corresponds to the Rankine-Hugoniot relation.
Conversely, the same calculation proves that if the Rankine-Hugoniot relation is
satisfied then u is a weak solution.
Example:

Consider the Burgers equation which corresponds to f (u) =

u2
2 :

u 1 u2
+
= 0,
t
2 x
with the initial condition


0 if x < 0,
u0 (x) =
1 if x 0.

1. Using the Rankine-Hugoniot relation let us see under what condition a piecewise
constant solution can be a weak solution of this problem. Because the characteristics of the Burgers equation are straight lines, it is natural to look for piecewise
smooth solutions which are separated by a straight line in the t x plane. A
straight line can be parametrized by (t, t). Then the candidate piecewise solution
writes

0 if x < t,
u(t, x) =
1 if x t.
As (t) = t, we have s = 0 (t) = and so the Rankine-Hugoniot condition
(u1 u2 )s = f (u1 ) f (u2 ) becomes
1
(0 1) = (0 1)
2
and thus the only weak solution consisting of two constant states is u corresponding
to = 12 . Such a discontinuous solution is called a shock.
2. Let us now define


0 if x < 0,

u(t, x) = x/t if 0 x < t,
1 if x t.

The two constant values are obviously classical solutions in their part of the dox
u2
2x
main. On the other hand u(t, x) = x/t verifies u
t (t, x) = t2 and x (t, x) = t2
so that
1 u2
x
1 2x
u
(t, x) +
(t, x) = 2 +
=0
t
2 x
t
2 t2
and u(t, x) = x/t is also a classical solution on its domain of definition. It is
straightforward to verify that this solution is continuous on the lines x = 0 and
32

x = t, so that the Rankine-Hugoniot conditions are automatically verified (no


jump) but not C 1 . This solution which goes continuously from one constant state
to the other, but is not C 1 on two lines is called a rarefaction wave.
So we proved that this is also a weak solution of the Burgers equation, which
means that the equation we are considering in this example has at least two weak
solutions. In practice it is possible to construct other weak solutions with several
constant states separated by a line of discontinuity (see [3] for examples).

4.2.3

Entropy solution

The last example shows us that the uniqueness of a weak solution is not guaranteed.
However the physical solution is unique. This means that there is a piece missing in our
theory that will enable us to characterise the physically correct weak solution. The idea
that will lead us there is that an exact conservation law is unlikely to be physical. There
is always a small amount of dissipation. Mathematically this can be modelled by adding
a small diffusion term to our conservation law: for a small positive  we consider
u f (u)
2u
+
 2 = 0.
t
x
x

(4.2)

Associated to a smooth initial condition this equation has a unique smooth solution. And
the physically correct solution of our conservation law, will be the unique limit of this
solution (which depends on ) when  0. This unique solution can be characterised
using the notion of entropy.
Definition 7 Let U, F C 2 (R) such that
(i) U is strictly convex,
(ii) F 0 = U 0 f 0
Then U is called an entropy and F an entropy flux associated to the conservation law
u f (u)
+
= 0.
t
x
f (u)
Theorem 3 Assume the conservation law u
t + x = 0 can be associated to an entropy
U with entropy flux F . Then if (u ) is a sequence of smooth solutions of (4.2) such
that ku kL is bounded uniformly in  and u u almost everywhere then u is a weak
solution of the conservation law and verifies the entropy condition

U (u) F (u)
+
0
t
x

(4.3)

in the weak sense, i.e.


Z

(U (u)
0

+ F (u) ) dt dx 0
t
x

Cc1 (]0, T [R), 0.

Definition 8 A weak solution that verifies the entropy condition (4.3) is called an entropy solution.

33

It has been proven by Kruzhkov (see [3] and references therein) that for a bounded
initial condition u0 and a smooth flux f (u), there exists a unique entropy solution of the
scalar conservation law.
In the case of a strictly convex flux (f 00 (u) > 0) which is verified by the Burgers
equation, there is a simple characterisation of the entropy condition for a shock, which
is called the Lax entropy condition:
f 0 (uL ) > s > f 0 (uR ),

with s =

f (uR ) f (uL )
.
uR uL

This is the Rankine-Hugoniot setting and s is the shock speed. This tells us that the
shock is the entropy solution provided the characteristics lines cross at the shock line
x = st. Else the entropy solution is a rarefaction wave.
Remark 7 Because f is strictly convex f 0 is increasing and the entropy condition can
be satisfied only if uL > uR and in this case because of the strict convexity, s which is the
slope of the line joining (uL , f (uL )) and (uR , f (uR )) on the graph of f lies necessarily
between f 0 (uR ) and f 0 (uL ).
2

1 u
Example: Consider again the Burgers equation u
t + 2 x = 0 with piecewise constant
initial conditions u0 (x) = uL for x < 0 and u0 (x) = uR for x 0. We compute previously
the characteristics for Burgers equation which are given by X(t; x) = x + tu0 (x). So in
our case we have two families of characteristics, those issued from the negative values of
x at t = 0 which become X(t; x) = x + tuL and those issued from the positive values of x
at t = 0 which become X(t; x) = x + tuR . The flux of the Burgers equation f (u) = 21 u2
is strictly convex, so the Lax criterion applies and the Lax entropy condition can be
satisfied only if uL > uR which is the case when the characteristics cross.

Remark 8 In the case of a linear flux f (u) = au for some constant a. If U is an


u
entropy, the associated entropy flux is F = au and we have, as u satisfies u
t + a x = 0,
U
F
U
U
u
u
+
=
+a
= U 0 (t)(
+ a ) = 0.
t
x
t
x
t
x
So the entropy condition is always verified in this case.
This is a simple case, which also can appear for some of the waves in non linear
systems. In this case two different values are just propagated side by side without interfering. This case is called a contact discontinuity.
We will restrict in this lecture to the two simplest cases of scalar conservation laws,
namely the cases of a strictly convex flux and of a linear flux. These often occur in
applications, but not always. A more complex study is needed in other cases. We refer
to the book of Leveque [4] for details.

4.3

The Riemann problem

Let us here compute the exact entropy solution of the Riemann problem for a strictly
convex flux i.e. f 00 (u) > 0.
The Riemann problem for the 1D scalar non linear conservation laws reads
u f (u)
+
= 0,
t
x
34


u
u0 (x) = L
uR

if
if

x < 0,
x 0.

Case 1: uL < uR . We saw previously that in this case the shock does not satisfy the
Lax entropy condition and is not entropic. We thus look for a rarefaction wave which is
continuous across each of the characteristics issued from the left and right of 0. Between
these two characteristics, we look for an explicit solution of the form u(t, x) = v(x/t).
Plugging this function into our equation, we get

x 0 x
x 1 x
v ( ) + f 0 (v( )) v 0 ( ) = 0.
2
t
t
t t
t

Setting = x/t this becomes


(f 0 (v()) )v 0 () = 0.
So we get two kinds of solution, either v() = C (C constant), or v such that f 0 (v()) =
0. The constant solution is not possible as it would yield non entropic shocks. As f 00 > 0,
f 0 is strictly increasing and thus invertible. Hence the solution can be expressed as
v() = (f 0 )1 ().
Finally the entropy solution of the Riemann problem in this case is the rarefaction
wave defined by


uL
if x < f 0 (uL )t,
0 1
u(t, x) = (f ) (x/t) if f 0 (uL )t x < f 0 (uR )t,

uR
if x f 0 (uR )t.
Case 2: uL > uR . In this case we saw that the shock separating the two constant states
uL and uR and propagating at the speed defined by the Rankine-Hugoniot condition
s=

f (uL ) f (uR )
uL uR

is entropic and so is the unique physical solution. This is defined by



u
if x < st,
u(t, x) = L
uR if x st.

4.4

Numerical methods

4.4.1

The Godunov method

The idea of the Godunov method is that if one is looking for piecewise constant solutions
then as the solution is propagating at a finite speed which is known, the solution can be
computed exactly by solving a Riemann problem for each cell interface. If s = max |f 0 (u)|
denotes the fastest possible wave then the neighbouring Riemann problems do note
interact provided the time step t verifies st 21 x. In principle the solution at time
tn+1 could be computed by integrating over the solutions of all the Riemann problems at
this time to find the new constant value on each cell by averaging. However, this could

35

be quite complicated. The Godunov method can be expressed in a much simpler way
by integrating the starting equation in space and time over one cell and one time step:

Z tn+1 Z xi+1 
u f (u)
dx dt = 0.
+
t
x
tn
xi
R xi+1
1
Denoting by uni+ 1 = x
u(tn , x) dx the average of u over a cell at time tn = nt,
xi
2

this yields

x(un+1
i+ 21

uni+ 1 )
2

tn+1

(f (u(t, xi+1 )) f (u(t, xi ))) dt = 0.

+
tn

Now as u is assumed to be constant on each cell and equal to its average, the value of
the flux f (u(t, xi )) can be exactly computed by solving the Riemann problem between
the left value uni 1 and the right value uni+ 1 at the interface xi . Note that we only need
2

the solution of the Riemann problem exactly on the interface and this solution does not
depend on time as long as there is no interaction with the neighbouring problem on
the interface, which cannot occur if st x. Here because we are only interested in
the solution directly at the interface we gain a factor 12 as we do not care about other
interactions.
As we saw previously, assuming that the flux f (u) is strictly convex, the entropy
solution of the Riemann problem can be either a shock or a rarefaction wave, but u(t, xi )
will take a value different from uni 1 or uni+ 1 only if the rarefaction wave has an interaction
2

with x = xi which is the case only if f 0 (uL ) < 0 and f 0 (uR ) > 0. Let us clearly distinguish
the cases and express the solution of the local Riemann problem at the interface x = xi .
Recall that we consider the local Riemann problem with the initial condition (at time
tn )
n
u 1 if x < xi ,
i
u(tn , x) = n 2
ui+ 1 if x xi .
2

If uni 1 < uni+ 1 the entropy solution is a rarefaction wave so that


2
2

un 1
if f 0 (uni 1 ) > 0,

2
2
0 i
1
0 n
0 n
u(t, xi ) = (f ) (0) if f (ui 1 ) 0 f (ui+ 1 ),
2
2
un
0 (un ) < 0.
if
f

1
1
i+
i+
2

If uni 1 > uni+ 1 the entropy solution is a shock wave so that


2
2
n
u 1 if s > 0,
i
u(t, xi ) = n 2
ui+ 1 if s < 0,
2

(uR )
where s is the shock speed given by the Rankine-Hugoniot condition: s = f (uuLL)f
.
uR
n
n
n
Using the exact solution of the Riemann problem, the numerical flux gi = g(ui 1 , ui+ 1 )
2

is taken to be the flux associated to the exact solution of the Riemann problem, i.e.
f (u(t, xi )) for x = xi and t > tn . Noticing that for a strictly convex flux, the minimum
of f is reached at the points (f 0 )1 (0) and going through the different cases we find that


f (u) if ui 1 < ui+ 1 ,
u[u min
2
2
1 ,ui+ 1 ]
i

2
2
gin =
f (u) if ui 1 > ui+ 1 .
max

2
2
u[ui 1 ,ui+ 1 ]
2

36

With this, the flux at the interface can be computed and the Godunov scheme is
well defined (here in the case of a strictly convex flux). In the case of a linear flux, the
Godunov scheme amounts to the first order upwind scheme. For other smooth fluxes,
the Riemann problem can generally also be solved even though it is more complicated.

4.4.2

Approximate Riemann solvers

Computing the exact solution of the local Riemann problems is in principle the best
method. However it can be complicated or numerically expensive. Moreover, as there
are other sources of numerical errors anyway, it might be numerically as good or almost
as good to replace the solution of the exact Riemann problem by some approximation.
1. The Rusanov flux (also called Local Lax-Friedrichs flux). The idea behind this flux,
instead of approximating the exact Riemann solver, is to recall that the entropy
solution is the limit of viscous solutions and to take a centred flux to which some
viscosity (with the right sign) is added. This flux is given by
"
#
1
n
n
0
n
gi = g(ui 1 , ui+ 1 ) =
f (ui 1 ) + f (ui+ 1 )
|f (u)|(ui+ 1 ui 1 ) .
max
2
2
2
2
2
2
2
u[ui 1 ,ui+ 1 ]
2

Taking the viscosity parameter as the largest local wave speed guarantees the
stability of the scheme. On the other hand taking the largest local wave speed
(instead of globally) adds less viscosity in regions where the solution is smooth
and ui 1 is close to ui+ 1 .
2

2. The Roe flux (also called Murman or Murman-Roe in the scalar case): the idea here
is to linearise the flux f (u) around the cell interface and then use the upwind flux.
Assuming two constant values uL and uR on each side, this amounts to replacing
f (u) by a(uL , uR )u with a well chosen velocity a(uL , uR ) that will enable us to get
a flux which is close to the flux given by the solution of the exact Riemann problem.
Looking at the solution of the exact Riemann problem an approximation that looks
interesting for both rarefaction wave and shocks is to take the Rankine-Hugoniot
velocity
f (uL ) f (uR )
a(uL , uR ) =
,
uL uR
if there is a discontinuity i.e. uL 6= uR , and simply a(uL , uR ) = f 0 (uL ) = f 0 (uR ) if
uL = uR . Indeed if uL and uR are close, which will be the case in smooth regions
a(uL , uR ) defines a good approximation of both f 0 (uL ) and f 0 (uR ). Moreover if
both f 0 (uL ) and f 0 (uR ) have the same sign a(uL , uR ) will also have this sign, so
that the same value of the numerical flux will be obtained as in the exact Riemann
problem was solved. The only case where this does not yield the same solution as
the exact Riemann solver is the case of a rarefaction wave with f 0 (uL ) and f 0 (uR ) of
different signs. In this case, it is possible that the solution obtained by the scheme
is not the correct entropy solution, the non entropic shock being approximated
instead of the rarefaction wave. So a fix, known as entropy fix, is needed. For this
let us express the Murman-Roe flux as a centred flux with some added viscosity.

37

The Murman-Roe flux is an upwind flux for the linearised problem we just introduced. The case uL = uR being trivial, we consider only uL 6= uR . Then

f (uL ) if a(uL , uR ) > 0,
g(uL , uR ) =
f (uR ) if a(uL , uR ) 0.
If a(uL , uR ) =

f (uL )f (uR )
uL uR

1
f (uL ) =
2
and if a(uL , uR ) =

> 0 and uL 6= uR we have



f (uR ) f (uL )
f (uL ) + f (uR )
(uR uL ) ,
uR uL

f (uL )f (uR )
uL uR

1
f (uR ) =
2

< 0 and uL 6= uR we have



f (uR ) f (uL )
f (uL ) + f (uR ) +
(uR uL ) ,
uR uL

so that we can define the numerical flux in all cases by


g(uL , uR ) =

1
(f (uL ) + f (uR ) |a(uL , uR )|(uR uL )) .
2

Here we also see that the numerical viscosity vanishes when a(uL , uR ) 0 which
can happen close to the minimum of the convex function f . Then a non entropic
shock might be selected by the scheme. A simple fix, introduced by Harten, consists
in smoothing the graph of the absolute value close to 0 (see [4] for details). This
consists in replacing the absolute value in the formula defining the flux by


||
|| ,
() = 2
( + 2 )/(2) || < .
This ensures that ()  and that there is always some dissipation. This works
and yields the correct entropy solution provided  is well tuned to the problem at
hand.

4.4.3

Higher order methods

Higher order methods can be designed using the Finite Volume or the Discontinuous
Galerkin methodology. In both cases what needs to be defined is the numerical flux
constructed from the two values uL and uR on each side of the interface. The same
numerical fluxes as previously can be used in this case. However for the scheme to be
stable close to discontinuities a limiting procedure needs to be introduced so that the
scheme falls back to first order in the cells where a discontinuity is detected. There is a
vast literature on limiters that we shall not detail here.

4.4.4

Strong stability preserving (SSP) Runge-Kutta methods.

Fluxes are generally constructed so that the associated scheme has some stability properties when used with a first order explicit Euler time solver. When higher order methods
in x are used, it makes sense to use also higher order methods in time. A specific class
of Runge-Kutta methods has been developed to keep the strong stability properties of
the explicit Euler solver at each stage (See [5]).
38

4.5

Nonlinear systems of conservation laws

Going from the scalar case to systems in the non linear case, is similar to what is done
in the linear case. The hyperbolicity of the system is essential so that the system can
be locally diagonalised and the eigenvalues explicitely used in the definition of the flux.
The derivation of a Finite Volume or Discontinuous Galerkin scheme can be done
component by component and so reduces to the scalar case except for the definition
of the numerical flux which in general mixes the different components and needs to be
specific to the system at hand. We shall restrict in this lecture to the introduction of
two of the most used numerical fluxes, namely the Rusanov (or local Lax-Friedrichs)
flux and the Roe flux.

4.5.1

The Rusanov flux

As in the scalar case, the main idea here is to use a centred flux to which just enough
dissipation is added to ensure stability in all cases. In the scalar case the needed viscosity
was given by the largest local wave speed. A system of n components corresponds to
the superposition of n waves the local speed of each being given by the corresponding
eigenvalue. So taking the viscosity coefficient in the flux as the maximum over all
eigenvalues should do the job. This yields the Rusanov flux for systems, which is the
simplest stable flux. It is defined for a linear system of the form
U F(U)
+
= 0,
t
x
as

1
G(UL , UR ) =
2


F(UL ) + F(UR )

max

U [UL ,UR ]


|(F (U))|(UR UL ) ,
0

where maxU [UL ,UR ] |(F0 (U))| denotes the maximum modulus of the eigenvalues of the
Jacobian matrix F0 (U).

4.5.2

The Roe flux

Roes method consists in locally linearising the non linear flux with a well chosen procedure. The linearised matrix between two constant states UL and UR is denoted by
A(UL , UR ) and constructed such that the following properties are verified:
F(UR ) F(UL ) = A(UL , UR )(UR UL ).
A(U, U) = F0 (U).
A(UL , UR ) is diagonalisable, has real eigenvalues and a complete system of eigenvectors.
Such a matrix is not always easy to find, but there are procedures, described in [4] for
example, to construct them. Moreover classical Roe matrices are known for the most
usual systems [4].
Once the Roe matrix is defined, the flux can be computed by solving the corresponding linear Riemann problem that we treated in Chapter 3. Let us rewrite the formula,
so that we can also include the entropy fix, which as in the scalar case is needed for non
linear systems to make sure that the scheme always converges to the correct entropy
solution.
39

In the case of linear systems, the flux was defined as AU (t, 0) = A+ UL + A UR .


Defining the absolute value of a matrix as |A| = A+ A , the flux can also be expressed
as
1
AU (t, 0) = (AUL + AUR |A|(UR UL )) .
2
Using the properties of the Roe matrix in the non linear case, the same expression
will be used to define the Roe flux:
G(UL , UR ) =

1
(A(UL , UR )UL + A(UL , UR )UR |A(UL , UR )|(UR UL )) .
2

The same entropy fix consisting in replacing the absolute value by the function which
is bounded away from 0 can be applied in this formula.

40

Bibliography
[1] J.-P. Berrut, L.N. Trefethen: Barycentric Lagrange Interpolation, SIAM Rev. 46
(3), pp. 501517.
[2] G. Cohen, Higher-Order Numerical Methods for Transient Wave equation,
Springer-Verlag, 2001.
[3] E. Godlewski and P.A. Raviart: Numerical Approximation of Hyperbolic Systems
of Conservation Laws, Springer, 1996.
[4] R. J. Leveque: Finite Volume Methods for Hyperbolic Problems, Cambridge Texts
in Applied Mathematics, 2002.
[5] J. S. Hesthaven and T. Warburton: Nodal Discontinuous Galerkin methods,
Springer, 2008.
[6] D. A. Di Pietro and A. Ern: Mathematical Aspects of Discontinuous Galerkin
Methods, vol. 69 SMAI Mathematiques et Applications, Springer, 2012.
[7] Chi-Wang Shu: High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems. SIAM Rev. 51 (1), pp. 82126.

41

You might also like