Exam duration — 3 hours Reading time — 15 minutes
This paper is printed two-sided and consists of six (6) pages of examination.
Identical Examination Papers: There are none.
Authorized Materials: Pens, pencils, rubbers, rulers. Students may use any material from the subject website hosted on the Learning Management System (LMS), any MAST30028 les on the Lab Server and any part of the provided software system MATLAB. Students MAY NOT access their email accounts or any other material hosted on the LMS or any other websites.
Instructions to Invigilators: Each candidate should be issued with an examination booklet, and have access to a desktop computer equipped with MATLAB.
During reading time, students may read the paper, log in to a computer, delete all existing les from C:\…\MATLAB and move any relevant les from the server or the subject website to C:\…\MATLAB. They may start MATLAB but MAY NOT open any M- les.
Instructions to Students: This exam consists of (7) seven questions. Each question states the number of marks it is worth. The total number of marks available is 80.
At the end of the examination, zip up all M- les, gures and documents you wish to include into an archive (Right-click → Send to → Compressed (zipped) folder) labelled with your student number (e.g. NMSCexam_350468.zip) and upload via the exam page of the LMS (like you did with the assignments). The examination booklet will be collected by the invigilator.
This paper may be reproduced and lodged in the Baillieu Library.
Page 1 of 6
Question 1 [8 marks]
Explain, in terms of concepts covered in this course, the output of the
following MATLAB code:
a. format short e
x = realmax*2; x/2 x = realmax+2; x-2
b. format short e; x=4;k=0;
while 4+x ~= 4
x=x/2;k=k+1;
end x
k
Question 2 [5 marks]
Explain, in terms of concepts covered in this course, the output of the
following MATLAB code:
for n = [4 8 12 16]
A = pascal(n);
xTrue = ones(n,1)*10/3;
b = A*xTrue;
format long;
x = A\b
relerr = norm(x – xTrue)/norm(xTrue);
fprintf(’Relative error = %8.4e\n’,relerr)
pause
end
Feel free to modify the code to help you understand what is going on.
Question 3 [8 marks]
In lectures, the code testEuler was run to illustrate several features of
Euler’s Method. Explain in detail which three properties of Euler’s method the code is demonstrating, and how it does this.
You may use any code used in the Labs.
Page 2 of 6
Question 4 [14 marks]
The Simple Monte Carlo method can be used to approximate de nite inte-
grals such as I = b f (x) dx. By using points randomly placed in the interval a
[a, b], the integral can be estimated by
I^ = ( b − a ) f
where f is the sample mean of the function values at the random points
1 N
f = N
a. Write a MATLAB function that returns the estimate I^ as described
above, with de nition function I = simpleMC(f,a,b,N)
b. Since I haven’t told you how to estimate the statistical error for this method from a single simulation, do the following instead. Write a driver examMC.m that
(i) repeats the simple MC simulation M = 1000 times, each one with N = 10, 000 x values, to estimate the integral
^
(iii) and hence computes a 95% con dence interval for I, using the set of
estimates for the simple MC method. You may use mean and std. What properties of the histogram did you use to derive the 95% con –
dence interval?
c. Compare the halfwidth of the con dence interval for the simple MC
method with that obtained from hit-and-miss MC [you may adapt hitmiss.m].
d. How big would N have to be for the halfwidth to be < 0.001, for each method? Comment.
Page 3 of 6
f(xk)
I =
(ii) plots the estimates using a histogram plot
k=1
3 1
exp(−x) dx
Question 5 [16 marks]
The following data: t = 0.5, 1.0, 1.5 · · · 5.0, y =
7.1584e-02 8.5252e+00 1.0868e+01 2.1101e+01 3.1813e+01
4.9507e+01 9.9131e+01 1.6219e+02 2.8002e+02 4.8945e+02
is thought to show exponential behaviour: y ∼ c1ct2
a. What form of plot would con rm exponential behaviour? Brie y explain
why and con rm the behaviour with a suitable plot.
b. Inspired by your plot, transform the data to new variables in which the model appears linear in the parameters. Now t the transformed data to a line by constructing a suitable linear system and using \. Report the corresponding parameters c1 and c2.
c. Now t the original data directly to the exponential model, by using lsqcurvefit. Choose initial guesses for the parameters by using part (b). Compare the parameter values with those from part (b).
d. By plotting the data and both ts on the same plot, decide which method gives the best t for this data.
Question 6 [9 marks]
a. Use Newton's method, the method of bisection and fzero to nd the
root of
For Newton's method, use x0 = 1 and for the other methods use an
f (x) = (sin(x) − x) exp(x) initial interval of [−1, 1].
Use tolerances of 10−10 and a maximum of 40 iterations, where possible.
b. Plot the residuals of each iterate for Newton's method and bisection on a suitable plot and describe what you see. For fzero use the Display option set to iter.
Comment on the performance of all three methods. You may use any code used in the Labs.
Page 4 of 6
Question 7 [20 marks]
The following ODEs describe a model for the spread of an infectious disease
with lifelong immunity
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
s_ = −βsi
i_ = βsi−γi
r_ = γi
where s(t), i(t), r(t) are the fraction of the population susceptible to the disease, infected by the disease and recovered from the disease, respectively. There are 2 parameters β,γ, representing rates of transmission and recovery. Keep γ
xed at γ = 1/7, corresponding to an infectious period of 7 days. Suitable initial conditions are (s0, i0, 0), where s0 + i0 = 1.
a. Write a MATLAB function sir.m, suitable for input into ode23 etc., that describes the system given above, with the parameters set to β = 3/7, γ = 1/7.
b. Write a MATLAB driver that solves the system using ode23, with default tolerances. Test your driver by solving the problem with i0 = 0.01, over the interval [0, 100].
Use plot commands or edit the plot window so that your results look like Figure 1.
ExamQ3
s
r
i
0
0 10 20 30 40 50 60 70 80 90 100
t (days)
Figure 1: The epidemic model for β = 3/7, γ = 1/7.
Page 5 of 6
c. Since s, i, r are the fractions of the population in each class, they must add to 1:
s(t)+i(t)+r(t) = 1
Check how accurately your numerical results satisfy this condition.
Since population fractions can't be negative, use an option to odeset that ensures that s, i, r ≥ 0 and repeat your calculation.
d. Since r can always be found from s and i, we often only plot s (horizontal axis) versus i (vertical axis) (a phase plot). Plot your solution from (b) as a phase plot and use a Textarrow to indicate the initial condition.
The phase plot should look `chunky' i.e. not very smooth. By changing an option to ode23, produce a smoother phase plot.
e. If we add (equal) birth and death rates to this model, we get the equa- tions:
s_ = −βsi+μ−μs i_ = βsi−γi−μi
r_ = γi−μr
A realistic value for μ is μ = 1/(70 × 365), corresponding to a human
lifespan in days, but we will take μ = 1/(100 × 7).
Write a new function M- le sir_birth.m describing this model, includ-
ing the xed parameter μ.
Solve the new model for the case β = 1.5/7,i0 = 0.01 over [0,6000] and describe what new features in the solution are produced by adding birth/death rates.
f. Using tic/toc to measure the time taken to solve the Initial Value Prob- lem from part (e), compare the eciency of ode45,ode15s for this prob- lem. Do your results suggest this is a sti problem?
Reminder: You may use any code used in the Labs. Include your code and any relevant output and plots as evidence to support your answers. You may answer either in the booklet or in a document.
END OF THE EXAMINATION
Page 6 of 6