Numerical Methods & Matlab

ENG2093 – Numerical Methods coursework – 2017/18
This coursework is worth 25% of the ENG2093 Experimental and Numerical Methods
Instructions – IMPORTANT: please read carefully
The coursework is divided in two parts. You should already have started engaging with the first part [20%]
involving weekly readings and multiple-choice tests or coding assignments to complete (see “Weekly
assignments” folder on the ENG2093/Numerical Methods SurreyLearn pages). The second part [80%] is
presented in the following pages and it involves solving two engineering problems by applying numerical
methods and, in some cases, Matlab programming. Coding is only a minor part of the assessment (please see
the marking scheme for each problem). That means that it is more important to answer the questions correctly
by using a Matlab code that works, rather than producing the best possible computer programme. Some parts
of the coursework won’t depend on Matlab programming, but on mathematical derivations and written
explanations and comments. Matlab must only be used where specified to obtain numerical answers or plots.
Your coursework paper must include all the answers to the questions, including numerical values, plots and
comments. Answers not reported on the main coursework document (i.e. where only the Matlab code is
present) won’t be considered. The coursework paper must be submitted as an electronic document (either a
Word document, .docx or .doc, an OpenOffice/LibreOffice document, .odt or a PDF file .pdf). Mathematical
formulas can be entered using the Word built-in “Equation” tool (or OpenOffice/LibreOffice equivalent) or by
digitalising (picture or scan) a hand-written paper and including it into the main .docx/.doc/.odt/.pdf file. The
former method is recommended and preferred. Sketches can either be produced using an external software
(e.g. PowerPoint) or, again, hand-written and then scanned/photographed. Please submit only one document
containing the whole coursework paper. If you submit multiple documents I will only assess the latest
submission unless otherwise specified.
You must submit all your Matlab codes as m-files. Do not copy them onto the main document. If your
numerical answers and plots are not supported by the code you’ve submitted, then they won’t be considered
correct. You can submit multiple files on SurreyLearn, even at different points in time. As the University has a
site licence for Matlab, if you do not already have a copy, or have access to a copy, you should be able to get
one for personal University use. To obtain this go to https://www.surrey.ac.uk/fepsit/services/ and follow
instructions. Please contact IT services if you have a problem: I am afraid I cannot help with this.
Several questions will require using data values that are different depending on you URN, and in particular, the
last two digits. Make sure you use the right values from your URN. Double check them before submitting. If
you fail to use the correct URN values, I’m afraid I cannot help.
Any evidence of plagiarism or collusion with others will be dealt with severely under the University rules.
Please ensure that what you submit is your work and only your work; do not “assist” or seek assistance from
others by e.g. “sharing” parts of programmes, and do not under any circumstances copy code from
books/websites. Both the coursework paper and the code will be submitted to specialist plagiarism-detection
software to spot potential problems.
Submit your coursework paper and Matlab codes (m-files) to your “Numerical Methods Coursework”
assignment folder on SurreyLearn. The deadline is Tuesday 20th March 2018 at 16:00. Do allow plenty of
time for submission, as the process will likely take a while: do not try to submit at the very last minute. The
time of submission of the last file is the one that will be considered for the whole coursework for late
submission purposes. Please submit all the individual files separately, do not compress them in a single .zip
file.Data for use in coursework
Please include the following info at the beginning of you coursework document:
Name:
Surname:
URN:
Last two digits from the URN:
6 x x x x x x
aURN bURN
aURN =
bURN =
Please also include the COURSEWORK Cover Sheet
You can find it on SurreyLearn in the “MES UG Programme Information” pages, in the folder
“Coursework schedules and submission forms”Problem 1 [40 marks]
One of the classical laws of planetary motion due to Kepler says that a planet revolves around the
sun in an elliptic orbit as shown in the figure below:
Figure 1 – Problem 1
Suppose one needs to find the position (�. �) of the planet at time �. This can be determined by the
following formulae:
� = � cos(� − �)
� = �√1 − �2 sin(�) (1)
where � is the semi-major axis of the ellipse, � is the eccentricity of the ellipse and � is called
eccentric anomaly and it depends on time (it can vary from 0 to 2� radians – don’t worry too much
about the meaning of those terms. It isn’t very important as far as the calculations are concerned).
To determine the position (�, �) one must know the value of �, which can be computed from
Kepler’s equation of motion:
� = � − � sin(�) (2)
where � is the mean anomaly (again, it can vary with time from 0 to 2� radians and the actual
meaning is not very important here, however if you’re curious you can find more in the Wikipedia
page of the Kepler’s laws of planetary motion), and 0 < � < 1. This equation, thus, relates the
eccentric anomaly, �, to the mean anomaly, �. Thus to find � we can solve the nonlinear equation:
�(�) = � − � + � sin(�) = 0 (3)
The following table shows the eccentricity (�) and the semi-major axis (�) values for the orbits of
several planets and dwarf planets in the solar system.Planet / Dwarf planet Eccentricity � Semi-major axis � [106 km] ����
Mercury 0.2056 57.909 0
Venus 0.0069 108.208 1
Earth 0.0167 149.598 2
Mars 0.0934 227.939 3
Ceres 0.0758 414.010 4
Jupiter 0.0489 778.57 5
Saturn 0.0565 1433.53 6
Uranus 0.0464 2875.04 7
Neptune 0.0095 4500.00 8
Pluto 0.2488 5906.38 9
Questions for Problem 1
1) Show that the simple (one-point) iteration method converges when using the following
formula:
� = �(�) = � + � sin(�)
[4 marks]
2) Use the data for the planet/dwarf planet corresponding to your ���� value and � = 1 + ����
10
(in radians), where ���� and ���� are the last two digits from your URN. Choose an
appropriate initial approximation for open methods by plotting the function �(�) using
MATLAB. Please report your plot and choice for �0 in the coursework document and submit the
MATLAB code you’ve used to produce the plot.
[4 marks]
3) Explain the basis of the Newton-Raphson method, using a sketch if helpful.
[3 marks]
4) Using the same data as in question 2, we now want to solve equation (3) using the NewtonRaphson method. Show the calculations (by hand)

for the first two iterations, including the
approximate relative error at each iteration. Work to 5 decimal digits.
[5 marks]
5) Calculate the solution to question 4 above using MATLAB using a tolerance (calculated on the
approximate relative error) of 10−9. Report your answer in the coursework document and
submit the MATLAB code you’ve used to calculate the answer.
[5 marks]
6) Investigate the behaviour of the approximate relative error as the number of iterations
increases in both the simple (one-point) iteration and Newton-Raphson methods in solving the
above problem. Comment on this behaviour: how do the two methods compare? Are the errors
decreasing as you expected? Are they not? Why? Also, please discuss the relationship betweenthe errors you are calculating and the true

errors. Ideally, use MATLAB to analyse the errors
iteration-after-iteration (in this case, submit the MATLAB code you’ve used to do this analysis,
in addition to describing the methodology you’ve used). However, you can still answer this
question (and getting partial marks) by using calculations by hand.
[9 marks]
7) Referring to equations (1) and (3) and using the same data as in question 2, write a MATLAB
programme to plot the values of � and � (on the same plot) against non-dimensional time � = �

(where � is the orbital period of the planet/dwarf planet: i.e. if the orbital period is 365 days, a
� value of 0.5 means 182.5 days, a value of 1 means 365 days, and so on). Please report your
plot in the coursework document and submit the MATLAB code you’ve used. In order to solve
this question, you need to be aware that in equations (1), (2) and (3), E and M are timedependent. Specifically, M varies uniformly from 0

when � = 0 to 2� radians when � = 1. So
you’ll need to calculate the corresponding E values (equations (2) and (3)) for a range of M
values. These E values can be used to calculate the corresponding values for x and y in equation
(1). This problem can get a bit complicated. If you’re not able to complete it successfully, you
can still get partial marks by describing the procedure. I recommend you describe your
procedure anyway, even if you submit the code, as your MATLAB programme might still be
wrong.
[10 marks]
A note on marking the questions where MATLAB code is requested:
Marks will be given for correct answers, of course, but also for how general your code is. For
example, a neat MATLAB function that solves the problem using any value for the tolerance (e.g. by
passing it as an input argument) or even for any f function (or g function, or whatever
formula/variable you need for that particular method) will be given more marks than for a
programme that only solves your particular problem. It is up to you to judge how general your code
is. My recommendation for people who are not very skilled at MATLAB is to go for just solving the
problem, without paying too much attention to these details (you’ll still get at least 40 to 75% of the
marks, depending on the question, by solving the problem correctly). On the other hand, if you are
very skilled with MATLAB, what I’m describing will make more sense. In any case, try to walk before
you run. A correct basic code is way better than a wrong fancy code… Code commenting and
indentation will also be (a small) part of the assessment.Problem 2 [40 marks]
Figure 2 – Problem 2
If the velocity distribution of a fluid flowing through a pipe is known (see figure above), the flow rate
� (that is, the volume of water passing through the pipe per unit time) can be computed by:
� = ∫ � �� (4)
where � is the velocity and � is the pipe’s cross-sectional area (to grasp the meaning of this
relationship physically, recall the close connection between summation and integration). For a
circular pipe, � = ��2 and �� = 2�� ��. Therefore,
� = ∫ �(2��)��
�0
0
(5)
where � is the radial distance measured from the centre of the pipe. Let’s assume that the velocity
distribution (in m/s) is given by
� = 2 (1 − ��0)
16
(6)
where �0 is the total radius.
Questions for Problem 2
1) Two methods to solve numerically the integral in equation (5) are the trapezoidal rule and the
Simpson’s 1/3 rule.
a. Explain the basics of the single-application trapezoidal rule, using a sketch if helpful.
b. Using a Taylor series expansion, show that the truncation error for the singleapplication trapezoidal rule is �(ℎ3).c. Both the methods

above (trapezoidal and Simpson’s) include a “single-application”
version and a “multiple-application” version. Explain the differences between the
two versions (single- and multiple-application) using either the trapezoidal or the
Simpsons’s rule as an example. Also, describe and discuss about the behaviour of
the truncation error when going from the single-application version to the multipleapplication version of the method.
[9 marks]
2) Consider a value for �0 = 2 + ���� + ����
10
(in mm), where ���� and ���� are the last two
digits from your URN. Use MATLAB to plot the velocity distribution against � with 0 < � < �0.
Please report your plot in the coursework document and submit the MATLAB code you’ve used
to produce the plot.
[3 marks]
3) Using the �0 value as in question 2, we now want to solve equation (5) using a numerical
integration method. Using sub-intervals with step-size ℎ = �0/6 and working to 5 significant
figures solve the equation numerically by hand, using:
a. The multiple-application trapezoidal rule.
b. The multiple-application Simpson’s rule.
[5 marks]
4) Calculate the solution to question 3 above (using both the multi-application trapezoidal and
Simpson’s rules) using MATLAB and sub-intervals with step-size ℎ = �0/1000. Report your
answers and submit the MATLAB code you’ve used to calculate them.
[5 marks]
5) Assuming the analytical solution of the integral in equation (5) is, in mm3/s:

� = [− 24000 91 � (1 − ��0)1 6 (� − �0)(7� − 6�0)]0 �0
calculate the true relative error for the results for both methods (either from Q3 or Q4).
Ideally, use MATLAB to calculate the errors (in this case, submit the MATLAB code you’ve used
to do the calculations, in addition to describing the methodology you’ve used). However, you
can still answer this question (and getting partial marks) by using calculations by hand.
[6 marks]
6) Let’s consider a slightly different problem. The pressure gradient for laminar flow through a
constant radius tube is given by
��
�� = −
8��
��4
(7)
where � = pressure (Pa), � = streamwise distance along the pipe’s centre line (m), � = dynamic
viscosity (N s/m2), � = flow (m3/s), and � = radius (m).a. Using the method of your choice (justify the choice), determine the pressure

drop for 10
cm length tube for a viscous liquid (� = 0.005 N s/m2, density = � = 1 ⋅ 103 kg/m3)
with a flow rate of � = (10 + ����) ⋅ 10−6 m3/s and the following varying radii along
its length.
x [cm] 0 2 4 6 8 10
r [mm] 2 1.35 1.34 1.58 1.49 2
b. Compare your result with the pressure drop that would have occurred if the tube had a
constant radius equal to the average radius.
c. Determine the average Reynolds number for the tube to verify that the flow is truly
laminar (�� < 2100), where
�� =
���

(8)
� is the (average) velocity, and � is the diameter.
Ideally, use MATLAB to carry out the above calculations (in this case, submit the MATLAB code
you’ve used to do the calculations, in addition to describing the methodology you’ve used).
However, you can still answer this question (and getting partial marks) by using calculations by
hand.
[12 marks]
A note on marking the questions where MATLAB code is requested:
Marks will be given for correct answers, of course, but also for how general your code is. For
example, a neat MATLAB function that solves the problem using any value for the tolerance (e.g. by
passing it as an input argument) or even for any f function (or g function, or whatever
formula/variable you need for that particular method) will be given more marks than for a
programme that only solves your particular problem. It is up to you to judge how general your code
is. My recommendation for people who are not very skilled at MATLAB is to go for just solving the
problem, without paying too much attention to these details (you’ll still get at least 40 to 75% of the
marks, depending on the question, by solving the problem correctly). On the other hand, if you are
very skilled with MATLAB, what I’m describing will make more sense. In any case, try to walk before
you run. A correct basic code is way better than a wrong fancy code… Code commenting and
indentation will also be (a small) part of the assessment.
M. Carpentieri 27/2/2018