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