Financial Stochastic Partial Differential Equations : Wiener Processes and Ito's Lemma

Suppose that a random variable satisfies

,

.

PDE methods for
Pricing Derivative Securities

Diane Wilcox

University of Cape Town
.
2005
Contents

I Mathematical Theory: Parabolic PDE 4
1 Introduction to Partial Differential Equations and Diffusion processes 5
1.1 Some general definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Classification of 2nd order linear PDE . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Function spaces and operators . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Linear DE and the principle of superposition . . . . . . . . . . . . . . . . . . 10
1.6 Types of problems within the solution of a PDE . . . . . . . . . . . . . . . . 11
1.7 Random walk derivation of a diffusion equation . . . . . . . . . . . . . . . . . 12
1.8 Heat flux derivation of a diffusion equation . . . . . . . . . . . . . . . . . . . 14
1.9 On mathematical modelling and Solving PDE . . . . . . . . . . . . . . . . . . 16

2 Fourier series methods 19
2.1 Fourier's solution to the heat equation - method of separation of variables . . 19
2.2 Trigonometric Series and pointwise convergence . . . . . . . . . . . . . . . . . 21
2.3 Bases for infinite-dimensional vectors spaces . . . . . . . . . . . . . . . . . . . 29
2.4 Sturm-Liouville problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.5 Applications to Finance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3 Integral transform methods 36
3.1 A heuristic motivation for fourier integral transform . . . . . . . . . . . . . . 36
3.1.1 Some basic properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.2 Convolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 The Convolution Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 Applying Fourier transforms to solve the IVP . . . . . . . . . . . . . . . . . . 41
3.4 Higher dimensional transforms . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4.1 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.5 Fourier sine and cosine transforms . . . . . . . . . . . . . . . . . . . . . . . . 43
3.6 Comments on Integral transforms . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.7 Laplace transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.8 Laplace Transforms and Convolution . . . . . . . . . . . . . . . . . . . . . . . 48

4 Fundamental solutions and Greens functions 50
4.1 Scaling and Similarity solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Introduction to Greens function . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3 The Dirac Delta function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.4 Rapidly Decreasing Functions and Distributions . . . . . . . . . . . . . . . . . 56
4.5 Distributional Evaluation of Functions and Greens Functions . . . . . . . . . 58
4.6 Application to Diffusion problems . . . . . . . . . . . . . . . . . . . . . . . . . 64




1
II Application to Derivative Securities 67
5 Black-Scholes PDE and solutions 68
5.1 The BSM differential equation . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1.1 Derivation of the PDE . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.2 Truncated expectation of log-normally distributed random variable . . . . . . 70
5.3 Solving the Black-Scholes PDE . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.4 The price of a European vanilla call . . . . . . . . . . . . . . . . . . . . . . . 72

6 First extensions of the Black-Scholes model 75
6.1 The model for dividends and applications . . . . . . . . . . . . . . . . . . . . 75
6.2 Greeks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3 Compound Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3.1 The value of a European call on European call . . . . . . . . . . . . . 75

7 Valuing American options 77
7.1 Bounds and put-call parity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.1.1 Upper bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.1.2 Lower bounds on European option values . . . . . . . . . . . . . . . . 77
7.1.3 Lower bounds and early exercise for American options on non-dividend
paying stock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.1.4 Put-call symmetry for American options . . . . . . . . . . . . . . . . . 80
7.2 American options on discrete dividend paying stock . . . . . . . . . . . . . . 81
7.3 The perpetual American put . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.4 The optimal exercise boundary & constraints for VA . . . . . . . . . . . . . . 84
7.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.4.2 Continuity of delta and the smooth-pasting condition . . . . . . . . . 85
7.4.3 The Black-Scholes inequality . . . . . . . . . . . . . . . . . . . . . . . 86
7.5 Integral solutions for American options on stock paying continuous dividends 87
7.5.1 Further properties of the optimal exercise boundary . . . . . . . . . . 87
7.5.2 Non-homogeneous PDE for American options . . . . . . . . . . . . . . 88
7.5.3 Integral solution for the option value . . . . . . . . . . . . . . . . . . . 90
7.6 Linear Complimentarity and Variational Inequality formulations . . . . . . . 91
7.6.1 An obstacle problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.6.2 Linear Complementarity problems and formulation for an American put 92
7.6.3 Variational inequality formulation for an American put . . . . . . . . 93

A Glossary on Linear Algebra terms 96

B On flow, flux, divergence, curl and Green's theorem in 2D 100
B.1 Basics from vector calculus (calculus of several variables) . . . . . . . . . . . 100
B.2 Green's Theorem Theorem (in R2 ) . . . . . . . . . . . . . . . . . . . . . . . . 105

C Taylor series, derivative matrices, approximations and Jacobians 106
C.1 Taylor series for several variables . . . . . . . . . . . . . . . . . . . . . . . . . 106
C.2 Affine approximation, the derivative matrix . . . . . . . . . . . . . . . . . . . 107
C.3 The Jacobian I - geometry of a transformation . . . . . . . . . . . . . . . . . 108
C.4 The Jacobian II - the scaling factor in integrals . . . . . . . . . . . . . . . . . 109




2
D Quick review for solving 1st and 2nd order linear ODE 111
D.1 First and simplest methods for solving simple 1st order linear ODE . . . . . . 111
D.2 On existence and uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
D.3 Linear Independence, the Wronskian, fundamental and general solutions . . . 114
D.4 Differential operator methods for linear ODE . . . . . . . . . . . . . . . . . . 116
D.5 Solving 2nd order linear ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . 117




3
Part I

Mathematical Theory:
Parabolic PDE




4
Chapter 1

Introduction to Partial
Differential Equations and
Diffusion processes

Many phenomena are described by functions whose values at a given point (for example, in
time) depends on values at neighbouring points. The equation which determines such a
function usually contains derivatives of the function in order to capture information about
rates of change of the function with respect to the underlying variable(s) or parameter(s).
Such equations are referred to as differential equations .


1.1 Some general definitions
Formally, an ordinary differential equation (ODE) is an equation involving an unknown
function u of one independent variable x and derivatives of that function. By appropriate
rearrangement of terms, such an equation may be written as :

du d2 u d3 u
F (x, u, , , , . . .) = 0. (1.1)
dx dx2 dx3
The order of the equation is the order of the highest derivative and the equation is said to
be linear if u and all it's derivatives are of the 1st degree (there are no terms of the form
u3 or ( du )2 , etc ... ) and there are no products of u and it's derivatives.
dx

A general nth order linear ODE may be written as :
du dn u
a0 (x)u(x) + a1 (x) + . . . + an (x) n = f (x), (1.2)
dx dx
where ai (x), i = 1, . . . n are known functions of x. If f (x) = 0, then equation ( 1.2) is said
to be homogeneous.

A partial differential equation (PDE) is an equation involving an unknown function
u of two or more one independent variables and partial derivatives of that function. By
appropriate rearrangement of terms, such an equation in two independent variables x and
y may be written as

u u 2 u 2 u 2 u ku
F (x, u, , , 2, , 2 , . . . , r s , . . .) = 0, (1.3)
x y x xy y x y


5
where k = r + s.

As before, the order of the equation is the order of the highest derivative and the equation
is said to be linear if u and all it's derivatives are of the 1st degree (i.e. there are no terms
2
of the form u2 or ( u )3 , etc... ) and there are no products of u and it's derivatives.
x2

A 1st order linear PDE may be written as :
u u
a(x, y) + b(x, y) + c(x, y)u = g(x, y). (1.4)
x y

The general form of a 2nd order linear PDE is given by

auxx + 2buxy + cuyy + dux + euy + f u = g, (1.5)
u 2u
where a, b, c, d, e and f are known functions of x and y and ux := x , uxx := x2 and
2u
uxy := xy , etc ...

If g = 0 in equation ( 1.4) or equation ( 1.5) then the equation is said to be homogeneous.


1.2 Examples

The best known PDE in the mathematics of financial markets is, of course, the Black-Scholes
equation for valuing options:


V 1 2V V
+ 2 S 2 2 + rS - rV = 0,
t 2 S S

where S denotes stock price, t denotes the time, and V , the value of the option, is a function
of S and t. The equation itself is quite a general object and one cannot speak of a solution
without reference to additional constraints, such as possible values which the option may
attain when exercised. These conditions specify the valuation problem at hand.

Historically, three 2nd-order linear PDE have been of fundamental interest because they
exhibit distinct behaviour and have led to the clarification of general theories and methods:

Let u(x) : R2 R be a function of two variables.

· The Laplace equation originated out of problems in potential theory in physics :
2u 2u
+ 2 = 0.
x2 y


· The Wave Equation describes the amplitudes u(x, t) of vibrating membranes and
other wave functions :


2u 2u
= 2.
x2 t



6
· The Heat Equation describes the temperature distribution u(x, t) in the plane when
heat is allowed to flow from warm areas to cool ones :


2u u
2
= .
x t
The above examples are models of problems with 2 independent variables. More generally
we may let u(x) : Rn R be a function of x = (x1 , x2 , . . . , xn ) Rn and the above
examples may be formulated naturally in higher dimensions. Since methods of solution often
generalise to higher dimensions, it is often sufficient to understand examples for n = 1, 2 or
3. However, it is sometimes the case that a particular method for solving is not suited to
higher dimensions.


1.3 Classification of 2nd order linear PDE
The behavior of known solutions of the classic 2nd order PDE's indicates that the presence
of the higher order terms is significant. Thus, for a general theory of 2nd order PDE's, we
consider an equivalent form of equation ( 1.5):
2u 2u 2u
a(x, y) 2
+ 2b(x, y) + c(x, y) 2 = F (x, y, u, ux , uy ). (1.6)
x xy y
The discriminant1 of a 2nd order linear PDE is defined as :
D := ac - b2 . (1.9)
The PDE is said to be


elliptic if D > 0,
parabolic if D = 0,
hyperbolic if D < 0.


N.B.: If a, b or c are function of the independent variables, then the discriminant varies
with the values of these variables.

By this classification,

the heat equation is parabolic everywhere (a = 1, b = 0 and c = 0),

the Laplace equation is elliptic everywhere and

the wave equation is hyperbolic everywhere.

1 Alternatively, if we consider the general form ( 1.5):
2u 2u 2u
a(x, y) + b(x, y) + c(x, y) 2 = F (x, y, u, ux , uy ).. (1.7)
x2 xy y
then the discriminant is equivalently defined as :
D := 4ac - b2 . (1.8)
b2
Sometimes, it defined D := - 4ac. In this case, the inequalities for the classification of PDE as elliptic or
hyperbolic must be reversed.


7
Examples 1.3.1. For the PDE, 2xuxx - utt = 0, the discriminant is D = -2x. Thus, the
equation is

elliptic if x < 0,
parabolic if x = 0,
hyperbolic if x > 0.

We will see that the Black-Scholes equation is also a parabolic equation. In fact, by a
suitable change in variables it is possible to transform the BS equation into one which has
the form of the heat equation.

The flow of heat through a solid, and the continuous Brownian motion of particles in a
liquid are examples of diffusion processes. These processes can be modelled by the heat
equation. The forward Kolmogorov equation, also known as the Fokker-Planck equation, is
another parabolic equation derived in the modelling diffusion processes :
u u 1 2u
= -(c ) + D 2 ,
t x 2 x
where c and D are constants.


1.4 Function spaces and operators
A Vector space V is a non-empty set of vectors x, y, z, . . . endowed two operations referred
to as addition and scalar multiplication (i.e. multiplication by a scalar quantity) such that
the following axioms are satisfied:

Let K denote R or C (scalar field). For any x, y and z V and , K we have

Property Interpretation
V1 x + y is a vector in V V is closed under addition
V2 x+y =y+x Addition is commutative
V3 (x + y) + z = z + (y + z) Addition is associative
V4 There exists a unique vector 0 V V has an additive identity
such that x + 0 = 0 + x = x
V5 For every x V there exists a unique vector For every x V there exists an
-x V such that x + (-x) = 0 = (-x) + x additive inverse
V6 x V V is closed under scalar multiplication
V7 (x + y) = x + y Distributivity I
V8 ( + )x = x + x Distributivity II
V9 (x) = ()x Multiplication is associative
V10 If = 1, then x = x Scaling by unity


Remarks : Because of this linear structure, vector spaces are often referred to as linear
spaces. To establish that a set is a vector space, one need only verify that it is closed under
addition and scalar multiplication.


8
Examples 1.4.1.

1. The spaces Rn or Cn of n-tuples of real and complex numbers, respectively, are vector
spaces under usual addition and scalar multiplication.
2. The set Mn (R) of n × n matrices with real entries, is a vector space under matrix
addition.
3. The set Pn of polynomials of degree at most n is a vector space under the usual addition
and scalar multiplication of functions2 .
4. The set C([a, b]) := {f : [a, b] R | f is continuous on [a, b] } is a vector space under
the usual addition and scalar multiplication of functions.

5. If Rn then the set

C() := { u : R | u is continuous on }

of real-valued functions of n variables is a vector space. Similarly, for m 1,

C m () := {u : R | u and its partial derivatives to order m are continuous on }

of real-valued functions of n variables is a vector space.


Definition 1.4.2. Given vector spaces V and W , an operator T : V W is a mapping
from V into W which assigns a vector w = T v W for each vector v in the domain D(T )
of T .
Examples 1.4.3.

1. An r × s matrix is an operator from Rr into Rs .
df
2. Let R. Then D : C 1 () C() given by Df := dx is a differential operator.
3. The operator L : C n () C() given by
n
u
Lu := , 1 =
x1
1 . . . xn
n
i=1

is a partial differential operator.

Definition 1.4.4. A linear operator is an operator L : V W which satisfies
L(x + y) = L(x) + L(y)
for all x, y V and scalars and .


Examples 1.4.5. The operators in Examples 1.4.3 are all linear operators.
2 (f + g)(x) = f (x) + g(x) and (f )(x) = (f (x))




9
The following linear operators are basic to the study of PDE :

2
2
1) Laplace operator L = := x2 + . . . + x2
1 n

2) Diffusion operator L= t -
2
3) D'Alembert operator L= t2 -


The above operators are referred to as elliptic, parabolic and hyperbolic operators, respec-
tively, by the classification of the associated equation Lu = f . By a suitable change of
variables, an elliptic operator L may be transformed into the Laplacian in an alternative
co-ordinate system. Similarly, parabolic and hyperbolic operators may be transformed into
the diffusion and d'Alembert operators, respectively.


1.5 Linear DE and the principle of superposition

Recall that the operation of differentiation acts as a linear transformation on the set of
continuous real-valued functions which have 1st derivatives. In other words, if f and g are
continuous real-valued functions which have 1st derivatives, and and are scalars ( i.e.
, R are constants ) then
d d d
(f + g) = f + g.
dx dx dx

The same is true for partial differentiation. Thus, linear partial differential equations may be
written in the form Lu = f , where u is unknown and f is given function, and in general, a lin-
ear partial differential operator, L, transforms a function u of
x = (x1 , x2 , . . . xn ) Rn , into another function L(u) given by:

u u u
L(u) = a(x)u + b1 (x) + b2 (x) + . . . + bn (x)
x1 x2 xn
2u 2u 2u
+ c11 (x) 2 + c12 (x) + . . . + cnn (x) 2 + . . .
x1 x1 x2 xn
n n
u 2u
= a(x)u + bi (x) + cij (x) + ...,
i=1
xi i,j=1
xi xj

where the dots indicate higher-order terms, but it is understood that the sum contains only
finitely many terms. Clearly we may write L independently of the function u as follows:
n n
2
L := a(x) + bi (x) + cij (x) + ...
i=1
xi i,j=1
xi xj

If u1 , . . . uk are functions which are sufficiently smooth, i.e. u1 , . . . uk have enough deriva-
tives for L(ui ) to be defined for 1 i k, and if c1 , . . . ck are scalars (constants )
then
L(c1 u1 + . . . + ck uk ) = c1 L(u1 ) + . . . + ck L(uk ).
For infinite linear combinations such that the series

ck uk and ck L(uk )
k=1 k=1


10
both converge, we have

L[ ck uk ] = ck L(uk ).
k=1 k=1

For integrals, we generalise as follows: suppose u(x, ) is a function of x Rn and of the
parameter , where a < < b, and suppose g() is an integrable function of on (a, b),
then
b b

L[ g()u(x)d ] = g() L [u(x)] d.
a a

b
Here we may think of g()u(x)d as the continuous form of the series ck uk in
a k=1
which the coefficients ck are replaced with g() in the "continuous" sum.

The principle of superposition states :
If ui satisfy the a linear homogenous equation, then an arbitrary linear combi-
nation c1 u1 + . . . + ck uk satisfies the same homogeneous equation.
Symbolically, we may write this as follows :

L(uk ) = 0 f or all k L[ ck uk ] = ck L(uk ) = 0.
k=1 k=1

For the integral (continuous sum) version , the principle of superposition states:
b b

L[u(x, )] = 0 f or a < < b L[ g()u(x, )d ] = g()L[u(x, )]d = 0.
a a


1.6 Types of problems within the solution of a PDE
(1) Boundary Value Problems (BVP)

e.g. the Dirichelet problem on a domain Rn , where is open and connected3 , refers to
the problem of finding a function u which satisfies the Laplace equation within , i.e.
u(x) = 0 for x
and which obeys the additional condition that
u(x) = f for x .
where denotes the boundary of and f is a known function defined on .

(2) Initial (and Final) Value Problems (IVP)

e.g. the heat equation specified as:


ut - uxx = 0 f or - < x < , t > 0
u(x, 0) = f (x) f or - < x < , t = 0,
3 for now we may take these conditions to simply mean that is a nice and unbroken subset of Rn


11
where f denotes some initial temperature distribution and the value of u(x, t) evolves with
time.

e.g. the Black-Scholes PDE for which the payoff function is given. Here the payoff function
is a final condition and depends on ST , the price of the stock at some expiry (final) time T .


V 1 2V V
+ 2 S 2 2 + rS - rV = 0 f or 0 t < T
t 2 S S
V (ST , T ) = payoff at t = T,

where S denotes stock price, t denotes the time, T is the expiry date and V , the value of
the option, is a function of S and t. The BS-PDE can be transformed in to an initial value
problem by letting := T - t. The final condition V (ST , T ) is now an initial condition in
the variable , i.e. V (S , 0) = payoff at = 0. The PDE is referred to as a backward PDE:

(3) Eigenvalue Value Problems (EVP)

Definition 1.6.1. Let L : V W be a linear operator. A non-zero vector v V is said
to be an eigenvector for L if there exists a scalar such that Lv = v. If V is a function
space, then an eigenvector is often referred to as an eigenfunction.
Consider the following ODE4

v = v, 0 < x <
v(0) = v() = = 0,
dv
where v : R R and v := dx . Clearly v = 0 is a solution. We are interested in non-zero
solutions, i.e. we want to find the eigenfunctions vn (x) and their corresponding eigenvalues
d2
for the problem Lv = v, where L = dx2 .


1.7 Random walk derivation of a diffusion equation
Let (, P) denote a probability space. A discrete valued random variable (r.v.), x, defined
on all of assumes values xi with probability pi ( pi = 1). Its expectation and variance
i
are

E(x) := xi pi , denoted x
i
2
var(x) := x- x = x2 - x 2 .

Consider an unrestricted 1-D random walk: A particle at the origin executes random
steps + (right) or - (left). Let xi denote the r.v. which assumes the values + or - at
the ith step; assume that each step is independent of the others. By construction, the xi are
independent, identically distributed (i.i.d.) r.v. The position of the particle after the nth
step is given by
n
Xn = xi .
i=1
4 such problems arise in the method of separation of variables and Sturm-Liouville problems which we

will meet in subsequent chapters.


12
Let

P(xi = +) = p
P(xi = -) = q,

where p + q = 1. It follows that

n n
E(Xn ) = E( xi ) = E(xi ) = (p - q)n (1.10)
i=1 i=1
E(x2 ) = (+) p + (-) q = 2
i
2 2

n n n
var(Xn ) = var( xi ) = var(xi ) = [ x2 - xi 2 ]
i
i=1 i=1 i=1
n
= [ 2 - (p - q)2 2 ] = 4p q 2 n. (1.11)
i=1


For the continuous case, consider the following mathematical model for physical Brownian
motion. Experimentally, the average displacement of the particle per unit time is c, and the
variance of the observed displacement is D > 0.

Suppose the particle executes r jumps per unit time. Then from ( 1.10) and ( 1.11) we have

(p - q)r c (1.12)
4p q 2 r D. (1.13)

The motion of particles appear continuous, so we examine behaviour under the limit 0
and r 0, constrained by the conditions that c and D remain fixed. If p = q and

lim (p - q) = 0 = lim (p - q),
0 r

then as 0, r jointly, we have
c
r - .
p-q
Hence, as 0, r jointly,

4p qc
4p q 2 r - - .
p-q

But, from ( 1.13),

4p q 2 r - D = 0.

Thus, lim(p - q) 0, with p + q = 1, so that lim p = lim q = 2 . Now p = q = 1 implies
1
2
that c = 0 and D = 0, which in turn implies that there is no variance and the process is
deterministic. Thus, ( 1.10) and ( 1.11) hold if we have
1
p = (1 + b)
2
1
q = (1 - b),
2


13
where b is a constant such that 0 p, q 1 and p + q = 1. In this case we have p - q = b.
It follows that as 0, r jointly, 2 r D, p 1 , q 2 and b = D .
2
1 c

1
Now there are r steps per unit time if and only if 1 step of length occurs in time := r
units.

Returning to the motion of our particle starting at x = 0, t = 0, we are interested in the
probability of the particle being at position x at time t after n steps, i.e. we wish to know
v(x, t) := P(Xn = x, n = t).
For the continuum limit, we construct difference equations for v(x, t) and show that it
converges to the a PDE. We have (the master equation) :

v(x, t + ) = p v(x - , t) + q v(x + , t). (1.14)
Expanding the Taylor series

v(x, t + ) = v(x, t) + vt (x, t) + O( 2 )
1
v(x ± , t) = v(x, t) ± vx (x, t) + 2 vxx (x, t) + O( 3 ), (1.15)
2
O(y k )
where O(y k ) means that lim yk
< 5 . Substituting ( 1.15) into ( 1.14), we obtain
y0

1 2 3
vt (x, t) = (q - p) vx (x, t) + ( )vxx (x, t) + O( ) + O( ). (1.16)
2
As 0 and 0, equation ( 1.16) converges to
1
vt = -cvx + Dvxx ,
2
and v(x, t) may interpreted as the probability density function for the continuous r.v. x at
time t.


1.8 Heat flux derivation of a diffusion equation
Suppose the problem is to model heat flowing through a medium, e.g. a gas or a liquid,
located in 3-D space.




Let x = (x1 , x2 , x3 ) denote the co-ordinates of a point in space. We want to find the
temperature distribution in the medium, u(x, t). The value of u(x, t) will depend on
5 i.e. the error due to ignoring higher order terms is smaller in absolute value than some constant times

y k if y is close enough to 0.


14
· the presence of heat in the body and
· the flow of heat into or out of the body through the surface.

There are two laws ( = equations ) which model the situation :

1. The balance or conservation equation for the conservation of energy.
2. The constitutive equation which describes how heat flows in the medium.


Assuming there are no other types of energy in the situation, (1) can be stated as follows:

rate of change heat generated by sources within the body
of = +
thermal energy in body flow of heat into body from outside.


Using the notation :
x= (x1 , x2 , x3 ) = position in space
t= time
= small volume element (= x1 x2 x3 )
R= region occupied by medium
S= surface of medium
c= heat capacity ( = amount of heat generated per unit mass and per unit rise
in temperature )
= mass density
u= temperature
f (x, t) = energy generated inside the body per unit volume and per unit time
q(x, t) = (normal component of) heat flux through surface (which may be negative)

and assuming that properties of the medium are identical throughout the medium ( the
medium is said to be homogeneous in this case), we have :

total thermal energy in region = c(x)(x)u(x)dV
R
( c.(. ).u ) ,


total energy from source = f (x, t)dV,
R

and
total flux through surface = - q(x, t).ndA
S
= - div q.dV,
R


where n is the outward normal vector to S and the last equality follows from the divergence
theorem. Thus, from the balance equation, we have :

d
cpu dV = f dV - div q dV (1.17)
dt
R R R



15
or
du
(cp - f + div q)dV = 0.
dt
R

Thus, since the volume is arbitary and the integrand is smooth,

du
cp - f + div q = 0.
dt
The constitutive equation is given by Fourier's law :


heat flux = q = -k u,


where k is referred to as the thermal conductivity function and depends on the material of
the medium6 . Thus ( 1.17) becomes

u f
- div(k u) = Q, Q= , (1.18)
t cp
which is referred to as the heat equation.

To solve this equation we need information about u on S and u at time t = 0.

The following properties may be deduced from this equation :

1. heat energy is conserved
2. heat flows in a manner that distributes evenly throughout the medium (from hot to
cold and faster when there is greater initial temperature differences).



1.9 On mathematical modelling and Solving PDE
Very broadly, the process of solving a problem by mathematical methods involves three
components, namely:

(1) formulating the problem mathematically,
(2) solving the mathematical problem, and
(3) interpretting the solution.

Once a model has been formulated, one needs to know :
6 Here


· k < 0 implies heat flows from hot to cold

· u = 0 implies temperature is constant and there is no heat flow

· u = large implies q is large, i.e. the greater the temp. difference, the greater the flow.




16
· does a solution exists?

This is sort of question is referred to as the existence problem. For example, there are some
mathematical formulations of minimisation problems which do not have solutions. In this
case, the problem probably does not lie with the math, but with the formulation. 7 We will
be fortunate to work with problems for which solutions are known to exist. The simplest
way to prove existence is to construct a solution that meets all the conditions originally
imposed. One also is interested:
· is the solution unique ?

8
· is the solution stable ?

These questions lead to the following definition:

Definition 1.9.1. A solution to a PDE is said to depend continuously on data if small
changes in the data (boundary or initial conditions) produce small changes in the solution.
A problem is said to be well-posed if

(1) a solution to the problem exists,
(2) the solution is unique, and
(3) the solution depends continuously on data


In general, the number of initial / boundary conditions which are required to solve a PDE
uniquely, depends on several factors - in general, it is complicated to determine the appro-
priate form for data needed for a solution. The notion of well-posedness was proposed by
Hadamard to serve as set of a guidelines for formulating a problem and determine which
data are necessary for a solution to be found 9 . The last condition is important for many
applications - it is preferred that the solution changes only a little when conditions which
specify the problem change a little. More can be said on the well-posedness of a problem
with the aid of deeper analysis of geometry of the classification of 2nd order PDE's - we will
not investigate that here.

Another question is :
· what method was used to construct the solution ?

This is particularly relevant when numerical methods are used to solve the problem: one
wants to know how good an approximation is and is concerned about the error in the
solution, how long an iteration takes to converge, etc... Rigour is always important when
there may be exceptional cases to a seemingly logical argument. Later we will see that when
using a series expansion to represent a solution, one must be sure that the series converges!

7 Problems of proving the existence of solutions have been important in the development of mathematics.
The Dirichelet problem, for example, had a significant early role in the development of functional analysis.
8 Here we are asking whether the solution changes significantly if a part of the input information / data

is changed only very slightly.
9 When the problem is set in a homogeneous medium, when the shape of the boundary is "nice" and when

the boundary data is simple, then solutions to PDE exist. In particular, there are always solutions for PDE
in this course!


17
And another is:
· how differentiable is the solution found ?

This refers to the problem if regularity of the solution. Depending on the model, it may
be required that the solution be very regular, say k-times continuously differentiable. Such
problems may be really hard to solve and proofs would have to include verifications that solu-
tions are smooth enough. An alternative strategy is to consider the existence and smoothness
problems separately. The idea is to consider the PDE for a very wide class of functions and
to find a weak solution. Since not much is being expected of such a solution, it is anticipated
that questions on existence, uniqueness and stability are easier to solve. Once it is estab-
lised that the problem is well-posed for a broad class of generalised solutions, the question of
regularity of solutions may be addressed - for some problems, it may be that weak solutions
are smooth enough to qualify as a meaningful solution10 .

In applications, we are not interested in general solutions - we usually seek some specific
solution to the given problem, where the unique solution is specified by additional data, i.e.
boundary conditions. In the case where one of the variables is a time parameter, we may
have some initial or final condition for the problem.

In the case of the diffusion equation, which contains a time dependence, we have will find it
appropriate to specify the function u at initial time t = 0, i.e. the initial value of the density
distribution is given as part of the problem . For the wave equation, which also contains
time dependency, both u(x, t) and ut (x, t) are usually specified, since the problem usually
involves a description of the vibrations of the system as a result of some initial impulse.

When the values of x are constrained to a bounded or partially bounded region, then u(x, t)
and/or ux (x, t) (or a linear combination of these) must be given on the boundaries for
all t > 0. Even when the region is not constrained, there are usually explicit or implicit
conditions for the behaviour of the solution at infinity.

Boundary value problems for parabolic and hyperbolic equations may not well-posed. So-
lutions to these equations evolve in time, and their behaviour at later times is given by
previous states. Thus, a boundary value problem which arbitrarily specifies the values of
the solution at two or more distinct times, is usually not reasonable.

Examples 1.9.2.
Consider the hyperbolic equation uxy = 0 on the square [0, 1] × [0, 1]. Since uxy = 0, it follows
that ux is constant as a function of the parameter y, and in particular, ux (x, 0) = ux (x, 1). Now
if ux (x, 0) = f (x) and ux (x, 1) = g (x) are given by the boundary conditions u(x, 0) = f (x) and
u(x, 1) = g(x), then the problem cannot be solved for arbitrary f and g; a solution exists only if
f (x) = g (x).




10 We will see more of this notion when we consider Green's functions and solutions to PDE


18
Chapter 2

Fourier series methods

2.1 Fourier's solution to the heat equation - method of
separation of variables
The problem which Joseph Fourier (1768 - 1830) considered in his paper On the Propagation
of Heat in Solid Bodies1 : was
find a solution for the diffusion of heat in a 1-dimensional rod of length l
given by the equation with boundary conditions and initial condition:

u 2u
- k 2 = 0, (2.1)
t x
u(0, t) = u(l, t) = 0 f or t > 0, , (2.2)
u(x, 0) = f (x) f or 0 < x < l. (2.3)

His solution can be summarised as follows:

Let u(x, t) = (x)(t). Then, substituting into equation ( 2.1), we get
xx (x) kt (t)
= . (2.4)
(x) (t)
Since equation ( 2.4) holds for all x such that 0 x l, and for all t > 0, the LHS and
RHS must be equal to some constant, say -. Thus, ( 2.4) is given by:


xx (x) + (x) = 0, and (2.5)
t (x) + k(t) = 0. (2.6)

Now the boundary conditions give at ( 2.2) imply that (0) = (l) = 0. Thus, the solution
to the eigenvalue problem ( 2.5) can be written:

nx
(x) = cn . sin( ), (2.7)
l
1 His work on the topic began in 1804 and was completed in 1807. Under review by Lagrange, Laplace,

Monge and Lacroix, the work was highly regarded but criticized, mainly for the controversial trigonometric
expansions for the representation functions. In this regard, his arguments were considered incomplete and
not general.


19
with corresponding eigenvalues
n 2
n = ( ) .
l
The solution to equation ( 2.6) is then given by

(t) = e-n kt . (2.8)

Thus, combining ( 2.7) and ( 2.8), it follows that we have solutions
nx -( n )2 kt
un (x, t) = cn sin( )e l , n N. (2.9)
l
Since ( 2.1) is a linear equation, it follows from ( 2.9) that our general solution can be written

u(x, t) = un (x, t)
n=1

nx -( n )2 kt
= cn sin( )e l (2.10)
n=1
l

with the initial condition that u(x, 0) = f (x). The problem is solved with :
provided

nx
f (x) = cn sin( ). (2.11)
n=1
l

This last statement leads to the question: when can a function be represented as a trigono-
metric series?

A rigorous answer for pointwise representation was first given by Dirichet in 1829. The
contemporary form of a much more powerful theorem on the uniform convergence of a series
representation for square-integrable functions was developed only after the theory of Hilbert
spaces was fleshed out.




20
2.2 Trigonometric Series and pointwise convergence
We continue with the question : which functions admit trigonometric representations of the
form ( 2.11)? Since these are sums of periodic functions, a natural place to begin is with
the class of periodic functions (since we asking about the convergence of a series made up
of periodic functions ).

Definition 2.2.1. A function f defined on the real line is said to be periodic with period
P if f (x) = f (x + P ) for all x R.

Recall that the sine function satisfies sin(x) = - sin(-x) and that the cosine function satisfies
cos(x) = cos(-x) and we have the following general definitions:

Definition 2.2.2. A function f defined on an interval [a, b] is said to be even if it satisfies
f (x) = f (-x) for all x [a, b], it is said to be odd if it satisfies f (x) = -f (-x) for all
x [a, b].
We note that the definition is vacuous when a 0 or when b 0.

Any function can be written as the sum of an odd function, f = f0 + f1 , where f0 is an the
function and f1 is even function defined as follows:
1 1
f0 (x) := ( f (x) - f (-x) ) and f1 (x) := ( f (x) + f (-x) ).
2 2

We could rephrase our question: which functions on the interval [-, ] admit trigonometric
representations of the form :


a0
f (x) = + [ an cos(nx) + bn sin(nx) ] ? (2.12)
2 n=1


The restriction to functions which are 2 periodic on the interval [-, ] is just
a mathematical convenience. The results which are developed hold equivalently
for functions of period P by making suitable adjustments (by rescaling shift-
ing of co-ordinate axes) for the problem considered. For an arbitrary function
f on a bounded interval [a, b], one may pass to its periodic extension2 - any
~ ~
mathematical conclusions for f are valid for the restriction of f to f defined on
[a, b].

Hence, we will proceed to tackle the problem: which functions on the interval [-l, l] admit
trigonometric representations of the form :


a0 nx nx
f (x) = + [ an cos( ) + bn sin( )] ? (2.13)
2 n=1
l l

Equivalently, we could develop the theory for functions on in the interval [0, l].
2 The periodic extension of a function f defined on [a, b] is given by

f (x)
for x [a, b],
~
f = f (x - k(b - a)) for x > b, x = y + k(b - a) for some y [a, b] and k N+

f (x + k(b - a)) for x < a, x = y - k(b - a) for some y [a, b] and k N+
~
Clearly f is (b - a)-periodic and defined on the entire R


21
Orthogonality of sines and cosines and calculating coefficients for
Fourier series

For m, n N, we have
l l
nx mx nx mx
sin( ) cos( )dx = sin( ) cos( )dx = 0,
l l l l
-l 0

and
l l
1 mx nx 2 mx nx 0 if m = n,
sin( ) sin( )dx = sin( ) sin( )dx =
l l l l l l 1 if m = n = 0,
-l 0

and finally
l l

1 mx nx 2 mx nx 0 if m = n,
cos( ) cos( )dx = cos( ) cos( )dx = 1 if m = n = 0,
l l l l l l
-l 0 2 if m = n = 0.


Now suppose

nx
f (x) = bn sin( ). (2.14)
n=1
l

1
Then multiplying the LHS and RHS of ( 2.14) by l sin( kx ) and integrating over [-l, l], we
l
have

l l
1 kx 1 nx kx
f (x) sin( )dx = bn sin( ) sin( )dx
l l l l l
-l -l n=1

l
bn nx kx
= sin( ) sin( )dx
n=1
l l l
-l
= ck

Exercises 2.2.3. Derive the formula for coefficients for a function defined on [0, l].




22
Definition 2.2.4. If f is 2l-periodic and Riemann integrable on [-l, l], then the coefficients
ak and bk are given by

l
1 kx
ak := f (x) cos( )dx
l l
-l
l
1 kx
bk := f (x) sin( )dx
l l
-l

are called the Fourier coefficients of f and the corresponding series,

a0 nx nx
+ [ an cos( ) + bn sin( )]
2 n=1
l l

is called the Fourier series of f .


NB : We have NOT shown or assumed that we have equality ( 2.13). We have
merely given a name to the series expansion appearing on the RHS of ( 2.13).


Exercises 2.2.5.

1. Consider Definition 2.2.4. Give the corresponding definition if f is l-periodic and
Riemann integrable on [0, l].
2. If f is periodic with period P then show that the following integral is independent of a:
a+P

f (x)dx
a


3. Show that

a 0
if f is odd

f (x)dx = a
2
f (x)dx if f is even
-a
0


4. Let f be the 2 periodic function given by f (x) = |x| for - x . Sketch the
function and find its Fourier series.




23
Complex Fourier series
We recall the properties of the complex exponential function and how it is related to the
cosine and sine functions of a variable x R:

eix= cos x + i sin x
1 ix
cos x = (e + e-ix )
2
1 ix
sin x = (e - e-ix )
2i
The advantages of working with the sine and cosine functions are that they are odd and
even functions, respectively. The advantage of the working with the exponential function is
that we have
d ix
(eix ) = (e ) = ieix , and
dx
ei(x1 +x2 ) = eix1 + eix2 .

Using the complex exponential function we can rewrite equation ( 2.13) by :

f (x) = cn einx/l (2.15)
-

where
a0
c0 = ,
2
1
cn = (an - ibn ), and
2
1
c-n = (an + ibn )
2
for n N. Equivalently we may start with equation ( 2.15) and derive equation ( 2.13).
The odd and even properties of sine and cosine are used when n is negative in equation
( 2.15)), and the coefficients are given by

a0 = 2c0 ,
an = cn - c-n , and
bn = i(cn - c-n ),

for n N.

Assuming that a function f defined on [-l, l] is given by the series representation of ( 2.15)
(or equivalently ( 2.13), then we may find any coefficients cm (or, equivalently, am and bm )
by multiplying the LHS and RHS of ( 2.15) by e-imx/l , integrating from -l to l, and then
solving for cm .
Exercises 2.2.6. Show that
l
1 0