MENG21712
10 Credits
Modelling with computers (TB1)
Dr. A.G.W. Lawrie
B Stream function descriptions of incompressible ﬂow
C Derivation of thin membrane mechanics
D Numerical convergence
1
Introduction
The aim of Modelling 2 is to introduce students to the ‘art’ of using a range of
modelling techniques to understand and ﬁnd solutions to engineering problems.
The part taken in Teaching Block 1 (TB1) this year is a short project where
the task is to model three closely related physical problems using Matlab. The
followon component in Teaching Block 2 involves the design/build/test of an
amphibious vehicle. The following course notes pertain to TB1.
1
.1 Computing for Mechanical Engineers
In ComputerBased Modelling (MENG11511) in Year 1 you have been trained
in the use of Matlab to encode a list of tasks into a procedural algorithm that
can be executed on computer, using the software package Matlab. Matlab is
a convenient ‘Swiss army knife’ spanning everything from calculations you might
perform on a pocket calculator to complex numerical problemsolving. The TB1
component of this course (MENG21712) will teach the principles of solving real
world physical problems whose structure is governed by partial diﬀerential equa
tions. Our aim will be to develop the skills, tools and programming techniques
to express engineering problems in a form that a computer can implement for us,
then you will design/build/test a method that with small adjustments can solve
a broad range of such physical problems.
2
1
.2 Educational Aims
Nowadays especially, enabling a computer to work for you is a core skill in any
practicing Engineer, and more broadly across all numerate, technical and problem
solving professions. Some students retain apprehensions regarding the use of com
puters as a tool under their control to solve their problems. This is likely to be
detrimental to their careers, whether they be in Engineering or elsewhere. This
course aims to give you experience of using computers to produce credible Engi
neering solutions to physical problems, and should allow you to grow comfortable
with the tools, algorithms, implementation strategies and some level of technical
background. Engaging with this course should help you overcome any remaining
misgivings.
Any form of systematic procedure to solve a problem is called an algorithm,
and it can be broken down into several small steps. Individually each step may
be simple and easily understood, but as the length of the algorithm increases
then the potential for unexpected behaviour quickly grows. Programming is a
skill that develops through practice because it takes repeated exposure in order
to restructure your own thinking so that it matches the rigorously systematic
execution of a computer. Although devising an implementation of an algorithm is
often tricky, understanding errors discovered when testing it is the real challenge
in programming. With practice and experience, both things this course aims to
provide, your errors will become less serious and less frequent, and your conﬁdence
will grow.
On this course we aim to provide a supportive environment for the development
of programming skills. To manage your workload, a stagegate in week 5 will
ensure that you have produced code that meets the passstandard of the course
by the halfway point. During the ﬁnal 5 weeks you will be encouraged to improve
your code’s performance, and pointed towards some techniques and algorithms
that can help you achieve that. The ﬁnal submission (in week 10 to avoid collision
with other deadlines) will comprise your code and a report justifying your use of
algorithms, implementation choices and presenting the results of your performance
optimisation.
The university assesses coursework on a 21point scale that matches marks to
descriptors of achievement, and ﬁrstclass (¿70%) marks are awarded for going
beyond what is explicitly taught in the course. To allow you to demonstrate
this achievement, in the second half of the course support will still be available
but you will be given guidance pointing you towards ways of improving your code
rather than oﬀering explicit teaching. As a further incentive, a prize is awarded
for the bestperforming code!
3
1
.3 Computing Facilities
The teaching computer resources for the Faculty of Engineering are located in the
Merchant Venturers Building (MVB) on Woodland Rd. Before the start of this
course all students should have a login and password for the Faculty computers.
If not you will need to obtain one from IT services. Laboratory sessions will take
place in both of the PC labs, rooms 1.07 and 1.08 in the MVB, and all
students enrolled on the course will be given a 2hour slot in one of these rooms.
In the laboratory there will be teaching assistants who can assist you wherever
you have diﬃculty in completing your tasks. You must attend the laboratory slot
you are given in your timetable because we have calculated your allocation to
either group 1 or group 2 on the available capacity of each room, and distributed
teaching assistants accordingly. Swapping session without permission from
the unit director is not permitted. The University also provides nocost
oﬀcampus access to Matlab in a variety of ways: ﬁrstly the software itself can
be downloaded onto your personal laptop using a licence key purchased by the
University and available through IT services, (www.bristol.ac.uk/software) and
secondly the University provides Matlab as one of the installed packages on
the Remote Desktop service (www.bristol.ac.uk/studentdesktop). The various
editions of Matlab all operate in slightly diﬀerent ways and code developed on
one version cannot be guaranteed to work on another. A common reason for this
is that certain Matlab toolboxes may be installed in the version you have used to
develop your code, but not on the version used to assess your code. This becomes
a problem should you have called a function that is contained only in this toolbox.
To avoid such problems, a list of admissible Matlab functions has been provided
in §5.7. Before submitting your ﬁnal work, it is your responsibility to
ensure that your code only uses functions from the prescribed lists.
1
.4 Attendance and Study
Attendance at laboratory sessions is the best means of obtaining assistance and
feedback from teaching assistants, and is also compulsory. Attendance lists will
be recorded for each laboratory session, and in total you will be expected to
spend approximately 70 hours on this part of the course, including contact time
and independent study. The introductory lectures will cover the mathematical
background you need to understand the physical problem and the approach we
will use to convert the problem into a form the computer can easily execute. The
laboratory sessions will be performed in groups of 2 or 3 students that you will
form in week 1. Each group will be required to hand in a group report.
4
1
.5 TB1 Assessment
Although the work is performed in groups, submissions are made per group, and
each package of submissions is assessed as a unit, marks are awarded to indi
vidual students. In line with uuniversity policy on group work, conﬁdential
peer assessment is used, where appropriate, to redistribute marks, though it
is anticipated that in most cases there will be no need to do so.
There is a stagegate that must be reached by your laboratory session on
Monday 28th October 2019, and this ensures that your group’s code meets
the passstandard for this element of the course well in advance of the ﬁnal dead
line. The assessment will involve demonstrating your group’s code to teaching
assistant during the laboratory and they may also ask a few questions to test your
understanding. For the purposes of this assessment, your code should produce
simple graphical output, showing a ﬁnal solution and a splot showing the path
taken to converge to it. If satisfactory, the teaching assistant will sign each in
dividual’s stagegate assessment form and assign your group a unique identiﬁer
number. It is your responsiblity to return the completed form to the unit director.
For the ﬁnal deadline, 5pm, Friday 6th December 2019, your submission
must include a code (a single Matlab .m ﬁle containing ASCII text only) and a
report (a single .pdf ﬁle of no more than 10 pages produced using the L^{A}T_{E}X doc
ument preparation system), submitted electronically on Blackboard. Both will
pass through a sequence of manual and electronic honestychecks to guard against
plagiarism. We may also perform random spotchecks, and request constituent
components of your report, eg. ﬁgures, .tex .bib ﬁles and so forth. Your report
will be marked conventionally against the 21point scale, and your code will be
tested exhaustively in an automated script to assess its performance on a range
of testcases. For the purposes of this assessment your code must not produce
graphical output or require user input.
Marks for each part of the course will be accumulated and an individual total
compiled at the end of the year. The work in TB1 contributes in total 60% of
the mark for this course, and this is subdivided into a 10% weighting towards the
week 5 stagegate, 50% towards your code and 40% to the report.
There is no timed examination for this course. Note that no additional credit
is awarded for the sophistication of graphical output: graphical user interfaces
are elegant but timeconsuming to develop, and in this course the focus is on
algorithmic eﬃciency rather than aesthetics of output. In recent years Matlab
has tried to improve its underlying execution speed (which is slow relative to other
coding platforms) and now oﬀers a variety of tools to mitigate its disadvantages.
Compilation to binary of any form is strictly forbidden (since there is no way
to verify authorship or quality of content), and the only admissible submission
format is a single .m ASCII text ﬁle.
5
1
.6 Rubric and intended learning outcomes
At the end of TB1 students should be able to:
•
•
represent realworld physical problem in mathematical notation using partial
derivatives
transform partial diﬀerential equations into a sequence of instructions a
computer can execute
•
•
•
•
•
devise systematic ways of testing instruction sequences (code) for errors
repackage code to conform to a speciﬁcation
graphically present results obtained by using your code
use the L^{A}T_{E}X document preparation system
summarise their work in a written report
A ﬁrstclass student will also be able to:
•
•
•
•
•
optimise the implementation of an algorithm for eﬃciency
optimise the choice of algorithm to improve scaling performance
critically evaluate their choices, articulate their reasoning in some depth
show evidence of independent study of alternative algorithms
evaluate their code’s performance over a suitably chosen range of testcases
6
Your code will primarily be evaluated by automated script that will perform
a timing test on a range of test cases. This will test your code for
•
•
•
•
successful execution
returning correct results
correct error handling
scaling performance
In the event of unsuccessful execution, incorrect results, or incorrect error han
dling, the automated script will continue on to the next case in the test suite. In
the event of successful, correct results, a time threshold will decide if the scaling
performance exhibited by your code warrants running a more complex test case.
For consistency, codes from all groups will be tested on a common machine. All
code will be evaluated against a benchmark implementation and graded accord
ingly.
The report will be marked according to the University’s 21point scale for
coursework, and the relevant interpretation of the scale for this task is tabled
beside the descriptor in ﬁgure 1.
There is no L^{A}T_{E}Xtemplate provided for this course, you should use your own
choice of formatting. In general terms, your marker will be seeking to see:
•
•
•
•
•
•
•
clarity of structure
lucidity of prose: brevity without castration
sections and paragraphs used as intended
correct typesetting in L^{A}T_{E}X
ﬁgures with correctly labelled axes
ﬁgures with captions
a relevant bibliography, cited in the body of text
Outline the physics of the problems, but focus on the algorithms developed in
your code, key optimisations, and evaluate why and how well they work. One
thing I strongly encourage is to interleave results and discussion together
in a single section. Undiscussed results have little meaning, discussion detached
from results involve a lot of tedious scrolling…
Bear in mind that your ﬁnal mark for TB1 has three contributions, 10% from
the week 5 stagegate, 50% from the ﬁnal code submission and 40% from the
report, so we are placing emphasis on quality and behaviour of your code over
other elements of assessment.
7
1820: Report reaches publication standard in
terms of content for a reputable academic journal,
with little other than format adjustment to ﬁt the
style requirements of the journal.
1517: Clear evidence of originality, ingenuity and
depth of understanding apparent from report content,
with lucid written communication, and with the
message intended of diagrams and ﬁgures
unambiguous at a ﬁrst reading.
1214: Clear evidence in the report of successful
execution of the required tasks, expressed with
written clarity and good support from diagrams and
ﬁgures, but lacking in originality when seeking
solutions to problems.
911: Evidence of some thought on the problems at
hand from a thorough reading of the report, but the
message is delivered without adequate clarity in text,
diagrams or ﬁgures, and there is little or no evidence
of originality or ingenuity.
68: Evidence of negligible effort in execution of the
project and/or preparation of the report, widespread
lack of clarity in communicating a message and clear
evidence of insufﬁcient thought on the problems at
hand.
05: Partial, invalid or nonsubmission of report,
clear evidence of dysfunctional execution of the
project and/or preparation of the report.
Figure 1: Interpreting the University 21point marking scale for coursework.
8
1
.7 Reference material
The main reference text which applies to this course is [1]. For the programming
aspects of the course, see also [2] and [3]. In addition there are many good web
based resources which include Matlab tutorials. These can be found by searching
the web using a good search engine, as well as searching within the help function
of Matlab. For the modelling aspects of the course, a further useful reference is
[
4]. All necessary materials are provided in this set of course notes, but reading
more broadly is one indicator of going beyond what is explicitly taught. In your
report references in L^{A}T_{E}X should be cited using \cite{Hahn} in the following
manner,
\
\
begin{thebibliography}{99}
bibitem{Hahn} B.D.~Hahn, 2007. {\it Essential MATLAB
for engineers and scientists.} 3rd Edition, Elsevier.
\
{
\
bibitem{AUSCH} M.~Austin and D.~Chancogne, 1999.
\it Engineering programming: C, Matlab, Java}, John Wiley.
bibitem{HUNT} B.R.~Hunt~{\it et al}, 2001. {\it A guide to
MATLAB : for beginners and experienced users}, Cambridge University Press.
bibitem{FAUS} L.V.~Fausett, 1999. {\it Applied numerical
analysis using MATLAB}, Prentice Hall.
\
\
end{thebibliography}
resulting in formatted output as follows:
References
[
[
[
[
1] B.D. Hahn, 2007. Essential MATLAB for engineers and scientists. 3rd Edi
tion, Elsevier.
2] M. Austin and D. Chancogne, 1999. Engineering programming: C, Matlab,
Java, John Wiley.
3] B.R. Hunt et al, 2001. A guide to MATLAB : for beginners and experienced
users, Cambridge University Press.
4] L.V. Fausett, 1999. Applied numerical analysis using MATLAB, Prentice
Hall.
1
.8 Document organisation
In this part of the course we will look at three physical problems that arise in
diﬀerent areas of engineering: Thermodynamics, Fluid dynamics and Structural
9
Mechanics, but are nonetheless related. The objective is to solve each of the given
problems using the same solution technique and develop a single common code able
to solve all three problems. This is possible because the diﬀerences between them
are simply in the physical context, the underyling partial diﬀerential equation
shares the same form and so they can all be solved using a common method.
The remainder of this document is organised as follows: The section on so
lution techniques, §2 walks you through an approach to computermodelling of
the problems posed. We demonstrate in §2.1 how to take an equation derived
from physical principles and solve it on a simple onedimensional problem using
analytical methods. Then in §2.2 we show how to discretise the equation, and in
§
2.3 illustrate how a numerical solution strategy can be implemented in Matlab.
Thereafter, in §3 the speciﬁc remit for each problem is discussed in detail. In
section §4 there are suggestions for ways to optimise your code for greater eﬃ
ciency, and ways to maintain eﬃciency as the problem resolution grows. Finally
§
5 contains important information about the ﬁnal submission speciﬁcation and
format. In the appendices, we provide background information on how several
diﬀerent physical principles result in the same form of partial diﬀerential equa
tion, and additional technical details on why numerical methods get less eﬃcient
as the problem size grows. The information in the appendices will inform your
report writing and ﬁnalisation of the code towards the end of your project, but it
is not necessary to read this material to begin work on your code.
The ﬁrst step to make progress towards a solution is to copy the code for the
rope bridge (demonstrated in lectures) into a text ﬁle that Matlab can run. A
good approach in computer programming is to keep the programming running but
slowly migrate it slowly from something simple that you already have working, to
something more complex that does more of what you want. A sensible approach
would be to modify the code in some simple ways so that instead of solving a 1D
UDL (the rope bridge problem) it solves the analogous problem of a 2D UDL on
a thin membrane. This fairly simple step opens the door to tackling all of the
three projects listed the section on problem descriptions. The next step is then
to modify this 2D code so that it addresses each specﬁc problem.
2
Solution techniques
This section introduces the common approach to solving all the problems posed.
Most engineering problems can be decomposed into simple rules that govern be
haviour. In many familiar cases these rules are based on Newton’s laws coupled
with a ‘constitutive relation’ that comes from physics, for example, Hooke’s law
to relate stress and strain, speciﬁc heat capacity to relate heat energy and tem
perature, or the perfect gas law to relate pressures and densities in a gas.
1
0
Figure 2: Uniformly distributed load on a rope.
Here, we’ll look at an inuitive problem to get started. A rope bridge is a
centuriesold engineering solution to crossing a gorge, long before our local hero,
Isambard Kingdom Brunel, designed the ﬁrst suspension bridge in the modern
style to span the Avon gorge at Clifton. Supported only at each end, and only
able to resist axial stress (and does not resist stress due to bending) a rope bridge
hangs under its own weight, as shown in diagram 2. The rule governing the
vertical displacement, φ, of the system is as follows:
^{2}φ
k^{∂}_{∂ }
_{x}_{2 }= mg,
(1)
where m is the mass per unit length, g is gravity and x is the ordinate in the
spanwise direction.
2
.1 Analytical solution
Equation (1) is a diﬀerential equation we can solve directly by integrating twice.
The ﬁrst integration gives,
^{∂}_{∂}^{φ}_{x }=
mg
_{k }x + C_{1},
(2)
where C_{1} is a constant of integration. Integrating again we have,
mg
φ(x) =^{ 1}_{2}_{ k} x + C_{1}x + C_{2}.
2
(3)
We can set the constants C_{1} and C_{2} by applying boundary conditions. In the
trampoline problem, and as is common in many beam deformation problems, the
boundaries have ‘Dirichelet’ boundary conditions, the prescribed displacement
φ = 0. We can solve equation (3) at both boundaries x = 0 and x = L:
1
2
= C_{2 }
mg
k
2
x = 0, φ = 0 : 0 =
0 + 0C_{1} + C_{2 }
(4)
(5)
0
Similarly, and using the fact that C_{2} = 0,
x = L, φ = 0 : 0 =
1
2
mg
2
L + C_{1}L + 0
(6)
(7)
k
C_{1 }
= −^{1 }
^{m}_{k}^{g}L.
2
1
1
So the deformation can be described as,
1
2
mg
k
1 mg_{L}_{x }_{+ }_{0 }
2 k
2
φ(x) =
x −
(8)
(9)
mg
=
_{2}_{k }x (x − L) ,
which is consistent with the parabolic sag shown in diagram 2.
2
.2 Numerical solution
As engineering problems become more realistic they become more complex and it
becomes progressively more diﬃcult to ﬁnd analytical solutions. Instead we try
to approach problems in a way that can handle more general cases, and we resort
to numerical methods of the type we will study here. Analytical calculations
still have a place in our understanding, particularly here in checking that our
computer software is providing satisfactory answers. To use a computer to model
an engineering rule like equation (1), we can discretise gradients using ‘ﬁnite
diﬀerences’, essentially using the original deﬁnition on which Calculus is based,
f^{0}(x) =^{ f}^{(}^{x }^{+}^{ h}^{)}^{ −}^{ f}^{(}^{x}^{) }
,
(10)
h
but instead of allowing h → 0 we let h → ∆x, some small nonzero value. If
we discretise the x axis into small segments, like links in a chain (rather than a
continuous rope), then the point where two links hook together we can in general
call position i, the previous join is position i−1 and the next join is position i+1.
Thus a gradient at a location halfway between i and i + 1 is given by,
ꢀ
ꢀ
ꢀ
ꢀ
∂
∂
φ
x
_{≈ }φ_{i}_{+}_{1} − φ_{i}_{ . }
∆x
(11)
i+ ^{1}_{2 }
Often, as in this case, the rules of engineering problems involve second derivatives,
the ‘derivative of the derivative’. Thus,
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
∂φ
∂φ
∂x i−_{ 2 }
ꢀ
∂^{2 }
x^{2 }
φ
∂x i+^{ 1}_{2 }
−
∆x
1
≈
.
(12)
∂
i
ꢀ
ꢀ
ꢀ
∂
∂
φ
∂φ
, but we still need
ꢀ
.
We already have the value for
x i+^{ 1}_{2 }
1
2
. Shifting the
∂x i−
indices i, i + 1 backwards to i − 1, i we have,
∂
∂
φ
x
_{≈ }φ_{i} − φ
∆x
^{i}^{−}^{1 }.
(13)
1
i− _{2 }
1
2
q
(x_{i}_{ï}_{1},q_{i}_{ï}_{1})
(x_{i}_{+}_{1},q_{i}_{+}_{1})
(x_{i },q_{i})
x
6
x
1
2
1
2
x_{i}_{ï}_{1 }
xiï ï
x_{i }
xi+ ï
x
i+1
Figure 3: A typical outofplane displacement φ in space.
as can be inferred from the typical representation of displacement φ against x
shown in diagram 3. So now the second derivative is given by,
ꢀ
φi+1−φi
φi−φi−1
∆x
∂^{2 }
φ
ꢀ
ꢀ
ꢀ
−
∆
x
≈
.
(14)
(15)
(16)
∂
x^{2 }
∆x
i
Simplifying, we have,
ꢀ
ꢀ
ꢀ
ꢀ
∂^{2 }
x^{2 }
φ
(φ_{i}_{+}_{1 }− φ_{i}) − (_{2}φ_{i} − φ_{i}_{−}_{1 }
(∆x)
_{i}_{+}_{1 }− 2φ_{i} + φ
≈
^{) }= ^{φ }
^{i}^{−}^{1 }.
2
∂
(∆x)
i
So the discretisation of our rule, equation (1), is as follows:
k^{φ }
_{i}_{+}_{1 }− 2φ_{i} + φ
^{i}^{−}^{1 }= mg.
2
(
∆x)
Written like this, it is expressed as a statement of fact: if we know all the correct
values on the lefthandside then the expression must balance the righthandside.
However in the absence of knowing the correct values, we are stuck. Instead
we seek a procedure that takes some initial guess at what the values should be
and update them until these guesses progressively become more acccurate. This
process is known as convergence on the correct values. To implement such an
approach we need a procedure that moves us from one guess to another, a ‘number
machine’ that takes in our ﬁrst guess and returns our second guess, takes in this
second guess and returns our third guess, and so on, subject to some parameters
that come from the problem. The information we need is encoded in equation (16),
but we need to turn this equation into an update procedure, a formula acting on
the initial guess.
1
3
As a ﬁrst step we rewrite our discrete equation (16), rearranging in terms of
the value of φ at our general position i and obtain,
2
∆x) mg
(
φ_{i}_{+}_{1 }− 2φ_{i} + φ_{i}_{−}_{1 }
=
=
(17)
(18)
k
φ_{i}_{+}_{1 }+ φ
2
(∆x) mg_{. }
φ_{i }
^{i}^{−}^{1 }−
2
2k
This statement now explicitly deﬁnes a value of φ_{i} that is dependent on its neigh
bours φ_{i}_{−}_{1} and φ_{i}_{+}_{1} and details of the problem, in fact we take the average of
surrounding values and account for the local contribution of the uniformly dis
2
(
∆x) mg
tributed load as an additional constant term
and decreases with stiﬀer ropes.
that increases with the load
2
k
Equation (18) holds true when the elastic displacement responses φ_{i}_{+}_{1}, φ_{i} and
φ_{i}_{−}_{1 }exactly balance the applied force, resulting in a static equilibrium. However it
can also form the basis of a numerical method that starts marching from an initial
guess and converges towards this exact equilibrium solution. If the marching index
is n, then,
n
i+1
n
2
φ
+ φ
(∆x) mg_{. }
n+1
^{i}^{−}^{1 }−
φ_{i }
=
(19)
2
2k
This update formula only exactly corresponds to the original problem statement
n+1
n
(
n
16) when there is zero diﬀerence between φ
and φ and so the update indices
and ^{n}^{+}^{1} are interchangeable and therefore can be safely ignored. We take this
to be the deﬁnition of ‘convergence’. Prior to convergence, the current values of
φ will not balance equation (16). By careful construction of the update formula
(
and behind this lies an area of mathematics known as numerical analysis) the de
gree of imbalance of the equation is guaranteed always to reduce until it becomes
asymptotically close to zero. In practice, we have ﬁnite patience and tend to trun
cate our marching index when the degree of imbalance is within some acceptable
(
userdeﬁned and problemdependent) limit.
More detail on how and why these methods work is provided below, but in the
meantime one could code this update formula in Matlab.
2
.3 Algorithm and implementation
The numerical method could be described in words as follows:
1
2
3
. set physical constants, eg. gravity, spatial dimensions
. set userspeciﬁed parameters for gridspacing, tolerances
. create arrays to hold old and new values of the solution variable φ
1
4
4
5
6
. set ‘old’ array to starting guess
. set values on the various boundaries according to the speciﬁed problem
. calculate the updated value of φ at all points in the domain, using equation
(
19)
7
8
9
. estimate how far we are from convergence
. copy updated ‘new’ solution to ‘old’ array in preparation for next iteration
. if we are further from convergence than given tolerance, return to step 5,
else continue
1
0. plot the function φ(x) vs. x
This stepbystep algorithm can be transformed into the following Matlab code.
set constants
%
nx=100;
mass=10;
gravity=9.8;
stiff=0.1;
length=1;
phinew=zeros(nx);
phiold=zeros(nx);
dx=length/nx;
force=0.5*dx*dx*mass*gravity/stiff;
%set up while loop
errc=1;
erro=errc;
tol=1e6;
it=0;
while errc/erro>tol
%enforce boundary displacement
phiold(1)=0;
phiold(nx)=0;
%calculate updated displacements
for ix=2:nx1
phinew(ix)=0.5*phiold(ix+1)+0.5*phiold(ix1)force;
end
%calculate gradient of the error
errc=0;
for ix=1:nx
1
5
errc=errc+abs(phinew(ix)phiold(ix));
end
errc=errc/nx;
if it==1
erro=errc;
end
%
swap old and new for next iteration after calculating error
for ix=1:nx
phiold(ix)=phinew(ix);
end
it=it+1;
end
plot(0:dx:lengthdx,phinew);
Note in particular how the convergence tolerance is enforced. The tolerance
tol is a userdeﬁned constant, which in all the exercises is set to 1e6. The error
between the current guess and the exact solution is never computed (since we don’t
a priori know the exact answer). Instead we compute the rate of change of the
error with respect to the iteration count. We do this by taking the absolute value
in each grid cell of the diﬀerence between the current iteration and the previous
one, then summing the contributions from each cell. To obtain a measure that is
independent of the grid resolution nx, then we normalise this total error by the
number of cells so that the error is expressed as a spatial average. We want to
provide a relative measure of convergence that is common to all problems. The
ﬁnal ‘converged’ value is measured with respect to the value at the end of the ﬁrst
iteration, given by the ratio errc/erro. We are thus deﬁning convergence to be
obtained when the rate of change of error per iteration drops from its initial value
to a value at least 6 orders of magnitude smaller.
2
.4 Two dimensions
In more realistic engineering problems we often encounter more than just one
spatial dimension, and constraints that cannot be resolved analytically. In the
following, we consider problems in two dimensions where, like the rope bridge
above, the physics requires discretisation of second derivatives. We generalise the
idea of discretising like links in a chain, to discretising like a ﬁshing net with
square holes. We have to calculate gradients in both the x and y directions, so
we introduce a second index j to identify the grid position in this new direction.
A rectangular mesh with spacing ∆x, ∆y is shown in ﬁgure 4 In many cases the
physics established for a one dimensional problem can be very simply translated
1
6
^{q}i,j+1
^{q}iï1,j q_{i}_{,}_{j }
q_{i}_{+}_{1}_{,}_{j }
^{q}i,jï1
Figure 4: Stencil on a rectangular mesh.
2
∂
φ
into two dimensions. Wherever derivatives_{ ∂}_{ 2} appear, these translate to
x
ꢁ
ꢂ
∂^{2 }
∂^{2}φ
∂y^{2 }
^{φ}_{2 }+
x
.
(20)
∂
The equation is balanced by a term that represents the physics of the particular
problem, and in general we call it a ‘forcing’ term f(x, y). In the above one
mg
dimensional example this forcing would be f(x) =_{ k} . The resulting equation is
known as ‘Poisson’s equation’,
ꢁ
ꢂ
∂^{2 }
∂^{2}φ
∂y^{2 }
^{φ}_{2 }+
x
= f(x, y).
(21)
2
∇
φ =
∂
You are required to solve an equation of this form for each physical problem,
but initially we will consider the special case of f(x, y) = 0, which is known as
‘
Laplace’s equation’. We already know how the second derivative in x is discre
tised:
ꢀ
ꢀ
ꢀ
ꢀ
∂^{2 }
x^{2 }
φ
φ_{i}_{+}_{1 }− 2φ_{i} + φ
≈
^{i}^{−}^{1 },
(22)
2
∂
(∆x)
i
and the y derivative follows suit:
ꢀ
ꢀ
ꢀ
ꢀ
∂^{2 }
y^{2 }
φ
φ_{j}_{+}_{1 }− 2φ_{j} + φ
≈
^{j}^{−}^{1 }.
(23)
2
∂
(∆y)
j
When we put them together, we need both i and j indices simultaneously to
clarify where exactly in space we’re retrieving values of φ. Our approximation to
Laplace’s equation becomes,
∂^{2 }
∂ φ
2
^{φ}i+1,j − 2φ_{i}_{,}_{j} + φ
^{φ}i,j+1 − 2φ_{i}_{,}_{j} + φ
^{φ}_{2 }+_{ ∂}_{y }
_{2 }= 0 ≈
^{i}^{−}^{1}^{,}^{j }+
^{i}^{,}^{j}^{−}^{1 }. (24)
2
∇
φ =_{ ∂ }
2
2
x
∆x
∆y
1
7
Laplace’s equation is a statement about the physics of a problem, not in itself
a solution technique. Written in this form we require to already know values of
φ(x, y) and we can then verify that the answer is zero, whereas we really want
to establish a numerical method for ﬁnding those values, based on some initial
guess. Following the suggested approach from the onedimensional example above,
we rearrange the physical statement into a procedure that we can implement
iteratively and believe will converge progressively from our initial guess to the
correct answer. Our ﬁrst step is to note that the spatially central value φ_{i}_{,}_{j }
appears more often than any other, it is said to be ‘dominant’ in the solution, so
we should perhaps treat it as special. We can rearrange our statement of physics
so that it is described by a function g of neighbouring values:
φ_{i}_{,}_{j }= g (φ_{i}_{−}_{1}_{,}_{j}, φ_{i}_{+}_{1}_{,}_{j}, φ_{i}_{,}_{j}_{−}_{1}, φ_{i}_{,}_{j}_{+}_{1}, ∆x, ∆y) .
(25)
and then repeatedly compute new values of φ_{i}_{,}_{j} until they converge. Multiplying
2
2
equation (24) by ∆x ∆y to get a common denominator, we have,
2
2
0
≈ ∆y (φ_{i}_{+}_{1}_{,}_{j} − 2φ_{i}_{,}_{j} + φ_{i}_{−}_{1}_{,}_{j}) + ∆x (φ_{i}_{,}_{j}_{+}_{1} − 2φ_{i}_{,}_{j} + φ_{i}_{,}_{j}_{−}_{1})
(26)
(27)
and teasing out the φ_{i}_{,}_{j} terms, we have,
2 2 2 2
∆y φ_{i}_{,}_{j} + 2∆x φ_{i}_{,}_{j} ≈ ∆y (φ_{i}_{+}_{1}_{,}_{j} + φ_{i}_{−}_{1}_{,}_{j}) + ∆x (φ_{i}_{,}_{j}_{+}_{1} + φ_{i}_{,}_{j}_{−}_{1}) .
2
We really want the equation in terms of φ_{i}_{,}_{j} alone, so we pull it out as a common
factor and divide both sides by the other factor, ie.
φ_{i}_{,}_{j }≈^{ ∆}y^{2} (φ_{i}_{+}_{1}_{,}_{j} + φ_{i}_{−}_{1}_{,}_{j}) + ∆x^{2} (φ_{i}_{,}_{j}_{+}_{1} + φ_{i}_{,}_{j}_{−}_{1})_{. }
(28)
2
∆x^{2 }+ 2∆y^{2 }
Of course this becomes a much neater expression if ∆y = ∆x, as on a regular
square grid, and in this case would reduce to calculating the average of four
numbers:
_{φ}_{i}_{,}_{j }_{≈} φ_{i}_{+}_{1}_{,}_{j }+ φ_{i}_{−}_{1}_{,}_{j} + φ_{i}_{,}_{j}_{+}_{1}^{ +}^{ φ}^{i}^{,}^{j}^{−}^{1} .
4
(29)
Coding this twodimensional update scheme equation in Matlab is a straight
forward extension of the script provided above. Further details are provided in
sections speciﬁc to each problem.
2
.5 Underlying principles
So what does this equation mean? It says that the centre value φ_{i}_{,}_{j} is the average
of the surrounding values in x ad y. It turns out that this averaging process
mirrors the behaviour of the physics of problems that contain second derivatives,
1
8
such as heat ﬂow, ﬂuid ﬂow, and stress distribution. It implies that there can be
no sudden jumps in the value of φ, since to raise one value egregiously with respect
to its neighbours would mean that it was not the average of its surroundings. So
how can we use it? Well, we’re aiming to develop an iterative scheme to walk
towards the correct answer from some arbitrary guess. Until we get to the correct
answer, or at least ‘close enough’, this equation probably won’t balance. The
amount by which it is unbalanced is called the ‘residual’ – the residue of error that
prevents the physical law from being satisﬁed. This equation forms the basis of a
technique called relaxation that works by taking an initial guess and repeatedly
smoothing, averaging and ‘relaxing’ values until the equation balances. So how
does this work? We take equation (29) and although it solves a steadystate
problem, we treat it in the same way we’d treat a timeintegration in some sort of
‘
pseudotime’ τ as we march towards convergence. We say the RHS is evaluated
at the ‘old’ time level n, say, and the LHS is the ‘new’ updated value at timelevel
n + 1, giving,
n
i+1,j
n
i−1,j
n
i,j+1
n
^{i}^{,}^{j}^{−}^{1 }.
φ
+ φ
+ φ
+ φ
n+1
i,j
φ
≈
(30)
4
Of course, when φ^{n}^{+}^{1} = φ^{n} then we’ve converged and therefore the equation
must be in balance and the steadystate correctly represented. Until convergence,
it describes a marching scheme, which can be written in the following form for
∆
τ = 1,
n
i+1,j
n
i−1,j
n
n
i,j−1
n
− 4φ
i,j
φ
+ φ
+ φ
+ φ
n+1
i,j
n
i,j
i,j+1
φ
≈ φ + ∆τ
.
(31)
4
Rearranging, it is easier to see that for an appropriate timescale τ (hidden in the
change of timescale from t to τ is an implied constitutive relation of some sort)
the update scheme is eﬀectively a discrete version of a time derivative.
n+1
i,j
n
i,j
n
n
i−1,j
n
n
i,j−1
− 4φ^{n}_{c }
φ
− φ
φ
+ φ
+ φ
+ φ
i+1,j
i,j+1
≈
(32)
(33)
∆
τ
4
∂
φ
τ
∂^{2}φ_{2 }_{+} ∂^{2}φ
=
.
∂y^{2 }
∂
∂x
So eﬀectively our relaxation technique obtains a steadystate solution by evolving
our arbitrary initial guess towards a steady state through a contrived sort of time.
This observation is often used to appropriately modify the update formula to
speed up convergence to a steady state. If we take larger steps through time,
then we could reach a steady state faster. Of course as we converge, then rates
of change of values reduce and by deﬁnition these time or pseudotime derivatives
tend towards zero. When they reach zero, we satisfy the rule that was obtained
from physics.
1
9
Problem Descriptions
3
.1 A. Thermodynamics
This problem concerns the calculation of temperatures inside a plate of size 0.8m×
.4m as shown in ﬁgure 5. Temperature distributions T(x, y) satisfy Laplace’s
equation in stationary equilibrium,
0
ꢁ
ꢂ
∂^{2 }
∂^{2}T
∂y^{2 }
^{T}_{2 }+
x
= 0.
(34)
2
∇
T =
∂
The plate is heated along parts of its top and left edges and cooled along parts
of its lower and right edges. The faces and remainder of the edges are insulated.
There is also a circular hole in the plate. The hole contains a highly conducting
ﬂuid, such that all points around the boundary of the hole have equal temperature.
Your Matlab program will meet the following minimum requirements for this
problem:
•
•
Uses a relaxation method to calculate the steadystate temperature distri
bution in the plate speciﬁed below.
The user should be able to specify x and y resolution (number of grid points),
a temperature convergence criterion, and be able to change the boundary
temperatures and the length of the boundary held hot or cold, and position
of the circular hole.
•
For the stagegate assessment and preparing your report, the code should
be able to provide clear graphical representation of a converged temperature
distribution T(x, y) and a convergence plot of ∆T vs. iteration count. This
should be turned oﬀ for the code performance assessment.
edg^{T}e ^{h}o^{e}f ^{p}th^{l}^{a}e^{t}^{e}h^{s}e^{d}at^{i}^{m}ed^{e}^{n}a^{s}n^{i}^{o}d^{n}c^{s}o^{a}o^{r}l^{e}ed^{0}^{.}^{8}se^{m}ct^{×}io^{0}n^{.}^{4}s^{m}ar^{.}e^{B}0^{y}.2^{d}m^{e}^{f}^{a}a^{u}n^{l}^{t}d^{t}0^{h}.^{e}1m^{l}^{e}^{n}a^{g}l^{t}o^{h}n^{a}g^{l}^{o}t^{n}h^{g}e^{ t}h^{h}o^{e}r^{v}iz^{e}o^{r}n^{t}^{i}t^{c}a^{a}l,^{l }
and maintained respectively at 100K and 0K. By default the circular hole has
radius 0.05m and is positioned at 0.3m to the right and 0.3m above the lower left
corner of the plate.
In locations where boundaries have a prescribed temperature, the value at
points on the outer edge is held ﬁxed. Boundaries that are described as ‘insulating’
indicate that there is no heat ﬂux through that part of the surface. Where there
is no heat ﬂux, there can be no temperature gradient driving one, so the edge
value and the nextinnermost value must be held to be identical (though with no
constraint on the speciﬁc value). On the edge of the circular hole, the values of
temperature must all hold a common value, though again there is no constraint
on the speciﬁc value, simply that they must hold the same value.
2
0
Figure 5: Boundary conditions for heated plate. Figure not to scale.
3
.2 B. Fluid dynamics
This problem concerns the calculation of steadystate incompressible ﬂow patterns
around obstacles in a twodimensional channel. Speciﬁcally, the purpose is to
estimate the directions of streamlines and local velocities in a 4m long section of
a 2m wide channel partially obstructed by several cylinders (circles), as depicted
below. In steady state, the velocities U and V can be described by a single variable
called a ‘stream function’ Ψ(x, y) according to the following relationships,
∂
∂
Ψ
x
=
=
−V
(
35)
∂
∂
Ψ
y
U.
Under appropriate conditions (detailed in appendix B), this stream function sat
isﬁes Laplace’s equation,
ꢁ
ꢂ
∂^{2 }
∂^{2}Ψ
∂y^{2 }
^{Ψ}_{2 }+
= 0.
(36)
2
∇
Ψ =
∂
x
In addition to the requirements for problem A, your Matlab program must also
satisfy the following minimum requirements for problem B:
•
Uses a relaxation method to calculate the steadystate velocity ﬁeld in the
2D channel speciﬁed below.
2
1
Figure 6: Typical arrangement of Streamlines around cylinders. Figure not to
scale.
•
•
The user should be able to specify x and y resolution (number of grid points),
a convergence criterion for stream function, and be able to change the inﬂow
velocity, and the position/radius of the cylinders.
For the stagegate assessment and preparing your report, the code should
be able to provide clear graphical representation of a converged velocity
ꢃ
ꢄ
U(x, y)
V (x, y)
a convergence plot of ∆Ψ vs. iteration count. This should be turned oﬀ for
the code performance assessment.
ﬁeld
, plotted as streamlines (contours of stream function), and
The quantities and dimensions required to construct the geometry are tabu
lated for your convenience:
Quantity
Symbol Dimension
Channel length
Channel width
Upstream ﬂow speed
Upstream ﬂow speed
(x)
(y)
(U_{∞})
(V_{∞})
4m
2m
1ms
0ms^{−}^{1 }
−
1
Cylinder locations and radii for the default problem are:
2
2
Cyl x_{c} (m) y_{c} (m) radius
1
2
3
4
5
6
0.5
1.0
1.5
1.9
2.5
2.3
0.7
1.4
0.5
1.1
0.4
1.6
0.2
0.3
0.3
0.2
0.2
0.2
On bounding walls and on cylinders in the interior, there is no penetration of
ﬂow through these boundaries. In such cases, the stream function Ψ (which is an
indirect measure of mass ﬂux per unit distance) takes a constant value. There
is no direct way of specifying what particular value that should be, just that all
values on the same boundary should have a common value. To provide a baseline
value, set Ψ = 0 in the bottom left corner. Thus all of the bottom wall will have
Ψ = 0. To be clear, the top wall and the cylinders will all have diﬀerent, but
individually constant values. For inﬂow we use the formula,
Ψ(y) = U_{∞}y
(37)
ie. a linear variation of Ψ with y whose gradient is the prescribed upstream ﬂow
speed. It follows that the stream function value on the wall at y = 2m has a
constant value Ψ = 2U_{∞}. At the outﬂow the value of Ψ must vary smoothly from
Ψ = 0 to Ψ = 2U_{∞}, however disturbances from the cylinders may still aﬀect the
distribution of velocity in the wake. To account for this, we do not enforce a
speciﬁc velocity proﬁle, but simply that there is no V component to the outﬂow
∂
Ψ
velocity. From equation (35) this implies_{ x} = 0 and we implement this by holding
∂
the edge value and the nextinnermost value identical (though with no constraint
on the speciﬁc value).
3
.3 C. Structural Mechanics
This problem concerns the displacement of a trampoline when subject to the
weight of a gymnast standing oﬀcentre as depicted in ﬁgure 7. Stretched mem
branes can only resist inplane stress, and do not support bending moment. There
is, unfortunately, a rip in the trampoline skin. Stress cannot be borne across the
rip; stresses can only be supported in the direction parallel to the rip.
In steady state, the displacement u satisﬁes Poisson’s equation,
ꢁ
ꢂ
∂^{2 }
∂^{2}u
∂y^{2 }
^{u}_{2 }+
x
= f(x, y),
(38)
2
∇
u =
∂
where the forcing is provided by the weight of the gymnast distributed over two
2
3
Figure 7: Skin of a trampoline with numerical grid superimposed, showing typical
gymnast contact patches, a typical rip in the membrane and the zerodisplacement
boundary around the edge. Figure not to scale.
small circles around their feet. In those two circles, the forcing is given by
f(x, y) =^{ m}_{A}_{k}^{g }
(39)
where m is the mass of the gymnast, g is gravity, A is the total area of contact
of the gymnast’s feet, and k is the stiﬀness of the elastic material. In addition to
the requirements for problems A and B, your Matlab program must also satisfy
the following minimum requirements for problem C:
•
•
Uses a relaxation method to calculate the steadystate displacement of the
trampoline skin speciﬁed below.
The user should be able to specify x and y resolution (number of grid points),
a convergence criterion for displacement, and be able to change the position
of the rip, and the position/weight of the gymnast’s feet, and the extensional
stiﬀness of the material.
•
For the stagegate assessment and preparing your report the code should be
able to provide clear graphical representation of a converged displacement
u(x, y) and a convergence plot of ∆u vs. iteration count. This should be
turned oﬀ for the code performance assessment.
The trampoline has a length of 4m and a width of 2m and is ﬁxed around
the edge. Taking the bottomleft corner as an origin, the default location of the
2
4
rip runs from (1.5m, 0.5m) to (1.5m, 0.8m). The default weight of the gymnast is
6kg, the contact area of their feet is approximated by two circles of radius 20cm.
5
By default the gymnast’s feet lie at (3m, 0.75m) and (3m, 1.25m). The stiﬀness
value you shoudl choose to scale your computation is largely determined by the
precise details of how you deﬁne it (whether by length, area, normalised strain etc).
for our purposes here it is not important how it is deﬁned, so choose an appropriate
value for the material’s extensional stiﬀness such that the maximum displacement
^{i}^{s }^{a}A^{r}^{o}l^{u}l^{n}o^{d}f t^{4}h^{5}e^{−}ed^{5}g^{0}e^{c}s^{m}o^{,}f^{t}t^{y}h^{p}e^{i}^{c}t^{a}r^{l}a^{o}m^{f}p^{t}o^{h}l^{e}in^{d}e^{e}^{ﬂ}ar^{e}e^{c}^{t}ﬁ^{i}^{o}x^{n}e^{s}d,^{f}^{o}s^{u}o^{n}t^{d}he^{o}y^{n} h^{a}a^{r}v^{e}e^{a}d^{l} i^{t}s^{r}p^{a}l^{m}ac^{p}e^{o}m^{l}^{i}e^{n}n^{e}t^{.} u = 0
everywhere on the edge. The rip is not able to support stress in the direction
across the rip (ie. in the xdirection), but it can support stress in the direction
parallel to the rip just as well as everywhere else. Since outofplane displacement
gradients are an elastic response to inplane stress, then there can be no gradient
of displacement in the xdirection anywhere adjacent to the rip – on either side.
Thus values on the very edge of the rip must be equal to those on the nextinner
most row of grid points, and the equivalent must be true on the opposite side of
the rip. Note that the displacements at both ends of the rip must be equal, so the
average diﬀerence in stress between the end points must be equal. In this case
we can make the convenient modelling approximation that the rip lies exactly
coincident with a mesh boundary, but the vertical displacements on either side of
the rip may take diﬀerent values.
4
Optimising your code
4
.1 The rationale
Now that you have passed the Week 5 stagegate you are expected to have pro
grams that solve each of the three problems, thermodynamics, ﬂuid dynamics,
and structural mechanics. You will have noticed that the time taken to converge
to a solution is dependent both on the problem you’re solving and how you write
the code that solves it. Optimisation is the process of carefully rewriting a work
ing code to make it faster. Good programming practice dictates that you ﬁrst
get a program working before you attempt to optimise it: only once you have
full oversight of what the code needs to do in order to work should you consider
taking ‘shortcuts’ that may (or may not) make the code more eﬃcient. The Week
5
stagegate has now given you this oversight and only now should you proceed
to optimise your code. Code performance on the sorts of physical problems you
study in this unit is not simply a matter of timing code on a particular case:
the aim is to have a general purpose tool that solves a range of related physical
problems to a userspeciﬁed accuracy in terms of convergence tolerance and grid
resolution.
2
5
4
.2 Eﬃciency and scaling
Timing information from one particular case can be misleading because there
are, broadly speaking, two components to the evaluation of a Matlab command:
parsing and execution. A command is parsed from a string of characters that
you ask Matlab to evaluate, and only then can it be executed as an operation
acting on data held in memory. The parsing step has a considerable overhead
cost, which is dependent on the length and complexity of the character string
that it receives. The execution cost is the bit that gets the job done, and as the
problem size increases, the cost of parsing doesn’t change but the relative fraction
of total time spent doing useful work on data held in memory increases. Thus at
very low grid resolution the cost of parsing dominates, but the overhead decreases
proportionally as grid resolution increases.
Matlab have tried very hard to reduce the overhead costs, and the vectorised
notation : for scanning over an array helps with this, since operations that would
naturally be expressed over several lines as a loop operation can be expressed in
a single line:
phinew(1,:)=0.0;
has the same eﬀect as
for iy=1:ny
phinew(1,iy)=0.0;
end
but requires fewer parsing operations to interpret. Internally the operation is
executed by a for loop in a compiled language (usually C or Fortran) in which
the parsing step is performed in advance rather than ’ontheﬂy’ and the for
operation is encoded as lowlevel binary instructions that execute much more
rapidly. Execution speed in Matlab has long been a drawback of the package,
especially with computationally intensive calculations.
The use of for and while shoudl be minimised in Matlab scripts. Careful
coding of array operations can yield equivalent outcomes to explicitly coded loops,
but are more eﬃciently executed. Consider the following code:
x = 0:dx:grid_length;
y = 0:dy:grid_height;
nx = length(x); ny = length(y);
[x_arr, y_arr] = meshgrid(x, y);
phi = zeros(ny, nx);
This generates arrays x arr and y arr having appropriate entries (note that a
sensible array ordering and standard matrix row/column ordering are not always
synonymous). For instance, this allows statements such as:
2
6
dirichelet1 = find((y_arr == grid_height) & (x_arr <= grid_length – bcportion));
dirichelet2 = find((y_arr >= grid_height – bcportion) & (x_arr == 0));
diricheletboundary = union(dirichelet1,dirichelet2);
where ‘Dirichelet’ is the type of boundary where the value is ﬁxed. These state
ments identify the grid points that can then be set to the desired value:
phi(diricheletboundary) = boundaryvalue;
Similar constructions can be used to identify other boundaries. Naturally there
are other methods to achieve the same results, but they are generally less easy
to interpret than the above. As an example, in the special case ∆x = ∆y, the
update step can be performed extremely eﬃciently:
inner_x = 2:nx1;
inner_y = 2:ny1;
x_diff = phi(inner_y, inner_x+1) + phi(inner_y, inner_x1);
y_diff = phi(inner_y+1, inner_x) + phi(inner_y1, inner_x);
old_inner = temp(inner_y, inner_x);
phi(inner_y, inner_x) = (x_diff + y_diff)/4;
though it should be noted that in the general case ∆x = ∆y and your code should
have the capability to handle the general case
There are other ways of improving the eﬃciency of your program; the gap
between your starting guess and the ﬁnal answer plays a big role in the time
taken, so how would you improve your starting guess? You could start by solving
a simpler relevant problem, and using the answer to this simpler problem as a
starting guess. If the target resolution is, say, an 80 × 40 grid, you could instead
start with a 10 × 5 grid, a much larger ∆x, ∆y, and a much shorter time to
converge on this approximate problem. Then interpolate this coarsely resolved
solution as a starting guess of a ﬁner mesh, maybe 20×10, solve this, then repeat
on 40 × 20, then ﬁnally on an 80 × 40. This method is a form of ‘multigrid’,
sonamed for obvious reasons, and it one of the fastest and most ﬂexible ways of
solving Poissonequation problems. The copying and interpolation process from
a coarse grid to a ﬁner one takes a little care to get the indexing correct, because
only every second value in each direction will have a directly supplied value. You
will need to perform a simple linear interpolation to ﬁll in the gaps. One hint
for designing your multigrid: reﬁnement by factors of 2 require interpolations
that are straightforward to implement, so it makes sense to choose the number of
gridpoints in x and y to be a power of two. You may well be asked to perform the
calculation on a grid with a nonpowerof2 on either side, so you should consider
carefully your strategy for such cases.
2
7
In all cases you must check that your modiﬁcations actually make your code
more eﬃcient… From Matlab 2015a onwards there is a new button ‘Run and
time’ on the top bar, which will give you performance statistics about your script.
In general multigrid will actually run less eﬃciently on very small problems, but
will get progressively more eﬃcient as the number of x and y grid points increases.
This is a favourable scaling property. The basic iterative method gets very slow
as the problem size increases, but this may be mitigated using a wellcrafted
multigrid method. You should plot the performance of your code at a variety of
resolutions. It is often instructive to consider the time taken to converge divided
by the total number of grid points as a measure of eﬃciency. If the resulting
graph were horizontal, then the time taken scales linearly with the problem size,
and this is the ideal case for an eﬃcient code.
The second way in which examining just one particular case can be misleading
is that time to convergence doesn’t scale linearly with grid resolution. You might
expect that the cost per grid point is a constant, independent of the size of the
problem, but in fact the cost per grid point increases. The reason for this is that
physical situations that give rise to Poisson or Laplace equations have one key
thing in common. A change of property at any location in space aﬀects every
other point, so a loose analogy might be the handshakes problem (10 people in
a room, 45 handshakes needed) every point is connected to every other and so
it takes more iterations before everything equilibrates and the current solution
balances the original equation to within an acceptable tolerance. Some of the
linear algebra behind this is detailed in the appendix, for those interested.
5
Final Submission
5
.1 Report submission guidance
Only one person from your group should submit a copy of the report. The only
acceptable format is .pdf (do not try to upload .zip ﬁles or any other format). The
ﬁlename of the uploaded .pdf should be group###.pdf where ### is your group
number expressed as three digits, and lowercase g is used. For example group 7
would submit a ﬁle group007.pdf. No other ﬁlenames will be accepted.
The only other stipulation on format is that the .pdf contains less than 10
pages, and note that there is no L^{A}T_{E}Xtemplate: you have a free format within
those 10 pages, but do take note of the marking rubric given in §1.6, which suggests
how you might organise your submission.
2
8
5
.2 Automated marking
The need to test how your code performs over a range of resolutions is imperative,
so we need to invoke your program multiple times for each problem at diﬀerent res
olutions and convergence tolerances. We will then plot the performance statistics
we need to evaluate correctness, error handling and scaling performance of your
code and assign a mark accordingly. It is unnecessary to do this datagathering
by hand, provided a common format is used to invoke your program, and the
solution is returned in a common format. The speciﬁcation in §5.3 details exactly
the arguments we will supply to invoke your code, and the solution we expect in
return. You will need to build a small wrapper around your code to accept our
arguments and return the data in the correct format. To ensure you can prepare
your code accordingly, a simpliﬁed assessment script is provided in §5.4, and af
ter the week 5 stage gate a representative exemplar submission solving the 2D
rope bridge problem will be issued, showing how to integrate your code with the
automated script.
5
.3 Code submission speciﬁcation
The code submission is 50% of the computational part of Modelling 2. You should
submit a single .m ﬁle containing all your code that you wish to be assessed.
The ﬁlename of the uploaded .m should be group###.m where ### is your group
number expressed as three digits, and lowercase g is used. For example group 7
would submit a ﬁle group007.m. Inside this ﬁle the name of the ﬁrst function
must correspond directly to the ﬁlename. It must also accept input and output
arguments in the following format, illustrated for group 7:
function [phiend,converge,st1,st2,st3,error]=…
group007(phistart,tol,prob,probargs);
%
put your code in here
end
The input arguments are speciﬁed as follows:
•
•
phistart is the array on which you are expected to return a solution. You
will need both the values and its its dimensions (number of elements in each
row and column). Note carefully that the grid will be interpreted as vertex
centred, so elements on the last row or column of the array represent values
that are on the boundary.
tol represents the convergence tolerance, and this is deﬁned to behave in
the same way as the template provided: it measures the gradient of error at
the beginning (comparing the original [phistart with the values after one
2
9
iteration) with respect to the gradient of error at the current iteration (com
paring current with immediately previous), and your iterative loop should
ﬁnish when the ratio of starting to current gradient is less than the speciﬁed
value.
•
•
prob is a single character speciﬁer for the problem you will solve and set
ting it to ’a’ should solve the thermodynamics problem, ’b’ should solve
the ﬂuid dynamics problem and ’c’ should solve the structural mechanics
problem. In the case shown below, ’d’ is used for the rope bridge exemplar.
probargs is a problemspeciﬁc parameter array, whose format is detailed
next.
The input argument probargs will in each case just be an array of numbers cod
ifying the setup of the problem. For the thermodynamics case ’a’, the argument
variables will determine only the location and size of the circle, and are speciﬁed
in the following order:
probargs(1)=x; %x coordinate of centre
probargs(2)=y; %y coordinate of centre
probargs(3)=r; %radius of circle
For the ﬂuid dynamics case ’b’, the argument variables will determine only the
location and sizes of the circles, and are speciﬁed in the following order:
probargs(1)=x1; %x coordinate of centre of circle 1
probargs(2)=y1; %y coordinate of centre of circle 1
probargs(3)=r1; %radius of circle 1
probargs(4)=x2; %x coordinate of centre of circle 2
probargs(5)=y2; %y coordinate of centre of circle 2
probargs(6)=r2; %radius of circle 2
.
..
and there may be more than 6 circles speciﬁed as input so you will have to check
the size of the incoming array to be sure how many there will be. For the structural
mechanics case ’c’, the argument variables will determine the position of the rip
and the position of the circles but not their loading (which should be preset along
with an appropriate stiﬀness for the membrane to achieve the speciﬁed maximum
displacement).
probargs(1)=x1; %x coordinate of centre of circle 1
probargs(2)=y1; %y coordinate of centre of circle 1
probargs(3)=r1; %radius of circle 1
3
0
probargs(4)=x2; %x coordinate of centre of circle 2
probargs(5)=y2; %y coordinate of centre of circle 2
probargs(6)=r2; %radius of circle 2
probargs(7)=xrip; %x coordinate of the rip
probargs(8)=ybot; %y coordinate of the bottom of the rip
probargs(9)=ytop; %y coordinate of the top of the rip
Note that the rip will always be oriented as per the description surrounding ﬁgure
7
. All dimensions and coordinates will be in the units of the original problem
speciﬁcation.
All information not provided by these arguments must be hard
coded into the function.
The output arguments are speciﬁed as follows:
•
phiend is the array which contains the solution. This must have exactly
as many elements in each row and column as phistart. Note carefully
that the grid will be interpreted as vertex centred, so elements on the
last row or column of the array represent values that are on the boundary.
You should comment in your report about any consequences for accuracy
that this vertexcentred representation introduces. Be careful to ensure that
numel(size(phiend) is identical to numel(size(phistart). The safest
approach is to set phiend=phistart and copy values into an unchanged
array.
•
•
converge is a two column data array supplying data for a 2D plot of con
vergence gradient ratio vs iteration count. The orientation of the array that
the calling script expects is given by converge=zeros(nrows,2);.
st1, st2 and st3 are character strings for the ﬁrst part of the email addresses
for each member of the group, and these should be set as follows:
st1=’jb007’;
st2=’cm16039’;
st3=’null’; %only 2 in this group
Failing to provide these in the format speciﬁed here will mean your code
will not contribute to your overall mark.
•
error should return 0 if your code has reached convergence and you are
conﬁdent that the solution is a good response to the input parameters.
Instead return 1 if your code recognises that there is an unsuitable input
speciﬁcation, or 2 if a failure to converge has occurred during execution.
3
1
5
.4 Calling Script
%
%
modelling 2 201920
simplifed script for automated assessement
clear all
close all
ngroups=61;
timeout=30;
nx=100;
ny=100;
%
%
default settings for problems A. B. and C.
thermodynamics 0.8m x 0.4m domain
proba=zeros(3,1);
proba(1,1)=0.3; %x coordinate of centre
proba(2,1)=0.3; %y coordinate of centre
proba(3,1)=0.05; %radius of circle
%
fluid dynamics 4m x 2m domain
probb=zeros(18,1);
probb(1,1)=0.5; %circle 1
probb(2,1)=0.7;
probb(3,1)=0.2;
probb(4,1)=1.0; %circle 2
probb(5,1)=1.4;
probb(6,1)=0.3;
probb(7,1)=1.5; %circle 3
probb(8,1)=0.5;
probb(9,1)=0.3;
probb(10,1)=1.9; %circle 4
probb(11,1)=1.1;
probb(12,1)=0.2;
probb(13,1)=2.5; %circle 5
probb(14,1)=0.4;
probb(15,1)=0.2;
probb(16,1)=2.3; %circle 6
probb(17,1)=1.6;
probb(18,1)=0.2;
%
structural mechanics 4m x 2m domain
3
2
probc=zeros(9,1);
probc(1,1)=3.0; %circle 1
probc(2,1)=0.9;
probc(3,1)=0.2;
probc(4,1)=3.0; %circle 2
probc(5,1)=1.1;
probc(6,1)=0.2;
probc(7,1)=1.5; %rip xcoordinate
probc(8,1)=0.5; %rip ystart
probc(9,1)=0.8; %rip yfinish
%rope bridge
probd=zeros(5,1);
probd(1)=10; %mass
probd(2)=9.8; %gravity
probd(3)=0.1; %stiff
probd(4)=1; %length
%
set up space for output statistics
durationA=zeros(ngroups,1);
durationB=zeros(ngroups,1);
durationC=zeros(ngroups,1);
durationD=zeros(ngroups,1);
users=zeros(ngroups,3,10);
crash=zeros(ngroups,3);
figure(1);
for ig=7:7 %1:ngroups
group=sprintf(’group%03d’,ig);
disp(group);
%test sequence
probnum=’d’;
probcur=probd;
tol=1e6;
phistart=zeros(nx,ny);
try
tstart=cputime;
%the line which calls your code is below:
3
3
[
phiend,converge,st1,st2,st3,error]=…
feval(group,phistart,tol,probnum,probcur);
durationD(ig,1)=cputimetstart;
users(ig,1,1:numel(st1))=st1(:);
users(ig,2,1:numel(st2))=st2(:);
users(ig,3,1:numel(st3))=st3(:);
if error==0 %expecting no error on this testcase
err=sum(sum(abs(phiendphistart)));
if err<1e4
disp(’test 1 failed to return a valid solution’);
crash(ig,1)=1; %flag unrecognised internal error
else
%
%
%
verify solution convergence by further iteration…
%
display outputs
disp(durationD(ig,1));
subplot(3,2,1);
plot(converge(:,1),converge(:,2));
subplot(3,2,2);
surf(phiend);
%consider recursively doubling resolution
if durationA(ig,1)<timeout
%
%
%
%
%
%
close all
clear all
phistart=zeros(nx*2,ny*2);
try…
catch…
end
end
end
else
disp(’test 1 failed to recognise incorrect input’);
crash(ig,1)=2;
end
catch
disp(’test 1 crashed internally’);
crash(ig,1)=3; %flag unrecoverable internal error
end
3
4
end
5
.5 RopeBridge/WeightedMembrane exemplar
function [phiend,converge,st1,st2,st3,error]=…
group007(phistart,tol,prob,probargs)
%prepare userids
st1=’jb007’;
st2=’cm16039’;
st3=’null’; %only 2 in this group
if prob==’d’
%obtain resolution and initial guess from phistart
sz=size(phistart);
nx=sz(1);
ny=sz(2);
phinew=phistart;
phiold=phistart;
phiend=phistart; %ensure output array has right size
%create some space for convergence plot
converge=zeros(30000,2);
%set constants
mass=probargs(1);
gravity=probargs(2);
stiff=probargs(3);
length=probargs(4);
dx=length/nx;
force=0.5*dx*dx*mass*gravity/stiff;
%set up while loop
errc=1;
erro=errc;
it=1;
while errc/erro>tol
%
enforce boundary displacement
phiold(1,:)=0;
phiold(:,1)=0;
3
5
phiold(nx,:)=0;
phiold(:,ny)=0;
%calculate updated displacements
for iy=2:ny1
for ix=2:nx1
phinew(ix,iy)=0.25*phiold(ix+1,iy)+0.25*phiold(ix1,iy)…
+
0.25*phiold(ix,iy1)+0.25*phiold(ix,iy+1)force;
end
end
%calculate error gradient
errc=0;
for iy=1:ny
for ix=1:nx
errc=errc+abs(phinew(ix,iy)phiold(ix,iy));
end
end
errc=errc/(nx*ny);
if it==1
erro=errc;
end
%
swap old and new for next iteration after calculating error
for iy=1:ny
for ix=1:nx
phiold(ix,iy)=phinew(ix,iy);
end
end
converge(it,1)=it;
converge(it,2)=errc;
it=it+1;
end
phiend(1:nx,1:ny)=phinew(1:nx,1:ny); %copy into output array
error=0; %rash assumption: don’t do this at home!
end
end
3
6
5
.6 Final checklist
Never trust user input. You have no control over what the automated script will
ask your code to do, you can’t be certain that you’ve interpreted the units in the
same way as the user and it is your responsibility to ensure that the code returns
gracefully. Don’t assume that because the automated script will have lots of code
to mark that it will not attempt to pass inconsistent or corrupted data. It will.
Matlab has two keywords, try and catch that enable us to test your code in a
safe ‘playpen’ and if it crashes the marking script will simply continue to the next
invocation attempt.
Ensure your usernames are provided within the ﬁle and returned as argument
•
•
Ensure your code returns a nonzero error if the combination of input
arguments don’t make sense, or your code can’t solve the problem posed by
these input arguments.
Ensure your code measures the initial rate of change of error from the phis
tart array exactly as given to you as an input argument: the input may
not be a zero initial condition. Measure the ratio of your current rate of
change of error to this initial rate of change of error to decide whether you
have met the convergence tolerance before returning phiend as an argument.
This calculation must be done at the resolution given by the phistart input
argument.
•
•
Ensure that the array size of your return argument phiend exactly matches
the array size of the input argument phistart, and that the requested
convergence tolerance tol as deﬁned above is met at this grid resolution,
not just on any other internal grids you may use.
Ensure that the physical problem (the boundary conditions, geometric di
mensions, positions of circles etc.) you solve is independent (excepting dis
cretisation errors) of the dimensions of the phistart array; there is no guar
antee that the array dimensions will match the aspect ratio of the physical
dimensions.
Ensure that your code does not require user intervention (GUIs, waiting for
•
Ensure that your ﬁnal submission does not produce any plots.
3
7
•
Ensure that your code only uses functions in the prescribed lists given in
5.7.
§
5
.7 Matlab function list
The following lists show several available functions and arithemtic/algebraic oper
ators that are commonly used in Matlab. Your program, when ﬁnally submitted,
is likely to make use of commands from the ﬁrst two lists: Programming, logic
and arithmetic and Common Functions. There is unlikely to be need to use
commands from the subsequent lists in your ﬁnal submission, however the plotting
commands will be especially useful for debugging and prepring plots to present
in your report. To obtain a good mark for your code, it is essential that you
remove any commands that will require user input, pauses execution
for any reason, or displays information to screen, because this will slow
your code down and possibly halt it altogether, resulting in a performance penalty
that you can easily avoid.
Programming, logic and arithmetic operators
:
+
+
/
.
&

<
>
>
<
=
Perform operation on an entire array row/column
Arithmetic on scalars
Algebraic operations on matrices
Gaussian elimination
Elementwise operations on arrays
Logical operator for and
Logical operator for or
Relational operator for less than
Relational operator for greater than
Relational operator for greater than or equal
Relational operator for less than or equal
Relational operator for equal
Relational operator for not
– * /
– *
* ./
=
=
=
f^{∼}or
Open a for loop
if, else, elseif, end Conditional statements
while Open a while loop
3
8
Common Functions
abs
Absolute value
ceil
cos
max
Round up to nearest integer values
Return cosine of input variable
Maximum value in a row/column of an array
Mean value of a given array
mean
min
Minimum value in a row/column of an array
round Round value(s)to nearest integer
sin
sign
size
Return sine of input variable
Returns sign of variable
Return size of an array
numel Number of elements in an array
pi
tan
3.141… to machine precision
Return tangent of input variable
zeros Array of zeros
ones Array of ones
Plotting commands
bar
Plot a bar chart
compass Compass plot
get
grid
List attributes of plots etc.
Show grid
gtext
hist
hold
Put text on the plot using mouse
Plot a histogram
Maintain current plot while others are plotted
Mesh surface plot
mesh
plot
2d plot function
plot3
set
3d plot
Change attributes of plots etc.
Stairstep plot
stairs
stem
surf
Plot a stem plot
3d coloured surface plot
Put text on the plot
Plot title
Label for the x axis
Label for the y axis
text
title
xlabel
ylabel
3
9
Matrix algebra, complex number and statistical commands
inv
Matrix inverse
pinv
tril
triu
det
Matrix pseudoinverse
Lower triangular part of a matrix
Upper triangular part of a matrix
Determinant of a matrix
Identity matrix
eye
mldivide Left matrix divide
conj
real
imag
angle
rand
randn
cov
Complex conjugate
Real part of complex number
Imaginary part of complex number
Phase angle or complex argument
Generate a random number
Normally distributed random numbers
Covariance
std
Standard deviation of a given array
4
0
Miscellaneous commands
cd
Change directory
clear
clock
date
Clear all variables
Current date and time as date vector
Display current date
diary
dir, ls
edit
Start recording a diary
List all ﬁles in the directory
Start the editor
eps
format
help
Floating point relative accuracy
Set output format for displaying numbers
Display Matlab help related to command
Display the help window
Load ﬁle into Matlab workspace
Look for commands or .m ﬁles
Get/set search path
helpwin
load
lookfor
path
quit
Exit from Matlab
save
who
whos
Save workspace variables to disk
List current variables
List current variables showing type and size
Prompt user for input
Wait for user response of length of time
Display input variable
input
pause
disp
clc
fprintf
num2str
Clear command window
Write data output to a ﬁle
Convert a number to a string
meshgrid Create an indexed array
A Physics of heat conduction
Fourier’s Law of thermal conduction says that the local heat ﬂux density, q, is
equal to the product of material’s thermal conductivity, k, and the negative local
temperature gradient, −∇T. The heat ﬂux density is the amount of energy that
ﬂows through a unit area per unit time. The reason for the negative sign is that
the heat ﬂows away from hot things towards cold things, ‘downhill’ and therefore
against the temperature gradient as its vector direction is deﬁned. Thus the heat
ﬂux is given by,
q = −k∇T,
(40)
or in one dimension, simply,
q = −k^{∂}^{T}_{x} .
∂
(41)
4
1
To describe how temperature (or any other quantity) diﬀuses in time through a
medium, we invoke Fick’s 2nd Law of diﬀusion, which in one dimension says that
^{∂}_{∂}^{T}_{t }=_{ ∂}^{∂}_{x}^{q }
.
(42)
To see why this makes physical sense, we integrate the LHS in space over an arbi
R
L
∂
T
trary length 2L.
dx can’t be solved immediately, but we recognise that in
−
L ∂t
most physical situations we can swap the order of integration and diﬀerentiation,
R
L
∂
t
and equivalently write_{ ∂ }
Tdx, which is now much easier to interpret as the
−L
time rate of change of the total ’amount’ of T in the length 2L. Doing the same
integration on the RHS, we can say that
R
L
∂q
L ∂x
L
−L
dx = q = q(L) − q(−L) is the
−
diﬀerence in the amount of q at both extremes of length 2L at this timeinstant.
Remember, q is a ﬂux of heat, and so heat could be ﬂowing very quickly through
the length 2L without the total amount contained in the length 2L changing at
all. The temperature T will only increase in time if there is more heat ﬂowing in
through, say, the −L boundary than, say, is ﬂowing out through the L boundary
of our arbitrary length.
Substituting the heat ﬂux deﬁnition (41) into equation (42), we realise we have
a mixed partial diﬀerential equation in one dimension,
^{ꢁ}−k
ꢂ
.
(43)
^{∂}_{∂}^{T}_{t }=
∂
∂T
∂x
∂x
Often the thermal conductivity k is not a spatial variable, so it can be pulled
∂
outside the
operator. Also, we often – as in this problem – consider only
∂T
∂
∂
x
steadystate behaviour, where the time derivatives become negligible:_{ t} → 0. In
such cases, we have the much simpler onedimensional form
∂^{2 }
^{0 }= −k_{ ∂ }
T
x^{2 }
(44)
In the more general case of two dimensions, this becomes,
= −k^{ ꢁ}^{ T}_{2} +
ꢂ
,
(45)
∂^{2 }
∂^{2}T
∂y^{2 }
0
∂
x
which appears in a wide range of physical contexts. Here, it describes the varia
tion in temperature T in a plate with no internal generation of heat. Since the
conductivity k is normally constant and nonzero, then the second term must be
zero. The vector form of this simpler equation is generally known as Laplace’s
equation,
∂^{2 }
∂ T
2
^{T}_{2 }+_{ ∂}_{y }
_{2 }= 0.
(46)
2
∇
T =_{ ∂ }
x
4
2
B Stream function descriptions of incompress
ible ﬂow
Incompressible ﬂuid is a ﬂuid whose volume is a conserved quantity, you cannot
squeeze the ﬂuid into a smaller space (think ‘water’!). This means that in one
dimension (think ‘pipe’!) you can measure the ﬂux of volume (same as the mass
ﬂux for a density of 1) at one end of an (inelastic) pipe, and it must exactly match
the ﬂux of volume at the other. Flux of volume Q and average ﬂow rate U^{¯} are
connected, Q = U^{¯}A, where A is the crosssectional area of the pipe. If our pipe
has length 2L between −L and L, then we can say that Q(L) − Q(−L) = 0. Now
R
L
∂Q
dx. So what we’re really saying
think about the solution to this integral:
−
L ∂x
by an incompressible onedimensional ﬂow is that,
Z
L
^{∂}^{Q}_{x }dx = 0.
(47)
_{−}L ^{∂ }
In two dimensions it is only a little more complicated, but instead of thinking
of a pipe, we think of a square box where ﬂow can only pass through its faces
at rightangles to them. Again incompressibility means that the total amount of
‘
stuﬀ’ in the box is constant, but for every new bit that ﬂows in, an old bit must
ﬂow out. This time, however, there is no constraint on the direction that ﬂuid can
ﬂow in or ﬂow out, just that the total volume mustn’t change. We can extend the
integral description above to two dimensions by noting that we need to deﬁne a
vector version of Q. We could use U and V for the x and y velocities respectively,
and we should clarify whether we are talking about face areas in the vertical or
the horizontal, so Q = U^{¯}A_{x} + V^{¯ }A_{y}. In this context we can choose the shape of
our box, so let’s make it square, A_{x} = A_{y} = A = const., and with sides 2L as
before. Thus we could write,
Z
Z
L
L
^{∂}^{Q}_{x }dx +
^{∂}_{∂}^{Q}_{y }dy = 0.
(48)
_{−}L ^{∂ }
−L
There’s not much we can do with this equation until we realise that with the ﬂow
∂
∂
Q
∂Q
we’ve deﬁned through our box,_{ x} doesn’t vary in y and_{ ∂}_{y} doesn’t vary in x, so
R
L
∂Q
L ∂x
adding terms like
dy, simply adds zero to the equation. However it helps
us group the terms meaningfully,
−
Z
ꢁ
ꢂ
Z
_{ꢁ}_{∂}_{∂}_{Q}_{x }_{+}_{ ∂}_{∂}_{Q}_{y} ꢂ
L
L
∂
^{Q}_{x }+^{ ∂}_{∂}^{Q}_{y} dx +
dy = 0.
(49)
∂
−
L
−L
Following from the fundamental theorem of calculus, if we let our box tend to zero
size without changing its shape, then we could disregard the integrals altogether,
4
3
Figure 8: Illustrative plot of stream function and corresponding streamlines in a
uniform ﬂow from the left.
and write down a diﬀerential form,
^{∂}_{∂}^{Q}_{x }+^{ ∂}_{∂}^{Q}_{y} = 0.
(50)
∂
∂
Q
∂Q
∂y
From the way we’ve deﬁned the ﬂow in our box,_{ x} = AU^{¯} and
= AV^{¯}. The
area A is a common factor which tends towards but never reaches zero, so we
could say that in two dimensions we satisfy the following equation,
∂
^{U}_{x }+^{ ∂}^{V }
∂y ^{=}^{ 0}^{, }
(51)
∂
which is called the ‘Continuity Equation’ for incompressible ﬂow.
We need to introduce another concept – that of the stream function. Steady
incompressible ﬂow in two dimensions has a special property that its velocity,
a vector ﬁeld quantity, can be described completely by a separate scalar ﬁeld, a
single number whose value varies throughout the domain. Returning to our square
box, this single number is a relative measure of the volume ﬂux through a surface
–
or a bit of a surface. It needs a datum somewhere, usually a bottom lefthand
corner is set to zero. For every unit of distance up the left side of our box from
the bottom corner to a variable point Y , it contributes a small amount δΨ to the
stream function Ψ:
Z
Y
δΨ(Y ) =
Udy
L
(52)
−
If we let Y = L, then we recover the volume ﬂux through the whole left side of the
R
L
box, _{L} Udy = U^{¯}A_{x} as introduced earlier. The accumulation of stream function
−
for a uniform ﬂow is illustrated in ﬁgure 8. Another contribution comes if we now
4
4
move from top left to top right, and to establish a consistent sign convention, we
say postive volume ﬂux is V > 0 ﬂow out the positive yaxis direction (just as
postive volume ﬂux is U > 0 ﬂow in the postive xaxis direction). So we have,
Z
Z
L
L
Ψ =_{ L} Udy −
V dx.
(53)
−
−L
One of the interesting features of the stream function is that provided the ﬂow is
incompressible, it doesn’t matter how one moves from top left to bottom right,
it can be any curve so long as we do this without forming a loop or a knot; the
volume ﬂux is still how much ‘stuﬀ’ moves from one side of the line to the other
and the diﬀerence in values of Ψ at bottom left and top right reﬂect this. These
values are not dependent on the curve drawn between them.
Letting our box once more tend to zero size without changing shape, we obtain
the diﬀerential form
δΨ = Uδy − V δx.
(54)
If we diﬀerentiate the stream function, then since xdiﬀerentiation has no eﬀect
on y and viceversa, we have,
∂
Ψ
=
=
−V
(55)
(56)
∂
∂
x
Ψ
U,
∂
y
thus a single ﬁeld Ψ can describe both components of velocity. But what if
we don’t like the sign convention? We could invent another single ﬁeld that is
ꢃ
ꢄ
U
V
everywhere oriented at 90 degrees to Ψ: the vector
is always perpendicular
ꢃ
ꢄ
−
V
U
to
. Let’s call our new ﬁeld φ, deﬁned as follows:
∂
φ
=
=
U
(57)
(58)
∂
∂
x
φ
V,
∂
y
So what constraints are there on φ? We need to ensure incompressibility, so we
substitute expressions for U and V into the continuity equation (51),
∂
U
+ ^{∂}^{V }
= 0
= 0
= 0.
(59)
(60)
(61)
^{ꢂ }−^{ ∂}^{x }ꢁ ∂y
ꢁ
ꢂ
∂
∂φ
∂x
∂
∂φ
∂
x
∂y ∂y
∂^{2 }
φ
∂^{2}φ
_{2 }+_{ ∂}_{y}_{2 }
x
∂
4
The vector form of this simpler equation is generally known as Laplace’s equation,
2
∇
φ = 0,
(62)
that crops up very often in all sorts of physical problems. Here, φ is in fact the
symbol traditionally used for this quanitity, and is known as the velocity potential.
We can show that the stream function Ψ also satisﬁes incompressibility, by
substituting as follows, and recognising that for smooth ﬁelds the order of the
derivatives in x and y can be swapped:
ꢁ
^{ꢂ }+_{ ∂}^{∂}_{y }
^{ꢁ}−
ꢂ
^{∂}_{∂}^{U}_{x }+^{ ∂}_{∂}^{V}_{y }
∂
∂Ψ
∂Ψ
= 0
(63)
^{= }∂x ∂y
∂x
If we substitute U and V for Ψ, then,
ꢁ
_{∂}_{x }^{ꢂ }+_{ ∂}^{∂}_{y}^{ ꢁ }
ꢂ
∂
Ψ =_{ ∂ }
x
∂Ψ
∂Ψ
∂y ^{=} ∂^{∂}x (−V ) +_{ ∂}^{∂}_{y }^{(}^{U}^{)}^{ , }
2
∇
(64)
and by comparing with the vorticity vector in two dimensions, we have the fol
lowing:
so we can reduce this to the following simple equation
2
∇
Ψ = −~ω
(66)
For the condition ∇^{2}Ψ = 0 (the equation we are actually solving) to be true,
we require that ω = 0 everywhere as an additional constraint. This is called
the condition of ‘irrotationality’ and only holds exactly true for straight ﬂow in
free space. However we tell ourselves that this holds approximately everywhere
because it makes obtaining solutions a bit easier, and we cross our ﬁngers and
2
hope that the error isn’t too large. In this problem we will aim to solve ∇ Ψ = 0,
but be aware that the ﬂows it produces are not always physically meaningful. You
should identify in your report where you believe the ﬂows are representative, and
where they are not.
C Derivation of thin membrane mechanics
To a good approximation, thin elastic membranes can be modelled as supporting
tensile stress but have no resistance to compression and consequently no nega
tive strain, and also cannot resist bending moment or shear. If any one part of
4
6
Figure 9: Discrete analog of an elastic string.
the elastic material locally undergoes only small displacements as a result of in
plane tensile stress, then a linear elastic model is appropriate, even if globally,
displacements from the unstressed state are large. If initially we consider a one
dimensional elastic material (think ‘elastic band’!) as a slightly simpler analogue,
then we can derive an equation of motion relating displacements u to axial forces
F, and material stiﬀness k. To make it easier to conceptualise, we pretend that
the mass of the elastic is lumped at regular intervals, and in between are exten
sional springs, as depicted in ﬁgure 9. We derive the equation of motion for the
mass at a general position x + h. Newton’s second law states that,
ΣF = mu¨,
(67)
and in the absence of local applied forcing, the force applied to the mass is due to
tensile stress in the elastic springs, ie.,
mu¨ = ΣF = F_{r}_{i}_{g}_{h}_{t} + F_{l}_{e}_{f}_{t }
(68)
(69)
(70)
=
=
k [u(x + 2h) − u(x + h)] + k [u(x) − u(x + h)]
k [u(x + 2h) − 2u(x + h) + u(x)] ,
noting that extension of the right spring leads to a rightward force on the mass at
x+h, but extension of the left spring leads to a leftward (negative) force. Now we
need to account for contributions from all the point masses, and the cumulative
stiﬀness and length of all the springs put together. Thus we say the total length L
is given by L = Nh, where N is the total number of point masses. The total mass
4
7
M = Nm, and the total stiﬀness K = k/N (because as we introduce more and
more springs, any displacement is distributed amongst the springs, so the ratio
of applied forcing to strain response, ie. the stiﬀness, is reduced). We know that
the local dynamical system for the mass at x + h will oscillate in simple harmonic
motion with a natural frequency given by,
ω =_{ m}^{k }
,
(71)
2
but we need to derive a new frequency to describe the whole system. By rear
ranging the expressions for M, L and K, we have,
ꢁ
ꢂ
2
_{m}^{k }=
(NK)
M
= N_{ M}^{K} =
L
h
K
.
(72)
2
M
N
Thus by dividing equation 70 by m we can write,
_{u}_{¨ }_{= }_{L}_{2} K [u(x + 2h) − 2u(x + h) + u(x)]
(73)
2
M
h
and following the fundamental theorem of calculus we can let h → 0. We know
that the numerator takes the form of a second derivative of u evaluated at location
q
K
M
x + h, so if c = L
we can say,
∂^{2 }
_{2 }∂^{2}u
∂x^{2 }
^{u}_{2 }=
t
c
,
(74)
∂
for onedimensional inplane elastic motion with small strain. The quantity c has
q
K
M
dimensions of velocity and is interpreted as a wave speed, with
a timescale
that represents the natural frequency of the system as a whole. This equation is
also valid for displacements in the direction normal to the plane that result in small
angular deviation of the elastic (eg. a drum skin, or in this case a trampoline),
since the local behaviour is still dominated by small strain. In the more general
case of two dimensions, this becomes
ꢁ
ꢂ
∂^{2 }
∂^{2}u_{2 }_{+} ∂^{2}u
^{u}_{2 }=
c
,
(75)
2
∂
t
∂x
∂y^{2 }
a form of wave equation that appears in a wide range of physical contexts. Wave
equations of this sort describe the vibration of membranes such as drum skins,
or in this case, of a trampoline. In this problem we will be interested in the
2
∂
∂
u
equilibrium stress on the torn trampoline, when
= 0 in a steady state. In
t^{2 }
this case, then since c is not normally zero, equation (75) reduces to Laplace’s
equation,
∂^{2 }
∂ u
2
^{u}_{2 }+_{ ∂}_{y }
_{2 }= 0.
(76)
2
∇
u =_{ ∂ }
x
4
8
D Numerical convergence
To understand why larger grid sizes lead to slower convergence, we need to express
the numerical method as operations in linear algebra. Taking the 1D equation for
a rope bridge for simplicity,
∂^{2 }
mg_{, }
^{φ}_{2 }=
x
(77)
∂
k
then we treat the unknown values of φ on the grid as a vector φ, one element
for each grid point. We also treat the loading as a vec_{2}tor of the same length,
∂
x
where in this case every element is equal. The operator_{ ∂}_{ 2 }can be discretised and
represented as a square matrix with the same number of rows and columns as the
vector has elements. Using our normal discretisation,
∂^{2 }
φ
_{2 }_{≈} φ_{i}_{+}_{1} − 2φ_{i} + φ_{i}_{−}_{1 }
(78)
∂
x
(∆x)^{2 }
we can separate out the coeﬃcients in front of φ from φ itself, and put those
coeﬃcients in the matrix that is otherwise zero. Eﬀectively we are saying that the
action of performing a doublederivative on any vector of values can be written
down as a matrix in the following way:
Thus, the linear algebra problem that is approximately equivalent to the physical
statement for the rope bridge is given by,
which we might also write in a more compact form as Aφ = f. Direct inversion,
φ = A f is in principle possible (provided we check that A has no zero eigen
values and is therefore invertible), but using a direct method to compute A ,
then for every N grid points it takes O(N ) arithmetic operations to compute the
−
1
−
1
3
exact solution. Faster methods tend to use iteration to convergence and are less
4
9
exact, but the very best methods can reduce the cost to O(N), so the beneﬁt is
worthwhile on larger problems relative to the loss of accuracy. The pathway to
these fast methods is described in §4.2. Below we discuss the principle behind all
iterative approaches to solving this (and many other) problems.
The iteration process can also be interpreted as a matrix operation on vectors.
This time, we need to think of the matrix as a tool for moving a vector from point
ing in one particular direction towards another, which is what matrices do when
operating on vectors. We start with an initial guess, and move this initial guess
around in space until it converges. The elementwise operation being performed
during an iteration of the rope bridge case is the following,
old
i+1
old
φ
+ φ
2 ^{m}^{g }
^{i}^{−}^{1 }− (∆x)
_{2}_{k },
new
φ_{i }
=
(81)
2
and expressed for all elements we have,
which can be thought of as a sequential scaling and rotation of the vector φ^{ˆ} (the
modiﬁed vector composed of φ with a 1 added on the end) by the matrix that
we might call B. We need the extra column on the matrix and the extra element
mg
2
on the vector to account for the loading −_{ 2}_{k} (∆x) . Repeatedly applying B to a
vector has the following structure,
1
0
φ = Bφ^{ˆ }
ˆ
2
1
φ = Bφ^{ˆ }
ˆ
3
2
φ = Bφ^{ˆ }
ˆ
(
83)
4
3
φ = Bφ^{ˆ }
ˆ
5
4
φ = Bφ^{ˆ }
ˆ
.
.
.
k+1
k
The general case is φ^{ˆ }
= Bφ^{ˆ} , and when convergence is reached then the
≈ φ^{ˆ} . This
k+1
k
solution changes by an imperceptible amount, so by deﬁnition φ^{ˆ }
5
0
is just a special case of an eigenvalue problem where the eigenvalue happens to
be 1:
1
.φ^{ˆ }= Bφ^{ˆ}.
(84)
so the converged solution φ^{ˆ}^{∞} is simply an eigenvector associated with the matrix
B constructed from the numerical method. In fact the only way to compute
eigenvalues and eigenvectors in large matrices (with 6 or greater rows and columns)
is to iteratively change a vector from a starting location to convergence. The only
additional trick is that when the eigenvalue is not 1, the vector length needs to
be rescaled, and the amount by which it needs to be rescaled is the eigenvalue
associated with that eigenvector. This is called a power iteration. The choice of
initial vector is purely arbitrary, provided no eigenvalues of B are larger than 1
then the numerical method will converge to the principal eigenvector, and this is
∞
the solution φ^{ˆ} of the problem.
What speeds up or slows down a numerical method is how quickly the principal
eigenvector can be isolated from the others, and this is related to the relation
ship between its eigenvalue and all the other eigenvalues. The largest eigenvalue
tends to dominate, so when the ratio of λ_{m}_{a}_{x }
2
the same tolerance grows as O(N ) in the number of points, and not O(N) as one
might have originally expected.
5
1