UNIVERSITY OF MELBOURNE
SEMESTER 2 ASSESSMENT, NOVEMBER 2017
SCHOOL OF MATHEMATICS AND STATISTICS MAST30028 NUMERICAL METHODS AND SCIENTIFIC COMPUTING Exam duration — 3 hours
Reading time — 15 minutes
This paper is printed two-sided and consists of seven (7) 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 files 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 files from C:\…\MATLAB and move any relevant files from the server or the subject website to C:\…\MATLAB. They may start MATLAB but MAY NOT open any M-files.
Instructions to Students: This exam consists of (6) six 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-files, figures 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 7
Question 1 [8 marks] The following code
close all;
f=@(t,y) -4*pi^2*(y-cos(t))-sin(t);
y0=1; tspan=[0 2];
n=1000;
[tvals,yvals] = FixedEuler(f,tspan,y0,n);
plot(tvals,yvals,’o’);
xlabel(’t’);ylabel(’y(t)’);
hold on;
n=100;
[tvals,yvals] = FixedEuler(f,tspan,y0,n);
plot(tvals,yvals,’x’);
n=10;
[tvals,yvals] = FixedEuler(f,tspan,y0,n);
plot(tvals,yvals,’r’);
produces the following plot.
5 x 10
1
0
−1
−2
−3
−4
−5
−6
−7
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
t
Which solutions do you think are accurate? Any? Explain in general terms what the code is doing (not a line-by-line description) and, in terms of concepts covered in this subject, explain what is happening in this problem.
You may use any code used in the Labs.
Page 2 of 7
y(t)
Question 2 [8 marks]
One way to estimate π is to count the proportion of points thrown randomly
at the unit square [0, 1] × [0, 1] that lie inside the unit quarter-circle.
a. Explain how the proportion can be used to estimate π/4 (and hence π).
b. Write a MATLAB function examQ2.m that estimates π/4 using a Monte Carlo method as above. Your MATLAB function should return your estimate and the half-width of the 95% confidence interval, given the number of points as an input argument.
c. Determine by simulation how many points are required for the statistical error (measured by the half-width) to be be less than 1.25×10−3 i.e. the estimate for π should have 2 decimal place accuracy.
d. Estimate how many points you would need for 3 decimal place accuracy. Explain your answer.
Question 3 [16 marks]
The following data: t = 0.5, 1.0, 1.5 · · · 5.0, y =
2.3564e+00 2.2169e+00 8.0175e+00 1.6653e+01
5.4348e+01 8.5735e+01 1.2883e+02 1.8279e+02
3.1486e+01
2.5020e+02
is thought to show power law behaviour: y ∼ c1tc2
a. What form of plot would confirm power law behaviour? Briefly explain
why and confirm 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 fit the transformed data to a line by constructing a suitable linear system and using \. Report the corresponding parameters c1 and c2.
c. Now fit the original data directly to the power law 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 fits on the same plot, decide which method gives the best fit for this data.
Page 3 of 7
Question 4 [12 marks]
a. By timing the solution of a system of 4000 random linear equations, compare the LU and QR factorizations in terms of efficiency. What benefits, if any, does the QR factorization have over the LU factorization?
Include the factorization step in your timing and use explicit calls to lu and qr.
b. Show that the solution of a square linear system can be written as x = Vy
where
Σy = UT b
and U, Σ, V come from the Singular Value Decomposition of A. Hence time the solution of the same random system from part (a) using svd and comment.
c. Explain 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
Page 4 of 7
Question 5 [14 marks]
a. One variant of Newton’s method, called the chord method, is just New- ton’s method but with the derivative function f′ being evaluated only at x0 and re-used, rather than being re-computed at each new iterate.
Modify Newton.m to produce a function M-file Chord.m that implements this method.
b. Another variant consists of using the forward difference approximation f(x+δ)−f(x)
δ
to approximate f′(x). What is the advantage in such a variant? Explain
why a sensible default value of δ might be δ = ε1/2. M
Modify Newton.m to produce a function M-file NewtonFD.m, with δ passed as a second parameter in the vector parms and with the above default value.
c. Compare the performance of Newton’s method with that of Chord.m and NewtonFD.m on the equation
sin(x) − cos(2×2) = 0
with x0 = 1 and a relative tolerance (only) of 10−11. For NewtonFD try both δ = 0.01 and the default value.
Suitable plots may be useful. Based on your results, describe the be- haviour of the chord method and NewtonFD compared to Newton’s method.
You may use any code used in the Labs.
Page 5 of 7
Question 6 [22 marks] The following ODE
y′′ +δy′ −y+y3 =γcos(t)
describes a mass moving in a certain potential, subject to periodic forcing. The parameter δ is the damping, and γ gives the amount of periodic forcing with period 2π. Here we investigate the solutions for the initial conditions
y(0) = y′(0) = 0.5
a. Write a MATLAB function examQ6.m, suitable for input into ode23 etc., that describes the ODE given above, with the parameters set to δ = γ = 0 inside the function.
b. Write a MATLAB driver Q6driver.m that solves the system using ode45, with default tolerances. Test your driver by solving the problem over the interval [0, 20].
Use plot commands or edit the plot window so that your results look like Figure 1.
Figure1: Asolutionforδ=γ=0.
Since the solution looks periodic, a phase plot should show a closed orbit.
Produce a phase plot for the solution above.
Page 6 of 7
c. Repeat the calculation from part (b) but with the solution from 2 initial conditions on the same plot: y(0) = 0.5, y′(0) = 0.4, 0.5. Show the results using two suitable plots and describe what you see.
d. Now add damping by setting δ to 0.3 (by editing examQ6.m). Using two suitable plots, explain how the system behaves in this case.
e. If we add 2π-periodic forcing by setting γ = 0.5, we get more complicated behaviour.
Repeat the calculation from part (d) with damping and 2 different initial conditions, but over the interval [0,100] and with γ = 0.5. Show the results using two suitable plots and describe what you see.
f. Since the forcing is periodic with period 2π, it is useful to obtain output of the solution every 2π units of time and plot the values as a phase plot. Use a feature of ode45 to find the solution with δ = 0.3, γ = 0.5 every 2π units of time over the interval [0, 4000π], starting from y(0) = 0.5, y′(0) = 0.5.
Plot the solution as a phase plot, using a point marker with MarkerSize set to 1 and no line. You should obtain an intriguing pattern signalling chaotic solutions.
g. Using tic/toc to measure the time taken to solve the Initial Value Prob- lem from part (f), compare the efficiency of ode45,ode15s for this prob- lem. Do your results suggest this is a stiff 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 7 of 7