8.2 Finite Difference, Finite Element and Finite Volume Methods For Partial Differential Equations

8.
2
FINITE DIFFERENCE, FINITE ELEMENT
AND FINITE VOLUME METHODS
FOR PARTIAL DIFFERENTIAL
EQUATIONS
Joaquim Peiró and Spencer Sherwin

Department of Aeronautics, Imperial College, London, UK
There are three important steps in the computational modelling of any

physical process: (i) problem definition, (ii) mathematical model, and
(iii) computer simulation.
The first natural step is to define an idealization of our problem of
interest in terms of a set of relevant quantities which we would like to mea-
sure. In defining this idealization we expect to obtain a well-posed problem,
this is one that has a unique solution for a given set of parameters. It might not
always be possible to guarantee the fidelity of the idealization since, in some
instances, the physical process is not totally understood. An example is the
complex environment within a nuclear reactor where obtaining measurements
is difficult.
The second step of the modeling process is to represent our idealization of
the physical reality by a mathematical model: the governing equations of the
problem. These are available for many physical phenomena. For example, in
fluid dynamics the Navier–Stokes equations are considered to be an accurate
representation of the fluid motion. Analogously, the equations of elasticity in
structural mechanics govern the deformation of a solid object due to applied
external forces. These are complex general equations that are very difficult to
solve both analytically and computationally. Therefore, we need to introduce
simplifying assumptions to reduce the complexity of the mathematical model
and make it amenable to either exact or numerical solution. For example, the
irrotational (without vorticity) flow of an incompressible fluid is accurately
represented by the Navier–Stokes equations but, if the effects of fluid viscos-
ity are small, then Laplace’s equation of potential flow is a far more efficient
description of the problem.
1
S. Yip (ed.),
Handbook of Materials Modeling. Volume I: Methods and Models, 1–32.
c 2005 Springer. Printed in the Netherlands.
2 J. Peiró and S. Sherwin
After the selection of an appropriate mathematical model, together with

suitable boundary and initial conditions, we can proceed to its solution. In this
chapter we will consider the numerical solution of mathematical problems
which are described by partial differential equations (PDEs). The three classical
choices for the numerical solution of PDEs are the finite difference method
(FDM), the finite element method (FEM) and the finite volume method (FVM).
The FDM is the oldest and is based upon the application of a local Taylor
expansion to approximate the differential equations. The FDM uses a topo-
logically square network of lines to construct the discretization of the PDE.
This is a potential bottleneck of the method when handling complex geome-
tries in multiple dimensions. This issue motivated the use of an integral form
of the PDEs and subsequently the development of the finite element and finite
volume techniques.
To provide a short introduction to these techniques we shall consider each
type of discretization as applied to one-dimensional PDEs. This will not allow
us to illustrate the geometric flexibility of the FEM and the FVM to their full
extent, but we will be able to demonstrate some of the similarities between the
methods and thereby highlight some of the relative advantages and disadvan-
tages of each approach. For a more detailed understanding of the approaches
we refer the reader to the section on suggested reading at the end of the chapter.
The section is structured as follows. We start by introducing the concept of
conservation laws and their differential representation as PDEs and the alter-
native integral forms. We next discusses the classification of partial differential
equations: elliptic, parabolic, and hyperbolic. This classification is important
since the type of PDE dictates the form of boundary and initial conditions
required for the problem to be well-posed. It also, permits in some cases, e.g.,
in hyperbolic equations, to identify suitable schemes to discretise the differen-
tial operators. The three types of discretisation: FDM, FEM and FVM are then
discussed and applied to different types of PDEs. We then end our overview by
discussing the numerical difficulties which can arise in the numerical solution
of the different types of PDEs using the FDM and provides an introduction to
the assessment of the stability of numerical schemes using a Fourier or Von
Neumann analysis.
Finally we note that, given the scientific background of the authors, the
presentation has a bias towards fluid dynamics. However, we stress that the
fundamental concepts presented in this chapter are generally applicable to
continuum mechanics, both solids and fluids.
1. Conservation Laws: Integral and Differential Forms

The governing equations of continuum mechanics representing the kine-
matic and mechanical behaviour of general bodies are commonly referred
Finite methods for partial differential equations 3
to as conservation laws. These are derived by invoking the conservation of

mass and energy and the momentum equation (Newton’s law). Whilst they are
equally applicable to solids and fluids, their differing behaviour is accounted
for through the use of a different constitutive equation.
The general principle behind the derivation of conservation laws is that the
rate of change of u(x, t) within a volume V plus the flux of u through the
boundary A is equal to the rate of production of u denoted by S(u, x, t). This
can be written as

∂
u(x, t) dV + f(u) · n dA − S(u, x, t) dV = 0 (1)
∂t
V A V
which is referred to as the integral form of the conservation law. For a fixed
(independent of t) volume and, under suitable conditions of smoothness of the
intervening quantities, we can apply Gauss’ theorem

∇ · f dV = f · n dA
V A
to obtain

∂u
+ ∇ · f (u) − S dV = 0. (2)
∂t
V
For the integral expression to be zero for any volume V , the integrand must be
zero. This results in the strong or differential form of the equation
∂u
+ ∇ · f (u) − S = 0. (3)
∂t
An alternative integral form can be obtained by the method of weighted
residuals. Multiplying Eq. (3) by a weight function w(x) and integrating over
the volume V we obtain

∂u
+ ∇ · f (u) − S w(x) dV = 0. (4)
∂t
V
If Eq. (4) is satisfied for any weight function w(x), then Eq. (4) is equivalent
to the differential form (3). The smoothness requirements on f can be relaxed
by applying the Gauss’ theorem to Eq. (4) to obtain

∂u
− S w(x) − f (u) · ∇w(x) dV + f · n w(x) dA = 0.
∂t
V A
(5)
This is known as the weak form of the conservation law.

Although the above formulation is more commonly used in fluid mechan-

ics, similar formulations are also applied in structural mechanics. For instance,
the well-known principle of virtual work for the static equilibrium of a body
[1], is given by

δW = (∇ + f ) · δ dV = 0
V
where δW denotes the virtual work done by an arbitrary virtual velocity δ ,

is the stress tensor and f denotes the body force. The similarity with the
method of weighted residuals (4) is evident.
2. Model Equations and their Classification

In the following we will restrict ourselves to the analysis of one-dimensional
conservation laws representing the transport of a scalar variable u(x, t) defined
in the domain = {x, t : 0 ≤ x ≤ 1, 0 ≤ t ≤ T }. The convection–diffusion-
reaction equation is given by

∂u ∂ ∂u
L(u) = + au − b −r u =s (6)
∂t ∂x ∂x
together with appropriate boundary conditions at x = 0 and 1 to make the prob-
lem well-posed. In the above equation L(u) simply represents a linear differen-
tial operator. This equation can be recast in the form (3) with f (u) = au − ∂u/∂ x
and S(u) = s + r u. It is linear if the coefficient a, b, r and s are functions of x
and t, and non-linear if any of them depends on the solution, u.
In what follows, we will use for convenience the convention that the pres-
ence of a subscript x or t under a expression indicates a derivative or partial
derivative with respect to this variable, for example
du ∂u ∂ 2u
u x (x) = (x); u t (x, t) = (x, t); u x x (x, t) = (x, t).
dx ∂t ∂x2
Using this notation, Eq. (6) is re-written as
u t + (au − bu x )x − r u = s.
2.1. Elliptic Equations

The steady-state solution of Eq. (6) when advection and source terms are
neglected, i.e., a=0 and s =0, is a function of x only and satisfies the Helmholtz
equation
(bu x )x + r u = 0. (7)
This equation is elliptic and its solution depends on two families of integration
constants that are fixed by prescribing boundary conditions at the ends of the
domain. One can either prescribe Dirichlet boundary conditions at both ends,
e.g., u(0) = α0 and u(1) = α1 , or substitute one of them (or both if r =/ 0) by a
Neumann boundary condition, e.g., u x (0) = g. Here α0 , α1 and g are known
constant values. We note that if we introduce a perturbation into a Dirichlet
boundary condition, e.g., u(0) = α0 + , we will observe an instantaneous
modification to the solution throughout the domain. This is indicative of the
elliptic nature of the problem.
2.2. Parabolic Equations

Taking a = 0, r = 0 and s = 0 in our model, Eq. (6) leads to the heat or
diffusion equation
u t − (b u x )x = 0, (8)
which is parabolic. In addition to appropriate boundary conditions of the form

used for elliptic equations, we also require an initial condition at t = 0 of the
form u(x, 0) = u 0 (x) where u 0 is a given function.
If b is constant, this equation admits solutions of the form u(x, t) = Aeβt
sin kx if β + k 2 b = 0. A notable feature of the solution is that it decays when
b is positive as the exponent β < 0. The rate of decay is a function of b. The
more diffusive the equation (i.e., larger b) the faster the decay of the solution
is. In general the solution can be made up of many sine waves of different
frequencies, i.e., a Fourier expansion of the form

u(x, t) = Am eβm t sin km x,
m
where Am and km represent the amplitude and the frequency of a Fourier mode,
respectively. The decay of the solution depends on the Fourier contents of the
initial data since βm = −km2 b. High frequencies decay at a faster rate than the
low frequencies which physically means that the solution is being smoothed.
This is illustrated in Fig. 1 which shows the time evolution of u(x, t) for
an initial condition u 0 (x) = 20x for 0 ≤ x ≤ 1/2 and u 0 (x) = 20(1 − x) for
1/2 ≤ x ≤ 1. The solution shows a rapid smoothing of the slope disconti-
nuity of the initial condition at x = 1/2. The presence of a positive diffusion
(b > 0) physically results in a smoothing of the solution which stabilizes it. On
the other hand, negative diffusion (b < 0) is de-stabilizing but most physical
problems have positive diffusion.
11
u(x)
10
t0
9
t T
8
t 2T
7
t 3T
6
t 4T
5
t 5T
4
t 6T
3
2
1
0
0.0 0.5 1.0
x
Figure 1. Rate of decay of the solution to the diffusion equation.
2.3. Hyperbolic Equations

A classic example of hyperbolic equation is the linear advection equation
u t + a u x = 0, (9)
where a represents a constant velocity. The above equation is also clearly

equivalent to Eq. (6) with b = r = s = 0. This hyperbolic equation also re-
quires an initial condition, u(x, 0) = u 0 (x). The question of what boundary
conditions are appropriate for this equation can be more easily be answered
after considering its solution. It is easy to verify by substitution in (9) that the
solution is given by u(x, t) = u 0 (x − at). This describes the propagation of
the quantity u(x, t) moving with speed “a” in the x-direction as depicted in
Fig. 2. The solution is constant along the characteristic line x − at = c with
u(x, t) = u 0 (c).
From the knowledge of the solution, we can appreciate that for a > 0 a
boundary condition should be prescribed at x = 0, (e.g., u(0) = α0 ) where in-
formation is being fed into the solution domain. The value of the solution at
x = 1 is determined by the initial conditions or the boundary condition at x = 0
and cannot, therefore, be prescribed. This simple argument shows that, in a hy-
perbolic problem, the selection of appropriate conditions at a boundary point
depends on the solution at that point. If the velocity is negative, the previous
treatment of the boundary conditions is reversed.
u (x,t ) Characteristic
x at c
t
u (x,0 ) x
Figure 2. Solution of the linear advection equation.
The propagation velocity can also be a function of space, i.e., a = a(x) or

even the same as the quantity being propagated, i.e., a = u(x, t). The choice
a = u(x, t) leads to the non-linear inviscid Burgers’ equation
u t + u u x = 0. (10)
An analogous analysis to that used for the advection equation shows that
u(x, t) is constant if we are moving with a local velocity also given by u(x, t).
This means that some regions of the solution advance faster than other re-
gions leading to the formation of sharp gradients. This is illustrated in Fig. 3.
The initial velocity is represented by a triangular “zig-zag” wave. Peaks and
troughs in the solution will advance, in opposite directions, with maximum
speed. This will eventually lead to an overlap as depicted by the dotted line
in Fig. 3. This results in a non-uniqueness of the solution which is obviously
non-physical and to resolve this problem we must allow for the formation and
propagation of discontinuities when two characteristics intersect (see Ref. [2]
for further details).
3. Numerical Schemes
There are many situations where obtaining an exact solution of a PDE is
not possible and we have to resort to approximations in which the infinite set
of values in the continuous solution is represented by a finite set of values
referred to as the discrete solution.
For simplicity we consider first the case of a function of one variable u(x).
Given a set of points xi ; i = 1, . . . , N in the domain of definition of u(x), as
t3
t2
u
u>0 t1
x
t 0
u 0
Figure 3. Formation of discontinuities in the Burgers’ equation.
ui
ui 1
ui 1
x
x1 xi 1 xi xi 1 xn
Ωi
xi 12 xi 12
Figure 4. Discretization of the domain.
shown in Fig. 4, the numerical solution that we are seeking is represented by a

discrete set of function values {u 1 , . . . , u N } that approximate u at these points,
i.e., u i ≈ u(xi ); i = 1, . . . , N .
In what follows, and unless otherwise stated, we will assume that the points
are equally spaced along the domain with a constant distance x = xi+1 − xi ;
i = 1, . . . , N − 1. This way we will write u i+1 ≈ u(xi+1 ) = u(xi + x). This
partition of the domain into smaller subdomains is referred to as a mesh or
grid.
3.1. The Finite Difference Method (FDM)

This method is used to obtain numerical approximations of PDEs written
in the strong form (3). The derivative of u(x) with respect to x can be defined
as
u(xi + x) − u(xi )
u x |i = u x (xi ) = lim
x→0 x
u(xi ) − u(xi − x)
= lim (11)
x→0 x
u(xi + x) − u(xi − x)
= lim .
x→0 2x
All these expressions are mathematically equivalent, i.e., the approximation
converges to the derivative as x → 0. If x is small but finite, the various
terms in Eq. (11) can be used to obtain approximations of the derivate u x of
the form
u i+1 − u i
u x |i ≈ (12)
x
u i − u i−1
u x |i ≈ (13)
x
u i+1 − u i−1
u x |i ≈ . (14)
2x
The expressions (12)–(14) are referred to as forward, backward and centred
finite difference approximations of u x |i , respectively. Obviously these approx-
imations of the derivative are different.
3.1.1. Errors in the FDM
The analysis of these approximations is performed by using Taylor expan-

sions around the point xi . For instance an approximation to u i+1 using n + 1
terms of a Taylor expansion around xi is given by

x 2 dn u x n
u i+1 = u i + u x |i x + u x x |i + ··· +
2 dx n i n!
dn+1 u ∗ x n+1
+ (x ) . (15)
dx n+1 (n + 1)!
The underlined term is called the remainder with xi ≤ x ∗ ≤ xi+1 , and repre-
sents the error in the approximation if only the first n terms in the expansion
are kept. Although the expression (15) is exact, the position x ∗ is unknown.
To illustrate how this can be used to analyse finite difference approxima-

tions, consider the case of the forward difference approximation (12) and use
the expansion (15) with n = 1 (two terms) to obtain
u i+1 − u i x
= u x |i + u x x (x ∗ ). (16)
x 2
We can now write the approximation of the derivative as
u i+1 − u i
u x |i = + T (17)
x
where T is given by
x
T = − u x x (x ∗ ). (18)
2
The term T is referred to as the truncation error and is defined as the
difference between the exact value and its numerical approximation. This term
depends on x but also on u and its derivatives. For instance, if u(x) is a linear
function then the finite difference approximation is exact and T = 0 since the
second derivative is zero in (18).
The order of a finite difference approximation is defined as the power p
such that limx→0 (T /x p ) = γ /=0, where γ is a finite value. This is often
written as T = O(x p ). For instance, for the forward difference approxima-
tion (12), we have T = O(x) and it is said to be first-order accurate ( p = 1).
If we apply this method to the backward and centred finite difference
approximations (13) and (14), respectively, we find that, for constant x, their
errors are
u i − u i−1 x
u x |i = + u x x (x ∗ ) ⇒ T = O(x) (19)
x 2
u i+1 − u i−1 x 2
u x |i = − u x x x (x ) ⇒ T = O(x 2 ) (20)
2x 12
with xi−1 ≤ x ∗ ≤ xi and xi−1 ≤ x ≤ xi+1 for Eqs. (19) and (20), respectively.
This analysis is confirmed by the numerical results presented in Fig. 5 that
displays, in logarithmic axes, the exact and truncation errors against x for the
backward and the centred finite differences. Their respective truncation errors
T are given by (19) and (20) calculated here, for lack of a better value, with
x ∗ = x = xi . The exact error is calculated as the difference between the exact
value of the derivative and its finite difference approximation.
The slope of the lines are consistent with the order of the truncation error,
i.e., 1:1 for the backward difference and 1:2 for the centred difference. The dis-
crepancies between the exact and the numerical results for the smallest values
of x are due to the use of finite precision computer arithmetic or round-off
error. This issue and its implications are discussed in more detail in numerical
analysis textbooks as in Ref. [3].
1e 00
Backward FD Total Error
Backward FD Truncation Error
Centred FD Total Error
1e 02 Centred FD Truncation Error
1e 04
1
1e 06
1
ε 2
1e 08
1
1e 10
1e 12
1e 14
1e 00 1e 02 1e 04 1e 06 1e 08 1e 10 1e 12
∆x
Figure 5. Truncation and rounding errors in the finite difference approximation of derivatives.
3.1.2. Derivation of approximations using Taylor expansions
The procedure described in the previous section can be easily transformed

into a general method for deriving finite difference schemes. In general, we can
obtain approximations to higher order derivatives by selecting an appropriate
number of interpolation points that permits us to eliminate the highest term
of the truncation error from the Taylor expansions. We will illustrate this with
some examples. A more general description of this derivation can be found in
Hirsch (1988).
A second-order accurate finite difference approximation of the derivative
at xi can be derived by considering the values of u at three points: xi−1 , xi and
xi+1 . The approximation is constructed as a weighted average of these values
{u i−1 , u i , u i+1 } such as
αu i+1 + βu i + γ u i−1
u x |i ≈ . (21)
x
Using Taylor expansions around xi we can write
x 2 x 3
u i+1 = u i + x u x |i + u x x |i + u x x x |i + · · · (22)
2 6
x 2 x 3
u i−1 = u i − x u x |i + u x x |i − u x x x |i + · · · (23)
2 6
Putting (22), (23) into (21) we get
αu i+1 + βu i + γ u i−1 1
= (α + β + γ ) u i + (α − γ ) u x |i
x x
x x 2
+ (α + γ ) u x x |i + (α − γ ) u x x x |i
2 6
x 3
+ (α + γ ) u x x x x |i + O(x 4 ) (24)
12
We require three independent conditions to calculate the three unknowns α,

β and γ . To determine these we impose that the expression (24) is consistent
with increasing orders of accuracy. If the solution is constant, the left-hand side
of (24) should be zero. This requires the coefficient of (1/x)u i to be zero and
therefore α+β +γ = 0. If the solution is linear, we must have α−γ =1 to match
u x |i . Finally whilst the first two conditions are necessary for consistency of
the approximation in this case we are free to choose the third condition. We
can therefore select the coefficient of (x/2) u x x |i to be zero to improve the
accuracy, which means α + γ = 0.
Solving these three equations we find the values α = 1/2, β = 0 and γ =
−(1/2) and recover the second-order accurate centred formula
u i+1 − u i−1
u x |i = + O(x 2 ).
2x
Other approximations can be obtained by selecting a different set of points,

for instance, we could have also chosen three points on the side of xi , e.g.,
u i , u i−1 , u i−2 . The corresponding approximation is known as a one-sided for-
mula. This is sometimes useful to impose boundary conditions on u x at the
ends of the mesh.
3.1.3. Higher-order derivatives
In general, we can derive an approximation of the second derivative using

the Taylor expansion
αu i+1 + βu i + γ u i−1 1 1
= (α + β + γ ) 2 u i + (α − γ ) u x |i
x 2 x x
1 x
+ (α + γ ) u x x |i + (α − γ ) u x x x |i
2 6
x 2
+ (α + γ ) u x x x x |i + O(x 4 ). (25)
12
Using similar arguments to those of the previous section we impose


α + β + γ = 0
α−γ =0 ⇒ α = γ = 1, β = −2. (26)

α+γ =2
The first and second conditions require that there are no u or u x terms on the
right-hand side of Eq. (25) whilst the third conditon ensures that the right-
hand side approximates the left-hand side as x tens to zero. The solution of
Eq. (26) lead us to the second-order centred approximation
u i+1 − 2u i + u i−1
u x x |i = + O(x 2 ). (27)
x 2
The last term in the Taylor expansion (α − γ )xu x x x |i /6 has the same coeffi-
cient as the u x terms and cancels out to make the approximation second-order
accurate. This cancellation does not occur if the points in the mesh are not
equally spaced. The derivation of a general three point finite difference ap-
proximation with unevenly spaced points can also be obtained through Taylor
series. We leave this as an exercise for the reader and proceed in the next
section to derive a general form using an alternative method.
3.1.4. Finite differences through polynomial interpolation
In this section we seek to approximate the values of u(x) and its derivatives
by a polynomial P(x) at a given point xi . As way of an example we will
derive similar expressions to the centred differences presented previously by
considering an approximation involving the set of points {xi−1 , xi , xi+1 } and
the corresponding values {u i−1 , u i , u i+1 }. The polynomial of minimum degree
that satisfies P(xi−1 ) = u i−1 , P(xi ) = u i and P(xi+1 ) = u i+1 is the quadratic
Lagrange polynomial
(x − xi )(x − xi+1 ) (x − xi−1 )(x − xi+1 )
P(x) = u i−1 + ui
(xi−1 − xi )(xi−1 − xi+1 ) (xi − xi−1 )(xi − xi+1 )
(x − xi−1 )(x − xi )
+ u i+1 . (28)
(xi+1 − xi−1 )(xi+1 − xi )
We can now obtain an approximation of the derivative, u x |i ≈ Px (xi ) as
(xi − xi+1 ) (xi − xi−1 ) + (xi − xi+1 )
Px (xi ) = u i−1 + ui
(xi − xi−1 )
+ u i+1 . (29)
(xi+1 − xi−1 )(xi+1 − xi )
If we take xi − xi−1 = xi+1 − xi = x, we recover the second-order accu-
rate finite difference approximation (14) which is consistent with a quadratic
interpolation. Similarly, for the second derivative we have

2u i−1 2u i
Px x (xi ) = +
2u i+1
+ (30)
(xi+1 − xi−1 )(xi+1 − xi )
and, again, this approximation leads to the second-order centred finite differ-
ence (27) for a constant x.
This result is general and the approximation via finite differences can be
interpreted as a form of Lagrangian polynomial interpolation. The order of the
interpolated polynomial is also the order of accuracy of the finite diference
approximation using the same set of points. This is also consistent with the
interpretation of a Taylor expansion as an interpolating polynomial.
3.1.5. Finite difference solution of PDEs
We consider the FDM approximation to the solution of the elliptic equation

u x x = s(x) in the region = {x : 0 ≤ x ≤ 1}. Discretizing the region using N
points with constant mesh spacing x = (1/N − 1) or xi = (i − 1/N − 1), we
consider two cases with different sets of boundary conditions:
1. u(0) = α1 and u(1) = α2 , and
2. u(0) = α1 and u x (1) = g.
In both cases we adopt a centred finite approximation in the interior points
of the form
u i+1 − 2u i + u i−1
= si ; i = 2, . . . , N − 1. (31)
x 2
The treatment of the first case is straightforward as the boundary conditions
are easily specified as u 1 = α1 and u N = α2 . These two conditions together with
the N − 2 equations (31) result in the linear system of N equations with N
unknowns represented by
    
1 0 ... 0 u1 α1
 1 −2 1 0 ... 0  u2   x 2 s2 
    
 0 1 −2 1 0 ...    x 2 s3 
 0  u3   
 .. .. ..  ..   .. 
 . . .  = .
  .   . 
    


0 ... 0 1 −2 1 0   u
  N −2
  x 2 s N −2
 


 0 ... 0 1 −2 1  u N −1
   x 2 s N −1 
0 ... 0 1 uN α2
This matrix system can be written in abridged form as Au = s. The matrix

A is non-singular and admits a unique solution u. This is the case for most
discretization of well-posed elliptic equations.
In the second case the boundary condition u(0) = α1 is treated in the same
way by setting u 1 = α1 . The treatment of the Neumann boundary condition
u x (1) = g requires a more careful consideration. One possibility is to use a
one-sided approximation of u x (1) to obtain
u N − u N −1
u x (1) ≈ = g. (32)
x
This expression is only first-order accurate and thus inconsistent with the
approximation used at the interior points. Given that the PDE is elliptic, this
error could potentially reduce the global accuracy of the solution. The alterna-
tive is to use a second-order centred approximation
u N +1 − u N −1
u x (1) ≈ = g. (33)
x
Here the value u N +1 is not available since it is not part of our discrete set of
values but we could use the finite difference approximation at x N given by
u N +1 − 2u N + u N −1
= sN
x 2
and include the Neumann boundary condition (33) to obtain
1
u N − u N −1 = (gx − s N x 2 ). (34)
2
It is easy to verify that the introduction of either of the Neumann boundary
conditions (32) or (34) leads to non-symmetric matrices.
3.2. Time Integration
In this section we address the problem of solving time-dependent PDEs

in which the solution is a function of space and time u(x, t). Consider for
instance the heat equation
u t − bu x x = s(x) in = {x, t : 0 ≤ x ≤ 1, 0 ≤ t ≤ T }
with an initial condition u(x, 0) = u 0 (x) and time-dependent boundary condi-

tions u(0, t) = α1 (t) and u(1, t) = α2 (t), where α1 and α2 are known
functions of t. Assume, as before, a mesh or spatial discretization of the

domain {x1 , . . . , x N }.
3.2.1. Method of lines
In this technique we assign to our mesh a set of values that are functions
of time u i (t) = u(xi , t); i = 1, . . . , N . Applying a centred discretization to
the spatial derivative of u leads to a system of ordinary differential equations
(ODEs) in the variable t given by
du i b
= 2 {u i−1 (t) − 2u i (t) + u i+1 (t)} + si ; i = 2, . . . , N − 1
dt x
with u 1 = α1 (t) and u N = α2 (t). This can be written as
 
bα1 (t)
     s2 +
u2 −2 1 u2 
 x 2 

 u3   1 −2 1     
d   b    u3   s3 
 ..   .. .. ..   .   .. 
 . =  . . .   ..  +  . 

dt   x 2     
u N −2   1 −2 1   u N −2  
 s N −2 

u N −1 1 −2 u  bα2 (t) 
N −1
s N −1 +
x 2
or in matrix form as
du
(t) = A u(t) + s(t). (35)
dt
Equation (35) is referred to as the semi-discrete form or the method of lines.
This system can be solved by any method for the integration of initial-value
problems [3]. The numerical stability of time integration schemes depends on
the eigenvalues of the matrix A which results from the space discretization.
For this example, the eigenvalues vary between 0 and −(4α/x 2 ) and this
could make the system very stiff, i.e., with large differences in eigenvalues, as
x → 0.
3.2.2. Finite differences in time
The method of finite differences can be applied to time-dependent prob-

lems by considering an independent discretization of the solution u(x, t) in
space and time. In addition to the spatial discretization {x1 , . . . , x N }, the dis-
cretization in time is represented by a sequence of times t 0 = 0 < · · · < t n <
· · · < T . For simplicity we will assume constant intervals x and t in space
and time, respectively. The discrete solution at a point will be represented by
u in ≈ u(xi , t n ) and the finite difference approximation of the time derivative

follows the procedures previously described. For example, the forward differ-
ence in time is given by
u(x, t n+1 ) − u(x, t n )
u t (x, t n ) ≈
t
and the backward difference in time is
u(x, t n+1 ) − u(x, t n )
u t (x, t n+1 ) ≈
t
both of which are first-order accurate, i.e., T = O(t).
Returning to the heat equation u t − bu x x = 0 and using a centred approx-
imation in space, different schemes can be devised depending on the time at
which the equation is discretized. For instance, the use of forward differences
in time leads to
u in+1 − u in b n
= u i−1 − 2u in + u i+1
n
. (36)
t x 2
This scheme is explicit as the values of the solution at time t n+1 are obtained
directly from the corresponding (known) values at time t n . If we use backward
differences in time, the resulting scheme is
u in+1 − u in b n+1
= u − 2u n+1
+ u n+1
. (37)
t x 2 i−1 i i+1
Here to obtain the values at t n+1 we must solve a tri-diagonal system of equa-
tions. This type of schemes are referred to as implicit schemes.
The higher cost of the implicit schemes is compensated by a greater numer-
ical stability with respect to the explicit schemes which are stable in general
only for some combinations of x and t.
3.3. Discretizations Based on the Integral Form
The FDM uses the strong or differential form of the governing equations.
In the following, we introduce two alternative methods that use their integral
form counterparts: the finite element and the finite volume methods. The use
of integral formulations is advantageous as it provides a more natural treat-
ment of Neumann boundary conditions as well as that of discontinuous source
terms due to their reduced requirements on the regularity or smoothness of the
solution. Moreover, they are better suited than the FDM to deal with complex
geometries in multi-dimensional problems as the integral formulations do not
rely in any special mesh structure.
These methods use the integral form of the equation as the starting point
of the discretization process. For example, if the strong form of the PDE is
L(u) = s, the integral from is given by
1 1
L(u)w(x) dx = sw(x) dx (38)
0 0
where the choice of the weight function w(x) defines the type of scheme.
3.3.1. The finite element method (FEM)
Here we discretize the region of interest = {x : 0 ≤ x ≤ 1} into N − 1

subdomains or elements i = {x : xi−1 ≤ x ≤ xi } and assume that the approx-
imate solution is represented by

N
u δ (x, t) = u i (t)Ni (x)
i=1
where the set of functions Ni (x) is known as the expansion basis. Its support
is defined as the set of points where Ni (x) =/ 0. If the support of Ni (x) is the
whole interval, the method is called a spectral method. In the following we will
use expansion bases with compact support which are piecewise continuous
polynomials within each element as shown in Fig. 6.
The global shape functions Ni (x) can be split within an element into two
local contributions of the form shown in Fig. 7. These individual functions are
referred to as the shape functions or trial functions.
3.3.2. Galerkin FEM
In the Galerkin FEM method we set the weight function w(x) in Eq. (38)
to be the same as the basis function Ni (x), i.e., w(x) = Ni (x).
Consider again the elliptic equation L(u) = u x x = s(x) in the region with
boundary conditions u(0) = α and u x (1) = g. Equation (38) becomes
1 1
w(x)u x x dx = w(x)s(x) dx.
0 0
At this stage, it is convenient to integrate the left-hand side by parts to get the
weak form
1 1
− wx u x dx + w(1) u x (1) − w(0) u x (0) = w(x) s(x) dx. (39)
0 0
ui 1 ui
u1 ui 1
Ωi uN
x
x1 xi 1 xi xi 1 xN
u1 x 1
x
..
.
Ni (x)
ui x 1
x
..
.
uN x 1
x
N
Figure 6. A piecewise linear approximation u δ (x, t) = i=1 u i (t)Ni (x).
ui
ui 1 Ni Ni 1
Ωi ui 1 ui 1 1
xi x i 1 xi x i 1 x x i 1
Figure 7. Finite element expansion bases.
This is a common technique in the FEM because it reduces the smoothness

requirements on u and it also makes the matrix of the discretized system sym-
metric. In two and three dimensions we would use Gauss’ divergence theorem
to obtain a similar result.
The application of the boundary conditions in the FEM deserves attention.
The imposition of the Neumann boundary condition u x (1) = g is straightfor-
ward, we simply substitute the value in Eq. (39). This is a very natural way
of imposing Neumann boundary conditions which also leads to symmetric
matrices, unlike the FDM. The Dirichlet boundary condition u(0) = α can be
applied by imposing u 1 = α and requiring that w(0) = 0. In general, we will
impose that the weight functions

w(x) are zero at the Dirichlet boundaries.
Letting u(x) ≈ u δ (x) = Nj=1 u j N j (x) and w(x) = Ni (x) then Eq. (39) be-
comes
1 N 1
dNi dN j
− (x) uj (x) dx = Ni (x) s(x) dx (40)
dx j =1
dx
0 0
for i =2, . . . , N . This represents a linear system of N − 1 equations with N − 1

unknowns: {u 2 , . . . , u N }. Let us proceed to calculate the integral terms corre-
sponding to the ith equation. We calculate the integrals in Eq. (40) as sums of
integrals over the elements i . The basis functions have compact support, as
shown in Fig. 6. Their value and their derivatives are different from zero only
on the elements containing the node i, i.e.,



x − xi−1

 xi−1 < x < xi
xi−1
Ni (x) =



xi+1 − x
 xi < x < xi+1
xi



1

 xi−1 < x < xi
dNi (x) xi−1
=
dx 
 −1

 xi < x < xi+1
xi
with xi−1 = xi − xi−1 and xi = xi+1 − xi . This means that the only integrals
different from zero in (40) are
xi
xi+1
dNi dNi−1 dNi dNi dNi dNi+1
− u i−1 + ui − ui + u i+1 dx
dx dx dx dx dx dx
xi−1 xi
xi
xi+1
= Ni s dx + Ni s dx (41)
xi−1 xi
The right-hand side of this equation expressed as

xi
xi+1
x − xi−1 xi+1 − x
F= s(x) dx + s(x) dx
xi−1 xi
xi−1 xi
can be evaluated using a simple integration rule like the trapezium rule

xi+1
g(xi ) + g(xi+1 )
g(x) dx ≈ xi
2
xi
and it becomes

xi−1 xi
F= + si .
2 2
Performing the required operations in the left-hand side of Eq. (41) and includ-
ing the calculated valued of F leads to the FEM discrete form of the equation
as
u i − u i−1 u i+1 − u i xi−1 + xi
− + = si .
xi−1 xi 2
Here if we assume that xi−1 = xi = x then the equispaced approximation

becomes
u i+1 − 2u i + u i−1
= x si
x
which is identical to the finite difference formula. We note, however, that the
general FE formulation did not require the assumption of an equispaced mesh.
In general the evaluation of the integral terms in this formulation are more
efficiently implemented by considering most operations in a standard element
st = {−1 ≤ x ≤ 1} where a mapping is applied from the element i to the
standard element st . For more details on the general formulation see Ref. [4].
3.3.3. Finite volume method (FVM)
The integral form of the one-dimensional linear advection equation is given

by Eq. (1) with f (u) = au and S = 0. Here the region of integration is taken to
be a control volume i , associated with the point of coordinate xi , represented
by xi−(1/2) ≤ x ≤ xi+(1/2) , following the notation of Fig. 4, and the integral
form is written as

xi+(1/2)
xi+(1/2)
u t dx + f x (u) dx = 0. (42)
xi−(1/2) xi−(1/2)
This expression could also been obtained from the weighted residual form (4)
by selecting a weight w(x) such that w(x) = 1 for xi−(1/2) ≤ x ≤ xi+(1/2) and
w(x) = 0 elsewhere. The last term in Eq. (42) can be evaluated analytically to
obtain

xi+(1/2)

f x (u) dx = f u i+(1/2) − f u i−(1/2)
xi−(1/2)
and if we approximate the first integral using the mid-point rule we finally
have the semi-discrete form

u t |i xi+(1/2) − xi−(1/2) + f u i+(1/2) − f u i−(1/2) = 0.
This approach produces a conservative scheme if the flux on the boundary
of one cell equals the flux on the boundary of the adjacent cell. Conservative
schemes are popular for the discretization of hyperbolic equations since, if
they converge, they can be proven (Lax-Wendroff theorem) to converge to a
weak solution of the conservation law.
3.3.4. Comparison of FVM and FDM
To complete our comparison of the different techniques we consider the

FVM discretization of the elliptic equation u x x = s. The FVM integral form of
this equation over a control volume i = {xi−(1/2) ≤ x ≤ xi+(1/2) } is

xi+(1/2)
xi+(1/2)
u x x dx = s dx.
xi−(1/2) xi−(1/2)
Evaluating exactly the left-hand side and approximating the right-hand side by
the mid-point rule we obtain

u x xi+(1/2) − u x xi−(1/2) = xi+(1/2) − xi−(1/2) si . (43)
If we approximate u(x) as a linear function between the mesh points i − 1 and
i, we have
u i − u i−1 u i+1 − u i
u x |i−(1/2) ≈ , u x |i+(1/2) ≈ ,
xi − xi−1 xi+1 − xi
and introducing these approximations into Eq. (43) we now have
u i+1 − u i u i − u i−1
− = (xi+(1/2) − xi−(1/2) ) si .
xi+1 − xi xi − xi−1
If the mesh is equispaced then this equation reduces to
u i+1 − 2u i + u i−1
= x si ,
x
which is the same as the FDM and FEM on an equispaced mesh.
Once again we see the similarities that exist between these methods
although some assumptions in the construction of the FVM have been made.
FEM and FVM allow a more general approach to non-equispaced meshes
(although this can also be done in the FDM). In two and three dimensions,
curvature is more naturally dealt with in the FVM and FEM due to the integral
nature of the equations used.
4. High Order Discretizations: Spectral Element/ p-Type

Finite Elements
All of the approximations methods we have discussed this far have dealt
with what is typically known as the h-type approximation. If h = x denotes
the size of a finite difference spacing or finite elemental regions then conver-
gence of the discrete approximation to the PDE is achieved by letting h → 0.
An alternative method is to leave the mesh spacing fixed but to increase the
polynomial order of the local approximation which is typically denoted by p
or the p-type extension.
We have already seen that higher order finite difference approximations
can be derived by fitting polynomials through more grid points. The draw-
back of this approach is that the finite difference stencil gets larger as the
order of the polynomial approximation increases. This can lead to difficulties
when enforcing boundary conditions particularly in multiple dimensions. An
alternative approach to deriving high order finite differences is to use com-
pact finite differences where a Padé approximation is used to approximate the
derivatives.
When using the finite element method in an integral formulation, it is
possible to develop a compact high-order discretization by applying higher
order polynomial expansions within every elemental region. So instead of us-
ing just a linear element in each piecewise approximation of Fig. 6 we can
use a polynomial of order p. This technique is commonly known as p-type
finite element in structural mechanics or the spectral element method in fluid
mechanics. The choice of the polynomial has a strong influence on the nu-
merical conditioning of the approximation and we note that the choice of an
equi-spaced Lagrange polynomial is particularly bad for p > 5. The two most
commonly used polynomial expansions are Lagrange polynomial based on the
Gauss–Lobatto–Legendre quadratures points or the integral of the Legendre
polynomials in combination with the linear finite element expansion. These
two polynomial expansions are shown in Fig. 8. Although this method is more
(a) (b)
1 1 1 1
0 0 0 0
1 1 1 1
1 1 1 1
0 0 0 0
1 1 1 1
1 1 1 1
0 0 0 0
1 1 1 1
Figure 8. Shape of the fifth order ( p = 5) polynomial expansions typically used in (a) spectral
element and (b) p-type finite element methods.
involved to implement, the advantage is that for a smooth problem (i.e., one
where the derivatives of the solution are well behaved) the computational cost
increases algebraically whilst the error decreases exponentially fast. Further
details on these methods can be found in Refs. [5, 6].
5. Numerical Difficulties
The discretization of linear elliptic equations with either FD, FE or FV

methods leads to non-singular systems of equations that can easily solved by
standard methods of solution. This is not the case for time-dependent problems
where numerical errors may grow unbounded for some discretization. This is
perhaps better illustrated with some examples.
Consider the parabolic problem represented by the diffusion equation u t −
u x x = 0 with boundary conditions u(0) = u(1) = 0 solved using the scheme
(36) with b = 1 and x = 0.1. The results obtained with t = 0.004 and 0.008
are depicted in Figs. 9(a) and (b), respectively. The numerical solution (b)
corresponding to t = 0.008 is clearly unstable.
A similar situation occurs in hyperbolic problems. Consider the one-
dimensional linear advection equation u t + au x = 0; with a > 0 and various
explicit approximations, for instance the backward in space, or upwind,
scheme is
u in+1 − u in u n − u i−1
n
+a i = 0 ⇒ u in+1 = (1 − σ )u in + σ u i−1
n
, (44)
t x
the forward in space, or downwind, scheme is
u in+1 − u in u n − u in
+ a i+1 =0 ⇒ u in+1 = (1 + σ )u in − σ u i+1
n
, (45)
t x
(a) (b)
0.3 0.3
t0.20 t0.20
t0.24 t0.24
t0.28 t0.28
t0.32 t0.32
0.2 0.2
0.1 0.1
u(x,t)
u(x,t)
0 0
0.1 0.1
0.2 0.2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
Figure 9. Solution to the diffusion equation u t + u x x = 0 using a forward in time and centred
in space finite difference discretization with x = 0.1 and (a) t = 0.004, and (b) t = 0.008.
The numerical solution in (b) is clearly unstable.
u(x,t)
 a
0
 x ≤ −0.2 1.0

1 + 5x −0.2 ≤ x ≤ 0
u(x, 0) =

 1 − 5x 0 ≤ x ≤ 0.2
 0.0
0 x ≥ 0.2 x
0.2 0.2
Figure 10. A triangular wave as initial condition for the advection equation.
and, finally, the centred in space is given by
u in+1 − u in u n − u i−1
n
σ n
+ a i+1 =0 ⇒ u in+1 = u in − (u − u i−1
n
)
t 2x 2 i+1
(46)
where σ = (at/x) is known as the Courant number. We will see later that
this number plays an important role in the stability of hyperbolic equations.
Let us obtain the solution of u t + au x = 0 for all these schemes with the initial
condition given in Fig. 10.
As also indicated in Fig. 10, the exact solution is the propagation of this
wave form to the right at a velocity a. Now we consider the solution of the
three schemes at two different Courant numbers given by σ = 0.5 and 1.5. The
results are presented in Fig. 11 and we observe that only the upwinded scheme
when σ ≤ 1 gives a stable, although diffusive, solution. The centred scheme
when σ = 0.5 appears almost stable but the oscillations grow in time leading
to an unstable solution.
6. Analysis of Numerical Schemes

We have seen that different parameters, such as the Courant number, can
effect the stability of a numerical scheme. We would now like to set up a
more rigorous framework to analyse a numerical scheme and we introduce the
concepts of consistency, stability and Convergence of a numerical scheme.
6.1. Consistency
A numerical scheme is consistent if the discrete numerical equation tends
to the exact differential equation as the mesh size (represented by x and t)
tends to zero.
1 3
0.9
2
0.8
0.7
1
0.6
u(x,t)
u(x,t)
0.5 0
0.4
1
0.3
0.2
2
0.1
0 3
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
σ 0.5 σ 1.5
1 30
20
2
10
1
0
u(x,t)
0
10
u(x,t)
1
20
2 30
3 40
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
x x
σ 0.5 σ 1.5
1.2 3
1
2
0.8
1
0.6
0.4 0
u(x,t)
u(x,t)
0.2
1
0
2
0.2
3
0.4
0.6 4
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
x x
σ 0.5 σ 1.5
Figure 11. Numerical solution of the advection equation u t + au x = 0. Dashed lines: initial
condition. Dotted lines: exact solution. Solid line: numerical solution.
Consider the centred in space and forward in time finite diference approxi-
mation to the linear advection equation u t + au x = 0 given by Eq. (46). Let us
consider Taylor expansions of u in+1 , u i+1
n n
and u i−1 around (xi , t n ) as
t 2
u in+1 = u in + t u t |in + u tt |in + · · ·
2
x 2 x 3
n
u i+1 = u in + x u x |in + u x x |in + u x x x |in + · · ·
2 6
x 2 x 3
n
u i−1 = u in − x u x |in + u x x |in − u x x x |in + · · ·
2 6
Substituting these expansions into Eq. (46) and suitably re-arranging the terms
we find that
u in+1 − u in u n − u i−1n
+ a i+1 − (u t + au x )|in = T (47)
t 2x
where T is known as the truncation error of the approximation and is
given by
t x 2
T = u tt |in + au x x x |in + O(t 2 , x 4 ).
2 6
The left-hand side of this equation will tend to zero as t and x tend to zero.
This means that the numerical scheme (46) tends to the exact equation at point
xi and time level t n and therefore this approximation is consistent.
6.2. Stability
We have seen in the previous numerical examples that errors in numeri-

cal solutions can grow uncontrolled and render the solution meaningless. It
is therefore sensible to require that the solution is stable, this is that the dif-
ference between the computed solution and the exact solution of the discrete
equation should remain bounded as n → ∞ for a given x.
6.2.1. The Courant–Friedrichs–Lewy (CFL) condition
This is a necessary condition for stability of explicit schemes devised by

Courant, Friedrichs and Lewy in 1928.
Recalling the theory of characteristics for hyperbolic systems, the domain
of dependence of a PDE is the portion of the domain that influences the so-
lution at a given point. For a scalar conservation law, it is the characteristic
passing through the point, for instance, the line P Q in Fig. 12. The domain
of dependence of a FD scheme is the set of points that affect the approximate
solution at a given point. For the upwind scheme, the numerical domain of
dependence is shown as a shaded region in Fig. 12.
The CFL criterion states that a necessary condition for an explicit FD
scheme to solve a hyperbolic PDE to be stable is that, for each mesh point,
the domain of dependence of the FD approximation contains the domain of
dependence of the PDE.
(a) (b)
t ∆x t ∆x
Characteristic a∆t
P P
a∆t
∆t
∆t
x x
Q Q
Figure 12. Solution of the advection equation by the upwind scheme. Physical and numerical
domains of dependence: (a) σ = (at/x) > 1, (b) σ ≤ 1.
For a Courant number σ = (at/x) greater than 1, changes at Q will

affect values at P but the FD approximation cannot account for this.
The CFL condition is necessary for stability of explicit schemes but it is
not sufficient. For instance, in the previous schemes we have that the upwind
FD scheme is stable if the CFL condition σ ≤ 1 is imposed. The downwind
FD scheme does not satisfy the CFL condition and is unstable. However, the
centred FD scheme is unstable even if σ ≤ 1.
6.2.2. Von Neumann (or Fourier) analysis of stability
The stability of FD schemes for hyperbolic and parabolic PDEs can be

analysed by the von Neumann or Fourier method. The idea behind the method
is the following. As discussed previously the analytical solutions of the model
diffusion equation u t − b u x x = 0 can be found in the form
∞

u(x, t) = eβm t e I km x
m=−∞
if βm + b km2 = 0. This solution involves a Fourier series in space and an expo-

nential decay in time since βm ≤ 0 for b > 0. Here we have included the√ com-
plex version of the Fourier series, e I km x = cos km x + I sin km x with I = −1,
because this simplifies considerably later algebraic manipulations.
To analyze the growth of different Fourier modes as they evolve under the
numerical scheme we can consider each frequency separately, namely
u(x, t) = eβm t e I km x .
A discrete version of this equation is u in = u(xi , t n ) = eβm t e I km xi . We can take,

n
without loss of generality, xi = ix and t n = nt to obtain

n
u in = eβm nt e I km ix = eβm t e I km ix .
The term e I km ix = cos(km ix) + I sin(km ix) is bounded and, therefore, any
growth in the numerical solution will arise from the term G = eβm t , known
as the amplification factor. Therefore the numerical method will be stable, or
the numerical solution u in bounded as n → ∞, if |G| ≤ 1 for solutions of the
form
u in = G n e I km ix .
We will now proceed to analyse, using the von Neummann method, the stabil-
ity of some of the schemes discussed in the previous sections.
Example 1 Consider the explicit scheme (36) for the diffusion equation
u t − bu x x = 0 expressed here as
bt
u in+1 = λu i−1
n
+ (1 − 2λ)u in + λu i+1
n
; λ= .
x 2
We assume u in = G n e I km ix and substitute in the equation to get
G = 1 + 2λ [cos(km x) − 1] .
Stability requires |G| ≤ 1. Using −2 ≤ cos(km x) − 1 ≤ 0 we get 1 − 4λ ≤

G ≤ 1 and to satisfy the left inequality we impose
1
−1 ≤ 1 − 4λ ≤ G =⇒ λ≤ .
2
This means that for a given grid size x the maximum allowable timestep is
t = (x 2 /2b).
Example 2 Consider the implicit scheme (37) for the diffusion equation
u t − bu x x = 0 expressed here as
bt
λu i−1
n+1
+ −(1 + 2λ)u in+1 + λu i+1
n+1
= −u in ; λ= .
x 2
The amplification factor is now
1
G=
1 + λ(2 − cos βm )
and we have |G| < 1 for any βm if λ > 0. This scheme is therefore uncondi-
tionally stable for any x and t. This is obtained at the expense of solving
a linear system of equations. However, there will still be restrictions on x
and t based on the accuracy of the solution. The choice between an explicit
or an implicit method is not always obvious and should be done based on the
computer cost for achieving the required accuracy in a given problem.
Example 3 Consider the upwind scheme for the linear advection equa-
tion u t + au x = 0 with a > 0 given by
at
u in+1 = (1 − σ )u in + σ u i−1
n
; σ= .
x
Let us denote βm = km x and introduce the discrete Fourier expression in the
upwind scheme to obtain
G = (1 − σ ) + σ e−Iβm
The stability condition requires |G| ≤ 1. Recall that G is a complex number

G = ξ + I η so
ξ = 1 − σ + σ cos βm ; η = −σ sin βm
This represents a circle of radius σ centred at 1 − σ . The stability condition

requires the locus of the points (ξ, η) to be interior to a unit circle ξ 2 + η2 ≤ 1.
If σ < 0 the origin is outside the unit circle, 1 − σ > 1, and the scheme is
unstable. If σ > 1 the back of the locus is outside the unit circle 1 − 2σ < 1 and
it is also unstable. Therefore, for stability we require 0 ≤ σ ≤ 1, see Fig. 13.
Example 4 The forward in time, centred in space scheme for the advec-
tion equation is given by
σ n at
u in+1 = u in − (u − u i−1
n
); σ= .
2 i+1 x
1 λ
1
G
λ ξ
Figure 13. Stability region of the upwind scheme.

The introduction of the discrete Fourier solution leads to

σ
G = 1 − (e Iβm − e−Iβm ) = 1 − I σ sin βm
2
Here we have |G|2 = 1 + σ 2 sin2 βm > 1 always for σ/=0 and it is therefore
unstable. We will require a different time integration scheme to make it stable.
6.3. Convergence: Lax Equivalence Theorem
A scheme is said to be convergent if the difference between the computed

solution and the exact solution of the PDE, i.e., the error E in = u in − u(xi , t n ),
vanishes as the mesh size is decreased. This is written as
lim |E in | = 0
x,t→0
for fixed values of xi and t n . This is the fundamental property to be sought

from a numerical scheme but it is difficult to verify directly. On the other hand,
consistency and stability are easily checked as shown in the previous sections.
The main result that permits the assessment of the convergence of a scheme
from the requirements of consistency and stability is the equivalence theorem
of Lax stated here without proof:
Stability is the necessary and sufficient condition for a consistent linear FD
approximation to a well-posed linear initial-value problem to be convergent.
7. Suggestions for Further Reading
The basics of the FDM are presented a very accessible form in Ref. [7].
More modern references are Refs. [8, 9].
An elementary introduction to the FVM can be consulted in the book by
Versteeg and Malalasekera [10]. An in-depth treatment of the topic with an
emphasis on hyperbolic problems can be found in the book by Leveque [2].
Two well established general references for the FEM are the books of
Hughes [4] and Zienkiewicz and Taylor [11]. A presentation from the point
of view of structural analysis can be consulted in Cook et al. [11]
The application of p-type finite element for structural mechanics is dealt
with in book of Szabo and Babus̆ka [5]. The treatment of both p-type and spec-
tral element methods in fluid mechanics can be found in book by Karniadakis
and Sherwin [6].
A comprehensive reference covering both FDM, FVM and FEM for fluid
dynamics is the book by Hirsch [13]. These topics are also presented using a
more mathematical perspective in the classical book by Quarteroni and Valli
[14].
References
[1] J. Bonet and R. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis.
Cambridge University Press, 1997.
[2] R. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University
Press, 2002.
[3] W. Cheney and D. Kincaid, Numerical Mathematics and Computing, 4th edn.,
Brooks/Cole Publishing Co., 1999.
[4] T. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element
Analysis, Dover Publishers, 2000.
[5] B. Szabo and I. Babus̆ka, Finite Element Analysis, Wiley, 1991.
[6] G.E. Karniadakis and S. Sherwin, Spectral/hp Element Methods for CFD, Oxford
University Press, 1999.
[7] G. Smith, Numerical Solution of Partial Differential Equations: Finite Diference
Methods, Oxford University Press, 1985.
[8] K. Morton and D. Mayers, Numerical Solution of Partial Differential Equations,
Cambridge University Press, 1994.
[9] J. Thomas, Numerical Partial Differential Equations: Finite Difference Methods,
Springer-Verlag, 1995.
[10] H. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynam-
ics. The Finite Volume Method, Longman Scientific & Technical, 1995.
[11] O. Zienkiewicz and R. Taylor, The Finite Element Method: The Basis, vol. 1, Butter-
worth and Heinemann, 2000.
[12] R. Cook, D. Malkus, and M. Plesha, Concepts and Applications of Finite Element
Analysis, Wiley, 2001.
[13] C. Hirsch, Numerical Computation of Internal and External Flows, vol. 1, Wiley,
1988.
[14] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equa-
tions, Springer-Verlag, 1994.

8.2 Finite Difference, Finite Element and Finite Volume Methods For Partial Differential Equations

Uploaded by

Copyright:

Available Formats

8.2 Finite Difference, Finite Element and Finite Volume Methods For Partial Differential Equations

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

8.2 Finite Difference, Finite Element and Finite Volume Methods For Partial Differential Equations

Uploaded by

Copyright:

Available Formats

8.

Joaquim Peiró and Spencer Sherwin

There are three important steps in the computational modelling of any

After the selection of an appropriate mathematical model, together with

1. Conservation Laws: Integral and Differential Forms

to as conservation laws. These are derived by invoking the conservation of

This is known as the weak form of the conservation law.

Although the above formulation is more commonly used in fluid mechan-

where δW denotes the virtual work done by an arbitrary virtual velocity δ ,

2. Model Equations and their Classification

2.1. Elliptic Equations

2.2. Parabolic Equations

which is parabolic. In addition to appropriate boundary conditions of the form

Figure 1. Rate of decay of the solution to the diffusion equation.

2.3. Hyperbolic Equations

where a represents a constant velocity. The above equation is also clearly

Figure 2. Solution of the linear advection equation.

The propagation velocity can also be a function of space, i.e., a = a(x) or

Figure 3. Formation of discontinuities in the Burgers’ equation.

Figure 4. Discretization of the domain.

shown in Fig. 4, the numerical solution that we are seeking is represented by a

3.1. The Finite Difference Method (FDM)

3.1.1. Errors in the FDM

The analysis of these approximations is performed by using Taylor expan-

To illustrate how this can be used to analyse finite difference approxima-

3.1.2. Derivation of approximations using Taylor expansions

The procedure described in the previous section can be easily transformed

Putting (22), (23) into (21) we get

We require three independent conditions to calculate the three unknowns α,

Other approximations can be obtained by selecting a different set of points,

3.1.3. Higher-order derivatives

In general, we can derive an approximation of the second derivative using

Using similar arguments to those of the previous section we impose

3.1.4. Finite differences through polynomial interpolation

interpolation. Similarly, for the second derivative we have

3.1.5. Finite difference solution of PDEs

We consider the FDM approximation to the solution of the elliptic equation

This matrix system can be written in abridged form as Au = s. The matrix

and include the Neumann boundary condition (33) to obtain

3.2. Time Integration

In this section we address the problem of solving time-dependent PDEs

with an initial condition u(x, 0) = u 0 (x) and time-dependent boundary condi-

functions of t. Assume, as before, a mesh or spatial discretization of the

3.2.1. Method of lines

3.2.2. Finite differences in time

The method of finite differences can be applied to time-dependent prob-

u in ≈ u(xi , t n ) and the finite difference approximation of the time derivative

3.3. Discretizations Based on the Integral Form

3.3.1. The finite element method (FEM)

Here we discretize the region of interest = {x : 0 ≤ x ≤ 1} into N − 1

3.3.2. Galerkin FEM

Figure 7. Finite element expansion bases.

This is a common technique in the FEM because it reduces the smoothness

for i =2, . . . , N . This represents a linear system of N − 1 equations with N − 1

The right-hand side of this equation expressed as

Here if we assume that xi−1 = xi = x then the equispaced approximation

3.3.3. Finite volume method (FVM)

The integral form of the one-dimensional linear advection equation is given

3.3.4. Comparison of FVM and FDM

Here if we assume that xi−1 = xi = x then the equispaced approximation

For a Courant number σ = (at/x) greater than 1, changes at Q will

without loss of generality, xi = ix and t n = nt to obtain

Stability requires |G| ≤ 1. Using −2 ≤ cos(km x) − 1 ≤ 0 we get 1 − 4λ ≤