Chapter 9

Finite Volume Discretization

9.1 The One-dimensional Setting: Basic Ideas

Consider the Cauchy problem of finding u = u(x, t), such that
ut + (f (u))x = 0, in R ⇥ (0, 1) (9.1)
u(x, 0) = u0 (x) x 2 R. (9.2)
The general idea behind finite volume schemes is to define a numerical method for volume integrals
(or volume averages), taken over a number of non-overlapping control volumes, whose union makes
up the computational domain. (Here: the real line.)
S To make this more precise, partition the real line into non-overlapping intervals, such that R =
1
I
i2Z i , where I i = [x 1
i 2 , x 1
i+ 2 ]. Denote the centroid of an interval x i := 2 (x i 2 + xi+ 12 ), and
1

integrate (9.1) over Ii to obtain

Z
ut + (f (u))x dx = 0. (9.3)
Ii

The volume of an interval is defined

Z xi+ 1
2
hi := dx = xi+ 12 xi 1
2
.
xi 1
2

In general, h is di↵erent in each interval. For simplicity we assume that hi = h is constant. Next we
introduce the volume average Z
1
ui (t) := u(x, t) dx,
h Ii
and write (9.3) as
Z
dui
h + (f (u))x dx =0, (Definition of volume average)
dt Ii
dui
, h + f (u(xi+ 12 , t)) f (u(xi 12 , t)) =0, (Fundamental theorem of calculus) (9.4)
dt
Note that no approximation has been used so far. We have simply derived an equation for the volume
averages, i.e., if u satisfies Eq. (9.1), then the volume averages satisfy Eq. (9.4) for all intervales Ii .
We now define a numerical scheme by introducing as the unknowns the volume averages ui . To
get a pointwise definition of the numerical solution, we could then use a step function
X
uh := ui i , (9.5)
i2Z

where (Ii ) is the characteristic function defined on the interval Ii :

⇢
1, x 2 Ii
i (x) = (9.6)
0, otherwise .

91
92 CHAPTER 9. FINITE VOLUME DISCRETIZATION

A Taylor series argument reveals that an approximation of the type (9.5) cannot be better than first-
order accurate in the point-wise sense. Alternatively, we could attempt a higher order polynomial
reconstruction from the volume averages. Of course, whether this makes sense, depends on how
well we approximate the volume averages.1 For now, it will suffice to work with the volume averages
directly. We may think of these volume averages as being arranged in a sequence of values u := (ui )i2Z
Formulating a scheme for the volume averages is certainly made simpler by the fact that we can
easily derive an exact equation for the volume averages (Eq. (9.4)), which we can use as a starting
point. Inspection of Eq. (9.4) reveals that we merely need some kind of method at the interval
interfaces to approximate the flux f (u(xi+ 12 , t)). For example, we may express the flux function there
in terms of the volume averages in adjacent intervals. To this end, define a function fˆ, such that

fˆi+ 12 = fˆ(ui , ui+1 ) ⇡ f (u(xi+ 12 , t)).

We shall give specific examples for such an approximation below. At this point, we require only that
the numerical flux be consistent, which means fˆ(u, u) = f (u). An approximation to Eq. (9.4) is then
written
dui
h + fˆi+ 12 fˆi 1 = 0, (t > 0, i 2 Z). (9.7)
dt 2

For now, initial conditions are set by simply starting with the exact volume averages of the initial
conditions:
Z
1 xi+ 12
ui (0) = u0 (x) dx, (i 2 Z). (9.8)
h x 1
i
2

We need not assume that the numerical flux is di↵erentiable, but merely that it is locally Lipschitz
continuous. It is now easy enough to give a fully discrete version of the scheme by simply discretizing
the time derivative.

Scheme 9.1 (1D Finite Volume Scheme) Let

Z
1 xi+ 12
u0i = u0 (x) dx, (i 2 Z). (9.9)
h x 1
i
2

For n = 0, 1, 2, . . . , let

un+1 uni 1 n
i
+ (fˆi+ 1 fˆin 1 ) = 0, (⌧ > 0, i 2 Z), (9.10)
⌧ h 2 2

where the superscript n indicates, as usual, evaluation at discrete time instance tn = n⌧ . In

particular, fˆi+
n ˆ n n
1 = f (ui , ui+1 ).
2

For the schemes we discuss here, R xit 1will become apparent that, instead of using exact integration,
the midpoint rule u0i = u0 (xi ) ⇡ h1 x i+12 u0 (x)dx may be used to set initial conditions. (The midpoint
i
2
rule is second order accurate. The schemes we discuss below, are only first order accurate.)

Example 9.1 Consider the linear advection equation

ut + aux = 0.
1 One has to carefully distinguish between the accuracy of u , as defined by (9.5), and the accuracy of the volume
h
averages ui . In this chapter, we restrict ourselves to first order accurate schemes, i.e. the scheme for the volume
averages ui is first order consistent (see below). If the volume averages are known to higher accuracy, then a pointwise
approximation of higher order can be generated using interpolation.
9.1. THE ONE-DIMENSIONAL SETTING: BASIC IDEAS 93

Setting f (u) := au, we can write the advection equation as a conservation law

ut + (f (u))x = 0.

Recall the upwind scheme (6.10). We may write that scheme in the form (9.10) if we set

a n |a| n
fˆi+
n
1 = (u + uni ) (u uni ). (9.11)
2 2 i+1 2 i+1
It should be noted that whether we view the upwind scheme for this problem as a finite di↵erence
scheme or as a finite volume scheme is a matter of interpretation. Also, recall that the upwind
scheme is only first order consistent, so clearly one can evaluate the initial conditions using the
midpoint rule.

Even if we treat a more general equation, the interpretation of the 1D finite volume method as a
finite di↵erence method is straightforward, as the next example shows.

Example 9.2 Consider the general conservation law

ut + (f (u))x = 0.

Approximate this using the finite volume nomenclature (9.10) with a flux function inspired by the
linear advection case (9.11),
1 ↵ n
fˆi+
n
1 = (f (uni+1 ) + f (uni )) (u uni ), (↵ > 0), (9.12)
2 2 2 i+1
where ↵ is now some ”di↵usion” coefficient, which remains to be determined. Note that the flux
is consistent, i.e., fˆ(u, u) = f (u). We can write, for n = 0, 1, . . . ,

un+1 uni 1 ↵ uni 1 2uni + uni+1

i
+ (f (uni+1 ) f (uni 1 )) = 0, (i 2 Z),
⌧ 2h 2 h
where in the finite di↵erence parlance we interpret the second term as a central di↵erence approx-
imation of (f (u))x , and the last term is of the kind

↵ u(x h, t) 2u(x, t) + u(x + h, t) h↵

= u(x, t)xx + O(h3 ).
2 h 2
We have used di↵usive terms before to stabilize central di↵erence approximations for hyperbolic
equations, and we may take the general form of the flux function (9.12) as a blueprint for the design
of finite volume schemes. (Of course, we need to subject this to some more rigorous analysis later.
In particular, the stabilization parameter ↵ needs to be determined.)
It is easy to see that the scheme with flux (9.12) is overall first order consistent for smooth
solutions, as the di↵usive term is already scaled with h. (The central di↵erence of the flux is second
order accurate for smooth enough f and u.)

So far the di↵erence between finite volume methods and finite di↵erence methods has appeared
to be largely a matter of interpretation, and we may ask why we should even bother with finite
volume schemes. The strength of the finite volume approach is seen in the multidimensional case.
The extension of the 1D method to the multidimensional case operating on rather general meshes
is straightforward in the finite-volume framework. The finite di↵erence formulation is not so easily
extended to the multidimensional case on general meshes2 .
2 The phrase ”general meshes” is important here. Extending the finite di↵erence method to simple two-dimensional
94 CHAPTER 9. FINITE VOLUME DISCRETIZATION

9.2 Multidimensional Extension

The advantage of the finite volume method is the easy extension to multiple dimensions. Consider
the d-dimensional scalar conservation law
d
X
ut + (fi (u))xi = 0, in Rd ⇥ (0, 1) (9.13)
i=1

For a proper discretization in multiple dimensions, we need a more general concept of a grid. To
avoid dealing with boundary conditions, we formally consider a partition of Rd into non-overlapping
polyhedral elements i . For purposes of illustration we often use d = 2, and assume the elements to
be triangles.
Consider the volume average Z
1
ui (t) := u(x, t) dx,
|i | i
where Z
|i | := dx
i

is the volume of the element, and integrate (9.13) over i . We obtain, using the divergence theorem,
Z
dui
|i | + f (u) · n ds = 0 . (9.14)
dt @i

Here we have defined f = (f1 , . . . , fd ). Again, it is relatively straightforward to derive an equation for
the volume averages in each domain i . So far, no approximation has been made.
For polyhedral mesh elements, the boundary of i consists of, say bi , planar faces. (For general
polyhedra, the number of boundary faces may vary from cell to cell. For triangles the boundary is
given by three straight line segments, called edges.) Let the faces comprising @i be denoted ek for
k = 1, . . . bi . We have
bi Z
dui X
|i | + f (u) · nk ds = 0,
dt ek
k=1
where nk is the outward pointing normal on face ek .
We re-write the boundary integral as a sum over the interfaces between i and neighboring ele-
ments,
dui XZ
|i | + f (u) · nij ds = 0,
dt Sij j2Ni

where Ni := {j : j is a neighbor of i }. See Fig. 9.1 for an illustration for the case of triangles.3 ) The
interface shared by elements i and j is denoted by Sij , while nij is the outward pointing normal (i.e.
pointing from i to j ). If an interface Sij between two elements is a complete face of both elements,
then Sij is simply one of the ek , defined above. This is the case in Fig. 9.1a, but not in Fig. 9.1b.
As in the one-dimensional case, we define a numerical scheme by considering the volume averages as
the unknowns, and take the numerical solution as a piecewise constant function (i.e. constant on each
element.) The values on the interfaces Sij are thus double-valued, which necessitates again a suitable
approximation of the flux. We may approximate the integral over the face with a midpoint quadrature
rule and express f (u) · nij ⇡ fˆ(ui , uj ; nij ) as a function of the volume average of neighboring elements.
We obtain the scheme
dui X
|i | + fˆ(ui , uj ; nij )|Sij | = 0 , i = 1, 2 . . . , (9.15)
dt
j2Ni

meshes in space that look like fig. 3.1 is not hard. However, for many problems of practical interest with complex
geometry meshes composed of, say, triangles are used. Here the finite volume framework is more convenient.
3 It is worth pointing out that while for the purposes of this section there is no need to specify in more detail exactly

what kind meshes we consider, in later sections, we shall see that some numerical schemes do not allow configurations
as in 9.1b. (Instead we will stipulate that if two triangles are neighbors, the shared interface must be an edge of both
triangles, as in 9.1a. Furthermore, when we carry out convergence analysis, we need to specify more clearly in what
way one may ”let the mesh size go to zero”.
9.2. MULTIDIMENSIONAL EXTENSION 95

i
i
j

j

(a) Triangles i and j are neighbors (b) Triangles i and j are neighbors

i

j

(c) Triangles i and j are not neigh-

bors

Figure 9.1: Illustration of di↵erent configurations of triangles with non-empty intersection, and non-
overlapping interior. A neighbor of i is a triangle that shares an edge, or a segment of an edge, with
i . If a triangle just shares an isolated point (a vertex) with i it is not a neighbor.

Here, |Sij | denotes the area of the interface Sij . This is the classical finite-volume scheme in semi-
discrete form. The numerical flux may be approximated by a variety of di↵erent methods. For
example, we may use
1 ↵ij
fˆ(ui , uj ; nij ) := (f (uj ) + f (ui )) · nij (uj ui ).
2 2
which is inspired by the flux we used as a blueprint in the one-dimensional case (cf. eq. (9.12)). It
will be shown below that this is indeed a reasonable choice.
Equation (9.15) is a system of nonlinear ordinary di↵erential equations (ODE). The simplest fully
discrete finite-volume scheme is obtained by using a one-step forward discretization in time, i.e.

dui un+1 uni

= i + O(⌧ ),
dt ⌧
where uni denotes the numerical approximation of the cell avarage at time tn = n⌧ , which leads to
⌧ X ˆ n n
un+1
i = uni f (ui , uj ; nij )|Sij |. (9.16)
|i |
j2Ni
96 CHAPTER 9. FINITE VOLUME DISCRETIZATION
Chapter 10

Stability Analysis for Finite Volume

Schemes

10.1 Stability in Maximum Norm

The properties of finite volumes schemes are essentially determined by the choice of numerical flux
function. We do not assume that the numerical flux functions are smooth. But we will assume that
they are at least locally Lipschitz continuous (i.e. Lipschitz on closed and bounded sets). In addition,
we note a few basic definitions that will be useful in the analysis below.

Definition 10.1 A numerical flux function fˆ is called conservative if it satisfies

fˆ(u, v; n) = fˆ(v, u; n) . (10.1)

Definition 10.1 ensures that flux contributions from interior edges cancel in the summation of
Eq. (9.15) over an arbitrary number of connected mesh elements. This is a discrete analog of the
property of the integral conservation laws that the rate of change of conserved quantities, for arbitrary
control volumes, be given by the flux across the boundary of the control volume1 .
The numerical flux was introduced as an approximation to f (u) · n at interfaces between two mesh
elements. Such an approximation, using the values from both adjacent elements, is necessary, because
the numerical solution is double-valued at these interfaces. It makes sense that fˆ should give back the
exact flux function f (u) · n if the solution is single-valued on the interface. This leads to the following
notion of consistency:

Definition 10.2 We call a numerical flux function fˆ consistent if it satisfies

fˆ(u, u; n) = f (u) · n (10.2)

Note that, in particular this holds if u is the exact solution. It is easily seen that scheme (9.15) is
consistent, if the numerical flux is consistent.
In the following we always assume that the numerical flux functions which we use are conservative
and consistent. These ingredients are essential. A third property is very helpful for the analysis of
finite-volume schemes, although it is not always verified in practice.

1 Recall the derivation of the compressible Euler equations in Chapter 8.1.

97
98 CHAPTER 10. STABILITY ANALYSIS FOR FINITE VOLUME SCHEMES

Definition 10.3 A numerical flux function fˆ(u, v; n) is called monotone if it is monotonically

increasing / decreasing in the first and second argument, respectively. In the case of di↵erentiability
this means
@ fˆ @ fˆ
0, 0. (10.3)
@u @v

Before going to the fully discrete scheme, we note the consequence of the monotonicity definition.
Begin by writing scheme (9.15) as

dui 1 X ˆ
= f (ui , uk ; nik )|Sik | (Definition of scheme)
dt |i |
k2Ni
1 X ⇣ˆ ⌘
= f (ui , uk ; nik ) f (ui ) · nik |Sik | (See Eq. (10.5))
|i |
k2Ni
1 X ⇣ˆ ⌘
= f (ui , uk ; nik ) fˆ(ui , ui ; nik ) |Sik | (Consistency) (10.4)
|i |
k2Ni

where we have used

! ✓Z ◆
X X
f (ui ) · nik |Sik | = f (ui ) · nik |Sik | = f (ui ) · nds = 0, (10.5)
k2Ni k2Ni @Ti

as we integrate over a closed surface. (Also recall that the boundary is composed of planar faces. This
explains the second-to-last equality in Eq. (10.5).) Note that if ui = uk , the corresponding term in
the sum on the right-hand side of (10.4) is zero. Assume now that ui 6= uk . Then

dui 1 X fˆ(ui , uk ; nik ) fˆ(ui , ui ; nik )

= (uk ui )|Sik |
dt |i | uk ui
k2Ni
1 X
= Cik (u)(uk ui ) , (10.6)
|i |
k2Ni

where we use u = (ui )i2N to emphasize the dependence of the coefficients on the solution. Finally,
the monotonicity property, Definition 10.3, ensures that

fˆ(ui , uk ; nik ) fˆ(ui , ui ; nik )

Cik = |Sik | 0 (10.7)
uk ui
| {z }
0

Since the numerical flux is Lipschitz continuous, if the numerical solution remains bounded, the
coefficients are bounded from above as well. 2
Of course, we are mainly interested in deriving a stability condition for the fully discrete case (9.16).
This can be written X
un+1 = Ceik (un )un , (10.8)
i k
k2Ni0

where
Ni0 = Ni [ {i},
2 Indeed, it is easy to see that the solution verifies a local maximum principle, i.e. has the property that local

maxima/minima cannot increase/decrease. The right-hand side of (10.6), noting that Cik 0, directly implies that,
if ui is a local maximum, then du dt
i
 0. So the local maximum cannot increase. Likewise, a local minimum cannot
decrease. This can be used to show the well-posedness of the ODE problem with bounded initial data. However, we
will be more interested in the fully discrete case.
10.1. STABILITY IN MAXIMUM NORM 99

eik (un ) =: C
and the coefficients C e n are given by
ik
(
⌧
Cn k 6= i
eik = |i | ik P
C n
(10.9)
1 |⌧i | l2Ni Ciln k=i

with Cik defined as in (10.7). Note that clearly

X
e n = 1,
C (10.10)
ik
k2Ni0

We have already noted in Section 4.3 for the linear 1D case that any consistent scheme written in the
form of (10.8) will satisfy (10.10). We have also seen in Section 4.3 that positivity of the coefficients
plays a major role in establishing stability in maximum norm. It turns out that this holds for the
nonlinear case as well.

Proposition 10.1 Consider the scheme given by Eq. (10.8), and assume that Eq. (10.10) holds.
Furthermore, assume that
max |u0i | < 1.
i

Then, one has

max |un+1
i |  max |uni | n = 0, 1, 2, . . . , (10.11)
i i

en
if and only if C 0, for all i, k.
ik

P en
Proof We show that k |C ik |  1, for all i is necessary and sufficient for (10.11) to hold. If (10.10)
holds, this condition is equivalent to C en 0.
ik
P eik |  1 implies (10.11). Indeed, for any i we have
Step 1: We show k |C

X
|un+1 |= eik
C n n
uk (Definition of scheme)
i
k2Ni0
X
 eik
|C n
| |unk | (Triangle inequality)
k2Ni0
X
 max |uni | eik
|C n
| (Trivial : |unk |  max |uni |, 8k)
i i
k2Ni0

P e n |  1 holds for all i, then (10.11) follows.

If k |C ik
P en
Step 2: Next, we show |Cik |  1 is necessary. We merely have to find one counterexample.
e n ) for k 2 N 0 , and un = 0 for
Consider the linear case f (u) = au. For fixed i, let unk = sign(C ik i k
0 n
k 62 Ni . Note that this implies maxi |ui | = 1. We have
X
|un+1 | = en |
|C
i ik
k2Ni0

P en
if |Cik | > 1 for any i, clearly maxi |un+1
i | > 1. ⌅

Proposition 10.1 is a rather general result, and applies to any scheme that can be written in the
e n are given
canonical form (10.8). In particular, for our finite-volume scheme, where the coefficients C ik
by (10.9), we can derive a time-step restriction that ensures positivity.
100 CHAPTER 10. STABILITY ANALYSIS FOR FINITE VOLUME SCHEMES

Proposition 10.2 Assume we are given a discrete finite volume scheme of the form (9.16) with
consistent and monotone fluxes. Assume that the time step restriction
⌧n X
1 Cik (un ) 0, 8i, (10.12)
|i |
k2Ni

holds for each n , where the coefficients Cik (un ) are defined in (10.7). Then the scheme is stable
in maximum norm, and assuming that the initial data is bounded, the following estimate holds:

inf u(x, 0)  uni  sup u(x, 0) , 8i, n . (10.13)

x2Rd x2Rd

Proof Using the constraint (10.12) in the fully discrete form, written as in (10.8), we have non-
negative coefficients (cf. eq. (10.9)). Recall that monotonicity of the numerical flux ensures that
the coefficients Cik are non-negative.) Stability then follows from Proposition 10.1. Upon applying
maxi |un+1
i |  maxi |uni | recursively up to the initial conditions, one immediately finds (10.13) by
noting that maxi2Z u0i  supx2R u(x, 0) and mini2Z u0i inf x2R u(x, 0) if either the volume average
or the midpoint integration rule is used to define u0i from u(x, 0). ⌅

Very few ingredients are thus needed in order to establish nonlinear stability of the finite volume
scheme on general grids. The most succinct summary of the machinery which we have just established
is that for consistent and monotone fluxes, we can always ensure stability of the discrete scheme (9.16)
by choosing the time step as in (10.12). There are, however, many open points: Firstly, due to the
nonlinearity of the equations, establishing convergence to the entropy weak solution and pertaining
convergence rates is much more difficult. It has to be stated that no general equivalence theorem
between consistency and stability on the one hand, and convergence on the other hand applies in the
nonlinear case. Secondly, one may attempt to design higher order consistent schemes. (The present
formulation is restricted to first order consistency.) This necessitates a careful interpolation of the
volume averages to produce, for example, a piecewise linear representation of the solution. We will
not consider this here.

10.2 Constructing Schemes

In this section we investigate how, in practice, we can design a stable finite volume scheme. From the
analysis in the preceding section it is clear that this can be done by constructing monotone fluxes.
Consider the scheme
dui 1 X ˆ
= f (ui , uk ; nik )|Sik | (10.14)
dt |i |
k2Ni
with the flux function
1 ↵ik
fˆ(ui , uk ; nik ) := (f (uk ) + f (ui )) · nik (uk ui ).
2 2
We claim the scheme is stable for
↵ik |aik |, (10.15)
where (
(f (uk ) f (ui ))·nik
u k ui 6 ui
uk =
aik = (10.16)
f 0 (ui ) · nik uk = ui
With this definition, the flux is consistent and conservative. Assume ui 6= uk , and write
1 ↵ik
fˆ(ui , uk ; nik ) = f (ui ) · nik + (f (uk ) f (ui )) · nik (uk ui )
2 2
1
= f (ui ) · nik + (aik ↵ik ) (uk ui ) (10.17)
2
10.3. FURTHER READING 101

Note that (f (uk ) f (ui )) · nik = aik (uk ui ) holds by definition (10.16). Alternatively, isolating
f (uk ), we may write the flux as
1 ↵ik
fˆ(ui , uk ; nik ) = f (uk ) · nik (f (uk ) f (ui )) · nik (uk ui )
2 2
1
= f (uk ) · nik (aik + ↵ik ) (uk ui ) (10.18)
2

It is now easy to see that fˆ is monotonically decreasing in its second argument, and monotonically
increasing in its first argument. Indeed, using consistency of the flux we have

fˆ(ui , uk ; nik ) fˆ(ui , ui ; nik ) 1

(10.17) ) = (aik ↵ik )  0
uk ui 2
ˆ ˆ
f (ui , uk ; nik ) f (uk , uk ; nik ) 1
(10.18) ) = (aik + ↵ik ) 0
ui uk 2

where the inequalities follow from (10.15). Hence, the flux is monotone in the sense of Definition 10.3.
We can write the scheme in the from (10.6), with non-negative coefficients. To make this explicit,
substitute (10.17) into (10.14) to obtain

dui 1 X 1
= (↵ik aik ) (uk ui )|Sik |.
dt |i | 2
k2Ni

(Note that the first term in (10.17) vanishes in the summation over a closed surface.) Clearly the
coefficients Cik = 1/2 (↵ik aik ) |Sik | are non-negative under the constraint (10.15). We can conclude
that, for ui a local minimum (respectively maximum) du dt
i
0 (respectively du
dt  0). The fully discrete
i

scheme can be made stable under the simple time step restriction as in Theorem 10.2.

10.3 Further Reading

Barth and Ohlberger [2] discuss many results covered here. (Theorems are generally not proved
in [2]. But it is nevertheless a good reference, because the most relevant results are summarized, and
references are given where proofs can be found.)
102 CHAPTER 10. STABILITY ANALYSIS FOR FINITE VOLUME SCHEMES
Chapter 11

Convergence

For multidimensional nonlinear systems of hyperbolic conservation laws, the well-posedness theory is
highly deficient. Naturally then, there is no complete convergence theory for numerical schemes either.
For scalar equations (which have received most of our attention so far), there is both a complete well-
posedness theory, and a convergence theory for finite-volume schemes. However, a full account of this
theory is beyond the scope of these introductory notes. We restrict ourselves to the one-dimensional
case, and give a useful theorem, which states that a convergent conservative scheme must necessarily
converge to a proper weak solution. This is the content of the Lax-Wendro↵ theorem.

11.1 The Lax-Wendro↵ Theorem

In the following we extend (as beore) the numerical solution to a function defined for all (x, t) 2
R ⇥ (0, 1) in a piecewise constant manner:
uh (x, t) := uni , for xi 1
2
< x < xi+ 12 , tn < t < tn+1

Theorem 11.1 (Lax-Wendro↵ ) Assume that the one-dimensional scheme (Scheme 9.1) with
forward di↵erence time stepping (eq. (9.10)), is consistent with the conservation law (9.1), and let
⌧ = O(h). Assume that there exists a constant C, such that

sup |uni | = sup |uh (x, t)|  C

i2Z (x,t)2R⇥(0,1)
n 0

uniformly in h, and that for some function u,

Z
lim |uh u|dx = 0 (11.1)
h!0 K

for all compact subsets K ⇢ R ⇥ (0, 1). Then u is a weak solution of (9.1), i.e. u satisfies
Z 1Z Z
uvt + f (u)vx dxdt + u0 (x)v(x, 0) dx = 0, 8v 2 C01 (R ⇥ [0, 1)) (11.2)
0 R R

Proof (Sketch) Choose an arbitrary test function v 2 C01 (R ⇥ [0, 1)), and set vin := v(xi , tn ).
Multiply (9.10) by h ⌧ vin , and sum over all i, n, to obtain
1
XX 1
XX ⇣ ⌘
h vin un+1
i uni + ⌧ vin fˆi+
n
1 fˆin 1 = 0.
2 2
i2Z n=0 i2Z n=0

Note that the infinite sums do not present a problem here, due to the compact support of the test
functions: The sums may be truncated to finite sums that cover the support of v. Rearranging the
103
104 CHAPTER 11. CONVERGENCE

sums, and multipying by ”-1”, we have

1
XX 1
XX X
h un+1
i vin+1 vin + ⌧ fˆi+
n
1
n
vi+1 vin + h u0i vi0 = 0.
2
i2Z n=0 i2Z n=0 i2Z

Define vh as a piecewise constant function in the same manner as uh , and define fˆh , such that

fˆh = fˆi+ 12 , for xi < x < xi+1 , tn < t < tn+1 .

(Note that, because the fluxes are associated with the interfaces, the piecewise constant flux function
has an o↵set of h/2 compared to the uh and vh .) Now, noting that for any of the piecewise constant
functions defined thus far, e.g. for uh , we have
X Z 1Z
h⌧ uni = uh dxdt,
i,n 0 R

and also noting the initial conditions

Z xi+ 1
1 2
u0i = u0 (x) dx,
h xi 1
2

we have
Z 1 Z ✓ ◆
vh (x, t + ⌧ ) vh (x, t)
uh (x, t + ⌧ ) dxdt
0 R ⌧
Z Z !
1
vh (x + h2 , t) vh (x h
2 , t)
+ fh (x, t) dxdt
0 R h
Z
+ u0 vh (x, 0) dx = 0. (11.3)
R

If uh converges to u as in (11.1), from the preceding equation (11.3) one can proof, using standard
theorems from integration theory, that the limit function u satisfies (11.2) in the limit h ! 0. (Note
that the terms in the brackets are consistent di↵erence approximations to vt and vx .) ⌅

11.2 Further Reading

The Lax-Wendro↵ Theorem is discussed and proved in Godlewski and Raviart [6]. Again, Barth and
Ohlberger [2] discuss many results related to finite volume schemes including convergence and further
extensions, such as higher order methods.
Chapter 12

Further Remarks on Finite-Volume

Schemes

12.1 Boundary Conditions

We have not so far discussed boundary conditions. For hyperbolic problems it is important to set
boundary conditions consistent with the propagation of information established by the characteristics.
Recall first the simplest case, the scalar advection equation, on a finite domain, for example,
ut + aux = 0, x 2 (0, 1), t > 0, a 2 R.
From the analysis in section 5.3 we recall that for a > 0 we have to set boundary conditions at the
”left” boundary x = 0, but we must not specify the solution at x = 1. (If a < 0, the situation is
opposite.) For linear systems (cf. section 7), the situation is similar. In fact, after a system, such as
eq. (7.1) has been transformed to characteristic variables (e.g. eq. 7.3), the problem decomposes into
independent advection equations, see (7.5). Consequently, we can set boundary conditions for each of
the characteristic variables, just like for the scalar linear advection equation.
For nonlinear scalar problems, i.e. (8.15), recall that the characteristics are still straight lines, and
thus one must still prescribe boundary conditions whenever characteristics enter the domain. There
may be situations, where boundary conditions must be prescribed at both boundaries, or no boundary
conditions may be prescribed at all. There may also be cases, where boundary data must be prescribed
on some portions on the boundary, but not on others.

Example 12.1 Consider the Burgers’ equation

✓ 2◆
u
ut + = 0, x 2 (0, 1), t > 0.
2 x
u(x, 0) = u0 (x), x 2 [0, 1]

Convince yourself that for the initial conditions

1
u0 (x) = x
2
all characteristics are outgoing, and we must not prescribe boundary conditions either at x = 0,
or at x = 1. Conversely, for the initial conditions
1
u0 (x) = x
2
all characteristics run towards the interior of the domain, and will eventually form a shock wave.
At least initially, for some t > 0, we must prescribe boundary conditions at both boundaries. What
happens for larger t then depends on what those boundary conditions are.

105
106 CHAPTER 12. FURTHER REMARKS ON FINITE-VOLUME SCHEMES

For nonlinear systems the picture is more complicated, and we will not discuss the issue here.

Implementation of the boundary conditions is straightforward, but di↵erent from the finite-di↵erence
schemes. The latter have been extensively discussed in previous chapters. For finite-volume schemes
we implement boundary conditions as follows. Consider again a scalar one-dimensional problem. Let
us assume we need to set boundary conditions at the left boundary, and boundary data is given as
u = (t) there. Let us consider a finite-volume scheme, and assume that the first cell has index i = 1.
Then we have for that cell
⌧ ⇣ˆ n n ⌘
un+1
1 = un1 f (u1 , u2 ) f ( (tn )) ,
h
meaning we evaluate the flux function with the boundary conditions instead of the numerical solution.

12.2 Nonlinear Systems of Conservation Laws

Recall first the linear system
ut + Aux = 0,
where A 2 R m⇥m
is a constant matrix, and recall the upwind scheme for linear systems (cf. (7.7) in
Chapter 7), which is obtained by applying the upwind scheme to each advection equation that results
from spectral decomposition of the linear system. Recall that such a decomposition is always possible
by definition of a linear hyperbolic system. We write this scheme as

⌧A n ⌧ |A| n
un+1 = uni u uni + ui 2uni + uni+1 .
i
2h i+1 1
2h 1

(See Chapter 7 for pertaining information.) Setting f (u) = Au, and defining a numerical flux

1 |A| n
f̂i+ 12 = f (uni+1 ) + f (uni ) ui+1 uni , (12.1)
2 2
we can write the scheme as
⌧ ⇣ ⌘
un+1
i = uni f̂i+ 12 f̂i 12 .
h
Thus, we note that the di↵erence scheme has a straightforward interpretation as a finite-volume
scheme. The upwind method is hardly the only way of defining a reasonable scheme for systems.
Instead of using |A| to scale the di↵usion terms we may consider a scalar di↵usion coefficient ↵, such
that
1 ↵ n
f̂i+ 12 = f (uni+1 ) + f (uni ) u uni
2 2 i+1
where we can define ↵ in such a way that the scheme has at least as much di↵usion as the upwind
scheme by setting ↵ = max1pm | p |, where p is the pth eigenvalue of the matrix A.
For a nonlinear system
ut + f (u)x = 0,
we consider as a general flux function
1 Bn n
f̂i+ 12 = f (uni+1 ) + f (uni ) ui+1 uni
2 2
where a choice similar to (12.1) leads to

B n = |A(uni+1 , uni )|,

and A(u, v) is a suitable linearization of the flux function, i.e. a consistent approximation to the
Jacobian, such that A(u, u) = fu . Concrete Examples will be discussed in the tutorial. Alternatively,
a scalar di↵usion coefficient may be written

B n = ↵(uni+1 , uni )I
12.3. FURTHER READING 107

where I is the identity matrix, and ↵(u, v) approximates the spectral radius (i.e. the eigenvalue with
maximum absolute value) of the flux Jacobian.
In multiple dimensions, a similar method may be used. Consider the conservation law
d
X
ut + (f (u))xi = 0,
i=1

and the corresponding finite volume scheme of Section 9.2,

⌧ X
un+1
i = uni f̂ (uni , unk ; nik )|Sik |.
|i |
k2Ni

We may consider the general flux function

1 B(u, v; n)
f̂ (u, v; n) = (f (u) · n + f (v) · n) (v u) ,
2 2
where now
B = |A(u, v; n)|
is a consistent linearization of the normal flux function, i.e. A(u, u; n) = fu · n. Alternatively we may
use B = ↵(u, v; n)I, where ↵ is an approximation of the spectral radius of fu · n.

12.3 Further Reading

The theory for systems is both challenging and incomplete. A good reference is [7]. This book covers
both theory and numerical methods.