BEST代写-线上编程学术专家

Best代写-最专业靠谱代写IT | CS | 留学生作业 | 编程代写Java | Python |C/C++ | PHP | Matlab | Assignment Project Homework代写

Matlab代写 | MENG21712 Modelling with computers

Matlab代写 | MENG21712 Modelling with computers

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 h0 we let hx, 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 i1 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  
i1 .  
(13)  
1
i2  
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φi1  
x  
2  
φ
x
.
(14)  
(15)  
(16)  
x2  
x  
i
Simplifying, we have,  
2  
x2  
φ
(φi+1 φi) (2φi φi1  
(∆x)  
i+1 2φi + φ  
) = φ  
i1 .  
2
(∆x)  
i
So the discretisation of our rule, equation (1), is as follows:  
kφ  
i+1 2φi + φ  
i1 = 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 + φi1  
=
=
(17)  
(18)  
k
φi+1 + φ  
2
(∆x) mg.  
φi  
i1 −  
2
2k  
This statement now explicitly defines a value of φi that is dependent on its neigh-  
bours φi1 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  
φi1 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  
i1 −  
φ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 + φ  
i1 ,  
(22)  
2
(∆x)  
i
and the y derivative follows suit:  
2  
y2  
φ
φj+1 2φj + φ  
j1 .  
(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 ≈  
i1,j +  
i,j1 . (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 (φi1,j, φi+1,j, φi,j1, φi,j+1,x,y) .  
(25)  
and then repeatedly compute new values of φi,j until they converge. Multiplying  
2
2
equation (24) by ∆xy to get a common denominator, we have,  
2
2
0
y (φi+1,j 2φi,j + φi1,j) +x (φi,j+1 2φi,j + φi,j1)  
(26)  
(27)  
and teasing out the φi,j terms, we have,  
2 2 2 2  
y φi,j + 2∆x φi,jy (φi+1,j + φi1,j) +x (φi,j+1 + φi,j1) .  
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 + φi1,j) +x2 (φi,j+1 + φi,j1).  
(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 + φi1,j + φi,j+1 + φi,j1 .  
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
i1,j  
n
i,j+1  
n
i,j1 .  
φ
+ φ  
+ φ  
+ φ  
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
i1,j  
n
n
i,j1  
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
i1,j  
n
n
i,j1  
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  
0ms1  
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) = Uy  
(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 t4h5eed5g0ecsmo,fttyhpeictarlaomfptohleindeeareectioxnesd,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  
– * /  
– *  
* ./  
=
=
=
for  
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 = kT,  
(40)  
or in one dimension, simply,  
q =kTx .  
(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 andy 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 h0. 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 + φi1  
(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 = 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  
i1 (∆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
φ =ˆ  
ˆ
2
1
φ =ˆ  
ˆ
3
2
φ =ˆ  
ˆ
(
83)  
4
3
φ =ˆ  
ˆ
5
4
φ =ˆ  
ˆ
.
.
.
k+1  
k
The general case is φˆ  
= ˆ , 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
.φˆ = ˆ.  
(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
bestdaixie

评论已关闭。