IB Methods (Full)

Part IB — Methods

Based on lectures by D. B. Skinner

Notes taken by Dexter Chua

Michaelmas 2015

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.

Self-adjoint ODEs

Periodic functions. Fourier series: definition and simple properties; Parseval’s theorem.

Equations of second order. Self-adjoint differential operators. The Sturm-Liouville

equation; eigenfunctions and eigenvalues; reality of eigenvalues and orthogonality of

eigenfunctions; eigenfunction expansions (Fourier series as prototype), approximation

in mean square, statement of completeness. [5]

PDEs on bounded domains: separation of variables

Physical basis of Laplace’s equation, the wave equation and the diffusion equation.

General method of separation of variables in Cartesian, cylindrical and spherical

co ordinates. Legendre’s equation: derivation, solutions including explicit forms of

and

, orthogonality. Bessel’s equation of integer order as an example of a

self-adjoint eigenvalue problem with non-trivial weight.

Examples including potentials on rectangular and circular domains and on a spherical

domain (axisymmetric case only), waves on a finite string and heat flow down a

semi-infinite rod. [5]

Inhomogeneous ODEs: Green’s functions

Properties of the Dirac delta function. Initial value problems and forced problems with

two fixed end points; solution using Green’s functions. Eigenfunction expansions of the

delta f unction and Green’s functions. [4]

Fourier transforms

Fourier transforms: definition and simple properties; inversion and convolution theorems.

The discrete Fourier transform. Examples of application to linear systems. Relationship

of transfer function to Green’s function for initial value problems. [4]

PDEs on unbounded domains

Classification of PDEs in two independent variables. Well posedness. Solution by

the method of characteristics. Green’s functions for PDEs in 1, 2 and 3 independent

variables; fundamental solutions of the wave equation, Laplace’s equation and the

diffusion equation. The method of images. Application to the forced wave equation,

Poisson’s equation and forced diffusion equation. Transient solutions of diffusion

problems: the error function. [6]

Contents

0 Introduction

1 Vector spaces

2 Fourier series

2.1 Fourier series

2.2 Convergence of Fourier series

2.3 Differentiability and Fourier series

3 Sturm-Liouville Theory

3.1 Sturm-Liouville operators

3.2 Least squares approximation

4 Partial differential equations

4.1 Laplace’s equation

4.2 Laplace’s equation in the unit disk in R

4.3 Separation of variables

4.4 Laplace’s equation in spherical polar coordinates

4.4.1 Laplace’s equation in spherical polar coordinates

4.4.2 Legendre Polynomials

4.4.3 Solution to radial part

4.5 Multipole expansions for Laplace’s equation

4.6 Laplace’s equation in cylindrical coordinates

4.7 The heat equation

4.8 The wave equation

5 Distributions

5.1 Distributions

5.2 Green’s functions

5.3 Green’s functions for initial value problems

6 Fourier transforms

6.1 The Fourier transform

6.2 The Fourier inversion theorem

6.3 Parseval’s theorem for Fourier transforms

6.4 A word of warning

6.5 Fourier transformation of distributions

6.6 Linear systems and response functions

6.7 General form of transfer function

6.8 The discrete Fourier transform

6.9 The fast Fourier transform*

7 More partial differential equations

7.1 Well-posedness

7.2 Method of characteristics

7.3 Characteristics for 2nd order partial differential equations

7.4 Green’s functions for PDEs on R

7.5 Poisson’s equation

0 Introduction

In the previous courses, the (partial) differential equations we have seen are

mostly linear. For example, we have Laplace’s equation:

∂

∂x

∂φ

∂y

= 0,

and the heat equation:

∂φ

∂t

= κ



∂

∂x

∂

∂y



The Schr¨odinger’ equation in quantum mechanics is also linear:

∂Φ

∂t

= −

∂

∂x

+ V (x)Φ(x).

By being linear, these equations have the property that if

, φ

are solutions,

then so are λ

+ λ

(for any constants λ

Why are all these linear? In general, if we just randomly write down a

differential equation, most likely it is not going to be linear. So where did the

linearity of equations of physics come from?

The answer is that the real world is not linear in general. However, often we

are not looking for a completely accurate and precise description of the universe.

When we have low energy/speed/whatever, we can often quite accurately approx-

imate reality by a linear equation. For example, the equation of general relativity

is very complicated and nowhere near being linear, but for small masses and

velocities, they reduce to Newton’s law of gravitation, which is linear.

The only exception to this seems to be Schr¨odinger’s equation. While there

are many theories and equations that superseded the Schr¨odinger equation, these

are all still linear in nature. It seems that linearity is the thing that underpins

quantum mechanics.

Due to the prevalence of linear equations, it is rather important that we

understand these equations well, and this is our primary objective of the course.

1 Vector spaces

When dealing with functions and differential equations, we will often think of

the space of functions as a vector space. In many cases, we will try to find a

“basis” for our space of functions, and expand our functions in terms of the basis.

Under different situations, we would want to use a different basis for our space.

For example, when dealing with periodic functions, we will want to pick basis

elements that are themselves periodic. In other situations, these basis elements

would not be that helpful.

A familiar example would be the Taylor series, where we try to approximate

a function f by

f(x) =

∞

n=0

(n)

(0)

Here we are thinking of

n ∈ N}

as the basis of our space, and trying to

approximate an arbitrary function as a sum of the basis elements. When writing

the function

as a sum like this, it is of course important to consider whether

the sum converges, and when it does, whether it actually converges back to f.

Another issue of concern is if we have a general set of basis functions

}

how can we find the coefficients

such that

(

) =

(

)? This is the

bit where linear algebra comes in. Finding these coefficients is something we

understand well in linear algebra, and we will attempt to borrow the results and

apply them to our space of functions.

Another concept we would want to borrow is eigenvalues and eigenfunctions,

as well as self-adjoint (“Hermitian”) operators. As we go along the course, we

will see some close connections between functions and vector spaces, and we can

often get inspirations from linear algebra.

Of course, there is no guarantee that the results from linear algebra would

apply directly, since most of our linear algebra results was about finite basis and

finite linear sums. However, it is often a good starting point, and usually works

when dealing with sufficiently nice functions.

We start with some preliminary definitions, which should be familiar from

IA Vectors and Matrices and/or IB Linear Algebra.

Definition

(Vector space)

A vector space over

(or

) is a set

with an

operation + which obeys

(i) u + v = v + u (commutativity)

(ii) (u + v) + w = u + (v + w) (associativity)

(iii) There is some 0 ∈ V such that 0 + u = u for all u (identity)

We can also multiply vectors by a scalars λ ∈ C, which satisfies

(i) λ(µv) = (λµ)v (associativity)

(ii) λ(u + v) = λu + λv (distributivity in V )

(iii) (λ + µ)u = λu + λv (distributivity in C)

(iv) 1v = v (identity)

Often, we wouldn’t have just a vector space. We usually give them some

additional structure, such as an inner product.

Definition

(Inner product)

An inner product on

is a map (

·, ·

) :

V ×V → C

that satisfies

(i) (u, λv) = λ(u, v) (linearity in second argument)

(ii) (u, v + w) = (u, v) + (u, w) (additivity)

(iii) (u, v) = (v, u)

∗

(conjugate symmetry)

(iv) (u, u) ≥ 0, with equality iff u = 0 (positivity)

Note that the positivity condition makes sense since conjugate symmetry entails

that (u, u) ∈ R.

The inner product in turn defines a norm

kuk

(u, u)

that provides the

notion of length and distance.

It is important to note that we only have linearity in the second argument.

For the first argument, we have (λu, v) = (v, λu)

∗

= λ

∗

(v, u)

∗

= λ

∗

(u, v).

Definition

(Basis)

A set of vectors

, v

, ··· , v

}

form a basis of

iff any

u ∈ V can be uniquely written as a linear combination

u =

i=1

for some scalars

. The dimension of a vector space is the number of basis

vectors in its basis.

A basis is orthogonal (with respect to the inner product) if (

, v

) = 0

whenever i 6= j.

A basis is orthonormal (with respect to the inner product) if it is orthogonal

and (v

, v

) = 1 for all i.

Orthonormal bases are the nice bases, and these are what we want to work

with.

Given an orthonormal basis, we can use the inner product to find the

expansion of any u ∈ V in terms of the basis, for if

u =

taking the inner product with v

gives

, u) =

, v

) = λ

using additivity and linearity. Hence we get the general formula

= (v

, u).

We have seen all these so far in IA Vectors and Matrices, where a vector is a list

of finitely many numbers. However, functions can also be thought of as elements

of an (infinite dimensional) vector space.

Suppose we have

f, g

: Ω

→ C

. Then we can define the sum

(f + g)(x) = f(x) + g(x). Given scalar λ, we can also define (λf )(x) = λf (x).

This also makes intuitive sense. We can simply view a functions as a list

of numbers, where we list out the values of

at each point. The list could be

infinite, but a list nonetheless.

Most of the time, we don’t want to look at the set of all functions. That

would be too huge and uninteresting. A natural class of functions to consider

would be the set of solutions to some particular differential solution. However,

this doesn’t always work. For this class to actually be a vector space, the sum

of two solutions and the scalar multiple of a solution must also be a solution.

This is exactly the requirement that the differential equation is linear. Hence,

the set of solutions to a linear differential equation would form a vector space.

Linearity pops up again.

Now what about the inner product? A natural definition is

(f, g) =

f(x)

∗

g(x) dµ,

where

is some measure. For example, we could integrate d

, or d

. This

measure specifies how much weighting we give to each point x.

Why does this definition make sense? Recall that the usual inner product

on finite-dimensional vector spaces is

∗

. Here we are just summing the

different components of

and

. We just said we can think of the function

a list of all its values, and this integral is just the sum of all components of

and g.

Example. Let Σ = [a, b]. Then we could take

(f, g) =

f(x)

∗

g(x) dx.

Alternatively, let Σ = D

⊆ R

be the unit disk. Then we could have

(f, g) =

2π

f(r, θ)

∗

g(r, θ) dθ r dr

Note that we were careful and said that d

is “some measure”. Here we are

integrating against d

θ r

. We will later see cases where this can be even more

complicated.

If Σ has a boundary, we will often want to restrict our functions to take

particular values on the boundary, known as boundary conditions. Often, we

want the boundary conditions to preserve linearity. We call these nice boundary

conditions homogeneous conditions.

Definition

(Homogeneous boundary conditions)

A boundary condition is

homogeneous if whenever

and

satisfy the boundary conditions, then so does

λf + µg for any λ, µ ∈ C (or R).

Example.

Let Σ = [

a, b

]. We could require that

(

) + 7

(

) = 0, or maybe

(

) + 3

(

) = 0. These are examples of homogeneous boundary conditions.

On the other hand, the requirement f(a) = 1 is not homogeneous.

2 Fourier series

The first type of functions we will consider is periodic functions.

Definition

(Periodic function)

A function

is periodic if there is some fixed

R such that f(x + R) = f (x) for all x.

However, it is often much more convenient to think of this as a function

→ C

from unit circle to

, and parametrize our function by an angle

ranging from 0 to 2π.

Now why do we care about periodic functions? Apart from the fact that

genuine periodic functions exist, we can also use them to model functions on a

compact domain. For example, if we have a function defined on [0

1], we can

pretend it is a periodic function on

by making infinitely many copies of the

function to the intervals [1, 2], [2, 3] etc.

x x

2.1 Fourier series

So. As mentioned in the previous chapter, we want to find a set of “basis

functions” for periodic functions. We could go with the simplest case of periodic

functions we know of — the exponential function

inθ

. These functions have a

period of 2

, and are rather easy to work with. We all know how to integrate

and differentiate the exponential function.

More importantly, this set of basis functions is orthogonal. We have

imθ

, e

inθ

) =

−π

−imθ

inθ

dθ =

−π

i(n−m)θ

dθ =

(

2π n = m

0 n 6= m

= 2πδ

We can normalize these to get a set of orthonormal functions

√

2π

inθ

: n ∈ Z

Fourier’s idea was to use this as a basis for any periodic function. Fourier

claimed that any f : S

→ C can be expanded in this basis:

f(θ) =

n∈Z

inθ

where

2π

inθ

, f) =

2π

−π

−inθ

f(θ) dθ.

These really should be defined as

(

) =

inθ

√

2π

with



inθ

√

2π

, f



, but for

convenience reasons, we move all the constant factors to the

coefficients.

We can consider the special case where

→ R

is a real function. We

might want to make our expansion look a bit more “real”. We get

(

)

∗



2π

−π

−inθ

f(θ) dθ



∗

2π

−π

inθ

f(θ) dθ =

−n

So we can replace our Fourier series by

f(θ) =

∞

n=1



inθ

−n

−inθ



∞

n=1



inθ

∗

−inθ



Setting

+ ib

, we can write this as

f(θ) =

∞

n=1

cos nθ + b

sin nθ)

∞

n=1

cos nθ + b

sin nθ).

Here the coefficients are

−π

cos nθf(θ) dθ, b

−π

sin nθf(θ) dθ.

This is an alternative formulation of the Fourier series in terms of sin and cos.

So when given a real function, which expansion should we use? It depends. If

our function is odd (or even), it would be useful to pick the sine/cosine expansion,

since the cosine (or sine) terms will simply disappear. On the other hand, if we

want to stick our function into a differential equation, exponential functions are

usually more helpful.

2.2 Convergence of Fourier series

When Fourier first proposed the idea of a Fourier series, people didn’t really

believe in him. How can we be sure that the infinite series actually converges?

It turns out that in many cases, they don’t.

To investigate the convergence of the series, we define the partial Fourier

sum as

f =

m=−n

imθ

The question we want to answer is whether

“converges” to

. Here we

have to be careful with what we mean by convergence. As we (might) have seen

in Analysis, there are many ways of defining convergence of functions. If we

have a “norm” on the space of functions, we can define convergence to mean

lim

n→∞

f − fk = 0. Our “norm” can be defined as

f − fk =

2π

−π

f(θ) − f(θ)|

dθ.

However, this doesn’t really work if

and

can be arbitrary functions, as in

this does not necessarily have to be a norm. Indeed, if

lim

n→∞

differs from

on finitely or countably many points, the integral will still be zero. In particular,

lim S

and

can differ at all rational points, but this definition will say that

f converges to f.

Hence another possible definition of convergence is to require

lim

n→∞

f(θ) − f(θ) = 0

for all

. This is known as pointwise convergence. However, this is often too

weak a notion. We can ask for more, and require that the rate of convergent is

independent of θ. This is known as uniform convergence, and is defined by

lim

n→∞

sup

f(θ) − f(θ)| = 0.

Of course, with different definitions of convergence, we can get different answers

to whether it converges. Unfortunately, even if we manage to get our definition

of convergence right, it is still difficult to come up with a criterion to decide if a

function has a convergent Fourier Series.

It is not difficult to prove a general criterion for convergence in norm, but

we will not do that, since that is analysis. If interested, one should take the IID

Analysis of Functions course. In this course, instead of trying to come up with

something general, let’s look at an example instead.

Example. Consider the sawtooth function f(θ) = θ.

−3π

−π

3π

Note that this function is discontinuous at odd multiples of π.

The Fourier coefficients (for n 6= 0) are

2π

−π

−inθ

θ dθ =



−

2πin

−inθ



−π

2πin

−π

inθ

dθ =

(−1)

n+1

We also have

= 0.

Hence we have

θ =

n6=0

(−1)

n+1

inθ

It turns out that this series converges to the sawtooth for all

θ 6

= (2

+ 1)

, i.e.

everywhere that the sawtooth is continuous.

Let’s look explicitly at the case where

. Each term of the partial Fourier

series is zero. So we can say that the Fourier series converges to 0. This is the

average value of lim

ε→0

f(π ± ε).

This is typical. At an isolated discontinuity, the Fourier series is the average

of the limiting values of the original function as we approach from either side.

2.3 Differentiability and Fourier series

Integration is a smoothening operator. If we take, say, the step function

Θ(x) =

(

1 x > 0

0 x < 0

Then the integral is given by

Θ(x) dx =

(

x x > 0

0 x < 0

This is now a continuous function. If we integrate it again, we make the positive

side look like a quadratic curve, and it is now differentiable.

On the other hand, when we differentiate a function, we generally make it worse.

For example, when we differentiate the continuous function

Θ(

) d

, we obtain

the discontinuous function Θ(

). If we attempt to differentiate this again, we

get the δ-(non)-function.

Hence, it is often helpful to characterize the “smoothness” of a function by

how many times we can differentiate it. It turns out this is rather relevant to

the behaviour of the Fourier series.

Suppose that we have a function that is itself continuous and whose first

m −

1 derivatives are also continuous, but

(m)

has isolated discontinuities at

{θ

, θ

, ··· , θ

We can look at the kth Fourier coefficient:

2π

−π

−ikθ

f(θ) dθ



−

2πik

−ikθ

f(θ)



−π

2πik

−π

−ikθ

(θ) dθ

The first term vanishes since

(

) is continuous and takes the same value at

and −π. So we are left with

2πik

−π

−ikθ

(θ) dθ

= ···

(ik)

2π

−π

ikθ

(m)

(θ) dθ.

Now we have to be careful, since

(m)

is no longer continuous. However provided

that

(m)

is everywhere finite, we can approximate this by removing small strips

(

−ε, θ

) from our domain of integration, and take the limit

ε →

0. We can

write this as

= lim

ε→0

(ik)

2π

−ε

−π

−ε

+ε

+ ··· +

+ε

ikθ

(m)

(θ) dθ

(ik)

m+1

2π

s=1

(m)

(θ

) − f

(m)

(θ

−

)) +

(−π,π)\θ

−ikθ

(m+1)

(θ) dθ

We now have to stop. So

decays like





m+1

if our function and its (

m −

derivatives are continuous. This means that if a function is more differentiable,

then the coefficients decay more quickly.

This makes sense. If we have a rather smooth function, then we would expect

the first few Fourier terms (with low frequency) to account for most of the

variation of f. Hence the coefficients decay really quickly.

However, if the function is jiggly and bumps around all the time, we would

expect to need some higher frequency terms to account for the minute variation.

Hence the terms would not decay away that quickly. So in general, if we can

differentiate it more times, then the terms should decay quicker.

An important result, which is in some sense a Linear Algebra result, is

Parseval’s theorem.

Theorem (Parseval’s theorem).

(f, f) =

−π

|f(θ)|

dθ = 2π

n∈Z

Proof.

(f, f) =

−π

|f(θ)|

dθ

−π

m∈Z

∗

−imθ

n∈Z

inθ

dθ

m,n∈Z

∗

−π

i(n−m)θ

dθ

= 2π

m,n∈Z

∗

= 2π

n∈Z

Note that this is not a fully rigorous proof, since we assumed not only that

the Fourier series converges to the function, but also that we could commute the

infinite sums. However, for the purposes of an applied course, this is sufficient.

Last time, we computed that the sawtooth function

(

) =

has Fourier

coefficients

= 0,

i(−1)

for n 6= 0.

But why do we care? It turns out this has some applications in number theory.

You might have heard of the Riemann ζ-function, defined by

ζ(s) =

∞

n=1

We will show that this obeys the property that for any

) =

for

some

q ∈ Q

. This may not be obvious at first sight. So let’s apply Parseval’s

theorem for the sawtooth defined by

(

) =

. By direct computation, we know

that

(f, f) =

−π

dθ =

2π

However, by Parseval’s theorem, we know that

(f, f) = 2π

n∈Z

= 4π

∞

n=1

Putting these together, we learn that

∞

n=1

= ζ(2) =

We have just done it for the case where

= 1. But if we integrate the sawtooth

function repeatedly, then we can get the general result for all m.

At this point, we might ask, why are we choosing these

imθ

as our basis?

Surely there are a lot of different sets of basis we can use. For example, in finite

dimensions, we can just arbitrary choose random vectors (that are not linearly

dependent) to get a set of basis vectors. However, in practice, we don’t pick

them randomly. We often choose a basis that can reveal the symmetry of a

system. For example, if we have a rotationally symmetric system, we would like

to use polar coordinates. Similarly, if we have periodic functions, then

imθ

often a good choice of basis.

3 Sturm-Liouville Theory

3.1 Sturm-Liouville operators

In finite dimensions, we often consider linear maps

V → W

. If

}

is a

basis for

and

}

is a basis for

, then we can represent the map by a

matrix with entries

= (w

, Mv

A map

V → V

is called self-adjoint if

†

as matrices. However, it is

not obvious how we can extend this notion to arbitrary maps between arbitrary

vector spaces (with an inner product) when they cannot be represented by a

matrix.

Instead, we make the following definitions:

Definition

(Adjoint and self-adjoint)

The adjoint

of a map

V → V

is a

map such that

(Bu, v) = (u, Av)

for all vectors u, v ∈ V . A map is then self-adjoint if

(Mu, v) = (u, M v).

Self-adjoint matrices come with a natural basis. Recall that the eigenvalues

of a matrix are the roots of

det

(

M − λI

) = 0. The eigenvector corresponding to

an eigenvalue λ is defined by Mv

= λ

In general, eigenvalues can be any complex number. However, self-adjoint

maps have real eigenvalues. Suppose

= λ

Then we have

, v

) = (v

, Mv

) = (Mv

, v

) = λ

∗

, v

So λ

= λ

∗

Furthermore, eigenvectors with distinct eigenvalues are orthogonal with

respect to the inner product. Suppose that

= λ

, M v

= λ

Then

, v

) = (v

, Mv

) = (Mv

, v

) = λ

, v

Since λ

6= λ

, we must have (v

, v

) = 0.

Knowing eigenvalues and eigenvalues gives a neat way so solve linear equations

of the form

Mu = f.

Here we are given

and

, and want to find

. Of course, the answer is

u = M

−1

f. However, if we expand in terms of eigenvectors, we obtain

Mu = M

Hence we have

Taking the inner product with v

, we know that

So far, these are all things from IA Vectors and Matrices. Sturm-Liouville theory

is the infinite-dimensional analogue.

In our vector space of differentiable functions, our “matrices” would be linear

differential operators L. For example, we could have

L = A

(x)

+ A

p−1

(x)

p−1

+ ··· + A

(x)

+ A

(x).

It is an easy check that this is in fact linear.

We say L has order p if the highest derivative that appears is

In most applications, we will be interested in the case

= 2. When will our

L be self-adjoint?

In the p = 2 case, we have

Ly = P

+ R

− Qy

= P



−



= P



−





−



Let p = exp





. Then we can write this as

= P p

−1







−



We further define

. We also drop a factor of

P p

−1

. Then we are left with

L =



p(x)



− q(x).

This is the Sturm-Liouville form of the operator. Now let’s compute (

f, Lg

). We

integrate by parts numerous times to obtain

(f, Lg) =

∗





− qg



= [f

∗

]

−



∗

+ f

∗



= [f

∗

− f

0∗

pg]



∗



− qf

∗



g dx

= [(f

∗

− f

0∗

g)p]

+ (Lf, g),

assuming that p, q are real.

So 2nd order linear differential operators are self-adjoint with respect to this

norm if

p, q

are real and the boundary terms vanish. When do the boundary

terms vanish? One possibility is when

is periodic (with the right period), or if

we constrain f and g to be periodic.

Example. We can consider a simple case, where

L =

Here we have p = 1, q = 0. If we ask for functions to be periodic on [a, b], then

∗

dx =

∗

g dx.

Note that it is important that we have a second-order differential operator. If it

is first-order, then we would have a negative sign, since we integrated by parts

once.

Just as in finite dimensions, self-adjoint operators have eigenfunctions and

eigenvalues with special properties. First, we define a more sophisticated inner

product.

Definition

(Inner product with weight)

An inner product with weight

written ( ·, ·)

, is defined by

(f, g)

∗

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

where w is real, non-negative, and has only finitely many zeroes.

Why do we want a weight

(

)? In the future, we might want to work with

the unit disk, instead of a square in

. When we want to use polar coordinates,

we will have to integrate with

, instead of just d

. Hence we need the

weight of

. Also, we allow it to have finitely many zeroes, so that the radius

can be 0 at the origin.

Why can’t we have more zeroes? We want the inner product to keep the

property that (

f, f

)

= 0 iff

= 0 (for continuous

). If

is zero at too many

places, then the inner product could be zero without f being zero.

We now define what it means to be an eigenfunction.

Definition

(Eigenfunction with weight)

An eigenfunction with weight

is a function y : [a, b] → C obeying the differential equation

Ly = λwy,

where λ ∈ C is the eigenvalue.

This might be strange at first sight. It seems like we can take any nonsense

, apply

, to get some nonsense

. But then it is fine, since we can write it

as some nonsense

times our original

. So any function is an eigenfunction?

No! There are many constraints

has to satisfy, like being positive, real and

having finitely many zeroes. It turns out this severely restraints what values

can take, so not everything will be an eigenfunction. In fact we can develop

this theory without the weight function

. However, weight functions are much

more convenient when, say, dealing with the unit disk.

Proposition. The eigenvalues of a Sturm-Liouville operator are real.

Proof. Suppose Ly

= λ

. Then

, y

)

= λ

, wy

) = (y

, Ly

) = (Ly

, y

) = (λ

, y

) = λ

∗

, y

)

Since (y

, y

)

6= 0, we have λ

= λ

∗

Note that the first and last terms use the weighted inner product, but the

middle terms use the unweighted inner product.

Proposition.

Eigenfunctions with different eigenvalues (but same weight) are

orthogonal.

Proof. Let Ly

= λ

and Ly

= λ

. Then

, y

)

= (y

, Ly

) = (Ly

, y

) = λ

, y

)

Since λ

6= λ

, we must have (y

, y

)

= 0.

Those were pretty straightforward manipulations. However, the main results

of Sturm–Liouville theory are significantly harder, and we will not prove them.

We shall just state them and explore some examples.

Theorem. On a compact domain, the eigenvalues λ

, λ

, ··· form a countably

infinite sequence and are discrete.

This will be a rather helpful result in quantum mechanics, since in quantum

mechanics, the possible values of, say, the energy are the eigenvalues of the

Hamiltonian operator. Then this result says that the possible values of the

energy are discrete and form an infinite sequence.

Note here the word compact. In quantum mechanics, if we restrict a particle

in a well [0

1], then it will have quantized energy level since the domain is

compact. However, if the particle is free, then it can have any energy at all since

we no longer have a compact domain. Similarly, angular momentum is quantized,

since it describe rotations, which takes values in S

, which is compact.

Theorem.

The eigenfunctions are complete: any function

: [

a, b

]

→ C

(obeying

appropriate boundary conditions) can be expanded as

f(x) =

(x),

where

∗

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

Example.

Let [

a, b

] = [

−L, L

= 1, restricting to periodic functions

Then our eigenfunction obeys

= λ

(x),

Then our eigenfunctions are

(x) = exp



inπx



with eigenvalues

= −

for n ∈ Z. This is just the Fourier series!

Example

(Hermite polynomials)

We are going to cheat a little bit and pick

our domain to be R. We want to study the differential equation

− xH

= λH,

with H : R → C. We want to put this in Sturm-Liouville form. We have

p(x) = exp



−

2t dt



= e

−x

ignoring constant factors. Then q(x) = 0. We can rewrite this as



−x



= −2λe

−x

H(x).

So we take our weight function to be w(x) = e

−x

We now ask that

(

) grows at most polynomially as

|x| → ∞

. In particular,

we want

−x

(

)

→

0. This ensures that the boundary terms from integration

by parts vanish at the infinite boundary, so that our Sturm-Liouville operator is

self-adjoint.

The eigenfunctions turn out to be

(x) = (−1)



−x



These are known as the Hermite polynomials. Note that these are indeed

polynomials. When we differentiate the

−x

term many times, we get a lot

of things from the product rule, but they will always keep an

−x

, which will

ultimately cancel with e

Just as for matrices, we can use the eigenfunction expansion to solve forced

differential equations. For example, if might want to solve

Lg = f(x),

where f (x) is a forcing term. We can write this as

Lg = w(x)F (x).

We expand our g as

g(x) =

n∈Z

ˆg

(x).

Then by linearity,

Lg =

n∈Z

ˆg

(x) =

n∈Z

ˆg

w(x)y

(x).

We can also expand our forcing function as

w(x)F(x) = w(x)

n∈Z

(x)

Taking the (regular) inner product with

(

) (and noting orthogonality of

eigenfunctions), we obtain

w(x)ˆg

= w(x)

This tells us that

ˆg

So we have

g(x) =

n∈Z

(x),

provided all λ

are non-zero.

This is a systematic way of solving forced differential equations. We used to

solve these by “being smart”. We just looked at the forcing term and tried to

guess what would work. Unsurprisingly, this approach does not succeed all the

time. Thus it is helpful to have a systematic way of solving the equations.

It is often helpful to rewrite this into another form, using the fact that

= (y

, F )

. So we have

g(x) =

n∈Z

, F )

(x) =

n∈Z

∗

(t)y

(x)w(t)F(t) dt.

Note that we swapped the sum and the integral, which is in general a dangerous

thing to do, but we don’t really care because this is an applied course. We can

further write the above as

g(x) =

G(x, t)F (t)w(t) dt,

where G(x, t) is the infinite sum

G(x, t) =

n∈Z

∗

(t)y

(x).

We call this the Green’s function. Note that this depends on

and

only.

It depends on the differential operator

, but not the forcing term

. We can

think of this as something like the “inverse matrix”, which we can use to solve

the forced differential equation for any forcing term.

Recall that for a matrix, the inverse exists if the determinant is non-zero,

which is true if the eigenvalues are all non-zero. Similarly, here a necessary

condition for the Green’s function to exist is that all the eigenvalues are non-zero.

We now have a second version of Parseval’s theorem.

Theorem (Parseval’s theorem II).

(f, f)

n∈Z

Proof. We have

(f, f)

Ω

∗

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

n,m∈Z

Ω

∗

(x)

(x)w(x) dx

n,m∈Z

∗

, y

)

n∈Z

3.2 Least squares approximation

So far, we have expanded functions in terms of infinite series. However, in real

life, when we ask a computer to do this for us, it is incapable of storing and

calculating an infinite series. So we need to truncate it at some point.

Suppose we approximate some function

: Ω

→ C

by a finite set of eigen-

functions y

(x). Suppose we write the approximation g as

g(x) =

k=1

(x).

The objective here is to figure out what values of the coefficients

are “the

best”, i.e. can make g represent f as closely as possible.

Firstly, we have to make it explicit what we mean by “as closely as possible”.

Here we take this to mean that we want to minimize the

-norm (

f −g, f −g

)

By linearity of the norm, we know that

(f − g, f − g)

= (f, f)

+ (g, g)

− (f, g)

− (g, f)

To minimize this norm, we want

0 =

∂

∂c

(f − g, f − g)

∂

∂c

[(f, f)

+ (g, g)

− (f, g)

− (g, f)

We know that the (

f, f

)

term vanishes since it does not depend on

. Expanding

our definitions of g, we can get

0 =

∂

∂c

∞

i=1

−

k=1

∗

−

k=1

∗

= c

∗

−

∗

Note that here we are treating

∗

and

as distinct quantities. So when we vary

∗

is unchanged. To formally justify this treatment, we can vary the real and

imaginary parts separately.

Hence, the extremum is achieved at

∗

. Similarly, varying with respect

to c

∗

, we get that c

To check that this is indeed an minimum, we can look at the second-derivatives

∂

∂c

(f − g, f − g)

∂

∂c

∗

(f − g, f − g)

= 0,

while

∂

∂c

∗

∂c

(f − g, f − g)

= δ

≥ 0.

Hence this is indeed a minimum.

Thus we know that (f − g, f − g)

is minimized over all g(x) when

= (y

, f)

These are exactly the coefficients in our infinite expansion. Hence if we truncate

our series at an arbitrary point, it is still the best approximation we can get.

4 Partial differential equations

So far, we have come up with a lot of general theory about functions and stuff.

We are going to look into some applications of these. In this section, we will

develop the general technique of separation of variables, and apply this to many

many problems.

4.1 Laplace’s equation

Definition

(Laplace’s equation)

Laplace’s equation on

says that a (twice-

differentiable) equation φ : R

→ C obeys

∇

φ = 0,

where

∇

i=1

∂

∂x

This equation arises in many situations. For example, if we have a conservative

force

−∇φ

, and we are in a region where the force obeys the source-free

Gauss’ law

∇·F

= 0 (e.g. when there is no electric charge in an electromagnetic

field), then we get Laplace’s equation.

It also arises as a special case of the heat equation

∂φ

∂t

κ∇

, and wave

equation

∂

∂t

= c

∇

φ, when there is no time dependence.

Definition

(Harmonic functions)

Functions that obey Laplace’s equation are

often called harmonic functions.

Proposition.

Let Ω be a compact domain, and

∂

Ω be its boundary. Let

f : ∂Ω → R. Then there is a unique solution to ∇

φ = 0 on Ω with φ|

∂Ω

= f.

Proof.

Suppose

and

are both solutions such that

∂Ω

. Then

let Φ = φ

− φ

. So Φ = 0 on the boundary. So we have

0 =

Ω

Φ∇

Φ d

x = −

Ω

(∇Φ) · (∇Φ) dx +

∂Ω

Φ∇Φ · n d

n−1

We note that the second term vanishes since Φ = 0 on the boundary. So we have

0 = −

Ω

(∇Φ) · (∇Φ) dx.

However, we know that (

∇

Φ)

(

∇

Φ)

≥

0 with equality iff

∇

Φ = 0. Hence

Φ is constant throughout Ω. Since at the boundary, Φ = 0, we have Φ = 0

throughout, i.e. φ

= φ

Asking that

φ|

∂Ω

(

) is a Dirichlet boundary condition. We can instead

ask

n ·∇φ|

∂Ω

, a Neumann condition. In this case, we have a unique solution

up to a constant factor.

These we’ve all seen in IA Vector Calculus, but this time in arbitrary dimen-

sions.

4.2 Laplace’s equation in the unit disk in R

Let Ω = {(x, y) ∈ R

: x

+ y

≤ 1}. Then Laplace’s equation becomes

0 = ∇

φ =

∂

∂x

∂

∂y

∂

∂z∂ ¯z

where we define z = x + iy, ¯z = x − iy. Then the general solution is

φ(z, ¯z) = ψ(z) + χ(¯z)

for some ψ, χ.

Suppose we want a solution obeying the Dirichlet boundary condition

(

z, ¯z

)

∂Ω

(

), where the boundary

∂

Ω is the unit circle. Since the do-

main is the unit circle S

, we can Fourier-expand it! We can write

f(θ) =

n∈Z

inθ

∞

n=1

inθ

∞

n=1

−n

−inθ

However, we know that

iθ

and

¯z

iθ

. So on the boundary, we know

that z = e

iθ

and ¯z = e

−iθ

. So we can write

f(θ) =

∞

n=1

(

−n

¯z

)

This is defined for |z| = 1. Now we can define φ on the unit disk by

φ(z, ¯z) =

∞

n=1

∞

n=1

−n

¯z

It is clear that this obeys the boundary conditions by construction. Also,

(

z, ¯z

)

certainly converges when

|z| <

1 if the series for

(

) converged on

∂

Ω. Also,

(

z, ¯z

) is of the form

(

) +

(

¯z

). So it is a solution to Laplace’s equation on

the unit disk. Since the solution is unique, we’re done!

4.3 Separation of variables

Unfortunately, the use of complex variables is very special to the case where

Ω ⊆ R

. In higher dimensions, we proceed differently.

We let Ω =

{

(

x, y, z

)

∈ R

: 0

≤ x ≤ a,

≤ y ≤ b, z ≥

}

. We want the

following boundary conditions:

ψ(0, y, z) = ψ(a, y, z) = 0

ψ(x, 0, z) = ψ(x, b, z) = 0

lim

z→∞

ψ(x, y, z) = 0

ψ(x, y, 0) = f(x, y),

where

is a given function. In other words, we want our function to be

when

z = 0 and vanish at the boundaries.

The first step is to look for a solution of ∇

ψ(x, y, z) = 0 of the form

ψ(x, y, z) = X(x)Y (y)Z(z).

Then we have

0 = ∇

ψ = Y ZX

+ XZY

+ XY Z

= XY Z





As long as ψ 6= 0, we can divide through by ψ and obtain

= 0.

The key point is each term depends on only one of the variables (

x, y, z

). If we

vary, say,

, then

does not change. For the total sum to be 0,

must

be constant. So each term has to be separately constant. We can thus write

= −λX, Y

= −µY, Z

= (λ + µ)Z.

The signs before λ and µ are there just to make our life easier later on.

We can solve these to obtain

X = a sin

√

λx + b cos

√

λx,

Y = c sin

√

λy + d cos

√

λx,

Z = g exp(

λ + µz) + h exp(−

λ + µz).

We now impose the homogeneous boundary conditions, i.e. the conditions that

ψ vanishes at the walls and at infinity.

– At x = 0, we need ψ(0, y, z) = 0. So b = 0.

– At x = a, we need ψ(a, y, z) = 0. So λ =



nπ



– At y = 0, we need ψ(x, 0, z) = 0. So d = 0.

– At y = b, we need ψ(x, b, z) = 0. So µ =



mπ



– As z → ∞, ψ(x, y, z) → 0. So g = 0.

Here we have n, m ∈ N.

So we found a solution

ψ(x, y, z) = A

n,m

sin



nπx



sin



mπy



exp(−s

n,m

z),

where A

n,m

is an arbitrary constant, and

n,m





This obeys the homogeneous boundary conditions for any

n, m ∈ N

but not the

inhomogeneous condition at z = 0.

By linearity, the general solution obeying the homogeneous boundary condi-

tions is

ψ(x, y, z) =

∞

n,m=1

n,m

sin



nπx



sin



mπy



exp(−s

n,m

The final step is to impose the inhomogeneous boundary condition

(

x, y,

0) =

f(x, y). In other words, we need

∞

n,m=1

n,m

sin



nπx



sin



mπy



= f(x, y). (†)

The objective is thus to find the A

n,m

coefficients.

We can use the orthogonality relation we’ve previously had:

sin



kπx



sin



nπx



dx = δ

k,n

So we multiply (†) by sin



kπx



and integrate:

∞

n,m=1

n,m

sin



kπx



sin



nπx



dx sin



mπy



sin



kπx



f(x, y) dx.

Using the orthogonality relation, we have

∞

m=1

k,m

sin



mπy



sin



kπx



f(x, y) dx.

We can perform the same trick again, and obtain

k,j

[0,a]×[0,b]

sin



kπx



sin



jπy



f(x, y) dx dy. (∗)

So we found the solution

ψ(x, y, z) =

∞

n,m=1

n,m

sin



nπx



sin



mπy



exp(−s

n,m

where

m,n

is given by (

∗

). Since we have shown that there is a unique solution

to Laplace’s equation obeying Dirichlet boundary conditions, we’re done.

Note that if we imposed a boundary condition at finite

, say 0

≤ z ≤ c

then both sets of exponential

exp

(

√

λ + µz

) would have contributed. Similarly,

does not vanish at the other boundaries, then the

cos

terms would also

contribute.

In general, to actually find our

, we have to do the horrible integral for

m,n

, and this is not always easy (since integrals are hard ). We will just look at

a simple example.

Example. Suppose f(x, y) = 1. Then

m,n

sin



nπx



sin



mπy



dy =

(

n, m both odd

0 otherwise

Hence we have

ψ(x, y, z) =

n,m odd

sin



nπx



sin



mπy



exp(−s

m,n

z).

4.4 Laplace’s equation in spherical polar coordinates

4.4.1 Laplace’s equation in spherical polar coordinates

We’ll often also be interested in taking

Ω = {(s, y, z) ∈ R

+ y

+ z

≤ a}.

Since our domain is spherical, it makes sense to use some coordinate system

with spherical symmetry. We use spherical coordinates (r, θ, φ), where

x = r sin θ cos φ, y = r sin θ sin φ, z = r cos θ.

The Laplacian is

∇

∂

∂r



∂

∂r



sin θ

∂

∂θ



sin θ

∂

∂θ



sin

∂

∂φ

Similarly, the volume element is

dV = dx dy dz = r

sin θ dr dθ dφ.

In this course, we’ll only look at axisymmetric solutions, where

(

r, θ, φ

) =

(

r, θ

)

does not depend on φ.

We perform separation of variables, where we look for any set of solution of

∇

ψ = 0 inside Ω where ψ(r, θ) = R(r)Θ(θ).

Laplace’s equation then becomes



) +

Θ sin θ

dθ



sin θ

dΘ

dθ



= 0.

Similarly, since each term depends only on one variable, each must separately

be constant. So we have





= λR,

dθ



sin θ

dΘ

dθ



= −λ sin θΘ.

Note that both these equations are eigenfunction equations for some Sturm-

Liouville operator and the Θ(

) equation has a non-trivial weight function

w(θ) = sin θ.

4.4.2 Legendre Polynomials

The angular equation is known as the Legendre’s equation. It’s standard to set

cos θ

(which is not the Cartesian coordinate, but just a new variable). Then

we have

dθ

= −sin θ

So the Legendre’s equation becomes

−sin θ



sin θ(−sin θ)

dΘ



+ λ sin θΘ = 0.

Equivalently, in terms of x, we have



(1 − x

)

dΘ



= −λΘ.

Note that since 0 ≤ θ ≤ π, the domain of x is −1 ≤ x ≤ 1.

This operator is a Sturm–Liouville operator with

p(x) = 1 − x

, q(x) = 0.

For the Sturm–Liouville operator to be self-adjoint, we had

(g, Lf ) = (Lg, f) + [p(x)(g

∗

− g

0∗

f)]

−1

We want the boundary term to vanish. Since

(

) = 1

− x

vanishes at our

boundary

1, the Sturm-Liouville operator is self-adjoint provided our

function Θ(

) is regular at

1. Hence we look for a set of Legendre’s

equations inside (

−

1) that remains regular including at

1. We can try a

power series

Θ(x) =

∞

n=0

The Legendre’s equation becomes

(1 − x

)

∞

n=0

n(n − 1)x

n−2

− 2

∞

n=0

+ λ

∞

n=0

= 0.

For this to hold for all

x ∈

(

−

1), the equation must hold for each coefficient

of x separately. So

(λ − n(n + 1)) + a

n−2

(n + 2)(n + 1) = 0.

This requires that

n+2

n(n + 1) −λ

(n + 2)(n + 1)

This equation relates

n+2

. So we can choose

and

freely. So we get

two linearly independents solutions Θ

(x) and Θ

(x), where they satisfy

(−x) = Θ

(x), Θ

(−x) = −Θ

(x).

In particular, we can expand our recurrence formula to obtain

(x) = a



1 −

−

(6 − λ)λ

+ ···



(x) = a



x +

(2 − λ)

(12 − λ)(2 −λ)

+ ···



We now consider the boundary conditions. We know that Θ(

) must be regular

(i.e. finite) at x = ±1. This seems innocent, However, we have

lim

n→∞

n+2

= 1.

This is fine if we are inside (

−

1), since the power series will still converge.

However, at

1, the ratio test is not conclusive and does not guarantee conver-

gence. In fact, more sophisticated convergence tests show that the infinite series

would diverge on the boundary! To avoid this problem, we need to choose

such

that the series truncates. That is, if we set

(

+ 1) for some

` ∈ N

, then

our power series will truncate. Then since we are left with a finite polynomial,

it of course converges.

In this case, the finiteness boundary condition fixes our possible values of

eigenvalues. This is how quantization occurs in quantum mechanics. In this case,

this process is why angular momentum is quantized.

The resulting polynomial solutions

(

) are called Legendre polynomials.

For example, we have

(x) = 1

(x) = x

(x) =

(3x

− 1)

(x) =

(5x

− 3x),

where the overall normalization is chosen to fix P

(1) = 1.

It turns out that

(x) =

− 1)

The constants in front are just to ensure normalization. We now check that this

indeed gives the desired normalization:

(1) =

[(x − 1)

(x + 1)

]



x=1



`!(x + 1)

+ (x − 1)(stuff)



x=1

= 1

We have previously shown that the eigenfunctions of Sturm-Liouville operators

are orthogonal. Let’s see that directly now for this operator. The weight function

is w(x) = 1.

We first prove a short lemma: for 0 ≤ k ≤ `, we have

− 1)

= Q

`,k

− 1)

`−k

for some degree k polynomial Q

`,k

(x).

We show this by induction. This is trivially true when k = 0. We have

k+1

+ 1)



`,k

(x)(x

− 1)

`−k



= Q

`,k

− 1)

`−k

+ 2(` − k)Q

`,k

x(x

− 1)

`−k−1

= (x

− 1)

`−k−1

`,k

− 1) − 2(` − k)Q

`,k

Since

`,k

has degree

`,k

has degree

k −

1. So the right bunch of stuff has

degree k + 1.

Now we can show orthogonality. We have

, P

) =

−1

(x)P

(x) dx

−1

− 1)

(x) dx



`−1

− 1)

(x)



−

−1

`−1

− 1)

= −

−1

`−1

− 1)

dx.

Note that the first term disappears since the (

` −

derivative of (

−

still

has a factor of

−

1. So integration by parts allows us to transfer the derivative

from x

− 1 to P

Now if

` 6

, we can wlog assume

< `

. Then we can integrate by parts

times until we get the `

0th

derivative of P

, which is zero.

In fact, we can show that

, P

) =

2δ

2` + 1

hence the P

(x) are orthogonal polynomials.

By the fundamental theorem of algebra, we know that

(

) has

roots. In

fact, they are always real, and lie in (−1, 1).

Suppose not. Suppose only

m < `

roots lie in (

−

1). Then let

(

) =

r=1

(

x−x

), where

, x

, ··· , x

}

are these

roots. Consider the polynomial

(

)

(

). If we factorize this, we get

r=m+1

(

x−r

)

r=1

(

x−x

)

. The first

few terms have roots outside (

−

1), and hence do not change sign in (

−

1).

The latter terms are always non-negative. Hence for some appropriate sign, we

have

−1

(x)Q

(x) dx > 0.

However, we can expand

(x) =

r=1

(x),

but (P

, P

) = 0 for all r < `. This is a contradiction.

4.4.3 Solution to radial part

In our original equation

∇

= 0, we now set Θ(

) =

(

cos θ

). The radial

equation becomes

)

= `(` + 1)R.

Trying a solution of the form

(

) =

, we find that

(

+ 1) =

(

+ 1). So

the solution is α = ` or −(` + 1). Hence

φ(r, θ) =



`+1



(cos θ)

is the general solution to the Laplace equation

∇

= 0 on the domain. In we

want regularity at r = 0, then we need b

= 0 for all `.

The a

are fixed if we require φ(a, θ) = f(θ) on ∂Ω, i.e. when r = 1.

We expand

f(θ) =

∞

`=0

(cos θ),

where

= (P

, f) =

(cos θ)f(θ) d(cos θ).

φ(r, θ) =

∞

`=0





(cos θ).

4.5 Multipole expansions for Laplace’s eq uation

One can quickly check that

φ(r) =

|r − r

solves Laplace’s equation

∇

= 0 for all

r ∈ R

\ r

. For example, if

where

k is the unit vector in the z direction, then

|r −

√

+ 1 − 2r cos θ

∞

`=0

(cos θ).

To find these coefficients, we can employ a little trick. Since

(1) = 0, at

= 0,

we have

∞

`=0

1 − r

∞

`=0

So all the coefficients must be 1. This is valid for r < 1.

More generally, we have

|r − r

∞

`=0





(

r ·

This is called the multiple expansion, and is valid when r < r

. Thus

|r − r

r ·

+ ··· .

The first term

is known as the monopole, and the second term is known as

the dipole, in analogy to charges in electromagnetism. The monopole term is

what we get due to a single charge, and the second term is the result of the

interaction of two charges.

4.6 Laplace’s equation in cylindrical coordinates

We let

Ω = {(r, θ, z) ∈ R

: r ≤ a, z ≥ 0}.

In cylindrical polars, we have

∇

φ =

∂

∂r



∂φ

∂r



∂

∂θ

∂

∂z

= 0.

We look for a solution that is regular inside Ω and obeys the boundary conditions

φ(a, θ, z) = 0

φ(r, θ, 0) = f(r, θ)

lim

z→∞

φ(r, θ, z) = 0,

where f is some fixed function.

Again, we start by separation of variables. We try

φ(r, θ, z) = R(r)Θ(θ)Z(z).

Then we get

(rR

)

= 0.

So we immediately know that

= µZ

We now replace Z

/Z by µ, and multiply by r

to obtain

(rR

)

+ µr

= 0.

So we know that

= −λΘ.

Then we are left with

+ rR

+ (µr

− λ)R = 0.

Since we require periodicity in θ, we have

Θ(θ) = a

sin nθ + b

cos nθ, λ = n

, n ∈ N.

Since we want the solution to decay as z → ∞, we have

Z(z) = c

−

√

µz

The radial equation is of Sturm-Liouville type since it can be written as





−

R = −µrR.

Here we have

p(r) = r, q(r) = −

, w(r) = r.

Introducing x = r

√

µ, we can rewrite this as

+ x

+ (x

− n

)R = 0.

This is Bessel’s equation. Note that this is actually a whole family of differential

equations, one for each

. The

is not the eigenvalue we are trying to figure out

in this equation. It is already fixed by the Θ equation. The actual eigenvalue we

are finding out is µ, not n.

Since Bessel’s equations are second-order, there are two independent solution

for each n, namely J

(x) and Y

(x). These are called Bessel functions of order

of the first (

) or second (

) kind. We will not study these functions’

properties, but just state some useful properties of them.

The

’s are all regular at the origin. In particular, as

x →

(

)

∼ x

These look like decaying sine waves, but the zeroes are not regularly spaced.

On the other hand, the

’s are similar but singular at the origin. As

x →

(x) ∼ ln x, while Y

(x) ∼ x

−n

Now we can write our separable solution as

φ(r, θ, z) = (a

sin nθ + b

cos nθ)e

−

√

µz

µ,n

√

µ) + d

µ,n

(r,

√

µ)].

Since we want regularity at

= 0, we don’t want to have the

terms. So

µ,n

= 0.

We now impose our first boundary condition

(

a, θ, z

) = 0. This demands

√

µ) = 0. So we must have

√

µ =

where

is the

th root of

(

) (since there is not much pattern to these roots,

this is the best description we can give!).

So our general solution obeying the homogeneous conditions is

φ(r, θ, z) =

∞

n=0

i∈roots

sin nθ + B

cos nθ) exp



−







We finally impose the inhomogeneous boundary condition

(

r, θ,

0) =

(

r, θ

). To

do this, we use the orthogonality condition









r dr =

)]

which is proved in the example sheets. Note that we have a weight function

here. Also, remember this is the orthogonality relation for different roots of

Bessel’s functions of the same order. It does not relate Bessel’s functions of

different orders.

We now integrate our general solution with respect to cos mθ to obtain

−π

f(r, θ) cos mθ dθ =

i∈roots





Now we just have Bessel’s function for a single order

. So we can use the

orthogonality relation for the Bessel’s functions to obtain

)]

πa

−π

cos mθJ





f(r, θ)r dr dθ.

How can we possibly perform this integral? We don’t know how

looks like,

what the roots

are, and we multiply all these complicated things together

and integrate! Often, we just don’t. We ask our computers to do this numerically

for us. There are indeed some rare cases where we are indeed able to get an

explicit solution, but these are rare.

4.7 The heat equation

We choose Ω =

, ∞

). We treat

as the “space” variable, and [0

, ∞

) as

our time variable.

Definition (Heat equation). The heat equation for a function φ : Ω → C is

∂φ

∂t

= κ∇

φ,

where κ > 0 is the diffusion constant.

We can think of this as a function

representing the “temperature” at

different points in space and time, and this equation says the rate of change of

temperature is proportional to ∇

φ.

Our previous study of the Laplace’s equation comes in handy. When our

system is in equilibrium and does not change with time, then we are left with

Laplace’s equation.

Let’s first look at some interesting properties of the heat equation.

The first most important property of the heat equation is that the total

“heat” is conserved under evolution by the heat equation, i.e.

φ(x, t) d

x = 0,

since we can write this as

φ(x, t) d

x =

∂φ

∂t

x = κ

∇

φ d

x = 0,

provided ∇φ(x, t) → 0 as |x| → ∞.

In particular, on a compact space, i.e. Ω =

M ×

, ∞

) with

compact,

there is no boundary term and we have

φ dµ = 0.

The second useful property is if

(

x, t

) solves the heat equation, so do

(

x, t

) =

φ(x − x

, t − t

), and φ

(x, t) = Aφ(λx, λ

t).

Let’s try to choose A such that

(x, t) d

x =

φ(x, t) d

We have

(x, t) d

x = A

φ(λx, λ

t) d

x = Aλ

−n

φ(y, λ

t) d

where we substituted y = λx.

So if we set

, then the total heat in

at time

is the same as the

total heat in

at time

. But the total heat is constant in time. So the total

heat is the same all the time.

Note again that if

∂φ

∂t

κ∇

on Ω

, ∞

), then so does

(

λx, λ

). This

means that nothing is lost by letting

√

tκ

and consider solutions of the form

φ(x, t) =

(κt)

n/2



√

κt



(κt)

n/2

F (η),

where

η =

√

κt

For example, in 1 + 1 dimensions, we can look for a solution of the form

√

κt

(

We have

∂φ

∂t

∂

∂t



√

κt

F (η)



−1

√

κt

F (η) +

√

κt

dη

(η) =

−1

√

κt

[F + ηF

]

On the other side of the equation, we have

∂

∂x



√

κt

F (η)



√

κt

∂

∂x



∂η

∂x



√

κt

So the equation becomes

0 = 2F

+ ηF

+ F = (2F

+ ηF )

So we have

+ ηF = const,

and the boundary conditions require that the constant is zero. So we get

−η

We now see a solution

F (η) = a exp



−η



By convention, we now normalize our solution φ(x, t) by requiring

∞

−∞

φ(x, t) dx = 1.

This gives

a =

√

4π

and the full solution is

φ(x, t) =

√

4πκt

exp



−

4κt



More generally, on R

× [0, ∞), we have the solution

φ(x, t) =

(4π(t − t

))

n/2

exp



−(x − x

)

4κ(t − t

)



At any fixed time

t 6

, the solutions are Gaussians centered on

with standard

deviation

2κ(t − t

), and the height of the peak at x = x

4πκ(t−t

)

We see that the solution gets “flatter” as

increases. This is a general property

of the evolution by the heat equation. Suppose we have an eigenfunction

(

) of

the Laplacian on Ω, i.e.

∇

−λy

. We want to show that in certain cases,

must be positive. Consider

−λ

Ω

|y|

x =

Ω

∗

(x)∇

y(x) d

x =

∂Ω

∗

n · ∇y d

n−1

x −

Ω

|∇y|

If there is no boundary term, e.g. when Ω is closed and compact (e.g. a sphere),

then we get

λ =

|∇y|

Ω

|y|

Hence the eigenvalues are all positive. What does this mean for heat flow?

Suppose

(

x, t

) obeys the heat equation for

t >

0. At time

= 0, we have

φ(x, 0) = f(x). We can write the solution (somehow) as

φ(x, t) = e

t∇

φ(x, 0),

where we (for convenience) let κ = 1. But we can expand our initial condition

φ(x, 0) = f(x) =

(x)

in a complete set {y

(x)} of eigenfunctions for ∇

By linearity, we have

φ(x, t) = e

t∇

(x)

t∇

(x)

−λ

(x)

= c

(t)y

(x),

where

(

) =

−λ

. Since

are all positive, we know that

(

) decays

exponentially with time. In particular, the coefficients corresponding to the

largest eigenvalues

|λ

decays most rapidly. We also know that eigenfunctions

with larger eigenvalues are in general more “spiky”. So these smooth out rather

quickly.

When people first came up with the heat equation to describe heat loss, they

were slightly skeptical — this flattening out caused by the heat equation is not

time-reversible, while Newton’s laws of motions are all time reversible. How

could this smoothening out arise from reversible equations?

Einstein came up with an example where the heat equation can come out of

apparently naturally from seemingly reversible laws, namely Brownian motion.

The idea is that a particle will randomly jump around in space, where the

movement is independent of time and position.

Let the probability that a dust particle moves through a step

over a time

∆t be p(y, ∆t). For any fixed ∆t, we must have

∞

−∞

p(y, ∆t) dy = 1.

We also assume that

(

∆

) is independent of time, and of the location of the

dust particle. We also assume

(

∆

) is strongly peaked around

= 0 and

(

∆

) =

(

−y,

∆

). Now let

(

x, t

) be the probability that the dust particle

is located at x at time t. Then at time t + δt, we have

P (x, t + ∆t) =

∞

−∞

P (x − y, t)p(y, ∆t) dy.

For P (x −y, t) sufficiently smooth, we can write

P (x − y, t) ≈ P (x, t) − y

∂P

∂x

(x, t) +

∂

∂x

(x, t).

So we get

P (x, t + ∆t) ≈ P (x, t) −

∂P

∂x

(x, t)hyi+

∂

∂x

P (x, t) + ··· ,

where

i =

∞

−∞

p(y, ∆t) dy.

Since

is even, we expect

to vanish when

is odd. Also, since

is strongly

peaked at 0, we expect the higher order terms to be small. So we can write

P (x, t + ∆t) −P (x, t) =

∂

∂x

Suppose that as we take the limit ∆

t →

0, we get

2∆t

→ κ

for some

. Then

this becomes the heat equation

∂P

∂t

= κ

∂

∂x

Proposition.

Suppose

: Ω

, ∞

)

→ R

satisfies the heat equation

∂φ

∂t

κ∇

and obeys

– Initial conditions φ(x, 0) = f(x) for all x ∈ Ω

– Boundary condition φ(x, t)|

∂Ω

= g(x, t) for all t ∈ [0, ∞).

Then φ(x, t) is unique.

Proof. Suppose φ

and φ

are both solutions. Then define Φ = φ

− φ

and

E(t) =

Ω

dV.

Then we know that

(

)

≥

0. Since

, φ

both obey the heat equation, so does

Φ. Also, on the boundary and at t = 0, we know that Φ = 0. We have

Ω

dΦ

= κ

Ω

Φ∇

Φ dV

= κ

∂Ω

Φ∇Φ · dS −κ

Ω

(∇Φ)

= −κ

Ω

(∇Φ)

≤ 0.

So we know that

decreases with time but is always non-negative. We also

know that at time

= 0,

= Φ = 0. So

= 0 always. So Φ = 0 always. So

= φ

Example

(Heat conduction in uniform medium)

Suppose we are on Earth,

and the Sun heats up the surface of the Earth through sun light. So the sun will

maintain the soil at some fixed temperature. However, this temperature varies

with time as we move through the day-night cycle and seasons of the year.

We let

(

x, t

) be the temperature of the soil as a function of the depth

defined on R

≥0

× [0, ∞). Then it obeys the heat equation with conditions

(i) φ(0, t) = φ

+ A cos



2πt



+ B cos



2πt



(ii) φ(x, t) → const as x → ∞.

We know that

∂φ

∂t

= K

∂

∂x

We try the separation of variables

φ(x, t) = T (t)X(x).

Then we get the equations

= λT, X

From the boundary solution, we know that our things will be oscillatory. So we

let λ be imaginary, and set λ = iω. So we have

φ(x, t) = e

iωt



−

√

iω/Kx

+ b

√

iω/Kx



Note that we have

iω







(1 + i)

|ω|

ω > 0

(i − 1)

|ω|

ω < 0

Since

(

x, t

)

→

constant as

x → ∞

, we don’t want our

to blow up. So if

ω <

we need a

= 0. Otherwise, we need b

= 0.

To match up at x = 0, we just want terms with

|ω| = ω

2π

, |ω| = ω

2π

So we can write out solution as

φ(x, t) = φ

+ A exp



−



cos



t −



+ B exp



−



cos



t −



We can notice a few things. Firstly, as we go further down, the effect of the sun

decays, and the temperature is more stable. Also, the effect of the day-night

cycle decays more quickly than the annual cycle, which makes sense. We also

see that while the temperature does fluctuate with the same frequency as the

day-night/annual cycle, as we go down, there is a phase shift. This is helpful

since we can store things underground and make them cool in summer, warm in

winter.

Example

(Heat conduction in a finite rod)

We now have a finite rod of length

2L.

x = 0

x = −L x = L

We have the initial conditions

φ(x, 0) = Θ(x) =

(

1 0 < x < L

0 −L < x < 0

and the boundary conditions

φ(−L, t) = 0, φ(L, t) = 1

So we start with a step temperature, and then maintain the two ends at fixed

temperatures 0 and 1.

We are going to do separation of variables, but we note that all our boundary

and initial conditions are inhomogeneous. This is not helpful. So we use a

little trick. We first look for any solution satisfying the boundary conditions

(

−L, t

) = 0,

(

L, t

) = 1. For example, we can look for time-independent

solutions φ

(x, t) = φ

(x). Then we need

= 0. So we get

(x) =

x + L

By linearity,

(

x, t

) =

(

x, t

)

−φ

(

) obeys the heat equation with the conditions

ψ(−L, t) = ψ(L, t) = 0,

which is homogeneous! Our initial condition now becomes

ψ(x, 0) = Θ(x) −

x + L

We now perform separation of variables for

ψ(x, t) = X(x)T (t).

Then we obtain the equations

= −κλT, X

= −λX.

Then we have

ψ(x, t) =

a sin(

√

λx) + b cos(

√

λx)

−κλt

Since initial condition is odd, we can eliminate all

cos

terms. Our boundary

conditions also requires

λ =

, n = 1, 2, ··· .

So we have

φ(x, t) = φ

(x) +

∞

n=1

sin



nπx



exp



−

κn



where a

are the Fourier coefficients

−L



Θ(x) −

x + L



sin



nπx



dx =

nπ

Example

(Cooling of a uniform sphere)

Once upon a time, Charles Darwin

went around the Earth, looked at species, and decided that evolution happened.

When he came up with his theory, it worked quite well, except that there was

one worry. Was there enough time on Earth for evolution to occur?

This calculation was done by Kelvin. He knew well that the Earth started as

a ball of molten rock, and obviously life couldn’t have evolved when the world

was still molten rock. So he would want to know how long it took for the Earth

to cool down to its current temperature, and if that was sufficient for evolution

to occur.

We can, unsurprisingly, use the heat equation. We assume that at time

= 0,

the temperature is

, the melting temperature of rock. We also assume that

the space is cold, and we let the temperature on the boundary of Earth as 0. We

then solve the heat equation on a sphere (or ball), and see how much time it

takes for Earth to cool down to its present temperature.

We let Ω =

{

(

x, y, z

)

∈ R

, r ≤ R}

, and we want a solution

(

r, t

) of the heat

equation that is spherically symmetric and obeys the conditions

– φ(R, t) = 0

– φ(r, 0) = φ

We can write the heat equation as

∂φ

∂t

= κ∇

φ =

∂

∂r



∂φ

∂r



Again, we do the separation of variables.

φ(r, t) = R(r)T (t).

So we get

= −λ

κT,





= −λ

We can simplify this a bit by letting

(

) =

S(r)

, then our radial equation just

becomes

= −λ

We can solve this to get

R(r) = A

sin λr

+ B

cos λr

We want a regular solution at

= 0. So we need

= 0. Also, the boundary

condition φ(R, t) = 0 gives

λ =

nπ

, n = 1, 2, ···

So we get

φ(r, t) =

n∈Z

sin



nπr



exp



−κn



where our coefficients are

= (−1)

n+1

nπ

We know that the Earth isn’t just a cold piece of rock. There are still volcanoes.

So we know many terms still contribute to the sum nowadays. This is rather

difficult to work with. So we instead look at the temperature gradient

∂φ

∂r

n∈Z

(−1)

n+1

cos



nπr



exp



−κn



+ sin stuff.

We evaluate this at the surface of the Earth, R = r. So we get the gradient

∂φ

∂r



= −

n∈Z

exp



−

κn



≈

∞

−∞

exp



−

κπy



dy = φ

πκt

So the age of the Earth is approximately

t ≈





4πK

where

V =

∂φ

∂r



We can plug the numbers in, and get that the Earth is 100 million years. This is

not enough time for evolution.

Later on, people discovered that fission of metals inside the core of Earth

produce heat, and provides an alternative source of heat. So problem solved!

The current estimate of the age of the universe is around 4 billion years, and

evolution did have enough time to occur.

4.8 The wave equation

Consider a string x ∈ [0, L] undergoing small oscillations described by φ(x, t).

φ(x, t)

Consider two points

separated by a small distance

δx

. Let

(

) be

the outward pointing tension tangent to the string at

(

). Since there is no

sideways (x) motion, there is no net horizontal force. So

cos θ

= T

cos θ

= T. (∗)

If the string has mass per unit length

, then in the vertical direction, Newton’s

second law gives

µδx

∂

∂t

= T

sin θ

− T

sin θ

We now divide everything by T , noting the relation (∗), and get

δx

∂

∂t

sin θ

cos θ

−

sin θ

cos θ

= tan θ

− tan θ

∂φ

∂x



−

∂φ

∂x



≈ δx

∂

∂x

Taking the limit

δx →

0 and setting

T/µ

, we have that

(

x, t

) obeys the

wave equation

∂

∂t

∂

∂x

From IA Differential Equations, we’ve learnt that the solution can be written as

(

x − ct

) +

(

). However, this does not extend to higher dimensions, but

separation of variables does. So let’s do that.

Assume that the string is fixed at both ends. Then

, t

) =

(

L, t

) = 0 for

all

. The we can perform separation of variables, and the general solution can

then be written

φ(x, t) =

∞

n=1

sin



nπx





cos



nπct



+ B

sin



nπct



The coefficients

are fixed by the initial profile

(

0) of the string, while

the coefficients

are fixed by the initial string velocity

∂φ

∂t

(

0). Note that we

need two sets of initial conditions, since the wave equation is second-order in

time.

Energetics and uniqueness

An oscillating string contains has some sort of energy. The kinetic energy of a

small element δx of the string is

µδx



∂φ

∂t



The total kinetic energy of the string is hence the integral

K(t) =



∂φ

∂t



dx.

The string also has potential energy due to tension. The potential energy of a

small element δx of the string is

T × extension = T (δs − δx)

= T (

δx

+ δφ

− δx)

≈ T δx

1 +



δφ

δx



+ ···

− T δx

δx



δφ

δx



Hence the total potential energy of the string is

V (t) =



∂φ

∂x



dx,

using the definition of c

Using our series expansion for φ, we get

K(t) =

µπ

∞

n=1



sin



nπct



− B

cos



nπct



V (t) =

µπ

∞

n=1



cos



nπct



+ B

sin



nπct



The total energy is

E(t) =

nπ

∞

n=1

+ B

What can we see here? Our solution is essentially an (infinite) sum of independent

harmonic oscillators, one for each

. The period of the fundamental mode (

= 1)

2π

= 2

π ·

πc

. Thus, averaging over a period, the average kinetic energy

K =

2L/c

K(t) dt =

V =

2L/c

V (t) dt =

Hence we have an equipartition of the energy between the kinetic and potential

energy.

The energy also allows us to prove a uniqueness theorem. First we show that

energy is conserved in general.

Proposition.

Suppose

: Ω

, ∞

)

→ R

obeys the wave equation

∂

∂t

∇

inside Ω × (0, ∞), and is fixed at the boundary. Then E is constant.

Proof. By definition, we have

Ω

∂

∂t

∂ψ

∂t

+ c

∇



∂φ

∂t



· ∇φ dV.

We integrate by parts in the second term to obtain

Ω

dφ



∂

∂t

− c

∇



dV + c

∂Ω

∂φ

∂t

∇φ · dS.

Since

∂

∂t

∇

by the wave equation, and

is constant on

∂

Ω, we know

that

= 0.

Proposition.

Suppose

: Ω

, ∞

)

→ R

obeys the wave equation

∂

∂t

∇

inside Ω × (0, ∞), and obeys, for some f, g, h,

(i) φ(x, 0) = f (x);

(ii)

∂φ

∂t

(x, 0) = g(x); and

(iii) φ|

∂Ω×[0,∞)

= h(x).

Then φ is unique.

Proof.

Suppose

and

are two such solutions. Then

− φ

obeys the

wave equation

∂

∂t

= c

∇

ψ,

and

ψ|

∂Ω×[0,∞)

= ψ|

Ω×{0}

∂ψ

∂t



Ω×{0}

= 0.

We let

(t) =

Ω



∂ψ

∂t



+ c

∇ψ · ∇ψ

dV.

Then since

obeys the wave equation with fixed boundary conditions, we know

is constant.

Initially, at

= 0, we know that

∂ψ

∂t

= 0. So

(0) = 0. At time

, we

have

Ω



∂ψ

∂t



+ c

(∇ψ) · (∇ψ) dV = 0.

Hence we must have

∂ψ

∂t

= 0. So

is constant. Since it is 0 at the beginning, it

is always 0.

Example. Consider Ω = {(x, y) ∈ R

, x

+ y

≤ 1}, and let φ(r, θ, t) solve

∂

∂t

= ∇

φ =



∂φ

∂r



∂

∂θ

with the boundary condition

φ|

∂Ω

= 0. We can imagine this as a drum, where

the membrane can freely oscillate with the boundary fixed.

Separating variables with φ(r, θ, t) = T (t)R(r)Θ(θ), we get

= −c

λT, Θ

= −µΘ, r(R

)

+ (r

λ − µ)R = 0.

Then as before,

and Θ are both sine and cosine waves. Since we are in polars

coordinates, we need

(

t, r, θ

+ 2

) =

(

t, r, θ

). So we must have

for

some m ∈ N. Then the radial equation becomes

r(rR

)

+ (r

λ − m

)R = 0,

which is Bessel’s equation of order m. So we have

R(r) = a

(

√

λr) + b

(

√

λr).

Since we want regularity at

= 0, we need

= 0 for all

. To satisfy the

boundary condition

φ|

∂Ω

= 0, we must choose

√

, where

is the

root of J

Hence the general solution is

φ(t, r, θ) =

∞

i=0

sin(k

ct) + B

cos(k

ct)]J

∞

m=1

∞

i=0

cos(mθ) + B

sin(mθ)] sin k

ctJ

∞

m=1

∞

i=0

cos(mθ) + D

sin(mθ)] cos k

ctJ

For example, suppose we have the initial conditions

, r, θ

) = 0,

∂

, r, θ

) =

(

). So we start with a flat surface and suddenly hit it with some force. By

symmetry, we must have

, B

, C

, D

= 0 for

m 6

= 0. If this does not

seem obvious, we can perform some integrals with the orthogonality relations to

prove this.

At t = 0, we need φ = 0. So we must have B

= 0. So we are left with

φ =

∞

i=0

sin(k

ct)J

r).

We can differentiate this, multiply with J

r)r to obtain

∞

i=0

r)J

r)r dr =

g(r)J

r)r dr.

Using the orthogonality relations for J

from the example sheet, we get

)]

g(r)J

r)r dr.

Note that the frequencies come from the roots of the Bessel’s function, and are

not evenly spaced. This is different from, say, string instruments, where the

frequencies are evenly spaced. So drums sound differently from strings.

So far, we have used separation of variables to solve our differential equations.

This is quite good, as it worked in our examples, but there are a few issues with

it. Of course, we have the problem of whether it converges or not, but there is a

more fundamental problem.

To perform separation of variables, we need to pick a good coordinate system,

such that the boundary conditions come in a nice form. We were able to

perform separation of variables because we can find some coordinates that fit our

boundary conditions well. However, in real life, most things are not nice. Our

domain might have a weird shape, and we cannot easily find good coordinates

for it.

Hence, instead of solving things directly, we want to ask some more general

questions about them. In particular, Kac asked the question “can you hear the

shape of a drum?” — suppose we know all the frequencies of the modes of

oscillation on some domain Ω. Can we know what Ω is like?

The answer is no, and we can explicitly construct two distinct drums that

sound the same. Fortunately, we get an affirmative answer if we restrict ourselves

a bit. If we require Ω to be convex, and has a real analytic boundary, then yes!

In fact, we have the following result: let

(

) be the number of eigenvalues

less than λ

. Then we can show that

4π

lim

→∞

N(λ

)

= Area(Ω).

5 Distributions

5.1 Distributions

When performing separation of variables, we also have the problem of convergence

as well. To perform separation of variables, we first find some particular solutions

of the form, say,

(

)

(

)

(

). We know that these solve, say, the wave equation.

However, what we do next is take an infinite sum of these functions. First of

all, how can we be sure that this converges at all? Even if it did, how do we

know that the sum satisfies the differential equation? As we have seen in Fourier

series, an infinite sum of continuous functions can be discontinuous. If it is not

even continuous, how can we say it is a solution of a differential equation?

Hence, at first people were rather skeptical of this method. They thought

that while these methods worked, they only work on a small, restricted domain

of problems. However, later people realized that this method is indeed rather

general, as long as we allow some more “generalized” functions that are not

functions in the usual sense. These are known as distributions.

To define a distribution, we first pick a class of “nice” test functions, where

“nice” means we can do whatever we want to them (e.g. differentiate, integrate

etc.) A main example is the set of infinitely smooth functions with compact

support on some set

K ⊆

Ω, written

∞

cpt

(Ω) or

(Ω). For example, we can

have the bump function defined by

φ(x) =

(

−1/(1−x

)

|x| < 1

0 otherwise.

We now define a distribution

to be a linear map

(Ω)

→ R

. For those

who are doing IB Linear Algebra, this is the dual space of the space of (“nice”)

real-valued functions. For φ ∈ D(Ω), we often write its image under T as T [φ].

Example.

The simplest example is just an ordinary function that is integrable

over any compact region. Then we define the distribution T

[φ] =

Ω

f(x)φ(x) dx.

Note that this is a linear map since integration is linear (and multiplication

is commutative and distributes over addition). Hence every function “is” a

distribution.

Of course, this would not be interesting if we only had ordinary functions.

The most important example of a distribution (that is not an ordinary function)

is the Dirac delta “function”.

Definition (Dirac-delta). The Dirac-delta is a distribution defined by

δ[φ] = φ(0).

By analogy with the first case, we often abuse notation and write

δ[φ] =

Ω

δ(x)φ(x) dx,

and pretend

(

) is an actual function. Of course, it is not really a function, i.e.

it is not a map Ω

→ R

. If it were, then we must have

(

) = 0 whenever

x 6

= 0,

since

[

] =

(0) only depends on what happens at 0. But then this integral

will just give 0 if

(0)

∈ R

. Some people like to think of this as a function that

is zero anywhere and “infinitely large” everywhere else. Formally, though, we

have to think of it as a distribution.

Although distributions can be arbitrarily singular and insane, we can nonethe-

less define all their derivatives, defined by

[φ] = −T [φ

This is motivated by the case of regular functions, where we get, after integrating

by parts,

Ω

(x)φ(x) dx = −

Ω

f(x)φ

(x) dx,

with no boundary terms since we have a compact support. Since

is infinitely

differentiable, we can take arbitrary derivatives of distributions.

So even though distributions can be crazy and singular, everything can be

differentiated. This is good.

Generalized functions can occur as limits of sequences of normal functions.

For example, the family of functions

(x) =

√

exp(−n

)

are smooth for any finite n, and G

[φ] → δ[φ] for any φ.

It thus makes sense to define

(φ) = −δ[φ

] = −φ

(0),

as this is the limit of the sequence

lim

n→∞

Ω

(x)φ(x) dx.

It is often convenient to think of

(

) as

lim

n→∞

(

), and

(

) =

lim

n→∞

(

)

etc, despite the fact that these limits do not exist as functions.

We can look at some properties of δ(x):

– Translation:

∞

−∞

δ(x − a)φ(x) dx =

∞

−∞

δ(y)φ(y + a) dx.

– Scaling:

∞

−∞

δ(cx)φ(x) dx =

∞

−∞

δ(y)φ





|c|

φ(0).

–

These are both special cases of the following: suppose

(

) is a continuously

differentiable function with isolated simple zeros at

. Then near any of

its zeros x

, we have f(x) ≈ (x − x

)

∂f

∂x



. Then

∞

−∞

δ(f(x))φ(x) dx =

i=1

∞

−∞

(x − x

)

∂f

∂x



φ(x) dx

i=1

φ(x

We can also expand the

-function in a basis of eigenfunctions. Suppose we live

in the interval [−L, L], and write a Fourier expansion

δ(x) =

n∈Z

inπx/L

with

−L

inπx/L

δ(x) dx =

So we have

δ(x) =

n∈Z

inπx/L

This does make sense as a distribution. We let

δ(x) =

n=−N

inπx/L

Then

lim

N→∞

−L

δ(x)φ(x) dx = lim

N→∞

n=−N

inπx/L

φ(x) dx

= lim

N→∞

n=−N

−L

inπx/L

φ(x) dx

= lim

N→∞

n=−N

−n

= lim

N→∞

n=−N

inπ0/L

= φ(0),

since the Fourier series of the smooth function

(

) does converge for all

x ∈

[−L, L].

We can equally well expand

(

) in terms of any other set of orthonormal

eigenfunctions. Let

(

)

}

be a complete set of eigenfunctions on [

a, b

] that are

orthogonal with respect to a weight function w(x). Then we can write

δ(x − ξ) =

(x),

with

∗

(x)δ(x − ξ)w(x) dx = y

∗

(ξ)w(ξ).

δ(x − ξ) = w(ξ)

∗

(ξ)y

(x) = w(x)

∗

(ξ)y

(x),

using the fact that

δ(x − ξ) =

w(ξ)

w(x)

δ(x − ξ).

5.2 Green’s functions

One of the main uses of the

function is the Green’s function. Suppose we wish

to solve the 2nd order ordinary differential equation

, where

(

) is a

bounded forcing term, and L is a bounded differential operator, say

L = α(x)

∂

∂x

+ β(x)

∂

∂x

+ γ(x).

As a first step, we try to find the Green’s function for

, which we shall write as

G(x, ξ), which obeys

LG(x, ξ) = δ(x − ξ).

Given G(x, ξ), we can then define

y(x) =

G(x, ξ)f(ξ) dξ.

Then we have

Ly =

LG(x, ξ)f(ξ) dξ =

δ(x − ξ)f(ξ) dξ = f(x).

We can formally say

−1

, and the Green’s function shows that

−1

means

an integral.

Hence if we can solve for Green’s function, then we have solved the differential

equation, at least in terms of an integral.

To find the Green’s function G(x, ξ) obeying the boundary conditions

G(a, ξ) = G(b, ξ) = 0,

first note that

= 0 for all

x ∈

[

a, ξ

)

∪

(

ξ, b

], i.e. everywhere except

itself.

Thus we must be able to expand it in a basis of solutions in these two regions.

Suppose

(

)

, y

(

)

}

are a basis of solutions to

= 0 everywhere on

[a, b], with boundary conditions y

(a) = 0, y

(b) = 0. Then we must have

G(x, ξ) =

(

A(ξ)y

(x) a ≤ x < ξ

B(ξ)y

(x) ξ < x ≤ b

So we have a whole family of solutions. To fix the coefficients, we must decide

how to join these solutions together over x = ξ.

(

x, ξ

) were discontinuous at

, then

∂

x=ξ

would involve a

function, while

∂

x=ξ

would involve the derivative of the

function. This is

not good, since nothing in

(

x − ξ

) would balance a

. So

(

x, ξ

) must

be everywhere continuous. Hence we require

A(ξ)y

(ξ) = B(ξ)y

(ξ). (∗)

Now integrate over a small region (ξ − ε, ξ + ε) surrounding ξ. Then we have

ξ+ε

ξ−ε



α(x)

+ β(x)

+ γ(x)G



dx =

ξ+ε

ξ−ε

δ(x − ξ) dx = 1.

By continuity of

, we know that the

γG

term does not contribute. While

discontinuous, it is still finite. So the

βG

term also does not contribute. So we

have

lim

ε→0

ξ+ε

ξ−ε

αG

dx = 1.

We now integrate by parts to obtain

lim

ε→0

[αG

]

ξ+ε

ξ−ε

ξ+ε

ξ−ε

dx = 1.

Again, by finiteness of G

, the integral does not contribute. So we know that

α(ξ)

∂G

∂x



−

∂G

∂x



−

= 1

Hence we obtain

B(ξ)y

(ξ) − A(ξ)y

(ξ) =

α(ξ)

Together with (∗), we know that

A(ξ) =

(ξ)

α(ξ)W (ξ)

, B(ξ) =

(ξ)

α(ξ)W (ξ)

where W is the Wronskian

W = y

− y

Hence, we know that

G(x, ξ) =

α(ξ)W (ξ)

(

(ξ)y

(x) a ≤ x ≤ ξ

(ξ)y

(x) ξ < x ≤ b.

Using the step function Θ, we can write this as

G(x, ξ) =

α(ξ)W (ξ)

[Θ(ξ − x)y

(ξ)y

(x) + Θ(x −ξ)y

(ξ)y

(x)].

So our general solution is

y(x) =

G(x, ξ)f(ξ) dξ

f(ξ)

α(ξ)W (ξ)

(ξ)y

(x) dξ +

f(ξ)

α(ξ)W (ξ)

(ξ)y

(x) dξ.

Example. Consider

Ly = −y

− y = f

for

x ∈

1) with

(0) =

(1) = 0. We choose our basis solution to satisfy

= −y as {sin x, sin(1 − x)}. Then we can compute the Wronskian

W (x) = −sin x cos(1 − x) −sin(1 − x) cos x = −sin 1.

Our Green’s function is

G(x, ξ) =

sin 1

[Θ(ξ − x) sin(1 −ξ) sin x + Θ(x −ξ) sin ξ sin(1 − x)].

Hence we get

y(x) = sin x

sin(1 − ξ)

sin 1

f(ξ) dξ + sin(1 − x)

sin ξ

sin 1

f(ξ) dξ.

Example.

Suppose we have a string with mass per unit length

(

) and held

under a tension T . In the presence of gravity, Newton’s second law gives

∂

∂t

= T

∂

∂x

+ µg.

We look for the steady state solution

˙y

= 0 shape of the string, assuming

y(0, t) = y(L, t) = 0. So we get

∂

∂x

= −

µ(x)

We look for a Green’s function obeying

∂

∂x

= −δ(x − ξ).

This can be interpreted as the contribution of a pointlike mass located at

0 L

The homogeneous equation y

= 0 gives

y = Ax + B(x − L).

So we get

G(x, ξ) =

(

A(ξx) 0 ≤ x < ξ

B(ξ)(x − L) ξ < x ≤ L.

Continuity at x = ξ gives

A(ξ)ξ = B(ξ)(ξ − L).

The jump condition on the derivative gives

A(ξ) − B(ξ) = 1.

We can solve these to get

A(ξ) =

ξ − L

, B(ξ) =

Hence the Green’s function is

G(x, ξ) =

ξ − L

Θ(ξ − x) +

(x − L)Θ(x −ξ).

Notice that

is always less that

. So the first term has a negative slope; while

ξ is always positive, so the second term has positive slope.

For the general case, we can just substitute this into the integral. We can

give this integral a physical interpretation. We can think that we have many

pointlike particles m

at small separations ∆x along the string. So we get

g(x) =

i=1

G(x, x

)

and in the limit ∆x → 0 and n → ∞, we get

g(x) →

G(x, ξ)µ(ξ) dξ.

5.3 Green’s functions for initial value problems

Consider

(

) subject to

(

) = 0 and

(

) = 0. So instead of

having two boundaries, we have one boundary and restrict both the value and

the derivative at this boundary.

As before, let

(

)

, y

(

) be any basis of solutions to

= 0. The Green’s

function obeys

L(G) = δ(t − τ ).

We can write our Green’s function as

G(t, τ) =

(

A(τ)y

(t) + B(τ )y

(t) t

≤ t < τ

C(τ)y

(t) + D(τ)y

(t) t > τ.

Our initial conditions require that



) y

)

) y

)









However, we know that the matrix is non-singular (by definition of

, y

being a

basis). So we must have

A, B

= 0 (which makes sense, since

= 0 is obviously

a solution for t

≤ t < τ ).

Our continuity and jump conditions now require

0 = C(τ)t

(t) + D(τ)y

(t)

α(τ)

= C(τ)y

(τ) + D(τ )y

(τ).

We can solve this to get C(τ) and D(τ). Then we get

y(t) =

∞

G(t, τ)f(τ) dτ =

G(t, τ)f(τ) dτ,

since we know that when

τ > t

, the Green’s function

(

t, τ

) = 0. Thus the

solution

(

) depends on the forcing term term

only for times

< t

. This

expresses causality!

Example.

For example, suppose we have

¨y

(

) with

(0) =

˙y

(0) = 0.

Then we have

G(t, τ) = Θ(t − τ )[C(τ) cos(t − τ ) + D(τ) sin(t − τ )]

for some

(

)

, D

(

). The continuity and jump conditions gives

(

) = 1

, C

(

) =

0. So we get

G(t, τ) = Θ(t − τ ) sin t(τ).

So we get

y(t) =

sin(t − τ )f(τ) dτ.

6 Fourier transforms

6.1 The Fourier transform

So far, we have considered writing a function

→ C

(or a periodic function)

as a Fourier sum

f(θ) =

n∈Z

inθ

2π

−π

−inθ

f(θ) dθ.

What if we don’t have a periodic function, but just an arbitrary function

f : R → C? Can we still do Fourier analysis? Yes!

To start with, instead of

2π

−π

−inθ

(

) d

, we define the Fourier

transform as follows:

Definition (Fourier transform). The Fourier transform of an (absolutely inte-

grable) function f : R → C is defined as

f(k) =

∞

−∞

−ikx

f(x) dx

for all k ∈ R. We will also write

f(k) = F[f(x)].

Note that for any k, we have

f(k)| =



∞

−∞

−ikx

f(x) dx



≤

∞

−∞

−ikx

f(x)| dx =

∞

−∞

|f(x)| dx.

Since we have assumed that our function is absolutely integrable, this is finite,

and the definition makes sense.

We can look at a few basic properties of the Fourier transform.

(i)

Linearity: if

f, g

R → C

are absolutely integrable and

, c

are constants,

then

F[c

f(x) + c

g(x)] = c

F[f(x)] + c

F[g(x)].

This follows directly from the linearity of the integral.

(ii) Translation:

F[f(x − a)] =

−ikx

f(x − a) dx

Let y = x − a. Then we have

−ik(y+a)

f(y) dy

= e

−ika

−iky

f(y) dy

= e

−ika

F[f(x)],

So a translation in x-space becomes a re-phasing in k-space.

(iii) Re-phasing:

F[e

−i`x

f(x)] =

∞

−∞

−ikx

−i`x

f(x) dx

∞

−∞

−i(k+`)x

f(x) dx

f(k + `),

where ` ∈ R and

f(x) = F[f(x)].

(iv) Scaling:

F[f(cx)] =

∞

−∞

−ikx

f(cx) dx

∞

−∞

−iky/c

f(y)

|c|





Note that we have to take the absolute value of

, since if we replace

y/c

and

is negative, then we will flip the bounds of the integral. The

extra minus sign then turns c into −c = |c|.

(v) Convolutions: for functions f, g : R → C, we define the convolution as

f ∗ g(x) =

∞

−∞

f(x − y)g(y) dy.

We then have

F[f ∗ g(x)] =

∞

−∞

−ikx



∞

−∞

f(x − y)g(y) dy



ik(x−y)

f(x − y)e

−iky

g(y) dy dx

−iku

f(u) du

−iky

g(y) dy

= F[f]F[g],

where

x − y

. So the Fourier transform of a convolution is the product

of individual Fourier transforms.

(vi)

The most useful property of the Fourier transform is that it “turns differ-

entiation into multiplication”. Integrating by parts, we have

F[f

(x)] =

∞

−∞

−ikx

∞

−∞

−

−ikx

)f(x) dx

Note that we don’t have any boundary terms since for the function to be

absolutely integrable, it has to decay to zero as we go to infinity. So we

have

= ik

∞

−∞

−ikx

f(x) dx

= ikF[f(x)]

Conversely,

F[xf(x)] =

∞

−∞

−ikx

xf(x) dx = i

∞

−∞

−ikx

f(x) dx = i

(k).

This are useful results, since if we have a complicated function to find the Fourier

transform of, we can use these to break it apart into smaller pieces and hopefully

make our lives easier.

Suppose we have a differential equation

L(∂)y = f,

where

L(∂) =

r=0

is a differential operator of pth order with constant coefficients.

Taking the Fourier transform of both sides of the equation, we find

F[L(∂)y] = F[f (x)] =

f(k).

The interesting part is the left hand side, since the Fourier transform turns

differentiation into multiplication. The left hand side becomes

˜y(k) + c

ik˜y(k) + c

(ik)

˜y(k) + ··· + c

(ik)

˜y(k) = L(ik)˜y(k).

Here

(

) is a polynomial in

. Thus taking the Fourier transform has changed

our ordinary differential equation into the algebraic equation

L(ik)˜y(k) =

f(k).

Since L(ik) is just multiplication by a polynomial, we can immediately get

˜y(k) =

f(k)

L(ik)

Example. For φ : R

→ C, suppose we have the equation

∇

φ − m

φ = ρ(x),

where

∇

i=1

is the Laplacian on R

We define the n dimensional Fourier transform by

φ(k) =

−ik·x

f(x) dx = F[φ(x)]

where now k ∈ R

is an n-dimensional vector. So we get

F[∇

φ − m

φ] = F[ρ] = ˜ρ(k).

The first term is

−ik·x

∇ · ∇φ d

x = −

∇(e

−ik·x

) · ∇φ d

∇

−ik·x

)φ d

= −k · k

−ik·x

φ(x) d

= −k · k

φ(k).

So our equation becomes

−k · k

φ(k) − m

φ(k) = ˜ρ(k).

So we get

φ(k) = −

˜ρ(k)

|k|

+ m

So differential equations are trivial in

space. The problem is, of course,

that we want our solution in

space, not

space. So we need to find a way to

restore our function back to x space.

6.2 The Fourier inversion theorem

We need to be able to express our original function

(

) in terms of its Fourier

transform

(

). Recall that in the periodic case where

(

) =

(

), we have

f(x) =

n∈Z

2inxπ/L

We can try to obtain a similar expression for a non-periodic function

R → C

by taking the limit L → ∞. Recall that our coefficients are defined by

L/2

−L/2

−2inπx/L

f(x) dx.

Now define

∆k =

2π

So we have

f(x) =

n∈Z

inx∆k

∆k

2π

L/2

−L/2

−inx∆k

f(u) du

As we take the limit L → ∞, we can approximate

L/2

−L/2

−ixn∆k

f(x) dx =

∞

−∞

−ix(n∆k)

f(x) dx =

f(n∆k).

Then for a non-periodic function, we have

f(x) = lim

∆k→∞

n∈Z

∆k

2π

in∆kx

f(n∆k) =

2π

∞

−∞

ikx

f(k) dk.

f(x) = F

−1

[f(k)] =

2π

∞

−∞

ikx

f(k) dk.

This is a really dodgy argument. We first took the limit as

L → ∞

to turn

the integral from

−L

∞

−∞

. Then we take the limit as ∆

k →

0. However,

we can’t really do this, since ∆

is defined to be 2

π/L

, and both limits have

to be taken at the same time. So we should really just take this as a heuristic

argument for why this works, instead of a formal proof.

Nevertheless, note that the inverse Fourier transform looks very similar to

the Fourier transform itself. The only differences are a factor of

2π

, and the sign

of the exponent. We can get rid of the first difference by evenly distributing the

factor among the Fourier transform and inverse Fourier transform. For example,

we can instead define

F[f(x)] =

√

2π

∞

−∞

−ikx

f(x) dx.

Then F

−1

looks just like F apart from having the wrong sign.

However, we will stick to our original definition so that we don’t have to deal

with messy square roots. Hence we have the duality

F[f(x)] = F

−1

[f(−x)](2π).

This is useful because it means we can use our knowledge of Fourier transforms

to compute inverse Fourier transform.

Note that this does not occur in the case of the Fourier series. In the Fourier

series, we obtain the coefficients by evaluating an integral, and restore the

original function by taking a discrete sum. These operations are not symmetric.

Example. Recall that we defined the convolution as

f ∗ g(x) =

∞

−∞

f(x − y)g(y) dy.

We then computed

F[f ∗ g(x)] =

f(k)

f(g).

It follows by definition of the inverse that

−1

[

f(k)˜g(k)] = f ∗ g(x).

Therefore, we know that

F[f(x)g(x)] = 2π

∞

−∞

f(k − `)˜g(`) d` = 2π

f ∗ ˜g(k).

Example. Last time we saw that if ∇

φ − m

φ = ρ(x), then

φ(k) = −

˜ρ(k)

|k|

+ m

Hence we can retrieve our φ as

φ(x) = F

−1

[

φ(k)] = −

(2π)

ik·x

˜ρ(k)

|k|

+ m

Note that we have (2

)

instead of 2

since we get a factor of 2

for each

dimension (and the negative sign was just brought down from the original

expression).

Since we are taking the inverse Fourier transform of a product, we know that

the result is just a convolution of F

−1

[˜ρ(k)] = ρ(x) and F

−1



|k|



Recall that when we worked with Green’s function, if we have a forcing term,

then the final solution is the convolution of our Green’s function with the forcing

term. Here we get a similar expression in higher dimensions.

6.3 Parseval’s theorem for Fourier transforms

Theorem

(Parseval’s theorem (again))

Suppose

f, g

R → C

are sufficiently

well-behaved that

and

˜g

exist and we indeed have

−1

[

] =

f, F

−1

[

˜g

] =

Then

(f, g) =

∗

(x)g(x) dx =

2π

(

f, ˜g).

In particular, if f = g, then

kfk

2π

So taking the Fourier transform preserves the

norm (up to a constant

factor of

2π

Proof.

(f, g) =

∗

(x)g(x) dx

∞

−∞

∗

(x)



2π

∞

−∞

ikx

˜g(x) dk



2π

∞

−∞



∞

−∞

∗

(x)e

ikx



˜g(k) dk

2π

∞

−∞



∞

−∞

f(x)e

−ikx



∗

˜g(k) dk

2π

∞

−∞

∗

(k)˜g(k) dk

2π

(

f, ˜g).

6.4 A word of warning

Suppose f(x) is defined by

f(x) =

(

1 |x| < 1

0 |x| ≥ 1.

This function looks rather innocent. Sure it has discontinuities at

1, but these

are not too absurd, and we are just integrating. This is certainly absolutely

integrable. We can easily compute the Fourier transform as

f(k) =

∞

−∞

−ikx

f(x) dx =

−1

−ikx

dx =

2 sin k

Is our original function equal to the integral F

−1

[

f(l)]? We have

−1



2 sin k



∞

−∞

ikx

sin k

dk.

This is hard to do. So let’s first see if this function is absolutely integrable. We

have

∞

−∞



ikx

sin k



dx =

∞

−∞



sin k



= 2

∞



sin k



≥ 2

n=0

(n+3/4)π

(n+1/4)π



sin k



dk.

The idea here is instead of looking at the integral over the whole real line, we

just pick out segments of it. In these small segments, we know that

|sin k| ≥

√

So we can bound this by

∞

−∞



ikx

sin k



dx ≥ 2

n=0

√

(n+3/4)π

(n+1/4)π

≥

√

n=0

(n+3/4)π

(n+1/4)π

(n + 1)π

√

n=0

n + 1

and this diverges as

N → ∞

. So we don’t have absolute integrability. So we

have to be careful when we are doing these things. Well, in theory. We actually

won’t. We’ll just ignore this issue.

6.5 Fourier transformation of distributions

We have

F[δ(x)] =

∞

−∞

ikx

δ(x) dx = 1.

Hence we have

−1

[1] =

2π

∞

−∞

ikx

dk = δ(x).

Of course, it is extremely hard to make sense of this integral. It quite obviously

diverges as a normal integral, but we can just have faith and believe this makes

sense as long as we are talking about distributions.

Similarly, from our rules of translations, we get

F[δ(x − a)] = e

−ika

and

F[e

−i`x

] = 2πδ(k −`),

Hence we get

F[cos(`x)] = F



i`x

+ e

−i`x

)



F[e

i`x

F[e

−i`x

] = π[δ(k−`)+δ(k+`)].

We see that highly localized functions in

-space have very spread-out behaviour

in k-space, and vice versa.

6.6 Linear systems and response f unctions

Suppose we have an amplifier that modifies an input signal

(

) to produce an

output

(

). Typically, amplifiers work by modifying the amplitudes and phases

of specific frequencies in the output. By Fourier’s inversion theorem, we know

I(t) =

2π

iωt

I(ω) dω.

This

I(ω) is the resolution of I(t).

We specify what the amplifier does by the transfer function

(

) such that

the output is given by

O(t) =

2π

∞

−∞

iωt

R(ω)

I(ω) dω.

Since this

(

)

(

) is a product, on computing

(

) =

−1

[

(

)

(

)], we

obtain a convolution

O(t) =

∞

−∞

R(t − u)I(u) du,

where

R(t) =

2π

∞

−∞

iωt

R(ω) dω

is the response function. By plugging it directly into the equation above, we see

that R(t) is the response to an input signal I(t) = δ(t).

Now note that causality implies that the amplifier cannot “respond” before

any input has been given. So we must have R(t) = 0 for all t < 0.

Assume that we only start providing input at t = 0. Then

O(t) =

∞

−∞

R(t − u)I(u) du =

R(t − u)I(u) du.

This is exactly the same form of solution as we found for initial value problems

with the response function R(t) playing the role of the Green’s function.

6.7 General form of transfer function

To model the situation, suppose the amplifier’s operation is described by the

ordinary differential equation

I(t) = L

O(t),

where

i=0

with

∈ C

. In other words, we have an

th order ordinary differential equation

with constant coefficients, and the input is the forcing function. Notice that

we are not saying that

(

) can be expressed as a function of derivatives of

(

). Instead,

(

) influences

(

) by being the forcing term the differential

equation has to satisfy. This makes sense because when we send an input into

the amplifier, this input wave “pushes” the amplifier, and the amplifier is forced

to react in some way to match the input.

Taking the Fourier transform and using F





= iω

O(ω), we have

I(ω) =

j=0

(iω)

O(ω).

So we get

O(ω) =

I(ω)

+ iωa

+ ··· + (iω)

Hence, the transfer function is just

R(ω) =

+ iωa

+ ··· + (iω)

The denominator is an

th order polynomial. By the fundamental theorem of

algebra, it has

roots, say

∈ C

for

= 1

, ··· , J

, where

has multiplicity

. Then we can write our transfer function as

R(ω) =

j=1

(iω − c

)

By repeated use of partial fractions, we can write this as

R(ω) =

j=1

r=1

(iω − c

)

for some constants Γ

∈ C.

By linearity of the (inverse) Fourier transform, we can find

(

) if we know

the inverse Fourier transform of all functions of the form

(iω − α)

To compute this, there are many possible approaches. We can try to evaluate

the integral directly, learn some fancy tricks from complex analysis, or cheat,

and let the lecturer tell you the answer. Consider the function

(t) =

(

αt

t > 0

0 otherwise

The Fourier transform is then

(ω) =

∞

−∞

−iωt

(t) dt =

∞

(α−iω)t

dt =

iω − α

provided Re(α) < 0 (so that e

(α−iω)t

→ 0 as t → ∞). Similarly, we can let

(t) =

(

αt

t > 0

0 otherwise

Then we get

(ω) = F[th

(t)] = i

dω

F[h

(t)] =

(iω − α)

Proceeding inductively, we can define

(t) =

(

αt

t > 0

0 otherwise

and this has Fourier transform

(ω) =

(iω − α)

p+1

again provided

(

)

0. So the response function is a linear combination

of these functions

(

) (if any of the roots

have non-negative real part,

then it turns out the system is unstable, and is better analysed by the Laplace

transform). We see that the response function does indeed vanish for all

t <

In fact, each term (except

) increases from zero at

= 0 to rise to some

maximum before eventually decaying exponentially.

Fourier transforms can also be used to solve ordinary differential equations

on the whole real line

, provided the functions are sufficiently well-behaved.

For example, suppose y : R → R solves

− A

y = −f(x)

with y and y

→ 0 as |x| → ∞. Taking the Fourier transform, we get

˜y(k) =

f(k)

+ A

Since this is a product, the inverse Fourier transform will be a convolution of

(

)

and an object whose Fourier transform is

. Again, we’ll cheat. Consider

h(x) =

2µ

−µ|x|

with µ > 0. Then

f(k) =

2µ

∞

−∞

−ikx

−µ|x|

∞

−(µ+ik)x

ik + µ

+ k

Therefore we get

y(x) =

∞

−∞

−A|x−u|

f(u) du.

So far, we have done Fourier analysis over some abelian groups. For example,

we’ve done it over

, which is an abelian group under multiplication, and

which is an abelian group under addition. We will next look at Fourier analysis

over another abelian group, Z

, known as the discrete Fourier transform.

6.8 The discrete Fourier transform

Recall that the Fourier transform is defined as

f(k) =

−ikx

f(x) dx.

To find the Fourier transform, we have to know how to perform the integral. If

we cannot do the integral, then we cannot do the Fourier transform. However,

it is usually difficult to perform integrals for even slightly more complicated

functions.

A more serious problem is that we usually don’t even have a closed form of

the function

for us to perform the integral. In real life,

might be some radio

signal, and we are just given some data of what value

takes at different points..

There is no hope that we can do the integral exactly.

Hence, we first give ourselves a simplifying assumption. We suppose that

is mostly concentrated in some finite region [

−R, S

]. More precisely,

(

)

| 

for x 6∈ [−R, S] for some R, S > 0. Then we can approximate our integral as

f(k) =

−R

−ikx

f(x) dx.

Afterwards, we perform the integral numerically. Suppose we “sample”

(

) at

x = x

= −R + j

R + S

for N a large integer and j = 0, 1, ··· , N − 1. Then

f(k) ≈

R + S

N−1

j=0

f(x

−ikx

This is just the Riemann sum.

Similarly, our computer can only store the result

(

) for some finite list of

k. Let’s choose these to be at

k = k

2πm

R + S

Then after some cancellation,

f(k

) ≈

R + S

2πimR

R+S

N−1

j=0

f(x

2πi

= (R + S)e

2πimR

R+S





N−1

j=0

f(x

)ω

−jm





where

ω = e

2πi

is an N th root of unity. We call the thing in the brackets

F (m) =

N−1

j=0

f(x

)ω

−jm

Of course, we’ve thrown away lots of information about our original function

(

), since we made approximations all over the place. For example, we have

already lost all knowledge of structures varying more rapidly than our sampling

interval

R+S

. Also,

(

) =

(

), since

= 1. So we have “forced” some

periodicity into our function F , while

f(k) itself was not necessarily periodic.

For the usual Fourier transform, we were able to re-construct our original

function from the

(

), but here we clearly cannot. However, if we know the

(

) for all

= 0

, ··· , N −

1, then we can reconstruct the exact values of

(

) for all

by just solving linear equations. To make this more precise, we

want to put what we’re doing into the linear algebra framework we’ve been using.

Let

{

, ω, ω

, ··· , ω

N−1

}

. For our purposes below, we can just treat this as

a discrete set, and

are just meaningless symbols that happen to visually look

like the powers of our Nth roots of unity.

Consider a function g : G → C defined by

g(ω

) = f(x

This is actually nothing but just a new notation for

(

). Then using this new

notation, we have

F (m) =

N−1

j=0

−jm

g(ω

The space of all functions

G → C

is a finite-dimensional vector space,

isomorphic to C

. This has an inner product

(f, g) =

N−1

j=0

∗

(ω

)g(ω

This inner product obeys

(g, f ) = (f, g)

∗

(f, c

+ c

) = c

(f, g

) + c

(f, g

)

(f, f) ≥ 0 with equality iff f(ω

) = 0

Now let e

: G → C be the function

(ω

) = ω

Then the set of functions

}

for

= 0

, ··· , N −

1 is an orthonormal basis

with respect to our inner product. To show this, we can compute

, e

) =

N−1

j=0

∗

(ω

)

N−1

j=0

−jm

N−1

j=0

= 1.

For n 6= m, we have

, e

) =

N−1

j=0

∗

(ω

)

N−1

j=0

j(m−n)

1 − ω

(m−n)N

1 − ω

m−n

Since

m−n

is an integer and

is an

th root of unity, we know that

(m−n)N

= 1.

So the numerator is 0. However, since

n 6

m −n 6

= 0. So the denominator is

non-zero. So we get 0. Hence

, e

) = 0.

We can now rewrite our F (m) as

F (m) =

N−1

j=0

−jm

f(x

) =

N−1

j=0

∗

(ω

)g(ω

) = (e

, g).

Hence we can expand our g as

g =

N−1

m=0

, g)e

N−1

m=0

F (m)e

Writing f instead of g, we recover the formula

f(x

) = g(ω

) =

N−1

m=0

F (m)e

(ω

If we forget about our

s and just look at the

, what we have effectively done is

take the Fourier transform of functions taking values on

{

, ω, ··· , ω

N−1

}

∼

This is exactly analogous to what we did for the Fourier transform, except

that everything is now discrete, and we don’t have to worry about convergence

since these are finite sums.

6.9 The fast Fourier transform*

What we’ve said so far is we’ve defined

F (m) =

N−1

j=0

−mj

g(ω

To compute this directly, even given the values of

−mj

for all

j, m

, this takes

N −

1 additions and

+ 1 multiplications. This is 2

operations for a single

value of m. To know F (m) for all m, this takes 2N

operations.

This is a problem. Historically, during the cold war, especially after the

Cuban missile crisis, people were in fear that the world will one day go into

a huge nuclear war, and something had to be done to prevent this. Thus the

countries decided to come up with a treaty to ensure people don’t perform

nuclear testings anymore.

However, the obvious problem is that we can’t easily check if other countries

are doing nuclear tests. Nuclear tests are easy to spot if the test is done above

ground, since there will be a huge explosion. However, if it is done underground,

it is much harder to test.

They then came up with the idea to put some seismometers all around

the world, and record vibrations in the ground. To distinguish normal crust

movements from nuclear tests, they wanted to take the Fourier transform and

analyze the frequencies of the vibrations. However, they had a large amount of

data, and the value of

is on the order of magnitude of a few million. So 2

will be a huge number that the computers at that time were not able to handle.

So they need to find a trick to perform Fourier transforms quickly. This is known

as the fast Fourier transform. Note that this is nothing new mathematically.

This is entirely a computational trick.

Now suppose N = 2M . We can write

F (m) =

2M−1

j=0

−jm

g(ω

)

M−1

k=0

−2km

g(ω

) + ω

−(2k+1)m

g(ω

2k+1

)

We now let η be an M th root of unity, and define

G(η

) = g(ω

), H(η

) = g(ω

2k+1

We then have

F (m) =

M−1

k=0

−km

G(η

) +

−m

M−1

k=0

−km

H(η

)

[

G(m) + ω

−m

H(m)].

Suppose we are given the values of

(

) and

(

) and

−m

for all

{

, ··· , N −

}

. Then we can compute

(

) using 3 operations per value of

So this takes 3 ×N = 6M operations for all M.

We can compute

−m

for all

using at most 2

operations, and suppose it

takes

operations to compute

(

) for all

. Then the number of operations

needed to compute F (m) is

= 2P

+ 6M + 2M.

Now let N = 2

. Then by iterating this, we can find

≤ 4N log

N  N

So by breaking the Fourier transform apart, we are able to greatly reduce the

number of operations needed to compute the Fourier transform. This is used

nowadays everywhere for everything.

7 More partial differential equations

7.1 Well-posedness

Recall that to find a unique solution of Laplace’s equation

∇

= 0 on a bounded

domain Ω

⊆ R

, we imposed the Dirichlet boundary condition

φ|

∂Ω

or the

Neumann boundary condition n · ∇φ|

∂Ω

= g.

For the heat equation

∂φ

∂t

κ∇

on Ω

, ∞

), we asked for

∂|

∂Ω×[0,∞)

and also φ|

Ω×{0}

= g.

For the wave equation

∂

∂t

∇

on Ω

, ∞

), we imposed

φ|

∂Ω×[0,∞)

φ|

Ω×{0}

= g and also ∂

φ|

Ω×{0}

= h.

All the boundary and initial conditions restrict the value of

on some co-

dimension 1 surface, i.e. a surface whose dimension is one less than the domain.

This is known as Cauchy data for our partial differential equation.

Definition

(Well-posed problem)

A partial differential equation problem is

said to be well-posed if its Cauchy data means

(i) A solution exists;

(ii) The solution is unique;

(iii)

A “small change” in the Cauchy data leads to a “small change” in the

solution.

To understand the last condition, we need to make it clear what we mean by

“small change”. To do this properly, we need to impose some topology on our

space of functions, which is some technicalities we will not go into. Instead, we

can look at a simple example.

Suppose we have the heat equation

∂

κ∇

. We know that whatever

starting condition we have, the solution quickly smooths out with time. Any

spikiness of the initial conditions get exponentially suppressed. Hence this is

a well-posed problem — changing the initial condition slightly will result in a

similar solution.

However, if we take the heat equation but run it backwards in time, we get

a non-well-posed problem. If we provide a tiny, small change in the “ending

condition”, as we go back in time, this perturbation grows exponentially, and

the result could vary wildly.

Another example is as follows: consider the Laplace’s equation

∇

on the

upper half plane (x, y) ∈ R ×R

≥0

subject to the boundary conditions

φ(x, 0) = 0, ∂

φ(x, 0) = g(x).

If we take g(x) = 0, then φ(x, y) = 0 is the unique solution, obviously.

However, if we take g(x) =

sin Ax

, then we get the unique solution

φ(x, y) =

sin(Ax) sinh(Ay)

So far so good. However, now consider the limit as A → ∞. Then

g(x) =

sin(Ax)

→ 0.

for all x ∈ R. However, at the special point φ



, y



, we get



, y



sinh(Ay)

→ A

−2

which is unbounded. So as we take the limit as our boundary conditions

(

)

→

we get an unbounded solution.

How can we make sure this doesn’t happen? We’re first going to look at first

order equations.

7.2 Method of characteristics

Curves in the plane

Suppose we have a curve on the plane x : R → R

given by s 7→ (x(s), y(s)).

Definition

(Tangent vector)

The tangent vector to a smooth curve

given

by x : R → R

with x(s) = (x(s), y(s)) is





If we have some quantity

→ R

, then its value along the curve

is just

φ(x(s), y(s)).

A vector field

(

x, y

) :

→ R

defines a family of curves. To imagine this,

suppose we are living in a river. The vector field tells us the how the water flows

at each point. We can then obtain a curve by starting a point and flow along

with the water. More precisely,

Definition

(Integral curve)

Let

(

x, y

) :

→ R

be a vector field. The

integral curves associated to

are curves whose tangent





is just

(

x, y

For sufficiently regular vector fields, we can fill the space with different curves.

We will parametrize which curve we are on by the parameter t. More precisely,

we have a curve

= (

(

)

, y

(

)) that is transverse (i.e. nowhere parallel) to our

family of curves, and we can label the members of our family by the value of

at which they intersect B.

B(t)

We can thus label our family of curves by (x(s, t), y(s, t)). If the Jacobian

J =

∂x

∂s

∂y

∂t

−

∂x

∂t

∂y

∂s

6= 0,

then we can invert this to find (

s, t

) as a function of (

x, y

), i.e. at any point, we

know which curve we are on, and how far along the curve we are.

This means we now have a new coordinate system (

s, t

) for our points in

. It turns out by picking the right vector field

, we can make differential

equations much easier to solve in this new coordinate system.

The method of characteristics

Suppose φ : R

→ R obeys

a(x, y)

∂φ

∂x

+ b(x, y)

∂φ

∂y

= 0.

We can define our vector field as

V(x, y) =



a(x, y)

b(x, y)



Then we can write the differential equation as

V ·∇φ = 0.

Along any particular integral curve of V, we have

∂φ

∂s

dx(s)

∂φ

∂x

dy(s)

∂φ

∂y

= V ·∇φ,

where the integral curves of V are determined by

∂x

∂s



= a(x, y),

∂y

∂s



= b(x, y).

This is known as the characteristic equation.

Hence our partial differential equation just becomes the equation

∂φ

∂s



= 0.

To get ourselves a well-posed problem, we have to specify our boundary data

along a transverse curve

. We pick our transverse curve as

= 0, and we

suppose we are given the initial data

φ(0, t) = h(t).

Since φ does not vary with s, our solution automatically is

φ(s, t) = h(t).

Then φ(x, y) is given by inverting the characteristic equations to find t(x, y).

We do a few examples to make this clear.

Example (Trivial). Suppose φ obeys

∂φ

∂x



= 0.

We are given the boundary condition

φ(0, y) = f(y).

The solution is obvious, since the differential equation tells us the function is

just a function of

, and then the solution is obviously

(

x, y

) =

(

). However,

we can try to use the method of characteristics to practice our skills. Our vector

field is given by

V =









Hence, we get

= 1,

= 0.

So we have

x = s + c, y = d.

Our initial curve is given by

t, x

= 0. Since we want this to be the curve

s = 0, we find our integral curves to be

x = s, y = t.

Now for each fixed t, we have to solve

∂φ

∂s



= 0, φ(s, t) = h(t) = f (t).

So we know that

φ(x, y) = f(y).

Example. Consider the equation

∂φ

∂x

∂φ

∂y

= 0.

with the boundary condition

φ(x, 0) = cosh x.

We find that the vector field is

V =





This gives

= e

= 1.

Thus we have

−x

= −s + c, y = s + d.

We pick our curve as x = t, y = 0. So our relation becomes

−x

= −s + e

−t

, y = s.

We thus have

∂φ

∂s



= 0, φ(s, t) = φ(0, t) = cosh(t) = cosh[−ln(y + e

−x

)]

So done.

Example.

Let

→ R

solve the inhomogeneous partial differential equation

∂

φ + 2∂

φ = ye

with φ(x, x) = sin x.

We can still use the method of characteristics. We have

u =





So the characteristic curves obey

= 1,

= 2.

This gives

x = s + t, y = 2s + t

so that x = y = t at s = 0. We can invert this to obtain the relations

s = y − x, t = 2x −y.

The partial differential equation now becomes

dφ



= u · ∇φ = ye

= (2s + t)e

s+t

Note that this is just an ordinary differential equation in

because

is held

constant. We have the boundary conditions

φ(s = 0, t) = sin t.

So we get

φ(x(s, t), y(s, t)) = (2 −t)e

(1 − e

) + sin t + 2se

s+t

Putting it in terms of x and y, we have

φ(x, y) = (2 − 2x + y)e

2x−y

+ sin(2x − y) + (y − 2)e

We see that our Cauchy data

(

x, x

) =

sin x

should be specified on a curve

B that intersects each characteristic curve exactly once.

If we tried to use a characteristic curve

that intersects the characteristics

multiple times, if we want a solution at all, we cannot specify data freely along

. its values at the points of intersection have to compatible with the partial

differential equation.

For example, if we have a homogeneous equation, we saw the solution will

be constant along the same characteristic. Hence if our Cauchy data requires

the solution to take different values on the same characteristic, we will have no

solution.

Moreover, even if it is compatible, if

does not intersect all characteristics,

we will not get a unique solution. So the solution is not fixed along such

characteristics.

On the other hand, if

intersects each characteristic curve transversely,

the problem is well-posed. So there exists a unique solution, at least in the

neighbourhood of

. Notice that data is never transferred between characteristic

curves.

7.3 Characteristics for 2nd order partial differential equa-

tions

Whether the method of characteristics works for a 2nd order partial differential

equation, or even higher order ones, depends on the “type” of the differential

equation.

To classify these differential equations, we need to define

Definition

(Symbol and principal part)

Let

be the general 2nd order differ-

ential operator on R

. We can write it as

L =

i,j=1

(x)

∂

∂x

i=1

(x)

∂

∂x

+ c(X),

where a

(x), b

(x), c(x) ∈ R and a

= a

(wlog).

We define the symbol σ(k, x) of L to be

σ(k, x) =

i,j=1

(x)k

i=1

(x)k

+ c(x).

So we just replace the derivatives by the variable k.

The principal part of the symbol is the leading term

(k, x) =

i,j=1

(x)k

Example. If L = ∇

, then

σ(k, x) = σ

(k, x) =

i=1

)

is the heat operator (where we think of the last coordinate

as the time,

and others as space), then the operator is given by

L =

∂

∂x

−

n−1

i=1

∂

∂x

The symbol is then

σ(k, x) = k

−

n−1

i=1

)

and the principal part is

(k, x) = −

n−1

i=1

)

Note that the symbol is closely related to the Fourier transform of the

differential operator, since both turn differentiation into multiplication. Indeed,

they are equal if the coefficients are constant. However, we define this symbol

for arbitrary differential operators with non-constant coefficients.

In general, for each x, we can write

(k, x) = k

A(x)k,

where

(

) has elements

(

). Recall that a real symmetric matrix (such as

A) has all real eigenvalues. So we define the following:

Definition

(Elliptic, hyperbolic, ultra-hyperbolic and parabolic differential

operators). Let L be a differential operator. We say L is

– elliptic at x if all eigenvalues of A(x) have the same sign. Equivalently, if

( ·, x) is a definite quadratic form;

– hyperbolic at x if all but one eigenvalues of A(x) have the same sign;

– ultra-hyperbolic at x if A(x) has more than one eigenvalues of each sign;

– parabolic at x if A(x) has a zero eigenvalue, i.e. σ

( ·, x) is degenerate.

We say L is elliptic if L is elliptic at all x, and similarly for the other terms.

These names are supposed to remind us of conic sections, and indeed if we

think of k

Ak as an equation in k, then we get a conic section.

Example. Let

L = a(x, y)

∂

∂x

+ 2b(x, y)

∂

∂x∂y

+ c(x, y)

∂

∂y

+ d(x, y)

∂

∂x

+ e(x, y)

∂

∂y

+ f(x, y).

Then the principal part of the symbol is

(k, x) =







a(x, y) b(x, y)

b(x, y) c(x, y)





Then

is elliptic at

− ac <

0; hyperbolic if

− ac >

0; and parabolic if

− ac = 0.

Note that since we only have two dimensions and hence only two eigenvalues,

we cannot possibly have an ultra-hyperbolic equation.

Characteristic surfaces

Definition (Characteristic surface). Given a differential operator L, let

f(x

, x

, ··· , x

) = 0

define a surface C ⊆ R

. We say C is characteristic if

i,j=1

(x)

∂f

∂x

∂f

∂x

= (∇f)

A(∇f) = σ

(∇f, x) = 0.

In the case where we only have two dimensions, a characteristic surface is just a

curve.

We see that the characteristic equation restricts what values

∇f

can take.

Recall that

∇f

is the normal to the surface

. So in general, at any point, we

can find what the normal of the surface should be, and stitch these together to

form the full characteristic surfaces.

For an elliptic operator, all the eigenvalues of

have the same sign. So there

are no non-trivial real solutions to this equation (

∇f

)

(

∇f

) = 0. Consequently,

elliptic operators have no real characteristics. So the method of characteristics

would not be of any use when studying, say, Laplace’s equation, at least if we

want to stay in the realm of real numbers.

is parabolic, we for simplicity assume that

has exactly one zero

eigenvector, say

, and the other eigenvalues have the same sign. This is the

case when, say, there are just two dimensions. So we have

. We

normalize n such that n · n = 1.

For any ∇f, we can always decompose it as

∇f = n(n · ∇f) + [∇f − n · (n · ∇f)].

This is a relation that is trivially true, since we just add and subtract the same

thing. Note, however, that the first term points along

, while the latter term is

orthogonal to n. To save some writing, we adopt the notation

∇

⊥

f = ∇f − n(n · ∇f ).

So we have

∇f = n(n · ∇f) + ∇f

⊥

Then we can compute

(∇f)

A(∇f) = [n(n · ∇f ) + ∇

⊥

A[n(n · ∇f ) + ∇

⊥

= (∇

⊥

A(∇

⊥

f).

Then by assumption, (

∇

⊥

)

(

∇

⊥

) is definite. So just as in the elliptic case,

there are no non-trivial solutions. Hence, if

defines a characteristic surface,

then

∇

⊥

= 0. In other words,

is parallel to

. So at any point, the normal to

a characteristic surface must be

, and there is only one possible characteristic.

is hyperbolic, we assume all but one eigenvalues are positive, and let

−λ

be the unique negative eigenvalue. We let

be the corresponding unit

eigenvector, where we normalize it such that

n ·n

= 1. We say

is characteristic

0 = (∇f)

A(∇f)

= [n(n · ∇f ) + ∇

⊥

A[n(n · ∇f ) + ∇

⊥

= −λ(n · ∇f )

+ (∇

⊥

A(∇

⊥

f).

Consequently, for this to be a characteristic, we need

n · ∇f = ±

(∇

⊥

A(∇

⊥

So there are two choices for

n · ∇f

, given any

∇

⊥

. So hyperbolic equations

have two characteristic surfaces through any point.

This is not too helpful in general. However, in the case where we have two

dimensions, we can find the characteristic curves explicitly. Suppose our curve

is given by

(

x, y

) = 0. We can write

(

). Then since

is constant along

each characteristic, by the chain rule, we know

0 =

∂f

∂x

∂f

∂y

Hence, we can compute

= −

−b ±

√

− ac

We now see explicitly how the type of the differential equation influences the

number of characteristics — if

−ac >

0, then we obtain two distinct differential

equations and obtain two solutions; if

−ac

= 0, then we only have one equation;

if b

− ac < 0, then there are no real characteristics.

Example. Consider

∂

φ − xy∂

φ = 0

. Then

−xy, b

= 0

, c

= 1. So

− ac

. So the type is elliptic if

xy < 0, hyperbolic if xy > 0, and parabolic if xy = 0.

In the regions where it is hyperbolic, we find

−b ±

√

− ac

= ±

√

Hence the two characteristics are given by

= ±

√

This has a solution

3/2

± x

1/2

= c.

We now let

u =

3/2

+ x

1/2

v =

3/2

− x

1/2

Then the equation becomes

∂

∂u∂v

+ lower order terms = 0.

Example. Consider the wave equation

∂

∂t

− c

∂

∂x

= 0

1,1

. Then the equation is hyperbolic everywhere, and the characteristic

curves are

x ± ct

= const. Let’s look for a solution to the wave equation that

obeys

φ(x, 0) = f(x), ∂

φ(x, 0) = g(x).

Now put u = x − ct, v = x + ct. Then the wave equation becomes

∂

∂u∂v

= 0.

So the general solution to this is

φ(x, t) = G(u) + H(v) = G(x − ct) + H(x + ct).

The initial conditions now fix these functions

f(x) = G(x) + H(x), g(x) = −cG

(x) + cH

(x).

Solving these, we find

φ(x, t) =

[f(x − ct) + f(x + ct)] +

x+ct

x−ct

g(y) dy.

This is d’Alembert’s solution to the 1 + 1 dimensional wave equation.

Note that the value of

at any point (

x, t

) is completely determined by

f, g

in the interval [

x −ct, x

]. This is known as the (past) domain of dependence

of the solution at (

x, t

), written

−

(

x, t

). Similarly, at any time

, the initial

data at (

0) only affects the solution within the region

− ct ≤ x ≤ x

This is the range of influence of data at x

, written D

−

(p)

(S)

We see that disturbances in the wave equation propagate with speed c.

7.4 Green’s functions for PDEs on R

Green’s functions for the heat equation

Suppose φ : R

× [0, ∞) → R solves the heat equation

∂

φ = D∇

φ,

where D is the diffusion constant, subject to

φ|

×{0}

= f.

To solve this, we take the Fourier transform of the heat equation in the spatial

variables. For simplicity of notation, we take

= 1. Writing

[

(

x, t

)] =

(

k, t

we get

∂

φ(k, t) = −Dk

φ(k, t)

with the boundary conditions

φ(k, 0) =

f(k).

Note that this is just an ordinary differential equation in

, and

is a fixed

parameter for each k. So we have

φ(k, t) =

f(k)e

−Dk

We now take the inverse Fourier transform and get

φ(x, t) =

2π

ikx

f(k)e

−Dk

This is the inverse Fourier transform of a product. This thus gives the convolution

φ(x, t) = f ∗ F

−1

−Dk

So we are done if we can find the inverse Fourier transform of the Gaussian.

This is an easy exercise on example sheet 3, where we find

F[e

−a

] =

√

−k

/4a

So setting

4Dt

So we get

−1

−Dk

] =

√

4πDt

exp



−

4Dt



We shall call this

(

x, t

), where the subscript 1 tells us we are in 1+1 dimensions.

This is known as the fundamental solution of the heat equation. We then get

φ(x, t) =

∞

−∞

f(y)S

(x − y, t) dy

Example. Suppose our initial data is

f(x) = φ

δ(x).

So we start with a really cold room, with a huge spike in temperature at the

middle of the room. Then we get

φ(x, t) =

√

4πDt

exp



−

4Dt



What this shows, as we’ve seen before, is that if we start with a delta function,

then as time evolves, we get a Gaussian that gets shorter and broader.

Now note that if we start with a delta function, at

= 0, everywhere outside

the origin is non-zero. However, after any infinitely small time

becomes

non-zero everywhere, instantly. Unlike the wave equation, information travels

instantly to all of space, i.e. heat propagates arbitrarily fast according to this

equation (of course in reality, it doesn’t). This fundamentally, is because the

heat equation is parabolic, and only has one characteristic.

Now suppose instead φ satisfies the inhomogeneous, forced, heat equation

∂

φ − D∇

φ = F (x, t),

but with homogeneous initial conditions

φ|

t=0

= 0. Physically, this represents

having an external heat source somewhere, starting at zero.

Note that if we can solve this, then we have completely solved the heat

equation. If we have an inhomogeneous equation and an inhomogeneous initial

condition, then we can solve this forced problem with homogeneous boundary

conditions to get

; and solve the unforced equation with homogeneous equation

to get φ

. Then the sum φ = φ

+ φ

solves the full equation.

As before, we take the Fourier transform of the forced equation with respect

to the spacial variables. As before, we will just do it in the cases with one spacial

dimension. We find

∂

φ(k, t) + Dk

φ(k, t) =

F (k, t),

with the initial condition

φ(k, 0) = 0.

As before, we have reduced this to a first-order ordinary differential equation in

t. Using an integrating factor, we can rewrite this as

∂

∂t

φ(k, t)] = e

F (k, t).

The solution is them

φ(k, t) = e

−Dk

F (k, u) du.

We define the Green’s function G(x, t; y, τ )) to be the solution to

[∂

− D∇

]G(x, t; y, τ ) = δ(x −y)δ(t −τ).

So the Fourier transform with respect to x gives

G(k, t, y, τ) = e

−Dk

iky

δ(t − τ ) du,

where e

iky

is just the Fourier transform of δ(x − y). This is equal to

G(k, t; y, τ) =

(

0 t < τ

−iky

−Dk

(t−τ)

t > τ

= Θ(t − τ )e

−iky

−Dk

(t−τ)

Reverting the Fourier transform, we get

G(x, t; y, τ ) =

Θ(t − τ )

2π

ik(x−y)

−Dk

(t−τ)

This integral is just the inverse Fourier transform of the Gaussian with a phase

shift. So we end up with

G(x, t; y, τ ) =

Θ(t − τ )

4πD(t −τ)

exp



−

(x − y)

4D(t − τ)



= Θ(t − τ )S

(x − y; t − τ).

The solution we seek is then

φ(x, t) =

F (y, τ )G(x, t; y, τ ) dy dτ.

It is interesting that the solution to the forced equation involves the same function

(

x, t

) as the homogeneous equation with inhomogeneous boundary conditions.

In general, S

(x, t) solves

∂S

∂t

− D∇

= 0

with boundary conditions S

(x, 0) = δ

(n)

(x − y), and we can find

(x, t) =

(4πDt)

n/2

exp



−

|x − y|

4Dt



Then in general, given an initial condition φ|

t=0

= f(x), the solution is

φ(x, t) =

f(y)S(x − y, t) d

Similarly, G

(x, t; y, t) solves

∂G

∂t

− D∇

= δ(t − τ )δ

(n)

(x − y),

with the boundary condition G

(x, 0; y, τ ) = 0. The solution is

G(x, t; y, τ ) = Θ(t − τ)S

(x − y, t −τ).

Given our Green’s function, the general solution to the forced heat equation

∂φ

∂t

− D∇

φ = F (x, t), φ(x, 0) = 0

is just

φ(x, t) =

∞

F (y, τ )G(x, t; y, τ ) d

y dτ

F (y, τ )G(x, t; y, τ ) d

y dτ.

Duhamel noticed that we can write this as

φ(x, t) =

(x, t; τ) dτ,

where

F (y, t)S

(x − y, t −τ) d

solves the homogeneous heat equation with φ

t=τ

= F (x, τ).

Hence in general, we can think of the forcing term as providing a whole

sequence of “initial conditions” for all

t >

0. We then integrate over the times

at which these conditions were imposed to find the full solution to the forced

problem. This interpretation is called Duhamel’s principle.

Green’s functions for the wave equation

Suppose φ : R

× [0, ∞) → C solves the inhomogeneous wave equation

∂

∂t

− c

∇

φ = F (x, t)

with

φ(x, 0) =

∂

∂t

φ(x, 0) = 0.

We look for a Green’s function G

(x, t; y, τ ) that solves

∂G

∂t

− c

∇

= δ(t − τ )δ

(n)

(x − y) (∗)

with the same initial conditions

(x, 0, y, τ ) =

∂

∂t

(x, 0, y, τ ) = 0.

Just as before, we take the Fourier transform of this equation with respect the

spacial variables x. We get

∂

∂t

+ c

|k|

= δ(t − τ )e

−ik·y

where

(k, t, y, τ ).

This is just an ordinary differential equation from the point of view of

, and

is of the same type of initial value problem that we studied earlier, and the

solution is

(k, t, y, τ ) = Θ(t − τ)e

−ik·y

sin |k|c(t − τ )

|k|c

To recover the Green’s function itself, we have to compute the inverse Fourier

transform, and find

(x, t; y, τ ) =

(2π)

ik·x

Θ(t − τ )e

−ik·y

sin |k|c(t − τ )

|k|c

Unlike the case of the heat equation, the form of the answer we get here does

depend on the number of spatial dimensions

. For definiteness, we look at the

case where

= 3, since our world (probably) has three dimensions. Then our

Green’s function is

G(x, t; y, τ ) =

Θ(t − τ )

(2π)

ik·(x−y)

sin |k|c(t − τ )

|k|

We use spherical polar coordinates with the

-axis in

-space aligned along the

direction of x − y. Hence k · (x − y) = kr cos θ, where r = |x − y| and k = |k|.

Note that nothing in our integral depends on

, so we can pull out a factor

of 2π, and get

G(x, t; y, τ ) =

Θ(t − τ )

(2π)

∞

ikr cos θ

sin kc(t − τ)

sin θ dθ dk

The next integral to do is the

integral, which is straightforward since it is an

exact differential. Setting α = c(t − τ), we get

Θ(t − τ )

(2π)

∞



ikr

− e

−ikr

ikr



sin kc(t − τ)

Θ(t − τ )

(2π)

icr



∞

ikr

sin kα dk −

∞

−ikr

sin kα dk



Θ(t − τ )

2πicr



2π

∞

−∞

ikr

sin kα dk



Θ(t − τ )

2πicr

−1

[sin kα],

Now recall F[δ(x − α)] = e

−ikα

. So

−1

[sin kα] = F

−1



ikα

− e

−ikα



[δ(x + α) −δ(x − α)].

Hence our Green’s function is

G(x, t; y, τ ) = −

Θ(t − τ )

4πc|x − y|



|x − y|+ c(t − τ)



− δ



|x − y|− c(t − τ)



Now we look at our delta functions. The step function is non-zero only if

t > τ

Hence

|x − y|

(

t − τ

) is always positive. So

(

|x − y|

(

t − τ

)) does not

contribute. On the other hand,

(

|x −y|−c

(

t −τ

)) is non-zero only if

t > τ

. So

Θ(

t − τ

) is always positive in this region. So we can write our Green’s function

G(x, t; y, τ ) =

4πc

|x − y|

δ(|x − y|− c(t − τ )).

As always, given our Green’s function, the general solution to the forced equation

∂

∂t

− c

∇

φ = F (x, t)

φ(x, t) =

∞

F (y, τ )

4πc|x − y|

δ(|x − y|− c(t − τ )) d

y dτ.

We can use the delta function to do one of the integrals. It is up to us which

integral we do, but we pick the time integral to do. Then we get

φ(x, t) =

4πc

F (y, t

ret

)

|x − y|

where

ret

= t −

|x − y|

This shows that the effect of the forcing term at some point

y ∈ R

affects

the solution

at some other point

not instantaneously, but only after time

|x − y|/c

has elapsed. This is just as we saw for characteristics. This, again,

tells us that information travels at speed c.

Also, we see the effect of the forcing term gets weaker and weaker as we

move further away from the source. This dispersion depends on the number of

dimensions of the space. As we spread out in a three-dimensional space, the

“energy” from the forcing term has to be distributed over a larger sphere, and

the effect diminishes. On the contrary, in one-dimensional space, there is no

spreading out to do, and we don’t have this reduction. In fact, in one dimensions,

we get

φ(x, t) =

F (y, τ )

Θ(c(t − τ ) − |x − y|)

dy dτ.

We see there is now no suppression factor at the bottom, as expected.

7.5 Poisson’s equation

Let φ : R

→ R satisfy the Poisson’s equation

∇

φ = −F,

where F (x) is a forcing term.

The fundamental solution to this equation is defined to be G

(x, y), where

∇

(x, y) = δ

(3)

(x − y).

By rotational symmetry, G

(x, y) = G

(|x − y|). Integrating over a ball

= {|x − y| ≤ r, x ∈ R

we have

1 =

∇

∂B

n · ∇G

sin θ dθ dφ

= 4πr

So we know

4πr

and hence

(x, y) = −

4π|x − y|

+ c.

We often set c = 0 such that

lim

|x|→∞

= 0.

Green’s identities

To make use of this fundamental solution in solving Poisson’s equation, we first

obtain some useful identities.

Suppose

φ, ψ

→ R

are both smooth everywhere in some region Ω

⊆ R

with boundary ∂Ω. Then

∂Ω

φn · ∇ψ dS =

Ω

∇ · (φ∇ψ) dV =

Ω

φ∇

ψ + (∇φ) ·(∇ψ) dV.

So we get

Proposition (Green’s first identity).

∂Ω

φn · ∇ψ dS =

Ω

φ∇

ψ + (∇φ) ·(∇ψ) dV.

Of course, this is just an easy consequence of the divergence theorem, but

when Green first came up with this, divergence theorem hasn’t existed yet.

Similarly, we obtain

∂Ω

ψn · ∇φ dS =

Ω

ψ∇

φ + (∇φ) ·(∇ψ) dV.

Subtracting these two equations gives

Proposition (Green’s second identity).

Ω

φ∇

ψ − ψ∇

φ dV =

∂Ω

φn · ∇ψ − ψn · ∇φ dS.

Why is this useful? On the left, we have things like

∇

and

∇

. These

are things we are given by Poisson’s or Laplace’s equation. On the right, we

have things on the boundary, and these are often the boundary conditions we

are given. So this can be rather helpful when we have to solve these equations.

Using the Green’s function

We wish to apply this result to the case

(

|x − y|

). However, recall that

when deriving Green’s identity, we assumed

and

are smooth everywhere in

our domain. However, our Green’s function is singular at

, but we want to

integrate over this region as well. So we need to do this carefully. Because of

the singularity, we want to take

Ω = B

− B

= {x ∈ R

: ε ≤ |x − y| ≤ R}.

In other words, we remove a small region of radius

centered on

from the

domain.

In this choice of Ω, it is completely safe to use Green’s identity, since our

Green’s function is certainly regular everywhere in this Ω. First note that since

∇

= 0 everywhere except at x = y, we get

Ω

φ∇

− G

∇

φ dV = −

Ω

∇

φ dV

Then Green’s second identity gives

−

Ω

∇

φ dV =

φ(n · ∇G

) − G

(n · φ) dS

φ(n · ∇G

) − G

(n · φ) dS

Note that on the inner boundary, we have n = −

r. Also, at S

, we have

= −

4πε



4πε

So the inner boundary terms are



φ(n · ∇G

) − G

(n · φ)



sin θ dθ dφ

= −

4πε

φ sin θ dθ dφ +

4πε

(n · ∇φ) sin θ dθ dφ

Now the final integral is bounded by the assumption that

is everywhere smooth.

So as we take the limit

ε →

0, the final term vanishes. In the first term, the

’s

cancel. So we are left with

= −

4π

φ dΩ

= −

→ −φ(y)

where

φ is the average value of φ on the sphere.

Now suppose ∇

φ = −F . Then this gives

Proposition (Green’s third identity).

φ(y) =

∂Ω

φ(n · ∇G

) − G

(n · ∇φ) dS −

Ω

(x, y)F (x) d

This expresses

at any point

in Ω in terms of the fundamental solution

, the forcing term

and boundary data. In particular, if the boundary values

of φ and n ·∇φ vanish as we take r → ∞, then we have

φ(y) = −

(x, y)F (x) d

So the fundamental solution is the Green’s function for Poisson’s equation on

However, there is a puzzle. Suppose

= 0. So

∇

= 0. Then Green’s

identity says

φ(y) =



− G

dφ



dS.

But we know there is a unique solution to Laplace’s equation on every bounded

domain once we specify the boundary value

φ|

∂Ω

, or a unique-up-to-constant

solution if we specify the boundary value of n ·∇φ|

∂Ω

However, to get

(

) using Green’s identity, we need to know

and and

n · ∇φ on the boundary. This is too much.

Green’s third identity is a valid relation obeyed by solutions to Poisson’s

equation, but it is not constructive. We cannot specify

and

n ·∇φ

freely. What

we would like is a formula of φ given, say, just the value of φ on the boundary.

Dirichlet Green’s function

To overcome this (in the Dirichlet case), we seek to modify G

via

→ G = G

+ H(x, y),

where

∇

= 0 everywhere in Ω,

is regular throughout Ω, and

∂Ω

= 0. In

other words, we find some

that does not affect the relations on

when acted

on by

∇

, but now our

will have boundary value 0. We will find this

later,

but given such an

, we replace

with

G − H

in Green’s third identity, and

see that all the H terms fall out, i.e. G also satisfies Green’s third identity. So

φ(y) =

∂Ω

[φn · ∇G −Gn · ∇φ] dS −

F G dV

∂Ω

φn · ∇G dS −

F G dV.

So as long as we find

, we can express the value of

(

) in terms of the values

on the boundary. Similarly, if we are given a Neumann condition, i.e. the

value of

n · ∇φ

on the boundary, we have to find an

that kills off

n · ∇G

the boundary, and get a similar result.

In general, finding a harmonic ∇

H = 0 with

∂Ω

4π|x − x



∂Ω

is a difficult problem. However, the method of images allows us to solve this in

some special cases with lots of symmetry.

Example. Suppose

Ω = {(x, y, z) ∈ R

: z ≥ 0}.

We wish to find a solution to

∇

−F

in Ω with

φ →

0 rapidly as

|x| → ∞

with boundary condition φ(x, y, 0) = g(x, y).

The fundamental solution

(x, x

) = −

4π

|x − x

obeys all the conditions we need except

z=0

= −

4π

[(x − x

)

+ (y − y

)

+ z

]

1/2

6= 0.

However, let

be the point (

, y

, −z

) . This is the reflection of

in the

boundary plane

= 0. Since the point

is outside our domain,

(

x, x

)

obeys

∇

(x, x

) = 0

for all x ∈ Ω, and also

(x, x

z=0

= G

(x, x

z=0

Hence we take

G(x, x

) = G

(x, x

) − G

(x, x

The outward pointing normal to Ω at z = 0 is n = −

z. Hence we have

n · ∇G|

z=0

4π



−(z − z

)

|x − x

−

−(z + z

)

|x − x



z=0

2π

[(x − x

)

+ (y − y

)

+ z

]

3/2

Therefore our solution is

φ(x

) =

4π

Ω



|x − x

−

|x − x



F (x) d

2π

g(x, y)

[(x − x

)

+ (y − y

)

+ z

]

3/2

dx dy.

What have we actually done here? The Green’s function

(

x, x

) in some sense

represents a “charge” at

. We can imagine that the term

(

x, x

) represents

the contribution to our solution from a point charge of opposite sign located at

. Then by symmetry, G is zero at z = 0.

z = 0

Our solution for

(

) is valid if we extend it to

∈ R

where

(

) is solved

by choosing two mirror charge distributions. But this is irrelevant for our present

purposes. We are only concerned with finding a solution in the region

z ≥

This is just an artificial device we have in order to solve the problem.

Example.

Suppose a chimney produces smoke such that the density

(

x, t

) of

smoke obeys

∂

φ − D∇

φ = F (x, t).

The left side is just the heat equation, modelling the diffusion of smoke, while

the right forcing term describes the production of smoke by the chimney.

If this were a problem for x ∈ R

, then the solution is

φ(x, t) =

F (y, τ )S

(x − y, t −τ) d

y dτ,

where

(x − y, t −τ) =

[4πD(t −τ)]

3/2

exp



−

|x − y|

4D(t − τ)



This is true only if the smoke can diffuse in all of

. However, this is not true

for our current circumstances, since smoke does not diffuse into the ground.

To account for this, we should find a Green’s function that obeys

n · ∇G|

z=0

= 0.

This says that there is no smoke diffusing in to the ground.

This is achieved by picking

G(x, t; y, τ ) = Θ(t − τ)[S

(x − y, t −τ) + S

(x − y

, t − τ )].

We can directly check that this obeys

∂

∇

G = δ(t − τ)δ

(x − y)

when x ∈ Ω, and also

n · ∇G|z

= 0.

Hence the smoke density is given by

φ(x, t) =

Ω

F (y, τ )[S

(x − y, t −τ) + S

(x − y

, t − τ )] d

We can think of the second term as the contribution from a “mirror chimney”.

Without a mirror chimney, we will have smoke flowing into the ground. With a

mirror chimney, we will have equal amounts of mirror smoke flowing up from

the ground, so there is no net flow. Of course, there are no mirror chimneys in

reality. These are just artifacts we use to find the solution we want.

Example.

Suppose we want to solve the wave equation in the region (

x, t

) such

that x > 0 with boundary conditions

φ(x, 0) = b(x), ∂

(x, 0) = 0, ∂

φ(0, t) = 0.

On R

1,1

d’Alembert’s solution gives

φ(x, t) =

[b(x − ct) + b(x + ct)]

This is not what we want, since eventually we will have a wave moving past the

x = 0 line.

x = 0

To compensate for this, we introduce a mirror wave moving in the opposite

direction, such that when as they pass through each other at

= 0, there is no

net flow across the boundary.

x = 0

More precisely, we include a mirror initial condition

(

0) =

(

) +

(

−x

where we set

(

) = 0 when

x <

0. In the region

x >

0 we are interested in, only

the

(

) term will contribute. In the

x <

0 region, only

x >

0 will contribute.

Then the general solution is

φ(x, t) =

[b(x + ct) + b(x + c) + b(−x − ct) + b(−x + ct)].