Year 4 — Practical Numerical Analysis

Based on lectures by Dr Kathryn Gillow

Notes taken by James Arthur

Michaelmas 2022

These notes are not endorsed by the lecturers, and I have modified them (often significantly) after lectures. They
are nowhere near accurate representations of what was actually lectured, and in particular, all errors are almost
surely mine (especially the typos!).

Contents
1 Introduction 2
1.1 Root Finding (and optimisation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Univariate Root Finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Regular Falsi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.3 Newton Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.4 Polynomial Root Finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.5 Multivariate Root Finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 IVPs 6
2.1 Euler Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Trapezium Rule / Crank Nicolson Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Generalisation – θ-method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Euler via Taylor Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.5 Truncation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Runge-Kutta Schemes 9
3.1 Derivation of Explicit Runge-Kutta Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 First Order PDEs 11

1 Introduction Practical Numerical Analysis

1 Introduction
Problem Sheet by Tuesday at 12 noon. Thursday on the problem sheets from week 2. TA is Robert McDonald.
We will study,
ˆ Rootfinding

ˆ ODEs

– Euler Schemes
– Runga Kutta
– Linear multistep
ˆ Parabolic PDEs - Heat Equation

1.1 Root Finding (and optimisation)

The idea is simple. Find some x such that f (x) = 0. We are being vague on purpose. The stating is simple,
the solving isn’t simple. The cubic formula isn’t so simple. Then for the quintic, then it’s only has closed
form when its Galois group is solvable. These closed forms, may not exist, be unstable, require evaluation of
special function. Hence we need numerical methods.

The relation to optimisation is quite nice. The max or min occurs at a turning point. That is,
∂g
fi := =0
∂xi
and then we have a root finding problem.

1.1.1 Univariate Root Finding

In the 1D case, given f : [a, b] → R then find c ∈ [a, b] such that f (c) = 0. The bisection algorithm uses the
intermediate value theorem (See Analysis / Topology). Thus to find a root of f (c) we find a and b such that
f (a) and f (b) have opposite signs. Then let c = (a + b)/2. Then compute f (c), check and repeat.

We need to decide when to terminate the algorithm, the normal are,

ˆ |b − a| < tol, or,

ˆ |f (c)| < tol.

For convergence, since c ∈ [a, b] the maximum error is b − a. Since the step halves, the error at step n,

b−a
|cn − c| ≤ .
2n
Thus to achieve an accuracy of tolerance requires,

log(b − a) − log(tol)
n≥
log(2)

steps of bisection algorithm.

Here is some pseudocode for this algorithm,

2 James Arthur
1 Introduction Practical Numerical Analysis

Algorithm 1 Bisection Method

Input: a, b, f (x), tol
Output: c
while |b − a| > tol do
c ← (a + b)/2
if f (a)f (c) < 0 then
b←c
else
a←c
end if
end while

1.1.2 Regular Falsi

This algorithm finds a clever guess. We are going to approximate f on [a, b],
x−a b−x
p1 (x) = f (b) + f (a)
b−a b−a
Then p1 (x) has a root at,
af (b) − bf (a)
c=
f (b) − f (a)
and then we use this in the bisection algorithm.
This fails for functions like f (x) = x10 − 1/10. There is a fix,

af (b)/2 − bf (a)
c=
f (b)/2 − f (a)
This is known as the illinois algorithm.

1.1.3 Newton Raphson

Suppose we approximate f (x) by g(x) which is the truncated Taylor series of f about the point xk ,

g(x) = f (xk ) + (x − xk )f 0 (xk )

Then g(x) has a root at,

f (xk )
c = xk −
f 0 (xk )
Thus the newton raphson method requires an intial guess and then iterates,
f (xk )
xk+1 = xk −
f 0 (xk )
until convergence is achieved. This is quicker, quadratic, if we have a guess sufficiently close.
We can variate this to secant method. The iterate is,
xk − xk+1
xk+1 = xk − f (xk )
f (xk ) − f (xk−1 )

This avoids the need for f 0 (xk ). We also have damped newtons,

f (xk )
xk+1 = xk − σk
f 0 (xk )

3 James Arthur
1 Introduction Practical Numerical Analysis

Algorithm 2 Newton-Raphson Method

Input: x0 , f (x), f 0 (x), tol
Output: x∗
k←0
while |f (xk )| > tol do
xk+1 ← xk − f (xk )/f 0 (xk )
k ←k+1
end while
x∗ ← xk

for some σk ∈ (0, 1]. A third is Halley’s method,

2f (xk )f 0 (xk )
xk+1 = xk −
2(f 0 (xk ))2 − f (xk )f 00 (xk )
This pworks if the roots x∗ satisfies f 0 (x∗ ) 6= 0 and is then Newtons method applied to the function
f (x)/ |f 0 (x)|.

1.1.4 Polynomial Root Finding

Suppose we want to find the root of a Polynomial,
n
X
p(z) = ai z i
i=0

where an = 1. Methods based on bisection will only find real roots, one at a time. Methods based on newtons
will find roots one at a time. With a complex initial guess such methods will find complex roots.

In order to find all the roots at once, we can define the companian matrix C ∈ Cn×n as,
 
0 −a0
.. 
1 . . .

. 
C=  
. . . 0 −a

n−2
 
1 −an−1
Then we can consider zI − c,
 
z a0
 .. .. 
−1 . . 
zI − C =  
 .. 
 . z an−2 
−1 z + an−1
and then we can see that det(zI − C) = p(z). Hence the roots of p(z) is equivalent to finding the eigenvalues
of matrix C. This is what matlab roots uses this approach. You may think this is circular, but we can find
eigenvalues in other ways, like the QR algorithm.

1.1.5 Multivariate Root Finding

Now we want to find x such that f (x) = 0 for f : Rn → Rn . Algorithms based on bisections aren’t extendable
to higher dimensions. Newton based algorithms are easier to extend. Again we approximate f (x) by the first
few terms,
f (x + δx) ≈ f (x) + J(x)δ(x)

4 James Arthur
1 Introduction Practical Numerical Analysis

where J is the jacobean. Then we need to find the newton iterate,

J(x)δx = −f (x)

such that,
J(xk )(xk+1 − xk ) = −f (xk )
We can also write the iterate as,
xk+1 = xk − (J(x))−1 f (xk )
thus damped can be written as,
xk+1 = xk − σk (J(x))−1 f (xk )
Here’s the kicker. We don’t really know how many solutions a root finding algorithm has. We need to look
the idea of newton fractals. For example, suppose we have f (z) = z 3 − 1. This
more closely. Hence we bring √
has three roots, z = 1, (−1 ± 3i)/2. By writing z = x + iy and equating real and imaginary parts. We get,

x3 − 3xy 2 − 1 = 0
3x2 y − y 3 = 0

which can be solved using Newtons Method. For each point for (x, y) ∈ [−2, 2] we compute a solution and
record what root it converges to. A newton fractal is what this converges to.

5 James Arthur
2 IVPs Practical Numerical Analysis

2 IVPs
Theorem 2.1 (Picard). Suppose that f (t, u) is a continuous function of t and u in a region Ω = [0, T ) ×
[u0 − α, u0 + α] of the (t, u) plane and that there exists an L > 0 such that,

|f (t, u) − f (t, v)| ≤ L|u − v|, ∀(t, u), (t, v) ∈ Ω

Suppose also that,

MT ≤ α
where M = maxΩ |f |. There is a unique continuously dufferentable function u(t) defined on [0, T ) satisying,

du
= f (u, t)
dt
u(0) = u0

Suppose we want to solve,

du
= f (u, t) t > 0
dt
u(0) = u0

In order to solve this numerically over [0, T ] we define a set of time points at which we wish to approximate
the solution. We set tn = n∆t for n = 0, 1, . . . , N where ∆t = T /N . Then we can integrate,
Z tn+1
u(tn+1 ) = u(tn ) + f (t, u(tn ))dt
tn

Using different approximatations of the integral lead to different numerical scheme

2.1 Euler Methods

Let Un be the numerical approximation to u at tn . For explicit (or forward) Euler we let,
Z tn+1
f (t, u(t))dt ≈ ∆tf (tn , u(tn ))
tn

This gives,
Un+1 = Un + ∆tf (t, Un )
or,
Un+1 − Un
= f (tn , Un ),
∆t
for n = 0, 1, . . . , N − 1 and U0 = u0 .

For implicit Euler, we use Z tn+1

f (t, u(t))dt ≈ ∆tf (tn+1 , u(tn+1 ))
tn

and this gives...

Explicit Euler is very simple, we can just computer the approximations. The implicit euler is more complica-
tion, as we get a non-linear equation. This then becomes a root finding method, but as the solutions are an
iterated solution to an ODE so the initial guess should be just the previous Un .

6 James Arthur
2 IVPs Practical Numerical Analysis

2.2 Trapezium Rule / Crank Nicolson Scheme

Another way to approximate is,
Z tn+1
∆t
f (t, u(t))dt ≈ (f (tn , u(tn )) + f (tn+1 , u(tn+1 )))
tn 2

This then gives a numerical scheme,

∆t
Un+1 = Un + (f (tn , Un ) + f (tn+1 , Un+1 ))
2

2.3 Generalisation – θ-method

Both the implict and explicit euler methods as well as Crank Nicolson are specific θ-methds,
Un+1 − Un
= θf (tn+1 , Un+1 ) + (1 − θ)f (tn , Un )
∆t
and our special cases are,
ˆ θ = 0, Explicit Euler
ˆ θ = 1, Implicit Euler
ˆ θ = 21 , Crank Nicolson.
For all non-zero θ, the method is implicit and a non-linear equation must be solved at each time-step.

2.4 Euler via Taylor Series

THe explicit and implicit euler scheme can be motivated using Taylor series expansions. Consider expanding
u(tn+1 ) about tn . We have,
u(tn+1 ) = u(tn ) + ∆tu0 (tn )
and then rearranging and fidding,

u(tn+1 − u(tn )
+ O(∆t) = f (tn , u(tn ))
∆tn

2.5 Truncation Error

We have just seen that the Euler methods can be derived by truncating taylor series. The truncation for the
θ-method is defined as,
un+1 − un
Tn = = θf (tn+1 , un+1 ) − (1 − θ)f (tn , un )
∆t
The truncation error can be computed by finding taylor series about some point. For θ = 0, the expansions
are usually about t = tn , while for q = 1 we usually go for t = tn+1 and for anything else, it doesnt matter so
we expand about tn+1/2 = tn +t2n+1 = tn + ∆t/2. For explicit we get,

1
Tn = u0 (tn ) − f (tn , u(tn )) + ∆tu00 (τn )
2
Thwn we remember u0 (t) = f (t, u(t)) and so,
1
Tn = ∆tu00 (τn ).
2

7 James Arthur
2 IVPs Practical Numerical Analysis

Then for the theta method, it gets worse. Welcome to Taylor series hell...
∆t
Tn = (1 − 2θ)u00 (tn+1/2 ) + O(∆t2 )
2
We see that the O(∆t2 ) never vanishes and so Tn will never be zero. Thus,
(
O(∆t) for θ 6= 1/2
.
O(∆t2 ) for θ = 1/2

More precisely, we can show that,

∆t 00 (1)

 2 u (τn )
 θ=0
Tn = − ∆t 000 (2)
12 u (τn ) θ = 1/2
 ∆ 00 (2)

− 2 u (τn ) θ=1

The order of a method is defined to be p where p is the largest integer such that Tn = O(∆tp )

8 James Arthur
3 Runge-Kutta Schemes Practical Numerical Analysis

3 Runge-Kutta Schemes
We continue the study of,

u0 (t) = f (t, u)
u(0) = u0

We say if θ 6= 0 we had to solve a non-linear equation, but we can avoid this with more complicated methods.
We also saw that Crank Nicholson was very good on error. We now go from Crank Nicholson, using Explicit
Euler. We call this improved euler,
Un+1 − Un 1
= (f (tn+1 , Un ) + ∆tf (tn , Un ) + f (tn , Un )).
∆t 2
We can write Runge-Kutta as,
s
Un+1 − Un X
= bi ki
∆t i=1

where, k1 = f (tn , Un ) and,  

i−1
X
ki = f tn + ci ∆t, Un + ∆t ai,j kj  ,
j=1

for i = 2, . . . , s. Then ki ’s are known as the stages of the method and the method is often referred to as an
s-stage method. We usually put the coefficients of explicit Runge-Kutta scheme in some Butcher Tableux.
The improved Euler Scheme can be written as,
Un+1 − Un 1
= (k1 + k2 )
∆t 2
where k1 = f (yn , Un ) and k2 = f (tn + ∆t, Un + ∆tk1 ). Another is the modified Euler-Scheme,

Un+1 − Un 1 1
= f (tn + ∆t, Un + ∆tf (tn , Un ))
∆t 2 2
The butcher table is,

Remark. We note that the b entries should sub to 1.

3.1 Derivation of Explicit Runge-Kutta Scheme

The coefficients are chosen to make the methods as high order as possible, we can do this using Taylor series
expansions of the truncation error. THe truncation error is,
s
u(tn+1 ) − u(tn ) X
Tn = = bi k̃i
∆t i=1

where k̃1 = f (tn , u(tn )) and so on. Now death by taylor series. (I am not typing this. )

Theorem 3.1. The order p of an explicit s-stage Runge-Kutta method is bounded by p ≤ s. Further, it is
possible to construct Runge-Kutta methods that achieve this maximal order.

9 James Arthur
3 Runge-Kutta Schemes Practical Numerical Analysis

3.2 Stability
u0 (t) = λu. If λ < 0, the exact solution limits of 0. The explicit euler scheme can be solved,

Un = (1 + λ∆t)n

For λ < 0 we require it to limit to zero. Hence,

−1 < 1 + λ∆t < 1

and so,
2
∆t < .
|λ|
We can do a similar to improved euler. We solve,
1
Un = (1 + λ∆t + (λ∆t)2 )
2
and we see again,
2
∆t < .
|λ|
The interval of absolute stability us the interval of values λ∆t such that limn→∞ Un = 0. Thus for both
explicit and improved Euler, λ∆t ∈ (−2, 0).

10 James Arthur
4 First Order PDEs Practical Numerical Analysis

4 First Order PDEs

We are going to consider a parabolic PDE,
∂u ∂2u
=
∂t ∂x2
with an initial condition u(x, 0) = u0 (x). Consider a uniform set of time steps tm = m∆t for m = 0, 1, 2 . . . .
We also consider xj = a + j∆x, where ∆x = (b − a)/N . We write um j = u(xj , tm ) and we seek to approximate
um m
j by Uj for j = 0, 1, 2 . . . , N and m = 0, 1, 2, . . . . We write,

∂2u u(xj+1 , t) − 2u(xj , t) + u(xj−1 , t)

2
(xj , t) = + O(∆x2 ).
∂x ∆x2
Similarly, we can do fortward and backward differences,

∂u u(x, tm+1 ) − u(x, tm )

(x, tm ) = + O(∆t),
∂t ∆t
and similar for a backward difference. Then forward,

Ujm+1 − Ujm m
Uj+1 − 2Ujm + Uj−1
m
= ,
∆t ∆x2
and similarly for the others. Then we can discretise the boundary conditions.

For ODEs, forward Euler was easy to implement. This kind of holds here. However, for implicit we have a
system of non-linear equations. Let’s consider forward euler first though,

Ujm+1 = Ujm + µ(Uj+1

m
− 2Ujm + Uj−1
m
)

Thus once we have the initial and boundary data we can find our starting points. Then we have all the Uj1
and then all the Uj2 and so on. We now let,

δx2 Uj = Uj+1 − 2Uj + Uj−1 .

Then the θ-method becomes,

Ujm+1 − Ujm δx2 Ujm+1 δx2 Ujm

=θ + (1 − θ)
∆t ∆x2 ∆x2
and it becomes a massive mess after some rearranging. We are going to set m = 0. We can write,

(I − µθA)U m+1 = (I 0 + µ(1 − θ)A)U m + g m+1

where U m is the vector of the Ujm ’s, then I 0 is the identity minus the top and bottom 1’s are zero, and
g m+1 = (ua (tm+1 ), 0, . . . , 0, un (tm+1 ))T .
If we have a BC of αu + β ∂u ∂t = γ. Then we can just take a forward difference given x0 = a,

u(x1 , t) − u(x0 , t)
+ O(∆x)
∆x
This means,
U1m+1 − U0m+1
αU0m+1 + β =γ
∆x
and then in Explicit Euler, we have something more complicated. It is harder to approximate the solution on
the boundary. We do something similar with the θ-method. Just rearrange and fiddle our equations to no

11 James Arthur
4 First Order PDEs Practical Numerical Analysis

longer have upper line of all zeros. However, this is not ∆x, so we want to get a second order finite difference
scheme for our BC. Hence,
∂u u(x1 , t) − u(x−1 , t)
(a, t) = + O(∆x2 )
∂t 2∆x
Hence we have,
U m+1 − U−1m+1
αU0m+1 + β 1 = γ,
2∆x
m+1
and we rearrange for U−1 and put in when appropriate.

12 James Arthur