AMME2000 & BMET2960/9960:Engineering Analysis
B. Thornber
June 17, 2021
1
Contents
1 Introduction to Numerical Methods 5
1.1 Why do we need numerical methods? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Interpolation, Numerical Integration and Differentiation . . . . . . . . . . . . . . . . . . . 6
1.2.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Taylor Series Derivation of the Central Difference Scheme . . . . . . . . . . . . . . . . . 8
2 Partial Differential Equations 10
3 Heat Equation 11
3.1 Introduction and Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Expected Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Use of the Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4 Simple Inhomogeneous Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.5 More Complex Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.5.1 Step 1: Steady State solution T ss(x) . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.5.2 Step 2: Formulation and Solution of the Homogeneous Problem . . . . . . . . . . 15
3.5.3 Step 3: Fourier Series for non-integer modes . . . . . . . . . . . . . . . . . . . . 16
3.5.4 Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4 Key Concepts in Numerical Methods 18
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Order Of Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4 Consistency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.5 Verification and Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.6 Quantifying the Errors with Error Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.7 Computing the Observed Order of Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . 19
5 Numerical Solution of the Heat Equation 21
5.1 Forward in Time, Central in Space (FTCS) . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.1.1 Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.1.2 von Neumann Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2
5.2 Backward in Time, Central in Space Scheme (BTCS) . . . . . . . . . . . . . . . . . . . . 24
5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6 The Wave Equation 26
6.1 Derivation of the Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.2 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.3 Physical Behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
7 Numerical Solution of the Wave Equation 29
7.1 Second Order Accurate Central Difference in Time and Space . . . . . . . . . . . . . . . 29
7.2 von Neumann Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
7.3 Advanced - Special Case of a Uniform Grid . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
8 Laplace and Poisson Equations 33
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
8.2 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
9 Numerical Solution of the Laplace and Poisson Equation 38
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
9.2 Second Order Central Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
9.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
10 Understanding PDEs 41
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
11 Fourier Integrals, Fourier Transform and Laplace Transform 43
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
11.2 Fourier Integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
11.2.1 Rationale and Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
11.3 Fourier Integrals in the Solution of PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
11.4 The Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
11.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
11.4.2 Formalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
11.4.3 Application to Signal Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3
11.4.4 Important Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
11.4.5 The Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
11.4.6 The Fourier Transform Solution to the Heat Equation . . . . . . . . . . . . . . . 51
11.4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
11.5 Laplace Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
11.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
11.5.2 A motivating example: application to ODEs . . . . . . . . . . . . . . . . . . . . . 53
11.6 Application to PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
11.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
12 Finite Element Analysis 58
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
12.2 Formalisation for General ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
12.3 Does this mean that the Finite Element Method has Zero Error? . . . . . . . . . . . . . . . 69
12.4 Specialisation of FEA to Structural Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 69
12.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
12.4.2 Predicting Stresses in a Simple Beam under Compression . . . . . . . . . . . . . 70
12.4.3 Predicting Stresses in a Prosthetic Leg . . . . . . . . . . . . . . . . . . . . . . . . 71
12.4.4 Predicting Stresses in a Helicopter Blade . . . . . . . . . . . . . . . . . . . . . . 74
12.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
13 Course Summary 79
A Derivation of the Heat Equation 81
B Derivation of the Fourier Sine Series Coefficients 82
4
Introduction
These outlines should be read in conjunction with the lecture slides and online videos for AMME2000 &
BMET2960.
1 Introduction to Numerical Methods
1.1 Why do we need numerical methods?
Numerical methods have been developed in order to give solutions to problems where there is no analytical
solution. An analytical solution is the answer to a given problem which when evaluated is exact. An
analytical solution may be very complex and require a computer to solve, but is exact up to the precision of
the method used to evaluate it - on a computer this is typically 8, 16 or 32 significant figures (single, double
or quadruple precision). An example of an analytical solution could be the derivative of f (x) = sin(x),
which is:
d f (x)
dx
= cos(x) (1)
There are a wide range of problems which do not have analytical solutions typically due to (i) geometric
complexity, (ii) nonlinearity in the governing equations which prevents analytical solutions or (iii) presence
of multiple physical effects which may not be easily separated. In this case we must resort to a numerical
solution. A numerical solution is an approximate solution to the problem. However, a carefully designed
numerical solution with sufficient computational effort applied can approach the level of accuracy required
to enable engineering analysis. Going back to our previous example, a numerical solution for the derivative
of f (x) = sin(x) may be computed from two evaluations of f (x) as follows:
d f (x)
dx
=
sin(x+∆x)− sin(x−∆x)
2∆x
+ ε (2)
where ε is the error made in the approximation, which would be a function of the interval ∆x between the
actual point of interest and the two sample points used to estimate the gradient. One can envisage that if ∆x
is made very small, the result would be very close to exact. Despite this the solution is not exact - it is only
exact in the limit of ∆x→ 0 which is never actually achieved. Despite this limitation, clearly such a small
value of ∆x could be chosen so that the numerical solution is sufficient for any engineering application.
Numerical methods are enabled by numerical approximations - these are approximate estimates of the real
quantities.
A word of caution
Numerical solutions are enabled by powerful computers. However, you must always be aware that the
possibility of an error when you write your first numerical solution is very high. It is very rare to write a
numerical solution, implement it, and for it to be perfect immediately. Treat your numerical results with
suspicion - they are wrong until proven to be correct. The best way to prove that they are correct is by
reproducing problems for which there is an analytical solution - demonstrating that the method matches
the exact solution. This is a role of huge importance for analytical solutions, such that in some fields a
substantial prize is offered for analytical solutions of certain governing equations (e.g. US$1m for the
Navier-Stokes equations). Make sure that any numerical solution you make is verified and validated.
5
1.2 Interpolation, Numerical Integration and Differentiation
Typical fundamental applications of numerical methods are to integrate or differentiate data which may
be presented only as a series point measurements f (xi) where xi are a series of discrete points. If f (xi)
represents the measured flow rate at different positions across a river, then the integral may represent the
flow rate for example. If f (ti) is a position measurement as a function of time, then the derivative is the
velocity.
Here we will summarise several approaches to interpolation, and extend this to numerical integration.
1.2.1 Interpolation
Lagrange Polynomials
Given a collection of points as shown on slide 21, the most straightforward method to integrate the data is
to simply assume that the values of the data are constant in the region x−∆x/2 < x < x+∆x/2, where ∆x
is the distance between the points. A second approach is to draw a straight line between each pair of points.
For an arbitrary number of points, a curve of arbitrary order can be defined which passes through all of the
points, called a Lagrange polynomial. This polynomial is defined as follows:
pN(x) = L0 f0+L1 f1+L2 f2+ ...LN fN =
N
∑
0
Li fi, (3)
where the number of individual points ranges from 0 to N, fi = f (xi) where xi is the position of point i. We
want the polynomial pN(x) to pass through all of the points, and this occurs if the function Li fi = fi at xi
and Li fi = 0 for all other values of i. With this as a constraint, the polynomials Li are given by:
Li = ∏
0≤m≤N,m6=i
(x− xm)
(xi− xm) (4)
where ∏ means that each term is multiplied by the next term. As an example, the quadratic Lagrange
polynomial:
p2(x) = L0 f0+L1 f1+L2 f2, (5)
L0 =
(x− x1)
(x0− x1)
(x− x2)
(x0− x2) (6)
L1 =
(x− x0)
(x1− x0)
(x− x2)
(x1− x2) (7)
L2 =
(x− x0)
(x2− x0)
(x− x1)
(x2− x1) (8)
Unfortunately, Carl Runge discovered a problem extending this method to very high order polynomials for
a large number of points - for some functions the error of the polynomial actually increases with the order of
the polynomial - i.e. the more points used in the Lagrange polynomial the error may increase! This counter
intuitive result led to a search for a better alternative. One possible option is to split the data up and fit it
with several lower order polynomials, however this has the disadvantages that it produces non-continuous
derivatives where the polynomials meet. This severely restricts the utility of this approach.
Splines
6
Splines were developed to address this. The main idea is to fit the data with curves which not only meet at
the same value at the joining points, but where the first or more derivatives of the function are also equal
at the joining points. The most common method is the cubic spline, where individual cubic curves are fit
between each pair of points, and at the polynomials are chosen to have the same first and second derivatives
at the junction between the two polynomials.
This leads to a more complicated fitting procedure than the Lagrange polynomials, in that each polynomial
depends on the subsequent one, and so on. However, it can be formulated in a very computationally efficient
manner. It can also be shown that the errors decrease smoothly as the number of points increase, thus
avoiding the problem with the Lagrange polynomials. This has found enormous application, particularly in
computer aided design.
Noisy Data
Some data is inherently noisy due to measurement uncertainty. A high order polynomial fit would make no
sense in that case, it would ’overfit’ the data, i.e. give a level of fit which would oscillate substantially and
not give a meaningful trend.
This is where Legendre’s Method of Least Squares is ideal. It is usually used to fit a lower order polynomial
to a given set of points, by minimising the mean square error between the polynomial value at each of
the points and the measurement data. Once more, this process can be implemented in a very efficient
computational algorithm for many thousands of points. This is usually a sound method for extracting
trends, but may be poor for derivatives.
1.2.2 Numerical Integration
A first application of interpolation is numerical integration. The first step in numerical integration is to
fill in the gaps between the discrete points. This can be done by any of the above approaches, although
obviously we would caution against using very high order Lagrange polynomials. The basic approach is to
fit the data with a curve, then integrate the area under each individual curve. For example, if a straight line
is fit between each pair of points, the numerical integration is the sum of the trapeziums formed by the area
under each pair of points. This simple approach is the Trapezoidal rule:
∫ x+∆x
x
f (x)dx≈ ∆x f (x+∆x)+ f (x)
2
(9)
This represents a single integration of the area under a single straight line joining points (a, f (a)), (b, f (b)).
This can then be extended to many points by a summation of many smaller areas, which gives:
∫ xN
x0
f (x)dx≈
N−1
∑
i=0
f (xi+1)+ f (xi)
2
∆xi, (10)
where ∆xi = xi+1− xi.
This is a second order accurate scheme, basically meaning if you increase the sampling points by a factor
of 2 (for the same integration length, so halving the spacing), error will decrease by 22 (we will talk about
what this means later). However, a very popular fifth order scheme is given by fitting a quadratic Lagrangian
polynomial to groups of three points, it is called Simpson’s rule:
∫ x+∆x
x−∆x
f (x)dx≈ ∆x
3
[ f (x−∆x)+4 f (x)+ f (x+∆x)] (11)
7
This method, when applied to multiple groups of 3 points, leads to a fifth order accurate result - i.e. if the
spacing between points decreases by a factor 2 then the error decreases by a factor of 25. This is a huge
improvement on the Trapezoidal rule without an enormous increase in effort to evaluate the integral.
A quick note on computational effort
Different operations take a different amount of time for a computer to evaluate. Additions, subtractions
and multiplications are very fast. Divisions are slow, typically ≈ 7 times slower than the other three. More
complex functions such as ln, cos, (.)r where r is a real number are even slower. A first strategy is to
rearrange and simplify the formulae to remove as many evaluations of the costly functions as possible. For
example, if you have to divide by ∆x multiple times, simply compute it once and store the value of 1/∆x
which you then multiply instead of divide (i.e. one divide and two multiplies is already faster than two
divides).
1.3 Taylor Series Derivation of the Central Difference Scheme
The definition of Taylor Series expansion at point x+∆x is:
f (x+∆x) = f (x)+
∂ f
∂x
|x=a(∆x)+ 12!
∂ 2 f
∂x2
|x=a(∆x)2+ 13!
∂ 3 f
∂x3
|x=a(∆x)3+ · · · · · ·+ 1n!
∂ n f
∂xn
|x=a(∆x)n (12)
Note that ∆x may be negative, in which case the odd powers will be negative. Assuming that central
difference approximation has the form:
∂ f (i)
∂x
|i= A f (i−1)+B f (i)+C f (i+1) (13)
The points used in the numeical estimation of the derivative is called ‘the stencil’. Here the stencil is
the collection of three points, points i−1, i and i+1.
Expanding each term in the right hand side at x = i point by Taylor Series to third order accuracy:
A f (i−1) = A f (i)−A∆x∂ f
∂x
|i+A(∆x)
2
2!
∂ 2 f
∂x2
|i−A(∆x)
3
3!
∂ 3 f
∂x3
|i+ · · · · · ·
B f (i) = B f (i)
C f (i+1) =C f (i)+C∆x
∂ f
∂x
|i+C (∆x)
2
2!
∂ 2 f
∂x2
|i+C (∆x)
3
3!
∂ 3 f
∂x3
|i+ · · · · · ·
(14)
Summing each term in the left and right hand side in the equation set (14), we now want to match the
left hand side of equation (13). We note that we have 3 coefficient, A, B and C, thus we can use these
coefficients to set at most three of the terms in the right hand side of equation (14) to zero. A neat property
of the Taylor series is that the first term in each expansion in equation (14) is the only term which contains
f (i), the second term is the only one with ∂ f/∂x etc. Thus, by inspection, we can match the sum of each
grouping of terms (columns in equation (14)) to the right hand side of equation (13) to give the following
set of simultaneous equations for A, B and C:
A+B+C = 0
(there are no instances of f (i) in the right hand side of equation (13)
8
−A∆x+C∆x = 1
(the right hand side of equation (13 is equal to ∂ f/∂x)
A+C = 0
(there are no instances of ∂ 2 f/∂x2 in the right hand side of equation (13)
The solution of equations are:
A =− 12∆x
B = 0
C = 12∆x
Substituting the value of A, B and C into equations set (16) and summing each term in the left and right
hand side, then obtain:
∂ f (i)
∂x
=
f (i+1)− f (i−1)
2∆x
− ∆x
2
6
∂ 3 f
∂x3
+ · · · · · · (15)
The order of accuracy is the leading order error term in ∆x of the approximation to the derivative, i.e. ∆xn
is nth order. Taking the first stencil and including the term of order:
∂ f (i)
∂x
=
f (i+1)− f (i−1)
2∆x
−O(∆x2) (16)
Thus central difference scheme is a second order accurate representation of the first derivative.
9
2 Partial Differential Equations
We are used to ordinary differential equations, such as d
2x
dt2 = −g, which describes position with respect
to time under acceleration by gravity. To solve this equation we simply integrate twice, and apply our
knowledge of the initial position and velocity to fix the constants of integration.
However, the vast bulk of interesting problems in engineering depend on more than one dimension, and
thus naturally lead us to have to solve differential equations which have derivatives with respect to many
variables - these are partial differential equations (PDEs). The notation is slightly different, in that a partial
derivative uses a ∂/∂x instead of d/dx to represent a derivative.
Similar to the solution of the ordinary differential equation (ODE), additional constraints are required to
determine the solution. However, unlike ODEs, the functional form of the solution on the variables is not
determined beforehand. For example, in the previous ODE the solution has the functional form x(t) =
−gt2
/ 2+Bt+C, with B and C constants determined by the initial conditions. For a PDE, both the functional
form (dependence on e.g. x, t) and the constants need to be determined from boundary conditions for steady
state problems, with the additional need for initial conditions in unsteady problems.
Thus PDEs are often denoted ‘Initial Value Problems’ or ‘Boundary Value Problems’, or both IVP and BVP.
There are many different types of equations, yet having a qualitative understanding of the behaviour of
a given equation is incredible useful. This is done via classification into one of three classes covered in
this course - and when broadened out in the rest of your degree these classifications cover a huge range of
different problems in Engineering and Physics.
The first are hyperbolic equations - these typically govern the behaviour of waves. These waves can travel
in one direction (e.g. a concentration of a pollutant being convected downstream in a river), or multiple
directions such as a wave expanding when a stone is dropped into a lake. These equations are typically
dependent on space and time.
The second class is parabolic - these equations represent the diffusion of a quantity such as heat, species in
space as time proceeds.
The final class is Elliptic - these equations typically represent a steady state (t → ∞, ∂/∂ t → 0) diffusion
problem. As it is steady state, it does not matter what the initial conditions were and so Elliptic problems are
determined entirely by their boundary conditions, hence being BVPs. Parabolic and Hyperbolic equations
are usually IVP and BVPs.
Boundary conditions come in two flavours in this course, Dirichlet boundary conditions take a fixed value
of the dependent variable, such as a specified Temperature in a heat diffusion problem, at the boundary.
Neumann boundary conditions are a fixed derivative of the dependent variable, such as a fixed ∂T/∂x in a
heat diffusion problem. Different boundary conditions lead to different solutions of the PDE.
10
3 Heat Equation
3.1 Introduction and Analytical Solution
The generic form of the heat equation is:
ut = c2uxx (17)
The heat equation is a Parabolic equation which describes the diffusion of a quantity in space and time.
Diffusion decreases the maxima and increases the minima of the quantity being diffused. As it acts to
reduce the variance of the quantity being diffused, the solution does not oscillate in time. As the diffusion
flux is proportional to the gradient of the quantity being diffused (as modelled by Fourier’s law), it acts
quickly to smooth very sharp discontinuities. This means that any feature with a high spatial frequency,
which we will call ‘high mode number’ or ‘high wavenumber modes’ will be dissipated quickly.
For a full derivation of the heat equation see the Appendix, Section A.
A key tool to solve the linear PDEs in this course is the use of ‘Separation of Variables’. The key steps are:
1. Write out a clear definition of the problem. The equation, the initial conditions, the boundary condi-
tions and the parameters needed to define the specific problem (e.g. diffusion coefficient).
2. Assume that the dependence of the solution on each dimension can be separated, such that u(x, t) =
F(x)G(t). This means that the solution to the PDE can be split into multiple ODEs which can be
directly solved.
3. Solve the ODE for F(x), which in most cases will result in a Fourier Series solution of the form
Fn(x) = An cos(px)+Bn sin(px). At this stage you can most often determine the value of An or Bn
based on the boundary condition at x = 0, and p may then be determined from the second boundary
condition by ensuring that the wavelength of the mode gives the desired answer at x= L. See Section
3.3)
4. Solve the ODE for G(t), given the value for p which has been determined in the previous step. For
the heat equation this gives Gn(t) = B′ne−λ
2
n t where B′n is a constant and λn = cp.
5. Combine the two to form the full solution for a single mode (single value of n), e.g. un(x, t) =
(An cos(px)+Bn sin(px))e−λ
2
n t , where the constants have been combined, and finally note that any
value of n is a solution thus we may represent all solutions as:
u(x, t) = ∑
n=0,∞
un(x, t) = ∑
n=0,∞
(An cos(px)+Bn sin(px))e−λ
2
n t (18)
6. Finally, determine the remaining coefficients An or Bn by setting the time t = 0 and matching the
summation with the initial conditions, i.e.:
u(x,0) = ∑
n=0,∞
un(x, t) = ∑
n=0,∞
(An cos(px)+Bn sin(px)) (19)
which shows that the An and Bn are coefficients of the Fourier Series representing the initial condition
u(x,0).
The Fourier Series is an incredibly useful tool. A derivation of the Fourier Coefficients for the Sine series is
given in Appendix B. The key message is that any initial condition for a problem with either defined spatial
11
extent (more on this later) or which is periodic can be represented using a Fourier Series. The Fourier series
then gives an exact match to the initial condition, however only when and infinite number of terms is taken.
For initial conditions with particularly steep features, such as square waves, the number of modes (terms,
i.e. n=1, 2..) required to give a certain accuracy may be substantial (i.e. thousands). Thanks to the physical
behaviour of the wave equation, you may only need to retain a few terms to gain an accurate solution at late
time, when the high mode numbers have largely dissipated.
3.2 Expected Physics
Thanks to the Fourier Series representation, we can now discuss the behaviour of the solution in terms of
the fundamental modes.
If a mode is not present in the initial condition, it will never appear in the solution. For example, if only a
single mode n = 2 is present in the initial condition, you will never see modes n = 0,1,3,4,5, ... appear in
the solution. This is a result of the equation being linear, i.e. all derivatives contain only u, not u2, u3 etc.
This is a very useful property.
The time decay is proportional to e−λ 2n t , where λn ∝ n, i.e. each mode decays in amplitude at a rate ∝ e−n2t.
Thus high mode numbers will decay much more rapidly than low mode numbers. In fact, at very late times
you may be able to get a very good answer to a problem by using only a few of the lowest mode numbers
in the Fourier Series. For example, given a rod of length pi with c = 1 and Dirichlet boundary conditions,
e−λ 2n t = e−n2t . Thus mode 10 would reach 1/10th of it’s initial amplitude in a time of 0.23s, mode 100 would
reach 1/10th of it’s initial amplitude in 0.23ms. The time taken to decay to a fixed proportion of the initial
amplitude of the mode is proportional to 1/n2 (i.e. t = ln(Acurrent/Ainitial)/n2 where A is the amplitude, in
this example).
Finally, the spatial terms are fixed in time. If you wanted, these could be computed first and stored, thus
saving computational time.
The solution may be described in terms of eigenvalues and eigenfunctions, which are the individual modes
in the problem. The Eigenvalues are λn, and the Eigenfunctions are sometimes taken to be the spatial
shape which corresponds to the Eigenvalue, thus An cos(px)+Bn sin(px). It should be noted here that the
recommended text does not follow the strict definition of an Eigenfunction, it should include both the spatial
and temporal terms thus being un(x, t) = (An cos(px)+Bn sin(px))e−λ
2
n t . In this course we will accept either
interpretation. In many engineering analysis the ‘mode shape’ is also important, and this is the spatial shape
formed by a given eigenfunction, which may then decay in time (heat equation) or oscillate in the case of
vibrations.
3.3 Use of the Boundary Conditions
An important function of Fourier Series is the representation of Boundary Conditions (BCs) and Initial
Conditions (ICs). Here we take a slightly different approach to illustrating how to constrain the solution
using the boundary and initial conditions. From separation of variables you gain a solution for u(x,t) in the
form of:
u(x, t) =
∞
∑
n=0
(An cos(px)+Bn sin(px))e−λ
2
n t (20)
The first step is the application of the x = 0 boundary condition. The boundary conditions must apply at for
any position in time, thus we set t = 0, giving:
12
u(x,0) =
∞
∑
n=0
(An cos(px)+Bn sin(px)) (21)
Next we have a boundary condition at x = 0 and x = L. The boundary condition at x = 0 may be Dirichlet
or Neumann, and is used to set either An or Bn to be zero if both sides have the same boundary condition
type and they are homogeneous (0 value). The boundary condition at x = L then determines p.
Thus with two boundary conditions we have determined two parameters. We are then left with something
which is either:
u(x,0) =
∞
∑
n=0
An cos(px) (22)
or
u(x,0) =
∞
∑
n=1
An sin(px) (23)
where p is known. This can then be equated to the initial condition u(x,0), and a Fourier Sine or Cosine
series used to determine the coefficients which match the initial condition u(x,0).
3.4 Simple Inhomogeneous Problems
A note about inhomogeneous problems. I often have questions about converting the initial condition from
the full problem to the homogeneous problem. Taking the Week 3 tutorial about diffusion of Nitrogen
down a corridor, the initial conditions were Y (x,0) = 0.78 and the boundary conditions Y (0, t) = 0.78 and
Y (L, t) = 1.0. The initial condition for the homogeneous problem is determined from:
Y (x, t) = YH(x, t)+YS(x) (24)
thus
YH(x,0) = Y (x,0)−YS(x) (25)
Also, the boundary conditions are computed in a similar fashion, i.e.:
YH(0, t) = Y (0, t)−YS(0) (26)
Which should both be zero. So, what shape does this give us for the initial condition? The solution is the
sawtooth wave given by the following formulae:
IC:
YH(x,0) =−0.22x/L (27)
BC:
YH(0, t) = YH(L, t) = 0 (28)
Note that a confusion which arises is what value does the initial condition take at x = L, given that the
boundary condition YH(L, t) = 0 seems to contradict the IC which gives Yh(x,0) = −0.22. In all cases the
13
boundary condition has ‘priority’ over the initial condition, i.e. the value of YH(L, t) = 0 always in this
example. That is an essential requirement of it being a homogeneous problem.
3.5 More Complex Boundary Conditions
In reality boundary conditions may be a mixture of homogeneous and inhomogeneous, or a mixture of
Dirichlet and Neumann. Let’s start with an example of a heat conduction through a plate, where one side
has a specified heat flux, the other side a fixed temperature.
The model problem. The governing equation is Tt−c2Txx = 0, and the boundary conditions are c2Tx(0, t) =
H, where H is the specified heat flux, and T (L, t) = T1 where T1 is the fixed temperature at x = L. Assume
that the initial temperature of the plate is T (x, t) = T1.
As a side note, T (x, t) can be chosen to be different to T1 and the methodology presented here will still
hold. This problem occurs in many places - for example in designing heat exchangers a Computational
Fluid Dynamics study is used to provide the heat flux, then the necessary cooling load determined.
There are three new complexities in this problem:
1. Inhomogeneous boundary conditions, thus requiring a steady state solution for which a functional
form must be assumed
2. Mixed boundary conditions lead to an eigenvalue which uses non-integer modes - i.e. (n−1/2)pi/L
instead of the more familiar npi/L.
3. The Fourier Series representation of the initial condition must use these non-integer modes
Let’s tackle these one at a time for the given problem. Firstly, we split the solution into a homogeneous and
steady state, such that T (x, t) = T ss(x)+T H(x, t)
3.5.1 Step 1: Steady State solution T ss(x)
Given the above boundary conditions, the first step is to transform from the heat flow to an equivalent
temperature gradient. Given c2, the temperature gradient at the boundary must be Tx(0, t) = T 0x = H/c
2.
Thus this is the actual boundary condition we will impose on the problem.
We now assume that the steady state solution is given by T ss(x) = Ax+B, with A and B being coefficients
to determine. This satisfies the steady state heat conduction equation Txx = 0. We now calibrate A and B to
satisfy the given boundary conditions of T (L, t) = T1 and Tx(0, t) = T 0x .
Given T ssx = A, then A = T
0
x . Then, T
ss(L, t) = T1 = T 0x L+B, thus B = T1−T 0x L. Thus the steady state
solution is
T ss(x) = T 0x (x−L)+T1 (29)
This raises a question - what steady state solutions do we have for other boundary condition choices? Well,
the options are:
• Both Dirichlet: T (0, t) = φ0, T (L, t) = φ1, T ss = φ0+(φ1−φ0)x/L
• Mixed BC Neumann/Dirichlet: Tx(0, t) = φ0, T (L, t) = φ1, T ss = φ0(x−L)+φ1
14
• Mixed BC Dirichlet/Neumann: T (0, t) = φ0, Tx(L, t) = φ1, T ss = φ0+φ1x
• Both Neumann: Tx(0, t) = φ0, Tx(L, t) = φ1, T ss = φ0x+(φ1−φ0) x22L . note that this last one is a little
different - it requires an additional ‘sink’ term in the equation as otherwise the temperatures
in the system will just increase without bound.
Note that this can be expanded to time dependent boundary conditions but it is beyond the scope of this
course.
3.5.2 Step 2: Formulation and Solution of the Homogeneous Problem
Given T H(x, t) = T (x, t)−T ss(x) we can now calculate the initial and boundary conditions for the homo-
geneous problem. As expected, the boundary conditions can be shown to be T Hx (0, t) = 0 and T
H(L, t) = 0.
The initial conditions become T H(x,0) = T (x,0)−(T 0x (x−L)+T1) = T 0x (L−x)+T1−T1, thus T H(x,0) =
T 0x (L− x).
Our problem is now completely defined, we now just need to solve T Ht − c2T Hxx = 0 using separation of
variables and applying the boundary conditions. Here we’ll shortcut the derivation and go straight to the
spatial Eigenfunction:
T Hn (x) = An cos(px)+Bn sin(px) (30)
which has the first spatial derivative:
T Hn,x(x) =−pAn sin(px)+ pBn cos(px) (31)
Noting that T Hx (0, t) = 0 tells us that Bn = 0. We are left with:
T Hn (x) = An cos(px) (32)
To satisfy the second boundary condition, T H(L, t) = 0, we need:
T Hn (L) = An cos(pL) = 0 (33)
By sketching out a cosine you can see that this is satisfied only if p = (n−1/2)pi/L, where n is an integer
n≥ 1.
Solving the separated ODE for the temporal behaviour gives the same solution as the Dirichlet-Dirichlet
problem considered in the first lecture, thus each mode will be described by:
T Hn (x, t) = An cos(px)e
−λ 2n t , (34)
with λn = cp where An has absorbed a second possible constant from the solution of the temporal ODE.
The full solution for T H is the summation over an infinite number of modes:
T H(x, t) = ∑
n=1,∞
An cos(px)e−λ
2
n t , (35)
Now we need to satisfy the initial conditions by determining An. Setting t = 0 gives:
15
T H(x,0) = ∑
n=1,∞
An cos(px). (36)
3.5.3 Step 3: Fourier Series for non-integer modes
This looks like a straightforward application of the Fourier Series. Although a sharp eyed observed may
question if it is still valid as p is no longer an integer multiple of pi/L as it has been in all previous cases.
It turns out that this isn’t a problem - the key is that the definition of the Fourier Series relies on each mode
being separated by a multiple of pi/L, and this is still the case here. For example, for n = 1, p = pi/2L
and for n = 2, p = 3pi/2L and so on. This means that the modes are still orthogonal and so satisfy the
definition of the Fourier Series. The derivation of the Fourier Series coefficients then follows as before, and
the coefficients are given by the same formulae, for example:
An =
2
L
∫ L
0
f (x)cos(px)dx, (37)
For the cosine modes we have here. Let’s derive the initial condition now. We need to have T H(x,0) =
∑n=1,∞An cos(px) = T 0x (L− x), thus:
An =
2
L
∫ L
0
T 0x (L− x)cos(px)dx =
2T 0x
L
∫ L
0
(Lcos(px)− xcos(px))dx (38)
The first term in the integral can be dealt with easily, the second term by application of integration by parts∫
uv′dx = u′v− ∫ u′vdx where u = x and v′ = cos(px):
An =
2T 0x
L
[
L/psin(px)− x/psin(px)− 1
p2
cos(px)
]L
0
(39)
The first two terms in square brackets are zero when x = 0. When x = L, they cancel. The final term in
square brackets is −1/p2 when x = 0 and due to the definition of p, cos(pL) = 0. The An thus simplify to:
An =
2T 0x
Lp2
(40)
Given the final solution to our mixed boundary homogeneous problem as:
T H(x, t) = ∑
n=1,∞
2T 0x
Lp2
cos(px)e−λ
2
n t . (41)
and the final solution for T (x, t) = T H(x, t)+T ss(x):
T (x, t) = T 0x (x−L)+T1+ ∑
n=1,∞
2T 0x
Lp2
cos(px)e−λ
2
n t . (42)
Plot this up in Matlab to check that it makes sense! A check of whether you have got the derivation right is
in implementing it and sanity checking against the initial conditions and boundary conditions. Also check
that if the length of the domain is change from 1, to a value which isn’t 1, is the solution still correct? That
helps find where you may have accidently dropped a factor of L in the derivation.
16
3.5.4 Periodic Boundary Conditions
Some cases may end up with periodic boundary conditions - an obvious one being conduction around a
thin ring - this can be unrolled to be a one dimensional problem where x = rθ . In this case there is a new
boundary condition, where both T (−L, t) = T (L, t) and Tx(−L, t) = Tx(L, t). Note that L = rpi , and the
choice of the location of x = 0 and the periodic boundary is deliberate to simplify the maths below. How
do we deal with this?
Following separation of variables, we are left with the spatial solution:
Fn(x) = An cos(px)+Bn sin(px) (43)
Applying the periodic boundary condition on T we note that:
An cos(−pL)+Bn sin(−pL) = An cos(pL)+Bn sin(pL) (44)
using the knowledge of odd/even behaviour of sin and cos about x = 0:
An cos(pL)−Bn sin(pL) = An cos(pL)+Bn sin(pL) (45)
thus 2Bn sin(pL) = 0. This could be satisfied with either Bn = 0 or p = npi/L. Applying the periodic
boundary condition on Tx we gain:
−pAn sin(−pL)+ pBn cos(−pL) =−pAn sin(pL)+ pBn cos(pL) (46)
i.e.
pAn sin(pL)+ pBn cos(pL) =−pAn sin(pL)+ pBn cos(pL) (47)
thus 2pAn sin(pL) = 0.
Thus there are two possible solutions, a trivial one which is not very useful (An = 0, Bn = 0) and the useful
solution which is p = npi/L where n = 0,1, ..∞.
The solution in the periodic case is then:
T (x, t) = ∑
n=0,∞
(An cos(px)+Bn sin(px))e−λ
2
n t . (48)
where An and Bn are given by the Fourier coefficients representing the initial condition.
17
4 Key Concepts in Numerical Methods
4.1 Introduction
This section will outline the application of our simple method of constructing approximation of derivatives
to solve governing equations. It will also introduce important concepts in more detail, including Order of
Accuracy, Stability and Verification and Validation.
4.2 Order Of Accuracy
The order of accuracy of a numerical method tells you how quickly the numerical solution approaches the
exact solution as the sample spacing ∆x decreases. Here we will use interchangeably the terminology of
decreasing the mesh spacing, or increasing the number of points as meaning the same thing - effectively
assuming that we are solving a problem for a fixed domain size. This is illustrated for the Central Difference
scheme in Section 1.3, where the leading order truncation error can be written as:
∂ f (i)
∂x
=
f (i+1)− f (i−1)
2∆x
− ∆x
2
6
∂ 3 f
∂x3
+ · · · · · · . (49)
This scheme is second order accurate, i.e. the error is proportional to ∆x2. Note that the full error also
includes the third derivative of f - hence the order of accuracy doesn’t actually tell you the magnitude of
the error for all solutions, only its dependence on the spacing between the points. You must rely on careful
Verification and Validation to ensure that actual error of the numerical scheme is within that required for
your given application. Note that the order of the derivative does not change the order of accuracy of the
scheme, as it is assumed that as the number of grid points increases, the derivatives become constant. Back
to the error, for a second order scheme, as you half the grid spacing you will decrease the error by a factor
of approximately 4. This is approximate, as the error term in Eqn. (49) is what is termed the leading order
error - there are other higher order error terms which also add to the error. It just happens that as the spacing
between the points becomes smaller, the other terms become negligible compared to the leading order.
In the lectures we introduced a second approximation - namely:
fx =
f (i+1)− f (i)
∆x
+
∆x
2
∂ 2 f
∂x2
+ · · · · · · . (50)
We noted that this is first order accurate, i.e. the error is proportional to ∆x1, where the 1 is included for
emphasis. Here if the grid spacing is halved, the error decreases by a factor of two, not four as for the
central difference scheme. This means that for roughly the same computational effort, as the grid is refined,
the second order scheme becomes substantially more accurate than the first order scheme - hence if there
were no other considerations we would always choose the second order scheme over the first.
4.3 Stability
A numerical scheme must be stable, otherwise it is no use at all. Unstable schemes do not just give a bad
answer, it is often in violation of physics and results in the simulation crashing. It is possible to have a very
high order method which looks beautifully constructed, being absolutely useless in practise as it is unstable.
An important concept in stability is that if a scheme is unstable, errors will commonly grow exponentially
thus rapidly rendering a solution useless. We cannot yet look at this, as to determine stability you must
analyse the numerical method in full as applied to the chosen governing equation - it is not typically applied
18
to e.g. a finite difference representation of the derivatives. We will cover this in more detail for the heat in
Section 5.1.2.
4.4 Consistency
A third requirement is consistency. This means that as the grid spacing tends to zero, the approximation
will become exact. This is clearly true in the above two approximations, where setting ∆x = 0 will make
the error term 0.
4.5 Verification and Validation
Verification is ensuring that your numerical method is derived and implemented correctly. It involves
checking that the numerical method converges to the actual solution of the governing equations, and con-
verges at the expected order of accuracy.
Validation checks if the governing equation is the correct choice to model the actual physics of the problem.
This usually means validation against experimental data, which of course has the ’full physics’ of the
problem. It checks that the modelling assumption in the derivation of the governing equations are sound.
In this course we will assume that the governing equations are sound, and will use the exact solution to
undertake both verification and validation in one step.
4.6 Quantifying the Errors with Error Norms
To validate, we need to quantify the errors. The typical choice is to use error norms, denoted Lp where p is
an integer ranging from 1,2...∞, where the most common choices are the L1, L2 and L∞ norms respectively.
They are calculated from:
Lp =
(
∆x
nx
∑
i=1
| fi− f (xi)|p
)1/p
(51)
where fi is the approximation at point i, f (xi) the exact solution, i is an individual point ranging from the
first point to point nx, the last. Careful inspection will tell you that the L1 error is the mean absolute error,
L2 the rms, and L∞ the maximum error. You should calculate this for your solution on multiple grid sizes.
Note that for problems which have two dimensions in space, such as in solving the Laplace or Poisson equa-
tions in (x,y), the additional dimension means that the error norm definition must be changed as follows:
Lp =
(
∆x∆y
nx
∑
i=1
ny
∑
j=1
| fi, j− f (xi,y j)|p
)1/p
(52)
4.7 Computing the Observed Order of Accuracy
The second step is to compute the observed order of accuracy. This is the actual order of accuracy of your
implementation of the algorithm, and so tests for bugs in the algorithm derivation (e.g. errors in the Taylor
Series approximations used to derive the approximations to the derivatives) or just straightforward coding
errors.
19
Table 1: An example of an actual presentation of errors and convergence rates (B Thornber, M Groom, D
Youngs, A five-equation model for the simulation of miscible and viscous compressible fluids, Journal of
Computational Physics 372, 256-280, 2018)
N L1 L2 L∞ O(L1) O(L2) O(L∞)
64 5.2677e-02 7.2627e-02 1.7744e-01 - - -
128 1.4011e-02 1.8492e-02 3.7163e-02 1.9106 1.9736 2.2554
256 3.4566e-03 4.5593e-03 9.2933e-03 2.0191 2.0200 1.9996
512 8.7023e-04 1.1486e-03 2.3790e-03 1.9899 1.9890 1.9659
1024 2.1800e-04 2.8778e-04 5.9729e-04 1.9971 1.9968 1.9938
2048 5.4519e-05 7.1936e-05 1.4825e-04 1.9995 2.0002 2.0104
4096 1.3724e-05 1.8020e-05 3.5851e-05 1.9900 1.9971 2.0480
We compute the order of accuracy based on the error norms. We assume that the behaviour of the error
norm is dominated by the leading order term. Thus L = C∆xn, with n the order of accuracy. Run two
calculations, with two substantially different grid sizes (usually varying by a factor of 2 from each other).
If the first mesh spacing ∆xA gives error LA, and the second mesh spacing ∆xB gives error LB, then:
LA = C∆xnA (53)
LB = C∆xnB (54)
Assuming that the problem is well resolved then C should be approximately a constant (it would be pro-
portional to derivatives of the function being solved for). There are two unknowns and two constants, so
solving gives:
n =
log
(
LA
LB
)
log
(
∆xA
∆xB
) (55)
Note that this is a model of the error - the real error consists of many terms with many different derivatives.
So this will only give the exact convergence rate as the number of points becomes very high. Table 1 shows
an actual convergence study published by Thornber, Groom and Youngs in the Journal of Computational
Physics in 2018 for a system of equations for compressible turbulent mixing. Important points to note are:
• The grid resolutions (number of points) vary by a factor of 2 from one calculation to the next - it
is very important to make a large change in grid size when estimating the error. A factor of 2 is
recommended.
• Three error norms are computed
• The convergence rate based on all three are calculated. Note how the error is not exactly 2 at any
point, in fact at very low grid resolutions as shown in the lecture slides the observed value of n can
differ dramatically from the expected. This is because the solution is not well resolved, meaning that
(i) the gradient term multiplying the factor of ∆xn changes as the number of points is changed and (ii)
the non-leading error terms contribute significantly to the error.
• At high grid resolutions the error asymptotes towards second order, which is expected for the numer-
ical scheme used
All convergence studies undertake in this course should be presented in exactly this manner.
20
5 Numerical Solution of the Heat Equation
5.1 Forward in Time, Central in Space (FTCS)
Rationale for Discretisation
In the computation of a given problem, we usually have initial and/or boundary conditions. Thus, at time
t = 0, we typically know the value of the function to be evolved, say u(x,0). So, our initial condition may
consist of known values of u(xi,0), where xi indicates the spatial location of i discrete points in a domain.
For example, if the problem is on the domain 0 ≤ x ≤ 2 and the total number of point is 5, then x1 = 0,
x2 = 0.5, x3 = 1.0, x4 = 1.5 and x5 = 2.0 and ∆x = 0.5. Assuming that we know the values of u at every
single point at time zero, we can use these known values to evolve the solution forward in time.
As an example, we know that ut = c2uxx. Thus, if you can estimate uxx and you know the value of c2, you
can work out ut directly from the heat equation itself. With ut known, you can predict the value in the future
using, for example, a first order in time Taylor Series expansion u(t+∆t) = u(t)+∆tut .
5.1.1 Discretisation
We need a numerical approximation for uxx(xi, tn) and then to evolve in time we need one for ut(xi, tn).
Here i is the integer labelling points in space, n is the integer labelling points in time, for example, if the
calculation starts at t = 0, then t = n∆t.
To gain the approximation for uxx, we can actually apply our central difference method twice, using
uxx(xi, tn) =
ux(xi+1/2, tn)−ux(xi−1/2, tn)
∆x
. (56)
Note that here tn is the time at time step n, not t raised to the power of n. Here the first derivative ux(xi+1/2, t)
is calculated at the midpoint of xi and xi+1. These first derivatives can be further computed from central
differences themselves:
ux(xi+1/2, t
n) =
u(xi+1, tn)−ux(xi, tn)
∆x
. (57)
ux(xi−1/2, tn) =
u(xi, tn)−ux(xi−1, tn)
∆x
. (58)
Combining Eqn. (58) and (56) gives:
uxx(xi, tn) =
u(xi+1, tn)−2u(xi, tn)+u(xi−1, tn)
∆x2
=
uni+1−2uni +uni−1
∆x2
. (59)
You may confirm yourself that this is a second order accurate approximation.
The time derivative is approximated using a first order, forward in time method called Euler time integration.
This is:
ut(xi, tn) =
u(xi, tn+1)−u(xi, tn)
∆t
=
un+1i −uni
∆t
. (60)
21
Although we advise in lectures that most scientific computing requires at least second order accuracy for
all terms, the heat equation is an exception to this. The reason for this is that although the error in the above
approximation for ut is first order in ∆t, we will see in the next section that to ensure stability of the resultant
scheme we will have to constrain ∆t ∝ ∆x2. Thus the first order in time error will reduce at a second order
rate (∝ ∆x2) because of the dependence of the time step size on the grid spacing ∆x.
Bringing the two numerical approximations of the derivatives into the heat equation ut = c2uxx becomes:
un+1i −uni
∆t
= c2
uni+1−2uni +uni−1
∆x2
(61)
Rearranging:
un+1i = u
n
i +
c2∆t
∆x2
(
uni+1−2uni +uni−1
)
(62)
Equation (62) is the Forward in Time, Central in space (FTCS) and is the most straightforward discretisation
of the heat equation. It is consistent, and has a formal error proportional to ∆t in time and ∆x2 in space.
It takes the known values of uni at the current time level n on the right hand side of the equation, and uses
them to predict the new values un+1i at the next time instant n+1 which is at time t
n+∆t.
Explicit Scheme: The use of the Forward in Time scheme makes this an Explicit scheme. Explicit schemes
directly calculate the value of un+1 using only values un at the current time step which are all known (the
RHS of Eqn. (62). We will explain the second major class of scheme shortly, which is Implicit.
This finishes the presentation of the scheme. However, we have one outstanding item (i) is the scheme
stable?
5.1.2 von Neumann Stability Analysis
von Neumann was a pioneer of mathematics, physics and computing. He developed a method of stability
analysis which may be applied easily to the discretised form of the governing equation/s. This requires the
following steps:
1. Assume a solution to the governing equation which can represent a range of the possible physical
behaviour
2. Substitute this into the discretised governing equations, and simplify
3. Calculate the amplification factor G describing the time dependent behaviour of each mode. If the
amplitude growth rate |G| > 1, the scheme is unstable, if |G| ≤ 1 it is stable. Rearrange to give an
expression for the time step size where G = 1, the limit of stability. Note that if G is composed
of a real and imaginary part then |G|2 = Re[G]2 + Im[G]2, e.g. if G = cosβm + I sinβm then |G|2 =
cos2βm+ sin2βm = 1..
This will tell you (i) is the scheme stable and (ii) under what conditions
Step 1: Assumed Solution
The basic approach is to assume a general form of the solution to the governing equation. Most solutions
to partial differential equations can be described as a Fourier Series in space and oscillatory, exponentially
decaying, or exponentially growing in time. In stability analysis, the complex exponential form of the
Fourier Series is employed:
22
f (x) =
∞
∑
k=−∞
ckeIknx (63)
where it should be noted that from Euler’s Formula:
eIknx = cos(knx)+ I sin(knx), e−Iknx = cos(knx)− I sin(knx) (64)
where ck is a complex wave amplitude and k = pinL is the wave number. By letting n be the temporal
index with increment ∆t and j be the spatial index with increment ∆x, the discrete version of the complex
exponential Fourier Series can be written as:
u(x) =
∞
∑
k=−∞
ckeIkni∆x (65)
Using this, the solution to the PDE is now written as:
u(x, t) =
∞
∑
k=−∞
GneIkni∆x (66)
where Gn(k, t) is the time dependent amplitude of the mode. This will capture the potential oscillatory,
exponentially decaying, or exponentially growing behaviour of the PDE solution in time. If |G| > 1 then
the amplitude will grow exponentially in time, if |G| = 1 the amplitude will not grow, and if |G| ≤ 1 the
amplitude will decrease. Thus a stable scheme has |G| ≤ 1, and the boundary of stability is |G|= 1.
For linear PDEs in this course, all modes demonstrate the same physical behaviour. This means that by
understanding the behaviour of one single mode, you can understand the behaviour of all the modes in the
solution. Thus in von Neumann analysis we assume that the solution consists of just one single mode, thus:
uni = G
neIkni∆x (67)
Step 2: Substitution into the Discretised Equation
Starting from Eqn. (62), substitute in Eqn. (67), giving:
Gn+1eIkni∆x = GneIkni∆x+
c2∆t
∆x2
(
GneIkni+1∆x−2GneIkni∆x+GneIkni−1∆x
)
(68)
Substantial simplification can be made by dividing through both sides by GneIkni∆x.
G = 1+
c2∆t
∆x2
(
eIkn∆x−2+ e−Ikn∆x
)
. (69)
Set σ = c
2∆t
∆x2 :
G = (1−2σ)+σ(e−Ik∆x+ eIk∆x) = 1−2σ(1− cos(k∆x)), (70)
where we have used Euler’s formula (Eqn. (64)) to express the sum of the two exponentials as a cosine
term.
23
Step 3: Determination of any Necessary Stability Criteria
To ensure that we have a stable scheme, we need to set |G| ≤ 1. This gives:
|1−2σ(1− cos(k∆x))| ≤ 1 (71)
Note that 0 ≤ (1− cosθ)≤ 2. Taking the lower limit where (1− cosθ) = 0, then |G|= 1 and the scheme
is stable. At the upper limit:
|1−4σ | ≤ 1, (72)
and it can be seen that in order to prevent the item in the brackets becoming less than −1, 0 ≤ σ ≤ 1/2.
Thus:
σ ≤ 1/2,
c2∆t
∆x2
≤ 1/2,
∆t ≤ ∆x
2
2c2
, (73)
Thus the FTCS scheme for the heat equation is stable only if the time step size is chosen to satisfy Eqn.
(73). Above this criteria the scheme will be unstable, as you will see in the tutorial in Week 4. An important
impact of this limitation is that if you have to run a calculation for a fixed time interval t, then the number
of time steps will be t/∆t ∝ 1/∆x2. So each time the grid spacing is halved, the required number of time
steps and hence the computational effort will increase by a factor of 4!
Thus, for short time calculations, the FTCS scheme is quite efficient. However, if you need to run a
calculation for a very long time, you may need to find a numerical method which does not have such a
stringent time step size restriction. This gives a strong motivation for the development of schemes which
are unconditionally stable, i.e. they are stable regardless of the choice of ∆t, thus ∆t can be chosen solely
based on accuracy requirements.
5.2 Backward in Time, Central in Space Scheme (BTCS)
Instead of using Eqn (59) to discretise uxx, how about the following:
uxx(xi, t) =
un+1i+1 −2un+1i +un+1i−1
∆x2
. (74)
where the derivative estimate is now calculated using values of u at the next time level, n+1. Thus:
un+1i = u
n
i +
c2∆t
∆x2
(
un+1i+1 −2un+1i +un+1i−1
)
(75)
This is the Backward in Time, Central in Space (BTCS) Scheme. If you repeat the von Neumann stability
analysis you will find that it is unconditionally stable, i.e. you can choose any value of ∆t and it will run
stably. However, how do we evolve this forward in time, since the value of u at every point at the end of the
time step depends on the value of u at every other point at the next time step, i.e.:
24
un+1i (1+2σ)−σ
(
un+1i+1 +u
n+1
i−1
)
= uni (76)
We have the same number of equations as unknowns, yet the only points we know the solution at time n+1
are the boundary conditions and the single term on the right hand side. For a full one dimensional domain,
we can write this in matrix form as:
Aun+1 = un, (77)
If the boundary conditions are Dirichlet, so un+1 = un at the boundary, then A is given by:
A =

1 0 0 0 0 0
σ (1−2σ) σ 0 0 0
0 σ (1−2σ) σ 0 0
0 0
. . . . . . . . . 0
0 0 0 σ (1−2σ) σ
0 0 0 0 0 1

(78)
The solution is then given by solving the linear system of equations in the form Ax= b, where the solution is
given by un+1 = A−1un. There are many different approaches to solve such a system of equations, where the
most efficient depends on the exact structure of the matrix A. One approach called the ‘Thomas Algorithm’
is highlighted in the lecture slides - this is useful for tridiagonal matrices.
This takes more computational expense than the Explicit FTCS scheme for one time step, but as the time
step may be much larger as it is unconditionally stable, it may be more computationally efficient when
integrating the equations over a long time period. Thus the choice of scheme depends on the problem that
you would like to solve.
5.3 Conclusions
This section has introduced the fundamentals of numerical methods, in the context of the heat equation.
Key items covered include the concepts of validation and verification, stability, accuracy, discretisation of a
PDE using explicit and implicit schemes. Explicit schemes are simple and fast, however the time step size
is limited by the stability criteria. Implicit schemes require the inversion of a potentially large matrix to
solve the problem, but have the advantage that the time step size may be arbitrarily large, so can be chosen
to give the desired accuracy.
25
6 The Wave Equation
6.1 Derivation of the Wave Equation
There is an excellent video covering the derivation of the wave equation on Canvas so I will not repeat that
here. Rather I want to emphasize the key underlying assumptions, and those are:
1. the mass/unit length is constant
2. the tension is much greater than the force due to gravity, i.e. T >> mg
3. the string performs small transverse motions so that e.g. sinθ ≈ θ etc.
These assumptions are the limitations of the application of the governing equation. These assumptions limit
the accuracy when you attempt to apply the equation to a real life case which does not satisfy the specified
limitations. In that case, it is necessary to validate your approach against experimental observations to
ensure that the error is tolerable for your specific application. Applying these assumptions, and examining
the force balance in the horizontal and vertical direction on a small segment of a long cable aligned in the
horizontal direction gives you the wave equation:
utt =
T
ρ
uxx (79)
in a more general form:
utt = c2uxx (80)
Examining the form of this equation shows that c1 = ρ/T > 0, c3 = −1 and so 4c1c3 < 0 and so the
governing equations is hyperbolic. A hyperbolic equation consists of waves travelling in space or time.
6.2 Analytical Solution
Solution is through application of the method of ‘Separation of Variables’, following an identical process
to the Heat Equation. The key steps are:
1. Write out a clear definition of the problem. The equation, the initial conditions, the boundary condi-
tions and the parameters needed to define the specific problem (e.g. tension, mass per unit length).
2. Assume that the dependence of the solution on each dimension can be separated, such that u(x, t) =
F(x)G(t).
3. Solve the ODE for F(x), which in most cases will result in a Fourier Series solution of the form
Fn(x) = An cos(px)+Bn sin(px). At this stage you can most often determine the value of An or Bn
based on the boundary condition at x = 0, and p may then be determined from the second boundary
condition by ensuring that the wavelength of the mode gives the desired answer at x = L.
4. Solve the ODE for G(t), given the value for p which has been determined in the previous step. Here
the solution differs from the heat equation, as the wave equation has a second derivative with respect
to time, thus we solve Gtt + c2 p2G = 0, or alternatively Gtt +λ 2n G = 0 where λn = cp. The solution
to this consists of cosine and sine modes. Thus for the wave equation:
26
Gn(t) = A∗n cos(λnt)+B
∗
n sin(λnt) (81)
This represents oscillations in time - an enormous difference between the physical solution of the
wave equation and the heat equation (which decays in time).
5. Combine the two to form the full solution for a single mode (single value of n),
un(x, t) = (An cos(px)+Bn sin(px))(A∗n cos(λnt)+B
∗
n sin(λnt)) (82)
At this point, alert readers will realise that we have an additional coefficient compared to the heat
equation as we have two terms with time dependence now, thus we need more information to solve
the equation. Looking only at the temporal variation, setting t = 0, we can understand the behaviour
of the temporal terms:
Gn(0) = A∗n,
dGn(t)
dt
= λnB∗n (83)
Thus the amplitude of the cosλnt term represents the amplitude of the initial condition, and the
sinλnt term represents the amplitude of the time rate of change of the initial condition - i.e. the initial
velocity. Hence, it turns out that the required information that we need is the initial velocity.
As with the heat equation, any value of n is a solution thus we may represent all solutions as:
u(x, t) = ∑
n=0,∞
un(x, t) = ∑
n=0,∞
(An cos(px)+Bn sin(px))(A∗n cos(λnt)+B
∗
n sin(λnt)) (84)
One way of writing this which may help with the physical interpretation is as follows:
u(x, t) = ∑
n=0,∞
un(x, t) = ∑
n=0,∞
(Ac1n cos(px)+B
c1
n sin(px))cos(λnt)︸ ︷︷ ︸
Fourier Series representing the initial amplitude distribution in space
+ (85)
∑
n=0,∞
(Ac2n cos(px)+B
c2
n sin(px))sin(λnt)︸ ︷︷ ︸
component determined from the initial velocity distribution in space
(86)
where I have combined the constant coefficients into new constant coefficients Ac1n ,B
c1
n ,A
c2
n ,B
c2
n .
In summary, for the wave equation we need to know both the spatial distribution of the initial
amplitude and the initial velocity.
6. Finally, determine the coefficients by:
(a) using the boundary condition at x = 0 to determine which of the cos(px) and sin(px) terms will
remain in the solution, and then the boundary condition at x = L to determine p.
(b) following on from the above step, set the time t = 0 and match the two summations with the
initial amplitude and velocity, i.e.:
u(x,0) = ∑
n=0,∞
(Ac1n cos(px)+B
c1
n sin(px)) (87)
ux(x,0) = ∑
n=0,∞
(Ac3n cos(px)+B
c3
n sin(px)) (88)
where Ac3n and B
c3
n are constant coefficients equal to λnAc2n and λnBc2n . The fact that these coef-
ficients include λn does not make the calculation any more complex - they are a still a constant
for each mode and a Fourier Coefficient.
27
6.3 Physical Behaviour
The solution to the wave equation consists of a Fourier Series in space, and oscillatory behaviour in time.
Each un is an eigenfunction. If we consider the guitar string example in the lecture slides, the eigenfunction
is:
un(x, t) = Bc1n sin(px)cos(λnt) (89)
The eigenfunction corresponds to the eigenvalue λn = cnpi/L. The eigenfunction may be interpreted as
the shape of a single mode in space and time. Inspecting Eqn. (89) it can be seen that the eigenfunction has
an amplitude which is determined by the amplitude of the mode in the initial condition and a wavelength
which ranges from mode n = 1,∞. As with the heat equation, modes not present in the initial condition
will never appear in the solution. Thus the spatial variation is that of a sine wave.
The temporal variation is represented by the cosine term. It starts at 1, and then oscillates in time, returning
to an amplitude of 1 in a time t = 2pi/λn. Note that the maximum and minimum amplitudes of the modes
in time are bounded, i.e. they oscillate between values fixed by the initial conditions. The period of the
oscillation is 2pi/λn, thus the frequency is:
F(Hz) =
λn
2pi
=
cn
2L
(90)
Thus the frequency can be increased by increasing c, decreasing L or examining a higher mode number
n. Given c =
√
T
ρ then c can be increased by either increasing the tension or decreasing the mass per
unit length. This can help guide engineering decisions - there are many systems where addition metal has
been added to a design with the pure goal of shifting the fundamental frequency away from a range where
resonance may occur. Resonance occurs when a time dependent force is applied to a system at close to an
eigenvalue of the system, and causes disastrous amplification of the mode amplitude.
Mode n = 1 is the fundamental mode. Others are denoted overtones or harmonics.
Eigenvalues and Eigenfunctions are properties of a specific system, dependent on the initial conditions,
boundary conditions and material properties.
6.4 Conclusions
In this section we covered:
• Physical phenomena represented by the wave equation
• Analytical solution of the wave equation using Fourier Series for an arbitrary initial condition
This has obvious applications in estimating solutions to more complex problems, sanity checks and verifi-
cation and validation of numerical schemes.
28
7 Numerical Solution of the Wave Equation
7.1 Second Order Accurate Central Difference in Time and Space
Here is where we start reaping the benefits of the hard work put in discretising the heat equation. For the
wave equation we have two second derivatives utt and uxx to approximate. Fortunately, we already know a
second order accurate scheme to approximate a second derivative with. This is:
uxx(xi, tn) =
uni+1−2uni +uni−1
∆x2
− ∆x
2
12
uxxxx, utt(xi, tn) =
un+1i −2uni +un−1i
∆t2
− ∆t
2
12
utttt , (91)
Given these two approximations, we have then have a discretised form of the governing equation utt = c2uxx
which when simplified results in:
un+1i = 2u
n
i −un−1i +
∆t2c2
∆x2
(
uni+1−2uni +uni−1
)
(92)
This equation combines two second order accurate approximations to make a final approximation which
is at least second order accurate. Here I note at least, because surprisingly if you can supply this method
with very accurate initial conditions (perhaps exact) then the above discretisation can be shown to be fourth
order accurate on grids with constant spacing (more on this later). The discretisation is consistent since it
reduces to the wave equation as ∆x→ 0 and ∆t→ 0.
An important difference here is the inclusion of un−1i , arising from the central discretisation of the time
derivative. This means that in the first time step we now need two values of ui, one at t = 0 and the other at
t = −∆t, i.e. we need to specify both u0i and u−1i . There are two ways to tackle this, one is basic but poor
accuracy, the other is much better but needs more information.
Basic but poor accuracy: Assume u−1i = u
0
i . It is important to relate this assumption to our knowledge of
interpolation methods to see that this is a first-order accurate approximation. One thought may be that
as this is just for one time step, it doesn’t matter. It actually turns out that when you run calculations the
error made in the first time step is actually bigger than the errors made in all other time steps put together!
This may be a huge surprise, given you may have tens of thousands of time steps, but it is true (test it and
see!).
Recommended 2nd order accurate approach: Let’s consider what this extra requirement represents. We
have seen that for the analytical solution we need two initial conditions, one is the amplitude, the other is
the initial velocity. It turns out that our requirement to specify u−1i is actually needed to represent the rate
of change of velocity at time zero (acceleration -which is what uxx means physically). Thus by getting this
wrong, we are setting up the wrong initial velocity. So, the better approach is to estimate our initial velocity
from known quantities, i.e. noting that the velocity at time t = 0 may be represented by:
ut(xi,0) = g(xi) =
u1−u−1
2∆t
, (93)
where g(xi) is the initial velocity at each point. Rearranging gives:
u−1 = u1−2∆tg(xi). (94)
We now insert this into the discretised form of the wave equation, Eqn. (92), to gain a new discretisation
which is only used in the first time step:
29
u1i = u
0
i +∆tg(xi)+
∆t2c2
2∆x2
(
u0i+1−2u0i +u0i−1
)
(95)
Now this equation can be directly solved, given g(xi). This equation should only be used for the first time
step, and then Eqn (92) used for all other steps.
7.2 von Neumann Stability Analysis
The remaining question is, is the scheme stable and what limits if any may there be on the choice of ∆t?
Here we follow the methodology laid out in Section 5.1.2 for the heat equation.
Step 1: Assumed Solution: As with the heat equation, we assume that the solution is a Fourier Series in
space, multiplied by Gn, where G is a complex coefficient, i.e. u(x, t) = ∑∞k=−∞GneIkni∆x. We look at only
one single mode, noting that all modes will behave in the same way, thus the assumed solution is:
uni = G
neIkni∆x (96)
Step 2: Substitution into the Discretised Equation. Substitute Eqn (96) in to Eqn (92), giving:
Gn+1eIkni∆x = 2GneIkni∆x−Gn−1eIkni∆x+ c
2∆t2
∆x2
(
GneIkni+1∆x−2GneIkni∆x+GneIkni−1∆x
)
(97)
Substantial simplification can be made by dividing through both sides by Gn−1eIkni∆x:
G2 = 2G−1+Gα2
(
eIkn∆x−2+ e−Ikn∆x
)
(98)
where α2 = c
2∆t2
∆x2 . Simplifying further and employing the Euler formulae to simplify e
Ikn∆x + e−Ikn∆x =
2cos(k∆x):
G2 = 2G(1−α2)−1+2Gα2 cos(k∆x) (99)
G2 = 2G(1−α2(1− cos(k∆x)))−1. (100)
Noting that 1− cos(k∆x) = 2sin2 ( k∆x2 ) gives:
G2−2βG+1 = 0, (101)
with β = 1−2α2 sin2 ( k∆x2 ).
Step 3: Determination of any Necessary Stability Criteria
So we are nearly there. The solution to Eqn (101) is
G = β ±
√
β 2−1. (102)
As before, we are looking to find the limit of stability of the scheme which is given by |G| ≤ 1. Inspecting
the above equation, if we limit |β | ≤ 1, then |G| ≤ 1.
30
Thus we require |β | ≤ 1 giving:
−1≤ 1−2α2 sin2
(
k∆x
2
)
≤ 1 (103)
As 0≤ sin2(x)≤ 1, then the upper bound is automatically satisfied, the above equation is always less than
1. For the lower bound we must take the maximum value of sin2(x) which is 1, which gives:
−1 ≤ 1−2α2 (104)
α2 ≤ 1 (105)
Recalling the definition of α2 gives us:
c2∆t2
∆x2
≤ 1 (106)
This is our stability requirement. In practise, it is used to determine the time step size as follows:
∆t ≤ ∆x
c
(107)
In some notations you may see the time step size ∆t = CFL∆xc where CFL stands for Courant-Friedrichs-
Lewy and must be set to be less than or equal to 1. This is also called the CFL criteria.
A quick note on the physical interpretation of this time step limit. The quantity ∆xc is the time it takes a wave
of velocity c to travel a distance ∆x, which is the grid spacing. Thus our largest time step size is the time
it takes the wave to travel one cell. This connects our choice of discretisation, which uses information one
cell at most away from the cell begin computed (ui+1 and ui−1), and the stable time step size. This makes
sense, because if the time step size was larger, the result in the cell would be influence by information from
outside the region xi+1 ≤ x≤ xi+1 and so the result could not possibly be accurate. It turns out it is also not
stable.
7.3 Advanced - Special Case of a Uniform Grid
Above we implied that the scheme is actually fourth order in space and time on a uniform grid. Why is
this? Well, it’s down to a quirk in how the truncation errors for the two discretisations interact. Including
the leading order discretisation errors, our approximations are:
un+1i −2uni +un−1i
∆t2
= utt +
∆t2
12
utttt +O(∆t4),and c2
uni+1−2uni +uni−1
∆x2
= c2uxx+ c2
∆x2
12
uxxxx+O(∆x4),
(108)
where O(∆x4),O(∆t4) is there to indicate that the next error term is of order ∆x4,∆t4. For our choice of
discretisation, we solve:
un+1i −2uni +un−1i
∆t2
= c2
uni+1−2uni +uni−1
∆x2
(109)
which substituting in Eqn (108) gives:
31
utt +
∆t2
12
utttt +O(∆t4) = c2
uni+1−2uni +uni−1
∆x2
+ c2
∆x2
12
uxxxx+O(∆x4) (110)
Gather the error terms on the RHS:
utt = c2uxx+ c2
∆x2
12
uxxxx− ∆t
2
12
utttt +O(∆x4,∆t4) (111)
Now, we have shown that to be stable, ∆t = ∆xc . Substituting this into the above equation gives:
utt = c2uxx+
∆x2
12
(
c2uxxxx− 1c2 utttt
)
+O(∆x4,∆t4) (112)
So, the CFL condition has allowed us to write the errors in ∆t2 as an error in ∆x2. Now we look at the
governing equation, utt = c2uxx. Take the second derivative of the wave equation with respect to space
gives uttxx = c2uxxxx. Take the second derivative of the wave equation with respect to time which gives
utttt = c2uxxtt . Combine the two and we see that:
utttt = c2uxxtt = c2(c2uxxxx) = c4uxxxx (113)
Replace utttt in Eqn. (112) with c4uxxxx given by Eqn. (113) gives:
utt = c2uxx+
∆x2
12
(
c2uxxxx− c2uxxxx
)
+O(∆x4,∆t4) (114)
You can now see that the second order error term in brackets will disappear! Thus the actual error for
constant grid spacing is:
utt = c2uxx+O(∆x4,∆t4) (115)
where if you repeat the process for the fourth order error you would find out that it doesn’t disappear in the
same way that the second order error does.
7.4 Conclusions
In this chapter we’ve covered discretisations of the wave equation which are second order accurate in space
and time, stable within the CFL criteria, consistent with the governing equation. Furthermore, we detailed
how to circumvent the numerical difficulty of producing a starting condition.
32
8 Laplace and Poisson Equations
8.1 Introduction
uxx+uyy = p(x,y) (116)
So this equation is a little different to some of the previous ones, in that it is not described by the physical
phenomenon it describes (e.g. diffusion, or waves). When p(x,y) = 0 the equation is called the Laplace
equation, when p(x,y) 6= 0 it is the Poisson equation.
There is a close relationship between this equation and the heat equation. Consider the heat equation,
extended to two dimensional problems in space:
ut = c2(uxx+uyy). (117)
Now, if that equation reaches a steady state where u is no longer a function of time, then ut = 0. The above
equation then becomes the Laplace equation, uxx+uyy = 0. Thus the equation describes steady state diffu-
sion problems. An important ramification of this is that the initial condition is no longer needed, since
there is no dependence on time and so the solution depends only on the specified boundary conditions.
This equation is important in many applications including the computation of stresses, as you will see in
the finite element part of this course, solving for the flow of an incompressible fluid, electric potentials
and even grid generation for complicated calculations. Calculations of steady state heat distributions may
require a source term in the domain to represent a source, or sink, of heat (could be a heater for example).
In that case the equation to solve becomes the Poisson equation uxx+uyy = p(x,y), where p(x,y) represents
any sources or sinks for heat or electric potential for example.
Regarding the class of equation, as c1 and c3 are positive then the Laplace and Poisson equations are
Elliptic.
8.2 Analytical Solution
T0
T1
T2
T3 =
T0
0
0
0 +
0
T1
0
0 +
0
0
T2
0 +
0
0
0
T3
Figure 1: Solution approach for the Laplace Equation with multiple Dirichlet boundary conditions
First a note on the solution methodology for general problems with varying Dirichlet boundary conditions
(we do not cover Neumann boundary conditions for the Laplace/Poisson equations in this UoS). Figure 1
illustrates how you can take advantage of the fact that you can add solutions linearly (a sum) to get a new
solution of the governing equation. This is a particularly useful property of linear partial differential equa-
tions. So, to solve a Laplace equation problem with four different Dirichlet boundary conditions T0, ...,T3,
you need to solve four simpler problems where only one of the four boundary conditions is non-zero, and
all of the others taken to be zero. When you add the four separate solution together, you note that you
then satisfy the boundary conditions for the full problem with four Dirichlet boundary conditions, thus it
is the solution. For example, the sum of the upper boundary conditions will be T0 + 0+ 0+ 0 = T0. This
methodology holds even when the boundary conditions are not constant in space, for example T0(x).
A further simplification can be made. The Laplace equation has the same derivative in the x and y direction.
Thus we need to be able to solve only two problems, the first is when the non-zero boundary conditions
33
is located at the left or bottom of the domain (x = 0 or y = 0), the second is when the non-zero boundary
condition is located at the right or upper boundary (x = L or y = L). Furthermore, if we can solve these two
problems in the x-direction, then the y-direction solution will be the same solution but swapping x for y and
y for x. This is because the equation behaves the same in both directions. So, we just need to solve the two
cases shown in Figure 2.
0 0
T0
0
0
0
T2
0
−−−−−−−−−−→
x
x
y
Figure 2: Solution approach for the Laplace Equation with multiple Dirichlet boundary conditions, the left
hand diagram illustrates BC1, the right BC2
We are now very well practised at the solution of PDEs, and this equation follows the well worn path of
separation of variables. Here we assume that the domain size is Lx×Ly, and (0,0) is located in the bottom
left hand corner of the domain. Let’s follow the steps through here for the two cases shown in Figure 2:
1. Write out a clear definition of the problem. The equation, the initial conditions, the boundary con-
ditions and the parameters needed to define the specific problem (here only the distribution of the
source term matters).
2. Assume that the dependence of the solution on each dimension can be separated, such that u(x,y) =
F(x)G(y). This results in two equations, Fxx+ p2F = 0 and Gyy− p2G = 0. Here p is the separation
constant, not the p(x,y) source term in the Poisson equation.
3. Solve the ODE for F(x), which will result in a Fourier Series solution of the form Fn(x)=Ac1n cos(px)+
Bc1n sin(px). At this stage you can most often determine the value of constants A
c1
n or B
c1
n based on
the boundary condition at x = 0, and p may then be determined from the second boundary condition
by ensuring that the wavelength of the mode gives the desired answer at x = Lx. Given that both
F(0) = 0 and F(Lx) = 0 then Fn(x) = Bc1n sin(px) where p = npi/L.
4. Solve the ODE for G(y), given the value for p which has been determined in the previous step. Here
the solution differs from the previous two equations as the ODE for G is different. The solution is
Gn(y) = Ac2n e
py+Bc2n e
−py, (118)
where Ac2n , B
c2
n are constants. The boundary condition combinations for G(y) were are interested in
are BC1:
G(0) = 0, G(Ly) = T0, (119)
or BC2,
G(0) = T2, G(Ly) = 0. (120)
Solution for G(y) for BC1: G(0) = 0, G(Ly) = T0: First look at y = 0, G(0) = Ac2n +Bc2n = 0.
Thus, Bc2n =−Ac2n . This gives us:
Gn(y) = Ac2n
(
epy− e−py) (121)
34
This may be simplified by noting that 2sinh(x) = ex− e−x, giving:
Gn(y) = 2Ac2n sinh
(
npiy
Lx
)
(122)
Solution for G(y) for BC2: G(0) = T2, G(Ly) = 0: The simplest solution of this requires one
small step which is a shift of the reference axes as shown in Figure 3. We shift the y axis by Ly,
by defining q = y−Ly. The axes we will solve this problem in are thus x,q, where 0 ≤ x ≤ Lx and
−Ly ≤ q ≤ 0. (x,q) = (0,0) is located at the upper left corner of the domain, and (x,q) = (Lx,−Ly)
is the lower right hand corner.
(0,0)
(Lx,Ly)
T2
−−−−−−−−−−→
x
x
y
(0,0)
(Lx,−Ly)T2
−−−−−−−−−−→
x
x
q
=⇒
Figure 3: Change of axes to simplify the solution of problem BC2, the corner coordinates are shown in
brackets for the original and transformed coordinate system
In this coordinate system, the solution of F(x) does not change as x has not been transformed. For
G(y), this becomes G(q) with boundary conditions G(0) = 0 and G(−Ly) = T2. Starting from Eqn.
(118) we now gain:
Gn(0) = Ac2n +B
c2
n = 0 (123)
and so Bc2n =−Ac2n as for BC1, and:
Gn(q) = Ac2n
(
epq− e−pq) . (124)
This may be simplified by noting that 2sinh(x) = ex− e−x, giving:
Gn(q) = 2Ac2n sinh
(
npiq
Lx
)
. (125)
We need to remember that our lower boundary condition which has not yet been satisfied must also
be written in the new reference frame, thus u(x,q =−Ly) = T2.
5. Combine the two to form the full solution for a single mode (single value of n),
un(x,y or q) =
A
c3
n sinh
(
npiy
Lx
)
sin
(
npix
Lx
)
, for BC1
Ac3n sinh
(
npiq
Lx
)
sin
(
npix
Lx
)
, for BC2
(126)
Where we have combined two coefficients into Ac3n = 2A
c2
n B
c1
n . Note that a very common ‘gotcha’ in
the above expression is not using Lx for both of the terms - it often feels intuitive to have y/Ly in the
sinh term but that is not correct.
As usual, the full solution is a summation of all modes, giving:
u(x,y or q) = ∑
n=0,∞
un(x,y or q) =
∑n=0,∞A
c3
n sinh
(
npiy
Lx
)
sin
(
npix
Lx
)
, for BC1
∑n=0,∞Ac3n sinh
(
npiq
Lx
)
sin
(
npix
Lx
)
, for BC2
(127)
35
So we have one remaining coefficient, Ac1n to determine to give a specific solution. This is done by
considering the remaining non-zero boundary condition.
For BC1:
u(x,Ly) = T0 = ∑
n=0,∞
Ac3n sinh
(
npiLy
Lx
)
sin
(
npix
Lx
)
(128)
Initially, this may not look useful. However, lets compare it to a Fourier Sine series in x:
∑
n=0,∞
Ac3n sinh
(
npiLy
Lx
)
sin
(
npix
Lx
)
= ∑
n=0,∞
Bn sin
(
npix
Lx
)
, (129)
we can see that the two are equivalent if Bn = Ac3n sinh
(
npiLy
Lx
)
. Furthermore, we know that if T0 =
∑n=0,∞Bn sin
(
npix
Lx
)
, then the coefficients Bn are the coefficients of the Fourier Series representing the
boundary condition at T (x,Ly) = T0(x). Given this distribution we can compute the Bn, and given the
Bn we can then compute Ac3n . Equivalently, we can just replace A
c3
n with Bn giving for BC1:
u(x,y) = ∑
n=0,∞
Bn
sinh
(
npiy
Lx
)
sinh
(
npiLy
Lx
) sin(npix
Lx
)
(130)
The lecture video has three examples of how to determine Bn for three different distributions of
u(x,Ly). Note in the video lengths a = Lx and b = Ly.
For BC2: Let’s follow the same process, except now we are looking at the boundary condition at
q =−Ly, giving:
u(x,−Ly) = T0 = ∑
n=0,∞
Ac3n sinh
(−npiLy
Lx
)
sin
(
npix
Lx
)
(131)
Equating with a Fourier Series once more gives:
Ac3n sinh
(−npiLy
Lx
)
sin
(
npix
Lx
)
= Bn sin
(
npix
Lx
)
, (132)
and so Bn = Ac3n sinh
(
−npiLyLx
)
, where Bn are the Fourier coefficients of the distribution of T (x,q =
−Ly). Rearranging to eliminate Ac3n :
u(x,q) = ∑
n=0,∞
Bn
sinh
(
npiq
Lx
)
sinh
(−npiLy
Lx
) sin(npix
Lx
)
(133)
Noting that q = y−Ly, we can return the equation back to the original (x,y) reference frame, giving
for BC2:
u(x,y) = ∑
n=0,∞
Bn
sinh
(
npi(y−Ly)
Lx
)
sinh
(−npiLy
Lx
) sin(npix
Lx
)
(134)
An example application of this is given at the end of the lecture recording. Note that the solutions
to the Week 8 tutorial give a methodology to solve BC2 problems without a change of frame of
reference, which ends in the same solution form.
36
Following the methodology illustrated in Figure 1 you are now able to solve any combination of Dirichlet
boundary conditions on the four boundaries of a two dimensional problem.
37
(i, j) = (0,0)
(0,1)
(0,2)
(0,3)
(0,4)
(0,5)
(1,0) (2,0) (3,0) (4,0) (5,0)
u(i, j)u(i−1, j) u(i+1, j)
u(i, j+1)
u(i, j−1)
Figure 4: Labelling of the grid in the solution of the Laplace or Poisson Equation, highlighting the stencil
of the discretisation in blue
9 Numerical Solution of the Laplace and Poisson Equation
9.1 Introduction
Here we will deal with the solution of the Laplace and Poisson equation only for steady state problems. If
you do need to look at the time dependent behaviour, the nature of the equation changes and it becomes a
two-dimensional heat equation - you can follow the discretisation method outlined in Section 5 in that case.
There are two common approaches to solving the Laplace/Poisson equations, the first is using finite differ-
ences following a similar route to our previous discretisations, borrowing elements of the representation of
the spatial derivative in the heat equation, and methods used to solve the Implicit form of the heat equation
discretisation outlined in Section 5.2. The second approach uses Finite Elements which we will get onto
later in the course.
9.2 Second Order Central Differences
The equation were are discretising is uxx + uyy = p(x,y), solved on a rectangular domain of size Lx×Ly.
The domain is split into nx cells in the x-direction such that x = inx where i is the index in the x-direction,
and split into ny points in the y-direction such that y = jny where j is the index in the y-direction. Thus
a given value u(xi,y j) may be written as ui, j. Grid spacing is assumed to be uniform in both directions,
∆x = xi+1, j− xi, j and ∆y = yi, j+1− yi, j. A schematic of the mesh is shown in Figure 4.
Applying second order central differences to the second derivatives gives
ui+1, j−2ui, j +ui−1, j
∆x2︸ ︷︷ ︸
uxx
+
ui, j+1−2ui, j +ui, j−1
∆y2︸ ︷︷ ︸
uyy
= p(xi,y j) (135)
In this course we will only deal with square grids, such that ∆x = ∆y= h. The above equation simplifies to:
38
ui+1, j−4ui, j +ui−1, j +ui, j+1+ui, j−1 = p(xi,y j)h2, (136)
or,
ui, j =
1
4
(
ui+1, j +ui−1, j +ui, j+1+ui, j−1− p(xi,y j)h2
)
. (137)
The computational points used in the algorithm, or the stencil, is a cross shape. ui, j cannot be determined
directly as it depends directly on the values of u at all other points. If you consider a full domain, then
eventually the stencils will include the boundaries, and hence the boundary conditions. Thus every single
point must be computed at once to gain the solution - i.e. you will form a system of simultaneous equations.
Just as with the implicit solution of the heat equation, the equations may be written in matrix form as
follows:
AU = P+B. (138)
Here, A is a matrix of the coefficients for all unknown points, U is the vector of unknowns, P is the vector
of source or sink terms and B the vector of known boundary values. For a mesh of 6×6 shown in Figure 4,
the unknown points which need to be solved for are 1≤ i≤ 4 and 1≤ j ≤ 4, since the i = 0 and i = 5 row
are the left and right boundary conditions and the j = 0 and j = 5 rows are the lower and upper boundary
conditions. To assist in constructing the matrices, the vector of unknowns U is constructed in a specific
order, it starts with u(1,1) and proceeds through i then j - i.e. it is listed from bottom left hand corner to
the top right hand corner in Figure 4.
In the case in Figure 4, the matrix A has size 16×16, and vectors U and P are given by:
A =

−4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0 0
0 1 −4 1 0 0 1 0 0 0 0 0 0 0 0 0
0 0 1 −4 0 0 0 1 0 0 0 0 0 0 0 0
1 0 0 0 −4 1 0 0 1 0 0 0 0 0 0 0
0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0 0
0 0 1 0 0 1 −4 1 0 0 1 0 0 0 0 0
0 0 0 1 0 0 1 −4 0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0 −4 1 0 0 1 0 0 0
0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0 0
0 0 0 0 0 0 1 0 0 1 −4 1 0 0 1 0
0 0 0 0 0 0 0 1 0 0 1 −4 0 0 0 1
0 0 0 0 0 0 0 0 1 0 0 0 −4 1 0 0
0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1 0
0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4 1
0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −4

(139)
39
U =

u(1,1)
u(2,1)
u(3,1)
u(4,1)
u(1,2)
u(2,2)
u(3,2)
u(4,2)
u(1,3)
u(2,3)
u(3,3)
u(4,3)
u(1,4)
u(2,4)
u(3,4)
u(4,4)

,P =

p(1,1)h2
p(2,1)h2
p(3,1)h2
p(4,1)h2
p(1,2)h2
p(2,2)h2
p(3,2)h2
p(4,2)h2
p(1,3)h2
p(2,3)h2
p(3,3)h2
p(4,3)h2
p(1,4)h2
p(2,4)h2
p(3,4)h2
p(4,4)h2

,B =−

u(0,1)+u(1,0)
u(2,0)
u(3,0)
u(4,0)+u(5,1)
u(0,2)
0
0
u(5,2)
u(0,3)
0
0
u(5,3)
u(0,4)+u(1,5)
u(2,5)
u(3,5)
u(4,5)+u(5,4)

. (140)
To solve for the vector of unknown variables U , the matrix A must be inverted. Thankfully, it has a regular
structure and as so efficient methods can be employed to invert the matrix, hence enabling the solution of
the linear algebra problem Eqn. (138). In Matlab this is as simple as a single line command - so most of
the effort in coding is actually in constructing the matrix in the above form. This can be automated with
some thought however. The vertical and horizontal lines in matrix A represent the start and end of each
row in the i direction, where the ‘gaps’ in the middle three diagonals (the pattern of −14− 1) are caused
by the start and end of each row respectively, where the stencil will refer to the left and right boundary
conditions (i+1 and i−1). The offset diagonals of value −1 represent the upper and lower elements in the
discretisation ( j+ 1 and j− 1) which in the first and last row may refer to boundary cells. All boundary
cells are represented as known values on the right hand side of the equation in the vector B.
Note that the source term p(x,y) is simply discretised by using the value of p(xi,y j).
9.3 Conclusions
This completes the description of a second order in space scheme to solve the Laplace or Poisson equations.
Here we do not expect you to memorise the form of the matrix, more that you should be familiar with the
method used to discretise the equations and construct the matrix. It is also a useful example of how to deal
with source terms in the discretisation of the governing equations.
40
10 Understanding PDEs
10.1 Introduction
We have already looked at a broad method to classify simple second order PDEs as elliptic, hyperbolic
or parabolic. This section details how you could go about assessing the behaviour of a set of governing
equations which you may have developed for a particular engineering problem of interest and so may not
be in a form which is amenable to that method. One example is the PDE ut + cux = 0, with c a constant, is
that system diffusive or oscillatory in time for example? At the moment you are not able to answer that in
a systematic sense. To provide an answer, we adapt our current tools to this problem in the following two
steps:
1. Assume a form of the solution to the governing equation. In undertaking von Neumann stability
analysis, we assumed that the solution to the govering equation was made up of a Fourier series
in space, and a variation in time - where the Fourier series in space was expressed as a complex
exponential. Here we adopt a similar approach, where it is assumed that the solution is a Fourier
series in space and a complex exponential in time, i.e.:
u(x, t) =
∞
∑
ω=−∞
uˆ(ω)eIωxeIβ t , (141)
where uˆ(ω) are the Fourier coefficients (complex) representing the amplitudes of modes in space and
time.
Firstly note that the spatial mode ω = pinL (as in Section 5.1.2), thus ω is real. However, β gives the
temporal behaviour of the governing equation and may be either real, or imaginary. If β is real, then
the exponential will have an imaginary power leading to a Fourier Series solution in time. If β is
complex, then the exponent will have a real power thus leading to either exponential growth in time,
or decay. If β = 0 then there is no change. Thus this assumed form can capture all of the potential
behaviours of the governing equation. We just need to determine β .
A quick note - examining a single mode initial condition (only mode ωn), u(x, t) = uˆ(ωn)eIωnxeIβnt .
At t = 0 we gain un(x,0) = uˆ(ωn)eIωnx, and this will be used to simplify our results.
2. Substitute into the governing equation, and determine β Substitute one mode from Eqn. (141)
into the governing equation, and simplify until you gain an expression relating the unknown β to the
known ω . As with our observations of the heat equation, wave equation and Laplace equation, if you
understand one mode, you understand all of the modes.
Here are two practical examples, more are given in the lectures.
Ex 1 The heat equation ut = c2uxx. First note that uxx =−ω2n uˆneIωnxeIβnt and ut = iβnuˆneIωnxeIβnt . Substitute
them into the governing equation, giving:
Iβ uˆneIωnxeIβnt =−c2ω2n uˆneIωnxeIβnt (142)
Simplifying:
βn =−c
2ω2n
i
= Ic2ω2n . (143)
Thus β is positive imaginary. Insert that back into the assumed solution Eqn. (141) giving:
41
un(x, t) = uˆneIωnxe−c
2ω2n t (144)
Noting that at t = 0, un(x,0) = uˆeIωnx, thus un(x, t) = un(x,0)e−c
2ω2n t . As c2ω2n is positive and real, then the
predicted behaviour is that each mode will decrease exponentially in amplitude as time progresses. This is
exactly what we expect for the heat equation, and so this equation is parabolic.
Ex 2 Now how about the new equation, ut +cux = 0? We note that ux = Iω uˆneIωnxeIβnt , ut = Iβ uˆneIωnxeIβnt .
Inserting this into the governing equation gives:
IβnuˆneIωnxeIβnt + cIωnuˆneIωnxeIβnt (145)
Simplifying:
βn =−cωn (146)
Giving
un(x, t) = uˆneIωnxe−Icωnt (147)
Once more, setting t = 0 gives un(x,0) = uˆneIωnx and so un(x, t) = un(x,0)e−Icωnt . As ω and c are real,
then the exponential is raised to a complex power which gives an oscillatory behaviour in time. Thus this
equation behaves like the wave equation and is hyperbolic.
There are some equations which are just badly behaved. For example, the equation ut =−c2uxx represents a
diffusive flow with negative diffusion. If you follow through the above methodology you should see that the
solution amplitude increases exponentially in time, thus any perturbation would rapidly grow - and there is
no bound to this growth (it would head quickly to infinity!).
42
11 Fourier Integrals, Fourier Transform and Laplace Transform
11.1 Introduction
We have seen already how to solve problems where there are boundary conditions placed at a well defined
position. However, a number of practical problems do not have easily definable boundary conditions.
Consider, for example, the diffusion of a pollutant (perhaps on oil spill) in the middle of the ocean. If
we needed to predict the concentration of the pollutant at a specific location, we would want to solve the
diffusion equation, ut = c2uxx. However, to determine the modes present and amplitudes, we would need to
specify boundary conditions. Where would these be? Ideally we would want to specify that u(x→∞, t) = 0
and u(x→−∞, t) = 0, where u(x, t) is the concentration. A second example is the simulation of pressure
waves in a long pipe - imagine if a valve located at x = 0 closes in a long pipe, sending a pressure wave
down the pipe. What boundary conditions would we use in the pipe to solve for the evolution of that wave?
Again, ideally we would want to specify that u(x= 0, t) = p(0, t) and u(x→∞, t) = p(x→∞, t) = 0, where
p(x, t) is the pressure fluctuation which may be represented by a hyperbolic equation.
These are examples of problems which may be said to be in an infinite domain (e.g. −∞≤ x≤ ∞) or semi-
infinte domain (e.g. 0≤ x≤∞). We cannot apply a Fourier Series as the Fourier Series can only be applied
to finite domains. Thus we need more mathematical tools to deal with these new boundary conditions.
These new mathematical tools are the Fourier Integral, Fourier Transform, and Laplace Transforms. The last
two are perhaps the most important mathematical transforms in Engineering. The first leads to extremely
useful mathematical solutions.
11.2 Fourier Integral
11.2.1 Rationale and Definition
As detailed above, the Fourier Series solutions only hold for finite length domains, for example the Fourier
coefficients are given by:
f (x) = ∑
n=0,∞
(an cos(ωnx)+bn sin(ωnx)) (148)
where ωn = npi/L and an and bn are given by:
a0 =
2
L
∫ L
−L
f (x)dx, an =
1
L
∫ L
−L
f (x)cos(ωnx)dx, bn =
1
L
∫ L
−L
f (x)sin(ωnx)dx (149)
We now explore what happens to the Fourier Coefficients a0, an and bn for a fixed initial distribution u(x,0)
which satisfies the requirement to have a zero value at infinity. Trends towards an infinite domain size can
be understood by looking at L→ ∞. In the lecture slides we assume a top hat distribution over a finite
range of x, so what happens to the amplitudes of each individual mode (called the amplitude spectrum) as L
increases? There are two main impacts, the first is that all of the above amplitudes vary as 1/L, the second
is that the individual modes move closer together. The latter effect can be seen more clearly if you look at a
difference between two modes, given ω = npi/L then the difference ωn+1−ωn = pi/L, thus this decreases
proportional to 1/L. If the amplitude of a given mode number is plotted against mode number, the shape of
the resultant curve of an or bn (dependent on whether it is an even or add function) can be seen to become
more and more smooth.
Now, lets use these observations to formalise a Fourier Integral. Let’s first insert Eqn. (149) into (148):
43
f (x) =
2
L
∫ L
−L
f (x)dx+
1
L ∑n=1,∞
(
cos(ωnx)
∫ L
−L
f (v)cos(ωnv)dv+ sin(ωnx)
∫ L
−L
f (v)sin(ωnv)dv
)
, (150)
where v has been used in the integrals in the Fourier coefficients instead of x so that the variable which is
being integrated is clear. Given ∆ω = ωn+1−ωn = pi/L, we can substitute in for 1/L giving:
f (x) =
2
L
∫ L
−L
f (x)dx+
1
pi ∑n=1,∞
(
cos(ωnx)∆ω
∫ L
−L
f (v)cos(ωnv)dv+ sin(ωnx)∆ω
∫ L
−L
f (v)sin(ωnv)dv
)
,
(151)
This holds for large L, but now let L→ ∞, and assuming that the function has a finite integral we see that
2
L
∫ L
−L f (x)dx→ 0. Next, we have two terms which are the sum of the spectral amplitudes multiplied by the
∆ω:
f (x) =
1
pi ∑n=1,∞
(
cos(ωnx)
∫ L
−L
f (v)cos(ωnv)dv+ sin(ωnx)
∫ L
−L
f (v)sin(ωnv)dv
)
∆ω, (152)
Given that the spectral amplitudes becomes progressively smoother we replace ∑n=1,∞(.)∆ω →
∫ ∞
0 (.)dω .
Mathematically there are some caveats in this procedure, but it is a valid assumption for most engineering
applications -particularly those tackled in this course. Using these two tools we gain:
f (x) =
1
pi
∫ ∞
0
(
cos(ωx)
∫ ∞
−∞
f (v)cos(ωv)dv+ sin(ωx)
∫ ∞
−∞
f (v)sin(ωv)dv
)
dω. (153)
Eqn. (153) represents the Fourier Integral representation of f (x), and may now be applied to non-periodic
functions where ω no longer represents integer modes. It could be seen as the ∆ω→ 0 limit of a numerical
integration of the spectral amplitudes which are plotted in the lecture slides.
Eqn. (153) can be simplified for even functions to be:
f (x) =
∫ ∞
0
A(ω)cos(ωx)dω, A(ω) =
2
pi
∫ ∞
0
f (v)cos(ωv)dv, (154)
and for odd functions,
f (x) =
∫ ∞
0
B(ω)sin(ωx)dω, B(ω) =
2
pi
∫ ∞
0
f (v)sin(ωv)dv. (155)
These are two very useful functions which allow the representation of non-periodic functions. A final useful
form is the following, gained through trigonometric relations:
f (x) =
1
pi
∫ ∞
0
[∫ ∞
−∞
f (v)cos(ωx−ωv)dv
]
dω. (156)
11.3 Fourier Integrals in the Solution of PDEs
We are now going to use this new mathematical tool to solve a heat diffusion problem in an infinite domain.
The solution steps are:
44
1. Solve the governing equations through the separation of variables method for a single mode, but do
not apply the boundary conditions (leaving it in general form)
2. Instead of summing the individual modes in the Fourier Series, replace that step with the Fourier
Integral representation of the initial condition
3. Use the initial condition to determine A(ω) and B(ω).
4. Evaluate the resultant integral to gain the solution
So lets undertake an example of this as applied to predict the temperature distribution as a function of time
in an infinitely long wire. The equation to solve is Tt = c2Txx, with an initial condition of T (x,0) = f (x) = T0
for |x|< 1, and T (|x|> 1,0) = 0. This forms a top hat centred around x = 0. The solution steps are:
Step 1 Separation of Variables Following separation of variables as done earlier for the heat equation the
solution for a single mode is:
Tn(x, t) = (An cos(px)+Bn sin(px))e−c
2 p2t (157)
Step 2 & 3 Use the Fourier Integral Instead of summing the modes, apply a Fourier Integral representation
(note that p=ω here) using Eqn. (156), which has already accounted for the initial conditions as the Fourier
Integral coefficients A(ω) and B(ω) were determined using them:
T (x, t) =
1
pi
∫ ∞
0
[∫ ∞
−∞
e−c
2 p2t f (v)cos(px− pv)dv
]
d p. (158)
Note that the exponential is inside the integral and cannot be moved outside as p features in the power of
the exponential. This can be rearranged into the following form:
T (x, t) =
1
pi
∫ ∞
−∞
f (v)
[∫ ∞
0
e−c
2 p2t cos(px− pv)d p
]
dv. (159)
This equation is very useful - it holds for any initial condition f (v) = f (x).
Step 4: Solve the integral We now employ the identity
∫ ∞
0 e
−s2 cos(2bs)ds =
√
pi
2 e
−b2 . The integral in the
square brackets above can be written in that form given s2 = c2 p2t and 2bs = px− pv, from which p = sc√t
and b = x−v2c√t . Note that the integral is in terms of p, thus we need to substitute ds = c
√
td p. Thus:
∫ ∞
0
e−c
2 p2t cos(px− pv)d p =
∫ ∞
0
e−s
2
cos(2bs)
ds
c
√
t
=
√
pi
2c
√
t
e−b
2
=
√
pi
2c
√
t
exp
(
−(x− v)
2
4c2t
)
(160)
Inserting this back into Eqn. (159) we now have:
T (x, t) =
1
2c
√
pit
∫ ∞
−∞
f (v)exp
(
−(x− v)
2
4c2t
)
dv, (161)
Given the initial condition T (x,0) = f (x) = T0 for |x|< 1, which is only non-zero between −1≤ x≤ 1, the
integration can be simplified to
45
T (x, t) =
T0
2c
√
pit
∫ 1
−1
exp
(
−(x− v)
2
4c2t
)
dv. (162)
Denoting z2 = (x−v)
2
4c2t , hence z=−
(x−v)
2c
√
t (one of two solutions - if you take the positive root you get the same
answer in the end). Thus v = x+2c
√
tz and hence dv = 2c
√
tdz. Inserting these into Eqn. 162 gives:
T (x, t) =
T0√
pi
∫ (1−x)
2c
√
t
− (x+1)2c√t
e−z
2
dz, (163)
Note that the limits of the integral were originally from v = −1 to v = 1, inserting these into z = − (x−v)2c√t
gives the new limits in the integral in z. It can be rewritten as:
T (x, t) =
T0√
pi
[∫ 0
− (x+1)2c√t
e−z
2
dz+
∫ (1−x)
2c
√
t
−0
e−z
2
dz
]
, (164)
We now swap the limits of the integration around in the first integral, giving:
∫ 0
− (x+1)2c√t
e−z
2
dz =−
∫ −(x+1)
2c
√
t
0
e−z
2
dz, (165)
Now the error function is defined as:
erf(x) =
2√
pi
∫ x
0
e−z
2
dz, (166)
Thus Eqn. 164 becomes
T (x, t) =
T0√
pi
[
−erf
(
−(x+1)
2c
√
t
)
+ erf
(
(1− x)
2c
√
t
)]√
pi
2
, (167)
Noting that the error function is an odd function, i.e. erf(x) = −erf(−x) means that we can simplify the
above equation. So our Eqn. 167 becomes finally:
T (x, t) =
T0
2
[
erf
(
(x+1)
2c
√
t
)
+ erf
(
(1− x)
2c
√
t
)]
, (168)
Remember to do a sanity check - at x = 0 and t = 0 this should be T0. As erf(∞) = 1, then the above
expression simplifies to 2T02 = T0.
The Error Function: As discussed in the lecture slides the Error Function is a fantastically useful function.
It not only allows you to go smoothly from one function to another, thus serving as an excellent blending
function, but also represents physically the diffusion process exactly. By combining two error functions as
in Eqn. (168) you can generate profiles which start initially as a sharp top hat, with none of the ringing
associated with undersampled Fourier Series, through to a late time very smooth diffusve peak. A very
useful function and it is well worth spending time in Matlab to play with the controlling parameters to
make sure that you understand how it behaves. Start with a simple plot of 12 (erf(c1x)+1) for −5 ≤ x ≤ 5
with c1 = 1, then explore what happens as you vary constant c1, how you make a sharper or more diffuse
looking profile.
46
11.4 The Fourier Transform
11.4.1 Introduction
In Section 11.2.1 we examined the amplitude spectrum, which is a compact method of understanding the
variation of the amplitude of each mode present in a solution as the mode number changes. This indicates
to us that there is a profoundly different and incredibly valuable way of understanding and interpreting
a function f (x) rather than simply examining it’s shape in physical space (x). This alternative way of
analysing a function is to calculate and plot the spectral amplitudes of the function. In doing that, it is
possible to identify the dominant frequencies or mode numbers in a given function or signal. The
word ‘signal’ is often employed since this approach is often used to understand the nature of time dependent
functions which have been recorded from experiments, or perhaps received from a broadcast, hence the
terminology ‘signals’. The Fourier Transform is a standardised method of computing and presenting
amplitude spectra for a given signal or function.
Why is this an important method of analysis?
Typically a given signal is made up of many different signals from many different sources. Each source has
a characteristic signal, and when breaking it down to individual components (say a single part in a machine,
such as a single valve in the heart, or a single fan blade on a jet engine) that component will have one or
many distinct characteristic frequencies. These characteristic frequencies may be, for example, the first and
harmonic modes identified in the wave equation. They may be the rotation rate of the device, or the pulse
rate, the vortex shedding frequency, or the broadcast signal for example. When viewed as a function of
space or time, all of the signals are mixed together to make a composite signal, as happens when summing
the Fourier Series. It is then very difficult to see the impact of individual components of a signal on the
composite signal f (x) or f (t). Furthermore, there may be substantial noise in the measurement which may
obscure the interesting part of the signal.
The amplitude spectrum fˆ (ω) of the composite signal allows the engineer to see exactly which frequencies
are present in the signal and which are the dominant or otherwise important frequencies. For the important
frequencies, the frequency/mode number may be used to identify the sources of that specific mode, enabling
mitigation, or further post-processing. Any noise in the signal in physical space (x or t) has a constant
amplitude spectrum and so does not contribute to any peaks.
A final reason for its importance is that plotting the amplitude spectrum can be a substantially more compact
method of visualising the data than the physical data itself. For example, in the aeroacoustics problem in
the lecture slides, the amplitude spectrum was computed for several million data points. Over that huge
number of points, there are some wavelengths which are several hundred thousand times the wavelength
of the shortest wavelengths in the signal. When viewing the data, it is impossible to simultaneously view
the short and long wavelengths together, and the long wavelengths are not visible when viewing the whole
signal because they are obscured by the short wavelengths.
An example - predictive maintenance
In high value industrial applications it is very important to minimise the down time of each system. When
running a fleet of helicopters, for example, you may be able to undertake regular scheduled maintenance,
but as you should appreciate this regular maintenance schedule is determined from the mean expected useful
lifetime of each component, not the actual lifetime of a single component. There can be quite a substantial
variation in how long a single component lasts dependent on uncertainty in manufacturing, or variations
in usage (consider sports car vs. long distance commute). To improve availability of systems, and cost-
effective maintenance, it may be better to replace difficult or high value components only when needed. Yet
- you don’t want to have a costly unscheduled breakdown.
The solution in may be made relatively straightforward with careful characterisation of the system of interest
47
and Fourier Transforms. For helicopter gearboxes, for example, an accelerometer is attached to the side of
the engine. This accelerometer measures the accelerations of the gearbox at the level of vibrations (not
acceleration of the helicopter itself). The vibrations are caused by the natural behaviour of the gearbox -
i.e. the rotation of various parts in the gearbox or opening and closing of valves, for example. The first step
is to measure the accelerations from a ’healthy’ undamaged gearbox and calculate the amplitude spectrum
using the Fourier Transform. Then, monitor the gearbox in use by calculating the amplitude spectrum at
periodic intervals. As the gearbox ages and wears, the amplitude spectrum will change. The change in each
individual peak in the spectrum identifies the wear of an individual component. Significant shifts in the
amplitude of location of the peak, or new peaks appearing in the amplitude spectrum indicate components
about to fail, or which have failed.
Given an understanding of which component produces which peak in the amplitude spectrum means that
the user knows exactly which component is broken or worn, and can then schedule a replacement. All
of this is done without even opening up the gearbox, simply by measuring vibrations and comparing the
amplitude spectra.
Such approaches are routinely used to diagnose medical problems, particularly with the heart or brain.
Changes in the frequency content of measured activity may point clearly to the underlying pathology.
11.4.2 Formalisation
The Fourier Transform of function f (x) is denoted by fˆ (ω). For even functions a cosine transform may be
used which is fˆc(ω) and for odd functions a sine transform may be used which is denoted fˆs(ω). These
transforms are proportional to the coefficients in the Fourier Integral A(ω) and B(ω) in Eqns. (154) and
(155). For even and odd functions respectively, the Fourier Cosine and Sine Transforms are:
fˆc(ω) =
√
pi
2
A(ω) =
√
2
pi
∫ ∞
0
f (x)cos(ωx)dx (169)
fˆs(ω) =
√
pi
2
B(ω) =
√
2
pi
∫ ∞
0
f (x)sin(ωx)dx (170)
From these definitions the Fourier Cosine and Sine Inverse Transforms are:
f (x) =
√
2
pi
∫ ∞
0
fˆc(ω)cos(ωx)dω, (171)
f (x) =
√
2
pi
∫ ∞
0
fˆs(ω)sin(ωx)dω. (172)
Notation may differ in some fields, for example fˆc(ω) may be denoted by Fc( f ). The combination of the
original and transformed signals are called ‘transform pairs’, for example f (x) and fˆc(ω) are Fourier Cosine
transform pairs. You should familiarise yourself with the transform of fundamental forms of f (x) such as
sine and cosine waves, delta functions, box/top hat and discrete pulses (comb).
The transform can also be written in terms of complex exponentials following the methodology outlined in
the lecture slides. This results in the definition of the Fourier Transform applicable to both odd and even
signals:
f (x) =
√
1
2pi
∫ ∞
−∞
fˆ (ω)expIωx dω, fˆ (ω) =
√
1
2pi
∫ ∞
−∞
f (x)exp−Iωx dx (173)
48
Note that as with the Fourier Integral, the signal must have a finite integral (it is not periodic).
11.4.3 Application to Signal Processing
In most practical cases we deal with a finite signal, so how is a Fourier Transform employed? There are
two tricks employed to deal with a finite length signal which is sampled (measured) at a time interval ∆t.
The first is to assume that it is periodic, i.e. if the signal is extended from the sampled range (e.g. 0≤ t ≤ T
to −∞ ≤ t ≤ ∞), the only modes present in the signal are those that fit fully in the the original sampling
range, i.e. the lowest mode present is 2pi/T . The second is to examine the signal in terms of it’s energy per
unit time or space, which is the power of the signal, meaning that we only want to compute the power in a
given temporal or spatial range (e.g. 0≤ t ≤ T ). This removes the problem of the integral above becoming
infinite. Given these observations, the complex form of the Fourier Series may be used to compute the
spectral amplitudes since the signal is now periodic and bounded in time or space, giving
f (t) =
∞
∑
n=−∞
cn expIωt , cn =
1
T
∫ T
0
f (t)exp−Iωt dt, (174)
where ω = n2pi/T and n is an integer −∞≤ n≤ ∞.
In the case of a signal which is measured for a time period T at a sampling rate ∆t, the modes are now
discrete and range from the lowest wavelength T and shortest wavelength 2∆t. Thus the fundamental
frequency is the first whole mode which fits in the sample range T , i.e. ω = 2pi/T . This is because the
signal is assumed to be periodic, so only full modes can be selected, as opposed to some previous cases
where mode 1 has been pi/T .
Given this, Eqn. 174 can be employed to compute the amplitude spectrum, restricted to modes n = 1 to
n = T/2∆t given ω = n2pi/T .
The power spectral density, often informally denoted the power spectrum, plots the power in a given signal
for each mode in the problem. It is defined so that the integration of the power spectrum gives the overall
power in the signal, i.e.:
∫ ∞
0
P(ω)dω = lim
T→∞
1
T
∫ T
0
f (t)2dt, (175)
which leads to a second definition of the Fourier Transform most commonly employed in computing power
spectra:
fˆ p(ω) =
√
1
T
∫ T
0
f (t)exp−Iωt dt (176)
The factor of
√
1
T arises to convert from an energy to power spectrum. The power spectrum is computed
via:
P(ω) = | fˆ p(ω)|2 = fˆ p(ω) fˆ p,∗(ω), (177)
where fˆ p,∗(ω) is the complex conjugate of fˆ p(ω). This is extremely useful as it allows you to see directly
which modes are contributing the most to the overall power in the signal, which is often the most important
49
quantity. This is termed the ‘Power Spectral Density’, which in many ways makes more sense as the units
are are actually power/(Units of ω).
In time dependent problems, engineers often refer to f (t) as the time domain representation and fˆ (ω) as
the frequency domain representation.
In space dependent problems, engineers often refer to f (x) as the space domain representation and fˆ (ω)
as the wavenumber domain representation.
11.4.4 Important Properties
There are several important properties of the Fourier Transform:
1. Linearity: F{a f +bg}= aF{ f}+bF{g}
2. Derivatives: F{ f ′}= IωF{ f},F{ f ′′}=−ω2F{ f}
3. Convolution theorem (for computing correlations): Given the convolution of f and g defined by
f ? g =
∫ ∞
−∞ f (x− p)g(p)d p, F{ f ? g} =
√
2piF{ f}F{g}. This is a very important result, as it
is much faster computationally to compute the Fourier Transforms, then multiply then undertake an
inverse Fourier Transform that to directly evaluate the integral in the correlation.
11.4.5 The Fast Fourier Transform
The most common method of computing the Fourier Transform is called the Fast Fourier Transform or
FFT. You will come across this terminology frequently during your degree. There are some important
observations to be made about the practical use of the FFT in signal processing:
1. Shannon’s sampling theory tells us that for a data sampling rate fS, the highest frequency we can
properly capture (without aliasing) is fNyq = fS/2, called the Nyquist frequency. This is equivalent
to a wavelength which has a period of 2∆t in time domain data, or 2∆x in space domain data, where
∆x or ∆t is the sampling interval. There is no point trying to represent frequencies above this because
our data sampling rate cannot support it.
2. The lowest frequency is mode 1 based on the full signal length, i.e. a frequency of 2pi/T or 2pi/L
where T or L is the temporal or spatial domain size. However, note that the spectrum also captures
mode zero, which is the mean. Often this is left out on plots since power spectra are typically plotted
on log-log axes.
3. Therefore, when assembling the frequency grid points for the FFT, most often the grid goes from
f = 2pi/T to f = fNyq.
4. The output of the FFT is complex-valued. To get the amplitude (real) spectrum we simply take the
absolute value of the output.
5. The FFT is also two-sided – that is, it outputs both the positive and negative frequency spectrum.
Often we just want the positive (single-sided) spectrum and wish to ignore the symmetric negative
frequencies. Therefore, we multiply the single-sided amplitude spectrum by 2 to normalise for this.
This is possible as most often we are interested in the power spectrum rather than the coefficients
themselves, and the power is the absolute value squared. The power spectrum is thus symmetric, and
so a simple doubling of the positive side of the spectrum accounts for the power on the negative side
(ω < 0).
50
In Matlab:
1. The FFT is computed using the Matlab function fft(y, NFFT), where y is the function you want to
transform and NFFT is the number of points being transformed.
2. NFFT must be a power of 2 to work most efficiently. So we have to specify the next power of 2 from
our number of data points. Notice that we can achieve this very easily using the Matlab nextpow2
function. It is possible to transform a non-power of 2 set of data but this is more expensive.
11.4.6 The Fourier Transform Solution to the Heat Equation
The Fourier Transform may be used to solve PDEs, however in this course there is no advantage of using
the Fourier Transform over the Fourier Integral solution, thus you may use either approach.
The first step is to apply the Fourier Transform to one dimension of the governing equation, which in this
case is Tt = c2Txx. Taking the transform of the heat equation with respect to x gives:
fˆ (Tt) =−c2ω2Tˆ . (178)
Interchanging the derivative with respect to time and the transform gives an ODE for Tˆ ,
Tˆt =−c2ω2Tˆ . (179)
The solution is Tˆ (ω, t) = c(ω)exp−c2ω2t . At time t = 0, T (x, t) = f (x), thus Tˆ (ω,0) = c(ω) = fˆ (ω).
Therefore Tˆ (ω, t) = fˆ (ω)exp−c2ω2t . Using the inverse transform (Eqn. 173) gives:
T (x, t) =
1√
2pi
∫ ∞
−∞
fˆ (w)e−c
2w2teiwxdw (180)
Note that the Fourier Transform of the initial condition f (x) is
fˆ (w) =
1√
2pi
∫ ∞
−∞
f (v)e−iwvdv (181)
Inserting this into Eqn. (180).
T (x, t) =
1
2pi
∫ ∞
−∞
f (v)
[∫ ∞
−∞
e−c
2w2teiw(x−v)dw
]
dv (182)
Using the Euler Formula:
e−c
2w2teiw(x−v) = e−c
2w2t cos(wx−wv)+ ie−c2w2t sin(wx−wv) (183)
where the cos term is real and even and sin is imaginary and odd. Inserting this into the integral from −∞
to ∞, note that the sin term will integrate to zero. As the integral of the cos term from −∞ to ∞ is twice the
integral from 0 to ∞ means that the final simplification is:
T (x, t) =
1
pi
∫ ∞
−∞
f (v)
[∫ ∞
0
e−c
2w2t cos(wx−wv)dw
]
dv (184)
51
This can be solved by following Step 4 of the solution from the Fourier Integral, since it is identical to Eqn.
(159).
11.4.7 Conclusions
The Fourier Transform is an incredibly powerful tool in Engineering. The primary utility is in the dissection
of signals in spectral space, through examination of the power spectrum. It may also be used to solve PDEs,
however for the PDEs dealt with in this course it provides no direct advantage over the Fourier Integral.
11.5 Laplace Transform
11.5.1 Introduction
The Laplace Transform is a very powerful tool permitting the solution of problems which are not straight-
forward to solve using our previous approaches - such as problems with time dependent forcing or boundary
conditions. It is also a very useful method of analysing the response of a system governed by a PDE to a
given input. The Laplace Transform is defined as:
L { f}(s) = F(s) =
∫ ∞
0
exp−st f (t)dt, s = σ + Iω, (185)
where σ and ω are both real. Here we have considered f (t) as the Laplace Transform is most often
undertaken on temporally varying data, but it may equally be applied to f (x). To give a few examples:
f (t) = 1 for t ≥ 0, L { f}(s) =
∫ ∞
0
exp−st dt =−1
s
exp−st |∞0 =
1
s
(186)
f (t) = expat for t ≥ 0, L { f}(s) =
∫ ∞
0
expat exp−st dt =− 1
a− s exp
(a−s)t |∞0 =
1
s−a (187)
f (t) = expiωt , L { f}(s) = 1
s− iω =
s
s2+ω2︸ ︷︷ ︸
cosωt
+
iω
s2+ω2︸ ︷︷ ︸
I sinωt
=L {cosωt+ I sinωt} (188)
The final line shows a useful means of deriving the Euler formula. The Laplace Transform satisfies the
following useful identities:
1. Linearity: L {a f +bg}= aL ( f )+bL (g)
2. Derivatives: L { f ′}= sL { f}− f (0),L { f ′′}= s2L { f}− s f (0)− f ′(0)
3. First shifting theorem: L {expat f (t)}= F(s−a) =L { f (t)}(s−a)
4. Integrals: L {∫ τ0 f (τ)dτ}= 1s F(s)
5. Scaling in time: L { f (ct)}= 1c F( sc)
The derivative identity arises from integration by parts, L { f ′} = ∫ ∞0 exp−st f ′(t)dt = [exp−st f (t)]∞0 +
s
∫ ∞
0 exp
−st f (t)dt =− f (0)+ sL { f}.
The second shifting theorem is worth a little more detail as it is very commonly seen in problems with
a step change in forcing - imagine a forcing term which is zero until a specific time, then turned on, then
52
turned off again for example. We start by taking a function f (t) and shift by a distance a in the t dimension,
giving f (t−a). Next, we want that function to ‘turn on’ at t = a, thus we multiply by the unit step function
u(t−a), where
u(t−a) =
{
0, if t < a
1, t ≥ a . (189)
The Laplace transform of this isL { f (t−a)u(t−a)}= exp−as F(s).
There is a very useful table of Laplace Transforms on the datasheet. Thanks to the linearity of the transform,
the inverse Laplace Transform of complicated functions can be easily achieved, as will be seen in the
proceeding example.
11.5.2 A motivating example: application to ODEs
Here we’ll use an example of the solution of an ODE to both introduce the terminology associated with the
solution using the Laplace Transform, and how powerful the approach can be. The process for solving an
ODE using the Laplace Transform is as follows:
1. Transform the Equation and gain the Subsidiary equation
2. Define a Transfer function and gain an expression for F(s) by rearranging the Subsidiary equation
3. Apply initial and/or boundary conditions, then take an inverse transform/s to gain the solution in
physical space
In our example, we will solve y′′+ay′+by = r(t) where r(t) is a forcing term as a function of time, which
is the input to the system. The response to the input is y(t), which may be seen as the output. The problem
statement is closed by boundary conditions y(0) = K0 and y′(0) = K1. To give a quantitative solution, the
constants a and b must be provided, but here it is not needed. Following the steps:
1. Transform the Equation and gain the Subsidiary equation: Defining Y = L{y(t)} and R =
L{r(t)}, and using the derivative relations above we gain the Subsidiary Equation:[
s2Y − sy(0)− y′(0)]+a [sY − y(0)]+bY = R(s) (190)
Collecting terms in Y :
(s2+as+b)Y = (s+a)y(0)+ y′(0)+R (191)
2. Define a Transfer function and gain an expression for F(s):
Define the Transfer Function Q(s) = 1s2Y+as+b , gives:
Y (s) = ((s+a)y(0)+ y′(0))Q(s)+R(s)Q(s) (192)
For a moment assume a system initially at rest where y(0) = y′(0) = 0, then the subsidiary equation
becomes Y (s) = R(s)Q(s), thus the Transfer function can be interpreted as the ratio of the output of
the system (Y (s)) divided by the input to the system (R(s)), Q = Y/R.
53
3. Apply ICs/BCs and gain the inverse transforms: At this point we must chose values of a and
b and a functional form of r(t) and then reduce Y (s) = ((s+ a)y(0) + y′(0))Q(s) + R(s)Q(s) via
partial fractions so that L −1{Y} may be found. As an example, given r(t) = t, a = 0, b = −1,
y(0) = y′(0) = 1, then we are solving y′′− y = t. ThusL {r(t)}= 1s2 , Q = 1s2−1 and:
Y (s) = (s+1)Q(s)+
1
s2
Q(s) =
s+1
s2−1 +
1
s2(s2− s) (193)
, =
1
s−1 +
1
s2−1 −
1
s2
, (194)
where the final line is gained using partial fractions. The solution temporal space is now gained by
comparing the individual terms with the transform tables (or by explicitly transforming each individ-
ual term). This is a useful result which follows from the linearity of the transform. Thus:
y(t) =L −1{Y}=L −1
(
1
s−1 +
1
s2−1 −
1
s2
)
= expt +sinh(t)− t (195)
The Laplace transform allowed us to reach quite a complex solution with several individually straightfor-
ward steps, following the model of (i) transforming the problem, (ii) solving in transformed space, then (iii)
inverting the transform to gain a solution in physical space.
11.6 Application to PDEs
The approach to solving PDEs is only made slightly more complex by introducing the need to solve an
ODE in the transformed space, in much the same manner as is done in the Fourier Transform solution to
PDEs outlined in Section 11.4.6. The steps are:
1. Take the Laplace Transform of the PDE with respect to one of the two variables (normally time)
2. Solve the resultant ODE to obtain the transformed solution, typically F(s,x)
3. Take the inverse transform to obtain the solution in physical space
Let’s examine this approach with the help of an example application to the wave equation ytt−c2yxx = 0
- however, unlike previously where the wave equation was solved on a finite domain with fixed boundary
conditions (homogeneous or inhomogeneous), we will now look at a semi-infinite problem 0 ≤ x ≤ ∞
where the end of the cable is given a forced motion, i.e. x(0, t) = f (t). At the other end of the cable
limx→∞ y(x, t) = 0, and the cable is at rest initially y(x,0) = 0 and yt(x,0) = 0. How can we solve this, and
how can we interpret the solution?
So, we have the governing equations, the boundary condition, and we’ll choose a specific form of f (t):
f (t) =
{
sin t, if 0 < t < 2pi
0, otherwise
. (196)
Let’s follow the steps to solve using the Laplace Transform:
54
1. Take the Laplace Transform of the PDE: Considering each term individually, we take the Laplace
transform with respect to time to gain an ODE in Y (x,s):
L {ytt(x, t)}= s2L {y}−sy(x,0)− yt(x,0)︸ ︷︷ ︸
=0 from BCs
= s2L {y}= s2Y (197)
c2L {yxx(x, t)}= c2
∫ ∞
0
yxx exp−st dt = c2
∂ 2
∂x2
∫ ∞
0
yexp−st dt = c2
∂ 2L {y}
∂x2
= c2Yxx (198)
Thus the transformed PDE becomes the ODE:
s2Y = c2Yxx (199)
or:
Y =
c2
s2
Yxx (200)
2. Solve the resultant ODE: The solution to this ODE is:
Y (x,s) = A(s)expsx/c+B(s)exp−sx/c (201)
We now need to apply boundary conditions and initial conditions to solve for A(s) and B(s). Defining
F(s) =L { f (t)}=Y (0,s) and limx→∞ y(x, t) = 0 becomes limx→∞Y (x,s) = 0. Applying these to Eqn
(201) and noting in particular what happens to the exponentials as x→ ∞ gives:
Y (x,s) = F(s)exp−sx/c (202)
3. Take the inverse transform to obtain the solution in physical space:
Recall the second shifting theorem,L { f (t−a)u(t−a)}= exp−as F(s), or equivalently f (t−a)u(t−
a) =L −1{exp−as F(s)}. Given a = x/c, we thus have:
L −1{Y (x,s)}= y(x, t) = f
(
t− x
c
)
u
(
t− x
c
)
, (203)
where:
u
(
t− x
c
)
=
{
0, if t < xc
1, t ≥ xc
, and f
(
t− x
c
)
=
{
sin
(
t− xc
)
, if 0≤ t− xc ≤ 2pi
0, otherwise
. (204)
Which is our final solution, in quite a surprisingly compact and straightforward form. A quick note
- you may be wondering how this satisfies the boundary condition that yt(0,0) = 0? Taking the first
derivative of the solution with respect to time we gain:
yt(x, t) = ft
(
t− x
c
)
u
(
t− x
c
)
,
Thus at (x, t) = (0,0), yt(0,0) = 0 as the unit step function, or the Heavidside function as it is also
known, takes the value u(0) = 0.
55
The physical interpretation of our solution is equally important as gaining the solution itself. What
physical behaviour does Eqn. (203) represent? It is dramatically different from our previous solutions. The
solution consists of a travelling wave propagating with speed c along the cable, unchanged in amplitude
and shape. This wave is the one generated by the disturbance at the boundary condition. Note that f (t) =
sin t, but the wave varies in space as sin
(
t− xc
)
, such that the shape of the wave is reversed as plotted in
space compared to it’s plot in time. This is a common ‘gotcha’ in this solution of the wave equation, which
arises as the first part of the wave in time propagates down the cable first, becoming the head of the wave
in space. As the wave travels at speed c the head of the wave passes position x at time t = x/c, and a single
wavelength passes that same location in time x/c+2pi .
The travelling wave solution illustrates exactly how this equation behaves - waves travel in x as t increases.
In the case where there are fixed boundary conditions at either side of a finite cable, the wave equation gives
solution which are standing waves, but these are in reality manifestations of waves travelling along the wire
- simply that they sum to make a standing wave. This is most apparent if you choose an asymmetric initial
condition for a cable with two zero Dirichlet boundary conditions at x = 0 and x = L.
The real beauty of the solution in Eqn. (203) is that up until the very final step we did not need to consider
the transform of the forcing term f (t), so it is a very general solution for any forcing f (t) for a semi-infinite
wave equation.
Our second example solves the heat equation on a semi-infinite domain with time dependent heating at
x = 0. The template problem is to calculate the time dependent heating of permafrost due to solar heating,
where x = 0 is the surface of the permafrost and x is positive downwards. The time dependent heating is
represented by a boundary condition T (o, t) = f (t), with the second boundary condition as T (x→∞, t) = 0.
We expect the heating to penetrate down through the permafrost as time increases.
The governing equation is Tt = c2Txx on a semi-infinite domain (0 ≤ x ≤ ∞) with initial temperature
T (x,0) = 0. Let’s follow the solution steps again:
1. Take the Laplace Transform of the PDE: Considering each term individually, we take the Laplace
transform with respect to time to gain an ODE inL {T (x, t)}= T (x,s):
L {Tt(x, t)}= sL {T} −T (x,0)︸ ︷︷ ︸
=0 from BCs
= sT , (205)
as with the wave equation,
c2L {Txx(x, t)}= c2T xx (206)
Thus the transformed PDE becomes the ODE:
sT = c2T xx (207)
or:
T =
c2
s
Yxx (208)
2. Solve the resultant ODE: The solution to this ODE is:
Y (x,s) = A(s)exp
√ s
c2
x
+B(s)exp−
√ s
c2
x (209)
We now need to apply boundary conditions and initial conditions to solve for A(s) and B(s). Defining
F(s) =L {T (t)} = T (0,s) and limx→∞T (x, t) = 0 becomes limx→∞T (x,s) = 0. Applying these to
Eqn (209) and noting in particular what happens to the exponentials as x→ ∞ gives:
56
T (x,s) = F(s)exp−
√ s
c2
x (210)
3. Take the inverse transform to obtain the solution in physical space:
Unlike the previous example of the wave equation, we must define the forcing term prior to taking
the inverse transform. Let’s choose an impulse f (t) = 1 for t ≥ 0, which when transformed yields
F(s) = 1s . Thus:
T (x,s) =
1
s
exp−
√ s
c2
x (211)
At this point we need to observe that:
L −1
{
exp−a
√
s
s
}
=
[
1− erf
(
a
2
√
t
)]
(212)
Setting a =−x/c gives:
T (x, t) =
[
1− erf
( −x
2c
√
t
)]
(213)
Which reassuringly is an error function solution as we may have expected from our experience with
the Fourier Integral solution to the heat equation.
11.7 Conclusions
In this section we have covered three incredibly important methods of analysis, the Fourier Integral, the
Fourier Transform and the Laplace Transform. These are very important in the solution of PDEs on infinite
and semi-infinite domains, for signal processing and computing power spectra. Laplace transforms are
particularly useful for the analysis of systems which are forced, either through the boundary condition or in
the domain itself (no covered explicitly in this course).
57
12 Finite Element Analysis
12.1 Introduction
Finite Element Analysis (FEA) is a numerical approach which has become the default method for solving a
plethora of engineering problems, but has found it’s most notable application in the computational analysis
of structures. Indeed, undertaking an FEA of a structure is synonymous with computing the deflections,
stresses and strains on a computational model of an engineering device. It has wider applications too, from
thermal analysis (heat equation) through to magnetic flux, seepage of groundwater flows and most recently
encroaching on fluid dynamics. Here we will first of all introduce the method as it would broadly apply to
a wider class of equation, then we will focus on the Poisson equation, which is the most basic model of a
structure which we have available and is perfect to help understand the complexities of the Finite Element
approach. We will, however, focus on ODEs since the extension to PDEs is a little more complex than
would be possible in the time allotted in this course.
So what is a finite element? It is a method of piecewise approximation in which the approximating function
is formed by connecting simple functions, each defined over a small region (element). A finite element is a
region in space in which a function is interpolated from nodal values on the boundary of the region in such
a way that inter-element continuity of the function tends to be maintained in the assemblage. Let’s compare
this to finite difference which may help clarify how FEA differs - finite difference assumes that we have a
single value of the quantity of interest (say Temperature) at a single point. Finite Element would instead
assume a distribution of temperature over an element - where an element may be envisaged to be a region
between point (or node) 1 and 2. Thus a finite element method represents not just the values of a property
at a point (as in Finite Differences) but also how that quantity varies over the element between the nodes
(points). This can be extended to two and three dimensions in a relatively straightforward manner.
However, that is focussing on a finite element - but FEA is more than a definition of an element, it is
encompasses an approach to solving a governing equation which has as it’s goal the minimisation of errors
at specific locations (nodes) within the domain. This is very different to the finite difference approach, and
appears almost like a magic trick the first time you look at it. The FEA approach is divided into two main
parts, the first is the establishment of the weak form of the governing equations - weak form in that it is
posed in integral rather than differential form. The second part is the solution of the weak form through
selection of the distribution of the quantity of interest within an element (basis and test function), the mesh
itself and then solving.
12.2 Formalisation for General ODEs
The solution of differential equations using the finite element method (FEM) is done by following the six
basic steps below:
1. Establish the strong formulation: this is the governing PDE
2. Establish the weak formulation by multiplying the error by an arbitrary field and integrating over the
whole domain
3. Select shape and weight functions and discretisation
4. Compute the stiffness matrix and load vector
5. Solve global system of equations for nodal values of the primary variables (displacements/temperature)
6. Compute temperature/stresses/strains etc. within the element using nodal values and the chosen shape
functions
58
In this process, steps 1 and 2 only have to be done once for each governing equation. Thus steps 3-6 are
the steps which are applied each time a new problem is solved. In a commercial software, steps 4 and 5 are
usually hidden from the user.
Complete Derivation to the Galerkin Approach to the Weak Form for the 1D Poisson Equa-
tion
Let’s got through each of the above steps in detail for the 1D Poisson equation, which is a model equation
which may represent both Steady Heat Equation with source/sink terms, and a stress analysis of a one
dimensional beam under loading. The equation is:
kTxx =−q(x) (214)
The motivating example we use here is thermal conduction in a rod of length L, thus the domain is [0,L],
T (0) = gD = 273K (Dirichlet) and the end of the rod is insulated, i.e. Tx(L) = 0 (Neumann). The rod is
also exposed to a constant forcing q.
To solve this problem we need to be a little more general than what you have seen so far in the lectures and
in doing so things become a little messy, so just hold on tight and try to follow all the steps.
Step 1: Establish the strong formulation
The strong (differential) form of the governing equation is:
kTxx =−q(x) (215)
Step 2: Establish the weak formulation
Starting with the Poisson equation given for steady-state temperature we will invoke the use of the Residual
theorem simply to introduce the variables necessary to solve the complete problem. Let the differential
operator L(T ) be such that;
L(T ) =
d2T
dx2
(216)
then the strong form is:
kL(T ) =−q(x) (217)
Now we introduce the approximate variables, such that if we solve L(T δ ) we have a Residual R given by;
R(T δ ) = kL(T δ )+q(x) (218)
If the approximate solution, T δ = T then we have the exact solution and the residual is zero. Thus the resid-
ual is the error in the discretised solution, and may be envisioned as the error introduced by the numerical
approximations as a function of space. For a typical solution using e.g. Taylor Series, this error would be
the truncation terms in the Taylor Series.
Now, the Weak Formulation is constructed by the method of weighted residuals. This approach introduces
the use of an arbitrary test function v(x) which multiplies the residual R(T δ ),
59
v(x)R(T δ ) = v(x)kL(T δ )+ v(x)q(x). (219)
The designation of this test function as ‘arbitrary’ is important, as we shall see. Now with an ideal discreti-
sation the error would be zero, and so we set
v(x)kL(T δ )+ v(x)q(x) = 0 = v(x)kL(T )+ v(x)q(x). (220)
We would like to set the global error to zero, i.e. error being zero throughout the whole domain. A global
error measure is the integral of the residual over the whole domain, however, there is a difficulty that if you
simply integrate the residual over the whole domain it is possible that negative errors may cancel positive
errors and so give a false measure of the true error. So, instead of integrating the Residual (error) over
the whole domain, we instead integrate the product of the Residual and the arbitrary test function, i.e.∫ L
0 R(x)v(x)dx. The only way that integral can be zero for any arbitrary v(x) is if R(x) is zero.
Following this rationale, the next step is to integrate Eqn. (220) over the whole domain, where we explicitly
write L(T ) =
d2T (x)
dx2
which here represents the exact derivative, but later on we will replace this with the
discretised form using the finite element method. This leads to the weak form of the governing equation:
∫ L
0
v(x)k
d2T (x)
dx2
dx+
∫ L
0
v(x)q(x)dx = 0, (221)
where k may be moved out of the integration if it is a constant in space. Using integration by parts, we can
expand out the first integral on the RHS to make a second, more useful weak form;
∫ L
0
v(x)q(x)dx+ k
[
v(x)
dT (x)
dx
]L
0
− k
∫ L
0
dv(x)
dx
dT (x)
dx
dx = 0. (222)
This form has some important advantages, most notably that the application of a heat flux (kTx) at the
boundaries is explicitly and algebraically incorporated into the equation. Note that in this example kTx is
the heat flux, however as applied to the Poisson equation for displacement u, the equivalent term is EAux
which is the force applied at that point.This natural inclusion of Neumann boundary conditions is a strength
of the FEM. A final point - we have not specified a functional form for v(x) - this is deliberate as in the FEM
a functional form does not need to be specified, we actually solve the system in such a way that the integral
of v(x)R(x) is zero for any v(x). As you will see shortly, this is done by putting the resultant approximated
equation into the form vδFδ (T δ ,qδ ) = 0 where Fδ is some function of the approximated distribution of T
(T δ ) and q (qδ ), vδ is the approximation to v(x), and then setting Fδ (T δ ,qδ ) = 0.
The next step covers the method of approximation of the derivatives and test function used in Eqn. (222).
Step 3: Select shape and weight functions and discretisation
Since we are interested in solving the weak form of the governing equation we need to approximate each
term continuously over the whole domain. A systematic approach to this is to represent the spatial variation
of each quantity (here T (x) and v(x)) using basis functions.
What are basis functions? Basis functions are simple building blocks which may be summed to form a
more complex function. An example of a basis function which we have already come across is the Fourier
Series, where the set of basis functions are cos(npix/L) and sin(npix/L) where n = 0,1, .....
60
1.0
0 21
M
ag
ni
tu
de
Node number
Element 1 Element 2
N1 N2
N0
0 21
Node number
Element 1 Element 2
T^
T^
T^
1
2
T
0
Figure 5: Linear basis functions Ni employed in a finite element analysis solved using two elements (left),
and an example distribution plotting T δ = ∑nei=0 TˆiNi (right)
We will use linear basis function in this Unit of Study. These basis functions are plotted in Fig. 5 for
a case where the domain is discretised into two finite elements, using three nodes, which results in three
basis function N0, N1, N2 where each basis function is centred on node i. They are defined by the following
equations:
N0 =

le− x
le
if 0≤ x≤ le
0 if x≥ le
(223)
N1 =

x/le if 0≤ x≤ le
(2le− x)/le if le ≤ x≤ 2le
0 if otherwise
(224)
and,
N2 =

x/le−1 if le ≤ x≤ 2le
(3le− x)/le if 2le ≤ x≤ 3le
0 if otherwise
(225)
where le is the element length. In the example shown in Figure 5, le = L/2 and the basis functions simplify
to be:
N0 =
{
1− 2x
L
if 0≤ x≤ L/2
0 if x≥ L/2
(226)
N1 =

2x/L if 0≤ x≤ L/2
2
(
1− x
L
)
if L/2≤ x≤ L
0 if otherwise
(227)
and,
N2 =
{ 2x
L
−1 if L/2≤ x≤ L
(0 if otherwise
(228)
Note that at each location in x, the ∑nei=0 Ni(x) = 1, and that each individual basis function Ni is only non-
zero in the two neighbouring elements i and i− 1. They are equal to 1 at the corresponding node, and 0
for every subsequent node. Consider now that we discretise the domain into ne elements of equal length
le. The basis function Ni(x) are centred at node i, and nodal values are represented by ˆ(.) (termed a basis
vector) where (.) is a quantity of interest. We can then approximate the variation of T (x) and q(x) over a
given range as T δ and qδ defined as:
61
T δ =
ne
∑
i=0
TˆiNi (229)
qδ =
ne
∑
i=0
qˆiNi (230)
When linear basis functions are chosen, this may be visualised as a piecewise linear variation of each
function over the domain, and is plotted in Figure 5(right). At a specific node i, the only basis function
which is non-zero is Ni. Thus T (xi) = Tˆi, and values in between are interpolations of the two neighbouring
nodal values, where the basis functions are the interpolation mechanism.
The final use of the basis functions is in the representation of the test function v(x). We will adopt the
classical Galerkin approach, where the test function is also approximated by the same basis functions as T
and q, giving:
vδ =
ne
∑
i=0
vˆiNi (231)
Representing functions using basis functions means that we now have a continuous representation which
can be integrated at all x locations and differentiated everywhere except at the nodes.
In the previous step we showed how Neumann boundary conditions are very simply represented by the
single term k [v(x)Tx]
L
0 . There is a final modification to how we approximate the function to enable the
easy integration of Dirichlet boundary conditions. We represent the unknown as a the sum of a Dirichlet
component T D, and internal nodes or non-Dirichlet boundary nodes forming a homogeneous component
T H . This gives the following form in our current example where the left hand boundary is Dirichlet, hence
Tˆ (0) = T (0) = gD:
T δ = T H +T D
where T H =
ne
∑
i=1
TˆiNi
and T D = gDN0
Furthermore, since the test function is arbitrary, and since errors are zero at a Dirichlet boundary, we can
set v = 0 at the Dirichlet boundary without any loss of generality of the method. Thus the test function
becomes;
vδ =
ne
∑
i=1
vˆiNi,
qδ =
ne
∑
i=0
qˆiNi,
noting that the approximation to the source/sink term qδ has remained the same and includes qˆ0N0.
Step 4: Compute the stiffness matrix and load vector
The problem is now fully defined, we will solve Eqn. (222) using the distributions defined in Eqn (12.2)
and (12.2) for the specific boundary conditions in this specific problem. Firstly, noting that we have split
the temperature field into Dirichlet T D(x) and homogeneous T H(x) components, we make that substitution
in the integrals in Eqn. (222) to give:
62
∫ L
0
v(x)q(x)dx+ k
[
v(x)
dT (x)
dx
]L
0
− k
∫ L
0
dv(x)
dx
dT H(x)
dx
dx− k
∫ L
0
dv(x)
dx
dT D(x)
dx
dx = 0. (232)
The easiest approach to solving Eqn. (232) is to consider each individual term one at a time. First we note
that the test functions are defined to be zero on the Dirichlet boundary, i.e. v(0) = 0 and that we have a
Neumann boundary condition Tx(L) = 0, thus;
k
[
v(x)
dT (x)
dx
]L
0
= kv(L)
dT (L)
dx
− kv(0)dT (0)
dx
(233)
= k · v(L) ·0− k ·0 · dT (0)
dx
(234)
= 0 (235)
Thus this term can be removed from Eqn. (232). Next, substitute the approximations vδ , qδ and T δ into the
remaining terms in Eqn. (232);
k
∫ L
0
dvδ
dx
dT H
dx
dx =
∫ L
0
vδqδdx− k
∫ L
0
dvδ
dx
dT D
dx
dx (236)
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx =
∫ L
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx− k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
gD
dN0
dx
dx (237)
In the above equation, the only dependence on x arises in the basis functions, which appears as the original
function and the first derivative. Also note that we dropped the (.)H notation on the Tˆi for brevity. As
detailed previously, the first shape function has the form;
N0 =

le− x
le
if 0≤ x≤ le
0 if x≥ le
(238)
Now the second node has to take into consideration approaches from either direction, i.e.;
N1 =

x/le if 0≤ x≤ le
(2le− x)/le if le ≤ x≤ 2le
0 if otherwise
(239)
and,
N2 =

x/le−1 if le ≤ x≤ 2le
(3le− x)/le if 2le ≤ x≤ 3le
0 if otherwise
(240)
This pattern will continue on for every shape function until the end of the rod (i.e. ne · le = L). Noticing that
the first time on the LHS and the last term on the RHS have derivatives of these shape functions, it is useful
to compute these terms as well now;
dN1
dx
=

1/le if 0≤ x≤ le
−1/le if le ≤ x≤ 2le
0 if otherwise
(241)
dN2
dx
=

1/le if le ≤ x≤ 2le
−1/le if 2le ≤ x≤ 3le
0 if otherwise
(242)
(243)
63
This is where things get a bit nicer. Let us consider that the integral over the entire domain can be subdivided
by the element discretisation, i.e.;∫ L
0
f (x)dx =
∫ le
0
f (x)dx+
∫ 2le
le
f (x)dx+ ...+
∫ L
(ne−1)le
f (x)dx (244)
We may use this relationship to compute each integral term of equation (236) in a block format which can
then be repeated for simplicity. Take the LHS of equation (236);
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx = k
∫ le
0
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx+ (245)
k
∫ 2le
le
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx+ ...
The first term is simple to solve since only half of the second shape function exists in this domain (refer to
Eqn. (242))). So this becomes;
k
∫ le
0
vˆ1
dN1
dx
Tˆ1
dN1
dx
dx = k
∫ le
0
vˆ1(1/le)Tˆ1(1/le)dx (246)
= k
[
x
l2e
vˆ1Tˆ1
]le
0
(247)
=
k
le
vˆ1Tˆ1 (248)
Now we go to the next term in the integral (note: you must here consider the bounds of the integral and
hence which shape functions are non-zero);
k
∫ 2le
le
(
vˆ1
dN1
dx
+ vˆ2
dN2
dx
)(
Tˆ1
dN1
dx
+ Tˆ2
dN2
dx
)
dx =
k
∫ 2le
le
(
−vˆ1 1le + vˆ2
1
le
)(
−Tˆ1 1le + Tˆ2
1
le
)
dx (249)
=
k
l2e
∫ 2le
le
(
vˆ1Tˆ1− vˆ1Tˆ2− vˆ2Tˆ1+ vˆ2Tˆ2
)
dx (250)
=
k
l2e
(
vˆ1Tˆ1− vˆ1Tˆ2− vˆ2Tˆ1+ vˆ2Tˆ2
)
[x]2lele (251)
=
k
le
(
vˆ1Tˆ1− vˆ1Tˆ2− vˆ2Tˆ1+ vˆ2Tˆ2
)
(252)
Now we can express this result as a matrix product since the basis vectors vˆi and Tˆi are independent, giving;
k
∫ 2le
le
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx =
k
le
[
vˆ1 vˆ2
][ 1 −1
−1 1
][
Tˆ1
Tˆ2
]
(253)
You should notice here that we have solved the stiffness matrix for a single element (element 2 to be exact),
but the formulation is exactly the same for each subsequent element due to the nature of the shape functions
(the pattern repeats identically for each element of identical length). For clarity, the next integration element
will have the result;
k
∫ 3le
2le
(
vˆ2
dN2
dx
+ vˆ3
dN3
dx
)(
Tˆ2
dN2
dx
+ Tˆ3
dN3
dx
)
dx =
k
le
[
vˆ2 vˆ3
][ 1 −1
−1 1
][
Tˆ2
Tˆ3
]
(254)
64
and;
k
∫ 4le
3le
(
vˆ3
dN3
dx
+ vˆ4
dN4
dx
)(
Tˆ3
dN3
dx
+ Tˆ4
dN4
dx
)
dx =
k
le
[
vˆ3 vˆ4
][ 1 −1
−1 1
][
Tˆ3
Tˆ4
]
(255)
and so on until the final element (Note; if you do not believe that this is the case you are strongly encouraged
to continue solving these elements to verify the results above). The first element only differed as it has a
Dirichlet boundary condition thus v(0)= 0 and it only includes the homogeneous component of the solution.
Now lets assemble the total solution based on these results;
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx =
k
le
vˆ1Tˆ1+
k
le
[
vˆ1 vˆ2
][ 1 −1
−1 1
][
Tˆ1
Tˆ2
]
+
k
le
[
vˆ2 vˆ3
][ 1 −1
−1 1
][
Tˆ2
Tˆ3
]
+ ...
k
le
[
vˆne−1 vˆne
][ 1 −1
−1 1
][
Tˆne−1
Tˆne
]
(256)
=
k
le
[
vˆ1 vˆ2 vˆ3 . . . vˆne−1 vˆne
]

2 −1 0 . . . 0 0
−1 2 −1 . . . 0 0
0 −1 2 . . . 0 0
...
. . . . . . . . .
...
...
0 0 . . . −1 2 −1
0 0 . . . 0 −1 1


Tˆ1
Tˆ2
Tˆ3
...
Tˆne−1
Tˆne

(257)
Now, we have completed one of the integrals from equation (236). Let us now look at the solution for the
first integral on the RHS;
∫ L
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx (258)
We know from the problem description that q(x) is a constant q. Thus if you were undertaking this deriva-
tion for a real problem or in an exam at this point you would recognise that ∑nej=1 qˆ jN j = q, and make this
substitution into the above equation to simplify the calculation. Here, we will evaluate the term as though
the body force is a function of x, so that you have an example of how to do this in the case of non-constant
q(x). We are going to break this function up in the same way we did with the LHS integral, therefore;∫ L
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx =
∫ le
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx+
∫ 2le
le
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx+ ... (259)
Similarly, we approach the first integral term as it is going to look slightly different to the other integrals;∫ le
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx =
∫ le
0
vˆ1N1 (qˆ0N0+ qˆ1N1)dx (260)
Note here that there are two terms for qδ in this range as we consider the both nodes, where as before we
65
know that vˆi and Tˆi must both start at the second element (i.e. i = 1). The integral may now be evaluated,∫ le
0
vˆ1N1 (qˆ0N0+ qˆ1N1)dx =
∫ le
0
vˆ1
x
le
(
qˆ0
(
1− x
le
)
+ qˆ1
x
le
)
dx (261)
=
∫ le
0
vˆ1
[
qˆ0
(
x
le
− x
2
l2e
)
+ qˆ1
(
x2
l2e
)]
dx (262)
= vˆ1
[
qˆ0
(
x2
2le
− x
3
3l2e
)
+ qˆ1
(
x3
3l2e
)]le
0
(263)
= vˆ1
(
le
6
qˆ0+
le
3
qˆ1
)
(264)
At this point we can say that if there is a uniform body force on the element, then qˆ0 = qˆ1 = q thus;∫ le
0
vˆ1N1 (qˆ0N0+ qˆ1N1)dx = vˆ1le
q
2
(265)
On to the next integral term (this is the most time consuming component of a derivation for varying q(x));∫ 2le
le
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx =
∫ 2le
le
(vˆ1N1+ vˆ2N2)(qˆ1N1+ qˆ2N2)dx (266)
=
∫ 2le
le
(
vˆ1qˆ1N21 + vˆ1qˆ2N1N2+ vˆ2qˆ1N1N2+ vˆ2qˆ2N
2
2
)
dx (267)
= vˆ1qˆ1
∫ 2le
le
N21 dx+ vˆ1qˆ2
∫ 2le
le
N1N2dx+ ...
vˆ2qˆ1
∫ 2le
le
N1N2dx+ vˆ2qˆ2
∫ 2le
le
N22 dx (268)
Now we can compute the shape functions and their squares and integrate directly since the vˆi and qi terms
are not dependent on the variable x;
N21 =
(
2− x
le
)2
= 4+
x2
l2e
− 4x
le
(269)
N22 =
(
x
le
−1
)2
= 1+
x2
l2e
− 2x
le
(270)
N1N2 =
(
2− x
le
)(
x
le
−1
)
=
3x
le
− x
2
l2e
−2 (271)
And the integrals of these functions are;∫ 2le
le
N21 dx =
∫ 2le
le
(
4+
x2
l2e
− 4x
le
)
dx =
[
4x+
x3
3l2e
− 2x
2
le
]2le
le
=
le
3
(272)
∫ 2le
le
N22 dx =
∫ 2le
le
(
1+
x2
l2e
− 2x
le
)
dx =
[
x+
x3
3l2e
− x
2
le
]2le
le
=
le
3
(273)
∫ 2le
le
N1N2dx =
∫ 2le
le
(
3x
le
− x
2
l2e
−2
)
dx =
[
3x2
2le
− x
3
3l2e
−2x
]2le
le
=
le
6
(274)
(275)
Combining these results with equation (268) gives us;∫ 2le
le
(vˆ1N1+ vˆ2N2)(qˆ1N1+ qˆ2N2)dx =
le
3
vˆ1qˆ1+
le
6
(vˆ1qˆ2+ vˆ2qˆ1)+
le
3
vˆ2qˆ2 (276)
66
Which in vector format appears as;
le
3
vˆ1qˆ1+
le
6
(vˆ1qˆ2+ vˆ2qˆ1)+
le
3
vˆ2qˆ2 =
[
vˆ1 vˆ2
]
le
[
qˆ1/3+ qˆ2/6
qˆ1/6+ qˆ2/3
]
(277)
Note for this problem, qˆ1 = qˆ2 = q then we get;∫ 2le
le
(vˆ1N1+ vˆ2N2)(qˆ1N1+ qˆ2N2)dx =
[
vˆ1 vˆ2
]
le
[
1/2
1/2
]
q (278)
As with the stiffness matrix terms, this pattern for the forces repeat themselves for each successive element,
and hence we can construct our force vector by simply adding terms;∫ L
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
qˆ jN jdx = vˆ1le
q
2
+
[
vˆ1 vˆ2
]
le
[
1/2
1/2
]
q (279)
+
[
vˆ2 vˆ3
]
le
[
1/2
1/2
]
q+ ...+
[
vˆne−1 vˆne
]
le
[
1/2
1/2
]
q (280)
= qle
[
vˆ1 vˆ2 vˆ3 . . . vˆne−1 vˆne
]

1
1
1
...
1
1/2

(281)
Now we have the complete (and required) forms for two of the three integrals! Fortunately, the last integral
contains our Dirichlet boundary and is nice and simple to compute;
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
gD
dN0
dx
dx (282)
You should notice immediately that here the N0 shape function is explicitly defined, which means we are
only interested in solving the integral in the region 0≤ x ≤ le since all other values for x are 0. So lets get
to it;
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
gD
dN0
dx
dx = k
∫ le
0
vˆ1
dN1
dx
gD
dN0
dx
dx (283)
= k
∫ le
0
vˆ1
(
1
le
)
gD
(−1
le
)
dx (284)
= vˆ1
kgD
l2e
[−x]le0 =−vˆ1
kgD
le
(285)
Now lets assemble the full expression from each of the integral parts we have formulated (Equations (257),
(281) and (285));
67
k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
ne
∑
j=1
Tˆj
dN j
dx
dx =
∫ L
0
ne
∑
i=1
vˆiNi
ne
∑
j=1
q jN jdx− k
∫ L
0
ne
∑
i=1
vˆi
dNi
dx
gD
dN0
dx
dx (286)
⇒ k
le
[vˆ]

2 −1 0 . . . 0 0
−1 2 −1 . . . 0 0
0 −1 2 . . . 0 0
...
. . . . . . . . .
...
...
0 0 . . . −1 2 −1
0 0 . . . 0 −1 1


Tˆ1
Tˆ2
Tˆ3
...
Tˆne−1
Tˆne

= qle [vˆ]

1
1
1
...
1
1/2

+ vˆ1
k
le
273 (287)
k
le
[vˆ]

2 −1 0 . . . 0 0
−1 2 −1 . . . 0 0
0 −1 2 . . . 0 0
...
. . . . . . . . .
...
...
0 0 . . . −1 2 −1
0 0 . . . 0 −1 1


Tˆ1
Tˆ2
Tˆ3
...
Tˆne−1
Tˆne

= le [vˆ]

q+273k/le
q
q
...
q
q/2

(288)
(289)
The final step here is merely one for mathematical consistency, but it is necessary to get rid of the test
functions, vˆi. If we move everything to the LHS we see;
[vˆ]

k
le

2 −1 0 . . . 0 0
−1 2 −1 . . . 0 0
0 −1 2 . . . 0 0
...
. . . . . . . . .
...
...
0 0 . . . −1 2 −1
0 0 . . . 0 −1 1


Tˆ1
Tˆ2
Tˆ3
...
Tˆne−1
Tˆne

− le

q+273k/le
q
q
...
q
q/2


= 0 (290)
Now this equation must be true for all test functions vˆ to be a valid solutions, and the only way that is
possible is if the matrix within the curly brackets equals zero. The final step in the process is to solve the
problem:
k
le

2 −1 0 . . . 0 0
−1 2 −1 . . . 0 0
0 −1 2 . . . 0 0
...
. . . . . . . . .
...
...
0 0 . . . −1 2 −1
0 0 . . . 0 −1 1


Tˆ1
Tˆ2
Tˆ3
...
Tˆne−1
Tˆne

= le

q+273k/le
q
q
...
q
q/2

. (291)
Step 5: Solve global system of equations
This requires most of the computational effort - the solution of the relatively straightforward problem (Eqn.
(291) in the generic formK Tˆ =R, where the matrixK is the global stiffness matrix andR is the global
load vector. This can be achieved using specialised routines, as has been discussed before in the sections
on implicit time stepping, or the discretisation of the Laplace & Poisson equations. Solving this will give
you the vector of unknown nodal values of Tˆi.
Step 6: Compute temperature/stresses/strains etc. within the element
Given that vector, the approximate distribution of T (x) which is T δ can be recovered by noting that:
68
T δ (x) = T H +T D (292)
where T H =
ne
∑
i=1
TˆiNi (293)
and T D = gDN0 (294)
which would give you the solution for T at any location along the bar. Any quantities which require
derivatives may be recovered easily.
General Aside: This whole process, although lengthy, has been straightforward, and hopefully you can
see how easy it would be to extend - particularly to an arbitrary number of elements. This procedure is
given here in full detail so that you, the student, can understand how the full matrix solution for an arbitrary
number of elements. It is important to understand the key steps, as you may encounter problems where
for example, the Neumann condition at the endpoint is non-zero, or your Dirichlet boundary is zero. You
should use this method as a guide to understand how to form the general element force and stiffness matri-
ces and how they combine to form the global matrix.
12.3 Does this mean that the Finite Element Method has Zero Error?
Have we somehow derived a solution with zero error? This seems like an impossibility - there is always
a loss of information between the real distribution and the approximated distribution, yet the whole basis
of the finite element method is setting the error to zero. It is very important to appreciate the following two
reasons why the error is not zero in most practical problems:
1. The foundation of the derivation was that
∫ L
0 v(x)R(T
δ )dx=
∫ L
0 v(x)
(
L(T δ )+q(x)
)
dx= 0. However,
this assumes that q(x) is approximated exactly, and that v(x) is continuous, whereas in application
both are approximated using the same linear basis functions as used for T (x). Thus, the residual
which is actually forced to zero is
∫ L
0 v
δ (x)
(
L(T δ )+qδ (x)
)
dx = 0. It turns out that the discretised
form only matches the actual form when the choice of basis functions is one order less than the exact
solution, for example a linear basis function can represent a quadratic solution exactly at the nodes.
2. When post processing the results, the basis functions are used. So, even if the nodal values computed
from the finite element method are exactly correct, unless the distribution of, e.g. T (x) is linear in x,
then the use of linear basis functions to determine values between the nodes will give an erroneous
answer. This is very clear when computing derived quantities such as stress, σ = Eux, where the
predicted stress is constant within each element when using linear basis functions - this is obviously
not exact for any solution which is quadratic or higher order.
In reality, both of these errors may be reduced by using higher order basis functions, such as quadratic,
cubic or quartic. In reality, of course, for structural analysis there would be additional uncertainties in point
load magnitudes and directions, material properties, joints etc which must be modelled.
12.4 Specialisation of FEA to Structural Analysis
12.4.1 Introduction
Here we are going to dive straight into example applications, this really is the best way to get a feel for the
process of approximating the unknown functions and forces with a chosen basis function, evaluating the
individual terms in the weak formulation, simplifying and solving for the unknowns.
69
12.4.2 Predicting Stresses in a Simple Beam under Compression
An engineering firm have asked you to compute the stress within a beam under compression. This can
be modelled as a one dimensional beam in x with constant cross-sectional area A, modulus of elasticity
E and total length L. The beam is fixed at x = 0, and at x = L it is compressed by a fixed displacement
of u(L) = −L/10. There are no body forces (q = 0) or further point loads applied. The displacement
distribution u(x) in the beam is given by Poisson’s equation ∂∂x (EAux) + q(x) = 0. Note that the weak
formulation of this equation is,
∫ L
0
vxEAuHx dx =
∫ L
0
vqdx+[vEAux]L0−
∫ L
0
EAvxuDx dx. (295)
We will now simplify the weak form of the governing equation and calculate the stiffness and load vector,
for a discretisation of the beam using two elements. The final step will be to calculate the displacement
at the middle node using linear shape functions. Thus the approximation of uδ follows approximately that
shown in Figure 5.
First of all simplify the governing equations to remove any terms which are zero. The weak form of the
governing equation is:
∫ L
0
vxEAuxdx =
∫ L
0
vqdx+[vEAux]L0−
∫ L
0
EAvxuDx dx (296)
The Dirichlet Boundary conditions in this problems are at x = 0 and x = L, such that the test functions
v(0) = 0 and v(L) = 0. Hence [vEAux]L0 = 0. We are given q(x) = 0 so
∫ L
0 vqdx = 0, leaving us with:
∫ L
0
vxEAuHx dx =−
∫ L
0
EAvxuDx dx (297)
Next step is to define the solution u(x) = uH +uD, where uH = uˆ1N1 and uD = uˆ0N0+ uˆ2N2, where uˆi is the
nodal displacement at node i and Ni is the corresponding shape function, defined for the 3 node bar as:
N0 =
{
1− 2xL 0≤ x < L2
0 L2 ≤ x≤ L
,N1 =
{
2x
L 0≤ x < L2
2
(
1− xL
) L
2 ≤ x≤ L
N2 =
{
0 0≤ x < L2
2x
L −1 L2 ≤ x≤ L
We can also specify the derivatives of the shape functions (as will be required in evaluating the RHS term
in Eqn. (297)):
N0,x =
{
− 2L 0≤ x < L2
0 L2 ≤ x≤ L
,N1,x =
{
2
L 0≤ x < L2
− 2L L2 ≤ x≤ L
N2,x =
{
0 0≤ x < L2
2
L
L
2 ≤ x≤ L
Now we can evaluate the first term in Eqn. (297):
70
∫ L
0
vxEAuHx dx =
∫ L/2
0
vxEAuHx dx+
∫ L
L/2
vxEAuHx dx (298)
= EA
∫ L/2
0
4
L2
uˆ1vˆ1dx+EA
∫ L
L/2
4
L2
uˆ1vˆ1dx (299)
= EA
4
L
uˆ1vˆ1 (300)
Where in the above it is noted that each integral is simply the integral of a constant.
Next for the right hand side of Eqn. (297). Noting that uD = uˆ0N0+ uˆ2N2:
−
∫ L
0
EAvxuDx dx = −
∫ L/2
0
EAvxuDx dx−
∫ L
L/2
EAvxuDx dx (301)
= −
∫ L/2
0
EA
−4
L2
vˆ1uˆ0dx−
∫ L
L/2
EA
−4
L2
vˆ1uˆ2dx (302)
= EA
2
L
vˆ1uˆ0+EA
2
L
vˆ1uˆ2 (303)
(304)
Equating the left hand and right hand side in Eqn. (297), and cancelling terms:
EA
4
L
uˆ1vˆ1 = EA
2
L
vˆ1uˆ0+EA
2
L
vˆ1uˆ2 (305)
2uˆ1 = uˆ0+ uˆ2 (306)
uˆ1 =
1
2
(uˆ0+ uˆ2) (307)
(308)
As the boundary values are known this gives the solution for only unknown uˆ1.
12.4.3 Predicting Stresses in a Prosthetic Leg
A prosthetic leg may be modelled as a vertical beam, attached to the ground (x = 0) with a load P applied
at the top (x = L). This can be modelled as a one dimensional beam in x with constant cross-sectional area
A, modulus of elasticity E and total length L. Body forces may be neglected (q = 0). The displacement
distribution u(x) in the beam is given by Poisson’s equation ∂∂x (EAux)+q(x) = 0.
Establish and analytical solution. We are given that E and A are constant and q(x) = 0. We integrate
twice to recover the displacement distribution:
∂
∂x
(EAux) = 0 → EAux = B → EAu = Bx+C
Noting we have the boundary conditions u(0) = 0 and EAux(L) =−P, we find that
u(0) = 0 → 0 = 0+C → C = 0
71
Figure 6: Discretisation of the prosthetic leg using two finite elements
, and
EAux(L) =−P → B =−P
. The analytical solution is:
u(x) =
−Px
EA
Determine the FEA approximation We will now simplify the weak form of the governing equation and
calculate the stiffness and load vector, under the assumption that we discretise the prosthetic using two
elements and using linear basis functions. The leg is discretised as in Figure 6.
The weak form of the governing equation is:
∫ L
0
vxEAuHx dx =
∫ L
0
vqdx+[vEAux]L0−
∫ L
0
EAvxuDx dx (309)
The Dirichlet Boundary condition in this problem is zero at x = 0 so
∫ L
0 EAvxu
D
x dx = 0. We are given
q(x) = 0 so
∫ L
0 vqdx = 0, leaving us with:
EA
∫ L
0
vxuxdx = EA[vux]L0 (310)
As area and Young’s modulus are constant. Now we set uδ = uH +uD, where uD = uˆ0N0 and uH = uˆ1N1+
uˆ2N2. Note that uˆ1 and uˆ2 are the unknowns we need to solve for. Note that the test function v= vˆ1N1+ vˆ2N2
as vˆ0 = 0 at the Dirichlet boundary condition. As previously, uˆi is the nodal displacement at node i and Ni
is the corresponding shape function, defined for the 2 element bar as given previously in Eqn. (12.4.2) and
derivatives in Eqn. (12.4.2) and derivatives in Eqn. (12.4.2), although note that in this problem N0 is not
needed. The notation here is simplified if we use vectors and matrices so:
uδ =
[
uˆ1 uˆ2
][N1
N2
]
= uˆT N
vδ =
[
vˆ1 vˆ2
][N1
N2
]
= vˆT N
With:
72
NT =
{
[2x/L 0] 0≤ x < L2
[2(1− x/L) 2x/L−1] L2 ≤ x≤ L
And:
NTx =
{
[2/L 0] 0≤ x < L2
[−2/L 2/L] L2 ≤ x≤ L
Now consider the first element on the LHS:
EA
∫ L/2
0
vxux dx =EA
∫ L/2
0
∂
∂x
(
vˆT N
) ∂
∂x
(
uˆT N
)
dx
=EA
∫ L/2
0
vˆT
(
NxNTx
)
uˆ dx
=EAvˆT
(
NxNTx
)
uˆ
∫ L/2
0
dx
=EAvˆT
[ 4
L2 0
0 0
]
uˆ · L
2
=vˆT
2EA
L
[
1 0
0 0
]
uˆ
Similarly, for the second element we have for the LHS:
EA
∫ L
L/2
vxux dx =EAvˆT
(
NxNTx
)
uˆ
∫ L
L/2
dx
=EAvˆT
[ 4
L2 − 4L2
− 4L2 4L2
]
uˆ · L
2
=vˆT
2EA
L
[
1 −1
−1 1
]
uˆ
Combining both for the complete LHS expression:
EA
∫ L
0
vxux dx = vˆT
2EA
L
[
2 −1
−1 1
]
uˆ
Now consider the RHS:
EA [vux]
L
0 = v
δ (L)EAux(L)− vδ (0)EAux(0)
The second term on the RHS goes to zero due to the zero Dirichlet BCs, the first term simplifies when we
apply the Neumann BC giving:
EA [vux]
L
0 =−vδP =−vˆT NP =−vˆT
[
0
1
]
P
Combining both sides we have:
73
vˆT
2EA
L
[
2 −1
−1 1
]
uˆ = −vˆT
[
0
1
]
P (311)
2EA
L
[
2 −1
−1 1
]
uˆ = −
[
0
1
]
P (312)
We can solve Eqn (312) explicitly for the nodal displacements either by inverting the matrix:
uˆ =− PL
2EA
[
2 −1
−1 1
]−1[0
1
]
=− PL
2EA
[
1 1
1 2
][
0
1
]
=− PL
2EA
[
1
2
]
(313)
or by noting from the upper row of Eqn. (312) that 2uˆ1 = uˆ2, and substituting that in the lower row to gain
−uˆ1+2uˆ1 =−PL/2EA.
Calculate the stress at a specific location and compare to the exact solution We will now use this
solution to compare with the exact solution at x = L/4 and x = 3L/4 as a function of EA. For the analytical
solution we have σ = Eux with ux = −PEA . As the strain is constant in this case, stress will also be constant,
ie:
σL/4 = σ3L/4 =
−EP
EA
=
−P
A
For the numerical solution we have from Eqn. (313) uˆ2 = 2uˆ1. Solving Eqn. (312) we have:
2EA
L
(−uˆ1+2uˆ1) =−P → uˆ1 = −PL2EA & uˆ2 =
−PL
EA
Now for the stresses we have:
σL/4 =
E(uˆ1− uˆ0)
L/2
=
E(−PL2EA )
L/2
=
−P
A
Similarly,
σ3L/4 =
−P
A
The numerical solution is identical to the analytical, since the analytical solution is linear and linear basis
functions can represent nodal values exactly for a quadratic solution, and inter-element values exactly for
linear solutions.
12.4.4 Predicting Stresses in a Helicopter Blade
A helicopter blade may be modelled as a one dimensional beam, fixed at the hub (r = 0) and free at the
rotor tip (r = L = 4m). Assuming constant cross-sectional area A = 1× 10−3m2, modulus of elasticity of
carbon fibre E = 140 GPa and total length L = 4m. The effects of centrifugal force are modelled as in the
tutorial in Week 13, where q(r) = mrω2 and the frequency of rotation is 5.4Hz.
The displacement distribution u(r) in the beam is given by Poisson’s equation ∂∂ r (EAur)+q(r) = 0.
An analytical solution to the problem.We are given that E and A are constant and q(r) = ρArω2. We
integrate twice to recover the displacement distribution:
∂
∂ r
(EAur)+ρArω2 = 0 → ∂∂ r (EAur) =−ρArω
2
74
Figure 7: Schematic of the discretisation of the helicopter blade
Integrating twice we have:
EAur =−ρAr
2ω2
2
+B → EAu(r) =−ρAr
3ω2
6
+Br+C
Noting we have the boundary conditions u(0) = 0 and EAur(L) = 0, we find that u(0) = 0 → 0 =
0+C → C= 0 and EAur(L) = 0 → 0=−ρAω
2L2
2 +B → B= ρAω
2L2
2 , leading to the analytical
solutions,
u(r) =
ρω2
2E
(
L2r− r
3
3
)
, ur(r) =
ρω2
2E
(
L2− r2)
Determine the FEA approximation. Next we simplify the weak form of the governing equation and
calculate the stiffness and load vector, discretising the solution using two elements and calculating the
displacement at all nodes using linear shape functions, as used in the previous two examples. The blade is
discretised as in Figure 7
The weak form of the governing equation is:
∫ L
0
vrEAurdr =
∫ L
0
vqdr+[vEAur]L0−
∫ L
0
EAvruDr dr (314)
As the Dirichlet and Neumann boundary conditions are zero, this simplifies to:
EA
∫ L
0
vrurdr =
∫ L
0
vqdr (315)
Now consider element 1. Let uδ (r) = uˆ0N0+ uˆ1N1, vδ (r) = vˆ0N0+ vˆ1N1 and qδ (r) = qˆ0N0+ qˆ1N1. For the
linear shape functions, we have:
N0 = 1− 2rL N1 =
2r
L
N0,r =−2L N1,r =
2
L
Substituting into the simplified weak form of Eqn. (315):
EA
∫ L/2
0
∂
∂ r
(vˆ0N0+ vˆ1N1)
∂
∂ r
(uˆ0N0+ uˆ1N1) dr =
∫ L/2
0
(vˆ0N0+ vˆ1N1)(qˆ0N0+ qˆ1N1) dr (316)
From the LHS:
75
EA
∫ L/2
0
(vˆ0N0,r + vˆ1N1,r)(uˆ0N0,r + uˆ1N1,r) dr
= EA
∫ L/2
0
(−2vˆ0
L
+
2vˆ1
L
)(−2uˆ0
L
+
2uˆ1
L
)
dr
=
4EA
L2
(−vˆ0+ vˆ1)(uˆ0+ uˆ1)
∫ L/2
0
dr
=
4EA
L2
(−vˆ0+ vˆ1)(uˆ0+ uˆ1) · L2
=
2EA
L
[
vˆ0 vˆ1
][ 1 −1
−1 1
][
uˆ0
uˆ1
]
Now for the RHS of Eqn. (315), and expressing terms in matrix notation:∫ L/2
0
vˆT NNT qˆ dr =
∫ L/2
0
[
vˆ0 vˆ1
][N0
N1
][
N0 N1
][qˆ0
qˆ1
]
dr
=
∫ L/2
0
[
vˆ0 vˆ1
][1− 2rL
2r
L
][
1− 2rL 2rL
][qˆ0
qˆ1
]
dr
=
∫ L/2
0
[
vˆ0 vˆ1
][1− 4rL + 4r2L2 2rL − 4r2L2
2r
L − 4r
2
L2
4r2
L2
][
qˆ0
qˆ1
]
dr
=
[
vˆ0 vˆ1
][r− 4r22L + 4r33L2 r2L − 4r33L2
r2
L − 4r
3
3L2
4r3
3L2
]L/2
0
[
qˆ0
qˆ1
]
=
[
vˆ0 vˆ1
][ L
6
L
12
L
12
L
6
][
qˆ0
qˆ1
]
=
[
vˆ0 vˆ1
] L
12
[
2 1
1 2
][
qˆ0
qˆ1
]
Combining the LHS and RHS, we have:
2EA
L
[
vˆ0 vˆ1
][ 1 −1
−1 1
][
uˆ0
uˆ1
]
=
[
vˆ0 vˆ1
] L
12
[
2 1
1 2
][
qˆ0
qˆ1
]
2EA
L
[
1 −1
−1 1
][
uˆ0
uˆ1
]
=
L
12
[
2 1
1 2
][
qˆ0
qˆ1
]
Applying the same procedure for element 2, we can then combine both elements to yield:
2EA
L
 1 −1 0−1 2 −1
0 −1 1
uˆ0uˆ1
uˆ2
= L
12
2 1 01 4 1
0 1 2
qˆ0qˆ1
qˆ2

76
Applying u(0) = 0, the first row and column of the stiffness matrix is no longer solved for, and the first row
of the load matrix is not needed giving:
2EA
L
[
2 −1
−1 1
][
uˆ1
uˆ2
]
=
L
12
[
1 4 1
0 1 2
]qˆ0qˆ1
qˆ2

Furthermore, as q = 0 at r = 0, then qˆ0 = 0, giving:
2EA
L
[
2 −1
−1 1
][
uˆ1
uˆ2
]
=
L
12
[
4 1
1 2
][
qˆ1
qˆ2
]
Solving for the nodal displacements:[
uˆ1
uˆ2
]
=
L2
24EA
[
2 −1
−1 1
]−1[4 1
1 2
][
qˆ1
qˆ2
]
=
L2
24EA
[
1 1
1 2
][
4 1
1 2
][
qˆ1
qˆ2
]
=
L2
24EA
[
5 3
6 5
][
qˆ1
qˆ2
]
ie uˆ1 = L
2
24EA (5qˆ1+3qˆ2) and uˆ2 =
L2
24EA (6qˆ1+5qˆ2).
Calculate the stress at a specific location and compare to the exact solution Let’s now look at the stress
at x = L/4 and x = 3L/4 as a function of EA, and compare this to the exact solution and comment on the
agreement.
Solution
Analytical solution:
σ = Eur
Where ur =
ρω2
2E (L
2− r2)
Now:
σ(L/4) = Eur(L/4) = E
ρω2
2E
(
L2− L
2
16
)
=
15ρω2L2
32
σ(3L/4) = Eur(3L/4) = E
ρω2
2E
(
L2− 9L
2
16
)
=
7ρω2L2
32
Numerical solution:
Calculate nodal displacements noting that q = ρArω2 (ie qˆ1 = ρAω
2L
2 and qˆ2 = ρAω
2L):
77
uˆ1 =
L2
24EA
(5qˆ1+3qˆ2) =
L2
24EA
(
5
ρAω2L
2
+3ρAω2L
)
=
11
48
ρω2L3
E
uˆ2 =
L2
24EA
(6qˆ1+5qˆ2) =
L2
24EA
(
6
ρAω2L
2
+5ρAω2L
)
=
1
3
ρω2L3
E
Now to calculate the stresses:
σ1 = E
uˆ1− uˆ0
L/2
=
2E
L
11
48
ρω2L3
E
=
11
24
ρω2L2
σ2 = E
uˆ2− uˆ1
L/2
=
2E
L
ρω2L3
E
(
1
3
− 11
48
)
=
5
48
ρω2L2
As the solution is quadratic, whereas the inter-element values are reconstructed from the linear profile, the
stress will not be exact since it is constant in an element. The nodal values of u are exact however, since the
solution is quadratic and the basis functions linear.
12.5 Conclusions
This chapter has outlined the foundations of the Finite Element Method and application in the analysis of
one dimensional, steady systems. We have seen how the methodology is founded on the weak form of
the governing equations, and has quite a different basis to the finite difference approach covered previously.
The remit of the FEM is to set the nodal errors in the approximated system to zero, the approximated system
being the original system multiplied by an arbitrary test function, integrated over the whole domain, and
then approximated using nodal values and basis functions. This results in a numerical method which can
deliver exact results at the node when compared to analytical solutions which are of order xp+1 where p is
the order of the basis functions (p = 1 for linear). The basis functions can be used to interpolate between
nodal results to give an approximation within an element, these approximations are exact for analytical
solutions up to the order of the basis function.
The Finite Element Method has become the default method to undertake stress analysis, and we have given
several simple one dimensional examples of how to apply the FEM to Finite Element Analysis (FEA) and
how to post-process the results to provide further quantities of physical interest such as stresses.
78
13 Course Summary
The goal of this course is to provide you with (i) an intuitive understanding of partial differential equations
which may be governed by time-dependent diffusion, waves or steady diffusion, (ii) an arsenal of mathe-
matical tools which you may employ to gain very useful analytical solutions and (iii) the underpinnings of
numerical methods along with their advantages and disadvantages. Importantly, we hope that you see the
intimate relationship which exists between theoretical and numerical solutions - in practical engineering
problems you want to make use of both of these. Analytical solutions guide design in that they provide
insight into the general behaviour of the often more complex real problem. The numerical solutions, once
verified and validated, may then provide quantitative predictions needed to finalise a given design. Theoret-
ical solutions can provide the necessary verification point, and in most practical cases experiments provide
final validation.
So what have we learned? In terms of the understanding of the behaviour of physical systems we have
covered:
• Qualitative classification of PDEs (diffusive, wave like, steady diffusion)
• Diffusion of heat in 1D (unsteady) and 2D (steady)
• Diffusion of species in 1D (unsteady) and 2D (steady)
• Vibrating beams: waves in vertical displacement, shear and torsional
• Displacement, shear and stresses in beams
• Electric fields and potential flows in 2D (Laplace)
You can now solve all of these problems! In order to solve these problems we have learned new analytical
approaches including:
• Classification: Parabolic, Hyperbolic, Elliptic
– Initial Value Problems
– Boundary Value Problems
• Exact Analytical Solution:
– Separation of Variables
– Fourier Series
– Fourier Integrals/Transforms (infinite domains)
– Laplace Transforms (semi-infinite domains, time varying forcing/source)
For more complicated problem you can now tackle these numerically with the following new numerical
skills and associated understanding:
• Taylor Series
• Interpolation
• Order of Accuracy
• Verification and Validation
79
• Discrete representation of the governing equation
• Determining Stability
• Implicit and Explicit methods
• Finite difference methods
• Finite element methods
• Much improved Matlab capabilities!
You now know a range of exact solutions to the most important PDEs in Engineering. You now know how
to solve those PDEs numerically to gain solutions where the analytical method will not work. You also
know how to verify the implementation and validate those algorithms. Most importantly you appreciate the
limitations of both approaches and understand where and when care must be taken.
Remember to always consider the simplest possible analytical solution to validate any numerical method
you employ, be it in:
• Stress analysis
• Fluid dynamics
• Control and Dynamics
Exploit large data effectively by representing it in frequency/wavenumber space. Make sure that during
assignments in your other UoS that you aim to develop an intuitive ‘feel’ for physical problems by making
use of simple governing models such as the ones presented here. Use your skills learned here to help
you complete complex assignments, which possibly require non-analytical solutions – can you make a
numerical solution? plot the solution, animate it, does the solution make sense with your understanding?
Does it match your expected boundary conditions?
As you move further in your degree program, we hope that this UoS will benefit you beyond the relatively
narrow scope of the material explicitly covered. Embrace the additional capabilities which come with
the use of computational methods - you are now able to solve problems easily which cannot be solved
analytically by the greatest minds on the planet.
80
Appendix A Derivation of the Heat Equation
FL =
∫ z0+δ z
z0
∫ y0+δy
y0 −kTx|x0dydz FR =
∫ z0+δ z
z0
∫ y0+δy
y0 −kTx|x0+δxdydz
Ht
(x0,y0,z0)
(x0+δx,y0+δy,z0+δ z)
Figure 8: Control volume for the derivation of the heat equation
Figure 8 shows the control volume used to derive the heat equation. The equation is a balance between net
flux of heat into the control volume and accumulation of heat within the control volume.
In this derivation it is assumed that the temperature profile only varies in the x direction and time, i.e. T (x, t),
and that the control volume has dimensions of δx0× δy0× δ z0 and has bottom left corner coordinates at
(x0,y0,z0) through to a top right coordinate (x0+δx0,y0+δy0,z0+δ z0).
The governing equation is derived by stating that:
(Rate of change of heat)=(heat flux at the left hand face)-(heat flux at the right hand face)
Ht = FL−FR (317)
The total heat within the control volume is defined as
H =
∫ z0+δ z
z0
∫ y0+δy
y0
∫ x0+δx
x0
ρσT dxdydz. (318)
The rate of change of heat is:
∂
∂ t
H =
∂
∂ t
[∫ z0+δ z
z0
∫ y0+δy
y0
∫ x0+δx
x0
ρσT dxdydz
]
. (319)
and swapping the order of the derivative and integral gives:
Ht =
∫ x0+δx
x0
ρσTtdxδyδ z. (320)
The left hand flux into the domain is given by Fourier’s law F =−kTx, integrated over the whole face:
FL =
∫ z0+δ z
z0
∫ y0+δy
y0
−kTx|x0dydz (321)
FL =−kTx|x0δyδ z (322)
For the right face:
FR =−kTx|x0+δxδyδ z (323)
Now we make the balance of the terms:
81
Ht = FL−FR (324)∫ x0+δx
x0
ρσTtdxδyδ z = −kTx|x0δyδ z+ kTx|x0+δxδyδ z (325)∫ x0+δx
x0
ρσTtdx = kTx|x0+δx− kTx|x0 (326)
Next we assume that the control volume becomes infinitesimally small, i.e. δx→ 0. The following limits
then apply:
∫ x0+δx
x0
ρσTtdx = ρσTtδx, (327)
kTx|x0+δx− kTx|x0 = k(Tx|x0+δx−Tx|x0) = kTxxδx (328)
Putting these back into Eqn. (326) gives:
Ttδx = kTxxδx, (329)
ρσTt = kTxx, (330)
Tt =
k
ρσ
Txx. (331)
This could alternatively be written as:
Tt = c2Txx,c2 =
k
ρσ
. (332)
Appendix B Derivation of the Fourier Sine Series Coefficients
The wave equation with Dirichlet boundary conditions has a general solution of the form:
u(x, t) =
∞
∑
n=1
un(x, t) =
∞
∑
n=1
(Bn cosλnt+B∗n cosλnt)sin
npi
L
x. (333)
and likewise the homogeneous heat equation with Dirichlet boundary conditions has a general solution of
the form:
u(x, t) =
∞
∑
n=1
Bn(sin
npi
L
x)e−λ
2
n t (334)
In order to satisfy the initial condition, we obtain from both (under the assumption of zero initial velocity
for the wave equation):
82
u(x,0) =
∞
∑
n=1
Bn sin
npi
L
x = f (x) (335)
Now we need to determine the Bn coefficient for general cases with an initial condition given by f (x).
Firstly, in order to use Fourier series, we need to extend our function (if it is not periodic) f (x) to be a
periodic function where in this case the function has period 2L, i.e.:
f (x+2L) = f (x). (336)
For example, in the guitar string case, even if we are only concerned about 0< x< L, we still have to extend
the initial ‘triangle wave’ to a periodic function (Section 5.2 of the Lecture slides).
The first step is to multiply both sides of Eqn. (335) by sin mpiL x (note that sin
mpi
L x is also a periodic function
in the range 0 < x < L), and integrate from −L to L. Following this we obtain:
∫ L
−L
f (x)sin(
mpi
L
x)dx =
∫ L
−L
[
∞
∑
n=1
Bn sin(
npi
L
x)
]
sin(
mpi
L
x)dx, (337)
where n= 1,2,3, · · · · · · and m= 1,2,3, · · · · · · . Integrating term by term, we see that the right side becomes:
∞
∑
n=1
Bn
∫ L
−L
sin(
npi
L
x)sin(
mpi
L
x)dx. (338)
Using trigonometric identities, we have:
sinα sinβ =
1
2
[cos(α−β )− cos(α+β )], (339)
then Eqn. (338) becomes:
∞
∑
n=1
Bn
[
1
2
∫ L
−L
cos[(n−m)(pi
L
x)]dx− 1
2
∫ L
−L
cos[(n+m)(
pi
L
x)]dx
]
. (340)
The second term equals zero since:
1
2
∫ L
−L
cos[(n+m)(
pi
L
x)]dx =
1
2
[
L
pi(n+m)
sin[
(n+m)pi
L
x]
]L
−L
= 0, (341)
thus Eqn. (340) is simplified to:
∞
∑
n=1
Bn
1
2
∫ L
−L
cos[(n−m)(pi
L
x)]dx. (342)
There are now two possible cases. For n 6= m, the integration in equation (342) is zero for the same reason
as Eqn. (341). For n = m, there is now only one single term in the summation, i.e. Eqn. (342) becomes:
Bm
1
2
∫ L
−L
cos[0]dx = Bm
1
2
[x]L−L = BmL, (343)
83
and substituting (343) into (342),then going back to (337), we obtain:
∫ L
−L
f (x)sin
mpi
L
xdx = BmL. (344)
Rearranging, we obtain:
Bm =
1
L
∫ L
−L
f (x)sin
mpi
L
xdx (345)
Since n and m are equivalent, from now on we shall use n, therefore obtaining:
Bn =
1
L
∫ L
−L
f (x)sin
npi
L
xdx. (346)
An analogous derivation can be used to derive the cosine coefficients.
Example 1: For the square wave initial condition in the tutorial in Week 3:
f (x) =
{
A 0≤ x < L,
−A −L≤ x < 0.
Bn =
1
L
∫ L
−L
f (x)sin(
npi
L
x) =
1
L
(∫ 0
−L
−Asin(npi
L
x)dx+
∫ L
0
Asin(
npi
L
x)dx
)
,
=
1
L
(
AL
npi
cos(
npi
L
x)
∣∣∣∣0
−L
− AL
npi
cos(
npi
L
x)
∣∣∣∣L
0
)
,
=
AL
npiL
[cos(0)− cos(−npi)− cos(npi)+ cos(0)],
=
2A
npi
(1− cos(npi)) n = 1,2,3, · · · · · ·
(347)
Example 2: For the saw-tooth shape initial condition:
f (x) =
Ax
L
−L < x≤ L, (348)
where A =−0.22. The coefficients Bn are:
84
Bn =
1
L
∫ L
−L
f (x)sin(
npi
L
x)dx =
2
L
∫ L
0
f (x)sin(
npi
L
x)dx
=
2
L
∫ L
0
A
x
L
sin(
npi
L
x)dx
=
2A
L2
∫ L
0
xsin(
npi
L
x)dx
=
2AL
L2npi
[[
xcos(
npi
L
x)
]L
0
−
∫ L
0
cos(
npi
L
x)dx
]
=
2AL
L2npi
(Lcos(npi)−0)
=
2A
npi