LectureNotesMATH2048 PDF

MATH2048
Mathematics for Engineering

and the Environment Part II
David Gammack
Kostas Skenderis
2016-17
Contents
1 Ordinary Differential Equations 4

1.1 Second order equations . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Linear Second Order ODEs . . . . . . . . . . . . . . . . 4
1.1.2 Inhomogeneous Linear Second Order ODEs . . . . . . . . 8
1.2 BVPs and IVPs . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1 Eigenvalues and functions . . . . . . . . . . . . . . . . . 12
1.2.2 Sturm-Liouville form . . . . . . . . . . . . . . . . . . . . 14
2 Fourier Series 19
2.1 Fourier series . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Orthogonality . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Fourier Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 Even and odd functions . . . . . . . . . . . . . . . . . . . 22
2.3.2 Even and odd extensions: half range series . . . . . . . . 27
2.4 Fourier’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 Calculus and Fourier Series . . . . . . . . . . . . . . . . . . . . . 31
2.5.1 Differentiation of Fourier Series . . . . . . . . . . . . . . 31
2.5.2 Integration of Fourier Series . . . . . . . . . . . . . . . . 33
2.6 Complex form of Fourier Series . . . . . . . . . . . . . . . . . . 34
2.7 Parseval’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.8 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . 36
3 Fourier transform 39
3.1 An example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.4 The response function . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
1
CONTENTS CONTENTS
4 Laplace transforms 50
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Definition and existence . . . . . . . . . . . . . . . . . . . . . . 51
4.2.1 Existence . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 General properties . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3.1 Transforming derivatives . . . . . . . . . . . . . . . . . . 53
4.4 Inverse Transform . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5 Special functions . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5.1 Heaviside function . . . . . . . . . . . . . . . . . . . . . 54
4.5.2 The Dirac δ-function . . . . . . . . . . . . . . . . . . . . 56
4.6 Convolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.7 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.7.1 Solving ODEs . . . . . . . . . . . . . . . . . . . . . . . 58
4.7.2 Systems of ODEs . . . . . . . . . . . . . . . . . . . . . . 59
4.7.3 PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Wave equation 62
5.1 Classification of PDEs . . . . . . . . . . . . . . . . . . . . . . . 62
5.2 Small amplitude vibrations . . . . . . . . . . . . . . . . . . . . . 63
5.2.1 Equations of motion . . . . . . . . . . . . . . . . . . . . 64
5.2.2 Energy of vibrating string . . . . . . . . . . . . . . . . . 65
5.3 D’Alembert’s solution . . . . . . . . . . . . . . . . . . . . . . . 66
5.3.1 Travelling waves . . . . . . . . . . . . . . . . . . . . . . 67
5.3.2 Wave reflection and transmission . . . . . . . . . . . . . . 67
5.4 Separation of variables . . . . . . . . . . . . . . . . . . . . . . . 69
5.4.1 Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.4.2 Solving the separated form: spatial dependence . . . . . . 70
5.4.3 Solving the separated form: time dependence . . . . . . . 72
5.4.4 Normal modes . . . . . . . . . . . . . . . . . . . . . . . 72
5.4.5 Superposition . . . . . . . . . . . . . . . . . . . . . . . . 73
5.4.6 Additional notes . . . . . . . . . . . . . . . . . . . . . . 74
6 Heat equation 76
6.1 Heat conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.2 Separation of Variables . . . . . . . . . . . . . . . . . . . . . . . 78
6.3 More complex material properties . . . . . . . . . . . . . . . . . 80
6.3.1 Separation of variables . . . . . . . . . . . . . . . . . . . 80
6.4 More complex boundary conditions . . . . . . . . . . . . . . . . 81
6.5 Inhomogeneous problems . . . . . . . . . . . . . . . . . . . . . . 83
6.5.1 Inhomogeneous equations . . . . . . . . . . . . . . . . . 83
6.5.2 Inhomogeneous boundary conditions . . . . . . . . . . . 84
2
CONTENTS CONTENTS
7 Laplace’s equation 87
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.1.1 Boundary conditions . . . . . . . . . . . . . . . . . . . . 88
7.2 Separable solutions . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.2.1 Cartesians . . . . . . . . . . . . . . . . . . . . . . . . . . 89
7.2.2 Other domains . . . . . . . . . . . . . . . . . . . . . . . 91
8 Vector Calculus 93
8.1 Vectors and scalars . . . . . . . . . . . . . . . . . . . . . . . . . 93
8.1.1 Curves and vector valued functions of 1-variable . . . . . 93
8.1.2 Scalar Fields . . . . . . . . . . . . . . . . . . . . . . . . 95
8.1.3 Vector fields . . . . . . . . . . . . . . . . . . . . . . . . 98
8.1.4 Directional derivatives . . . . . . . . . . . . . . . . . . . 99
8.1.5 Gradient and steepest ascent . . . . . . . . . . . . . . . . 101
8.1.6 Gradients and level surfaces . . . . . . . . . . . . . . . . 101
8.1.7 Properties of gradφ . . . . . . . . . . . . . . . . . . . . . 104
8.1.8 The divergence of a vector field . . . . . . . . . . . . . . 104
8.1.9 The Laplacian . . . . . . . . . . . . . . . . . . . . . . . . 106
8.1.10 The curl of a vector field . . . . . . . . . . . . . . . . . . 107
8.1.11 Vector Identities . . . . . . . . . . . . . . . . . . . . . . 110
8.2 Integrals of vector and scalar fields . . . . . . . . . . . . . . . . . 112
8.2.1 Line integrals . . . . . . . . . . . . . . . . . . . . . . . . 112
8.2.2 Conservative vector fields . . . . . . . . . . . . . . . . . 118
8.2.3 Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.2.4 The normal to a surface . . . . . . . . . . . . . . . . . . . 121
8.2.5 Surface integrals . . . . . . . . . . . . . . . . . . . . . . 122
8.2.6 Flux integrals . . . . . . . . . . . . . . . . . . . . . . . . 125
8.3 Integral Theorems . . . . . . . . . . . . . . . . . . . . . . . . . . 130
8.3.1 Stokes’s Theorem . . . . . . . . . . . . . . . . . . . . . . 130
8.3.2 Conservative Vector Fields . . . . . . . . . . . . . . . . . 134
8.3.3 Finding potentials for conservative vector fields . . . . . . 137
8.3.4 Stokes’s theorem in the plane . . . . . . . . . . . . . . . 138
8.3.5 The proof of Green’s theorem . . . . . . . . . . . . . . . 141
8.3.6 The divergence theorem . . . . . . . . . . . . . . . . . . 142
8.3.7 Vector potential . . . . . . . . . . . . . . . . . . . . . . . 146
3
Chapter 1
Ordinary Differential Equations
1.1 Second order equations

1.1.1 Linear Second Order ODEs
The general linear second order ODE can be written as
d2 y dy
A(x) + B(x) + C(x)y + D(x) = 0 (1.1)
dx2 dx
Taking the D(x) term onto the RHS and dividing by A(x) we may write such an
equation in the standard form
d2 y dy
2
+ p(x) + q(x)y = r(x). (1.2)
dx dx
where p(x) = B(x)/A(x), q(x) = C(x)/A(x) and r(x) = −D(x)/A(x).
Homogeneous Equations
In the case where r(x) ≡ 0 the equation becomes
d2 y dy
2
+ p(x) + q(x)y = 0. (1.3)
dx dx
and is called homogeneous. For such homogeneous equations the solution will be
of the form
y(x) = c1 y1 (x) + c2 y2 (x) (1.4)
where y1 (x) and y2 (x) are linearly independent solutions of the homogeneous
equation (1.3) and c1 and c2 are constants. The functions y1 (x) and y2 (x) are
linearly independent if
c1 y1 (x) + c2 y2 (x) ≡ 0 ⇔ c1 = 0 and c2 = 0. (1.5)
4
ORDINARY DIFFERENTIAL EQUATIONS SECOND ORDER EQUATIONS
The Wronskian determinant

The following section is NOT examinable. All non-examinable material will
be in red. Examainable material starts again with black text.
If we have two solutions y1 and y2 of the homogeneous linear ODE (1.3), how do
we know if they are linearly independent? One way to motivate this is to note that
if they were linearly dependent then
y1
= const. (1.6a)
y2
0
y1
⇒ =0 (1.6b)
y2
⇒ y1 y20 − y10 y2 = 0. (1.6c)
We rewrite this using the Wronskian determinant W [y1 , y2 ](x) which is defined
as follows:
W (x) = W [y1 , y2 ](x) (1.7a)

y1 y2
= 0
(1.7b)
y1 y20
= y1 y20 − y10 y2 . (1.7c)
Using the Wronskian one can show that two solutions y1 and y2 are linearly inde-
pendent if, and only if, W (x) 6= 0.
In addition, we can straightforwardly check by direct substitution into (1.3)
that the Wronskian obeys
dW
+ p(x)W = 0. (1.8)
dx
This follows from the explicit calculation below:
dW
= (y1 y200 + y10 y20 ) − (y100 y2 + y10 y20 ) (1.9a)
dx
= y1 y200 − y100 y2 (1.9b)
0 0
= −p(x) (y1 y2 − y2 y1 ) (1.9c)
where we have used that y1,2 are solutions of the homogeneous equation (1.3)
= −p(x)W. (1.9d)
As the simple equation (1.8) for the Wronskian is separable we can write the
solution as
Z Z Z
dW
= − p(x)dx ⇒ log W = − p(x)dx + C (1.10)
W
5
So that we may write r

W = Ke− p dx
(1.11)
C
where K = e is a constant of integration. From this result it immediately follows
that the Wronskian is either zero for all x in the range considered (if K = 0), or
never zero for any x in the range considered (if K 6= 0).
Homogeneous constant coefficient

After looking at the general theory we now turn to some examples. An impor-
tant class of second order ODEs are the simple linear, homogeneous, constant
coefficient equations
d2 y dy
a 2 + b + cy = 0, (1.12)
dx dx
where a, b and c are constants. We need to find linearly independent solutions. To
do this we assume that they take the form y ∝ eλx . Substituting this guess into
equation (1.12) gives the auxiliary equation
aλ2 + bλ + c = 0. (1.13)
Depending on the type of the roots of the auxiliary equation (real, complex,
repeated) there are different types of solution:
1. b2 − 4ac > 0: Real roots λ1 and λ2 giving independent solutions
y1 = eλ1 x , y2 = eλ2 x (1.14)
and general solution
y = c1 eλ1 x + c2 eλ2 x (1.15)
2. b2 − 4ac < 0: Complex conjugate roots λ = α ± jβ giving independent
(real) solutions
y1 = eαx cos(βx), (1.16a)
y2 = eαx sin(βx). (1.16b)
y(x) = (c1 sin(βx) + c2 cos(βx))) eαx (1.17)
3. b2 − 4ac = 0: Repeated root λ giving one solution eλx . Direct calculation

shows that a second solution can be computed to give
y1 = eλx , (1.18a)
λx
y2 = xe . (1.18b)
y(x) = (c1 + c2 x)eλx (1.19)
6
Euler equations
The Euler type of linear, homogeneous equation are also straightforward. These
have the form
d2 y dy
x2 2 + cx + d y = 0 (1.20)
dx dx
where c, d are constants. Again we need to find two independent solutions. We
guess that such solutions are y ∝ xk which, on substituting into equation (1.20),
gives the condition
k 2 + (c − 1)k + d = 0. (1.21)
If equation (1.21) has distinct real roots k1 and k2 then the general solution to
the Euler equation (1.20) can be written
y = c1 xk1 + c2 xk2 . (1.22)
If equation (1.21) has repeated or complex roots it is better to use the change
of variables t = log(x) to transform equation (1.20) to the constant coefficient
problem
d2 y dy
2
+ (c − 1) + d y = 0 (1.23)
dt dt
which can be solved using the techniques of section 1.1.1.
Reduction of order and the Wronskian

(Not examined)
We have seen that in general the solution of a second order equation requires
finding two independent solutions of the homogeneous equation. There are a num-
ber of situations where it is possible to find one solution of the equation (often a
simple polynomial) but it is not easy to find a second solution. We show below a
method of using one given solution to find a second solution. This makes use of
the fact that the Wronskian is given by (1.11).
Assume that we have somehow found one solution y1 of the general linear
homogeneous ODE given by equation (1.3). We assume that the second solution
y2 is related to y1 by
y2 (x) = v(x)y1 (x). (1.24)
This immediately implies that
y20 = v 0 y1 + vy10 (1.25)
7
and from this we can compute the Wronskian determinant (see section 1.1.1)

y1 vy1

W [y1 , y2 ](x) = 0 0 0
(1.26a)
y1 v y1 + vy1
= v 0 y12 (1.26b)
r
= Ke− p dx
(1.26c)
where the last equation follows from equation (1.11).
This allows us to directly compute v, the function relating the two solutions,
as
w e− r p dx
v=K dx + D (1.27)
y12
and hence the second solution is given by
y2 (x) = v(x)y1 (x) (1.28a)
r
w e− p dx
= Dy1 (x) + Ky1 (x) dx (1.28b)
y12
and, as we do not care about proportionality constants, we can write this as
w e− r p dx
y2 (x) = y1 (x) dx. (1.28c)
y12
(Example using Legendre polynomials)
1.1.2 Inhomogeneous Linear Second Order ODEs

If r(x) 6≡ 0 we say that the differential equation is inhomogeneous.
d2 y dy
2
+ p(x) + q(x)y = r(x) (1.29)
dx dx
We often call r(x) a source term for reasons that will be clear from the examples.
In order to find the general solution to an inhomogeneous equation we first
need to find a particular solution yp (x) (sometimes called a particular integral in
some textbooks). That is some (any) solution of the inhomogeneous equation. If
one then looks at the function ỹ(x) = y(x) − yp (x) then it follows that it satisfies
the homogeneous equation
d2 ỹ dỹ
2
+ p(x) + q(x)ỹ = 0. (1.30)
dx dx
We can therefore write the general solution in the form
ỹ(x) = c1 y1 (x) + c2 y2 (x) (1.31)
8
where y1 (x) and y2 (x) are solutions of the homogeneous equation

d2 y dy
2
+ p(x) + q(x)y = 0. (1.32)
dx dx
It follows that we can write the (general) solution of the inhomogeneous equa-
tion (1.2) in the form
y(x) = yp (x) + c1 y1 (x) + c2 y2 (x). (1.33)
So the only question remaining is how to find the particular integral solution of
the full inhomogeneous problem in equation (1.29).
Finding a particular solution: undetermined coefficients

All the techniques used above concentrated on homogeneous linear ODEs. We
now look at methods for finding particular solutions of the homogeneous equation.
The idea is to use a trial solution (i.e., make an informed guess for the solution)
which is based on the form of the inhomogeneous term in equation (1.29), r(x).
We leave a number of undetermined coefficients in this guess, and substitute the
result into the full equation (1.29). If this form is a potential solution there will be
a non-trivial solution for the coefficients.
Standard simple trial solutions are:
r(x) Trial function
Any polynomial degree n Most general polynomial degree n
eαx Ceαx
sin(αx) or cos(αx) C sin(αx) + D cos(αx)
Note that if the trial function can be written as a linear combination of the
solutions to the homogeneous equation then, of course, the trial function will only
be a solution of the homogeneous equation. In this case it is usual to multiply the
trial solution by the independent variable until the trial solution is independent of
the complementary functions.
Examples!
Variation of parameters
(Not examined)
The method of undetermined coefficients is simple and effective when the in-
homogeneous term r(x) is simple, but useless in general. A more general method
can be found by assuming that the particular integral yp (x) is related to the solu-
tions y1 and y2 of the homogeneous problem, by
yp (x) = v1 (x)y1 (x) + v2 (x)y2 (x). (1.34)
9
Here v1 and v2 are functions to be determined that allow us to compute the un-
known particular integral. Note that we have introduced two functions v1 and v2 in
order to find a single function yp , so that we have an additional degree of freedom.
We will use this later on to impose an additional condition on v1 and v2 which
will reduce things back to one degree of freedom in such a way as to simplify the
calculation.
To find v1 and v2 We calculate the first and second derivatives of yp and sub-
stitute into (1.2).
yp0 = (v1 y10 + v2 y20 ) + (v10 y1 + v20 y2 ). (1.35)
In order to simplify the expression for the second derivative we now use our ad-
ditional degree of freedom to impose the condition that the term in the second
bracket vanishes so that
v10 y1 + v20 y2 = 0. (1.36)
With this condition (1.35) becomes
yp0 = (v1 y10 + v2 y20 ). (1.37)
We can then differentiate again to obtain
yp00 = (v1 y100 + v2 y200 ) + (v10 y10 + v20 y20 ). (1.38)
Substituting into the inhomogeneous ODE (1.29) and remembering that y1 and y2
are solutions of the homogeneous ODE (1.32) we can see that
r(x) = yp00 + pyp0 + qyp (1.39a)

:0 :0
00 0 00 0 0 0 0 0
= v1 (y + py + qy1 ) +v (y
2 + py
+ qy2 ) +(v1 y1 + v2 y2 ). (1.39b)
| 1 {z1 | 2 {z2

} }
from equation (1.32) from equation (1.32)
Hence the unknowns v1 and v2 are completely determined from our assumption
in equation (1.36) and our condition in equation (1.39) which can be written as
v10 y1 + v20 y2 = 0, (1.40a)

v10 y10 + v20 y20 = r. (1.40b)
The unknown functions v1,2 can thus be determined by solving equation (1.40)
0
as a linear system for v1,2 and integrating the results. We find that
ry2
v10 = − (1.41a)
y1 y20
− y10 y2
r(x)y2 (x)
=− , (1.41b)
W [y1 , y2 ](x)
10
ORDINARY DIFFERENTIAL EQUATIONS BVPS AND IVPS
and similarly
r(x)y1 (x)
v20 = . (1.42)
W [y1 , y2 ](x)
This thus gives
w r(x)y2 (x)
v1 = − dx, (1.43a)
W [y1 , y2 ](x)
w r(x)y1 (x)
v2 = dx (1.43b)
W [y1 , y2 ](x)
and the particular integral of the original inhomogeneous ODE (1.29) follows
from equation (1.34).
1.2 Boundary and Initial Value Problems

We have already seen that the differential equation
d2 y dy
+ p(x) + q(x)y = r(x) (1.44)
dx2 dx
has general solution y(x) = yp (x) + c1 y1 (x) + c2 y2 (x) where yp (x) is a particular
solution of the equation and y1 (x) and y2 (x) are independent solutions of the
homogeneous equation and c1 and c2 are constants.
In order to obtain a unique solution to the differential equation we need some
way of fixing the constants c1 and c2 . One common way of doing this is to fix the
value of both y(x) and y 0 (x) at some point x0 . So that
y(x0 ) = C, (1.45a)
y 0 (x0 ) = D. (1.45b)
where C and D are given. Equation (1.44) together with the initial conditions
(1.45) is called an initial value problem.
Rather than give the value of the function y and its derivative y 0 at a single point
x0 one can give alternatively give the value of y at two point x0 and x1 , so that
the solution satisfies the differential equation (1.44) together with the boundary
conditions
y(x0 ) = α, (1.46a)
y(x1 ) = β. (1.46b)
11
This is an example of a boundary value problem. More generally we can specify

the value of some linear combination of the function and its derivatives at two pint
x0 and x1 and demand that
a0 y(x0 ) + b0 y 0 (x0 ) = α, (1.47a)

a1 y(x1 ) + b1 y 0 (x1 ) = β. (1.47b)
One can then show that the boundary value problem of solving the differential
equation (1.44) on the region x0 ≤ x ≤ x1 subject to the boundary conditions
(1.47) has a unique solution provided the associated homogeneous boundary prob-
lem (i.e. the case where r(x) ≡ 0, α = 0, β = 0) only has the trivial solution (this
result is known as the Fredholm Alternative)
1.2.1 Eigenvalues and functions

We now look at a class of BVPs where the point is to find the conditions under
which the problem has a solution. We start by looking at a simple example where
we try and find non-trivial solutions (i.e. solutions that are not zero everywhere)
of the simple BVP
d2 y
+ λy = 0, (1.48a)
dx2
y(0) = 0, (1.48b)
y(π) = 0. (1.48c)
The unknown parameter λ may take any value. However, it may not be true that
the BVP given by equation (1.48) has a solution for any value of λ. The task of
finding a solution of the BVP involves finding not just solutions y that satisfy the
differential equation (1.48a), but also of finding which values of λ allow such a
solution to satisfy the boundary conditions given by equations (1.48b-1.48c).
In this simple case we can solve the BVP by inspection. However, to proceed
more systematically we split the problem into three separate cases depending on
the sign of the parameter λ:
1. λ = −α2 < 0. We first find solutions of equation (1.48a). As the differential
equation is homogeneous, these combine to form √ the general solution. In
this case the auxiliary equation has real roots ± −λ = ±α. It follows that
the general solution is
y(x) = a1 eαx + a2 e−αx . (1.49)
The only solution compatible with the boundary conditions (1.48b-1.48c) is

the trivial solution y ≡ 0 when a1 = 0 = a2 .
12
2. λ = 0. In this case the equation (1.48a) obviously has the general solution
y(x) = b1 x + b2 . (1.50)
Again, the only solution compatible with the boundary conditions (1.48b-
1.48c) is the trivial solution y ≡ 0 when b1 = 0 = b2 .
3. λ √= µ2 > 0. In this case the auxiliary equation has imaginary roots

± −λ = ±jµ. It follows that the general solution can be written as
y(x) = c1 sin (µx) + c2 cos (µx) . (1.51)
From the left boundary condition, equation (1.48b), we must have that c2 =
0 and so this reduces to
y(x) = c1 sin (µx) . (1.52)
In order to satisfy the final boundary condition, equation (1.48c), we either

have c1 = 0 and hence the trivial solution y ≡ 0 again, or we have a
condition on µ that ensures that
sin (µπ) = 0. (1.53)
This condition is obviously that µ = n for any integer n1 ; given our previous
definitions, this gives the condition on λ that
λ = n2 , n = 1, 2, . . . (1.54)
Therefore the non-trivial solutions of the BVP (1.48) are
yn (x) = sin (nπx) . (1.55)
The values of λ for which a solution of the BVP can be obtained are denoted λn ,
in this case are
λn = n2 , n = 1, 2, . . . , (1.56)
and are called the eigenvalues of the problem. The associated solutions yn (x) are
called the eigenfunctions. Clearly the eigenvalues are uniquely determined but
the eigenfunctions are not; any constant multiple of an eigenfunction yn is also a
solution of the BVP, and is also an eigenfunction.
A number of points should be noted for later consideration.
1
The condition (1.53) is satisfied with µ = ±n, where n is a positive integer. Since
sin (−nπx) = − sin (nπx), we may restrict to µ = n > 0 without loss of generality.
13
1. The eigenvalues are all real.
2. The different eigenfunctions are orthogonal in this sense:

wπ
sin(mx) sin(nx) dx = 0, if m 6= n, (1.57a)
0
wπ
cos(mx) cos(nx) dx = 0, if m 6= n, (1.57b)
0
wπ
sin(mx) cos(nx) dx = 0, for all m, n, (1.57c)
0
3. The eigenvalues form an increasing sequence of positive numbers that ap-

proach ∞.
4. The nth eigenfunction (sin(nx)) has exactly n − 1 zeros inside the interval
[0, π] on which the BVP is defined. It also vanishes at the endpoints of the
interval.
1.2.2 Sturm-Liouville form

(Not Examined)
We consider the equation
d2 y dy
2
+ b(x) + c(x)y = λd(x)y. (1.58)
dx dx
This is a minor modification of previous second order ODEs that we have looked
at. We shall later specify boundary conditions in order to turn it in to a BVP.
r
A standard form is found by multiplying the whole equation by e b(x) dx to
rewrite equation (1.58) as
r 0 r r
e b(x) dx y 0 + c(x)e b(x) dx y(x) = λd(x)e b(x) dx y(x). (1.59)
This is called the self-adjoint form of the ODE. The conventional notation is
0
(p(x)y 0 ) − q(x)y + λw(x)y = 0. (1.60)
This is often written as

Ly = λw(x)y (1.61)
14
where w(x), the weight function, is positive (w ≥ 0) and non-trivial (w 6≡ 0), and
the differential operator L is given by
0
L : y → (p(x)y 0 ) − q(x)y. (1.62)
A solution yn of the regular Sturm-Liouville problem given by equation (1.60)

(with suitable boundary conditions) is called an eigenfunction of the problem, and
the associated value of λ, λn is called the eigenvalue. As the notation suggests,
we are expecting that more than one such solution exists.
Self-adjoint operators and boundary conditions

A key property is that the differential operator L is self-adjoint. This is true on x ∈
[x0 , x1 ] if, and only if, for all pairs of functions yA , yB satisfying the appropriate
boundary conditions,
wx1 wx1
yA LyB dx = yB LyA dx. (1.63)
x0 x0
For the Sturm-Liouville equation (1.60) we find that

wx1 wx1
0 0
− (pyB0 ) yA + qyA yB + (pyA0 ) yB − qyB yA dx

(yA LyB − yB LyA ) dx =
x0 x0
(1.64a)
wx1
x
= [−pyA yB0 + pyB yA0 ]x10 + [pyB0 yA0 − pyA0 yB0 ] dx
x0
(1.64b)
x
= [−pyA yB0 + pyB yA0 ]x10 (1.64c)
We want this to vanish in order to have a self-adjoint problem. We therefore have

a self-adjoint problem if, and only if,
x
[−pyA yB0 + pyB yA0 ]x10 = 0. (1.65)
This can be guaranteed if at x = x0 and x = x1 we insist on some appropriate

combination of the following boundary conditions:
1. y = 0: Dirichlet condition.
2. y 0 = 0: Neumann condition.
3. y + ky 0 = 0: Radiation condition.
15
4. p = 0: x = x0 , x1 are singular points of the ODE.
We note that the conditions at the boundaries x = x0 , x1 need not be the same.
There is one final possibility: periodic boundary conditions, where
y(x0 ) = y(x1 ), y 0 (x0 ) = y 0 (x1 ), p(x0 ) = p(x1 ). (1.66)
Any of these are Sturm-Liouville boundary conditions, and given these the differ-
ential operator L is self-adjoint.
Properties of solutions
Properties 1. The eigenvalues λn are real.
Proof. Suppose that the eigenvalue λ, and possibly also the eigenfunction y are
complex. By definition
Ly = λw(x)y. (1.67)
If we take the complex conjugate then, as both the operator L and the weight
function w(x) are real, we have
Lȳ = λ̄w(x)ȳ. (1.68)
Taking the difference between these two equations and integrating over the do-
main gives us
wx1 wx1
(ȳLy − yLȳ) dx = (λ − λ̄) w(x)y ȳ dx. (1.69)
x0 x0
As the operator is self-adjoint, we have from equation (1.63) that the left hand
side of equation (1.69) vanishes identically.
As the eigenfunction solution is non-trivial we can see that
y ȳ = |y|2 > 0, (1.70)
and by assumption we have that w(x) > 0. Therefore in order for the right hand
side to also vanish, we must have that λ = λ̄, and hence the eigenvalues must be
real. Without loss of generality, we can therefore also take the eigenfunctions to
be real.
Properties 2. The eigenfunctions yn are real.
16
Proof. This is essentially given by the previous argument, as noted above, but an
alternative is as follows.
Suppose that yn = fn (x) + jgn (x) where both f, g are real and the eigenfunc-
tion yn is non-trivial. As the differential equation (and hence operator L) are real
it follows that
L (fn + jgn ) = L (fn ) + jL (gn ) (1.71a)

= −λn w(x) (fn + jgn ) . (1.71b)
By equating real and imaginary parts in (1.71) it follows that
L (fn ) = −λw(x)fn , (1.72a)

L (gn ) = −λw(x)gn , (1.72b)
which, after checking the boundary conditions and using that the original BVP is
real, gives that both fn and gn must therefore be solutions of the original Sturm-
Liouville problem. Therefore we can take the non-zero real function of fn , gn to
be the real-valued eigenfunction with eigenvalue λn (at least one must be non-zero
as yn is non-trivial).
Properties 3. Eigenfunctions corresponding to distinct eigenvalues are orthogo-
nal.
Proof. Suppose that we have distinct eigenfunction solutions ym , yn with asso-
ciated eigenvalues λm 6= λn . Again we can use the self-adjoint nature of the
problem, given by equation (1.63), to show that
wx1
0= (ym Lyn − yn Lym ) dx (1.73a)
x0
wx1
= (λn − λm ) w(x)ym yn dx. (1.73b)
x0
Therefore, as the eigenvalues are distinct, we have that

wx1
w(x)ym yn dx = 0. (1.74)
x0
We say that the associated eigenfunctions ym , yn are orthogonal with respect to

the weight function w(x).
Properties 4. There are infinitely many eigenvalues λ1 < λ2 < . . . that approach
infinity, and the associated eigenfunction yn (x) has exactly (n − 1) zeros in the
interval (x0 , x1 ).
17
Properties 5. The eigenfunctions are complete in the sense that for any bounded,
piecewise continuous function F (x) on [x0 , x1 ] the error N
wx1
" N
#
X
N = w(x) F (x) − an yn (x) dx (1.75)
x0 n=1
becomes arbitrarily small as N → ∞.
18
Chapter 2
Fourier Series
2.1 Fourier series

Motivated by the results for boundary value problems discussed in section 1.2
we now look at how we can represent a particular function as a series of other
functions; in particular, as a series of trigonometric functions.
We start by considering periodic functions. Periodic functions occur fre-
quently in engineering problems and their representation in terms of simple pe-
riodic functions, which leads to Fourier series, is of great practical importance.
They are in a sense more general than Taylor series because they permit the rep-
resentation of discontinuous periodic functions which occur in many engineering
applications.
A Fourier series in simple terms is an expansion of a function in terms of cos
and sin functions. Suppose f (x) is a periodic function of period 2π so that
f (x + 2π) = f (x)
Then (provided f (x) is a “fairly nice” function that satisfies certain conditions) it
may be shown that it can be expanded in a Fourier series
∞
1 X
f (x) = a0 + [an cos nx + bn sin nx] (2.1)
2 n=1
where an , bn are constants called the Fourier Coefficients and are such that the
series on the right hand side converges for all x.
More generally suppose that f (x) is a periodic function of period 2L so that
f (x + 2L) = f (x)
Then by using the result above and rescaling the x-coordinate we may expand the
nπx
function in a Fourier series in terms of the 2L-periodic functions cos L and
19
FOURIER SERIES ORTHOGONALITY
nπx

sin L
so that we may write
∞ h
1 X nπx nπx i
f (x) = a0 + an cos + bn sin (2.2)
2 n=1
L L
The question then arises: given a 2L-periodic function f (x), how can we find the
Fourier coefficients an , bn in the Fourier Series?
2.2 Orthogonality
Throughout this section (and indeed this chapter) the following identities are re-
peatedly used; they are sufficiently important that you should know them.
sin n + 12 π = (−1)n ;

sin (nπ) = 0, (2.3a)
n 1

cos (nπ) = (−1) , cos n + 2 π = 0. (2.3b)
Consider the integral

wL mπx nπx
Imn = sin sin dx. (2.4)
L L
−L
We can straightforwardly use the formula cos(A+B) = cos A cos B −sin A sin B
and cos(A − B) = cos A cos B + sin A sin B to show that
2 sin A sin B = cos(A − B) − cos(A + B)
This enables us to write the integral as

wL mπx nπx
Imn = sin sin dx (2.5a)
L L
−L
1w
L
(m − n)πx (m + n)πx
= cos − cos dx (2.5b)
2 L L
−L
 (m−n)πx L
(m+n)πx
L  sin L
sin L
= −  (2.5c)
2π m−n m+n
−L
= 0 except when m = n. (2.5d)
In the special case when m = n it is straightforward to check (exercise) that
Inn = L. (2.6)
20
FOURIER SERIES FOURIER COEFFICIENTS
We can combine these results into one simple equation:

Imn = Lδmn (2.7)
where the symbol δmn represents
(
1 m = n,
δmn = . (2.8)
0 m=6 n
Similar use of the double angle formula give the results
wL mπx nπx
cos cos dx = Lδmn , (2.9)
L L
−L
wL mπx nπx
sin cos dx = 0. (2.10)
L L
−L
Note that the cos relation given by equation (2.9) does not hold for m = 0; in that
case the result is
wL nπx
cos dx = 2Lδ0n (2.11)
L
−L
Note: If we consider
wL
hf (x), g(x)i = f (x)g(x) dx (2.12)
−L
to be an inner product between the functions f and g, then the set of functions
n nπx nπx o
1, sin , cos (2.13)
L L
form an orthogonal basis for the space of 2L-periodic “nice” functions. In fact
we know from property 3 of the general theory of Sturm-Liouville boundary value
problems (even without doing the calculations above) that these functions are or-
thogonal. For Fourier series it is just as easy to do the explicit calculations but
for other types of expansion (e.g. Fourier-Bessel series which are used to describe
vibrating plates) it is much easier to use the general theory.
2.3 Fourier Coefficients

We remember the definition of the Fourier Series:
∞ h
1 X nπx nπx i
f (x) = a0 + an cos + bn sin . (2.2)
2 n=1
L L
21
If we multiply this equation by the appropriate sin term,

mπx
sin (2.14)
L
and integrate the function over its period 2L, then the orthogonality relations given
by equations (2.7, 2.10) give
wL mπx
Lbm = f (x) sin dx. (2.15)
L
−L
The orthogonality relations ensure that except for the single term bm , all other
terms vanish under integration. Similarly, using the appropriate cos term and the
orthogonality relations given by equations (2.9, 2.10) gives
wL mπx
Lam = f (x) cos dx, (2.16)
L
−L
and this equation holds even for m = 0 (due to the factor of 1/2 in the definition
of the Fourier Series, and using equation (2.11)).
We can summarise this to give the Euler relations:
1 w
L mπx
am = f (x) cos dx, (2.17a)
L L
−L
1 w
L mπx
bm = f (x) sin dx. (2.17b)
L L
−L
Note: Because the function f (x) is periodic, and the range of integration is a
single period, then one could equally use any 2L-periodic range in x; for example,
r 2L
one could equally well use 0 dx.
Note: Because the formulae for the Fourier coefficients an , bn are explicit, the
Fourier Series of a function is unique.
2.3.1 Even and odd functions

The function f (x) is an even function if
f (−x) = f (x)
22
Figure 2.1: Even function. The function is symmetric about the

y-axis.
for all x. Graphically this means that the function is symmetric about the y-axis.
The symmetrical integral of an even function can be expressed as follows
Z` Z`
f (x)dx = 2 f (x)dx
−` 0
The function f (x) is said to be an odd function if
f (−x) = −f (x)
for all x. Graphically this means that the function is symmetric under a reflection
about the y-axis and the x-axis.
The symmetrical integral of an odd function is zero i.e.
Z`
f (x)dx = 0
−`
Note, in addition, the following important properties
23
Figure 2.2: Odd function. The function is symmetric under a

reflection about the y-axis and the x-axis.
• If f (x), g(x) are even then f (x)g(x) is even since

f (−x)g(−x) = f (x)g(x)
• If f (x), g(x) are odd then f (x)g(x) is even since

f (−x)g(−x) = [−f (x)][−g(x)] = f (x)g(x)
• If f (x) is even and g(x) is odd then f (x)g(x) is odd since

f (−x)g(−x) = f (x)[−g(x)] = −f (x)g(x)
If f (x) is an even function then this gives rise to a Fourier cosine series only
since
Zπ Zπ
1 2
an = f (x) cos nxdx = f (x) cos nxdx
π π
−π 0
Zπ
1
bn = f (x) sin nxdx = 0
π
−π
and so
∞ Zπ
a0 X 2
f (x) = + an cos nx with an = f (x) cos nxdx
2 n=1
π
0
24
If f(x) is an odd function then this gives rise to a Fourier sine series only since
Zπ
1
an = f (x) cos nxdx = 0
π
−π
Zπ Zπ
1 2
bn = f (x) sin nxdx = f (x) sin nxdx
π π
−π 0
and so
∞ Zπ
X 2
f (x) = bn sin nx with bn = f (x) sin nxdx
n=1
π
0
Example.
Find the Fourier series of f (x) which is equal to x in the range −π < x < π
and periodic outside, see Fig. 2.3.
Since f (x) is an odd function, it has a Fourier sine series
∞
X
f (x) = bn sin nx
n=1
with
Zπ
2
bn = f (x) sin nxdx
π
0
Zπ
2
= x sin nxdx
π
0
Zπ
2 cos nx
= xd −
π n
0
Zπ
2 − cos nx 2
= x + cos nxdx
π n nπ
0 π
2h π n
i 2 sin nx
= − (−1) +
π n nπ n 0
2
= (−1)n−1
n
25
4
Sawtooth function
-1
-2
-3
-4
-8 -6 -4 -2 0 2 4 6 8
x
Figure 2.3: The function f (x) = x with range −π < x < π,

periodically extended to the real line (Sawtooth function).
Therefore
∞
X 2
f (x) = (−1)n−1 sin nx
n=1
n
1 1
= 2[sin x − sin 2x + sin 3x − · · · ]
2 3
If we consider the contribution of the first few terms to the sum we can gain some
idea of how the Fourier series approximates to the function in the interval 0 to π.
26
Then, in particular, at the following points

∞
X 2
x=0 : (−1)n−1 sin 0 = f (0) = 0
n=1
n
∞
X 2 1 1
x=π : (−1)n−1 sin nπ = [f (π − 0) + f (π + 0)] = [π + (−π)] = 0
n=1
n 2 2
∞
π X 2 nπ 1 1 π π
x= : (−1)n−1 sin = 2[1 − + − · · · ] = f ( ) =
2 n=1
n 2 3 5 2 2
This last result can be used to show that the sum of the series
1 1 1π π
1 − + − ··· = =
3 5 22 4
2.3.2 Even and odd extensions: half range series

Suppose that f (x) is given only on 0 ≤ x ≤ L; Then we can extend the function
f onto −L < x ≤ L by making
f (−x) = −f (x).
We then extend this to all values of x by making the function have period 2L, see
Fig. 2.4.
This is called an odd periodic extension of the function. Since the function is
odd it has a Fourier series which only contains the sine terms. Hence the function
is represented by
X∞ nπx
f (x) = bn sin , (2.18)
n=1
L
a Fourier sine series where
an = 0, (2.19a)
2 wL mπx
bn = f (x) sin dx, (2.19b)
L 0
L
where we have used the fact that f (x) is odd (and hence f (x) sin mπx

L
is even)
to write the expression for bn in terms of an integral over 0 < x < L (which is
half the period). For this reason such an Fourier sine series is called a half-range
Fourier series.
Alternatively given f (x) only on 0 ≤ x ≤ L we can extend the function f
onto −L < x ≤ L by making
f (−x) = f (x).
27
and then extend this to all values of x by making the function have period 2L, see
Fig. 2.4.
This is called an even periodic extension of the function. Since the function
is even it has a Fourier series which only contains the cosine terms. Hence the
function is represented by
∞
1 X nπx
f (x) = a0 + an cos , (2.20)
2 n=1
L
a Fourier cosine series where
bn = 0, (2.21a)
2w
L mπx
an = f (x) cos dx (2.21b)
L0 L
and we have again used the fact that the integrand is even to write the expression
for an as an integral over half the range.
Note: we could also extend the function by making f be periodic with period
L, which immediately implies that it is also periodic with period 2L, see Fig. 2.4.
However in applications to differential equations it is the odd and even extensions
that are most useful.
Example.
Find the Fourier series of (i) the odd periodic extension and (ii) the even peri-
odic extension of
0 < x < 2`

1,
f (x) = `
0, 2
<x<`
in the range (0, `).
1. In (0, `)
∞
X nπx
f (x) = bn sin
n=1
`
where
Z`
2 nπx
bn = f (x) sin dx
` `
0
`
Z2
2 nπx
= sin dx
` `
0
28
ff (x)
−L L t
ll
fu
f(x) fe (x)
even
x
L −L L t
od
d
fo (x)
interval of
periodicity −L L t
Figure 2.4: Periodic extensions of the function f (x) defined on

the half interval, 0 < x < L (left panel). Right panel, top: Ex-
tension as a periodic function of period L. Right panel, middle:
Extension as a even function. Right panel, bottom. Extension
as an odd function. The blue shading indicates the interval of
periodicity.
29
2 h nπx i 2`
= − cos
nπ ` 0
2 h nπ i
= 1 − cos
nπ 2
Therefore
∞
2X1h nπ i nπx
f (x) = 1 − cos sin
π n=1 n 2 `

2 πx 2πx 1 3πx 1 5πx
= sin + sin + sin + sin + ···
π ` ` 3 ` 5 `
2. In (0, `)
∞
a0 X nπx
f (x) = + an cos
2 n=1
`
where
Z`
2 nπx
an = f (x) cos dx
` `
0
`
Z2
2 nπx
= cos dx
` `
0
2 ` h nπx i 2`
= sin
` nπ ` 0
2 nπ
= sin
nπ 2
and
`
Z2
2
a0 = dx = 1
`
0
Therefore
∞
1 X 2 nπ nπx
f (x) = + sin cos
2 n=1 nπ 2 `

1 2 πx 1 3πx 1 5πx
= + cos − cos + cos − ···
2 π ` 3 ` 5 `
30
FOURIER SERIES FOURIER’S THEOREM
2.4 Fourier’s Theorem

We have said that Fourier Series work for “fairly nice” functions, without spec-
ifying what we mean. Here we note that sufficient conditions for the existence
of a Fourier Series representation, i.e. that the sum given in the definition of the
Fourier Series in equation (2.2) converges, is that
1. the function f (x) is bounded, and
2. the function has a finite number of maxima, minima and discontinuities, and
3. the function is 2L-periodic for some L.
Under these conditions we have the following:
Fourier Convergence Theorem. Let f (x) be a bounded, 2L-periodic function
with a finite number of maxima, minima and discontinuities. Then the sum
∞ h
1 X nπx nπx i
a0 + an cos + bn sin (2.22)
2 n=1
L L
converges for all finite x, and when an , bn are given by the Euler formulas (2.17)
the sum converges to f (x) if, at x, f is continuous. At points x where f is not
continuous the sum converges to 12 [f (x− ) + f (x+ )] where f (x− ) and f (x+ ) are
the values of the function as we approach from the left and right respectively.
Note: There may be “hidden” discontinuities at the end of the periodic range.
That is, if the function is 2L-periodic with f (x) = f (x + 2L) over the range
−L ≤ x < L, it may be that f (−L) 6= f (L).
2.5 Calculus and Fourier Series

2.5.1 Differentiation of Fourier Series
It is trivial to take the definition of a Fourier Series
∞ h
1 X nπx nπx i
f (x) = a0 + an cos + bn sin (2.23)
2 n=1
L L
and to differentiate or integrate the right hand side term by term. The question is,
under which conditions does the statement
∞ h
0
X nπ nπx nπ nπx i
f (x) = − an cos + bn sin (2.24)
n=1
L L L L
hold?
From the conditions given in section 2.4 we would expect that
31
FOURIER SERIES CALCULUS AND FOURIER SERIES
1. the function f (x) is bounded, and
2. the function has a finite number of maxima, minima and discontinuities, and
3. the function is 2L-periodic for some L,
to be a good start. Under this set of conditions the Fourier convergence theorem
of section 2.4 shows at least that the original definition of equation (2.23) holds
everywhere. If in addition we assume that
4. the function f (x) is continuous, and
5. its derivative f 0 (x) is piecewise smooth for all x
then application of the Fourier convergence theorem to f 0 gives that

∞ h
X nπx nπx i
0
f (x) = αn cos + βn sin (2.25)
n=1
L L
for some Fourier coefficients αn , βn . It is then straightforward to show that these

coefficients are exactly those given by equation (2.24).
Example 1.
In the range −π < x < π

1 cos 2x cos 3x
f (x) = x sin x = 1 − cos x − 2 − + ···
2 1.3 2.4
and f (−π) = f (π) = 0.

Since f (x) satisfies the three conditions above we may differentiate both sides
to give

1 2 sin 2x 3 sin 3x
sin x + x cos x = sin x + 2 − + ···
2 1.3 2.4
from which we get in addition

1 2 sin 2x 3 sin 3x
x cos x = − sin x + 2 − + ···
2 1.3 2.4
Note: It is important to note that although the conditions above are sufficient
and some of them can be relaxed, the example below shows that some of these
conditions are crucial.
Example 2.
32
FOURIER SERIES CALCULUS AND FOURIER SERIES
The periodic extension of f (x) = x in 0 < x < 2π is

∞
X sin nx
x=π−2 .
n=1
n
However
f (0) = 0 6= f (2π) = 2π
and therefore the series may not be differentiated term by term. If we do differen-
tiate the series term by term we obtain the series
∞
X
g(x) = −2 cos nx
n=1
which does not converge to f 0 (x) = 1 in 0 < x < 2π. Indeed the series is
divergent since the nnth term does not tend to zero as n → ∞.
2.5.2 Integration of Fourier Series

The Fourier series of f (x) in −` < x < `
∞
a0 X h nπx nπx i
f (x) = + an cos + bn sin
2 n=1
` `
may always be integrated term by term to give

 
Zx Zx ∞ Zx Zx
a0 X
an cos nπx dx + bn sin nπx dx
f (x)dx = dx +
2 n=1
` `
−` −` −` −`
∞
a0 X ` nπx nπ(−`)
= (x + `) + an sin − sin
2 n=1
nπ ` `
∞
X ` nπx nπ(−`)
+ bn − cos + cos
n=1
nπ ` `
∞
a0 `X1h nπx nπx n
i
= (x + `) + an sin − bn cos + (−1) bn .
2 π n=1 n ` `
The proof is straightforward but rather detailed. Of course, the constant term
in the new Fourier Series will need to be computed as it appears as the ‘constant
of integration’ when one integrates term by term.
Example.
33
FOURIER SERIES COMPLEX FORM OF FOURIER SERIES
In the interval −π < x < π, f (x) = x has odd periodic extension

∞
X sin nx
x=2 (−1)n−1 .
n=1
n
Integrating term by term from −π to x
2 x ∞
x X h cos nx ix
=2 (−1)n−1 −
2 −π n=1
n2 −π
∞
x2 π 2 X (−1)n
− =2 2
[cos nx − cos nπ]
2 2 n=1
n
∞
2 2
X (−1)n
x −π =4 [cos nx − (−1)n ]
n=1
n2
∞
2 a0 X (−1)n
x = +4 2
cos nx
2 n=1
n
where ∞
a0 X 1
= π2 − 4
2 n=1
n2
∞
X 1
We can determine a0 either by knowing the sum of the series 2
or directly
n=1
n
from the Fourier coefficient formula
Zπ π
1 2 1 x3 1 2π 3 2π 2
a0 = x dx = = . =
π π 3 −π π 3 3
−π
Hence we get
∞
π2 X (−1)n
x2 = +4 2
cos nx
3 n=1
n
2.6 Complex form of Fourier Series

We can express the Fourier series of a function more succinctly if we employ a
complex representation and use the exponential definitions for the trigonometric
functions,
nπx 1
ejnπx/L + e−jnπx/L ,

cos = (2.26a)
L 2
nπx 1 jnπx/L
− e−jnπx/L

sin = e (2.26b)
L 2j
34
FOURIER SERIES PARSEVAL’S THEOREM
and substitute these into the definition of the Fourier Series, equation (2.2), then
we will end up with
X∞
f (x) = cn ejnπx/L , (2.27)
n=−∞
where the complex coefficients cn are given by


1
 2 (an − jbn )
 n > 0,
cn = 21 a0 n = 0, , (2.28)

1
2
(a−n + jb−n ) n < 0
which implies, or can be written as,
1 w
L
cn = f (x)e−jnπx/L dx. (2.29)
2L
−L
Example.
Find the complex Fourier series for
f (x) = ex , −π < x < π with f (x + 2π) = f (x)
By above
Zπ
1
cn = ex e−jnx dx
2π
−π
π
1 e(1−jn)x

=
2π 1 − jn −π
1 1 (1−jn)π
− e−(1−jn)π

= . e
2π 1 − jn
(−1)n
= sinh π
π(1 − jn)
Therefore ∞
X (−1)n jnx
f (x) = e sinh π
−∞
π(1 − jn)
2.7 Parseval’s Theorem

(not examinable)
35
FOURIER SERIES SUMMARY OF RESULTS
Consider the real number

wL
{f (x)}2 dx. (2.30)
−L
If we expand out the Fourier series representation of f (x) given by (2.2) and use
the orthogonality relations given by equations (2.7, 2.9, 2.10) then we find
wL wL 1 nπx i 2
( ∞ h
)
X nπx
[f (x)]2 dx = a0 + an cos + bn sin dx
2 n=1
L L
−L −L
(2.31a)
wL
( ∞
)
1 2 Xh 2 nπx nπx i
= a0 + an cos2 + b2n sin2 dx
4 n=1
L L
−L
(2.31b)
" ∞
#
1 2 X
=L a + a2 + b2n . (2.31c)
2 0 n=1 n
This can be restated as
X
“kf k2 = squares of components of f on basis vectors”. (2.32)
The norm of f is often called the power contained in the function, and Parseval’s
Theorem or Parseval’s Identity, given by equation (2.31), relates the power to the
Fourier coefficients.
2.8 Summary of results

1. Fourier series of periodic and non-periodic functions
(a) Periodic function f (x) of period 2π.
The Fourier series is valid for all x.
(b) Non-periodic function f (x) defined over interval −` < x < `, say.
Periodic extension
∞
a0 X nπx nπx
f (x) = + an cos + bn sin
2 n=1
` `
where
Zl Zl
1 nπx 1 nπx
an = f (x) cos dx, bn = f (x) sin dx
` ` ` `
−` −`
36
f (x) is even −→ bn = 0
f (x) is odd −→ an = 0
(c) Non-periodic function f (x) defined over interval 0 < x < `, say.
i. Even periodic extension (half range expansion)
∞
a0 X nπx
f (x) = + an cos
2 n=1
`
where
Z`
2 nπx
an = f (x) cos dx
` `
0
ii. Odd periodic extension (half range expansion)

∞
X nπx
f (x) = bn sin
n=1
`
where
Z`
2 nπx
bn = f (x) sin dx
` `
0
2. Differentiation
The Fourier series for f (x) in −` < x < ` may be differentiated term by
term to give a Fourier series for f 0 (x) if
(a) f (x) is continuous in −` < x < `

(b) f 0 (x) has a Fourier series
(c) f (−`) = f (`)
3. Integration
The Fourier series for f (x) in −` < x < ` may always be integrated term
by term.
4. Complex Fourier series

The complex Fourier series for f (x) in −` < x < ` is
∞
X jnπx
f (x) = cn e `
−∞
37
where
Z`
1 jnπx
cn = f (x)e− ` dx
2`
−`
38
Chapter 3
Fourier transform
The Fourier series is a very powerful tool to analyse periodic functions or func-
tions defined in an interval, using their periodic extensions. However, we now
want to extend it to functions defined on the whole real line that are not periodic.
We do this by considering the case of a Fourier series of a simple function defined
on an interval of length 2L and then see what happens in the limit that L tends to
infinity.
3.1 An example
Consider the function

0 −L ≤ t < −1,

fL (t) = 1 −1 ≤ t < 1, with fL (t + 2L) = fL (t) ∀t ∈ R. (3.1)

0 1 ≤ t < L,

As L increases this function tends to the non-periodic function


0 t < −1,

f (t) = 1 −1 ≤ t < 1, (3.2)

0 t ≥ 1.

The function fL (t) has complex Fourier coefficients
1 w
+L
πn 1 πn
cn = fL (t)e−j L t dt = sin (3.3)
2L πn L
−L
and complex Fourier series

∞
X 1 nπ nπ
fL (t) = sin ej L t . (3.4)
n=−∞
nπ L
39
FOURIER TRANSFORM AN EXAMPLE
y fL(t)
−2L −L −1 +1 L 2L t
y
fL(t)
L increases
−L −1 +1 L t
y
f(t)
Limit of infinite L
−1 +1 t
Figure 3.1: The periodic function fL (t) defined in equation (3.1),

top and middle row, and its limit as L → ∞, f (t) defined by
equation (3.2), bottom row.
We want to see how equation (3.3) and (3.4) change as L tends to infinity. We
introduce the variable
nπ
ωn = (3.5)
L
and look at how this changes with n by considering
(n + 1)π nπ π
∆ωn = ωn+1 − ωn = − = . (3.6)
L L L
Note that as L increases ∆ωn tends to zero and the ωn variables becomes contin-
uous, rather than discrete. With this notation we also have that
1 ∆ωn
= . (3.7)
n ωn
We start with equation (3.4) and rewrite it as
∞
X 1 ∆ωn
fL (t) = sin(ωn )ejωn t . (3.8)
ω =−∞
π ωn
n
We now take the limit L → ∞ of this expression and obtain

f (t) ≡ lim fL (t) (3.9a)
L→∞
( ∞ )
X 1 ∆ωn
= lim sin(ωn )ejωn t (3.9b)
L→∞
ω =−∞
π ω n
n
1 w sin(ω) jωt
∞
= e dω . (3.9c)
π −∞ ω
40
FOURIER TRANSFORM DEFINITION
In this way we have expressed f (t) as the integral of the function sin(ω)/ω of the
continuous variable ω, multiplied by the complex exponential exp(jωt). We now
need to deal with equation (3.3). We define
c(ωn ) = Lcn . (3.10)
Using the same substitutions as for fL (t) we can write
1w
+L
c(ωn ) = fL (t)e−jωn t dt , (3.11)
2
−L
which, in the limit L → ∞ gives a function of the continuous variable ω,
1 w
∞
1 +1 −jωt sin(ω)
Z
−jωt
c(ω) = f (t)e dt = e = . (3.12)
2 −∞ 2 −1 ω
Note that c(ω) is the function that appears as integrand in the expression (3.9c)
of f (t). Equation (3.12) is the generalisation of (3.3) to a non-periodic function,
while equation (3.9c) is the generalisation of (3.4). We will see shortly that the
former is called the Fourier transform of f (t), while the latter is its inverse-Fourier
transform.
3.2 Definition
We are now in a position to define the Fourier transform of a function of a real
variable, f (t), with t ∈ R. We assume that the function f (t) satisfies the following
conditions:
w∞
1. |f (t)| dt < ∞, (note that this requires that |f (t)| → 0 as |t| → ∞),
−∞
2. f (t) has at most a finite number of maxima, minima and discontinuities in

any finite interval.
The function
1 w
∞
F (ω) = √ f (t)e−jωt dt . (3.13)
2π −∞
is called the Fourier transform of f (t) and we write
1 w
∞
F (ω) = F[f (t)](ω) := √ f (t)e−jωt dt. (3.14)
2π −∞
41
FOURIER TRANSFORM DEFINITION
Conversely, given F (ω) we can write f (t) as
1 w
∞
f (t) = √ F (ω)ejωt dω (3.15)
2π −∞
where this expression converges to
1. f (t) at all points where f (t) is continuous;

2. the average of right- and left-hand limits of f (t) at points where f (t) is
discontinuous.
In this respect the function f (t) is called the inverse Fourier transform of F (ω)
and we write
1 w
∞
−1
f (t) = F [F (ω)](t) := √ F (ω)ejωt dω (3.16)
2π −∞
w∞
Remark - The requirement that |f (t)| dt < ∞ is necessary to ensure that the
−∞
integrals (3.13) and (3.15) converge and the Fourier transform and its inverse are
defined.
Remark - There are a variety of definitions of Fourier and inverse-Fourier trans-
form: they all differ by the factor in front of the integral. The definition given here
is symmetric, the factors in√front of the Fourier and inverse-Fourier transform are
identical, both equal to 1/ 2π. The “definition” given in Section 3.1 has a factor
1/2 in front of the Fourier transform, eq. (3.12), and a factor 1/π in front of the
inverse-Fourier transform, eq. (3.9c). The important point is that the product of
the two factors must be equal to 1/2π.
Example - Find the Fourier transform of the function
(
sin(t) −π ≤ t ≤ π,
f (t) = (3.17)
0 Otherwise .
42
FOURIER TRANSFORM PROPERTIES
1 w 1 w
∞ +π
−jωt
F (ω) = √ f (t)e dt = √ sin(t)e−jωt dt
2π −∞ 2π −π
w
+π
1 j
− e−j(1+ω)t dt
j(1−ω)t
=√ − e
2π 2 −π
+π
e−j(1+ω)t
j(1−ω)t
j e
=− √ −
2 2π j(1 − ω) −j(1 + ω) −π

1 1 j(1−ω)π −j(1−ω)π
1 −j(1+ω)π j(1+ω)π

=− √ e −e + e −e
2 2π 1 − ω 1+ω

1 1 jωπ −jωπ
1 jωπ −jωπ

=− √ e −e + e −e
2 2π 1 − ω 1+ω

1 1 1 jωπ −jωπ

=− √ + e −e
2 2π 1−ω 1+ω
1 2 2j
=− √ 2
2j sin(ωπ) = √ sin(ωπ) .
2 2π 1 − ω 2
(ω − 1) 2π
(3.18)
3.3 Properties
We list here some of the basic properties of Fourier transforms. We indicate with
the symbol F the Fourier transform, in the sense that
1 w
∞
F[f (t)] = √ f (t)e−jωt dt = F (ω) . (3.19)
2π −∞
We indicate with the symbol F −1 the inverse Fourier transform.
Linearity
The Fourier transform is a linear operator:
F[αf (t) + βg(t)] = αF[f (t)] + βF[g(t)] . (3.20)
The same apply to the inverse Fourier transform.
43
FOURIER TRANSFORM PROPERTIES
Differentiation
If f (t) is continuous everywhere and its derivative has at most isolated jump dis-
continuities then
df
F = jωF[f (t)] . (3.21)
dt
The generalisation of this formula is
n
d f
F = (jω)n F[f (t)] . (3.22)
dtn
To prove (3.21) we compute the Fourier transform of the derivative of f (t):
1 w df −jωt
∞
df
F =√ e dt . (3.23)
dt 2π −∞ dt
The next step is to integrate by parts and use the fact that |f (t)| → 0 as |t| → ∞
(as required to ensure convergence of the Fourier transform). For simplicity we
consider here the case that the function is differentiable everywhere. In this case,
integration by parts gives
1 w
∞
df 1 −jωt t→+∞
f (t) e−jωt dt =jω F (ω)

F =√ f (t)e t→−∞
+ jω √
dt 2π 2π −∞
(3.24)
Equation (3.22) can be obtained by applying (3.21) n times.
Parseval’s theorem
(not examinable)
Starting from the definition of Fourier transform and using the δ-function it is
possible to prove the following equality, that goes under the name of Parseval’s
theorem1
w∞ w∞
2
|f (t)| dt = |F (ω)|2 dω (3.25)
−∞ −∞
This relation has a clear physical meaning. We have already mentioned when
studying Parseval’s theorem for Fourier series that the moduli of the coefficients
of the series, |an |, |bn | are a measure of the strength of the frequency component
ωn . A similar interpretation can be given for the Fourier transform: if f (t) is some
physical signal, e.g. the amplitude of a sound or electromagnetic wave measured
1
Marc-Antoine Parseval des Chênes, French mathematician (1755-1836)
44
FOURIER TRANSFORM THE RESPONSE FUNCTION
as a function of time, then the modulus square of its Fourier transform, |F (ω)|2 ,
is the energy density in frequency space, i.e. |F (ω)|2 δω is the energy of the signal
in a frequency band of width δω centred around the frequency ω (in appropriate
units). The total energy of the wave is the integral of this energy density, i.e.
w∞
Total energy of wave = |F (ω)|2 dω . (3.26)
−∞
On the other hand, the total energy of a wave of amplitude f (t) is given by inte-
grating |f (t)|2 over time (in appropriate units), i.e.
w∞
Total energy of wave = |f (t)|2 dt . (3.27)
−∞
These two expressions must be equal. In this interpretation, equation (3.25), is just
a statement that the total energy measured must be the same whether we integrate
over time, equation (3.27), or over frequency, equation (3.26).
Shift properties
The Fourier transform of a shifted signal is equal to the Fourier transform of the
original signal, multiplied by a phase factor proportional to the shift:
F[f (t − t0 )] = e−jωt0 F (ω). (3.28)
The converse of (3.28) is
F ejω0 t f (t) = F (ω − ω0 ) .

(3.29)
The signal processing interpretation of this result is that the multiplication of f (t)
by a constant frequency term, exp(jω0 t), is equivalent to shifting the entire fre-
quency spectrum of the signal by an amount ω0 . This is the mathematical basis
of the process of modulation, where a high frequency carrier at frequency ω0 is
modulated by a low frequency signal f (t).
3.4 Fourier transform and the response function

The Fourier transform can be used to study the response function of a linear
system to generic modulations. It is particular useful to analyse its frequency
response. We show this by considering the following example [based on Glyn
James, Modern Engineering Mathematics].
45
A dynamical system y(t) is modelled by the differential equation
d2 y dy du
2
+ 3 + 7y(t) = 3 + 2u(t), (3.30)
dt dt dt
where u(t) is a given signal. We can use the Fourier transform to study the re-
sponse of this system to the modulation due to the function u(t). We let Y (ω) =
F[y(t)] and U (ω) = F[u(t)] denote the Foureir transform of y(t) and u(t) respec-
tively. Taking the Fourier transform of both sides of (3.30) we obtain
(jω)2 Y (ω) + 3jωY (ω) + 7Y (ω) = 3jωU (ω) + 2U (ω) (3.31a)

2 + 3jω
=⇒ Y (ω) = U (ω) (3.31b)
7 − ω 2 + 3jω
= G(ω)U (ω), (3.31c)
where
2 + 3jω
G(ω) = . (3.32)
7 − ω 2 + 3jω
The function G(ω) is the transfer function of the system: it determines how the
system responds to the various frequency components of the forcing term u(t). In
particular, we have
|Y (ω)| = |G(ω)| |U (ω)| . (3.33)
In other words, the strength of the frequency component ω of the output is equal
to that of the input multiplied by the modulus of the transfer function. In the case
of the function G(ω) defined by (3.32) we have
r
9 ω2 + 4
|G(ω)| = (3.34)
49 − 5 ω 2 + ω 4
This quantity is plotted in Figure 3.2: from the graph we can deduce that the
system acts like a band pass filter with maximal response at
√
q
1
ω= −4 + 7 85 ' 2.59 . (3.35)
3
An example of the effect of the response function on an input signal is shown in
Figure 3.3. In this case the input signal consists in a Gaussian pulse in time (solid
line, Figure 3.3 right). Its Fourier power spectrum is also Gaussian (solid line,
Figure 3.3 left). The power spectrum of the output (dashed line, Figure 3.3 left)
is the product Y (ω) = G(ω)U (ω) where G(ω) is given by equation (3.32). The
time dependent output y(t) (dashed line, Figure 3.3 right) is the inverse Fourier
transform of Y (ω). A more detailed analysis of the response of a linear system can
46
Plot of the response function as a function of the frequency
1
|G(ω)|
0.5
0
0 5 10 15 20
ω
Figure 3.2: Plot of the modulus of the response function, equa-
tion (3.34), as a function of the frequency.
Frequency domain response Time domain response

1 0.15
Input: U(ω) Input: u(t)
Power spectrum − |U(ω)|2 and |Y(ω)|2
Output: Y(ω) Output: y(t)

0.8
0.1
u(t) and y(t)
0.6
0.05
0.4
0
0.2
0 −0.05
−60 −40 −20 0 20 40 60 0 1 2 3 4
Frequency (ω) Time [× π]
Figure 3.3: Input and output signal in the frequency (left) and
time (right) domain for equation (3.30). The relation between
the Fourier series of the input U (ω) and output Y (ω) signal is
Y (ω) = G(ω)U (ω), with the response function G(ω) given by
equation (3.32).
47
FOURIER TRANSFORM SUMMARY
be obtained by studying the singularities of the transfer function in the complex

plane.
In the next chapter we will look at a related integral transform called the
Laplace Transform. This is defined by
Z ∞
L[f (t)](s) = f (t)e−st dt
0
This is up to a factor just the Fourier transform where ω = js. It is used more
widely then the Fourier transform for solving differential equations as it is simpler
to use when one can claculate the inverse transform. However unlike the Fourier
transform there is no simple formula for calculating the inverse of a Lapalce trans-
form - one either has to look it up in a table (which does not always work) or else
regard it as a “complex Fourier transform” in which case one can calculate the
inverse using complex integration.
3.5 Summary
The Fourier transform is important because it we can use it to transform an ordi-
nary differential equation into an algebraic equation (or a PDE into an ODE). We
can then solve the algebraic equation and find the solution to the original ODE by
applying the inversse Fourier transform.
• Definition
The Fourier transform of a function f (t) is defined by
1 w
∞
F[f (t)](ω) = √ f (t)e−jωt dt = F (ω) .
2π −∞
• Inverse Fourier transform

The inverse Fourier transform of a function F (ω) is defined by
1 w
∞
−1
F [F (ω)](t) = √ F (ω)ejωt dω = f (t)
2π −∞
The key properties of the Fourier transform are the following:
• Linearity
F[αf (t) + βg(t)] = αF[f (t)] + βF[g(t)] .
48
FOURIER TRANSFORM SUMMARY
• Derivatives The Fourier transform of the derivative of a function f 0 (t) is

given by
F[f 0 (t)] = jωF[f (t)] .
More generally the Fourier transform of the n-th derivative of a function
f (n) (t) is given by
F[f (n) (t)] = (jω)n F[f (t)] .
• First Shift Theorem
F[f (t − t0 )] = e−jωt0 F (ω).
• Second Shift Theorem
F ejω0 t f (t) = F (ω − ω0 ) .

49
Chapter 4
Laplace transforms
4.1 Introduction
A standard method for solving differential equations is to tranform the problem
into an equivalent, simpler problem. For a broad range of differential equations
(especially linear differential equations) the integral transforms
wβ
ỹ(s) = K(s, x)y(x) dx (4.1)
α
convert the problem from a complex differential equation for y(x) to a (hopefully)
simpler equation for ỹ(s). Here the kernel K(s, x) is a given function which
together with the real integration limits α, β defines the transform.
The broad outline of the method is the same for all types of integral transform:
1. apply the transform to the original equation(s) to get a simpler problem for
ỹ;
2. solve the simpler problem for ỹ;
3. invert the transform to get y.
In certain cases it may only be formally possible to invert the transform; that is,
we may not be able to compute the solution in closed form, but only as an integral.
Here we will only study one type of integral transform, the Laplace transform.
50
LAPLACE TRANSFORMS DEFINITION AND EXISTENCE
4.2 Definition and existence

For a suitably well-behaved function f (x) we define the Laplace transform of
f (x) to be
w∞
f˜(s) = f (x)e−sx dx (4.2a)
0
= L [f (x)] . (4.2b)
We note that a wide variety of notation is used; the Laplace transform is also often
denoted f¯(s) or F (s).
4.2.1 Existence
Strictly the definition is not complete. For unrestricted functions f it is not at
all obvious that the integral in the definition (4.2) converges for any given value
of s. In fact, it is simple to check (using the strict definition of the unbounded
integral) that the integral need not converge for some values of s for certain simple
functions.
We look at ecx where c is some real nonzero constant. We therefore have
w∞ wA
cx
e dx = lim ecx dx (4.3a)
A→∞
0 0
A
ecx

= lim (4.3b)
A→∞ c 0
1 cA
= lim e −1 . (4.3c)
A→∞ c
We immediately see that the limit is finite, and hence the integral converges, if
c < 0, and that the integral diverges if c > 0. The case c = 0 is not covered by
this, but it is trivial to check that the integral diverges in this case as well.
This immediately motivates:
Theorem 4.2.1. If f (x) is a piecewise continuous function on 0 ≤ x ≤ A for
all A > 0, and |f (x)| ≤ Keax for all x ≥ M , where K, a, M are some real
constants and K, M are both positive, then the Laplace transform L [f (x)] = f˜(s)
as defined by (4.2) exists for s > a.
The proof is given, for example, in Boyce & DiPrima. We note that if the
conditions of the theorem hold then
lim f (x)e−sx = 0. (4.4)
x→∞
51
LAPLACE TRANSFORMS GENERAL PROPERTIES
4.3 General properties

There are a number of key basic properties that can easily be checked, and which
are essential in calculations involving Laplace transforms.
1. If L [f (x)] = f˜(s) and C is a constant, then
L [Cf (x)] = C f˜(s). (4.5)
2. If L [f1 (x)] = f˜1 (s) and L [f2 (x)] = f˜2 (s) then
L [f1 (x) + f2 (x)] = f˜1 (s) + f˜2 (s). (4.6)
Both of these properties follow immediately from the linearity of integra-

tion.
3. If L [f (x)] = f˜(s) then
L e−ax f (x) = f˜(s + a).

(4.7)
This is the first shift theorem and is a simple application of the defini-
tion (4.2). It is extremely useful when inverting transforms.
4. If L [f (x)] = f˜(s) then
df˜
L [xf (x)] = − (s). (4.8)
ds
This follows using the definition (4.2) and integration by parts. Less useful
in practice, it illustrates the difficulties that arise when transforming nonlin-
ear equations.
5. If L [f (x, a)] = f˜(s, a) then
∂ f˜

∂f
L = . (4.9)
∂a ∂a
This is an obvious consequence of the definition (4.2), and integration by

parts gives the additional result
w
"x #
1
L f (u) du = L [f (x)] . (4.10)
0
s
52
LAPLACE TRANSFORMS GENERAL PROPERTIES
4.3.1 Transforming derivatives

As our key application of Laplace transforms is going to be to differential equa-
tions, we are most interested in how a Laplace transform affects a derivative.
We directly apply the definition to find
w∞
df df −sx
L = e dx (4.11a)
dx 0
dx
w∞
−sx ∞
+ s f (x)e−sx dx.

= f (x)e 0
(4.11b)
0
We then recall from section 4.2.1 that in order for the Laplace transform to exist
we have
lim f (x)e−sx = 0. (4.4)
x→∞
We therefore get
df
L = sf˜(s) − f (0). (4.12)
dx
This allows us to transform first derivatives of f into algebraic expressions involv-
ing f˜ and the initial data f (0).
We can repeat this exercise for the second derivative, finding that
2
df d df
L 2
=L (4.13a)
dx dx dx

df
= sL − f 0 (0) (4.13b)
dx
h i
= s sf˜(s) − f (0) − f 0 (0) (4.13c)
= s2 f˜(s) − sf (0) − f 0 (0). (4.13d)
We note that in taking these steps we have assumed that the Laplace transform for
f 0 exists, and hence assumed that
lim f 0 (x)e−sx = 0. (4.14)
x→∞
Clearly we can extend this result to higher order by repeating these methods,
finding
n n−1
d f (x) n˜
X
L = s f (s) − sn−k−1 f (k) (0), (4.15)
dxn k=0
where we require
lim f (k) (x)e−sx = 0, k = 0, . . . , n. (4.16)
x→∞
53
LAPLACE TRANSFORMS INVERSE TRANSFORM
4.4 Inverse Transform

We first ask the key question of whether two different functions can have the same
Laplace transform. If the two functions f1 (x), f2 (x) both transform to f˜(s) then
we have that
w∞ w∞
−sx ˜
e f1 (x) dx = f (s) = e−sx f2 (x) dx, (4.17)
0 0
and hence
w∞
0= e−sx (f1 (x) − f2 (x)) dx (4.18a)
0
⇒ 0 = f1 (x) − f2 (x) (4.18b)
⇒ f1 (x) = f2 (x). (4.18c)
A truly rigorous proof of this requires various restrictions on the continuity of

f1 , f2 and on their behaviour as x → ∞, but continuity plus the restrictions given
in section 4.2.1 are sufficient.
Once we have that there is a one-to-one correspondence between a function
and its Laplace transform, we see that the simplest method to invert a Laplace
transform f˜ is simply to find or spot the function f that transforms to it. This is
usually done using tables given in books or in the Course Summary Notes. With
the addition of certain results, such as those given in section 4.3 and those to come
later, most simple transforms can rapidly be inverted.
There does exist a general formula giving the inverse transform:
1 w ˜
γ+j∞
f (x) = f (s)esx dx . (4.19)

2πj
γ−j∞
| {z }
Contour integral
This Bromwich inversion formula is beyond the scope of this course.
4.5 Special functions

4.5.1 Heaviside function
The Heaviside step function
(
0, x<0
H(x) = (4.20)
1, x>0
54
LAPLACE TRANSFORMS SPECIAL FUNCTIONS
is discontinuous at x = 0 and is useful in modelling problems where an effect

suddenly “switches on”. It also has a number of useful mathematical properties
when combined with Laplace transform.
We can use the definition to show
w∞
L [H(x − a)] = H(x − a)e−sx dx (4.21a)
0
wa w∞
−sx
= 0×e dx + 1 × e−sx dx (4.21b)
0 ∞ a
e−sx
= − (4.21c)
s a
−sa
e
= . (4.21d)
s
But more usefully we can multiply the Heaviside function by any function
f (x) that has a Laplace transform f˜(s), to find
w∞
L [H(x − a)f (x − a)] = H(x − a)f (x − a)e−sx dx (4.22a)
0
w∞
= f (x − a)e−sx dx (4.22b)
x=a
τw
=∞
= f (τ )e−s(τ +a) dτ (4.22c)
τ =0
using the change of variables τ = x − a, to get

w∞
−sa
=e f (τ )e−sτ dτ (4.22d)
0
=e −sa
f˜(s). (4.22e)
This is the second shift theorem, which should be compared to the first shift
theorem given in section 4.3.
These theorems are most useful to invert Laplace transforms, as they show that
h i
L−1 f˜(s + a) = e−sa f (x), (4.23a)
h i
L−1 e−sa f˜(s) = H(x − a)f (x − a). (4.23b)
55
LAPLACE TRANSFORMS SPECIAL FUNCTIONS
4.5.2 The Dirac δ-function

The Dirac δ-“function” is not a genuine function at all. It is used in modelling
instantaneous impulses, and as with the Heaviside function has a number of useful
mathematical properties as well. There are a number of ways of defining the δ-
function, but its essential properties are that
w∞
δ(x)f (x) dx = f (0), (4.24a)
−∞
wb
δ(x)f (x) dx = 0 0 6∈ [a, b]. (4.24b)
a
In particular, we note that to be strict the δ-function only makes sense when inte-
grated.
One way of thinking about the δ-function is that it takes the values
(
0, x 6= 0,
δ(x) = (4.25)
∞, x = 0.
This is obviously an unsatisfactory definition at x = 0. It is more usual to rely on
definitions based on limits, such as
δ(x) = lim D(x, ) (4.26a)
→0
where
(
1/ − 2 ≤ x ≤
2
D(x, ) = , (4.26b)
0 otherwise,
which gives

w∞ w2 1
D(x, ) dx = dx (4.26c)
−∞

2
=1 ∀. (4.26d)
Note that lim→0 D(0, ) = ∞.
We see that the δ-function, when multiplied by f and integrated, “picks out”
the value of f at the value of x where the argument of the δ-function vanishes.
Thus the Laplace transform has the form
w∞
L [δ(x − a)] = δ(x − a)e−sx dx (4.27a)
0
−as
=e , (4.27b)
56
LAPLACE TRANSFORMS CONVOLUTION
and in particular L [δ(x)] = 1.

There is a sense in which the δ function is the derivative of the Heaviside
function. This is certainly intuitively appealing and for certain arguments can be
made rigorous.
4.6 Convolution
(Not examined)
One point that we have not emphasized but which should be clear is that
L [f (x) × g(x)] 6= L [f (x)] × L [g(x)] . (4.28)
It should also be clear that

h i h i
L−1 f˜(s) × g̃(s) 6= L−1 f˜(s) × L−1 [g̃(s)] . (4.29)
As in practical examples we will frequently find the expression h̃ that we wish to

invert is a product f˜g̃, it is useful to find an expression for this case.
The method relies on the convolution theorem, which can be expressed
w
"x #
L f (u)g(x − u) du = f˜(s)g̃(s). (4.30)
0
In particular this implies

h i wx
L −1 ˜
f (s)g̃(s) = f (u)g(x − u) du. (4.31)
0
We define the convolution operation f ∗ g as

wx
(f ∗ g)(x) = f (u)g(x − u) du, (4.32)
0
and say that the Laplace transform of the convolution is the product of the Laplace
transforms.
Convolution Theorem.
L [(f ∗ g)(x)] = f˜(s)g̃(s). (4.30)
57
LAPLACE TRANSFORMS APPLICATIONS
Proof. We directly apply the definition of the Laplace transform to find

w
"x #
L [(f ∗ g)(x)] = L f (u)g(x − u) du (4.33a)
0
w∞ wx
!
= e−sx f (u)g(x − u) du dx (4.33b)
0 0
and, by interchanging the integration order, find

w∞ w∞
!
−sx
= f (u) e g(x − u) dx du (4.33c)
0 u
w∞ w∞
!
= f (u) e−s(x−u) e−su g(x − u) dx du (4.33d)
0 u
and, by introducing z = x − u
w∞ w∞
!
= f (u)e−su e−sz g(z) dz du (4.33e)
0 z=0
w∞ w∞
! !
= f (u)e−su du e−sz g(z) dz (4.33f)
0 0
= f˜(s)g̃(s) (4.33g)
as required.
Convolution is often used to write the final solution y(x) in terms of a formal
integral solution.
4.7 Applications
4.7.1 Solving ODEs
Laplace transforms may be used to solve ODEs. To do this we follow the follow-
ing steps.
1. We apply the Laplace transform to the equation we want to solve.
Using linearity this amounts to rewriting the equation in terms of the Laplace
transform of the function we want to obtain.
58
2. We use the formula for the Laplace transform of derivatives and insert the
initial condition.
3. This leads to an algebraic equation for the Laplace transform of the solution
we want to find, which we solve.
4. We determine the inverse Laplace transform to find the solution.
Example – Let us solve the following differential equation using Laplace

transform,
y 00 + y = 0, y(0) = 0, y 0 (0) = 1. (4.34)
The first step is to apply the Laplace transform to the equation:
0 = L [y 00 + y] = L [y 00 ] + L [y] (4.35)
In step 2 we use (4.13),
L [y 00 ] = s2 ỹ(s) − sy(0) − y 0 (0) = s2 ỹ(s) − 1 (4.36)
where in the last equality we used the boundary conditions. Inserting (4.36) in
(4.35) (step 3) yields an algebraic equation for ỹ(s), which we solve:
1
0 = s2 ỹ(s) − 1 + ỹ(s) ⇒ ỹ(s) = (4.37)
1 + s2
Finally (step 4), we need to determine the inverse Laplace transform of ỹ(s).
Looking at the Table of Laplace transforms in the Formula Sheet we find

−1 1
y(t) = L (t) = sin t (4.38)
1 + s2
4.7.2 Systems of ODEs

We know that we can transform a single linear ODE such as
dy
− ay = f (x) (4.39)
dx
to the algebraic equation in terms of its Laplace transform ỹ(s),
(s − a) ỹ(s) = f˜(s) + y(0). (4.40)
59
We could similarly transform the system of equations

dy
− Ay = f (x), (4.41)
dx
where A is a matrix compatible with y, f , to the system of algebraic equations
(sI − A) ỹ(s) = f˜(s) + y(0). (4.42)
Given the appropriate initial conditions for y we could then solve this vector equa-
tion for fixed s, as it is a standard linear algebra problem. In particular we note
that the solution
−1 ˜
ỹ(s) = (sI − A) f (s) + y(0) (4.43)
only exists when sI − A is non-singular, which is when s is not an eigenvalue of
A.
4.7.3 PDEs
As Laplace transforms are integral transforms they affect only one variable at a
time. We can therefore use a Laplace transform to convert a PDE to an ODE. For
example, if we consider the wave equation in spherical symmetry
1 ∂ 2y 1 ∂ 2 (ry)
(r, t) = (r, t), (4.44)
c2 ∂t2 r ∂r2
we can perform a Laplace transform with respect to t to find
1 ∂ 2 (rỹ)

1 2 ∂y
s ỹ(s, r) − sy(r, 0) − (r, 0) = (s, r). (4.45)
c2 ∂t r ∂r2
Once we have initial conditions for y and its time derivative then we will have
an ODE for ỹ. In particular if we make the simplest choice that y and its time
derivative vanish at t = 0 then we are left with the ODE
d2 (rỹ) s 2
− (rỹ) = 0, (4.46)
dr2 c
which obviously has the solution
sr sr
rỹ = De c + Ee− c . (4.47)
As we are assuming that c > 0, and r, s > 0 by definition, in order for the solution
to be regular as r → ∞ we must have D ≡ 0.
60
In order to compute the final solution we cannot directly invert the result that
follows from equation (4.47),
sr
e− c
ỹ(s) = E , (4.48)
r
as the “constant” E that follows from solving equation (4.46) may be E ≡ E(s).
Instead we have to use an appropriate boundary condition in order to fix E. If we
have a boundary condition such as

∂y
= g(t) (4.49)
∂r r=a
then we can Laplace transform this, finding

dỹ
(r, s) = g̃(s) (4.50a)
dr r=a
sa
e − c a2 s

⇒ −E 2 + 1 = g̃(s) (4.50b)
a c
a2 cg̃(s) sa
⇒ E=− 2 ec. (4.50c)
a s+c
We can then form our Laplace transformed solution,
a2 cg̃(s) s(a−r)
ỹ(s) = − e c , (4.51)
r (a2 s + c)
which can (in principle) be inverted to find the solution.
61
Chapter 5
Wave equation
5.1 Classification of PDEs

We consider second order, linear PDEs of the form
∂ 2u ∂ 2u ∂ 2u
a(x, y) + 2b(x, y) + c(x, y)
∂x2 ∂x∂y ∂y 2
(5.1)
∂u ∂u
+ d(x, y) + e(x, y) + f (x, y)u = 0.
∂x ∂y
We assume that c 6= 0 and perform the general change of variables
ξ = x + βy, η = x + δy, δ 6= β. (5.2)
This transforms the general form of equation (5.1) to

∂ 2u ∂ 2u
a + 2bβ + cβ 2 + 2 (a + b(δ + β) + cβδ) +
∂ξ 2 ∂ξ∂η
(5.3)
2

2 ∂ u ∂u ∂u
a + 2bδ + cδ = G u, , , ξ, η, F .
∂η 2 ∂ξ ∂η
We now use our freedom to choose β, δ to be roots of the quadratic equation
cλ2 + 2bλ + a = 0. (5.4)
This implies that the coefficients of uξξ and uηη in equation (5.3) vanish, so that it
simplifies to the canonical form
4 ∂ 2u
ac − b2 = G. (5.5)
c ∂ξ∂η
62
WAVE EQUATION SMALL AMPLITUDE VIBRATIONS
Figure 5.1: An elastic string is stretched between two fixed points

at x = 0, L. The tension is constant (T ) and the density (ρ) is
constant. The wave equation describes the motion of the string
through the behaviour of the (“small”) displacement y(x, t).
As we will see with the wave equation, when a PDE can be reduced to this form it
has simple solutions (at least locally). Clearly the type of the solution will depend
strongly on the qualitative form of the values β, δ used to change to the canonical
form. In particular if
• b2 > ac: the roots of equation (5.4), β and δ are real and distinct. The
equation is called hyperbolic.
• b2 = ac: the roots of equation (5.4), β and δ are real and repeated. The
equation is called parabolic.
• b2 < ac: the roots of equation (5.4), β and δ are complex conjugates. The
equation is called elliptic.
5.2 Small amplitude vibrations of a stretched string

Consider a purely elastic string having a uniform density ρ per unit length, stretched
at tension T , constant, between the fixed points x = 0 and x = L, subject to a
small deflection y(x, t).
Because of the fixed points, the boundary conditions in space are
y(0, t) = 0 = y(L, t). (5.6)
Given an initial shape and velocity, i.e. prescribed initial conditions (or boundary
conditions in time)
∂y
y(x, 0) = p(x), (x, 0) = q(x), (5.7)
∂t
what is the shape at subsequent times?
63
Figure 5.2: A small section of the string and the forces acting on
it. By examining the forces over this small section and insisting
that the string not break we find the wave equation holds for small
displacements.
5.2.1 Equations of motion

We consider motions with small slopes

∂y
1 ∀x, t, (5.8)
∂x
as for large slopes the methods used below will not work and the assumptions
made above do not hold. In the case of small slopes then we can use the small
angle expansion of the trigonometric functions to see that

∂y
α≈ . (5.9)
∂x x
We resolve the forces on the string in the x-direction to find
T1 cos α ≈ T2 cos β as δx → 0. (5.10)
But we can again use the small angle expansion to see that this implies that T1 ≈
T2 , i.e. the tension T is constant along the length of the string (as assumed above).
64
Resolving the forces in the y-direction gives

∂ 2y
ρ δx = T2 sin β − T1 sin α (5.11a)
∂t2
∂y ∂y
≈T − (5.11b)
∂x x+δx ∂x x
∂ 2y
≈T δx. (5.11c)
∂x2
This immediately gives
∂ 2y 2
2∂ y
= c , (5.12)
∂t2 ∂x2
the 1-dimensional wave equation, where the wave speed c is given by
T
c2 = . (5.13)
ρ
5.2.2 Energy of vibrating string

The total kinetic energy, that is the energy due to the motion of the string, is given
by
wL 1 ∂y 2
KE = ρ dx. (5.14)
0
2 ∂t
The elastic potential energy, that is the energy due to the tension caused by the
change in the length of the string, is given by
s 
wL 2
∂y
PE = T  1 + − 1 dx (5.15a)
0
∂x
wL 1 ∂y 2
≈T dx (5.15b)
0
2 ∂x
as the slopes, and hence the derivatives, are small (equation (5.8)).
Hence the total energy E is
wL 1 ∂y 2
" 2 #
∂y
E= ρ + c2 dx. (5.16)
0
2 ∂t ∂x
As there is no friction in the system, the total energy must be constant in time.
It is a straightforward exercise (check!) to prove from the wave equation and the
boundary conditions that the energy is constant.
65
WAVE EQUATION D’ALEMBERT’S SOLUTION
5.3 D’Alembert’s solution

We now look at the wave equation
∂ 2y 2
2∂ y
= c (5.17)
∂t2 ∂x2
in the absence of boundaries.
We introduce a new coordinate system, sometimes called characteristic coor-
dinates:
ξ = x + ct, (5.18a)
η = x − ct. (5.18b)
The motivation for these coordinates should be clear from section 5.1, but for the
moment we will just rewrite the wave equation (5.17) in terms of y = y(ξ, η) and
see where it takes us.
Firstly we note that the chain rule gives us
∂ ∂ ∂
=1· +1· , (5.19a)
∂x ∂ξ ∂η
∂ ∂ ∂
=c· −c· . (5.19b)
∂t ∂ξ ∂η
Thus the wave equation (5.17) becomes
2 2
∂ ∂ ∂ ∂
+ y= − y (5.20)
∂ξ ∂η ∂ξ ∂η
which gives
∂ 2y
= 0. (5.21)
∂ξ∂η
This can immediately be solved one coordinate at a time to give, e.g.,
∂y
= h(ξ) (5.22)
∂ξ
and hence
y = f (ξ) + g(η) (5.23)
where the functions f, g are arbitrary. Therefore the general solution of the wave
equation (5.17), ignoring initial and boundary conditions, can be written as
y(x, t) = f (x + ct) + g(x − ct). (5.24)
66
5.3.1 Travelling waves

The solution given by y = g(x − ct) is a travelling or progressive wave. It is
given this name as on the lines in the (x, t) plane where x − ct is constant, called
characteristics, the solution is constant. Therefore the wave retains its shape and
moves with constant speed c. Characteristics can be similarly defined for solutions
y = f (x + ct) which travel in the opposite direction. This property is useful when
the string is considered to be infinite so that the effects of the boundary conditions,
which are neglected here, are unimportant. However, it is generally true that for
the wave equation information propagates at speed c.
There is an even more special case solution, called a harmonic travelling wave
y = Aejω(t−x/c) . (5.25)
Here A is a complex constant, and it is understood that we always take the real
part of the right hand side to get the solution. If we write the complex constant A
as A = |A|ejφ then we see that
y = |A| cos (|{z}

ω (t − x/c) + φ ) . (5.26)
|{z} |{z}
amplitude frequency phase
The key quantity is typically the frequency of the wave, ω, which implicitly defines
its wavelength 2πc
ω
.
5.3.2 Wave reflection and transmission

Wave reflection and transmission arises where the “natural” properties, such as the
density of the string ρ, change discontinuously.
p This will lead to a discontinuous
change in the wave speed c = T /ρ. Typically a travelling wave is incident on
(i.e., hits) the discontinuity and is partially reflected and transmitted. We want to
predict the amount of reflection and transmission.
For simplicity we assume that the discontinuity arises in the density ρ only,
and we choose our coordinates such that it occurs at x = 0. Thus to the left,
x < 0, we have ρ = ρ− leading to the wave speed c− , whilst to the right, x > 0,
we have ρ = ρ+ leading to the wave speed c+ . We assume that the incident wave
is harmonic, i.e.,
y = Aejω(t−x/c− ) , (5.27)
and that the transmitted and reflected waves are given by travelling waves f, g
respectively, i.e.
x x
y →f t− +g t+ . (5.28)
c+ c−
67
From our assumption that the string does not break we have that y is contin-
uous at x = 0. This tells us that the sum of the incident wave and the reflected
wave is the transmitted wave at x = 0 for all times t.
We can also use the force balance at x = 0 and that the tension in the string is
constant to give
∂y ∂y
T =T , (5.29)
∂x x=0− ∂x x=0+
∂y
which implies that ∂x is continuous at x = 0 as well.
In order to get a solution we have to assume a specific form for the transmitted
and reflected waves. We look for harmonic travelling wave solutions for these as
well, where
f = Dejω(t−x/c+ ) , (5.30a)
g = Bejω(t+x/c− ) . (5.30b)
The continuity conditions at x = 0 give
y continuous: A+B =D (5.31a)
∂y jωA jωB jωD
continuous: − + =− . (5.31b)
∂x c− c− c+
We can solve these equations to get the unknown coefficients B, D in terms of the
incident wave defined by A to get
c+ − c−
B= A, (5.32a)
c+ + c−
2c+
D= A. (5.32b)
c+ + c−
Given that the problem is linear we could have taken A to be 1 from the outset,
which would give the amplitudes and phases of the reflected (B) and transmitted
(D) waves relative to the incident wave.
Note three special cases:
1. c+ = c− : In this case B = 0 and D = A, which is exactly what we would
expect. The material properties do not change and information is perfectly
transmitted with no reflection.
2. ρ+ /ρ− → ∞, or equivalently c+ /c− → 0: In this case D = 0 and B = −A.
This corresponds to no transmission or perfect reflection; the reflected wave
has the same amplitude but is out of phase by π. This is “reflection from a
solid wall”.
3. ρ+ /ρ− → 0: In this case D = 2A and B = A. This corresponds to the
string becoming “massless”.
68
WAVE EQUATION SEPARATION OF VARIABLES
5.4 Separation of variables

5.4.1 Ansatz
We wish to solve the 1-dimensional wave equation
∂ 2y 2
2∂ y
= c , 0 < x < L, t > 0, (5.33)
∂t2 ∂x2
with the boundary conditions
y(0, t) = 0 = y(L, t), t > 0, (5.34)
and initial conditions

∂y
y(x, 0) = p(x), (x, 0) = q(x) (5.35)
∂t
for given p, q.
We will firstly seek a solution in the very special form
y(x, t) = X(x)T (t), (5.36)
i.e. with variables x and t separated. We substitute this ansatz into the wave equa-
tion (5.33) finding
X T̈ = c2 X 00 T (5.37a)
or, on dividing by y = XT ,
1 T̈ X 00
= . (5.37b)
c2 T X
Now we note that the left hand side of this equation depends solely on T and its
derivatives with respect to time; that is, the left hand side is solely a function of t.
However, the right hand side depends solely on X and its derivatives with respect
to space; that is, the right hand side is solely a function of x.
How is it possible for a function solely of t to equal a function solely of x?
If both functions are trivial, i.e. equal to zero, then it is obviously true. If both
functions are constant then this gives us a constraint; the two functions must be
equal to the same constant. However, if either function varies then this equation
cannot hold, as (say) the left hand side will vary as t varies whilst the right hand
side stays constant. Therefore this variables separable form immediately implies
that both sides of the equation must be separately constant. In general, if the
equation can be written as the sum of functions of different variables then each
function must be separately constant.
69
In this case we have, from equation (5.37), that both sides must be equal to the
same constant. Let us call it λ. Then we can re-arrange our variables separable
form of equation (5.37) to give
X 00 − λX = 0, (5.38a)
T̈ − λT = 0. (5.38b)
So in this very special case we have reduced the partial DE defining the wave
equation, equation (5.33), into a pair of ordinary DEs which are coupled through
the separation constant λ.
5.4.2 Solving the separated form: spatial dependence

In order to solve the wave equation we have to solve the ODEs given by equa-
tion (5.38). Let us consider equation (5.38a) for the spatial behaviour. At present
it is not completely specified, because we have not given any boundary conditions.
We do not have any boundary conditions for X as it has no physical meaning. In-
stead we have to find boundary conditions for X that are compatible with the
boundary conditions for y given by equation (5.34). As these boundary condi-
tions force y to vanish at x = 0, L we could impose that X = 0 at x = 0, L. If we
consider what would happen if we assumed that X 6= 0 at one of these boundaries
we see that the only way for y = XT to vanish would be for T to vanish for all
time t, which would give the trivial solution y ≡ 0. Hence we need to solve
X 00 − λX = 0, X(0) = 0, X(L) = 0. (5.39)
This is a second order ODE with one additional unknown, the value of the
separation constant λ. So we do not expect to be able to completely solve the
problem with the information provided, as we only have two boundary conditions.
We do need to ensure that we find the value of the separation constant λ, however,
as it is this that couples the two ODEs that we are trying to solve in equation (5.38).
We find that it is best to consider the solution of the ODE given by equa-
tion (5.39) in terms of three separate cases.
1. λ > 0. In this case the general solution to the (constant coefficient, linear,
second order, homogeneous) ODE is
√ √
X(x) = Ae λx
+ Be− λx
. (5.40)
We need this to be compatible with the boundary conditions. By setting
70
x = 0, L we immediately get
x=0: A + B = 0, (5.41a)
√
2 λL
x=L: Ae + B = 0, (5.41b)
⇒ A = B = 0. (5.41c)
Hence we only have the solution X ≡ 0 which implies that y ≡ 0, the trivial
solution. This is of no use.
2. λ = 0. In this case the general solution to the (constant coefficient, linear,

X(x) = A + Bx. (5.42)
We need this to be compatible with the boundary conditions. By setting

x = 0, L we immediately get
x=0: A = 0, (5.43a)
x=L: A + BL = 0, (5.43b)
⇒ A = B = 0. (5.43c)
Hence we only have the solution X ≡ 0 which implies that y ≡ 0, the trivial
solution. This is of no use.
3. λ < 0. In this case the general solution to the (constant coefficient, linear,
√ √
X(x) = A cos −λx + B sin −λx . (5.44)
We need this to be compatible with the boundary conditions. If we impose

the boundary condition at x = 0 then we immediately see that A = 0. If we
impose the boundary condition at x = L then we see that
√
0 = B sin −λL . (5.45)
If this were to imply that B = 0 then once again we would only have the
trivial solution y = 0. However, there is another possibility. It could be
that the sin term vanishes instead. We know that sin vanishes whenever its
argument is an integer multiple of π. In other words, there is a non-trivial
solution to this equation if
√
−λL = nπ, n ∈ Z. (5.46)
71
Hence, by looking at all the possible cases, we see that the only non-trivial solu-
tions to the ODE for the spatial dependence, equation (5.39), are
nπx nπ 2
Xn (x) = An sin , λn = − , n ∈ Z. (5.47)
L L
Here An is an undetermined constant.
5.4.3 Solving the separated form: time dependence

We still have not found solutions to the full wave equation, as we have only solved
for the spatial dependence in our separated problem given by equation (5.38).
However, we have found the separation constant λ; it can take infinitely many dif-
ferent values, all discrete, depending on n as given by equation (5.47). So we can
use this value for the separation constant when solving for the time dependence of
the solution.
Substituting the value of λ into the ODE governing the time dependence, equa-
tion (5.38b), we find nπ 2
T̈ + T = 0. (5.48)
L
We can immediately write down the solution of this (linear, second order, constant
coefficient, homogeneous) ODE as the coefficients have definite sign, to get

nπct nπct
Tn (t) = Cn cos + Dn sin . (5.49)
L L
We have not yet used any conditions that could determine the value of the constant
coefficients Cn , Dn .
5.4.4 Normal modes

We now know, up to various free constants, the behaviour of the space (equa-
tion (5.47)) and time (equation (5.49)) dependence of the separated solution. There-
fore we can multiply the results together to get a solution (or set of solutions) to
the wave equation,

nπct nπct nπx
yn (x, t) = Cn cos + Dn sin sin . (5.50)
L L L
Here we have slightly abused notation by absorbing the free An coefficient in

the spatial dependence in equation (5.47) into the Cn , Dn coefficients in the time
dependence, but without changing the names of the coefficients.
72
The set of solutions given by equation (5.50) satisfy the wave equation itself,
and the boundary conditions, but not the initial conditions. They are a special set
of solutions that have the property of retaining a constant shape of the form of
a sin function with only the amplitude varying in time. These are called normal
modes.
For reasons that are related to Sturm-Liouville theory, the spatial shape sin nπx

L
are also called the eigenfunctions, whilst the separation constants λn are called the
eigenvalues.
5.4.5 Superposition
We still do not have the solution to the problem that we really want, which is a
complete solution of the wave equation (5.33) satisfying the boundary conditions
of equation (5.34) and the initial conditions of equation (5.35). The normal mode
solutions given by equation (5.50) solve the wave equation and the boundary con-
ditions, but say nothing about the initial conditions.
However, we can exploit the fact that the wave equation and the boundary
conditions are linear. Therefore if we have two solutions y1 , y2 that satisfy the
wave equation and the boundary conditions then any linear combination
y = A1 y1 + A2 y2 (5.51)
will also satisfy both the wave equation and the boundary conditions. Therefore
we can take a general solution of the wave equation to be a linear combination of
the normal mode solutions given by equation (5.50):
∞
X nπct nπct nπx
y(x, t) = Cn cos + Dn sin sin . (5.52)
n=1
L L L
This is a superposition of the normal mode solutions. This general solution satis-
fies the wave equation and the boundary conditions, but has sufficient freedom in
the choice of the Cn , Dn coefficients to also satisfy the initial condition given by
equation (5.35).
To actually use this in practise we simply set t = 0 in our general solution
which, when matched against the initial data of equation (5.35), gives
∞
X nπx
p(x) = Cn sin . (5.53)
n=1
L
This is just a Fourier sine series for the Cn coefficients. Similarly, by taking the
time derivative of our general solution given by equation (5.52) and then setting
73
t = 0, we find
∞
X nπc nπx
q(x) = Dn sin . (5.54)
n=1
L L
This is again just a Fourier sine series for the Dn coefficients.
Using the standard Euler formulas for the coefficients of Fourier series we can
explicitly write
2w
L nπx
Cn = p(x) sin dx, (5.55a)
L0 L
2 w
L nπx
Dn = q(x) sin dx. (5.55b)
nπc 0 L
5.4.6 Additional notes

Travelling waves
If we consider the case where the initial “velocity” of the wave as given by ∂y
∂t
vanishes, i.e. q ≡ 0 ≡ Dn , then our solution reduces to
∞
X nπct nπx
y(x, t) = Cn cos sin (5.56a)
n=1
L L
and using the compound angle formulas for sin we write this as
∞
1X nπ(x + ct) nπ(x − ct)
= Cn sin + sin (5.56b)
2 n=1 L L
1
= [p(x + ct) + p(x − ct)] , (5.56c)
2
where in the last step we have noted that
∞
X nπx
p(x) = Cn sin . (5.57)
n=1
L
We can interpret this by noting that p(x) was the initial “shape” of the wave,
i.e. y(x, 0) = p(x). Therefore p(x + ct) is just the original shape of the wave with
the origin of the x coordinate moved back from x = 0 to x = −ct. Therefore
the solution given by equation (5.56) represents the initial data splitting into two
travelling waves, one travelling to the left and one to the right, both being the
same “shape” as the initial data but with half the amplitude. Each wave moves
with speed c, giving the meaning of the wave speed.
74
It should be noted that this solution is not genuinely in the separation of vari-
ables form, however, this form, suitably re-interpreted, is d’Alembert’s solution
as seen in section 5.3.
Energy
We showed earlier in section 5.2.2 that the energy of the wave is
wL 1 ∂y 2
" 2 #
∂y
E= ρ + c2 dx. (5.58)
0
2 ∂t ∂x
If we substitute in our general solution given by equation (5.52) then we find,

on integration over the domain, that all the cross terms between sin and cos or
between sin and sin with different arguments, for example, will vanish. Hence we
find that
X
Total energy at time t = Energy in normal mode n at time t. (5.59)
normal
modes
It is a straightforward exercise to check, by using the equation for a single

normal mode, equation (5.50), in the equation for the energy, equation (5.58), that
1 nπc 2 2
Cn + Dn2 .

Energy in normal mode n = ρ (5.60)
4 L
This is independent of time. In particular, this shows that there is no exchange of
energy between different normal modes.
75
Chapter 6
Heat equation
6.1 Heat conduction in a rod

The model problem for all diffusive equations is that of heat transport in a solid
rod. The heat can only move along the axis, as we assume that it is perfectly
insulated in the other directions. The key assumption, which is borne out by ex-
periment, is that if we consider two parallel cross sections of area A separated by
a small distance d, an amount of heat (per unit time) will pass from the warmer
section to the cooler. In particular we will assume that the amount of heat trans-
ferred is proportional to the area A and the temperature difference |T2 − T1 | and
inversely proportional to the separation d, giving
A|T2 − T1 |
Amount of heat transferred per unit time = κ . (6.1)
d
The proportionality constant κ, called the thermal conductivity, is (to first approx-
imation) dependent only on the material of the rod. This is an empirical relation
which has its problems from the mathematical physics point of view, but leads to
an equation that produces accurate results in many cases.
We will consider the temperature as a function of x and t only, assuming that
it is constant in any cross section at constant x. We consider the section of the
bar such that x ∈ [x0 , x0 + δx], where x0 is arbitrary and δx small. From the
empirical law in equation (6.1) we have that the instantaneous rate of heat transfer
H(x0 , t) from left to right into the section under consideration through the surface
at x = x0 is given by
T x0 + d2 , t − T x0 − d2 , t

H(x0 , t) = − lim κA (6.2a)
d→0 d
∂T
= −κA (x0 , t). (6.2b)
∂x
76
HEAT EQUATION HEAT CONDUCTION
Figure 6.1: Schematic picture of heat flow in a section of a rod
This indicates that there will be positive heat flow into the section from the left if
the temperature is greater to the left of x = x0 , as we would expect. Repeating
the calculation for heat flow out of the section to the right through the surface at
x = x0 + δx gives
∂T
H(x0 + δx, t) = −κA (x0 + δx, t). (6.3)
∂x
Therefore the net rate at which heat flows into (or out of) the section of the bar
between x0 and x0 + δx is given by
Q = H(x0 , t) − H(x0 + δx, t) (6.4a)

∂T ∂T
= κA (x0 , t) − (x0 + δx, t) , (6.4b)
∂x ∂x
and hence the amount of heat entering the section in time δt is

∂T ∂T
Q δt = κA (x0 , t) − (x0 + δx, t) δt. (6.5)
∂x ∂x
This has to be balanced by the amount of heat actually absorbed by this section
of the bar. The average change in temperature δT in time δt is proportional to the
heat that flows into the section, Q δt and inversely proportional to the mass of the
section. As the mass is proportional to the volume of the section A δx, we write
this as
Q δt
δT = , (6.6)
sA δx
where s is related to (but is not exactly) the specific heat of the material of the bar.
Now, for a sufficiently small section of bar the average temperature change in the
section will be
∂T
δT ≈ δt (6.7a)
∂t
Q δt
= (6.7b)
sA δx
77
HEAT EQUATION SEPARATION OF VARIABLES
by using equation (6.6). Matching this against the expression for the heat flowing
in to the bar, equation (6.5), gives
κ ∂T (x0 , t) − ∂T

∂T ∂x ∂x
(x0 + δx, t) δt
δt ≈ . (6.8)
∂t s δx δt
Taking the limits δx → 0, δt → 0 gives
∂T ∂ 2T
= α2 2 . (6.9)
∂t ∂x
We note that the standard notation that we will use from now on is
∂y ∂ 2y
(x, t) = κ2 2 (x, t). (6.10)
∂t ∂x
Here κ is some constant related to the thermal conductivity of the material, or in
general the diffusivity of the problem that we are interested in.
6.2 Separation of Variables

We start with the simplest problem of the heat equation
∂y ∂ 2y
(x, t) = κ2 2 (x, t) (6.10)
∂t ∂x
with Dirichlet boundary conditions
y(0, t) = 0 = y(L, t), t > 0, (6.11)
and initial conditions

y(x, 0) = p(x) (6.12)
for given p.
The standard ansatz
y(x, t) = X(x)T (t) (6.13)
when substituted into the heat equation (6.10) gives
X Ṫ = κ2 X 00 T (6.14a)
or, on dividing by y = XT ,
1 Ṫ X 00
= . (6.14b)
κ2 T X
78
HEAT EQUATION SEPARATION OF VARIABLES
Again the two sides are separately constant, and by the same argument made
as for the wave equation in section 5.4.1 we can write the separated equations as
X 00 − λX = 0, (6.15a)
Ṫ − κ2 λT = 0. (6.15b)
Again, λ is some constant whose value is to be determined.
At this point we need boundary conditions in order to solve the ODEs give by
equation (6.15). Exactly as in the wave equation case the boundary conditions in
equation (6.11) imply Dirichlet boundary conditions for X at x = 0, L, giving the
fully-posed ODE problem
X 00 − λX = 0, X(0) = 0, X(L) = 0. (6.16)
This is identical to the problem in the wave equation case (equation (5.39) in
section 5.4.2). Hence we can write down the solution as
nπx nπ 2
Xn (x) = An sin , λn = − , n ∈ Z. (6.17)
L L
Here An is an undetermined constant.
We now need to solve for the time dependence, which is given by equa-
tion (6.15b). Once the eigenvalue λ is known, this is a simple first order ODE.
Given that the results for the spatial dependence, summarized in equation (6.17),
tell us that λ is negative, we can write down the solution to the time dependence
as
Tn (t) = Cn exp λn κ2 t

(6.18a)

nπ 2 2
= Cn exp − κt . (6.18b)
L
Again, Cn is an undetermined constant.
Combining our solutions and superposing gives the general solution
∞ nπx
X nπ 2 2
y(x, t) = Cn sin exp − κt . (6.19)
n=1
L L
Once again we have abused notation to absorb the free constants into the single
set of constants Cn .
Again we determine the free constants by noting that our initial conditions
given by equation (6.12) must match our general solution of equation (6.19) when
evaluated at t = 0, to give
X∞ nπx
p(x) = Cn sin . (6.20)
n=1
L
This is exactly the same Fourier Series problem as in the wave equation case.
79
HEAT EQUATION MORE COMPLEX MATERIAL PROPERTIES
6.3 More complex material properties

In the initial derivation of the heat equation given in section 6.1 we assumed that
all of the material properties of the rod, such as the conductivity κ, the cross-
sectional area A and the specific heat s were constant within the rod. It is a key
part of the derivation that they are constant within any cross-section. However, a
first generalization would allow them to depend on the x location within the rod.
Allowing the parameters to vary in this fashion obviously modifies the heat
flux used in equation (6.2) and later to give, for example,
∂T
H(x0 , t) = −κ(x0 )A(x0 ) (x0 , t). (6.21)
∂x
When these expressions are used to compute the net heat transfer Q and balanced
against the heat absorption we obtain a more general heat conduction equation

∂y ∂ ∂y
r(x) = p(x) . (6.22)
∂t ∂x ∂x
The spatial dependence of the various parameters has been collected into the func-
tions r and p which are positive on empirical grounds.
Additional complexity can be added by relaxing the assumption that the rod is
perfectly insulated and adding a source or sink term that adds or removes heat. To
make the problem tractable it is typically assumed that this source term G(x, t, T )
depends linearly on the temperature T and that the linear term is time independent,
i.e.
G(x, t, T ) = F (x, t) − q(x)T (x, t). (6.23)
When adding this source term to our increasingly complex heat conduction model
of equation (6.22) we obtain the generalized heat conduction equation

∂y ∂ ∂y
r(x) = p(x) − q(x)y + F (x, t). (6.24)
∂t ∂x ∂x
6.3.1 Separation of variables

We take our generalized heat conduction equation

∂y ∂ ∂y
r(x) = p(x) − q(x)y + F (x, t) (6.24)
∂t ∂x ∂x
and for the moment consider the homogeneous problem where F ≡ 0, giving

∂y ∂ ∂y
r(x) = p(x) − q(x)y. (6.25)
∂t ∂x ∂x
80
HEAT EQUATION MORE COMPLEX BOUNDARY CONDITIONS
We use our standard separation of variables ansatz, y = X(x)T (t) to obtain

0
r(x)X(x)Ṫ (t) = (p(x)X 0 (x)) T (t) − q(x)X(x)T (t). (6.26)
Dividing both sides by r(x)y(x, t) = r(x)X(x)T (t) gives
Ṫ (p(x)X 0 (x))0 q(x)

(t) = − . (6.27)
T r(x)X(x) r(x)
Once again we can note that both sides must be separately constant, leading to the
ODEs
Ṫ + λT = 0, (6.28a)
0 0
(p(x)X (x)) − q(x)X(x) + λr(x)X(x) = 0. (6.28b)
Once again we have an eigenvalue problem for the unknown separation con-
stant λ. However, the spatial dependence governed by equation (6.28b) is clearly
far more complex than before. However, we can note that equation (6.28b) is ex-
actly the form of a Sturm-Liouville problem as given in section 1.2.2; in particular,
equation (6.28b) should be compared to equation (1.60).
Particular results from Sturm-Liouville theory can immediately be used in the
solution of more complex boundary value problems such as this. For example,
the theory immediately tells us which boundary conditions guarantee a solution
(see section 1.2.2), that eigenvalues and eigenfunctions with qualitatively similar
properties to those seen earlier exist (see section 1.2.2), and that there are orthog-
onality conditions that can be used in place of the Fourier Series techniques to
determine the resulting free coefficients. Therefore once the solution to the spatial
dependence given by equation (6.28b) is found, all of the standard separation of
variables steps can be applied.
6.4 More complex boundary conditions

The results from Sturm-Liouville theory in section 1.2.2 suggest that we should
be able to solve the standard heat equation
∂y ∂ 2y
(x, t) = κ2 2 (x, t) (6.10)
∂t ∂x
with a broader range of boundary conditions than the simple Dirichlet conditions
y(0, t) = 0 = y(L, t), t > 0. (6.11)
81
HEAT EQUATION MORE COMPLEX BOUNDARY CONDITIONS
For the generalized (homogeneous) heat equation

∂y ∂ ∂y
r(x) = p(x) − q(x)y (6.25)
∂t ∂x ∂x
we can compare the results of the separation of variables ansatz in section 6.3.1
to those of Sturm-Liouville theory in section 1.2.2 to see that any appropriate
combination of
1. y = 0: Dirichlet condition;
2. y 0 = 0: Neumann condition;
3. y + ky 0 = 0: Radiation condition;
4. p = 0: x = x0 , x1 are singular points of the ODE;
5. y periodic;
fits the conditions necessary for separation of variables to be a practical technique.

For the standard heat equation (6.10) we have p = 1 so either the solution is
periodic or some appropriate combination of the Dirichlet, Neumann and Radia-
tion conditions must hold. We consider applying the general combined boundary
condition
α1 y(0) + β1 y 0 (0) = 0, α2 y(L) + β2 y 0 (L) = 0. (6.29)
Applied to the standard heat equation and using separation of variables as normal
gives the generalization of the ODE in equation (6.16)
X 00 − λX = 0, α1 X(0) + β1 X 0 (0) = 0, α2 X(L) + β2 X 0 (L) = 0. (6.30)
We look at the case λ = µ2 > 0 where the auxiliary equation has distinct
real roots ±µ. Rather than writing the general solution to the ODE in terms of
exponentials we instead write it as
X(x) = A sinh(µx) + B cosh(µx). (6.31)
It is straightforward to check that the boundary conditions imply that
0 = X(0) = α̃1 sinh(µ0) + β̃1 cosh(µ0), (6.32a)

0 = X(L) = α̃2 sinh(µL) + β̃2 cosh(µL). (6.32b)
It is straightforward to check that the only solution compatible with these bound-
ary conditions is the trivial solution. The same is true if we attempt to impose that
82
HEAT EQUATION INHOMOGENEOUS PROBLEMS
the solution y is periodic; this implies that X is periodic, which is not possible for
the hyperbolic trigonometric functions that make up this solution.
In contrast, the case λ = 0 may have a solution. The general solution of
equation (6.30) when λ = 0 is obviously
X(x) = Ax + B. (6.33)
In most cases we can still check that A = 0 (but this need not be true, for exam-
ple, for the radiation condition!), but the solution where B 6= 0 is compatible with
periodic boundary conditions, and with a range of Dirichlet and radiation condi-
tions. This leads to additional solutions for the time dependence that are similarly
of the form
T (t) = Ct + D. (6.34)
Typically the solution being bounded either in the past or the future implies that
C = 0. However, the constant solution is again frequently valid. It should be
noted that in more complex domains these additional solutions often lead to con-
siderably more complex behaviour that must be considered.
6.5 Inhomogeneous problems

6.5.1 Inhomogeneous equations
In this section we consider a small extension to the standard heat equation
∂y ∂ 2y
(x, t) = κ2 2 (x, t) (6.10)
∂t ∂x
suggested by the generalized heat equation (6.24); we add the source term F to
get
∂y ∂ 2y
(x, t) = κ2 2 (x, t) + F (x, t). (6.35)
∂t ∂x
For simplicity and for later use we will assume trivial Dirichlet boundary condi-
tions
y(0, t) = 0 = y(1, t) (6.36)
on the simple domain x ∈ [0, 1].
Clearly for general source term F we cannot use separation of variables, as our
ansatz y = XT does not lead to an equation containing only separably constant
terms. Instead we make the assumption that the spatial dependence of the solution
can still be described by a Fourier Series, and because of the simple boundary
83
conditions of equation (6.36), this implies a Fourier Sine Series. This gives us the
more general ansatz
X∞
y(x, t) = yn (t) sin(nπx). (6.37)
n=1
We still cannot start solving our problem without making additional assump-
tions. However, we can substitute our ansatz into the inhomogeneous heat equa-
tion (6.35) and see what results. We find
∞
X
ẏn (t) − (nπ)2 yn (t) sin(nπx) = F (x, t).

(6.38)
n=1
At any instant in time the left hand side is a Fourier Sine Series. Therefore in order
to be consistent the source term F must also be a Fourier Sine Series. Making that
assumption,
X∞
F (x, t) = Fn (t) sin(nπx), (6.39)
n=1
we can substitute this in to the inhomogeneous heat equation as well to find

∞
X
ẏn (t) − (nπ)2 yn (t) − Fn (t) sin(nπx) = 0.

(6.40)
n=1
By using the standard orthogonality relations we immediately see that
ẏn (t) − (nπ)2 yn (t) − Fn (t) = 0. (6.41)
This is a set of ODEs to solve for the unknown yn (t) in terms of the Fn , which
can be computed from the known source term F . In particular, this linear ODE
has the obvious solution
wt
yn (t) = Cn exp −(nπ) t + Fn (s) exp −(nπ)2 s ds.
2

(6.42)
0
Once Fn has been computed we can compute yn and our original ansatz in equa-
tion (6.37) gives the solution.
6.5.2 Inhomogeneous boundary conditions

In this section we consider the standard heat equation
∂y ∂ 2y
(x, t) = κ2 2 (x, t) (6.10)
∂t ∂x
84
with the inhomogeneous boundary conditions
y(0, t) = f (t), (6.43a)

y(1, t) = g(t). (6.43b)
We have assumed that the domain has length 1 for simplicity and that f, g are
given functions.
It should immediately be obvious that separation of variables fails. This is
because we can no longer solve the separated problem for the spatial dependence,
as the boundary conditions are functions of time. However, we can use a simple
trick inspired by techniques for solving inhomogeneous ODE problems to reduce
this inhomogeneous PDE to a form that can be solved by separation of variables.
First we note that the heat equation (6.10) is linear so we can superpose solu-
tions, even when we have complex boundary conditions such as equation (6.43).
Therefore if we had a particular solution yP (x, t) that satisfies the heat equation
and the boundary conditions (6.43), we could superpose it with another solution
yG (x, t) of the heat equation, and the combined solution
y(x, t) = yP (x, t) + yG (x, t) (6.44)
will satisfy the heat equation. More importantly, we note that if the solution
yG (x, t) vanishes at the boundaries x = 0, 1 then the combined solution y(x, t)
also satisfies the boundary conditions (6.43).
Unfortunately it is just as difficult to find a single solution to the full inho-
mogeneous problem as it is to find the general solution. Instead we suppose we
have a particular solution yP that satisfies the boundary conditions but not the
heat equation. When we form the combined solution using equation (6.44) we
find that the combined solution satisfies the boundary conditions but not the heat
equation (6.10). It will however satisfy the inhomogeneous heat equation (6.35),
where the source term F (x, t) can be computed directly from yP .
Thus we have reduced the problem of finding the general solution of the heat
equation (6.10) with complex boundary conditions (6.43) to the problem of find-
ing a simple function yP that satisfies the complex boundary conditions (6.43),
plus the problem of solving the inhomogeneous heat equation (6.35) with trivial
Dirichlet boundaries as studied in section 6.5.1.
We still have the problem of finding a sufficiently simple particular solution
yP . However, this can itself be simplified by assuming that it separates, yP (x, t) =
XP (x)TP (t). In particular, if we choose XP to be a linear polynomial, then it
trivially vanishes under the second derivative in the heat equation. We then choose
T to satisfy the boundary conditions. In particular the choice
yP (x, t) = xg(t) + (1 − x)f (t) (6.45)
85
has the correct behaviour at the boundaries and vanishes under the second spatial
derivative. Therefore we can see that this solution satisfies the equation
∂yP ∂ 2 yP
(x, t) = (x, t) + xġ(t) + (1 − x)f˙(t). (6.46)
∂t ∂x2
This gives the simple source term
F (x, t) = xġ(t) + (1 − x)f˙(t), (6.47)
meaning that the general solution of the heat equation with inhomogeneous bound-
ary conditions (6.43) reduces to
y(x, t) = yP (x, t) + yG (x, t) (6.48)
where yG is the general solution of the inhomogeneous heat equation (6.35) with
source term (6.47) and trivial Dirichlet boundary conditions.
86
Chapter 7
Laplace’s equation
7.1 Introduction
Laplace’s equation is simply
∇2 φ = 0. (7.1)
In R1 and Cartesian coordinates Laplace’s equation reduces to
d2 φ
= 0, (7.2)
dx2
where φ(x) represents, e.g., the static shape of a stretched string, or the steady
temperature in a conducting bar, or something similar.
2
∂2

∂
+ φ = 0, (7.3)
∂x2 ∂y 2
where φ(x, y) represents, e.g., the static shape of a stretched membrane such as a
soap film, or the steady temperature in a conducing sheet, or something similar.
2
∂2 ∂2

∂
+ + φ = 0, (7.4)
∂x2 ∂y 2 ∂z 2
where φ(x, y, z) represents, e.g., the steady temperature in a conducing solid, or

the velocity potential in an inviscid fluid, or the electrostatic potential in a charge-
free region, or the gravitational potential in a mass-free region, or something sim-
ilar.
In general we usually call φ the potential, and F = (−)∇φ is the correspond-
ing field.
87
LAPLACE’S EQUATION SEPARABLE SOLUTIONS
7.1.1 Boundary conditions

In order to pick out a unique solution it is necessary to specify boundary conditions
at the edges of the domain. The natural boundary conditions are to specify either
1. φ: Dirichlet boundary conditions, or
2. n̂ · ∇φ: Neumann boundary conditions,
at each point of the bounding surface. We note that this is two points in R1 , a
bounding curve in R2 , and a bounding surface in R3 . The (unit) vector n̂ is the
normal to the bounding surface.
7.2 Separable solutions of Laplace’s equation

We consider the two dimensional Laplace equation
2
∂2

∂
+ φ=0 (7.3)
∂x2 ∂y 2
on the unit square with trivial Dirichlet boundary conditions,

2
∂2

∂
+ φ = 0, 0 ≤ x, y ≤ 1,
∂x2 ∂y 2 (7.5)
φ(0, y) = φ(1, y) = φ(x, 0) = φ(x, 1) = 0.
We make the standard separation of variables assumption, φ(x, y) = X(x)Y (y).

As usual this gives two ODEs
X 00 − λX = 0, (7.6a)
Y 00 + λY = 0. (7.6b)
Both are BVPs with trivial boundary conditions. This is different to the situation
for the heat and wave equations where only one of the ODEs was a BVP.
We solve the eigenvalue problem for the x dependence as normal, finding that
the only non-trivial solution is
λn = −(nπ)2 , Xn (x) = sin (nπx) . (7.7)
We are left with the solution for the y dependence, and using this value of the
separation constant λn in the ODE of equation (7.6b) gives
Yn = An sinh (nπy) + Bn cosh (nπy) . (7.8)
88
The trivial Dirichlet boundary conditions immediately imply that the only solution
compatible with the boundary conditions is the trivial solution Yn ≡ 0 and hence
φ(x, y) ≡ 0.
This illustrates two points. One is the maximum principle: solutions of Laplace’s
equation always take their extreme values (both maxima and minima) at the bound-
aries of the domain. This holds in all dimensions and (reasonable) domains. Al-
though not proved or used in this course, it is easily motivated by considering the
stationary limit of the heat equation.
The second point is that for the vast majority of interesting problems using
Laplace’s equation we must expect them to be inhomogeneous, either through
their boundary conditions or a source term (Poisson’s equation). There is a straight-
forward method of dealing with these problems, exploiting the linearity of the
equation.
7.2.1 Cartesians
For simplicity we explain the technique using the two-dimensional Laplace equa-
tion 2
∂2

∂
+ φ=0 (7.3)
∂x2 ∂y 2
on the unit square (hence using Cartesian coordinates) with some set of boundary
conditions
φ(x, 1) = Φ1 (x), (7.9a)

φ(1, y) = Φ2 (y), (7.9b)
φ(x, 0) = Φ3 (x), (7.9c)
φ(0, y) = Φ4 (y), (7.9d)
which implicitly number the boundaries of the domain in a clockwise sense.

The technique has two steps.
1. Change variable from φ(x, y) to a function u(x, y) = φ(x, y) − v(x, y)
which satisfies the Laplace equation in the interior and also satisfies bound-
ary conditions Ui on boundary i that vanish in the corners of the domain.
2. Solve for u(x, y) by writing it as
u(x, y) = u1 (x, y) + u2 (x, y) + u3 (x, y) + u4 (x, y), (7.10)
where each ui satisfies Laplace’s equation in the interior, satisfies the bound-
ary condition Ui on boundary i and vanishes (i.e., satisfies a trivial Dirichlet
boundary condition) on all other boundaries.
89
Once we have found each of the ui , we can reconstruct the full solution to the
original problem using
φ(x, y) = v(x, y) + u1 (x, y) + u2 (x, y) + u3 (x, y) + u4 (x, y). (7.11)
Finding the difference function

The first required step is to find the difference function v(x, y). This is constructed
so that u has boundary conditions taking trivial values in the corners of the do-
main. In general it should be clear that as φ satisfies Laplace’s equation, and we
want u also to do so, v must satisfy Laplace’s equation. It should also be clear
that v must take the same value as the boundary conditions Φi in the corners of
the domain.
This is a considerably easier problem, as we can take any solution of Laplace’s
equation that is constrained at only four points. In Cartesians the simplest solu-
tions of Laplace’s equation are the low order polynomials 1, x, y, xy. Therefore
we can define our difference function to be
v(x, y) = α1 + α2 x + α3 y + α4 xy. (7.12)
The values of the α coefficients is set by the restrictions given at the corners of the
domain, so that
v(0, 0) = α1 = Φ3 (0) = Φ4 (0) (7.13a)

v(1, 0) = α1 + α2 = Φ2 (0) = Φ3 (1) (7.13b)
v(0, 1) = α1 + α3 = Φ1 (0) = Φ4 (1) (7.13c)
v(1, 1) = α1 + α2 + α3 + α4 = Φ1 (1) = Φ2 (1). (7.13d)
This simple linear system can be solved by forward substitution, provided the
boundary conditions are consistent in the corners.
Solving the sub-problems

We now have a Laplace equation problem for u subject to boundary conditions
that vanish in the corners. We have to solve the four sub-problems for ui ; for
example, the problem for u1 (x, y) is
2
∂2

∂
+ u1 = 0,
∂x2 ∂y 2 (7.14)
u1 (x, 1) = U1 (x), u1 (1, y) = u1 (x, 0) = u1 (0, y) = 0.
90
Using standard separation of variables techniques as in section 7.2 we can write

u1 = X(x)Y (y), and again compute the separation constant and x dependent
behaviour to find
λn = −(nπ)2 , Xn (x) = sin (nπx) (7.15)
as in equation (7.7). We then solve for the y dependence as in equation (7.8), and
the trivial boundary condition at y = 0 restricts the solution to
Yn = An sinh (nπy) . (7.16)
We can once again superpose solutions, finding the result
∞
X
u1 (x, y) = An sin (nπx) sinh (nπy) , (7.17)
n=1
where the value of the coefficients An is set by the boundary condition at y = 1.

This gives
∞
X
U1 (x) = Cn sin (nπx) (7.18)
n=1
where the new constant Cn is
Cn = An sinh (nπ) . (7.19)
Using standard Fourier Series techniques we can find Cn , from which we can find
An , and hence find u1 . Repeating this for all boundaries gives the solution to
each sub-problem ui , hence the solution for u, from which we can reconstruct the
solution to the original problem.
7.2.2 Other domains

The essential features of the problem remain as in Cartesian coordinates even
when the domain is not a simple Cartesian rectangle. However, it may be simpler
to solve. For example, if the problem is Laplace’s equation in two dimensions on
a circular disk, or on an annulus, then it is simpler to express Laplace’s equation
in polar coordinates as
1 ∂ 2φ

1 ∂ ∂φ
r + 2 2 = 0, (7.20)
r ∂r ∂r r ∂θ
where φ ≡ φ(r, θ). The domains for the disk D and annulus A are then simply
expressed as
D = {(r, θ) | 0 ≤ r < R, 0 ≤ θ < 2π} , (7.21a)
A = {(r, θ) | R1 ≤ r < R2 , 0 ≤ θ < 2π} . (7.21b)
91
The constant R defines the radius of the disk and the constants R1,2 the inner and
outer radii of the annulus.
There are two key points to note. The first is that the boundaries are smooth,
so there are no corners to worry about as in the Cartesian case. This means that the
step of finding the difference function as in section 7.2.1 is considerably simpler
– only one boundary needs considering in the case of the annulus, and the step is
completely unnecessary for the disk.
The second point is that there is no explicit boundary condition in the angular
direction, as there is no explicit boundary. However, in order to have a single-
valued solution we need φ to be 2π-periodic. In addition, in the case of the disk
we require the solution to be finite at r = 0. These combine with the standard
separation of variables techniques to give solutions of the form
∞
1 X
φ(r, θ) = a0 + rn (an cos(nθ) + bn sin(nθ)) (7.22)
2 n=1
for the disk problem, and imposing boundary conditions at r = R and using
orthogonality will give the free constants an , bn .
92
Chapter 8
Vector Calculus
8.1 Vector and scalar valued functions and their deriva-

tives
8.1.1 Curves and vector valued functions of 1-variable
A vector valued function is a vector r which depends on a real variable t say.
We write such a vector valued function as r(t). Sometimes t will take on all real
values and sometimes t will be in some interval, I. If we take a < t < b then I is
an open interval and we write I = (a, b) however if we include the end points so
that a ≤ t ≤ b then I is a closed interval and we write I = [a, b].
More formally we think of a vector valued function as a map
I → R3
t 7→ r(t)
For most of theses notes you may assume, unless specified otherwise, that the
map is smooth. This means that one can differentiate the map as often as one
likes. However in some of the theorems we need to be careful about the degree
of differentiability of functions. A continuous function is said to be C 0 , a differ-
entiable function with continuous first derivative is said to be C 1 , whilst a twice
differentiable function with continuous second derivative is sad to be C 2 , and so
on.
As t varies between a and b the vector r(t) traces out a curve.
We call this set of points in R3 a curve. However we refer to the map
I → R3
t 7→ r(t)
93
VECTOR CALCULUS VECTORS AND SCALARS
Figure 8.1: As t varies between a and b the vector r(t) traces out
a curve.
as a parameterised curve. This is because it is possible to obtain the same curve

in R3 with a different parameterisation. Let f : R → R be a monotonic increasing
(smooth) function. Let f (c) = a and f (d) = b and define I˜ = [c, d]. Then we
may define a new vector valued function r̃(s) = r(f (s)) for c ≤ s ≤ d. So that
I˜ → R3
s 7→ r̃(s)
is also a parameterised curve. However as s varies between c and d, r̃ traces out

the same set of points as r. Thus r and r̃ are simply different parameterisations
of the same curve.
The tangent to a parametrised curve r(t) is given by

r(t + h) − r(t)
v(t) = ṙ(t) = lim
h→0 h
If we give the components of r(t) as
r(t) = x(t)i + y(t)j + z(t)k
The ṙ(t) has components
ṙ(t) = ẋ(t)i + ẏ(y)j + ż(t)k
94
If we think of the parameter t representing time, then ṙ(t) is the velocity and
v(t) = |ṙ(t)| is the speed. This gives us the following formula for the length of
the curve between the point A with position vector r(t0 ) and the point B with
position vector r(t1 ):
Z t1 Z t1 Z t1
1/2
ẋ(t)2 + ẏ(t)2 + ż(t)2

L= v(t)dt = |ṙ(t)|dt dt
t0 t0 t0
What happens if we use some other parameterisation of the same curve? Let
s = f (t) and r̃(s) = r(f (s)) as above. Then by the chain rule
dr̃(s) dr(t) dt
=
ds dt ds
Let f (t0 ) = s0 and f (t1 ) = s1 . Then the length of the curve from A to B using
the s parameterisation is given by
Z s1
dr̃(s)
L̃ = ds ds

s
Z 0s1
dr(t) dt
= dt ds ds

s
Z 0s1
dr(t) dt
= dt ds ds

s0
Z t1
dr(t)
= dt dt

t0
=L
Hence the length of the curve is independent of the parameterisation. This feature
of the integral is called re-parameterisation invariance.
8.1.2 Scalar Fields

A scalar field is simply a scalar valued function that depend upon position. For
example the temperature in a room. In 2-dimensional space a scalar field is a map
φ : R2 → R
(x, y) 7→ φ(x, y)
An example of a scalar field in 2-dimensions is therefore given by

1 3 1 7
φ(x, y) = y − y − x2 +
12 4 2
95
Figure 8.2: Three dimensional plot of the function φ(x, y) =

1 3
12
y − y − 41 x2 + 72 .
To visualise this function we can plot the graph of the surface given by
z = φ(x, y)
using Maple, for example, see Fig. 8.2.

An alternative is to plot the contour lines or level surfaces given by z = c =constant.
Thus we plot the curves
1 3 1 7
y − y − x2 + = c
12 4 2
for a number of equally spaced contours. For the above function this gives the
picture in Fig. 8.3.
Example 2: For the scalar field given by f (x, y) = 4 − x2 − y 2 plot the graph
z = f (x, y) and the contour lines f (x, y) = c.
Often it is easier to draw the contour lines so we start with these. For this example
2 2
√given by x + y = 4 − c. These are a series of circles centre the origin
they are
radius 4 − c. Note that for real solutions we must have c ≤ 4, we also note that
the curves get closer together as the radius of the circles increases. This indicates
that the surface gets steeper as we move away from the origin.
96
y
4
0 x
-2
-4
-4 -2 0 2 4
1 3
Figure 8.3: Contour plot of the function φ(x, y) = 12
y −y−
1 2
4
x + 72 .
We now consider scalar functions in 3-dimensions. These are maps

φ : R3 → R
(x, y, z) 7→ φ(x, y, z)
The graphs of this function w = φ(x, y, z) are 3-dimensional hypersurfaces in

4-dimensional space - which are not so easy to draw! Instead we draw the level
surfaces given by
φ(x, y, z) = c
As with contour lines we can visualise the function by drawing a series of equally
spaced level surfaces. These are 2-dimensional surfaces in R3 so drawing them is
feasible.
Example: Draw the level surfaces of the function φ(x, y, z) = x + y + z. These
are simply given by x + y + z = c, which we know from chapter 1 consist of
planes with normal n = i + j + k.
97
8.1.3 Vector fields

A vector field is a vector valued function that depends upon position. In 2-
dimensions it is a map
F : R2 → R2
(x, y) 7→ F(x, y) = F1 (x, y)i + F2 (x, y)j
For example F(x, y) = yi − xj. How do we visualise a vector field?

At the point (x, y) we draw the vector F(x, y). For the above example this gives
the picture:
Example 2: Let r = xi + yj denote a typical point in R2 . Let the vector field

F(r) be defined by
r k
F(r) = −k 3 = − 2 r̂
|r| r
which represents an inverse square force law (such as gravitation or electromag-
netism). The picture we get of this vector field is as follows:
98
Note: The vector field is not defined at the origin.
8.1.4 Directional derivatives

Let φ(r) be a scalar field. What is the derivative of φ? It depends on what direction
you move in!
Let v be some fixed vector, then the directional derivative of φ in the v direction,
at the point r0 is given by

φ(r0 + tv) − φ(r0 )
∇v φ(r0 ) = lim
t→0 t
This measures the change in φ in the v direction.

This formal definition is not very useful when it comes to doing calculations so
we will convert it into something that we can use.
To do this we think of r0 and v as fixed and define the real valued function f (t) =
φ(r0 + tv). Then

φ(r0 + tv) − φ(r0 )
∇v φ(r0 ) = lim
t→0 t

f (t) − f (0)
= lim
t→0 t
0
= f (0)
99
We can write f (t) as:

f (t) = φ(x, y, z) (8.1)
where
x = r1 + tv1
y = r2 + tv2
x = r3 + tv3
and r0 = r1 i + r2 j + r3 k and v = v1 i + v2 j + v3 k.
Differentiating (8.1) with respect to t using the chain rule gives
df ∂φ dx ∂φ dy ∂φ dz
= + +
dt ∂x dt ∂y dt ∂z dt
∂φ ∂φ ∂φ
= v1 + v2 + v3
∂x ∂x ∂z
If we introduce the vector field gradφ (denoted ∇φ) defined by

∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
Then we see that
∇v φ = v·∇φ
Example: Let φ(x, y, z) = xy +z 2 . What is the rate of change of φ in the direction
of v̂, where v = 3i + 4k, at the point x = 1, y = 2, z = 1?
We first calculate ∇φ;
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
= yi + xj + 2zk
= 2i + j + 2k at the point (1, 2, 1)
Now |v| = 5 so that

3 4
v̂ = i + k
5 5
and hence
3 4 14
∇v φ(1, 2, 1) = i + k ·(2i + j + 2k) =
5 5 5
100
8.1.5 Gradient and steepest ascent

When calculating ∇φ using the formula above we need to know φ as a function of
the Cartesian coordinates x, y and z. However we sometimes use other coordinate
systems. In this section we give a geometrical definition of gradφ which enables
us to calculate it in any coordinate system.
Given a scalar field φ(r) in what direction is it increasing most? Let v̂ be a unit
vector then
∇v φ = v̂·∇φ
= |v̂||∇φ| cos θ
= |∇φ| cos θ
This is a maximum when θ = 0 and a minimum when θ = π, so that
−|∇φ| ≤ ∇v φ ≤ |∇φ|
We have therefore established the following theorem:

Theorem 1: The directional derivative ∇v φ is maximised when v̂ points in the
same direction ∇φ and minimised when it points in the opposite direction. Fur-
thermore the max and min values are given by ±|∇φ|.
8.1.6 Gradients and level surfaces

We want to look at ∇φ at the point P with coordinates (x0 , y0 , z0 ). Let φ(x0 , y0 , z0 ) =
c0 . Now consider the level surface through the point P which is given by the set
of points (x, y, z) such that φ(x, y, z) = c0 .
Now let v be a tangent vector at P (ie v is an element of the tangent plane at P ).

Then since v is tangent to the level surface φ(x, y, z) = c0 , the derivative of φ in
101
the v direction is zero. Thus
∇v φ = 0
⇒v·∇φ = 0
⇒the tangent vector v is perp. to ∇φ
⇒∇φis perpendicular to the tangent space at P
⇒∇φis normal to the level surface φ = c0
We have therefore established the following:

Theorem 2: The gradient of φ, ∇φ, is normal to the surfaces φ =const.
We now use the above two theorems to obtain a new formula for the gradient.
By Theorem 2, since ∇φ is normal to the surface φ=const. we know that
∇φ = k n̂ (8.2)
where k = |∇φ|.
On the other hand we know from Theorem 1 that the maximum value of ∇v φ is
|∇φ| and that it obtains this value when v points in the same direction as ∇φ, i.e.
when v̂ = n̂. Thus
|∇φ| = ∇n̂ φ (8.3)
Where we have taken n̂ to be the choice of normal that points in the direction of
increasing φ.
Combining formula (8.2) with formula (8.3) we obtain
∇φ = (∇n̂ φ) n̂
∂φ
To simplify notation we write to denote the directional derivative in the normal
∂n
direction ∇n̂ φ. Hence we can write
∂φ
∇φ = n̂
∂n
To summarise:
Gradient (geometric definition)
The gradient of the scalar field φ is the vector field ∇φ given by
∂φ
∇φ = n̂
∂n
102
∂φ
where n̂ is the unit normal to the surfaces φ = const and is the directional
∂n
derivative of φ in the n̂ direction.
Gradient (Cartesian definition)
If the scalar field φ is given in Cartesian coordinates by φ(x, y, z) then
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
Directional Derivative
The directional derivative

φ(r0 + tv) − φ(r0 )
∇v φ(r0 ) = lim
t→0 t
= v·∇φ
∂φ ∂φ ∂φ
= v1 + v2 + v3 (in Cartesian coordinates)
∂x ∂y ∂z
Example 1: Let a general point be denoted r = xi + yj + zk, and let r = |r|.

Consider the scalar field φ(r) = |r| = r. What is ∇φ?
Cartesian method:
φ(x, y, z) = (x2 + y 2 + z 2 )1/2 . So that
∂φ 1 x x
= 2x (x2 + y 2 + z 2 )−1/2 = 2 = .
∂x 2 (x + y 2 + z 2 )1/2 r
∂φ y ∂φ z
Similarly = , and = . Hence
∂y r ∂z r
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
x y z
= i+ j+ k
r r r
1 r
= (xi + yj + zk) = = r̂
r r
Geometric method:
The level surfaces of φ =const. are the surfaces r =const. which are simply
spheres centre the origin. Hence the normal to the surface points in the same
direction as the position vector. Thus n̂ = r̂. Hence
∂φ ∂φ
∇φ = n̂ = r̂ = 1.r̂
∂n ∂r
103
More generally if the scalar field φ only depends upon the distance r from the
origin and not the direction then we can write φ(r) = f (r). The level surfaces are
still spheres centre the origin and hence we still have n̂ = r̂. Thus
∂φ ∂φ df
∇φ = n̂ = r̂ = r̂ = f 0 (r)r̂
∂n ∂r dr
Example 2: Calculate the equation of the tangent plane to the surface S given by
x3 y − yz 2 + z 5 = 9 at the point P given by (3, −1, 2).
Let φ = x3 y − yz 2 + z 5 then at P , the value of φ = 9 so that S is the level surface
given by φ(x, y, z) = 9. Then
∇φ = 3x2 yi + (x3 − z 2 )j + (5z 4 − 2yz)k
So that at P , a normal to the surface is given by n = ∇φ = −27i + 23j + 84k.

The equation of the tangent plane at r0 is n·(r − r0 ) = 0 which is
(−27i + 23j + 84k)·((x − 3)i + (y + 1)j + (z − 2)k) = 0
or −27x + 23y + 84z = 64.
8.1.7 Properties of gradφ

(1.) Linearity: ∇(φ + ψ) = ∇φ + ∇ψ
∂ ∂φ ∂ψ
Proof: Use Cartesian coordinates and the fact that (φ + ψ) = +
∂x ∂x ∂x
(2.) Leibnitz rule: ∇(φψ) = ψ∇φ + φ∇ψ
∂ ∂φ ∂ψ
Proof: Use Cartesian coordinates and the fact that (φψ) = ψ +φ
∂x ∂x ∂x
8.1.8 The divergence of a vector field

In the previous section we saw how to calculate a vector field gradφ by taking the
derivative of a scalar field. In this section we show how to calculate a scalar field
by taking the derivative of a vector field F. We call this derived scalar field the
divergence of F and denote it divF or ∇·F.
For the special case where F is given by
F = φc
104
where φ is a scalar field and c is a constant vector, we define divF by the formula
divF = ∇φ·c
We also require that the divergence is a linear operator, so that if F and F̃ are two
vector fields then
div(F + F̃) = divF + divF̃
Using the fact that we may write a general vector field in a Cartesian basis as
F = F1 i + F2 j + F3 k
Then
divF = div(F1 i) + div(F2 j) + div(F3 k)
= ∇(F1 )·i + ∇(F2 )·j + ∇(F3 )·k
∂F1 ∂F2 ∂F3
= + +
∂x ∂y ∂z
Hence we have obtained the following result

Divergence in Cartesian coordinates Let F be a vector field with components
F = F1 i + F2 j + F3 k
Then the divergence of F is given by
∂F1 ∂F2 ∂F3
divF = + +
∂x ∂y ∂z
Note if we introduce the differential operator
∂ ∂ ∂
∇= + +
∂x ∂y ∂z
Then we may write the above as
∂F1 ∂F2 ∂F3
divF = ∇·F = + +
∂x ∂y ∂z
Example 1: Let F = y 2 zi + xzj − y 2 k Then
∂F1 ∂F2 ∂F3
∇·F = + +
∂x ∂y ∂z
∂y z ∂xz ∂ − y 2
2
= + +
∂x ∂y ∂z
=0+0+0=0
105
So that for F = y 2 zi + xzj − y 2 k then ∇·F = 0. This example shows that

divF = 0 does not imply that F is constant.
In fact the divergence measure the effect of the vector field on the volume of a
small sphere. If the flow lines are moving apart they diverge and the volume of
a small sphere moving with the flow increases, and the divergence is positive.
However if the flow lines are moving together they converge and the volume of
a small sphere moving with the flow decreases and the divergence is negative.
However it is possible for the flow lines to converge in some directions and diverge
in others so that there is no overall change in volume of a small sphere. In this
case the divergence vanishes (even if the sphere is distorted by the flow).
Proposition:
Let φ be a scalar field and F a vector field then
∇·(φF) = (∇φ)·F + φ∇·F
Proof: Use cartesian coordinates to calculate both sides.
8.1.9 The Laplacian

If we start with a scalar field φ we may first take the gradient to obtain ∇φ. This
is a vector field so we may take its divergence to obtain a scalar field, ∇·(∇φ).
We call this scalar field the “Laplacian” of φ.
Definition:
The Laplacian of a scalar field φ is denoted ∇2 φ and defined by
∇2 φ = div(gradφ) = ∇·(∇φ)
If we do the above calculation in Cartesian coordinates we have

∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
Taking the divergence of this we get

∂φ ∂φ ∂φ ∂φ ∂φ ∂φ
∇·(∇φ) = + +
∂x ∂x ∂y ∂y ∂z ∂z
∂ 2φ ∂ 2φ ∂ 2φ
= + + 2
∂x2 ∂y 2 ∂z
106
Hence in Cartesian coordinates

∂ 2φ ∂ 2φ ∂ 2φ
∇2 φ = + + 2
∂x2 ∂y 2 ∂z
Example:
Let φ(x, y, z) = ex y 3 sin z, calculate ∇2 φ.
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
= ex y 3 sin zi + 3ex y 2 sin zj + ex y 3 cos zk
∂F1 ∂F2 ∂F3

∇·(∇φ) = + +
∂x ∂y ∂z
∂(e y sin z) ∂(3ex y 2 sin z) ∂(ex y 3 cos z)
x 3
= + +
∂x ∂y ∂z
x 3 x x 3
= e y sin z + 6e y sin z − e y sin z
Alternatively
∂ 2φ ∂ 2φ ∂ 2φ
∇2 φ = + + 2
∂x2 ∂y 2 ∂z
∂ (e y sin z) ∂ 2 (ex y 3 sin z) ∂ 2 (ex y 3 sin z)
2 x 3
= + +
∂x2 ∂y 2 ∂z 2
x 3 x x 3
= e y sin z + 6e y sin z − e y sin z
8.1.10 The curl of a vector field

We have already considered a differential operator called the divergence that may
be symbolically written as ∇· in this section we consider a differential operator
called the “curl” or “rotation” which may be symbolically written as ∇×. As its
name suggests the “curl” measures the extent to which a small sphere rotates as it
moves along the flow lines of the vector field F.
For the special case where F is given by
F = φc
107
where φ is a scalar field and c is a constant vector, we define curlF by the formula
curlF = ∇φ×c
We also require that the curl is a linear operator, so that if F and F̃ are two vector
fields then
curl(F + F̃) = curlF + curlF̃
Using the fact that we may write a general vector field in a Cartesian basis as
F = F1 i + F2 j + F3 k
Then
curlF = curl(F1 i) + curl(F2 j) + curl(F3 k)
= ∇(F1 )×i + ∇(F2 )×j + ∇(F3 )×k

∂F1 ∂F1 ∂F2 ∂F2 ∂F3 ∂F3
= − k+ j + k− i + j− i
∂y ∂z ∂x ∂z ∂x ∂y

∂F3 ∂F2 ∂F3 ∂F1 ∂F2 ∂F1
= − i+ − j+ − k
∂y ∂z ∂x ∂z ∂x ∂y
Hence we have obtained the following result

Curl in Cartesian coordinates Let F be a vector field with components
F = F1 i + F2 j + F3 k
Then

∂F3 ∂F2 ∂F3 ∂F1 ∂F2 ∂F1
curlF = ∇×F = − i+ − j+ − k
∂y ∂z ∂x ∂z ∂x ∂y
This may be formally written using a determinant as

i j k
∂ ∂ ∂
∇×F = ∂x ∂y ∂z

F1 F2 F3
Example 1: Let F = xy 2 i + ez j + y sin xez k. Calculate curlF. Using the above

formula

i j k
∂ ∂ ∂

∇×F = ∂x ∂y = (sin xez − ez )i + (y cos xez )j + 2xyk
∂z
xy 2 ez y sin xez
108
Example 2: Let φ = x2 yz 3 and let F = ∇φ. Calculate ∇×F.
F = ∇φ = 2xyz 3 i + x2 z 3 j + 3x2 yz 2 k
Then using the above formula

i j k
∂ ∂ ∂

∇×F = ∂x

∂y ∂z

2xyz 3 x2 z 3 3x2 yz 2
= (3x2 z 2 − 3x2 z 2 )i + (−6xyz 2 + 6xyz 2 )j + (2xz 3 − 2xz 3 )k = 0
So in this example we find that ∇×∇φ = 0. In fact this is not just chance but is
a general result.
Theorem: Let F be a C 2 (ie twice differentiable with continuous second deriva-
tive) vector field then:
curl(gradφ) = ∇×∇φ = 0
Proof:
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
So that

i j k
∂ ∂ ∂
∇×∇φ = ∂x ∂y ∂z
∂φ ∂φ ∂φ
∂x ∂y ∂z
2
∂ 2φ
2
∂ 2φ
2
∂ 2φ

∂ φ ∂ φ ∂ φ
= − i+ − j+ − k
∂y∂z ∂z∂y ∂x∂z ∂z∂x ∂x∂y ∂y∂x
=0
∂ 2φ ∂ 2φ
where the C 2 condition is required to show that, for example, =
∂y∂z ∂z∂y
Example: Let F = x2 ey i + x ln zj + (x + z)k. Calculate G = ∇×F and also
∇·G.

i j k
∂ ∂ ∂

∇×F = ∂x = (0 − x/z)i + (0 − 1)j + (ln z − x2 ey )k
∂y ∂z
x2 ey x ln z (x + z)
109
x
Hence G = − i − j + (ln z − x2 ey )k. Calculating the divergence of this we
z
obtain
∂G1 ∂G2 ∂G3
∇·G = + +
∂x ∂y ∂z
1 1
=− +0+ =0
z z
Hence in this example ∇·(∇×F) = 0. Again this is not just chance but is a
general result:
Theorem: Let F be a C 2 vector field, then
div(curlF) = ∇·(∇×F) = 0
Proof:
Let F = F1 i + F2 j + F3 k, then

i j k
∂ ∂ ∂
G = ∇×F = ∂x ∂y

∂z
F1 F2 F3

∂F3 ∂F2 ∂F3 ∂F1 ∂F2 ∂F1
= − i+ − j+ − k
∂y ∂z ∂x ∂z ∂x ∂y
Now divG is given by

∂G1 ∂G2 ∂G3
∇·G = + +
∂x ∂y ∂z
2 2
∂ F3 ∂ F2 ∂ 2 F1 ∂ 2 F3 ∂ 2 F2 ∂ 2 F1
= − + − + −
∂x∂y ∂x∂z ∂y∂z ∂y∂x ∂z∂x ∂z∂y
=0
Where we have again used the C 2 condition to ensure that the mixed second
derivatives commute.
8.1.11 Vector Identities

In this section we collect together some useful vector identities:
(1.) ∇(φψ) = ψ∇φ + φ∇ψ
110
(2.) ∇·(φF) = ∇φ·F + φ∇·F

(3.) ∇×(φF) = ∇φ×F + φ∇×F
(4.) ∇·(F×G) = G·∇×F − F·∇×G
(5.) ∇×(∇φ) = 0
(6.) ∇·(∇×F) = 0
(7.) ∇·(∇φ) = ∇2 φ
(8.) ∇×(∇×F) = ∇(∇·F) − ∇2 F
where ∇2 F = ∇2 F1 i + ∇2 F2 j + ∇2 F3 k.
111
VECTOR CALCULUS INTEGRALS OF VECTOR AND SCALAR FIELDS
8.2 Integrals of vector and scalar fields

8.2.1 Line integrals
Let F(r) be a (smooth) vector field and let γ(t) for t0 ≤ t ≤ t1 be a curve
parameterised by t. Let r(t) denote the position vector of the point on the curve
given by γ(t).
In this section we give a definition of
Z
F·dr
γ
the line integral of F along the curve γ. We will see below that this integral
depends upon the curve γ, but not upon the parameterisation of the curve - only
upon the set of points on the curve.
In the usual way when considering integrals we can break the curve γ up into a
series of m small straight line segments.
Let δrk denote the line segment joining rk−1 to rk , and let Fk = F(rk ). Then we
may define the line integral by
Z Xm
F·dr = lim Fj ·δrj
γ m→∞
j=1
However this definition is not very useful when it comes to calculations. In oder
to calculate a line integral we convert it to an ordinary 1-dimensional integral.
112
We start be looking at the parameterised curve
γ : t 7→ r(t), a≤t≤b
Then
dr
dr = dt
dt
We also look at the value of the vector field along the curve
F = F(r(t))
Combining these the integrand becomes

dr
F·dr = F(r(t))· dt
dt
Finally the integral may be written:
Z Z b
dr
γ t=a) dt

dr
Notice that F(r(t))· is simply some function of t, so that the final integral
dt
Rb
has the form t=a f (t)dt of an ordinary integral.
In order to show that the integral only depends upon the curve γ and not upon the
parameterisation, we introduce a new parameterisation as in section 2.3
Let f : R → R be a monotonic increasing (smooth) map. Let f (c) = a and
f (d) = b. Then we may define a new vector valued function r̃(s) = r(f (s)) for
c ≤ s ≤ d. As s varies between c and d then r̃(s) traces out the same set of
points as r(t), so it is also a parameterisation of the same curve γ. In terms of this
parameterisation the line integral is given by
Z Z d
dr̃
F·dr = F(r̃(s))· ds
γ s=c) ds
However by the chain rule

dr̃ dr dt
=
ds dt ds
113
So changing variable in the integral we get

Z Z d
dr̃
F·dr = F(r̃(s))· ds
γ s=c) ds
Z d
dr dt
= F(r̃(s))· ds
s=c) dt ds
Z b
dr
= F(r(t))· dt
s=a) dt
which agrees with the previous formula. Hence the line integral is independent of
the parameterisation of the curve as claimed.
Example
Z 1a: Let F = 2xyi + (x2 − z 2 )j − 3xz 2 k. Evaluate the line integral
F·dr, where γ is the curve given by
γ
r(t) = 2ti + t3 j + 3t2 k, 0≤t≤1

going from the point (0, 0, 0) to the point (2, 1, 3).
(2,1,3)
z
g2 y
g1
Figure 8.4: Geometry for Examples 1a (γ1 ) and 1b (γ2 ).
From the formula for r(t) above we see that

x = 2t, y = t3 , and z = 3t2
substituting into the formula for F(r) this gives
F(r(t)) = 4t4 i + (4t2 − 9t4 )j − 54t5 k
114
Also differentiating the expression for r(t) we get
dr
= 2i + 3t2 j + 6tk
dt
Hence
Z Z 1
dr
γ t=0 dt
Z 1
= (4t4 i + (4t2 − 9t4 )j − 54t5 k)·(2i + 3t2 j + 6tk)dt
Zt=0
1
= (8t4 + 12t4 − 27t6 − 324t6 )dt
t=0
1
5 351
= 4t −
7 0
323
=−
7
Example 1b: Calculate the line integral for the bf same vector field F = 2xyi +
(x2 − z 2 )j − 3xz 2 k, but this time along the straight line connecting (0, 0, 0) to
(2, 1, 3).
This time the curve is given by
r(t) = 2ti + tj + 3tk, 0≤t≤1
Hence
x = 2t, y = t, and z = 3t
substituting into the formula for F(r) this gives
F(r(t)) = 4t2 i − 5t2 j − 54t3 k
Also differentiating the expression for r(t) we get
dr
= 2i + j + 3k
dt
115
Hence
Z Z 1
dr
γ t=0 dt
Z 1
= (4t2 i − 5t2 j − 54t3 k)·(2i + j + 3k)dt
Zt=0
1
= (6t8 + 36t8 + 12t8 )dt
t=0
1
54 9
= t
9 0
=6
This example shows the important fact that in general the value of the line integral
depends upon the path taken, not just the end points.
Example 2: Let F = ∇φ where φ = xy 2 z. Calculate the line integral of F for
both of the curves used in example 1.
F = ∇φ = y 2 zi + 2xyzj + xy 2 k
For part (a)
r(t) = 2ti + t3 j + 3t2 k, 0≤t≤1
So
F(r(t) = 3t8 i + 12t6 j + 2t7 k
and
dr
= 2i + 3t2 j + 6tk
dt
Hence
Z Z 1
dr
γ t=0 dt
Z 1
= (3t8 i + 12t6 j + 2t7 k)·(2i + 3t2 j + 6tk)dt
Zt=0
1
= (6t8 + 36t8 + 12t8 )dt
t=0
1
54 9
= t
9 0
=6
116
For part (b)

r(t) = 2ti + tj + 3tk, 0≤t≤1
Hence
F(r(t)) = 3t3 i + 12t3 j + 2t3 k
and
dr
= 2i + j + 3k
dt
Hence
Z Z 1
dr
γ t=0 dt
Z 1
= (3t3 i + 12t3 j + 2t3 k)·(2i + j + 3k)dt
Zt=0
1
= (6t3 + 12t3 + 6t3 )dt
t=0
1
24 4
= t
4 0
=6
So that in this example the two line integrals give the same answer. In fact this is
true for the line integral of any vector field that is the form gradφ.
Proposition:
1
Z vector field so that F = ∇φ, for some C scalar field φ, then
If F is a gradient
the line integral F·dr only depends upon the end points A and B of γ, not on
γ
the path taken. In fact
Z
(∇φ)·dr = φ(rB ) − φ(rA )
γ
where φ(rA ) and φ(rB ) are the values of φ at the start and end points of the curve
γ.
Proof: Let γ be parameterised by r(t), t0 ≤ t ≤ t1 then
Z Z t1
dr
(∇φ)·dr = ∇φ(r(t))· dt
γ t=t0 dt
Now
∂φ ∂φ ∂φ
∇φ = i+ j+ k
∂x ∂y ∂z
117
And if we write
r(t) = x(t)i + y(t)j + z(t)k, t0 ≤ t ≤ t1
then
dr dx dy dz
= i+ j+ k
dt dt dt dt
Hence
dr ∂φ dx ∂φ dy ∂φ dz
∇φ· = + +
dt ∂x dt ∂y dt ∂z dt
dφ(r(t))
=
dt
Thus
Z Z t1
dφ(r(t))
(∇φ)·dr = dt
γ t=t0 dt
= [φ(r(t))]tt=t
1
0
= φ(r(t1 )) − φ(r(t0 ))
= φ(rB ) − φ(rA )
8.2.2 Conservative vector fields

Why is the result of the previous section significant? One important application
of line integrals is to find out the work done by a force. For a constant force F
moving an from A to B, the work done is “the force times the distance moved in
the direction of the force”. The mathematical formula for this is
W = F·r, where r = AB
However if the force and path are changing, we have to divide the path into small
sections, and then take the limit
Xm Z
W = lim Fj δrj = F(r)·dr
m→∞ γ
j=1
in order to obtain the formula for the work done by the vector field F(r) in moving
from A to B along the path γ.
If we can find a scalar field such that
F = −∇φ NB note the minus sign
118
then we say that φ is a potential for the force F and the change in φ simply
represents the change in potential energy. In such a situation any loss of potential
energy results in a gain in kinetic energy (and vice-versa) so that the total energy
is conserved. For this reason a gradient vector field is often called a conservative
vector field, and we will use this name from now onwards in this course. In section
4.2 of these notes we will look at conditions which guarantee that a vector field is
conservative.
8.2.3 Surfaces
We start by giving a number of different definitions of a 2-dimensional surface in
three dimensional space R3 . There are essentially three different ways of specify-
ing a surface
(1.) As the graph of a function of 2 variables, z = f (x, y)
(2.) As a level surface of a function of 3 variables F (x, y, z) = c
(3). As a parameterised surface (s, t) 7→ (x(s, t), y(s, t), z(s, t))
For the purpose of calculating surface integrals we will mostly use definition (3)
and think of surfaces as (smooth) maps
S → R3
(s, t) 7→ r(s, t)
where S ⊂ R2 is some 2-dimensional set which specifies the region in which the
parameters s and t lie in. In many of the applications this will simply involve s0 ≤
s ≤ s1 and t0 ≤ t ≤ t1 so that S will be given by the rectangle [s0 , s1 ] × [t0 , t1 ].
If we draw a picture of this map
s1 s2
s3
s4 t5 t6
t4
s5 t3
t2
s6
t1
119
we see that the lines s =const. and t =const. are simply the coordinate lines on
the surface.
Example 1: A cylinder radius a and height h
r(s, t) = a sin si + a cos sj + tk, −π ≤ s < π, 0 ≤ t ≤ h
n̂
The lines t =const. are circles and the lines s =const. are straight vertical lines
on the surface of the cylinder.
Example 2: A sphere radius a centre the origin
r(s, t) = a cos s sin ti + a sin s sin tj + a cos tk, −π ≤ s < π, 0 ≤ t ≤ π
‫ݏ‬ൌͲ
‫ݐ‬ൌͲ
‫ ݐ‬ൌ ߨȀʹ
‫ݐ‬ൌߨ
120
The lines s =const. are the lines of longitude and the lines t =const. are the
lines of latitude on the surface of the sphere. Thus the s-coordinate is just the
φ-coordinate of spherical polar coordinates, while the t-coordinate is just the θ-
coordinate of spherical polar coordinates.
8.2.4 The normal to a surface

We now calculate an expression for the normal to a general surface.
∂r
The vector is tangent to the curve t =const.
∂s
∂r
While the vector is tangent to the curve s =const.
∂t
¶r
¶t
¶r
s1 ¶s
s2
t4
s3 t3
t2
s4
t1
A normal to the surface is therefore given by

∂r ∂r
n= ×
∂s ∂t
We calculate the normal to the surface in the two examples above:
Example 1: A cylinder radius a and height h
r(s, t) = a cos si + a sin sj + tk

∂r
= −a sin si + a cos sj + 0k
∂s
∂r
= 0i + 0j + k
∂t
121
Hence

i j k
n = −a sin s a cos s 0 = a cos si + a sin sj
0 0 1
Note this is simply the (i, j)-part of r(s, t) as one would expect from the picture
above.
Example 2: A sphere radius a centre the origin
r(s, t) = a cos s sin ti + a sin s sin tj + a cos tk
∂r
= −a sin s sin ti + a cos s sin tj + a cot tk
∂s
∂r
= a cos s cos ti + a sin s cos tj − a sin tk
∂t
Hence

i j k

n = −a sin s sin t a cos s sin t a cot t

a cos s cos t a sin s cos t −a sin t
= −a2 sin t(cos s sin ti + sin s sin tj + cos tk)
= −a2 sin tr
Thus n̂ = r̂ (as one would expect!).
8.2.5 Surface integrals

In order to calculate “surface integrals” we first need to find an expression for
the element of surface area dA. To do this we look again at the picture of the
coordinate lines on a surface.
dr2
dA
dr1
s1
s2
t4
s3 t3
t2
s4
t1
122
Suppose we start at the point r(s, t) and move to the nearby point r(s + δs, t).
Then the vector connecting these two points is
δr1 = r(s + δs, t) − r(s, t)
Now by Taylor’s theorem
∂r
r(s + δs, t) = r(s, t) + δs (s, t) + O(δs2 )
∂s
Hence
∂r
δr1 = δs (s, t) + O(δs2 )
∂s
Similarly if we move in the t direction to the nearby point r(s, t + δt) then the
vector measuring the displacement from r(s, t) is
δr2 = r(s, t + δt) − r(s, t)
Again using Taylor’s theorem we have
∂r
r(s, t + δt) = r(s, t) + δt (s, t) + O(δt2 )
∂t
and hence
∂r
δr2 = δt (s, t) + O(δt2 )
∂t
Thus the area of the parallelogram with sides δr1 and δr2 is
δA = |δr1 ×δr2 |

∂r ∂r
≈ × δsδt
∂s ∂s
So in the limit we get

∂r ∂r
dA = × dsdt
∂s ∂s
Example 1: we now use this result to calculate the area of a cylinder. The formula
for a parameterised cylinder of radius a and height h is given by
r(s, t) = a cos si + a sin sj + tk, −π ≤ s < π, 0 ≤ t ≤ h
So that

i j k
∂r ∂r
× = −a sin s a cos s 0
∂s ∂s
0 0 1
= a cos si + a sin sj
123
and hence

∂r ∂r
× = a(cos2 s + sin2 s)1/2
∂s ∂s
=a
Therefore by the above formula for the surface of a cylinder
dA = adsdt
Hence the area is given by

Z π Z h
A= adsdt = 2πah
s=−π t=0
Hence we have shown that the surface area of a cylinder of height h and radius a
is 2πah (in agreement with the standard formula).
Example 2: We now use the formula for dA to calculate the surface area of a
sphere. Recall that a parameterised sphere is given by
r(s, t) = a cos s sin ti + a sin s sin tj + a cos tk, −π ≤ s < π, 0 ≤ t ≤ π
So that

i j k
∂r ∂r
× = −a sin s sin t a cos s sin t a cot t
∂s ∂s
a cos s cos t a sin s cos t −a sin t
= −a2 sin t(cos s sin ti + sin s sin tj + cos tk)
and hence

∂r ∂r
× = a2 sin t cos2 s sin2 t + sin2 s sin2 t + cos2 t 1/2

∂s ∂s
1/2
= a2 sin t (cos2 s + sin2 s) sin2 t + cos2 t
= a2 sin t
Therefore we have shown that for a sphere
dA = a2 sin tdsdt
124
Note that a is the radius of the sphere and (s, t) correspond to the usual spherical
polar coordinates (φ, θ), so that the above formula agrees with the standard result
that
dA = r2 sin θdθdφ
We now use this to calculate the area of a sphere
Z π Z π
A= a2 sin tdsdt
s=−π t=0
Z π
2
= 2πa sin tdt
t=0
= 2πa [− cos t]πt=0
2
= 4πa2
Hence we have shown that the area of a sphere radius a is 4πa2 (in agreement with
the standard formula).
8.2.6 Flux integrals

A flux integral measures the “flow” of a vector field F through a surface S
ZZ
F·n̂dA
S
where n̂ is the unit normal to the surface.

If we think of F as giving the velocity of a fluid, then the flux integral measures
the net rate of flow of the fuid through the surface S.
As with line integrals the method of calculating flux integrals is to use the param-
terisation of the surface to convert them into ordinary integrals.
If we introduce the vector element of area
dS = n̂dA
then we may write the flux integral as

ZZ
F·dS
S
125
Then dS points in the same direction as n and has modulus |dS| = dA. From
this we see that
∂r ∂r
dS = × dsdt
∂s ∂t
This is because the above expression points in the same direction as n, since
∂r ∂r
n= ×
∂s ∂t
and has modulus dA since
∂r ∂r
dA = ×
∂s ∂t
Using this formula for dS we may write the flux integral as
ZZ ZZ
∂r ∂r
F·dS = F(r(s, t))· × dsdt
∂s ∂t
S S
where the integrand

∂r ∂r
F(r(s, t))· ×
∂s ∂t
is simply a function of s and t, so the flux integral has been converted into an
ordinary 2-dimensional integral.
Note: There is an important point to note about this result: the formula we use
to calculate the flux integral depends upon the paramterisation. However one can
show, as with the case of the line integral, that the final answer is independent of
the parameterisation chosen. We do not give the details of the calculation here, but
not that the essential point is that if we calculate dS using a different paramterisa-
tion, (s̃, t̃), it differs from the original one by being multiplied by the the Jacobian
determinant of the transformation from (s, t) to (s̃, t̃). This is exactly the factor
one needs when going from dsdt to ds̃dt̃.
Example: Let S be the parameterised surface given by
r(s, t) = s cos ti + s sin tj + tk, 0 ≤ s ≤ 1, 0 ≤ t ≤ 2π
126
Figure 8.5: Helicoid.
(a twisted ribbon or “helicoid”). Let F(r) ZbeZ the vector field given by F = xi +
yj + (z − 2y)k. Evaluate the flux integral F·dS. We first evaluate the vector
S
field on the surface
r(s, t) = s cos ti + s sin tj + tk
Now on the surface
x = s cos t, y = s sin t, and z = t
So that substituting into the formula F = xi + yj + (z − 2y)k we get
F(r(s, t)) = s cos ti + s sin tj + (t − 2s sin t)k
We now calculate dS. Differentiating the expression for r(s, t) gives
127
∂r
= cos ti + sin tj + 0k
∂s
∂r
= −s sin ti + s cos tj + k
∂t
Now

∂r ∂r
dS = × dsdt
∂s ∂t

i j k
= cos t sin t 0 dsdt
−s sin t s cos t 1
= sin ti + cos tj + (s cos2 t + s sin2 t)k dsdt

= (sin ti − cos tj + sk)dsdt
Hence
ZZ Z 1 Z 2π
F·dS = (s cos ti + s sin tj + (t − 2s sin t)k)·(sin ti − cos tj + sk)dsdt
s=0 t=0
S
Z 1 Z 2π
= (st − 2s2 sin t)dsdt
s=0 t=0
Z 1 2π
1 2 2
= st + 2s cos t ds
s=0 2 t=0
Z 1
= 2sπ 2 ds
s=0
2 2 1
= s π s=0
= π2
Hence ZZ
F·dS = π 2
S
We end this section by giving the expression for a flux integral in the special case
that the surface is the graph of a function z = f (x, y). To calculate the flux integral
in this case we take x and y themselves to be the paramters used to describe the
surface. Hence the parameterised surface is given by
r(x, y) = xi + yj + f (x, y)k, (x, y) ∈ S
128
Then
∂r ∂f
= i + 0j +
∂x ∂x
∂r ∂f
= 0i + j +
∂y ∂y
Hence

∂r ∂r
dS = × dxdy
∂x ∂y

i j k
1 0 ∂f

= ∂x dxdy

∂f
0 1 ∂y

∂f ∂f
= − i− j + k dxdy
∂x ∂y
Thus for the case of a suface which is the graph of a function we have the formula
ZZ ZZ
∂f ∂f
F·dS = F(x, y, f (x, y)· − i − j + k dxdy
∂x ∂y
S
As a final remark in this section we point out that it often saves time to think a bit
about the geometry of a flux integral before going ahead and mechanically doing
the calculation.
2 2 2 2
Example: Let F = xi+yj Z+zk,
Z and S the surface of the sphere x +y +z = a .
Calculate the flux integral F·dS.
S
In this example a little thought shows that F = r and that dS = n̂dA = r̂dA.
Hence
F·dS = r·r̂dA = r(r̂)·(r̂)dA = rdA
Hence ZZ ZZ ZZ
F·dS = rdA = a dA = a.4πa2 = 4πa3
S S S
(since r = a on S).
129
VECTOR CALCULUS INTEGRAL THEOREMS
8.3 Integral Theorems

In this section we look at the relationship between line integrals, flux integrals and
volume integrals using grad, div and curl.
8.3.1 Stokes’s Theorem

We start by looking at the relationship between the line integral round a closed
loop and the flux flowing through the loop.
Theorem (Stokes’s Theorem)
Let S be a bounded, piecewise smooth, orientable, surface in R3 , with boundary
∂S, which consists of a finite number of piecewise C 1 simple closed curves. Let
F be a C 1 vector field whose domain includes S. Then
ZZ I
(∇×F) ·dS = F·dr
∂S
S
I
Note 1: The symbol indicates that the line integral is round a closed curve.
This line integral round ∂S is computed in a clockwise direction relative to the
normal n to S used to define dS. (If the opposite normal is chosen both integrals
therefore change sign).
n̂
¶S
Note 2: In the picture above the boundary to S consists of just one component.
However the theorem allows for the boundary to consist of several components
(e.g. three as shown in the picture below).
130
n̂
¶S2
¶S3
¶S1
In such a case the line integral is the sum of the line integral for the various com-
ponents, so in the example above
I I I I
F·dr = F·dr + F·dr + F·dr
∂S ∂S1 ∂S2 ∂S3
where the direction of integration round all three loops ∂Si is clockwise relative
to the normal n to S.
Note 3: Stokes’s Theorem contains some fine print concerning the nature of the
surface “piecewise smooth, orientable”, the nature of the boundary “finite number
of piecewise C 1 simple closed curves” and concerning the nature of the vector
field F “C 1 on domain that includes S”. All these restrictions are important and
the theorem is not valid without them (or equivalent) conditions. For further ex-
planation on the precise meaning of these terms see the Appendix.
Example: Let S be the paraboloid
z = 9 − x2 − y 2 , z≥0
Verify Stokes’s theorem for the vector field
F = (2z − y)i + (x + z)j + (3x − 2y)k
(0,0,9)
z=9- x2-y2
(0,3,0)
(-3,0,0) (3,0,0)
(0,-3,0)
S: x2+y2=9
131
We first evaluate the flux integral of ∇×F.

i j k
∂ ∂ ∂

∇×F = ∂x
∂y ∂z

2z − y x + z 3x − 2y
= (−2 − 1)i + (2 − 3)j + (1 − −1)k
= −3i − j + 2k
S is the graph of z = f (x, y) = 9 − x2 − y 2 . Also from the diagram the condition

z ≥ 0 is equivalent to x2 + y 2 ≤ 9 so we parameterise S as:
r(x, y) = xi + yj + (9 − x2 − y 2 )k, x2 + y 2 ≤ 9
Then

∂r ∂r
dS = × dxdy
∂x ∂y

i j k
= 1 0 −2x dxdy
0 1 −2y
= (2xi + 2yj + k) dxdy
Hence
ZZ Z Z
(∇×F)·dS = (−3i − j + 2k)·(2xi + 2yj + k)dxdy
S x2 +y 2 ≤9
Z Z
= (−6x − 2y + 2)dxdy
x2 +y 2 ≤9
Z Z
= 2dxdy, (since the other terms cancel by symmetry)
x2 +y 2 ≤9
= 2 × area = 2.9π = 18π
To verify Stokes’s theorem we now calculate the line integral. From the diagram
we see that C = ∂S is a circle radius 3 centre the origin and lying in the xy-plane..
Since the normal n to S point upwards we parameterise C as:
r(t) = 3 cos ti + 3 sin tj + 0k, quad0 ≤ t < 2π
132
Then
F(r(t)) = −3 sin ti + 3 cos tj + (9 cos t − 6 sin t)k
and
dr
= −3 sin ti + 3 cos tj + 0k
dt
Thus
I Z 2π
dr(t)
∂S t=0 dt
Z 2π
= ((−3 sin ti + 3 cos tj + (9 cos t − 6 sin t)k)·(−3 sin ti + 3 cos tj + 0k)) dt
t=0
Z 2π
= (9 sin2 t + 9 cos2 t)dt
Zt=0
2π
= 9dt
t=0
= [9t]2π
t=0 = 18π
Hence in this example we have shown

ZZ I
(∇×F) ·dS = 18π = F·dr
∂S
S
ZZ
2 2
Example 2: Let S̃ be the disc x + y ≤ 9, z = 0. Calculate (∇×F) ·dS
S̃
where F is as before and n is the upward normal.
To do this we again uses Stokes’s theorem
ZZ I
(∇×F) ·dS = F·dr
∂ S̃
S̃
But S and S̃ both have the same boundary given by C, a circle radius 3 in the
xy-plane. So that ∂ S̃ = ∂S. Hence
ZZ I I
(∇×F) ·dS = F·dr = F·dr = 18π
∂ S̃ ∂S
S̃
The method used in this example is often useful so we state it in the general case
as a consequence of Stokes’s theorem.
133
Corollary to Stokes’s Theorem

Let S and S̃ be two bounded, piecewise smooth, orientable, surfaces in R3 , with
the same boundary (consisting of a finite number of piecewise C 1 simple closed
curves). Let F be a C 1 vector field whose domain includes both S and S̃. Then
ZZ ZZ
(∇×F) ·dS = (∇×F) ·dS
S S̃
8.3.2 Conservative Vector Fields

We start by defining a vector field with path-independent integral.
Definition: A continuous vector field F has path-independent line integral in the
region U if Z Z
F·dr = F·dr
γ1 γ2
for any 2 piecewise C 1 , oriented curves γ1 and γ2 which lie in the set U and both
start at the same point A, and end at both end at the same point B.
Proposition: Let F be a continuous vector field. Then F has path-independent
integrals in U if and only if I
F·dr = 0
γ
for all simple closed piecewise C 1 curves γ in U .
B
g1
g2
g3
g4
A
Proof: Suppose that F is a path-independent vector field and γ is a closed curve.

Then we may split it up into two curves γ1 and γ2 that start at A and end at B. Then
134
γ is just γ1 followed by −γ2 (ie γ2 traversed in the opposite direction). Hence

I Z Z
F·dr = F·dr − F·dr = 0
γ γ1 γ2
Conversely suppose the line integral round any closed loop vanishes and γ1 and
γ2 are two curves that start and finish at the same point. Then γ1 followed by −γ2
(ie γ2 traversed in the opposite direction) is just γ so that
Z Z I
F·dr − F·dr = F·dr = 0
γ1 γ2 γ
We now give a theorem that relates path-independent vector fields to conservative

(or gradient) vector fields.
Theorem:
Let F be a continuous vector field in the connected open region U of R3 . Then
there exists a C 1 scalar field φ such that F = ∇φ if and only if F has path
independent integral in U .
Proof: We have already shown that if F = ∇φ then
Z Z
F·dr = (∇φ)·dr = φ(rB ) − φ(rA )
γAB γAB
and hence the integral is independent of the path.

If F has path independent integral then we may define a scalar field φ by
Z
φ(r) = F·dr
γ
where γ is some path from some (arbitrary) fixed point A to the point P with po-
sition vector r. Note that because the vector field F has path independent integral
this is well defined without having to specify the path γ.
We now show by direct computation that
∂φ ∂φ ∂φ
= F1 , = F1 , = F1 , ie ∇φ = F.
∂x ∂x ∂x
Let Z
φ(x, y, z) = F·dr
γ
Then Z Z
φ(x + h, y, z) = F·dr + F·dr
γ γ1
135
where γ1 is a straight line parallel to the x-axis going from x to x + h. Hence

Z x+h
φ(x + h, y, z) − φ(x, y, z) = F·dr
x
Z x+h
= F1 dx(since dr = idx on γ1 )
x
= hF1 (x, y, z) + O(h2 )
Hence

∂φ φ(x + h, y, z) − φ(x, y, z)
= lim = lim (F1 (x, y, z) + O(h)) = F1 (x, y, z)
∂x h→0 h h→0
Similarly
∂φ ∂φ
= F2 , and = F3
∂y ∂z
Note that different choices of the arbitrary initial point A simply change φ by a
constant, but this has no effect on ∇φ which is unchanged.
The above theorem tells us that a conservative vector field is equivalent to one
with path-independent integral. Unfortunately this condition is not very easy to
check as we need to know it is true for any possible path! The following theorem
gives a more useful criterion, which is easy to check.
Theorem: Let F be a C 1 vector field defined on the simply connected region U
in R3 . Then
There exists some C 2 scalar field φ such that F = ∇ if and only if ∇×F = 0
Note: The fine print in this theorem is again important to note. In particular the
simply connected condition (which means we can continuously shrink any closed
path down to a point) is required - see the exercise sheets for a counterexample to
this theorem when U is not simply connected.
Proof:
We know from chapter 2 that ∇×(∇φ) = 0, so if F = ∇φ then we know that
∇×F = 0.
On the other hand if ∇×F = 0 then we know from Stokes’s theorem that
I ZZ
F·dr = (∇×F)·dS = 0
C
S
136
where S is some spanning surface with boundary C. (Note that such a spanning
surface S exists precisely because of the simply-connected condition.)
Hence by the first proposition in this section F has path-independent integral, and
thus by the previous theorem F = ∇φ, as required.
8.3.3 Finding potentials for conservative vector fields

In the previous section we showed how to test for a conservative vector field in
this section we show how one can actually find a potential for the vector field.
Example: F = (2xy+cos 2y)i+(x2 +2y−2x sin 2y)j, show that F is conservative
and find a scalar field φ such that ∇φ = F.
To show that F is conservative we compute ∇×F.

i j k
∂ ∂ ∂

∇×F = ∂x ∂y ∂z
2xy + cos 2y x2 + 2y − 2x sin 2y 0
= 0i + 0j + (2x − 2 sin 2y − 2y + 2 sin 2y)k
=0
Hence by the theorem of the previous section F is conservative. We can therefore

try and find a scalar field such that ∇φ = F. Since F only depends upon the
x and y coordinates and has no k terms, in this example we may assume that
φ = φ(x, y). Hence we want to find a scalar field φ such that
∂φ
= F1 = 2xy + cos 2y (8.1)
∂x
∂φ
= F2 = x2 + 2y − 2x sin 2y (8.2)
∂y
Integrating (1) we get
φ = x2 y + x cos 2y + f (y) (where f (y) is a “constant of integration”) (8.3)
Differentiating (3) with respect to y we get

∂φ
= x2 − 2x sin 2y = f 0 (y) (8.4)
∂y
Comparing (2) and (4) we see that
f 0 (y) = 2y ⇒ f (y) = y 2 + c
137
Substituting for f (y) in (3) gives
φ(x, y) = x2 y + x cos 2y + y 2 + c
We end this section by summarising the four equivalent ways of characterising a

C 1 conservative vector field (in a simply connected region U ).
(1.) Path-independent integral
Z Z
F·dr = F·dr
γ1 γ2
For all paths γi starting at A and ending at B

(2.) Vanishing integral round closed curve
I
F·dr = 0
γ
for all line integral round a simple closed curve γ

(3.) Gradient vector field
F = ∇φ
(4.) Curl free vector field

∇×F = 0
8.3.4 Stokes’s theorem in the plane

Stokes’s theorem states that
ZZ I
(∇×F) ·dS = F·dr
∂S
S
In this section we look at this in the case that the vector field F depends only
upon x and y and both it and the surface S lie in the xy-plane. Since F lies in the
xy-plane and only depends upon x and y, then
F = F1 (x, y)i + F2 (x, y)j
138
Then

i j k
∂ ∂ ∂

∇×F = ∂x ∂y ∂z
F1 (x, y) F2 (x, y) 0

∂F2 ∂F1
= 0i + 0j + − k
∂x ∂y

∂F2 ∂F1
= − k
∂x ∂y
Also for a surface lying in the xy-plane
dS = n̂dA = kdxdy
Hence ZZ ZZ
∂F2 ∂F1
(∇×F) ·dS = − dxdy
∂x ∂y
S S
On the other hand if we only consider points in the xy-plane then
r = xi + yj ⇒ dr = dxi + dyj
So that I I
F·dr = (F1 dx + F2 dy)
∂S ∂S
Hence in the xy-plane Stokes’s theorem becomes
ZZ I
∂F2 ∂F1
− dxdy = (F1 dx + F2 dy)
∂x ∂y ∂S
S
this result is known as Green’s theorem:

Green’s Theorem Let S be a bounded, piecewise smooth, orientable, surface in
R2 , with boundary ∂S, which consists of a finite number of piecewise C 1 simple
closed curves. Let F be a C 1 vector field in R2 whose domain includes S. Then
ZZ I
∂F2 ∂F1
∂x ∂y ∂S
S
We have seen that one can obtain Green’s theorem from Stokes’s theorem by
restricting things to be in a 2-dimensional space. However it is possible to derive
139
Stokes’ theorem from Green’s theorem. This follows from the fact that if one
parameterises the surface S in Stokes’s theorem in terms of s and t by
r(s, t), with (s, t) ∈ S̃
then one finds

ZZ ZZ
∂Ft ∂Fs
(∇×F) ·dS = − dudv
∂s ∂t
S S̃
∂r ∂r
where Fs = F· and Ft = F· . It then follows from Green’s theorem in the
∂s ∂t
st-plane that
ZZ I
∂Ft ∂Fs
− dudv = (Fs ds + Ft dt)
∂s ∂t ∂ S̃
S̃
On the other hand

∂r ∂r
dr = dr + ds
∂r ∂s
So that F·dr = Fr dr + Ft dt and hence
I I
(Fs ds + Ft dt) = F·dr
∂ S̃ ∂S
Putting all this together we get

ZZ ZZ
∂Ft ∂Fs
(∇×F) ·dS = − dudv
∂s ∂t
S
I S̃
= (Fs ds + Ft dt) (by Green’s theorem)
I ∂ S̃
= F·dr
∂S
Hence we have used Green’s theorem in the st-plane to prove Stokes’s theorem:
ZZ I
(∇×F) ·dS = F·dr
∂S
S
140
8.3.5 The proof of Green’s theorem

In this section we prove Green’s theorem (you do not require to know this proof
for the exam). As we have seen this also allows us to prove Stokes’s theorem.
Outline proof of Greens’ Theorem
ZZ I
∂F2 ∂F1
∂x ∂y ∂S
S
(x,y+)
C
+
¶S=C+-C-
C-
(x,y-)
Let C + be the upper curve and C − be the lower curve in the above diagram. Let
(x, y + (x)) denote points on the graph of the upper curve and (x, y − (x)) denote
points on the graph of the lower curve. Then if we traverse the closed curve in a
clockwise direction we move from left to right on the bottom curve and from right
to left on the upper curve. Hence
Z Z b
F1 dx = F1 (x, y − (x))dx
−
ZC Zx=a
a
F1 dx = F1 (x, y + (x))dx
C+ x=b
Z b
=− F1 (x, y + (x))dx
x=a
141
Thus
I Z Z
F1 dx = F1 dx + F1 dx
C C− C+
Z b
F1 (x, y + (x)) − F1 (x, y − (x)) dx

=−
x=a
Zb
y + (x)
=− [F1 (x, y)]y=y− (x) dx
x=a
Z (Z y+ (x) )
∂F1
=− dy dx
y=y − (x) ∂y
Hence I ZZ
∂F1
F1 dx = − dxdy
C ∂y
S
Similarly I ZZ
∂F2
F2 dy = dxdy
C ∂x
S
Hence I I ZZ
∂F2 ∂F1
F1 dx + F2 dy = − dxdy dxdy
C C ∂x ∂y
S
as required.
8.3.6 The divergence theorem

Stokes’s theorem relates a flux integral of the curl of a vector field over a surface
S to a line integral round the boundary ∂S. In this section we consider an integral
theorem which relates a volume integral of the divergence of a vector field over a
region V to a flux integral over the boundary surface ∂V . This theorem is called
Gauss’s theorem or the divergence theorem.
Theorem (Gauss’s theorem or the divergence theorem)
Let V be a bounded solid region in R3 with boundary surface ∂V which consists
of a finite number of piecewise smooth, closed orientable surfaces (with the ori-
entation chosen so that the normal point out of the surface). Let F be a C 1 vector
field then y {
(∇·F)dv = F·dS
V ∂V
142
{
where the symbol is used to indicate that the integral is over a closed surface.
S1 a
S3
h
S2
Example 1: Let F be the vector field F = xi + yj + zk. Let S be the closed

surface of the cylinder radius a and height h (with the centre of the base at the
origin). Calculate the flux integral
{
F = F·dS
S
The hard way is to calculate the flux integral directly. To do this we observe that
S consists of three surfaces: the top S1 , the bottom S2 and the curved surface S3
so that { x x x
F·dS = F·dS + F·dS + F·dS
S S1 S2 S3
On S1 we have dS = kdxdy so that F·dS = zdxdy

x x
F·dS = hdxdy = πa2 h
S1 S1
On S2 we have dS = −kdxdy so that F·dS = −zdxdy

x x
F·dS = 0dxdy = 0
S2 S2
On S3 we have dS = (xi + yj)(x2 + y 2 )−1/2

x x x
F·dS = (x2 + y 2 )1/2 dA == adA = a.Area of S3 = 2πa2 h
S3 S3 S3
143
Hence {
F = F·dS = πa2 h + 0 + 2πa2 h = 3πa2 h
S
The easy way is to use the divergence theorem
{ y
F = F·dS = (∇·F)dv
S V
Now
∂F1 ∂F2 ∂F3
∇·F = + +
∂x ∂y ∂z
=1+1+1=3
Hence
{ y y
F = F·dS = (∇·F)dv = 3dv = 3.volume of V = 3πa2 h
S V V
In fact this is a general result

Proposition
Let S be a closed surface then
{
r·dS = 3 × volume of region enclosed by S
S
Proof: ∇·r = 3 so bt the divergence theorem

{ y
r·dS = 3dv = 3 × volume of region enclosed by S
S V
x
Example 2: Use the divergence theorem to calculate F·dS where F = 2xy 2 i+
S
z 3 j − x2 yk and S is the hemisphere x2 + y 2 + z 2 = a2 , z ≥ 0.
(0,0,a)
(0,a,0)
~
S
(-a,0,0) (a,0,0)
(0,-a,0) ~
S : x2 + y2 £ a2; z = 0
S : x2 + y2 + z 2 = a2
144
Since S is not a closed surface we cannot immediately apply the divergence theo-
rem. However if we add on the disc S̃ given by x2 +y 2 ≤ a2 , z = 0 (see diagram),
then S ∪ S̃ is a closed surface. Hence
x x y
F·dS + F·dS = (∇·F)dv (∗)
S S̃ V
where V is the region enclosed by S and S̃.

We start by calculating the integral on the RHS.
∂F1 ∂F2 ∂F3
∇·F = + + = 2y 2 + 0 + 0 = 2y 2
∂x ∂y ∂z
Hence y y
(∇·F)dv = 2y 2 dydxdz
V V
If we do the integral with respect to x and z first we integrate over a semicircular
slice of radius (a2 − y 2 )1/2 , which has area 21 π(a2 − y 2 )
y Z a
1
2
2y dydxdz = 2y 2 π(a2 − y 2 )dy
−a 2
V
Z a
=π (a2 y 2 − y 4 )dy
−a
4π 5
= a
15
Now x x x
F·dS = F·kdA = x2 ydxdy = 0 (by symmetry)
S̃ S̃ S̃
Hence substituting into (*) gives
x x y
F·dS + F·dS = (∇·F)dv
S S̃ V
x 4π 5
F·dS + 0 = a
15
S
Hence x 4π 5
F·dS = a
15
S
145
8.3.7 Vector potential

For a C 2 scalar field we know that ∇×(∇φ) = 0 so that if F = ∇φ then
∇×F = 0. However we have also seen that in a simply connected region the
converse is true and if ∇×F = 0 there exists a scalar field φ such that F = ∇φ.
In this section we examine the consequences of the identity
∇·(∇×F) = 0
It follows from this that
G = ∇×F ⇒ ∇·G = 0
However in a simply connected region one can also show (although we do not
give the proof in this course) that
∇·G = 0 ⇒ there exists a vector field F such that G = ∇×F
We call such a vector field F a vector potential for G.

Note that if G is divergence free then it follows from the divergence theorem that
{ y y
G·dS = (∇·G)dv = 0dv = 0
S V V
Finally the above result shows that if G is divergence free then

x x
G·dS = G·dS
S1 S2
where S1 and S2 are two open surfaces with the same boundary C = ∂Si . Note,
that we have already come across this result when we considered Stokes’s theorem
in section 4.1.
146

LectureNotesMATH2048 PDF

Uploaded by

Copyright:

Available Formats

LectureNotesMATH2048 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

LectureNotesMATH2048 PDF

Uploaded by

Copyright:

Available Formats

MATH2048

Mathematics for Engineering

1 Ordinary Differential Equations 4

Ordinary Differential Equations

1.1 Second order equations

The Wronskian determinant

W (x) = W [y1 , y2 ](x) (1.7a)

So that we may write r

Homogeneous constant coefficient

3. b2 − 4ac = 0: Repeated root λ giving one solution eλx . Direct calculation

y = c1 xk1 + c2 xk2 . (1.22)

Reduction of order and the Wronskian

y20 = v 0 y1 + vy10 (1.25)

(Example using Legendre polynomials)

1.1.2 Inhomogeneous Linear Second Order ODEs

where y1 (x) and y2 (x) are solutions of the homogeneous equation

Finding a particular solution: undetermined coefficients

yp0 = (v1 y10 + v2 y20 ). (1.37)

We can then differentiate again to obtain

yp00 = (v1 y100 + v2 y200 ) + (v10 y10 + v20 y20 ). (1.38)

r(x) = yp00 + pyp0 + qyp (1.39a)

v10 y1 + v20 y2 = 0, (1.40a)

1.2 Boundary and Initial Value Problems

This is an example of a boundary value problem. More generally we can specify

a0 y(x0 ) + b0 y 0 (x0 ) = α, (1.47a)

1.2.1 Eigenvalues and functions

y(x) = a1 eαx + a2 e−αx . (1.49)

The only solution compatible with the boundary conditions (1.48b-1.48c) is

3. λ √= µ2 > 0. In this case the auxiliary equation has imaginary roots

y(x) = c1 sin (µx) + c2 cos (µx) . (1.51)

y(x) = c1 sin (µx) . (1.52)

In order to satisfy the final boundary condition, equation (1.48c), we either

sin (µπ) = 0. (1.53)

Therefore the non-trivial solutions of the BVP (1.48) are

yn (x) = sin (nπx) . (1.55)

1. The eigenvalues are all real.

2. The different eigenfunctions are orthogonal in this sense:

3. The eigenvalues form an increasing sequence of positive numbers that ap-

1.2.2 Sturm-Liouville form

This is often written as

A solution yn of the regular Sturm-Liouville problem given by equation (1.60)

Self-adjoint operators and boundary conditions

For the Sturm-Liouville equation (1.60) we find that

We want this to vanish in order to have a self-adjoint problem. We therefore have

This can be guaranteed if at x = x0 and x = x1 we insist on some appropriate

4. p = 0: x = x0 , x1 are singular points of the ODE.

y(x0 ) = y(x1 ), y 0 (x0 ) = y 0 (x1 ), p(x0 ) = p(x1 ). (1.66)

Lȳ = λ̄w(x)ȳ. (1.68)

y ȳ = |y|2 > 0, (1.70)

Properties 2. The eigenfunctions yn are real.

L (fn + jgn ) = L (fn ) + jL (gn ) (1.71a)

By equating real and imaginary parts in (1.71) it follows that

L (fn ) = −λw(x)fn , (1.72a)

Therefore, as the eigenvalues are distinct, we have that

We say that the associated eigenfunctions ym , yn are orthogonal with respect to

becomes arbitrarily small as N → ∞.

2.1 Fourier series

Consider the integral

2 sin A sin B = cos(A − B) − cos(A + B)