IB Numerical Analysis (Full)

Part IB — Numerical Analysis

Based on lectures by G. Moore

Notes taken by Dexter Chua

Lent 2016

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.

Polynomial approximation

Interpolation by polynomials. Divided differences of functions and relations to deriva-

tives. Orthogonal polynomials and their recurrence relations. Least squares approx-

imation by polynomials. Gaussian quadrature formulae. Peano kernel theorem and

applications. [6]

Computation of ordinary differential equations

Euler’s method and proof of convergence. Multistep methods, including order, the root

condition and the concept of convergence. Runge-Kutta schemes. Stiff equations and

A-stability. [5]

Systems of equations and least squares calculations

LU triangular factorization of matrices. Relation to Gaussian elimination. Column

pivoting. Factorizations of symmetric and band matrices. The Newton-Raphson

method for systems of non-linear algebraic equations. QR factorization of rectangular

matrices by Gram-Schmidt, Givens and Householder techniques. Application to linear

least squares calculations. [5]

Contents

0 Introduction

1 Polynomial interpolation

1.1 The interpolation problem

1.2 The Lagrange formula

1.3 The Newton formula

1.4 A useful property of divided differences

1.5 Error bounds for polynomial interpolation

2 Orthogonal polynomials

2.1 Scalar product

2.2 Orthogonal polynomials

2.3 Three-term recurrence relation

2.4 Examples

2.5 Least-squares polynomial approximation

3 Approximation of linear functionals

3.1 Linear functionals

3.2 Gaussian quadrature

4 Expressing errors in terms of derivatives

5 Ordinary differential equations

5.1 Introduction

5.2 One-step methods

5.3 Multi-step methods

5.4 Runge-Kutta methods

6 Stiff equations

6.1 Introduction

6.2 Linear stability

7 Implementation of ODE methods

7.1 Local error estimation

7.2 Solving for implicit methods

8 Numerical linear algebra

8.1 Triangular matrices

8.2 LU factorization

8.3 A = LU for special A

9 Linear least squares

0 Introduction

Numerical analysis is the study of algorithms. There are many problems we

would like algorithms to solve. In this course, we will tackle the problems of

polynomial approximation, solving ODEs and solving linear equations. These

are all problems that frequently arise when we do (applied) maths.

In general, there are two things we are concerned with — accuracy and speed.

Accuracy is particularly important in the cases of polynomial approximation and

solving ODEs, since we are trying approximate things. We would like to make

good approximations of the solution with relatively little work. On the other

hand, in the case of solving linear equations, we are more concerned with speed

— our solutions will be exact (up to numerical errors due to finite precision of

calculations), but we would like to solve it quickly. We might have to deal with

huge systems, and we don’t want the computation time to grow too quickly.

In the past, this was an important subject, since they had no computers.

The algorithms had to be implemented by hand. It was thus very important to

find some practical and efficient methods of computing things, or else it would

take them forever to calculate what they wanted. So we wanted quick algorithms

that give reasonably accurate results.

Nowadays, this is still an important subject. While we have computers that

are much faster at computation, we still want our programs to be fast. We would

also want to get really accurate results, since we might be using them to, say,

send our rocket to the Moon. Moreover, with more computational power, we

might sacrifice efficiency for some other desirable properties. For example, if we

are solving for the trajectory of a particle, we might want the solution to satisfy

the conservation of energy. This would require some much more complicated and

slower algorithms that no one would have considered in the past. Nowadays, with

computers, these algorithms become more feasible, and are becoming increasingly

more popular.

1 Polynomial interpolation

Polynomials are nice. Writing down a polynomial of degree

involves only

+ 1

numbers. They are easy to evaluate, integrate and differentiate. So it would be

nice if we can approximate things with polynomials. For simplicity, we will only

deal with real polynomials.

Notation.

We write

[

] for the real linear vector space of polynomials (with

real coefficients) having degree n or less.

It is easy to show that dim(P

[x]) = n + 1.

1.1 The interpolation problem

The idea of polynomial interpolation is that we are given

+ 1 distinct interpo-

lation points

}

i=0

⊆ R

, and

+ 1 data values

}

i=0

⊆ R

. The objective is

to find a p ∈ P

[x] such that

p(x

) = f

for i = 0, ··· , n.

In other words, we want to fit a polynomial through the points (x

, f

There are many situations where this may come up. For example, we may

be given

+ 1 actual data points, and we want to fit a polynomial through

the points. Alternatively, we might have a complicated function

, and want

to approximate it with a polynomial

such that

and

agree on at least that

n + 1 points.

The naive way of looking at this is that we try a polynomial

p(x) = a

+ a

n−1

+ ··· + a

and then solve the system of equations

= p(x

) = a

+ a

n−1

+ ··· + a

This is a perfectly respectable system of

+ 1 equations in

+ 1 unknowns.

From linear algebra, we know that in general, such a system is not guaranteed

to have a solution, and if the solution exists, it is not guaranteed to be unique.

That was not helpful. So our first goal is to show that in the case of polynomial

interpolation, the solution exists and is unique.

1.2 The Lagrange formula

It turns out the problem is not too hard. You can probably figure it out yourself

if you lock yourself in a room for a few days (or hours). In fact, we will come up

with two solutions to the problem.

The first is via the Lagrange cardinal polynomials. These are simple to

state, and it is obvious that these work very well. However, practically, this is

not the best way to solve the problem, and we will not talk about them much.

Instead, we use this solution as a proof of existence of polynomial interpolations.

We will then develop our second method on the assumption that polynomial

interpolations exist, and find a better way of computing them.

Definition

(Lagrange cardinal polynomials)

The Lagrange cardinal polynomials

with respect to the interpolation points {x

}

i=0

are, for k = 0, ··· , n,

(x) =

i=0,i6=k

x − x

− x

Note that these polynomials have degree exactly

. The significance of these

polynomials is we have

(

) = 0 for

i 6

, and

(

) = 1. In other words, we

have

) = δ

This is obvious from definition.

With these cardinal polynomials, we can immediately write down a solution

to the interpolation problem.

Theorem. The interpolation problem has exactly one solution.

Proof. We define p ∈ P

[x] by

p(x) =

k=0

(x).

Evaluating at x

gives

p(x

) =

k=0

) =

k=0

= f

So we get existence.

For uniqueness, suppose

p, q ∈ P

[

] are solutions. Then the difference

p − q ∈ P

[

] satisfies

(

) = 0 for all

, i.e. it has

+ 1 roots. However, a

non-zero polynomial of degree

can have at most

roots. So in fact

p − q

zero, i.e. p = q.

While this works, it is not ideal. If we one day decide we should add one more

interpolation point, we would have to recompute all the cardinal polynomials,

and that is not fun. Ideally, we would like some way to reuse our previous

computations when we have new interpolation points.

1.3 The Newton formula

The idea of Newton’s formula is as follows — for

= 0

, ··· , n

, we write

∈ P

[

]

for the polynomial that satisfies

) = f

for i = 0, ··· , k.

This is the unique degree-

polynomial that satisfies the first

+ 1 conditions,

whose existence (and uniqueness) is guaranteed by the previous section. Then

we can write

p(x) = p

(x) = p

(x) + (p

(x) − p

(x)) + ··· + (p

(x) − p

n−1

(x)).

Hence we are done if we have an efficient way of finding the differences

−p

k−1

We know that

and

k−1

agree on

, ··· , x

k−1

. So

−p

k−1

evaluates to

0 at those points, and we must have

(x) − p

k−1

(x) = A

k−1

i=0

(x − x

for some A

yet to be found out. Then we can write

p(x) = p

(x) = A

k=1

k−1

i=0

(x − x

This formula has the advantage that it is built up gradually from the interpo-

lation points one-by-one. If we stop the sum at any point, we have obtained

the polynomial that interpolates the data for the first

points (for some

Conversely, if we have a new data point, we just need to add a new term, instead

of re-computing everything.

All that remains is to find the coefficients

. For

= 0, we know

is the

unique constant polynomial that interpolates the point at x

, i.e. A

= f

For the others, we note that in the formula for

− p

k−1

, we find that

the leading coefficient of

. But

k−1

(

) has no degree

term. So

must

be the leading coefficient of

. So we have reduced our problem to finding the

leading coefficients of p

The solution to this is known as the Newton divided differences. We first

invent a new notation:

= f[x

, ··· , x

Note that these coefficients depend only on the first

interpolation points.

Moreover, since the labelling of the points

, ··· , x

is arbitrary, we don’t have

to start with x

. In general, the coefficient

f[x

, ··· , x

]

is the leading coefficient of the unique

q ∈ P

k−j

[

] such that

(

) =

for

i = j, ··· , k.

While we do not have an explicit formula for what these coefficients are, we

can come up with a recurrence relation for these coefficients.

Theorem

(Recurrence relation for Newton divided differences)

For 0

≤ j <

k ≤ n, we have

f[x

, ··· , x

] =

f[x

j+1

, ··· , x

] − f[x

, ··· , x

k−1

]

− x

Proof.

The key to proving this is to relate the interpolating polynomials. Let

, q

∈ P

k−j−1

[x] and q

∈ P

k−j

satisfy

) = f

i = j, ··· , k − 1

) = f

i = j + 1, ··· , k

) = f

i = j, ··· , k

We now claim that

(x) =

x − x

− x

(x) +

− x

(x).

We can check directly that the expression on the right correctly interpolates

the points

for

j, ··· , k

. By uniqueness, the two expressions agree. Since

[

, ··· , x

[

j+1

, ··· , x

] and

[

, ··· , x

k−1

] are the leading coefficients of

, q

respectively, the result follows.

Thus the famous Newton divided difference table can be constructed

f[∗, ∗] f[∗, ∗, ∗] ··· f[∗, ··· , ∗]

f[x

]

f[x

, x

]

f[x

] f[x

, x

]

f[x

, x

]

f[x

] f[x

, x

] ··· f[x

, x

, ··· , x

]

f[x

, x

]

f[x

]

f[x

]

From the first

columns, we can find the

+ 1th column using the recurrence

relation above. The values of

can then be found at the top diagonal, and

this is all we really need. However, to compute this diagonal, we will need to

compute everything in the table.

In practice, we often need not find the actual interpolating polynomial. If

we just want to evaluate

(

ˆx

) at some new point

ˆx

using the divided table, we

can use Horner’s scheme, given by

S <- f[x

,..., x

]

for k = n - 1,..., 0

S <- (^x - x

)S + f[x

,..., x

]

end

This only takes O(n) operations.

If an extra data point

n+1

, f

n+1

}

is added, then we only have to compute

an extra diagonal

[

, ··· , x

n+1

] for

n, ··· ,

0 in the divided difference table

to obtain the new coefficient, and the old results can be reused. This requires

O(n) operations. This is less straightforward for Lagrange’s method.

1.4 A useful property of divided differences

In the next couple of sections, we are interested in the error of polynomial inter-

polation. Suppose the data points come from

(

) for some complicated

f we want to approximate, and we interpolate f at n data points {x

}

i=1

with

. How does the error

(

) =

(

)

− p

(

) depend on

and the choice of

interpolation points?

At the interpolation point, the error is necessarily 0, by definition of interpo-

lation. However, this does not tell us anything about the error elsewhere.

Our ultimate objective is to bound the error

by suitable multiples of the

derivatives of

. We will do this in two steps. We first relate the derivatives to

the divided differences, and then relate the error to the divided differences.

The first part is an easy result based on the following purely calculus lemma.

Lemma.

Let

g ∈ C

[

a, b

] have a continuous

th derivative. Suppose

is zero

at m + ` distinct points. Then g

(m)

has at least ` distinct zeros in [a, b].

Proof.

This is a repeated application of Rolle’s theorem. We know that between

every two zeros of

, there is at least one zero of

∈ C

m−1

[

a, b

]. So by

differentiating once, we have lost at most 1 zeros. So after differentiating

times, g

(m)

has lost at most m zeros. So it still has at least ` zeros.

Theorem.

Let

}

i=0

∈

[

a, b

] and

f ∈ C

[

a, b

]. Then there exists some

ξ ∈ (a, b) such that

f[x

, ··· , x

] =

(n)

(ξ).

Proof.

Consider

f − p

∈ C

[

a, b

]. This has at least

+ 1 distinct zeros in

[

a, b

]. So by the lemma,

(n)

− p

(n)

must vanish at some

ξ ∈

(

a, b

). But

then p

(n)

= n!f[x

, ··· , x

] constantly. So the result follows.

1.5 Error bounds for polynomial interpolation

The actual error bound is not too difficult as well. It turns out the error

f − p

is “like the next term in the Newton’s formula”. This vague sentence

is made precise in the following theorem:

Theorem.

Assume

}

i=0

⊆

[

a, b

] and

f ∈ C

[

a, b

]. Let

¯x ∈

[

a, b

] be a non-

interpolation point. Then

(¯x) = f[x

, x

, ··· , x

, ¯x]ω(¯x),

where

ω(x) =

i=0

(x − x

Note that we forbid the case where

¯x

is an interpolation point, since it is

not clear what the expression

[

, x

, ··· , x

, ¯x

] means. However, if

¯x

is an

interpolation point, then both

(

¯x

) and

(

¯x

) are zero, so there isn’t much to

say.

Proof. We think of ¯x = x

n+1

as a new interpolation point so that

n+1

(x) − p

(x) = f[x

, ··· , x

, ¯x]ω(x)

for all

x ∈ R

. In particular, putting

¯x

, we have

n+1

(

¯x

) =

(

¯x

), and we get

the result.

Combining the two results, we find

Theorem.

If in addition

f ∈ C

n+1

[

a, b

], then for each

x ∈

[

a, b

], we can find

∈ (a, b) such that

(x) =

(n + 1)!

(n+1)

(ξ

)ω(x)

Proof.

The statement is trivial if

is an interpolation point — pick arbitrary

, and both sides are zero. Otherwise, this follows directly from the last two

theorems.

This is an exact result, which is not too useful, since there is no easy

constructive way of finding what

should be. Instead, we usually go for a

bound. We introduce the max norm

kgk

∞

= max

t∈[a,b]

|g(t)|.

This gives the more useful bound

Corollary. For all x ∈ [a, b], we have

|f(x) − p

(x)| ≤

(n + 1)!

(n+1)

∞

|ω(x)|

Assuming our function

is fixed, this error bound depends only on

(

which depends on our choice of interpolation points. So can we minimize

(

)

in some sense by picking some clever interpolation points ∆ =

}

i=0

? Here

we will have

fixed. So instead, we put ∆ as the subscript. We can write our

bound as

kf − p

∆

∞

≤

(n + 1)!

(n+1)

∞

kω

∆

∞

So the objective is to find a ∆ that minimizes kω

∆

∞

For the moment, we focus on the special case where the interval is [

−

1].

The general solution can be obtained by an easy change of variable.

For some magical reasons that hopefully will become clear soon, the optimal

choice of ∆ comes from the Chebyshev polynomials.

Definition

(Chebyshev polynomial)

The Chebyshev polynomial of degree

[−1, 1] is defined by

(x) = cos(nθ),

where x = cos θ with θ ∈ [0, π].

So given an

, we find the unique

that satisfies

cos θ

, and then find

cos

(

nθ

). This is in fact a polynomial in disguise, since from trigonometric

identities, we know

cos

(

nθ

) can be expanded as a polynomial in

cos θ

up to

degree n.

Two key properties of T

on [−1, 1] are

(i) The maximum absolute value is obtained at

= cos



πk



for k = 0, ··· , n with

) = (−1)

(ii) This has n distinct zeros at

= cos



2k − 1



for k = 1, ··· , n.

When plotted out, the polynomials look like this:

(x)

−1 1

−1

All that really matters about the Chebyshev polynomials is that the maximum

is obtained at

+ 1 distinct points with alternating sign. The exact form of the

polynomial is not really important.

Notice there is an intentional clash between the use of

as the zeros and

as the interpolation points — we will show these are indeed the optimal

interpolation points.

We first prove a convenient recurrence relation for the Chebyshev polynomials:

Lemma

(3-term recurrence relation)

The Chebyshev polynomials satisfy the

recurrence relations

n+1

(x) = 2xT

(x) − T

n−1

(x)

with initial conditions

(x) = 1, T

(x) = x.

Proof.

cos((n + 1)θ) + cos((n − 1)θ) = 2 cos θ cos(nθ).

This recurrence relation can be useful for many things, but for our purposes,

we only use it to show that the leading coefficient of T

is 2

n−1

(for n ≥ 1).

Theorem

(Minimal property for

n ≥

On [

−

1], among all polynomials

p ∈ P

[

] with leading coefficient 1,

n−1

minimizes

kpk

∞

. Thus, the

minimum value is

n−1

Proof.

We proceed by contradiction. Suppose there is a polynomial

∈ P

with leading coefficient 1 such that kq

∞

n−1

. Define a new polynomial

r =

n−1

− q

This is, by assumption, non-zero.

Since both the polynomials have leading coefficient 1, the difference must

have degree at most

n −

1, i.e.

r ∈ P

n−1

[

]. Since

n−1

(

) =

n−1

, and

(

)

| <

n−1

by assumption,

alternates in sign between these

+ 1 points.

But then by the intermediate value theorem,

has to have at least

zeros. This

is a contradiction, since r has degree n − 1, and cannot be zero.

Corollary. Consider

∆

i=0

(x − x

) ∈ P

n+1

[x]

for any distinct points ∆ = {x

}

i=0

⊆ [−1, 1]. Then

min

∆

kω

∆

∞

This minimum is achieved by picking the interpolation points to be the zeros of

n+1

, namely

= cos



2k + 1

2n + 2



, k = 0, ··· , n.

Theorem.

For

f ∈ C

n+1

[

−

1], the Chebyshev choice of interpolation points

gives

kf − p

∞

≤

(n + 1)!

(n+1)

∞

Suppose

has as many continuous derivatives as we want. Then as we

increase

, what happens to the error bounds? The coefficients involve dividing

by an exponential and a factorial. Hence as long as the higher derivatives of

don’t blow up too badly, in general, the error will tend to zero as

n → ∞

, which

makes sense.

The last two results can be easily generalized to arbitrary intervals [

a, b

], and

this is left as an exercise for the reader.

2 Orthogonal polynomials

It turns out the Chebyshev polynomials is just an example of a more general

class of polynomials, known as orthogonal polynomials. As in linear algebra, we

can define a scalar product on the space of polynomials, and then find a basis

of orthogonal polynomials of the vector space under this scalar product. We

shall show that each set of orthogonal polynomials has a three-term recurrence

relation, just like the Chebyshev polynomials.

2.1 Scalar product

The scalar products we are interested in would be generalization of the usual

scalar product on Euclidean space,

hx, yi =

i=1

We want to generalize this to vector spaces of functions and polynomials. We

will not provide a formal definition of vector spaces and scalar products on an

abstract vector space. Instead, we will just provide some examples of commonly

used ones.

Example.

(i)

Let

[

a, b

], where [

a, b

] is a finite interval and

s ≥

0. Pick a weight

function

(

)

∈ C

(

a, b

) such that

(

)

0 for all

x ∈

(

a, b

), and

integrable over [

a, b

]. In particular, we allow

to vanish at the end points,

or blow up mildly such that it is still integrable.

We then define the inner product to be

hf, gi =

w(x)f(x)d(x) dx.

(ii)

We can allow [

a, b

] to be infinite, e.g. [0

, ∞

) or even (

−∞, ∞

), but we have

to be more careful. We first define

hf, gi =

w(x)f(x)g(x) dx

as before, but we now need more conditions. We require that

(

)

to exist for all

n ≥

0, since we want to allow polynomials in our vector

space. For example,

(

) =

−x

on [0

, ∞

), works, or

(

) =

−x

(

−∞, ∞

). These are scalar products for

[

] for

n ≥

0, but we cannot

extend this definition to all smooth functions since they might blow up too

fast at infinity. We will not go into the technical details, since we are only

interested in polynomials, and knowing it works for polynomials suffices.

(iii) We can also have a discrete inner product, defined by

hf, gi =

j=1

f(ξ

)g(ξ

)

with

{ξ

}

j=1

distinct points and

}

j=1

0. Now we have to restrict

ourselves a lot. This is a scalar product for

m−1

[

], but not for

higher degrees, since a scalar product should satisfy

hf, f i >

0 for

f 6

= 0.

In particular, we cannot extend this to all smooth functions.

With an inner product, we can define orthogonality.

Definition

(Orthogonalilty)

Given a vector space

and an inner product

h·, ·i, two vectors f, g ∈ V are orthogonal if hf, gi = 0.

2.2 Orthogonal polynomials

Definition

(Orthogonal polynomial)

Given a vector space

of polynomials

and inner product

h·, ·i

, we say

∈ P

[

] is the

th orthogonal polynomial if

, qi = 0 for all q ∈ P

n−1

[x].

In particular, hp

, p

i = 0 for n 6= m.

We said the orthogonal polynomial, but we need to make sure such a polyno-

mial has to be unique. It is clearly not unique, since if

satisfies these relations,

then so does

λp

for all

λ 6

= 0. For uniqueness, we need to impose some scaling.

We usually do so by requiring the leading polynomial to be 1, i.e. it is monic.

Definition

(Monic polynomial)

A polynomial

p ∈ P

[

] is monic if the coeffi-

cient of x

is 1.

In practice, most famous traditional polynomials are not monic. They have

a different scaling imposed. Still, as long as we have some scaling, we will have

uniqueness.

We will not mess with other scalings, and stick to requiring them to be monic

since this is useful for proving things.

Theorem.

Given a vector space

of functions and an inner product

h·, ·i

there exists a unique monic orthogonal polynomial for each degree

n ≥

0. In

addition, {p

}

k=0

form a basis for P

[x].

Proof.

This is a big induction proof over both parts of the theorem. We induct

over

. For the base case, we pick

(

) = 1, which is the only degree-zero monic

polynomial.

Now suppose we already have {p

}

k=0

satisfying the induction hypothesis.

Now pick any monic

n+1

∈ P

n+1

[

], e.g.

n+1

. We now construct

n+1

from

n+1

by the Gram-Schmidt process. We define

n+1

= q

n+1

−

k=0

n+1

, p

This is again monic since q

n+1

is, and we have

n+1

, p

i = 0

for all m ≤ n, and hence hp

n+1

, pi = 0 for all p ∈ P

[x] = hp

, ··· , p

To obtain uniqueness, assume both

n+1

, ˆp

n+1

∈ P

n+1

[

] are both monic

orthogonal polynomials. Then r = p

n+1

− ˆp

n+1

∈ P

[x]. So

hr, ri = hr, p

n+1

− ˆp

n+1

i = hr, p

n+1

i − hr, ˆp

n+1

i = 0 − 0 = 0.

So r = 0. So p

n+1

= ˆp

n−1

Finally, we have to show that

, ··· , p

n+1

form a basis for

n+1

[

]. Now

note that every p ∈ P

n+1

[x] can be written uniquely as

p = cp

n+1

+ q,

where

q ∈ P

[

]. But

}

k=0

is a basis for

[

]. So

can be uniquely

decomposed as a linear combination of p

, ··· , p

Alternatively, this follows from the fact that any set of orthogonal vectors

must be linearly independent, and since there are

+ 2 of these vectors and

n+1

[x] has dimension n + 2, they must be a basis.

In practice, following the proof naively is not the best way of producing the

new

n+1

. Instead, we can reduce a lot of our work by making a clever choice of

n+1

2.3 Three-term recurrence relation

Recall that for the Chebyshev polynomials, we obtained a three-term recurrence

relation for them using special properties of the cosine. It turns out these

recurrence relations exist in general for orthogonal polynomials.

We start by picking

n+1

in the previous proof. We now use the fact

that

hxf, gi = hf, xgi.

This is not necessarily true for arbitrary inner products, but for most sensible

inner products we will meet in this course, this is true. In particular, it is clearly

true for inner products of the form

hf, gi =

w(x)f(x)g(x) dx.

Assuming this, we obtain the following theorem.

Theorem. Monic orthogonal polynomials are generated by

k+1

(x) = (x − α

(x) − β

k−1

(x)

with initial conditions

= 1, p

(x) = (x − α

where

hxp

, p

, β

, p

k−1

, p

k−1

Proof. By inspection, the p

given is monic and satisfies

, p

i = 0.

Using q

n+1

= xp

in the Gram-Schmidt process gives

n+1

= xp

−

k=0

hxp

, p

n+1

= xp

−

k=0

, xp

, p

We notice that

, xp

and vanishes whenever

has degree less than

. So

we are left with

= xp

−

hxp

, p

−

, xp

n−1

, p

n−1

= (x − α

−

, xp

n−1

, p

n−1

Now we notice that

n−1

is a monic polynomial of degree

so we can write

this as xp

n−1

= p

+ q. Thus

, xp

n−1

i = hp

, p

+ qi = hp

, p

Hence the coefficient of p

n−1

is indeed the β we defined.

2.4 Examples

The four famous examples are the Legendre polynomials, Chebyshev polynomials,

Laguerre polynomials and Hermite polynomials. We first look at how the

Chebyshev polynomials fit into this framework.

Chebyshev is based on the scalar product defined by

hf, gi =

−1

√

1 − x

f(x)g(x) dx.

Note that the weight function blows up mildly at the end, but this is fine since

it is still integrable.

This links up with

(x) = cos(nθ)

for x = cos θ via the usual trigonometric substitution. We have

, T

i =

√

1 − cos

cos(nθ) cos(mθ) sin θ dθ

cos(nθ) cos(mθ) dθ

= 0 if m 6= n.

The other orthogonal polynomials come from scalar products of the form

hf, gi =

w(x)f(x)g(x) dx,

as described in the table below:

Type Range Weight

Legendre [−1, 1] 1

Chebyshev [−1, 1]

√

1−x

Laguerre [0, ∞) e

−x

Hermite (−∞, ∞) e

−x

2.5 Least-squares polynomial approximation

If we want to approximate a function with a polynomial, polynomial interpolation

might not be the best idea, since all we do is make sure the polynomial agrees

with

at certain points, and then hope it is a good approximation elsewhere.

Instead, the idea is to choose a polynomial

[

] that “minimizes the error”.

What exactly do we mean by minimizing the error? The error is defined as

the function

f −p

. So given an appropriate inner product on the vector space of

continuous functions, we want to minimize

kf − pk

= hf − p, f − pi.

This is usually of the form

hf − p, f − pi =

w(x)[f(x) − p(x)]

dx,

but we can also use alternative inner products such as

hf − p, f − pi =

j=1

[f(ξ

) − p(ξ

)]

Unlike polynomial interpolation, there is no guarantee that the approximation

agrees with the function anywhere. Unlike polynomial interpolation, there is

some guarantee that the total error is small (or as small as we can get, by

definition). In particular, if

is continuous, then the Weierstrass approximation

theorem tells us the total error must eventually vanish.

Unsurprisingly, the solution involves the use of the orthogonal polynomials

with respect to the corresponding inner products.

Theorem.

}

k=0

are orthogonal polynomials with respect to

h·, ·i

, then

the choice of c

such that

p =

k=0

minimizes kf − pk

is given by

hf, p

and the formula for the error is

kf − pk

= kfk

−

k=0

hf, p

Note that the solution decouples, in the sense that

depends only on

and

. If we want to take one more term, we just need to compute an extra term,

and not redo our previous work.

Also, we notice that the formula for the error is a positive term

kfk

sub-

tracting a lot of squares. As we increase

, we subtract more squares, and the

error decreases. If we are lucky, the error tends to 0 as we take

n → ∞

. Even

though we might not know how many terms we need in order to get the error to

be sufficiently small, we can just keep adding terms until the computed error

small enough (which is something we have to do anyway even if we knew what

n to take).

Proof. We consider a general polynomial

p =

k=0

We substitute this in to obtain

hf − p, f − pi = hf, fi − 2

k=0

hf, p

i +

k=0

Note that there are no cross terms between the different coefficients. We minimize

this quadratic by setting the partial derivatives to zero:

0 =

∂

∂c

hf − p, f − pi = −2hf, p

i + 2c

To check this is indeed a minimum, note that the Hessian matrix is simply 2

which is positive definite. So this is really a minimum. So we get the formula for

the

’s as claimed, and putting the formula for

gives the error formula.

Note that our constructed

p ∈ P

[

] has a nice property: for

k ≤ n

, we have

hf − p, p

i = hf, p

i − hp, p

i = hf, p

i −

hf, p

, p

i = 0.

Thus for all q ∈ P

[x], we have

hf − p, qi = 0.

In particular, this is true when

, and tells us

hf, pi

hp, pi

. Using this to

expand hf − p, f − pi gives

kf − pk

+ kpk

= kfk

which is just a glorified Pythagoras theorem.

3 Approximation of linear functionals

3.1 Linear functionals

In this chapter, we are going to study approximations of linear functions. Before

we start, it is helpful to define what a linear functional is, and look at certain

examples of these.

Definition

(Linear functional)

A linear functional is a linear mapping

V →

R, where V is a real vector space of functions.

In generally, a linear functional is a linear mapping from a vector space to

its underlying field of scalars, but for the purposes of this course, we will restrict

to this special case.

We usually don’t put so much emphasis on the actual vector space

. Instead,

we provide a formula for

, and take

to be the vector space of functions for

which the formula makes sense.

Example.

(i) We can choose some fixed ξ ∈ R, and define a linear functional by

L(f) = f(ξ).

(ii) Alternatively, for fixed η ∈ R we can define our functional by

L(f) = f

(η).

In this case, we need to pick a vector space in which this makes sense, e.g.

the space of continuously differentiable functions.

(iii) We can define

L(f) =

f(x) dx.

The set of continuous (or even just integrable) functions defined on [

a, b

]

will be a sensible domain for this linear functional.

(iv)

Any linear combination of these linear functions are also linear functionals.

For example, we can pick some fixed α, β ∈ R, and define

L(f) = f(β) − f(α) −

β − α

(β) + f

(α)].

The objective of this chapter is to construct approximations to more compli-

cated linear functionals (usually integrals, possibly derivatives point values) in

terms of simpler linear functionals (usually point values of f itself).

For example, we might produce an approximation of the form

L(f) ≈

i=0

f(x

where V = C

[a, b], p ≥ 0, and {x

}

i=0

⊆ [a, b] are distinct points.

How can we choose the coefficients

and the points

so that our approxi-

mation is “good”?

We notice that most of our functionals can be easily evaluated exactly when

is a polynomial. So we might approximate our function

by a polynomial,

and then do it exactly for polynomials.

More precisely, we let

}

i=0

⊆

[

a, b

] be arbitrary points. Then using the

Lagrange cardinal polynomials `

, we have

f(x) ≈

i=0

f(x

(x).

Then using linearity, we can approximate

L(f) ≈ L

i=0

f(x

(x)

i=0

L(`

)f(x

So we can pick

= L(`

Similar to polynomial interpolation, this formula is exact for

f ∈ P

[

]. But we

could do better. If we can freely choose

}

i=0

and

}

i=0

, then since we now

have 2

+ 2 free parameters, we might expect to find an approximation that is

exact for

f ∈ P

2N+1

[

]. This is not always possible, but there are cases when

we can. The most famous example is Gaussian quadrature.

3.2 Gaussian quadrature

The objective of Gaussian quadrature is to approximate integrals of the form

L(f) =

w(x)f(x) dx,

where w(x) is a weight function that determines a scalar product.

Traditionally, we have a different set of notations used for Gaussian quadra-

ture. So just in this section, we will use some funny notation that is inconsistent

with the rest of the course.

We let

hf, gi =

w(x)f(x)g(x) dx

be a scalar product for

[

]. We will show that we can find weights, written

}

k=1

, and nodes, written {c

}

k=1

⊆ [a, b], such that the approximation

w(x)f(x) dx ≈

k=1

f(c

)

is exact for

f ∈ P

2ν−1

[

]. The nodes

}

k=1

will turn out to be the zeros of

the orthogonal polynomial

with respect to the scalar product. The aim of

this section is to work this thing out.

We start by showing that this is the best we can achieve.

Proposition.

There is no choice of

weights and nodes such that the approxi-

mation of

w(x)f(x) dx is exact for all f ∈ P

2ν

[x].

Proof. Define

q(x) =

k=1

(x − c

) ∈ P

[x].

Then we know

w(x)q

(x) dx > 0,

since q

is always non-negative and has finitely many zeros. However,

k=1

) = 0.

So this cannot be exact for q

Recall that we initially had the idea of doing the approximation by interpo-

lating f at some arbitrary points in [a, b]. We call this ordinary quadrature.

Theorem

(Ordinary quadrature)

For any distinct

}

k=1

⊆

[

a, b

], let

}

k=1

be the Lagrange cardinal polynomials with respect to

}

k=1

. Then by choosing

w(x)`

(x) dx,

the approximation

L(f) =

w(x)f(x) dx ≈

k=1

f(c

)

is exact for f ∈ P

ν−1

[x].

We call this method ordinary quadrature.

This simple idea is how we generate many classical numerical integration

techniques, such as the trapezoidal rule. But those are quite inaccurate. It turns

out a clever choice of

}

does much better — take them to be the zeros of

the orthogonal polynomials. However, to do this, we must make sure the roots

indeed lie in [

a, b

]. This is what we will prove now — given any inner product,

the roots of the orthogonal polynomials must lie in [a, b].

Theorem.

For

ν ≥

1, the zeros of the orthogonal polynomial

are real, distinct

and lie in (a, b).

We have in fact proved this for a particular case in IB Methods, and the

same argument applies.

Proof.

First we show there is at least one root. Notice that

= 1. Thus for

ν ≥ 1, by orthogonality, we know

w(x)p

(x)p

(x) dx =

w(x)p

(x) dx = 0.

So there is at least one sign change in (

a, b

). We have already got the result we

need for ν = 1, since we only need one zero in (a, b).

Now for

ν >

1, suppose

{ξ

}

j=1

are the places where the sign of

changes

in (a, b) (which is a subset of the roots of p

). We define

q(x) =

j=1

(x − ξ

) ∈ P

[x].

Since this changes sign at the same place as

, we know

maintains the same

hq, p

i =

w(x)q(x)p

(x) dx = 0,

which is impossible, since

does not change sign. Hence we must have

m = ν.

Theorem.

In the ordinary quadrature, if we pick

}

k=1

to be the roots of

(

), then get we exactness for

f ∈ P

2ν−1

[

]. In addition,

}

k=1

are all

positive.

Proof. Let f ∈ P

2ν−1

[x]. Then by polynomial division, we get

f = qp

+ r,

where

q, r

are polynomials of degree at most

ν −

1. We apply orthogonality to

get

w(x)f(x) dx =

w(x)(q(x)p

(x) + r(x)) dx =

w(x)r(x) dx.

Also, since each c

is a root of p

, we get

k=1

f(c

) =

k=1

(q(c

) + r(c

)) =

k=1

r(c

But

has degree at most

ν −

1, and this formula is exact for polynomials in

ν−1

[x]. Hence we know

w(x)f(x) dx =

w(x)r(x) dx =

k=1

r(c

) =

k=1

f(c

To show the weights are positive, we simply pick as special

. Consider

f ∈

}

k=1

⊆ P

2ν−2

[

], for

the Lagrange cardinal polynomials for

}

k=1

. Since

the quadrature is exact for these, we get

0 <

w(x)`

(x) dx =

j=1

) =

j=1

= b

Since this is true for all k = 1, ··· , ν, we get the desired result.

4 Expressing errors in terms of derivatives

As before, we approximate a linear functional L by

L(f) ≈

i=0

(f),

where

are some simpler linear functionals, and suppose this is exact for

f ∈ P

[x] for some k ≥ 0.

Hence we know the error

(f) = L(f) −

i=0

(f) = 0

whenever

f ∈ P

[

]. We say the error annihilates for polynomials of degree less

than k.

How can we use this property to generate formulae for the error and error

bounds? We first start with a rather simple example.

Example. Let L(f) = f(β). We decide to be silly and approximate L(f) by

L(f) ≈ f(α) +

β − α

(β) + f

(α)),

where α 6= β. This is clearly much easier to evaluate. The error is given by

(f) = f(β) − f(α) −

β − α

(β) + f

(α)),

and this vanishes for f ∈ P

[x].

How can we get a more useful error formula? We can’t just use the fact that

it annihilates polynomials of degree

. We need to introduce something beyond

this — the k + 1th derivative. We now assume f ∈ C

k+1

[a, b].

Note that so far, everything we’ve done works if the interval is infinite, as long

as the weight function vanishes sufficiently quickly as we go far away. However,

for this little bit, we will need to require [

a, b

] to be finite, since we want to make

sure we can take the supremum of our functions.

We now seek an exact error formula in terms of

(k+1)

, and bounds of the

form

(f)| ≤ c

(k+1)

∞

for some constant

. Moreover, we want to make

as small as possible. We

don’t want to give a constant of 10 million, while in reality we can just use 2.

Definition

(Sharp error bound)

The constant

is said to be sharp if for any

ε > 0, there is some f

∈ C

k+1

[a, b] such that

(f)| ≥ (c

− ε)kf

(k+1)

∞

This makes it precise what we mean by

“cannot be better”. This doesn’t

say anything about whether

can actually be achieved. This depends on the

particular form of the question.

To proceed, we need Taylor’s theorem with the integral remainder, i.e.

f(x) = f(a) + (x −a)f

(a) + ···+

(x − a)

(k)

(a) +

(x −θ)

(k+1)

(θ) dθ.

This is not really good, since there is an

in the upper limit of the integral.

Instead, we write the integral as

(x − θ)

(k+1)

(θ) dθ,

where we define (x − θ)

is a new function defined by

(x − θ)

(

(x − θ)

x ≥ θ

0 x < θ.

Then if λ is a linear functional that annihilates P

[x], then we have

λ(f) = λ

(x − θ)

(k+1)

(θ) dθ

for all f ∈ C

k+1

[a, b].

For our linear functionals, we can simplify by taking the

inside the integral

sign and obtain

λ(f) =

λ((x − θ)

(k+1)

(θ) dθ,

noting that

acts on (

x − θ

)

∈ C

k−1

[

a, b

] as a function of

, and think of

being held constant.

Of course, pure mathematicians will come up with linear functionals for

which we cannot move the

inside, but for our linear functionals (point values,

derivative point values, integrals etc.), this is valid, as we can verify directly.

Hence we arrive at

Theorem

(Peano kernel theorem)

annihilates polynomials of degree

less, then

λ(f) =

K(θ)f

(k+1)

(θ) dθ

for all f ∈ C

k+1

[a, b], where

Definition (Peano kernel). The Peano kernel is

K(θ) = λ((x − θ)

The important thing is that the kernel

is independent of

. Taking suprema

in different ways, we obtain different forms of bounds:

|λ(f)| ≤











|K(θ)| dθkf

(k+1)

∞

|K(θ)|

dθ

(k+1)

kK(θ)k

∞

(k+1)

Hence we can find the constant

for different choices of the norm. When

computing c

, don’t forget the factor of

By fiddling with functions a bit, we can show these bounds are indeed sharp.

Example. Consider our previous example where

(f) = f(β) − f(α) −

β − α

(β) + f

(α)),

with exactness up to polynomials of degree 2. We wlog assume α < β. Then

K(θ) = e

((x − θ)

) = (β − θ)

− (α − θ)

− (β − α)((β − θ)

+ (α − θ)

Hence we get

K(θ) =











0 a ≤ θ ≤ α

(α − θ)(β − θ) α ≤ θ ≤ β

0 β ≤ θ ≤ b.

Hence we know

(f) =

(α − θ)(β − θ)f

000

(θ) dθ

for all f ∈ C

[a, b].

Note that in this particular case, our function

(

) does not change sign on

[a, b]. Under this extra assumption, we can say a bit more.

First, we note that the bound

|λ(f)| ≤



K(θ) dθ



(k+1)

∞

can be achieved by

k+1

, since this has constant

+ 1th derivative. Also, we

can use the integral mean value theorem to get the bound

λ(f) =

K(θ) dθ

(k+1)

(ξ),

where ξ ∈ (a, b) depends on f. These are occasionally useful.

Example.

Continuing our previous example, we see that

(

)

≤

0 on [

a, b

and

K(θ) dθ = −

(β − α)

Hence we have the bound

(f)| ≤

(β − α)

000

∞

and this bound is achieved for x

. We also have

(f) = −

(β − α)

000

(ξ)

for some f-dependent value of some ξ ∈ (a, b).

Finally, note that Peano’s kernel theorem says if

(

) = 0 for all

f ∈ P

[

then we have

(f) =

K(θ)f

(k+1)

(θ) dθ

for all f ∈ C

k+1

[a, b].

But for any other fixed

= 0

, ··· , k −

1, we also have

(

) = 0 for all

f ∈ P

[x]. So we also know

(f) =

(θ)f

(j+1)

(θ) dθ

for all f ∈ C

j+1

[a, b]. Note that we have a different kernel.

In general, this might not be a good idea, since we are throwing information

away. Yet, this can be helpful if we get some less smooth functions that don’t

have k + 1 derivatives.

5 Ordinary differential equations

5.1 Introduction

Our next big goal is to solve ordinary differential equations numerically. We will

focus on differential equations of the form

(t) = f(t, y(t))

for 0 ≤ t ≤ T , with initial conditions

y(0) = y

The data we are provided is the function

R×R

→ R

, the ending time

T >

and the initial condition

∈ R

. What we seek is the function

: [0

, T

]

→ R

When solving the differential equation numerically, our goal would be to

make our numerical solution as close to the true solution as possible. This makes

sense only if a “true” solution actually exists, and is unique. From IB Analysis

II, we know a unique solution to the ODE exists if f is Lipschitz.

Definition

(Lipschitz function)

A function

R × R

→ R

is Lipschitz with

Lipschitz constant λ ≥ 0 if

kf(t, x) − f (t,

x)k ≤ λkx −

for all t ∈ [0, T ] and x,

x ∈ R

A function is Lipschitz if it is Lipschitz with Lipschitz constant

for some

It doesn’t really matter what norm we pick. It will just change the

. The

importance is the existence of a λ.

A special case is when

= 0, i.e.

does not depend on

. In this case, this

is just an integration problem, and is usually easy. This is a convenient test case

— if our numerical approximation does not even work for these easy problems,

then it’s pretty useless.

Being Lipschitz is sufficient for existence and uniqueness of a solution to

the differential equation, and hence we can ask if our solution converges to

this unique solution. An extra assumption we will often make is that

can

be expanded in a Taylor series to as many degrees as we want, since this is

convenient for our analysis.

What exactly does a numerical solution to the ODE consist of? We first

choose a small time step h > 0, and then construct approximations

≈ y(t

), n = 1, 2, ··· ,

with

. In particular,

− t

n−1

and is always constant. In practice,

we don’t fix the step size

− t

n−1

, and allow it to vary in each step. However,

this makes the analysis much more complicated, and we will not consider varying

time steps in this course.

If we make

smaller, then we will (probably) make better approximations.

However, this is more computationally demanding. So we want to study the

behaviour of numerical methods in order to figure out what h we should pick.

5.2 One-step methods

There are many ways we can classify numerical methods. One important classifi-

cation is one-step versus multi-step methods. In one-step methods, the value of

n+1

depends only on the previous iteration

and

. In multi-step methods,

we are allowed to look back further in time and use further results.

Definition

((Explicit) one-step method)

A numerical method is (explicit)

one-step if y

n+1

depends only on t

and y

, i.e.

n+1

= φ

, y

)

for some function φ

: R × R

→ R

We will later see what “explicit” means.

The simplest one-step method one can imagine is Euler’s method.

Definition (Euler’s method). Euler’s method uses the formula

n+1

= y

+ hf(t

, y

We want to show that this method “converges”. First of all, we need to make

precise the notion of “convergence”. The Lipschitz condition means there is a

unique solution to the differential equation. So we would want the numerical

solution to be able to approximate the actual solution to arbitrary accuracy as

long as we take a small enough h.

Definition

(Convergence of numerical method)

For each

h >

0, we can produce

a sequence of discrete values

for

= 0

, ··· ,

[

T/h

], where [

T/h

] is the integer

part of

T/h

. A method converges if, as

h →

0 and

nh → t

(hence

n → ∞

), we

get

→ y(t),

where

is the true solution to the differential equation. Moreover, we require

the convergence to be uniform in t.

We now prove that Euler’s method converges. We will only do this properly for

Euler’s method, since the algebra quickly becomes tedious and incomprehensible.

However, the proof strategy is sufficiently general that it can be adapted to most

other methods.

Theorem (Convergence of Euler’s method).

(i) For all t ∈ [0, T ], we have

lim

h→0

nh→t

− y(t) = 0.

(ii) Let λ be the Lipschitz constant of f. Then there exists a c ≥ 0 such that

k ≤ ch

λT

− 1

for all 0 ≤ n ≤ [T/h], where e

= y

− y(t

Note that the bound in the second part is uniform. So this immediately gives

the first part of the theorem.

Proof.

There are two parts to proving this. We first look at the local truncation

error. This is the error we would get at each step assuming we got the previous

steps right. More precisely, we write

y(t

n+1

) = y(t

) + h(f, t

, y(t

)) + R

and

is the local truncation error. For the Euler’s method, it is easy to get

, since

(

, y

(

)) =

(

), by definition. So this is just the Taylor series

expansion of

. We can write

as the integral remainder of the Taylor series,

n+1

− θ)y

(θ) dθ.

By some careful analysis, we get

∞

≤ ch

where

c =

∞

This is the easy part, and tends to go rather smoothly even for more complicated

methods.

Once we have bounded the local truncation error, we patch them together to

get the actual error. We can write

n+1

= y

n+1

− y(t

n+1

)

= y

+ hf(t

, y

) −



y(t

) + hf(t

, y(t

)) + R



= (y

− y(t

)) + h



f(t

, y

) − f(t

, y(t

))



− R

Taking the infinity norm, we get

n+1

∞

≤ ky

− y(t

∞

+ hkf(t

, y

) − f(t

, y(t

))k

∞

+ kR

∞

≤ ke

∞

+ hλke

∞

+ ch

= (1 + λh)ke

∞

+ ch

This is valid for all

n ≥

0. We also know

= 0. Doing some algebra, we get

∞

≤ ch

n−1

j=0

(1 + hλ)

≤

((1 + hλ)

− 1) .

Finally, we have

(1 + hλ) ≤ e

λh

since 1 +

λh

is the first two terms of the Taylor series, and the other terms are

positive. So

(1 + hλ)

≤ e

λhn

≤ e

λT

So we obtain the bound

∞

≤ ch

λT

− 1

Then this tends to 0 as we take h → 0. So the method converges.

This works as long as

λ 6

= 0. However,

= 0 is the easy case, since it is

just integration. We can either check this case directly, or use the fact that

λT

−1

→ T as λ → 0.

The same proof strategy works for most numerical methods, but the algebra

will be much messier. We will not do those in full. We, however, take note of

some useful terminology:

Definition

(Local truncation error)

For a general (multi-step) numerical

method

n+1

= φ(t

, y

, ··· , y

the local truncation error is

n+1

= y(t

n+1

) − φ

, y(t

), y(t

), ··· , y(t

)).

This is the error we will make at the (

+ 1)th step if we had accurate values

for the first n steps.

For Euler’s method, the local truncation error is just the Taylor series

remainder term.

Definition

(Order)

The order of a numerical method is the largest

p ≥

1 such

that η

n+1

= O(h

p+1

The Euler method has order 1. Notice that this is one less than the power

of the local truncation error, since when we look at the global error, we drop a

power, and only have e

∼ h.

Let’s try to get a little bit beyond Euler’s method.

Definition (θ-method). For θ ∈ [0, 1], the θ-method is

n+1

= y

+ h



θf(t

, y

) + (1 − θ)f(t

n+1

, y

n+1

)



If we put

= 1, then we just get Euler’s method. The other two most

common choices of θ are θ = 0 (backward Euler) and θ =

(trapezoidal rule).

Note that for

θ 6

= 1, we get an implicit method. This is since

n+1

doesn’t

just appear simply on the left hand side of the equality. Our formula for

n+1

involves

n+1

itself! This means, in general, unlike the Euler method, we can’t

just write down the value of

n+1

given the value of

. Instead, we have to

treat the formula as

(in general) non-linear equations, and solve them to find

n+1

In the past, people did not like to use this, because they didn’t have computers,

or computers were too slow. It is tedious to have to solve these equations in

every step of the method. Nowadays, these are becoming more and more popular

because it is getting easier to solve equations, and

-methods have some huge

theoretical advantages (which we do not have time to get into).

We now look at the error of the θ-method. We have

η = y(t

n+1

) − y(t

) − h



θy

) + (1 − θ)y

n+1

)



We expand all terms about t

with Taylor’s series to obtain



θ −



) +



θ −



000

) + O(h

We see that

gives us an order 2 method. Otherwise, we get an order 1

method.

5.3 Multi-step methods

We can try to make our methods more efficient by making use of previous values

instead of the most recent one. One common method is the AB2 method:

Definition

(2-step Adams-Bashforth method)

The 2-step Adams-Bashforth

(AB2) method has

n+2

= y

n+1

h (3f (t

n+1

, y

n+1

) − f(t

, y

)) .

This is a particular example of the general class of Adams-Bashforth methods.

In general, a multi-step method is defined as follows:

Definition (Multi-step method). A s-step numerical method is given by

`=0

n+`

= h

`=0

f(t

n+`

, y

n+`

This formula is used to find the value of y

n+s

given the others.

One point to note is that we get the same method if we multiply all the

constants

, σ

by a non-zero constant. By convention, we normalize this by

setting ρ

= 1. Then we can alternatively write this as

n+s

= h

`=0

f(t

n+`

, y

n+`

) −

s−1

`=0

n+`

This method is an implicit method if σ

6= 0. Otherwise, it is explicit.

Note that this method is linear in the sense that the coefficients

and

appear linearly, outside the

s. Later we will see more complicated numerical

methods where these coefficients appear inside the arguments of f .

For multi-step methods, we have a slight problem to solve. In a one-step

method, we are given y

, and this allows us to immediately apply the one-step

method to get higher values of

. However, for an

-step method, we need

to use other (possibly 1-step) method to obtain

, ··· , y

s−1

before we can get

started.

Fortunately, we only need to apply the one-step method a fixed, small number

of times, even as

h →

0. So the accuracy of the one-step method at the start

does not matter too much.

We now study the properties of a general multi-step method. The first thing

we can talk about is the order:

Theorem. An s-step method has order p (p ≥ 1) if and only if

`=0

= 0

and

`=0

= k

`=0

k−1

for k = 1, ··· , p, where 0

= 1.

This is just a rather technical result coming directly from definition.

Proof. The local truncation error is

`=0

y(t

n+`

) − h

`=0

n+`

We now expand the y and y

about t

, and obtain

`=0

y(t

) +

∞

k=1

`=0

− k

`=0

k−1

(k)

This is O(h

p+1

) under the given conditions.

Example

(AB2)

In the two-step Adams-Bashforth method, We see that the

conditions hold for p = 2 but not p = 3. So the order is 2.

Instead of dealing with the coefficients

and

directly, it is often convenient

to pack them together into two polynomials. This uses two polynomials associated

with the numerical method. Unsurprisingly, they are

ρ(w) =

`=0

, σ(w) =

`=0

We can then use this to restate the above theorem.

Theorem. A multi-step method has order p (with p ≥ 1) if and only if

ρ(e

) − xσ(e

) = O(x

p+1

)

as x → 0.

Proof. We expand

ρ(e

) − xσ(e

) =

`=0

− x

`=0

We now expand the e

in Taylor series about x = 0. This comes out as

`=0

∞

k=1

`=0

− k

`=0

k−1

So the result follows.

Note that

`=0

= 0, which is the condition required for the method to

even have an order at all, can be expressed as ρ(1) = 0.

Example (AB2). In the two-step Adams-Bashforth method, we get

ρ(w) = w

− w, σ(w) =

w −

We can immediately check that ρ(1) = 0. We also have

ρ(e

) − xσ(e

) =

+ O(x

So the order is 2.

We’ve sorted out the order of multi-step methods. The next thing to check

is convergence. This is where the difference between one-step and multi-step

methods come in. For one-step methods, we only needed the order to understand

convergence. It is a fact that a one step method converges whenever it has an

order p ≥ 1. For multi-step methods, we need an extra condition.

Definition

(Root condition)

We say

(

) satisfies the root condition if all its

zeros are bounded by 1 in size, i.e. all roots

satisfy

|w| ≤

1. Moreover any

zero with |w| = 1 must be simple.

We can imagine this as saying large roots are bad — they cannot get past 1,

and we cannot have too many with modulus 1.

We saw any sensible multi-step method must have

(1) = 0. So in particular,

1 must be a simple zero.

Theorem

(Dahlquist equivalence theorem)

A multi-step method is convergent

if and only if

(i) The order p is at least 1; and

(ii) The root condition holds.

The proof is too difficult to include in this course, or even the Part II version.

This is only done in Part III.

Example

(AB2)

Again consider the two-step Adams-Bashforth method. We

have seen it has order

= 2

≥

1. So we need to check the root condition. So

ρ(w) = w

− w = w(w − 1). So it satisfies the root condition.

Let’s now come up with a sensible strategy for constructing convergent

-step

methods:

(i) Choose a ρ so that ρ(1) = 0 and the root condition holds.

(ii) Choose σ to maximize the order, i.e.

σ =

ρ(w)

log w

(

O(|w − 1|

s+1

) if implicit

O(|w − 1|

) if explicit

We have the two different conditions since for implicit methods, we have

one more coefficient to fiddle with, so we can get a higher order.

Where does the

log w

come from? We try to substitute

(noting that

− 1 ∼ x). Then the formula says

σ(e

) =

ρ(e

) +

(

O(x

s+1

) if implicit

O(x

) if explicit

Rearranging gives

ρ(e

) − xσ(e

) =

(

O(x

s+2

) if implicit

O(x

s+1

) if explicit

which is our order condition. So given any

, there is only one sensible way to

pick σ. So the key is in picking a good enough ρ.

The root conditions is “best” satisfied if

(

) =

s−1

(

w −

1), i.e. all but one

roots are 0. Then we have

n+s

− y

n+s−1

= h

`=0

f(t

n+`

, y

n+`

where α is chosen to maximize order.

Definition

(Adams method)

An Adams method is a multi-step numerical

method with ρ(w) = w

s−1

(w − 1).

These can be either explicit or implicit. In different cases, we get different

names.

Definition

(Adams-Bashforth method)

An Adams-Bashforth method is an

explicit Adams method.

Definition

(Adams-Moulton method)

An Adams-Moulton method is an im-

plicit Adams method.

Example.

We look at the two-step third-order Adams-Moulton method. This

is given by

n+2

− y

n+1

= h



−

f(t

, y

) +

f(t

n+1

, y

n+1

) +

f(t

n+1

, y

n+2

)



Where do these coefficients come from? We have to first expand our

ρ(w)

log w

about

w − 1:

ρ(w)

log w

w(w − 1)

log w

= 1 +

(w − 1) +

(w − 1)

+ O(|w − 1|

These aren’t our coefficients of

, since what we need to do is to rearrange the

first three terms to be expressed in terms of w. So we have

ρ(w)

log w

= −

w +

+ O(|w − 1|

Another important class of multi-step methods is constructed in the opposite

way — we choose a particular σ, and then find the most optimal ρ.

Definition

(Backward differentiation method)

A backward differentiation

method has σ(w) = σ

, i.e.

`=0

n+`

= σ

f(t

n+s

, y

n+s

This is a generalization of the one-step backwards Euler method.

Given this

, we need to choose the appropriate

. Fortunately, this can be

done quite easily.

Lemma.

-step backward differentiation method of order

is obtained by

choosing

ρ(w) = σ

`=1

s−`

(w − 1)

with σ

chosen such that ρ

= 1, namely

`=1

−1

Proof. We need to construct ρ so that

ρ(w) = σ

log w + O(|w − 1|

s+1

This is easy, if we write

log w = −log





= −log



1 −

w − 1



∞

`=1



w − 1



Multiplying by σ

gives the desired result.

For this method to be convergent, we need to make sure it does satisfy the

root condition. It turns out the root condition is satisfied only for

s ≤

6. This is

not obvious by first sight, but we can certainly verify this manually.

5.4 Runge-Kutta methods

Finally, we look at Runge-Kutta methods. These methods are very complicated,

and are rather tedious to analyse. They have been largely ignored for a long

time, until more powerful computers came along and made these methods much

more practical. These are used quite a lot nowadays since they have many nice

properties.

Runge-Kutta methods can be motivated by Gaussian quadrature, but we

will not go into that connection here. Instead, we’ll go straight in and work with

the method.

Definition

(Runge-Kutta method)

General (implicit)

-stage Runge-Kutta

(RK) methods have the form

n+1

= y

+ h

`=1

where

= f





+ c

h, y

+ h

j=1





for ` = 1, ··· , ν.

There are a lot of parameters we have to choose. We need to pick

}

`=1

, {c}

`=1

, {a

}

`,j=1

Note that in general,

}

`=1

have to be solved for, since they are defined in

terms of one another. However, for certain choices of parameters, we can make

this an explicit method. This makes it easier to compute, but we would have

lost some accuracy and flexibility.

Unlike all the other methods we’ve seen so far, the parameters appear inside

. They appear non-linearly inside the functions. This makes the method much

more complicated and difficult to analyse using Taylor series. Yet, once we

manage to do this properly, these have lots of nice properties. Unfortunately, we

will not have time to go into what these properties actually are.

Notice this is a one-step method. So once we get order

p ≥

1, we will have

convergence. So what conditions do we need for a decent order?

This is in general very complicated. However, we can quickly obtain some

necessary conditions. We can consider the case where

is a constant. Then

is always that constant. So we must have

`=1

= 1.

It turns out we also need, for ` = 1, ··· , ν,

j=1

While these are necessary conditions, they are not sufficient. We need other

conditions as well, which we shall get to later. It is a fact that the best possible

order of a ν-stage Runge-Kutta method is 2ν.

To describe a Runge-Kutta method, a standard notation is to put the

coefficients in the Butcher table:

··· a

1ν

ν1

··· a

νν

··· b

We sometimes more concisely write it as

c A

This table allows for a general implicit method. Initially, explicit methods came

out first, since they are much easier to compute. In this case, the matrix

strictly lower triangular, i.e. a

whenever ` ≤ j.

Example.

The most famous explicit Runge-Kutta method is the 4-stage 4th

order one, often called the classical Runge-Kutta method. The formula can be

given explicitly by

n+1

= y

+ 2k

+ k

where

= f(x

, y

)

= f



h, y



= f



h, y



= f (x

+ h, y

+ hk

) .

We see that this is an explicit method. We don’t need to solve any equations.

Choosing the parameters for the Runge-Kutta method to maximize order

is hard. Consider the simplest case, the 2-stage explicit method. The general

formula is

n+1

= y

+ h(b

+ b

where

= f(x

, y

)

= f(x

+ c

h, y

+ r

To analyse this, we insert the true solution into the method. First, we need to

insert the true solution of the ODE into the k’s. We get

= y

)

= f(t

+ c

h, y(t

) + c

))

= y

) + c



∂f

∂t

.y(t

)) + ∇f(t

, y(t

))y

)



+ O(h

)

Fortunately, we notice that the thing inside the huge brackets is just

(

). So

this is

= y

) + c

) + O(h

Hence, our local truncation method for the Runge-Kutta method is

y(t

n+1

) − y(t

) − h(b

+ b

)

= (1 − b

− b

)hy

) +



− b



) + O(h

Now we see why Runge-Kutta methods are hard to analyse. The coefficients

appear non-linearly in this expression. It is still solvable in this case, in the

obvious way, but for higher stage methods, this becomes much more complicated.

In this case, we have a 1-parameter family of order 2 methods, satisfying

+ b

= 1, b

It is easy to check using a simple equation

λh

that it is not possible to get a

higher order method. So as long as our choice of

and

satisfy this equation,

we get a decent order 2 method.

As we have seen, Runge-Kutta methods are really complicated, even in the

simplest case. However, they have too many good properties that they are

becoming very popular nowadays.

6 Stiff equations

6.1 Introduction

Initially, when people were developing numerical methods, people focused mostly

on quantitative properties like order and accuracy. We then develop many

different methods like multi-step methods and Runge-Kutta methods.

More recently, people started to look at structural properties. Often, equations

come with some special properties. For example, a differential equation describing

the motion of a particle would most probably conserve energy. When we

approximate it numerically, we would like the numerical approximation to satisfy

conservation of energy as well. This is what recent developments are looking at —

we want to look at whether numerical methods preserve certain nice properties.

We are not going to look at conservation of energy — this is too complicated

for a first course. Instead, we look at the following problem. Suppose we have a

system of ODEs for 0 ≤ t ≤ T :

(t) = f(t, y(t))

y(0) = y

Suppose T > 0 is arbitrary, and

lim

t→∞

y(t) = 0.

What restriction on

is necessary for a numerical method to satisfy

lim

n→∞

This question is still too complicated for us to tackle. It can only be easily

solved for linear problems, namely ODEs of the form

(t) = Ay(t),

for A ∈ R

N×N

Firstly, for what

do we have

(

)

→

0 as

t → ∞

? By some basic linear

algebra, we know this holds only if

(

)

0 for all eigenvalues

. To simplify

further, we consider the case where

(t) = λy(t), Re(λ) < 0.

It should be clear that if

is diagonalizable, then it can be reduced to multiple

instances of this case. Otherwise, we need to do some more work, but we’ll not

do that in this course.

And that, at last, is enough simplification.

6.2 Linear stability

This is the easy part of the course, where we are just considering the problem

λy

. No matter how complicated are numerical method is, when applied to

this problem, it usually becomes really simple.

Definition (Linear stability domain). If we apply a numerical method to

(t) = λy(t)

with y(0) = 1, λ ∈ C, then its linear stability domain is

D =

z = hλ : lim

n→∞

= 0

Example. Consider the Euler method. The discrete solution is

= (1 + hλ)

Thus we get

D = {z ∈ C : |1 + z| < 1}.

We can visualize this on the complex plane as the open unit ball

−1

Example.

Consider instead the backward Euler method. This method is

implicit, but since our problem is simple, we can find what

th step is. It turns

out it is

= (1 − λh)

−n

Then we get

D = {z ∈ C : |1 − z| > 1}.

We can visualize this as the region:

We make the following definition:

Definition (A-stability). A numerical method is A-stable if

−

= {z ∈ C : Re(z) < 0} ⊆ D.

In particular, for

(

)

< λ

, A-stability means that

will tend to 0 regardless

of how large h is.

Hence the backwards Euler method is A-stable, but the Euler method is not.

A-stability is a very strong requirement. It is hard to get A-stability. In

particular, Dahlquist proved that no multi-step method of order

p ≥

3 is A-stable.

Moreover, no explicit Runge-Kutta method can be A-stable.

So let’s look at some other implicit methods.

Example

(Trapezoidal rule)

Again consider

(

) =

λy

, with the trapezoidal

rule. Then we can find



1 + hλ/2

1 − hλ/2



So the linear stability domain is

D =



z ∈ C :



2 + z

2 − z



< 1



What this says is that

has to be closer to

−

2 than to 2. In other words,

exactly C

−

When testing a numerical method for A-stability in general, complex analysis

is helpful. Usually, when applying a numerical method to the problem y

= λy,

we get

= [r(hλ)]

where r is some rational function. So

D = {z ∈ C : |r(z)| < 1}.

We want to know if

contains the left half-plane. For more complicated

expressions of

, like the case of the trapezoidal rule, this is not so obvious.

Fortunately, we have the maximum principle:

Theorem

(Maximum principle)

Let

be analytic and non-constant in an open

set Ω ⊆ C. Then |g| has no maximum in Ω.

Since

|g|

needs to have a maximum in the closure of Ω, the maximum must

occur on the boundary. So to show

|g| ≤

1 on the region Ω, we only need to

show the inequality holds on the boundary ∂Ω.

We try Ω =

−

. The trick is to first check that

is analytic in Ω, and

then check what happens at the boundary. This technique is made clear in the

following example:

Example. Consider

r(z) =

6 − 2z

6 − 4z + z

This is still pretty simple, but can illustrate how we can use the maximum

principle.

We first check if it is analytic. This certainly has some poles, but they are

2 ±

√

2i, and are in the right-half plane. So this is analytic in C

−

Next, what happens at the boundary of the left-half plane? Firstly, as

|z| → ∞

, we find

(

)

→

0, since we have a

at the denominator. The next

part is checking when

is on the imaginary axis, say

with

t ∈ R

. Then

we can check by some messy algebra that

|r(it)| ≤ 1

for

t ∈ R

. Therefore, by the maximum principle, we must have

(

)

| ≤

1 for all

z ∈ C

−

7 Implementation of ODE methods

We just did quite a lot of theory about numerical methods. To end this section,

we will look at some more practical sides of ODE methods.

7.1 Local error estimation

The first problem we want to tackle is what

we should pick. Usually, when

using numerical analysis software, you will be asked for an error tolerance, and

then the software will automatically compute the required

we need. How is

this done?

Milne’s device is a method for estimating the local truncation error of a

multi-step method, and hence changing the step-length

(there are similar

techniques for Runge-Kutta methods, but they are more complicated). This uses

two multi-step methods of the same order.

To keep things simple, we consider the two-step Adams-Bashforth method.

Recall that this is given by

n+1

= y

(2f(t

, y

) − f(t

n−1

, y

n−1

)).

This is an order 2 error with

n+1

000

) + O(h

The other multi-step method of order 2 is our old friend, the trapezoidal rule.

This is an implicit method

n+1

= y

(f(t

, y

) + f(t

n+1

, y

n+1

)).

This has a local truncation error of

n+1

= −

000

) + O(h

The key to Milne’s device is the coefficients of h

000

), namely

, c

= −

Since these are two different methods, we get different

n+1

’s. We distinguish

these by superscripts, and have

y(t

n+1

) − y

n+1

' c

000

)

y(t

n+1

) − y

n+1

' c

000

)

We can now eliminate some terms to obtain

y(t

n+1

) − y

n+1

−c

− c

n+1

− y

n+1

In this case, the constant we have is

. So we can estimate the local truncation

error for the trapezoidal rule, without knowing the value of

000

. We can then

use this to adjust h accordingly.

The extra work we need for this is to compute numerical approximations

with two methods. Usually, we will use a simple, explicit method such as the

Adams-Bashforth method as the second method, when we want to approximate

the error of a more complicated but nicer method, like the trapezoidal rule.

7.2 Solving for implicit methods

Implicit methods are often more likely to preserve nice properties like conservation

of energy. Since we have more computational power nowadays, it is often

preferable to use these more complicated methods. When using these implicit

methods, we have to come up with some way to solve the equations involved.

As an example, we consider the backward Euler method

n+1

= y

+ hf(t

n+1

, y

n+1

There are two ways to solve this equation for

n+1

. The simplest method is

functional iteration. As the name suggests, this method is iterative. So we use

superscripts to denote the iterates. In this case, we use the formula

(k+1)

n+1

= y

+ hf(t

n+1

, y

n+1

To do this, we need to find a

(0)

n+1

. Usually, we start

(0)

n+1

. Even better,

we can use some simpler explicit method to obtain our first guess of y

(0)

n+1

The question, of course, is whether this converges. Fortunately, this converges

to a locally unique solution if

λh

is sufficiently small, where

is the Lipschitz

constant of

. For the backward Euler, we will require

λh <

1. This relies on

the contraction mapping theorem, which you may have met in IB Analysis II.

Does this matter? Sometimes it does. Usually, we pick an h using accuracy

considerations, picking the largest possible

that still gives us the desired

accuracy. However, if we use this method, we might need to pick a much smaller

in order for this to work. This will require us to compute much more iterations,

and can take a lot of time.

An alternative is Newton’s method. This is given by the formula

(I − hJ

(k)

= y

(k)

n+1

− (y

+ hf(t

n+1

, y

(k)

n+1

))

(k+1)

n+1

= y

(k)

− z

(k)

where J

(k)

is the Jacobian matrix

(k)

= ∇f(t

n+1

, y

(k)

n+1

) ∈ R

N×N

This requires us to solve for

in the first equation, but this is a linear system,

which we have some efficient methods for solving.

There are several variants to Newton’s method. This is the full Newton’s

method, where we re-compute the Jacobian in every iteration. It is also possible

to just use the same Jacobian over and over again. There are some speed gains

in solving the equation, but then we will need more iterations before we can get

our y

n+1

8 Numerical linear algebra

In the last part of the course, we study numerical algorithms for solving certain

problems in linear algebra. Here we are not so much concerned about accuracy

— the solutions we obtained are theoretically exact. Instead, we will try to find

some efficient methods of solving these problems. However, we will occasionally

make some comments about accuracy problems that arise from the fact that we

only work with finite precision.

We start off with the simplest problem in linear algebra: given

A ∈ R

n×n

b ∈ R

, we want to find an x such that

Ax = b.

We all know about the theory — if

is non-singular, then there is a unique

solution for every possible b. Otherwise, there are no solutions for some b and

infinitely many solutions for other b’s.

Most of the time, we will only consider the case where

is non-singular.

However, we will sometimes comment on what happens when A is singular.

8.1 Triangular matrices

While matrices in general and big and scary, triangular matrices tend to be much

nicer.

Definition

(Triangular matrix)

A matrix

is upper triangular if

= 0

whenever

i > j

. It is lower triangular if

= 0 whenever

i < j

. We usually

denote upper triangular matrices as U and lower triangular matrices as L.

We can visualize these matrices as follows:

L U

Why are triangular matrices nice? First of all, it is easy to find the determinants

of triangular matrices. We have

det(L) =

i=1

, det(U) =

i=1

In particular, a triangular matrix is non-singular if and only if it has no zero

entries in its diagonal.

How do we solve equations involving triangular matrices? Suppose we have

a lower triangular matrix equation







0 ··· 0

··· 0

··· L































or, more visually, we have

There is nothing to solve in the first line. We can immediately write down the

value of

. Substituting this into the second line, we can then solve for

. In

general, we have





−

i−1

j=1





This is known as forward substitution.

For upper triangular matrices, we can do a similar thing, but we have to

solve from the bottom instead. We have





−

j=1+i





This is known as backwards substitution.

This can be used to find the inverse of matrices as well. The solution to

is the

th column of

−1

. It is then easy to see from this that the

inverse of a lower triangular matrix is also lower triangular.

Similarly, the columns of

−1

are given by solving

, and we see that

−1

is also upper triangular.

It is helpful to analyse how many operations it takes to compute this. If we

look carefully, solving Lx = b or U x = b this way requires O(n

) operations.

8.2 LU factorization

In general, we don’t always have triangular matrices. The idea is, for every

matrix A, find some lower and triangular matrices L, U such that A = LU .

If we can do this, then we can solve

in two steps — we first find a

such that Ly = b. Then we find an x such that Ux = y. Then

Ax = LUx = Ly = b.

So what we want is such a factorization. To guarantee uniqueness of the

factorization, we require that

is unit, i.e. the diagonals are all 1. Otherwise,

given any such factorization, we can divide

by some (non-zero) constant and

multiply L by the same constant to get another factorization.

Definition

(LU factorization)

. A

is an LU factorization if

is upper

triangular and L is unit lower triangular (i.e. the diagonals of L are all 1).

Note that since

has to be unit, it must be non-singular. However, we still

allow A and U to be singular. Note that

det(A) = det(L) det(U) = det(U).

So A is singular if and only if U is.

Unfortunately, even if

is non-singular, it may not have an LU factorization.

We take a simple example of



0 1

1 1



This is clearly non-singular, with determinant

−

1. However, we can manually

check that there is no LU factorization of this.

On the other hand, while we don’t really like singular matrices, singular

matrices can still have LU factorizations. For example,



0 0

0 1





1 0

0 1



0 0

0 1



is trivially an LU factorization of a singular matrix.

Fortunately, if

is non-singular and has an LU factorization, then this

factorization is unique (this is not necessarily true if

is singular. For example,

we know



0 0

0 1





1 0

a 1



0 0

0 1



for any real number a).

To understand when LU factorizations exist, we try to construct one, and

see when it could fail.

We write L in terms of its columns, and U in terms of its rows:

L =



··· l



, U =













Clearly, these rows and columns cannot be arbitrary, since

and

are triangular.

In particular, l

, u

must be zero in the first i − 1 entries.

Suppose this is an LU factorization of A. Then we can write

A = L · U = l

+ l

+ ··· + l

What do these matrices look like? For each

, we know

and

have the first

i −

1 entries zero. So the first

i −

1 rows and columns of

are zero. In

particular, the first row and columns only have contributions from

, the

second row/column only has contributions from l

and l

etc.

The plan is as follows:

(i)

Obtain

and

from the first row and column of

. Since the first entry

is 1,

is exactly the first row of

. We can then obtain

by taking

the first column of A and dividing by U

= A

(ii) Obtain l

and u

form the second row and column of A − l

similarly.

(iii) ···

(iv) Obtain l

and u

from the nth row and column of A −

n−1

i=1

We can turn this into an algorithm. We define the intermediate matrices, starting

with

(0)

= A.

For k = 1, ··· , n, we let

= A

(k−1)

j = k, ··· , n

(k−1)

i = k, ··· , n

(k)

= A

(k−1)

− L

i, j ≥ k

When

, we end up with a zero matrix, and then

and

are completely

filled.

We can now see when this will break down. A sufficient condition for

to exist is that

(k−1)

= 0 for all

. Since

(k−1)

, this sufficient condition

ensures

, and hence

is non-singular. Conversely, if

is non-singular and

an LU factorization exists, then this would always work, since we must have

(k−1)

= 0. Moreover, the LU factorization must be given by this

algorithm. So we get uniqueness.

The problem with this sufficient condition is that most of these coefficients

do not appear in the matrix

. They are constructed during the algorithm. We

don’t know easily what they are in terms of the coefficients of

. We will later

come up with an equivalent condition on our original A that is easier to check.

Note that as long as this method does not break down, we need

(

)

operations to perform this factorization. Recall we only needed

(

) operations

to solve the equation after factorization. So the bulk of the work in solving

Ax = b is in doing the LU factorization.

As before, this allows us to find the inverse of

if it is non-singular. In

particular, solving

gives the

th column of

−1

. Note that we are

solving the system for the same

for each

. So we only have to perform the

LU factorization once, and then solve

different equations. So in total we need

O(n

) operations.

However, we still have the problem that factorization is not always possible.

Requiring that we must factorize

is too restrictive. The idea is to factor

something closely related, but not exactly A. Instead, we want a factorization

P A = LU,

where

is a permutation matrix. Recall a permutation matrix acting on a

column vector just permutes the elements of the vector, and

acting on

would just permute the rows of

. So we want to factor

up to a permutation

of rows.

Our goal now is to extend the previous algorithm to allow permutations of

rows, and then we shall show that we will be able to perform this factorization

all the time.

Suppose our breakdown occurs at

= 1, i.e.

(0)

= 0. We find a

permutation matrix

and let it act via

(0)

. The idea is to look down the

first column of

, and find a row starting with a non-zero element, say

. Then

we use

to interchange rows 1 and

such that

(0)

has a non-zero top-most

entry. For simplicity, we assume we always need a

, and if

(0)

is non-zero in

the first place, we just take P

to be the identity.

After that, we can carry on. We construct

and

from

(0)

as before,

and set A

(1)

= P

(0)

− l

But what happens if the first column of

is completely zero? Then no

interchange will make the (1

1) entry non-zero. However, in this case, we don’t

actually have to do anything. We can immediately find our

and

, namely

set

(or anything) and let

be the first row of

(0)

. Then this already

works. Note however that this corresponds to

(and hence

) being singular,

and we are not too interested with these.

The later steps are exactly analogous. Suppose we have

(k−1)

= 0. Again

we find a

such that

(k−1)

has a non-zero (

k, k

) entry. We then construct

and u

from P

(k−1)

and set

(k)

= P

(n−1)

− l

Again, if the

th column of

(k−1)

is completely zero, we set

and

be the kth row of A

(k−1)

. But again this implies A and U will be singular.

However, as we do this, the permutation matrices appear all over the place

inside the algorithm. It is not immediately clear that we do get a factorization

of the form

P A

. Fortunately, keeping track of the interchanges, we do

have an LU factorization

P A =

LU,

where U is what we got from the algorithm,

P = P

n−1

···P

while

L is given by

L =



···



= P

n−1

···P

k−1

Note that in particular, we have

n−1

= l

n−1

= l

One problem we have not considered is the problem of inexact arithmetic. While

these formula are correct mathematically, when we actually implement things,

we do them on computers with finite precision. As we go through the algorithm,

errors will accumulate, and the error might be amplified to a significant amount

when we get the reach the end. We want an algorithm that is insensitive to

errors. In order to work safely in inexact arithmetic, we will put the element of

largest modulus in the (

k, k

)th position, not just an arbitrary non-zero one, as

this minimizes the error when dividing.

8.3 A = LU for special A

There is one thing we haven’t figured out yet. Can we just look at the matrix

and then determine whether it has an LU factorization (that does not involve

permutations)? To answer this, we start with some definitions.

Definition

(Leading principal submatrix)

The leading principal submatrices

∈ R

k×k

for k = 1, ··· , n of A ∈ R

n×n

are

)

= A

, i, j = 1, ··· , k.

In other words, we take the first k rows and columns of A.

Theorem.

A sufficient condition for the existence for both the existence and

uniqueness of A = LU is that det(A

) 6= 0 for k = 1, ··· , n − 1.

Note that we don’t need

to be non-singular. This is equivalent to the

restriction

(k−1)

= 0 for

= 1

, ··· , n −

1. Also, this is a sufficient condition,

not a necessary one.

Proof. Straightforward induction.

We extend this result a bit:

Theorem.

det

(

)

= 0 for all

= 1

, ··· , n

, then

A ∈ R

n×n

has a unique

factorization of the form

A = LD

where

is non-singular diagonal matrix, and both

and

are unit triangular.

Proof.

From the previous theorem,

exists. Since

is non-singular,

is non-singular. So we can write this as

U = D

where D consists of the diagonals of U and

U = D

−1

U is unit.

This is not hard, but is rather convenient. In particular, it allows us to

factorize symmetric matrices in a nice way.

Theorem.

Let

A ∈ R

n×n

be non-singular and

det

(

)

= 0 for all

= 1

, ··· , n

Then there is a unique “symmetric” factorization

A = LDL

with L unit lower triangular and D diagonal and non-singular.

Proof. From the previous theorem, we can factorize A uniquely as

A = LD

We take the transpose to obtain

A = A

This is a factorization of the form “unit lower”-“diagonal”-“unit upper”. By

uniqueness, we must have

U = L

. So done.

We just looked at a special type of matrices, the symmetric matrices. We

look at another type of matrices:

Definition

(Positive definite matrix)

A matrix

A ∈ R

n×n

is positive-definite if

Ax > 0

for x 6= 0 ∈ R

Theorem.

Let

A ∈ R

n×n

be a positive-definite matrix. Then

det

(

)

= 0 for

all k = 1, ··· , n.

Proof.

First consider

. To show

is non-singular, it suffices to show that

implies

. But we can multiply the equation by

to obtain

Ax = 0. By positive-definiteness, we must have x = 0. So done.

Now suppose

for

k < n

and

. Then

= 0. We

invent a new

x ∈ R

by taking

and pad it with zeros. Then

= 0. By

positive-definiteness, we know x = 0. Then in particular y = 0.

We are going to use this to prove the most important result in this section.

We are going to consider matrices that are symmetric and positive-definite.

Theorem.

A symmetric matrix

A ∈ R

n×n

is positive-definite iff we can factor

it as

A = LDL

where L is unit lower triangular, D is diagonal and D

> 0.

Proof. First suppose such a factorization exists, then

Ax = x

LDL

x = (L

D(L

x).

We let

. Note that

if and only if

, since

is invertible. So

Ax = y

Dy =

k=1

> 0

if y 6= 0.

Now if

is positive definite, it has an LU factorization, and since

symmetric, we can write it as

A = LDL

where

is unit lower triangular and

is diagonal. Now we have to show

0. We define

such that

, which exist, since

is invertible.

Then clearly y

6= 0. Then we have

= e

= y

LDL

= y

> 0.

So done.

This is a practical check for symmetric

being positive definite. We can

perform this LU factorization, and then check whether the diagonal has positive

entries.

Definition

(Cholesky factorization)

The Cholesky factorization of a symmetric

positive-definite matrix A is a factorization of the form

A = LDL

with L unit lower triangular and D a positive-definite diagonal matrix.

There is another way of doing this. We let

1/2

be the “square root” of

by taking the positive square root of the diagonal entries of D. Then we have

A = LDL

= LD

1/2

= (LD

1/2

)(LD

1/2

)

= GG

where

is lower triangular with

0. This is another way of presenting this

result.

Finally, we look at the LU factorization of band matrices.

Definition

(Band matrix)

A band matrix of band width

is a matrix

such

that A

6= 0 implies |i − j| ≤ r.

For example, a band matrix of band width 0 is a diagonal matrix; a band

matrix of band width 1 is a tridiagonal matrix.

The result is

Proposition.

If a band matrix

has band width

and an LU factorization

A = LU, then L and U are both band matrices of width r.

Proof. Straightforward verification.

9 Linear least squares

Finally, we consider a slightly different kind of problem, whose solution involves

a different kind of factorization.

The question we are interested in is how we can “solve”

Ax = b,

where A ∈ R

m×n

, with m > n, and b ∈ R

and x ∈ R

Of course, this is an over-determined system. In general, this has no solution.

So we want to find the “best” solution to this system of equations. More precisely,

we want to find an

∗

∈ R

that minimizes

kAx − bk

. Note that we use the

Euclidean norm. This is the least squares problem.

Why would one be interested in this? Often in, say, statistics, we have some

linear model that says the outcome

is related to the input variables

, ··· , A

B = A

+ A

+ ··· + A

+ ε,

where the

are some unknown parameters, and

is some random fluctuation.

The idea is to do lots of experiments, say

of them, collect the data in

A ∈ R

m×n

and

b ∈ R

, and then find the

’s that predict these results the most accurately,

i.e. find the x such that Ax is as close to b as possible.

So how can we find this

∗

? First of all, we find a necessary and sufficient

condition for x

∗

to be a solution to this.

Theorem. A vector x

∗

∈ R

minimizes kAx − bk

if and only if

(Ax

∗

− b) = 0.

Proof. A solution, by definition, minimizes

f(x) = hAx − b, Ax − bi

= x

AAx − 2x

b + b

Then as a function of x, the partial derivatives of this must vanish. We have

∇f(x) = 2A

(Ax − b).

So a necessary condition is

(Ax − b).

Now suppose our

∗

satisfies

(

∗

− b

) = 0. Then for all

x ∈ R

, we write

x = x

∗

+ y, and then we have

kAx − bk

= kA(x

∗

+ y) − bk

= kAx

∗

− bk + 2y

(Ax − b) + kAyk

= kAx

∗

− bk + kAyk

≥ kAx

∗

− bk.

So x

∗

must minimize the Euclidean norm.

Corollary.

A ∈ R

m×n

is a full-rank matrix, then there is a unique solution

to the least squares problem.

Proof. We know all minimizers are solutions to

A)x = A

The matrix

being full rank means

y 6

0 ∈ R

implies

Ay 6

0 ∈ R

. Hence

A ∈ R

n×n

is positive definite (and in particular non-singular), since

Ax = (Ax)

(Ax) = kAxk

> 0

for x 6= 0. So we can invert A

A and find a unique solution x.

Now to find the minimizing x

∗

, we need to solve the normal equations

Ax = A

has full rank, then

is non-singular, and there is a unique solution. If

not, then the general theory of linear equations tells us there are either infinitely

many solution or no solutions. But for this particular form of equations, it turns

out we can show this is always consistent. So there will always be a solution.

However, solving these equations is not the recommended algorithm for finding

∗

for accuracy reasons. Instead the solution of these is based on orthogonal

matrices.

Definition

(Orthogonal matrix)

A square matrix

Q ∈ R

n×n

is orthogonal if

Q = I.

A key property of the orthogonal matrices is that we have

kQxk = kxk

for all

x ∈ R

. This is helpful since what we are trying to do in the least-squares

problem involves minimizing the Euclidean norm of something.

The idea is that for any orthogonal matrix, minimizing

kAx −bk

is the same

as minimizing

kQAx − Qbk

. So

should we pick? Recall that a simple kind of

matrices is triangular matrices. So we might want to get something, say, upper

triangular.

Definition

(QR factorization)

A QR factorization of an

m × n

matrix

is a

factorization of the form

A = QR,

where

Q ∈ R

m×m

is an orthogonal matrix, and

R ∈ R

m×n

is “upper triangular”

matrix, where “upper triangular” means

R =







··· R

0 R

··· R

0 0 ··· R

0 0 ··· 0







Note that every

A ∈ R

m×n

has a QR factorization, as we will soon show,

but this is not unique (e.g. we can multiply both Q and R by −1).

Once we have the QR factorization, we can multiply

kAx −bk

−1

and get an equivalent problem of minimizing

kRx − Q

. We will not go into

details, but it should be clear that this is not too hard to solve.

There is an alternative “skinny” QR factorization, with

, where we

remove redundant rows and columns such that

Q ∈ R

m×n

and

R ∈ R

n×n

. This

is good enough to solve the least squares problem.

We will see that if

has full rank, then the skinny QR is unique up to a

sign, i.e. if we require

> 0 for k = 1, ··· , n.

We shall demonstrate three standard algorithms for QR factorization:

(i) Gram-Schmidt factorization

(ii) Givens rotations

(iii) Householder reflections

The latter two are based on some geometric ideas, while the first is based on the

familiar Gram-Schmidt process.

Gram-Schmidt factorization

This targets the skinny version. So we stop writing the tildes above, and just

write Q ∈ R

m×n

and R ∈ R

n×n

As we all know, the Gram-Schmidt process orthogonalizes vectors. What we

are going to orthogonalize are the columns of A. We write

A =



··· a



, Q =



··· q



By definition of the QR factorization, we need

i=1

This is just done in the usual Gram-Schmidt way.

(i) To construct column 1, if a

6= 0, then we set

, R

= ka

Note that the only non-unique possibility here is the sign — we can let

−ka

and

−

instead. But if we require

0, then

this is fixed.

In the degenerate case,

, then we can just set

= 0, and the pick

any q

∈ R

with kq

k = 1.

(ii) For columns 1 < k ≤ n, for i = 1, ··· , k − 1, we set

= hq

, a

and set

= a

−

k−1

i=1

, a

If d

6= 0, then we set

and

= kd

In the case where

, we again set

= 0, and pick

to be anything

orthonormal to q

, ··· , q

k−1

In practice, a slightly different algorithm (modified Gram-Schmidt process)

is used, which is (much) superior with inexact arithmetic. The modified Gram-

Schmidt process is in fact the same algorithm, but performed in a different order

in order to minimize errors.

However, this is often not an ideal algorithm for large matrices, since there

are many divisions and normalizations involved in computing the

, and the

accumulation of errors will cause the resulting matrix Q to lose orthogonality.

Givens rotations

This works with the full QR factorization.

Recall that in R

, a clockwise rotation of θ ∈ [−π, π] is performed by



cos θ sin θ

−sin θ cos θ







α cos θ + β sin θ

−α sin θ + β cos θ



By choosing θ such that

cos θ =

+ β

, sin θ =

+ β

this then becomes



+ β



Of course, by choosing a slightly different

, we can make the result zero in the

first component and

+ β

in the second.

Definition

(Givens rotation)

, where

m >

2, we define the Givens

rotation on 3 parameters 1 ≤ p < q ≤ m, θ ∈ [−π, π] by

Ω

[p,q]







cos θ sin θ

−sin θ cos θ







∈ R

m×m

where the sin and cos appear at the p, qth rows and columns.

Note that for

y ∈ R

, the effect of Ω

[p,q]

only alters the

and

components.

In general, for B ∈ R

m×n

, then Ω

[p,q]

B only alters the p and q rows of B.

Moreover, just like the

case, given a particular

z ∈ R

, we can choose

such that the qth component (Ω

[p,q]

= 0.

Hence,

A ∈ R

m×n

can be transformed into an “upper triangular” form by

applying

mn −

n(n+1)

Givens rotations, since we need to introduce

many

zeros. Then

···Q

A = R.

How exactly do we do this? Instead of writing down a general algorithm, we

illustrate it with a matrix A ∈ R

4×3

. The resultant R would look like

R =







0 R

0 0 R

0 0 0







We will apply the Givens rotations in the following order:

Ω

[3,4]

Ω

[3,4]

Ω

[2,3]

Ω

[1,4]

Ω

[1,3]

Ω

[1,2]

A = R.

The matrix A transforms as follows:

× × ×

Ω

[1,2]

× × ×

× ×

× × ×

Ω

[1,3]

× × ×

× ×

× × ×

Ω

[1,3]

× × ×

× ×

Ω

[2,3]

× × ×

× ×

0 0

× ×

Ω

[3,4]

× × ×

× ×

0 0

Ω

[3,4]

× × ×

× ×

0 0

0 0 0

Note that when applying, say, Ω

[2,3]

, the zeros of the first column get preserved,

since Ω

[2,3]

only mixes together things in row 2 and 3, both of which are zero in

the first column. So we are safe.

Note that this gives us something of the form

···Q

A = R.

We can obtain a QR factorization by inverting to get

A = Q

···Q

However, we don’t really need to do this if we just want to do use this to solve

the least squares problem, since to do so, we need to multiply by

, not

which is exactly Q

···Q

Householder reflections

Definition

(Householder reflection)

For

u 6

0 ∈ R

, we define the House-

holder reflection by

= I − 2

∈ R

m×m

Note that this is symmetric, and we can see

. So this is indeed

orthogonal.

To show this is a reflection, suppose we resolve

as the perpendicular and

parallel parts as x = αu + w ∈ R

, where

α =

, u

w = 0.

Then we have

x = −αu + w.

So this is a reflection in the m − 1-dimensional hyperplane u

y = 0.

What is the cost of computing H

z? This is evaluated as

z = z − 2

This only requires O(m) operations, which is nice.

Proposition.

A matrix

A ∈ R

m×n

can be transformed into upper-triangular

form by applying n Householder reflections, namely

···H

A = R,

where each

introduces zero into column

and leaves the other zeroes alone.

This is much less multiplications needed than the number needed for Givens

rotations, in general.

To show this, we first prove some lemmas.

Lemma.

Let

a, b ∈ R

, with

a 6

, but

kak

kbk

. Then if we pick

a −b

then

a = b.

This is obvious if we draw some pictures in low dimensions.

Proof. We just do it:

a = a −

2(kak − a

kak

− 2a

b + kbk

(a − b) = a − (a − b) = b,

where we used the fact that kak = kbk.

We will keep applying the lemma. Since we want to get an upper triangular

matrix, it makes sense to pick

to have many zeroes, and the best option would

be to pick b to be a unit vector.

So we begin our algorithm: we want to clear the first column of

. We let

be the first column of

, and assume

a ∈ R

is not already in the correct form,

i.e. a is not a multiple of e

. Then we define

u = a ∓ kake

where either choice of the sign is pure-mathematically valid. However, we will

later see that there is one choice that is better when we have to deal with inexact

arithmetic. Then we end up with

a = H

a = ±kake

Now we have

A =

× × ×

× ×

To do the next step, we need to be careful, since we don’t want to destroy the

previously created zeroes.

Lemma. If the first k − 1 components of u are zero, then

(i) For every x ∈ R

, H

x does not alter the first k − 1 components of x.

(ii) If the last (m − k + 1) components of y ∈ R

are zero, then H

y = y.

These are all obvious from definition. All these say is that reflections

don’t affect components perpendicular to

, and in particular fixes all vectors

perpendicular to u.

Lemma. Let a, b ∈ R

, with

























but

j=k

Suppose we pick

u = (0, 0, ··· , 0, a

− b

, ··· , a

− b

)

Then we have

a = (a

, ··· , a

k−1

, ··· , b

This is a generalization of what we’ve had before for

= 1. Again, the proof

is just straightforward verification.

Now we can mess with the second column of

. We let

be the second

column of

, and assume

, ··· , a

are not all zero, i.e. (0

, a

, ··· , a

)

not a multiple of e

. We choose

u = (0, a

∓ γ, a

, ··· , a

)

where

γ =

j=2

Then

a = (a

, ±γ, 0, ···).

Also, by the previous lemma, this does not affect anything in the first column

and the first row. Now we have

A =

× × ×

0 0

Suppose we have reached

k−1

···H

, where the first

k −

1 rows are of the

correct form, i.e.

k−1

···H

A =

k − 1

We consider the

th column of

k−1

···H

, say

. We assume

(0, ··· , 0, a

, ··· , a

)

is not a multiple of e

. Choosing

u = γ(0, ··· , 0, a

∓ γ, a

k+1

, ··· , a

)

, γ =

j=k

we find that

a = (a

, ··· , a

k−1

, ±γ, 0, ··· , 0)

Now we have

···H

A =

and H

did not alter the first k − 1 rows and columns of H

k−1

···H

There is one thing we have to decide on — which sign to pick. As mentioned,

these do not matter in pure mathematics, but with inexact arithmetic, we

should pick the sign in

∓ γ

such that

∓ γ

has maximum magnitude, i.e.

∓ γ

sgn

(

)

. It takes some analysis to justify why this is the right

choice, but it is not too surprising that there is some choice that is better.

So how do Householder and Givens compare? We notice that the Givens

method generates zeros at one entry at a time, while Householder does it column

by column. So in general, the Householder method is superior.

However, for certain matrices with special structures, we might need the

extra delicacy of introducing a zero one at a time. For example, if

is a band

matrix, then it might be beneficial to just remove the few non-zero entries one

at a time.