MENG21712
10 Credits
Modelling with computers (TB1)
Dr. A.G.W. Lawrie
B Stream function descriptions of incompressible flow
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 find 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
follow-on 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 Computer-Based 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 problem-solving. The TB1
component of this course (MENG21712) will teach the principles of solving real-
world physical problems whose structure is governed by partial differential 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 confidence
will grow.
On this course we aim to provide a supportive environment for the development
of programming skills. To manage your workload, a stage-gate in week 5 will
ensure that you have produced code that meets the pass-standard of the course
by the half-way point. During the final 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 final 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 21-point scale that matches marks to
descriptors of achievement, and first-class (¿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 offering explicit teaching. As a further incentive, a prize is awarded
for the best-performing 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 2-hour slot in one of these rooms.
In the laboratory there will be teaching assistants who can assist you wherever
you have difficulty 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 no-cost
off-campus access to Matlab in a variety of ways: firstly 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 different 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 final 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, confidential
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 stage-gate that must be reached by your laboratory session on
Monday 28th October 2019, and this ensures that your group’s code meets
the pass-standard for this element of the course well in advance of the final 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 final solution and a splot showing the path
taken to converge to it. If satisfactory, the teaching assistant will sign each in-
dividual’s stage-gate assessment form and assign your group a unique identifier
number. It is your responsiblity to return the completed form to the unit director.
For the final deadline, 5pm, Friday 6th December 2019, your submission
must include a code (a single Matlab .m file containing ASCII text only) and a
report (a single .pdf file of no more than 10 pages produced using the LATEX doc-
ument preparation system), submitted electronically on Blackboard. Both will
pass through a sequence of manual and electronic honesty-checks to guard against
plagiarism. We may also perform random spot-checks, and request constituent
components of your report, eg. figures, .tex .bib files and so forth. Your report
will be marked conventionally against the 21-point scale, and your code will be
tested exhaustively in an automated script to assess its performance on a range
of test-cases. 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 stage-gate, 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 time-consuming to develop, and in this course the focus is on
algorithmic efficiency 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 offers 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 file.
5
1
.6 Rubric and intended learning outcomes
At the end of TB1 students should be able to:
•
•
represent real-world physical problem in mathematical notation using partial
derivatives
transform partial differential 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 specification
graphically present results obtained by using your code
use the LATEX document preparation system
summarise their work in a written report
A first-class student will also be able to:
•
•
•
•
•
optimise the implementation of an algorithm for efficiency
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 test-cases
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 21-point scale for
coursework, and the relevant interpretation of the scale for this task is tabled
beside the descriptor in figure 1.
There is no LATEXtemplate 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 LATEX
figures with correctly labelled axes
figures 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 final mark for TB1 has three contributions, 10% from
the week 5 stage-gate, 50% from the final 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
18-20: Report reaches publication standard in
terms of content for a reputable academic journal,
with little other than format adjustment to fit the
style requirements of the journal.
15-17: 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 figures
unambiguous at a first reading.
12-14: Clear evidence in the report of successful
execution of the required tasks, expressed with
written clarity and good support from diagrams and
figures, but lacking in originality when seeking
solutions to problems.
9-11: 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 figures, and there is little or no evidence
of originality or ingenuity.
6-8: 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 insufficient thought on the problems at
hand.
0-5: Partial, invalid or non-submission of report,
clear evidence of dysfunctional execution of the
project and/or preparation of the report.
Figure 1: Interpreting the University 21-point 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 LATEX 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
different areas of engineering: Thermo-dynamics, 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 differences between them
are simply in the physical context, the underyling partial differential 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 computer-modelling of
the problems posed. We demonstrate in §2.1 how to take an equation derived
from physical principles and solve it on a simple one-dimensional 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 specific remit for each problem is discussed in detail. In
section §4 there are suggestions for ways to optimise your code for greater effi-
ciency, and ways to maintain efficiency as the problem resolution grows. Finally
§
5 contains important information about the final submission specification and
format. In the appendices, we provide background information on how several
different physical principles result in the same form of partial differential equa-
tion, and additional technical details on why numerical methods get less efficient
as the problem size grows. The information in the appendices will inform your
report writing and finalisation 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 first step to make progress towards a solution is to copy the code for the
rope bridge (demonstrated in lectures) into a text file 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 specfic 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, specific 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
centuries-old engineering solution to crossing a gorge, long before our local hero,
Isambard Kingdom Brunel, designed the first 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∂∂
x2 = mg,
(1)
where m is the mass per unit length, g is gravity and x is the ordinate in the
span-wise direction.
2
.1 Analytical solution
Equation (1) is a differential equation we can solve directly by integrating twice.
The first integration gives,
∂∂φx =
mg
k x + C1,
(2)
where C1 is a constant of integration. Integrating again we have,
mg
φ(x) = 12 k x + C1x + C2.
2
(3)
We can set the constants C1 and C2 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
= C2
mg
k
2
x = 0, φ = 0 : 0 =
0 + 0C1 + C2
(4)
(5)
0
Similarly, and using the fact that C2 = 0,
x = L, φ = 0 : 0 =
1
2
mg
2
L + C1L + 0
(6)
(7)
k
C1
= −1
mkgL.
2
1
1
So the deformation can be described as,
1
2
mg
k
1 mgLx + 0
2 k
2
φ(x) =
x −
(8)
(9)
mg
=
2k 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 difficult to find 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 ‘finite
differences’, essentially using the original definition on which Calculus is based,
f0(x) = f(x + h) − f(x)
,
(10)
h
but instead of allowing h → 0 we let h → ∆x, some small non-zero 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 half-way between i and i + 1 is given by,
ꢀ
ꢀ
ꢀ
ꢀ
∂
∂
φ
x
≈ φi+1 − φi .
∆x
(11)
i+ 12
Often, as in this case, the rules of engineering problems involve second derivatives,
the ‘derivative of the derivative’. Thus,
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
∂φ
∂φ
∂x i− 2
ꢀ
∂2
x2
φ
∂x i+ 12
−
∆x
1
≈
.
(12)
∂
i
ꢀ
ꢀ
ꢀ
∂
∂
φ
∂φ
, but we still need
ꢀ
.
We already have the value for
x i+ 12
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
(xiï1,qiï1)
(xi+1,qi+1)
(xi ,qi)
x
6
x
1
2
1
2
xiï1
xiï ï
xi
xi+ ï
x
i+1
Figure 3: A typical out-of-plane 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)
∂
x2
∆x
i
Simplifying, we have,
ꢀ
ꢀ
ꢀ
ꢀ
∂2
x2
φ
(φ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 left-hand-side then the expression must balance the right-hand-side.
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 first 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 first 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 defines 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 stiffer 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 difference between φ
and φ and so the update indices
and n+1 are interchangeable and therefore can be safely ignored. We take this
to be the definition 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 finite patience and tend to trun-
cate our marching index when the degree of imbalance is within some acceptable
(
user-defined and problem-dependent) 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 user-specified parameters for grid-spacing, 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 specified 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 step-by-step 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=1e-6;
it=0;
while errc/erro>tol
%enforce boundary displacement
phiold(1)=0;
phiold(nx)=0;
%calculate updated displacements
for ix=2:nx-1
phinew(ix)=0.5*phiold(ix+1)+0.5*phiold(ix-1)-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:length-dx,phinew);
Note in particular how the convergence tolerance is enforced. The tolerance
tol is a user-defined constant, which in all the exercises is set to 1e-6. 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 difference 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
final ‘converged’ value is measured with respect to the value at the end of the first
iteration, given by the ratio errc/erro. We are thus defining 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 fishing 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 figure 4 In many cases the
physics established for a one dimensional problem can be very simply translated
1
6
qi,j+1
qiï1,j qi,j
qi+1,j
qi,jï1
Figure 4: Stencil on a rectangular mesh.
2
∂
φ
into two dimensions. Wherever derivatives ∂ 2 appear, these translate to
x
ꢁ
ꢂ
∂2
∂2φ
∂y2
φ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φ
∂y2
φ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
x2
φ
φi+1 − 2φi + φ
≈
i−1 ,
(22)
2
∂
(∆x)
i
and the y derivative follows suit:
ꢀ
ꢀ
ꢀ
ꢀ
∂2
y2
φ
φ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 finding those values, based on some initial
guess. Following the suggested approach from the one-dimensional 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 first 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 ≈ ∆y2 (φi+1,j + φi−1,j) + ∆x2 (φi,j+1 + φi,j−1).
(28)
2
∆x2 + 2∆y2
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 two-dimensional update scheme equation in Matlab is a straight-
forward extension of the script provided above. Further details are provided in
sections specific 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 flow, fluid flow, 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 satisfied. 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 steady-state
problem, we treat it in the same way we’d treat a time-integration in some sort of
‘
pseudo-time’ τ 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 time-level
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 steady-state 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 time-scale τ (hidden in the
change of time-scale from t to τ is an implied constitutive relation of some sort)
the update scheme is effectively 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φnc
φ
− φ
φ
+ φ
+ φ
+ φ
i+1,j
i,j+1
≈
(32)
(33)
∆
τ
4
∂
φ
τ
∂2φ2 + ∂2φ
=
.
∂y2
∂
∂x
So effectively our relaxation technique obtains a steady-state 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 definition these time or pseudo-time derivatives
tend towards zero. When they reach zero, we satisfy the rule that was obtained
from physics.
1
9
Problem Descriptions
3
.1 A. Thermo-dynamics
This problem concerns the calculation of temperatures inside a plate of size 0.8m×
.4m as shown in figure 5. Temperature distributions T(x, y) satisfy Laplace’s
equation in stationary equilibrium,
0
ꢁ
ꢂ
∂2
∂2T
∂y2
T2 +
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
fluid, 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 steady-state temperature distri-
bution in the plate specified 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 stage-gate 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 off for the code performance assessment.
edgTe hoef pthlaetehsedatimedenasniodncsoaorleed0.8semct×io0n.4smar.eB0y.2dmefaaunltdt0h.e1mlenagltohnaglotnhge thhoervizeorntitcaal,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 fixed. Boundaries that are described as ‘insulating’
indicate that there is no heat flux through that part of the surface. Where there
is no heat flux, there can be no temperature gradient driving one, so the edge
value and the next-inner-most value must be held to be identical (though with no
constraint on the specific 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 specific 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 steady-state incompressible flow patterns
around obstacles in a two-dimensional channel. Specifically, 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-
isfies Laplace’s equation,
ꢁ
ꢂ
∂2
∂2Ψ
∂y2
Ψ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 steady-state velocity field in the
2D channel specified 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 inflow
velocity, and the position/radius of the cylinders.
For the stage-gate 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 off for
the code performance assessment.
field
, 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 flow speed
Upstream flow speed
(x)
(y)
(U∞)
(V∞)
4m
2m
1ms
0ms−1
−
1
Cylinder locations and radii for the default problem are:
2
2
Cyl xc (m) yc (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
flow through these boundaries. In such cases, the stream function Ψ (which is an
indirect measure of mass flux 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 different, but
individually constant values. For inflow we use the formula,
Ψ(y) = U∞y
(37)
ie. a linear variation of Ψ with y whose gradient is the prescribed upstream flow
speed. It follows that the stream function value on the wall at y = 2m has a
constant value Ψ = 2U∞. At the outflow the value of Ψ must vary smoothly from
Ψ = 0 to Ψ = 2U∞, however disturbances from the cylinders may still affect the
distribution of velocity in the wake. To account for this, we do not enforce a
specific velocity profile, but simply that there is no V component to the outflow
∂
Ψ
velocity. From equation (35) this implies x = 0 and we implement this by holding
∂
the edge value and the next-inner-most value identical (though with no constraint
on the specific value).
3
.3 C. Structural Mechanics
This problem concerns the displacement of a trampoline when subject to the
weight of a gymnast standing off-centre as depicted in figure 7. Stretched mem-
branes can only resist in-plane 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 satisfies Poisson’s equation,
ꢁ
ꢂ
∂2
∂2u
∂y2
u2 +
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 zero-displacement
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) = mAkg
(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 stiffness 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 steady-state displacement of the
trampoline skin specified 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
stiffness of the material.
•
For the stage-gate 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 off for the code performance assessment.
The trampoline has a length of 4m and a width of 2m and is fixed around
the edge. Taking the bottom-left 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 stiffness
value you shoudl choose to scale your computation is largely determined by the
precise details of how you define it (whether by length, area, normalised strain etc).
for our purposes here it is not important how it is defined, so choose an appropriate
value for the material’s extensional stiffness such that the maximum displacement
is aArolulnodf t4h5e−ed5g0ecsmo,fttyhpeictarlaomfptohleindeeflareectfiioxnesd,fosuontdheoyn haarveeadl itsrpalmacpeomliennet. u = 0
everywhere on the edge. The rip is not able to support stress in the direction
across the rip (ie. in the x-direction), but it can support stress in the direction
parallel to the rip just as well as everywhere else. Since out-of-plane displacement
gradients are an elastic response to in-plane stress, then there can be no gradient
of displacement in the x-direction anywhere adjacent to the rip – on either side.
Thus values on the very edge of the rip must be equal to those on the next-inner-
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 difference 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 different values.
4
Optimising your code
4
.1 The rationale
Now that you have passed the Week 5 stage-gate you are expected to have pro-
grams that solve each of the three problems, thermodynamics, fluid 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 re-writing a work-
ing code to make it faster. Good programming practice dictates that you first
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 efficient. The Week
5
stage-gate 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 user-specified accuracy in terms of convergence tolerance and grid
resolution.
2
5
4
.2 Efficiency 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 effect 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 ’on-the-fly’ and the for
operation is encoded as low-level 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 efficiently 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 fixed. 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 efficiently:
inner_x = 2:nx-1;
inner_y = 2:ny-1;
x_diff = phi(inner_y, inner_x+1) + phi(inner_y, inner_x-1);
y_diff = phi(inner_y+1, inner_x) + phi(inner_y-1, 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 efficiency of your program; the gap
between your starting guess and the final 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 finer mesh, maybe 20×10, solve this, then repeat
on 40 × 20, then finally on an 80 × 40. This method is a form of ‘multigrid’,
so-named for obvious reasons, and it one of the fastest and most flexible ways of
solving Poisson-equation problems. The copying and interpolation process from
a coarse grid to a finer 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 fill in the gaps. One hint
for designing your multigrid: refinement 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 non-power-of-2 on either side, so you should consider
carefully your strategy for such cases.
2
7
In all cases you must check that your modifications actually make your code
more efficient… 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 efficiently on very small problems, but
will get progressively more efficient 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 well-crafted
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 efficiency. If the resulting
graph were horizontal, then the time taken scales linearly with the problem size,
and this is the ideal case for an efficient 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 affects every
other point, so a loose analogy might be the handshakes problem (10 people in
a room, 45 hand-shakes 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 files or any other format). The
filename 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 file group007.pdf. No other filenames will be accepted.
The only other stipulation on format is that the .pdf contains less than 10
pages, and note that there is no LATEXtemplate: 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 different 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 data-gathering
by hand, provided a common format is used to invoke your program, and the
solution is returned in a common format. The specification 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 simplified 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 specification
The code submission is 50% of the computational part of Modelling 2. You should
submit a single .m file containing all your code that you wish to be assessed.
The filename 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 file group007.m. Inside this file the name of the first function
must correspond directly to the filename. 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 specified 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 defined 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
finish when the ratio of starting to current gradient is less than the specified
value.
•
•
prob is a single character specifier for the problem you will solve and set-
ting it to ’a’ should solve the thermodynamics problem, ’b’ should solve
the fluid 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 problem-specific 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 specified
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 fluid dynamics case ’b’, the argument variables will determine only the
location and sizes of the circles, and are specified 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 specified 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 pre-set along
with an appropriate stiffness for the membrane to achieve the specified 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 figure
7
. All dimensions and coordinates will be in the units of the original problem
specification.
All information not provided by these arguments must be hard-
coded into the function.
The output arguments are specified 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 vertex-centred 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 first 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 specified 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
confident that the solution is a good response to the input parameters.
Instead return 1 if your code recognises that there is an unsuitable input
specification, or 2 if a failure to converge has occurred during execution.
3
1
5
.4 Calling Script
%
%
modelling 2 2019-20
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=1e-6;
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)=cputime-tstart;
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(phiend-phistart)));
if err<1e-4
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 Rope-Bridge/Weighted-Membrane 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:ny-1
for ix=2:nx-1
phinew(ix,iy)=0.25*phiold(ix+1,iy)+0.25*phiold(ix-1,iy)…
+
0.25*phiold(ix,iy-1)+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 file and returned as argument
•
•
Ensure your code returns a non-zero 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 defined 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 final 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 finally submitted,
is likely to make use of commands from the first two lists: Programming, logic
and arithmetic and Common Functions. There is unlikely to be need to use
commands from the subsequent lists in your final 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
Element-wise 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 pseudo-inverse
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 files 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 file into Matlab workspace
Look for commands or .m files
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 file
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 flux density, q, is
equal to the product of material’s thermal conductivity, k, and the negative local
temperature gradient, −∇T. The heat flux density is the amount of energy that
flows through a unit area per unit time. The reason for the negative sign is that
the heat flows away from hot things towards cold things, ‘downhill’ and therefore
against the temperature gradient as its vector direction is defined. Thus the heat
flux is given by,
q = −k∇T,
(40)
or in one dimension, simply,
q = −k∂Tx .
∂
(41)
4
1
To describe how temperature (or any other quantity) diffuses in time through a
medium, we invoke Fick’s 2nd Law of diffusion, which in one dimension says that
∂∂Tt = ∂∂xq
.
(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 differentiation,
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
−
difference in the amount of q at both extremes of length 2L at this time-instant.
Remember, q is a flux of heat, and so heat could be flowing 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 flowing in
through, say, the −L boundary than, say, is flowing out through the L boundary
of our arbitrary length.
Substituting the heat flux definition (41) into equation (42), we realise we have
a mixed partial differential equation in one dimension,
ꢁ−k
ꢂ
.
(43)
∂∂Tt =
∂
∂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
steady-state behaviour, where the time derivatives become negligible: t → 0. In
such cases, we have the much simpler one-dimensional form
∂2
0 = −k ∂
T
x2
(44)
In the more general case of two dimensions, this becomes,
= −k ꢁ T2 +
ꢂ
,
(45)
∂2
∂2T
∂y2
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 non-zero, then the second term must be
zero. The vector form of this simpler equation is generally known as Laplace’s
equation,
∂2
∂ T
2
T2 + ∂y
2 = 0.
(46)
2
∇
T = ∂
x
4
2
B Stream function descriptions of incompress-
ible flow
Incompressible fluid is a fluid whose volume is a conserved quantity, you cannot
squeeze the fluid into a smaller space (think ‘water’!). This means that in one
dimension (think ‘pipe’!) you can measure the flux of volume (same as the mass
flux for a density of 1) at one end of an (inelastic) pipe, and it must exactly match
the flux of volume at the other. Flux of volume Q and average flow rate U¯ are
connected, Q = U¯A, where A is the cross-sectional 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 one-dimensional flow is that,
Z
L
∂Qx 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 flow can only pass through its faces
at right-angles to them. Again incompressibility means that the total amount of
‘
stuff’ in the box is constant, but for every new bit that flows in, an old bit must
flow out. This time, however, there is no constraint on the direction that fluid can
flow in or flow 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 define 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¯Ax + V¯ Ay. In this context we can choose the shape of
our box, so let’s make it square, Ax = Ay = A = const., and with sides 2L as
before. Thus we could write,
Z
Z
L
L
∂Qx dx +
∂∂Qy dy = 0.
(48)
−L ∂
−L
There’s not much we can do with this equation until we realise that with the flow
∂
∂
Q
∂Q
we’ve defined 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
ꢁ∂∂Qx + ∂∂Qy ꢂ
L
L
∂
Qx + ∂∂Qy 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 flow from the left.
and write down a differential form,
∂∂Qx + ∂∂Qy = 0.
(50)
∂
∂
Q
∂Q
∂y
From the way we’ve defined the flow 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,
∂
Ux + ∂V
∂y = 0,
(51)
∂
which is called the ‘Continuity Equation’ for incompressible flow.
We need to introduce another concept – that of the stream function. Steady
incompressible flow in two dimensions has a special property that its velocity,
a vector field quantity, can be described completely by a separate scalar field, a
single number whose value varies throughout the domain. Returning to our square
box, this single number is a relative measure of the volume flux through a surface
–
or a bit of a surface. It needs a datum somewhere, usually a bottom left-hand
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 flux through the whole left side of the
R
L
box, L Udy = U¯Ax as introduced earlier. The accumulation of stream function
−
for a uniform flow is illustrated in figure 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 flux is V > 0 flow out the positive y-axis direction (just as
postive volume flux is U > 0 flow in the postive x-axis 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 flow 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 flux is still how much ‘stuff’ moves from one side of the line to the other
and the difference in values of Ψ at bottom left and top right reflect 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 differential form
δΨ = Uδy − V δx.
(54)
If we differentiate the stream function, then since x-differentiation has no effect
on y and vice-versa, we have,
∂
Ψ
=
=
−V
(55)
(56)
∂
∂
x
Ψ
U,
∂
y
thus a single field Ψ can describe both components of velocity. But what if
we don’t like the sign convention? We could invent another single field that is
ꢃ
ꢄ
U
V
everywhere oriented at 90 degrees to Ψ: the vector
is always perpendicular
ꢃ
ꢄ
−
V
U
to
. Let’s call our new field φ, defined 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 + ∂y2
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 satisfies incompressibility, by
substituting as follows, and recognising that for smooth fields the order of the
derivatives in x and y can be swapped:
ꢁ
ꢂ + ∂∂y
ꢁ−
ꢂ
∂∂Ux + ∂∂Vy
∂
∂Ψ
∂Ψ
= 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 flow in
free space. However we tell ourselves that this holds approximately everywhere
because it makes obtaining solutions a bit easier, and we cross our fingers and
2
hope that the error isn’t too large. In this problem we will aim to solve ∇ Ψ = 0,
but be aware that the flows it produces are not always physically meaningful. You
should identify in your report where you believe the flows 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 stiffness 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 figure 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 = Fright + Fleft
(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
stiffness 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 stiffness 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 stiffness, 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,
ω = mk
,
(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
mk =
(NK)
M
= N MK =
L
h
K
.
(72)
2
M
N
Thus by dividing equation 70 by m we can write,
u¨ = L2 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 ∂2u
∂x2
u2 =
t
c
,
(74)
∂
for one-dimensional in-plane elastic motion with small strain. The quantity c has
q
K
M
dimensions of velocity and is interpreted as a wave speed, with
a time-scale
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
∂2u2 + ∂2u
u2 =
c
,
(75)
2
∂
t
∂x
∂y2
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
t2
this case, then since c is not normally zero, equation (75) reduces to Laplace’s
equation,
∂2
∂ u
2
u2 + ∂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 vec2tor 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 coefficients in front of φ from φ itself, and put those
coefficients in the matrix that is otherwise zero. Effectively we are saying that the
action of performing a double-derivative 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 benefit 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 element-wise operation being performed
during an iteration of the rope bridge case is the following,
old
i+1
old
φ
+ φ
2 mg
i−1 − (∆x)
2k ,
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
modified 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 − 2k (∆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 definition φˆ
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 |λmax
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
