Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
AMME2000 & BMET2960/9960:Engineering Analysis
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 .
- 经济Economics
- 金融Finance
- 会计Accounting
- 工商管理Business Management
- 计算机computer science
- 商业分析business analysis
- 法律Law
- 管理Management
- 市场营销Marketing
- 哲学Philosophy
- 心理学Psychology
- 物理Physical
- 统计学statistics
- 电气工程electric engineering
- 化学chemistry
- 生物biology
- 数学mathematics
- 语言学linguistics
- 土木工程civil engineering
- 工程力学engineering mechanics
- 艺术art
PROFESSIONAL HELP!
在线咨询
立即下单
Wechat：THEend8_ 
Email：