School of Mathematics and Statistics
MAST30028 Numerical Methods & Scientific Computing
2022
Assignment 1: Simulation & errors
Due:17:00 p.m. September 1
Note: This assignment is worth 20% of the total assessment in MAST30028. When you submit the assignment,
you should submit two files: one pdf file and one zip file. The pdf file contains the answer you write, the
numerical results you generated including figures or tables(0 mark if numerical results are not included in the
pdf file), the comment or the explanation of the results, and the printout of your code (you may use Matlab
publish or Matlab live scripts). The zip file should only contains all the m-files.
To aid marking, please start Questions 3 and 4 on a new page in your write-up.
1 Traffic [20 marks]
Here is a simple stochastic traffic model that incorporates a few simple rules for drivers to follow, yet predicts
features seen in real traffic situations.
The model is discrete in time and space and the car velocities are also discrete. We have N cars distributed
along a circular single-lane road L spatial units in length. Each car takes up 1 spatial unit. Cars are allowed to
travel a maximum speed vmax. In what follows, take vmax = 7.
a. Initially, distribute the cars randomly on the road by sampling the L possible positions without replacement
(there is a command in MATLAB to do this). Set the velocity of each car to zero. Record the distances
between each car and the car in front [remember the road is circular].
b. At each subsequent timestep, perform the following updates on all cars in this order:
(i) accelerate: each car increases its velocity by 1 unit if it is not travelling at the maximum speed.
(ii) avoid collisions: each car checks the distance d to the next car. If its velocity v ≡ d then it reduces
its speed to d? 1 to avoid a possible collision in the next timestep.
(iii) random braking: if the car is moving, it reduces its speed by 1 unit with probability p. This is the
only stochastic part of the model.
You should seed the random number generator with the integer 46 if the logical variable seed is true.
This is to aid marking.
(iv) update positions (and distances): move all positions forward by v units [remember the road is circular]
and update all distances.
c. The model is run for a number Tr timesteps to relax from the initial state, then for a further Te timesteps
in which the positions of each car are recorded at each timestep. It is convenient to record the position of
the cars as rows of a matrix trafficflow (initially full of zeroes) using a command such as
trafficflow(time,position) = 1;
where position is a row vector recording the position of all the cars at t = time.
1
d. Write a MATLAB function to simulate this model, with the first line
function trafficflow = trafficSim(L,N,p,vmax,Tr,Te,seed)
You can＊t avoid loops over time but I was able to write the updates in section b above without any loops
over cars, by using MATLAB commands operating on vectors and logical indexing. More marks will be
given if you can avoid update loops (and the code is correct).
e. Test your code by:
(i) writing a driver to simulate 1 car on a road of length 150, with p = 0.3, and seed = true. Set
Tr = 350 and Te = 150. You should get clear traffic
(ii) edit your driver to simulate 55 cars under the same conditions. You should obtain,
using imshow(~trafficflow), the image below.
f. Finally run your code for a road of length 1500 with p = 1/3, and seed = false. Set Tr = 3500 and Te =
1500. Perform simulations for various values of traffic density 老 = N/L, including 老 = 0.1, 0.15, 0.2, 0.25.
Describe what you see and interpret it in terms of common traffic patterns.
2
2 Hyperspheres [12 marks]
In lectures I mentioned we could estimate 羽 by using hit-and-miss Monte Carlo estimation on a unit circle.
Here we study the n-dimensional analogue 〞 unit hyperspheres i.e. unit balls in Rn. This allows us to test the
abilities of MC methods on higher-dimensional integrals.
We want to estimate the volume of hyperspheres with radius r = 2 in dimensions n = 2, 5, 8, 11, 14. To do
this, we distribute points in the hypercube [?2, 2]n randomly, and estimate the probability of landing inside the
unit hypersphere, p. Then the volume of a hypersphere Sn is given by
Sn = 4
np
since 4n is the volume of the hypercube with side length 4.
a. Write a MATLAB function with calling sequence
function [estimate, halfWidth] = ballMC(Nreps,n)
that uses hit-and-miss Monte Carlo with Nreps random points in Rn to estimate p using the sample
proportion P? . It returns the estimate p? and the half-width of the 95% confidence interval for p.
b. Write a driver that, for each dimension in n = 2, 5, 8, 11, 14 calls ballMC with Nreps taking values dis-
tributed logarithmically between 101 and 106. The driver should plot
(i) the estimated Sn versus Nreps on a suitable plot
(ii) the CI halfwidth (as a measure of statistical error) for Sn versus Nreps on a suitable plot
i.e. 10 plots in all. Include suitable labels, and titles including the value of n.
c. The volume Sn is known exactly:
Sn =
2n羽n/2
忙(1 + n/2)
Modify your driver so that the plots of statistical error also include the actual error of each estimate. I
suggest you use at least 100 values of Nreps for each n.
d. Based on your plots, answer the following:
do the Monte Carlo estimates stabilize as Nreps increases?
what happens to the estimated Sn as n increases?
does the CI halfwidth capture the actual error behaviour?
what happens to the CI halfwidth for Sn as n increases?
is there anything unusual in your plots for high n and low Nreps?
This example illustrates why one may need to use variance reduction techniques in high dimensions.
3
3 A bound on roundoff error [4 marks]
To analyse the roundoff error from n floating point multiplies, the following quantity arises:
旭ni=1(1 + 汛i) √ 1 + 牟n
where |汛i| ≒ u and u is unit roundoff i.e. 牟n is the relative error after n floating point multiplies.
Use mathematical induction to prove the bound
|牟n| < nu
1? nu
assuming nu < 1. The same bound holds (with ≒) if we include divides as well, but you do not need to prove
this.
4 A better derivative. [14 marks]
a. Use Taylor series to derive the truncation error of the approximation
f ∩∩(x) ＞ f(x+ h)? 2f(x) + f(x? h)
h2
assuming f ﹋ C4.
b. Explain why dividing by h produces no roundoff error if h = 2?k, k ﹋ Z.
c. Assuming that built-in functions return exact answers but that addition/subtraction produce roundoff
errors and that h = 2?k, k ﹋ Z, show that the roundoff error RE in computing this expression satisfies
the bound
|RE| ≒ (K1 +K2h)u/h2
where the constants K1,K2 depend on f and x.
d. Hence estimate, in terms of u, the optimal choice for h. You may assume hu? u.
e. Write a script file called SecondDerivative.m to test your estimates for f(x) = exp(x) and x = 1,
h = 0.5, 0.25, . . . 2?52 and plot the errors, as in ForwardDifference.m. Comment on your results.