Math 257: Finite difference methods

1 Finite Differences
Remember the definition of a derivative
f (x + ∆x) − f (x)
f 0 (x) = lim (1)
∆x→0 ∆x
Also recall Taylor’s formula:
∆x2 00 ∆x3 (3)
f (x + ∆x) = f (x) + ∆xf 0 (x) + f (x) + f (x) + . . . (2)
2! 3!
or, with −∆x instead of +∆x:
∆x2 00
0 ∆x3 (3)
f (x − ∆x) = f (x) − ∆xf (x) + f (x) − f (x) + . . . (3)
2! 3!

Forward difference
On a computer, derivatives are approximated by finite difference expressions; rearrang-
ing (2) gives the forward difference approximation
f (x + ∆x) − f (x)
= f 0 (x) + O(∆x), (4)
∆x
where O(∆x) means ‘terms of order ∆x’, ie. terms which have size similar to or smaller
than ∆x when ∆x is small.1 So the expression on the left approximates the derivative of
f at x, and has an error of size ∆x; the approximation is said to be ‘first order accurate’.
1
The strict mathematical definition of O(δ) is to say that it represents terms y such that
y
lim is finite.
δ→0 δ

So δ 2 and δ are O(δ), but δ 1/2 or 1, for instance, are not. Another commonly used notation is o(δ),
which represents terms y for which y
lim = 0.
δ→0 δ

In other words, o(δ) represents terms that are strictly smaller than δ as δ → 0. So δ 2 is o(δ), but δ is
not.

1
Backward difference
Rearranging (3) similarly gives the backward difference approximation
f (x) − f (x − ∆x)
= f 0 (x) + O(∆x), (5)
∆x
which is also first order accurate, since the error is of order ∆x.

Centered difference
Combining (2) and (3) gives the centered difference approximation
f (x + ∆x) − f (x − ∆x)
= f 0 (x) + O(∆x2 ), (6)
2∆x
which is ‘second order accurate’, because the error this time is of order ∆x2 .

Second derivative, centered difference

Adding (2) and (3) gives
∆x4 (4)
f (x + ∆x) + f (x − ∆x) = 2f (x) + ∆x2 f 00 (x) + f (x) + . . . (7)
12
Rearranging this therefore gives the centered difference approximation to the second
derivative:
f (x + ∆x) − 2f (x) + f (x − ∆x)
= f 00 (x) + O(∆x2 ), (8)
∆x2
which is second order accurate.

2 Heat Equation
Dirichlet boundary conditions
To find a numerical solution to the heat equation
∂u ∂ 2u
= α2 2 , 0 < x < L, (9)
∂t ∂x
u(0, t) = A, u(L, t) = B, u(x, 0) = f (x), (10)
approximate the time derivative using forward differences, and the spatial derivative
using centered differences;
u(x, t + ∆t) − u(x, t) u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t)
= α2 + O(∆t, ∆x2 ). (11)
∆t ∆x2
2
 t
unk+1

∆t ukn−1 ukn+1
ukn

u2−1 u20 u2N u2N +1

u1−1 u10 u1N u1N +1

u00 u01 u02 u0N

x0 x1 x2 xN
x=0 ∆x x=L

Figure 1: Mesh points and finite difference stencil for the heat equation. Blue points are
prescribed the initial condition, red points are prescribed by the boundary conditions.
For Neumann boundary conditions, fictional points at x = −∆x and x = L + ∆x can
be used to facilitate the method.

This approximation is second order accurate in space and first order accurate in time.
The use of the forward difference means the method is explicit, because it gives an
explicit formula for u(x, t + ∆t) depending only on the values of u at time t.
Divide the interval 0 < x < L into N + 1 evenly spaced points, with spacing ∆x;
ie. xn = n∆x, for n = 0, 1, . . . , N , and divide time into discrete levels tk = k∆t for
k = 0, 1, . . .. Then seek the solution by finding the discrete values

ukn = u(xn , tk ). (12)

From (11), these satisfy the equations

uk+1 − ukn uk − 2ukn + ukn−1

n
= α2 n+1 , (13)
∆t ∆x2
or, rearranging,
α2 ∆t k
uk+1 = ukn + k k


n u n+1 − 2un + u n−1 . (14)
∆x2
If the values of un at timestep k are known, this formula gives all the values at timestep
k + 1, and it can then be iterated again and again to march forward in time.
The initial condition gives
u0n = f (xn ), (15)

3
for all n, and the boundary conditions require

uk0 = A, ukN = B, (16)

for all values of k > 0 (equation (14) only has to be solved for 1 ≤ n ≤ N − 1).

Stability
Note, this method will only be stable, provided the condition

α2 ∆t 1
2
≤ , (17)
∆x 2
is satisfied; otherwise it will not work.

Neumann boundary conditions

To apply
∂u
(0, t) = C, (18)
∂x
instead of (10), notice that the centered difference version of this boundary condition
would be
u(0 + ∆x, t) − u(0 − ∆x, t)
= C, (19)
2∆x
and therefore
u(−∆x, t) = u(∆x, t) − 2∆xC. (20)
x = −∆x is outside the domain of interest, but knowing the value of u(−∆x, t) there
allows the discretised equation (14) to be used also for x = 0 (n = 0), and it becomes

α2 ∆t
uk+1 = uk0 + k k


0 2u 1 − 2u 0 − 2∆xC . (21)
∆x2
So in this case we must solve (14) for 1 ≤ n ≤ N − 1, and (21) for n = 0 at each
timestep. If there is a derivative condition at x = L the same procedure is followed and
an equation similar to (21) must be solved for n = N too.
A handy way to implement this type of boundary condition, which enables the same
formula (14) to be used for all points, is to introduce ‘fictional’ mesh points for n = −1
and n = N + 1, and to prescribe the value uk−1 given by (20) (or the equivalent for ukN +1 )
at those points.

4
3 Wave equation
For the wave equation,
∂ 2u 2
2∂ u
= c , 0 < x < L, (22)
∂t2 ∂x2
u(0, t) = 0, u(L, t) = 0, (23)
∂u
u(x, 0) = f (x), (x, 0) = g(x), (24)
∂t
discretise x into N + 1 evenly spaced mesh points xn = n∆x, discretise time into evenly
spaced time levels tk = k∆t, and seek the solution at these mesh points; ukn = u(xn , tk ).
Using centered difference approximations to the two second derivatives, the discrete
equation is
k k k
uk+1
n − 2ukn + uk−1
n 2 un+1 − 2un + nn−1
= c + O(∆x2 , ∆t2 ). (25)
∆t2 ∆x2
This can be rearranged to give

c2 ∆t2 k
uk+1 = 2ukn − uk−1 k k


n n + un+1 − 2un + nn−1 , (26)
∆x2
which gives an explicit method to calculate uk+1
n in terms of the values at the previous
two timesteps k and k − 1. This method is second order accurate in space and time - it
is sometimes referred to as the ‘leap-frog’ method.

Boundary conditions
To apply Dirichlet boundary conditions (23), the values of uk+1 0 and uk+1
N are simply
prescribed to be 0; there is no need to solve an equation for these end points.
If the boundary conditions are Neumann, on the other hand:
∂u ∂u
(0, t) = 0, (L, t) = 0, (27)
∂x ∂x
then (26) can still be solved for n = 0 and n = N if we use the centered difference
approximations to the boundary conditions in order to determine the ‘fictional’ values
uk−1 and ukN +1 respectively. At x = 0, for instance, the discretised condition (27) is

uk−1 − uk1
= 0, (28)
2∆x
and therefore uk−1 = uk1 . Similarly ukN +1 = ukN −1 .

5
 t
unk+1

∆t ukn−1 ukn+1
ukn

unk−1

u2−1 u20 u2N u2N +1

u1−1 u10 u1N u1N +1

u00 u01 u02 u0N

x0 x1 x2 xN
x=0 ∆x x=L

Figure 2: Mesh points and finite difference stencil for the wave equation.

Initial conditions
To apply the initial conditions, note that to use (26) for the first timestep k = 1, we
need to know the values of u0n and also u−1 0
n . The values of un follow from the initial
condition on u(x, 0);
u0n = f (xn ). (29)
The values of u−1
n come from considering the centered difference approximation to the
derivative condition (24);

u(x, 0 + ∆t) − u(x, 0 − ∆t)

= g(x), (30)
2∆t
from which
u−1 1
n = un − 2∆tg(xn ). (31)
Substituting this into the discrete equation (26), and rearranging, gives the formula

1 c2 ∆t2 0
u1n = u0n + un+1 − 2u0n + n0n−1 + ∆tg(xn ),


2
(32)
2 ∆x
for the first timestep. Once the solution has been initialised in this way, all subsequent
timesteps can be made using (26).

6
Stability
This method is stable provided
c∆t
≤ 1, (33)
∆x
which is often called the Courant-Friedrichs-Levy (or CFL) condition.

4 Laplace’s equation
Laplace’s equation on a square domain is
∂ 2u ∂ 2u
+ = 0, 0 < x < L, 0 < y < L, (34)
∂x2 ∂y 2
with boundary conditions, for instance,

u(0, y) = 0, u(L, y) = 0, u(x, 0) = f (x), u(x, L) = 0. (35)

Choose a mesh with spacing ∆x in the x direction, ∆y in the y direction (often it

is sensible to choose ∆x = ∆y), so grid points are xn = n∆x for n = 0, 1, . . . N , and
ym = m∆y for m = 0, 1, . . . M . Then seek solution values un m = u(xn , ym ). Using
second order accurate centered differences for the two derivatives, the discrete equation
is
un+1 m − 2un m + un−1 m un m+1 − 2un m + un m−1
+ = O(∆x2 , ∆y 2 ), (36)
∆x2 ∆y 2
and if ∆x = ∆y this simplifies to
1
un m = (un+1 m + un−1 m + un m+1 + un m−1 ) (37)
4
Note this equation shows that the solution has the property that the value at each point
(xn , ym ) is the average of the values at its four neighbouring points. This is an important
property of Laplace’s equation.

Boundary conditions
The boundary conditions (35) determine the values un m on each of the four boundaries,
so (37) does not have to be solved at these mesh points:

u0 m = 0, uN m = 0, un 0 = f (xn ), un M = 0. (38)

Neumann boundary conditions can be incorporated by calculating values for u−1 m (for
instance), as for the heat equation.

7
 u0 M uN M
y = L yM

un m+1

∆y un−1 m un m un+1 m

un m−1

u0 1
y1

u0 0 u1 0 u2 0 uN 0
y = 0 y0
x0 x1 x2 xN
x=0 ∆x x=L

Figure 3: Mesh points and finite difference stencil for Laplace’s equation.

Jacobi Iteration
Note (37) is not an explicit formula, since the solution at each point depends on the
unknown values at other points, and it is therefore harder to solve than the previous
explicit methods. One way to do it, however, is to take a guess at the solution (eg.
un m = 0 everywhere), and then go through each of the mesh points updating the
solution according to (37) by taking the values of the previous guess on the right hand
side. Iterating this process many times, the successive approximations will hopefully
change by less and less, and the values they converge to provide the solution. Given an
(0) (1) (k) (k+1)
initial guess un m , the successive iterations un m , . . . , un m , un m , . . ., are given by
1  (k) (k) (k) (k)

u(k+1) = u + u + u + u n m−1 . (39)
nm
4 n+1 m n−1 m n m+1

Here the superscript (k) denotes the values for the kth iteration. The process is continued
until the change between successive iterations is less than some desired tolerance.