Main Examination period 2019/20
MTH739U / MTH739P: Topics in Scientific Computing
2019/20 – Coursework 2
Assignment date: Monday 25/11/2019
Submission deadline: 20/01/2020 at 23:55
The coursework is due by Monday, 20th January 2020 at 23:55. Please submit a report
(preferably in pdf) containing the answers to the questions, complete with written explana-
tions and printed figures. All figures should contain a title, axes labels, and a legend (if more
than one curve are in the same figure). It is necessary to provide a sufficiently detailed expla-
nation of all the original algorithms used in the solution of the questions (except for material
explicitly discussed in the lectures, e.g. the Runge-Kutta method). You normally need to
show that your program works using suitable examples. All the code produced to answer
each question should be submitted in a single zip file, aside with the report. It might be useful
to organise the code in different directories, one for each question. Only material submit-
ted through QMPlus will be accepted. Late submissions will be treated in accordance
to current College regulations. Plagiarism is an assessment offence and carries severe
penalties. See the accompanying guidelines for additional information.
Part A: Coursework [45 marks]
Question 1. [10 marks]
Generating random numbers.
For the probability distributions detailed below, construct functions to obtain random numbers
from these distributions, and test these functions by creating histograms from a sufficient
number of samples and plotting the histograms together with the given distribution functions.
(a) Uniform distribution over the interval [−2pi, pi]. [2]
(b) Uniform distribution over the union of the three intervals [1, 2] ∪ [3, 4] ∪ [5, 6]. [2]
(c) Gaussian distribution with a given mean value µ and variance σ2. Use either the
Box-Muller method or the Marsaglia’s polar method, and test your function by sampling
from a Gaussian with µ = 12.5 and σ = 3. [2]
(d) Continuous distribution with probability density function:
f(x;λ) = λe−λx
for x ≥ 0, where λ is a parameter provided by the user. Test your function by sampling
from the three distributions obtained by setting λ respectively equal to 0.7, 1.5, and 3.5. [2]
1
(e) Continuous distribution with cumulative density function:
F (x) =
1
6
(
x2 + x
)
for x ∈ [0, 2]. [2]
Question 2. [10 marks] [10 marks]
Acceptance-Rejection Sampling.
We want to sample random variables from the probability density function:
f(x) =
1
1000
(
100 +
4
5
x3 − 6x2
)
x ∈ [0, 10] (1)
using Acceptance-Rejection sampling. Hence, we need to find a function g(x) such that
g(x) ≥ f(x) ∀x ∈ [0, 10] and we are able to draw random samples from g(x)/A where
A =
∫ 10
0
g(x)dx.
(a) Construct a function which returns random samples distributed according to f(x) using
the acceptance-rejection method with g(x) = 0.3 ∀x ∈ [0, 10]. [3]
(b) Use the function you constructed in point (a) to obtain N = 10000 samples from f(x).
Produce a histogram of the samples, and compare it with the expected distribution. [3]
(c) Construct a function to obtain random samples distributed according to f(x) using the
acceptance-rejection method with
g(x) = a+ b(x− 5)2
where a and b are two constants whose value should be appropriately determined by the
student in order to reduce the rejection rate as much as possible. [2]
(d) Use the function you constructed in point (c) to obtain N = 10000 samples from f(x).
Compare the distribution of the samples with the expected one, and report the value of
the rejection rate. [2]
Question 3. [10 marks] Importance sampling.
Consider the integral
I =
∫ 1
0
x3(1− x)1/2dx = 32
315
which is the value of the Euler Beta function:
B(a, b) =
∫ 1
0
xa−1(1− x)b−1dx
for a = 4 and b = 3/2.
2
(a) Construct a function to determine a Monte-Carlo estimate IU of I using N = 1000
uniform variates in [0, 1], and provide the absolute value of the estimate and the
corresponding Mean Squared Error. [3]
(b) Determine a Monte-Carlo estimate of I using importance sampling. Choose N = 1000
random variables with a probability density function f(x) = 5x4, x ∈ [0, 1] as sampling
points. Provide the absolute error of the estimate. How much does this estimate
improves/deteriorates over the one obained in point (a) through simple sampling? [3]
(c) Now repeat point (b) by trying the three functions f1(x) = 4x3, f2(x) = 3x2 and
f3(x) = 2x as the sampling distributions for importance sampling, considering N = 1000
samples in each case. Compare the results with those obtained in point (a) and (b). Which
of the four estimates is the most accurate? Can you provide an explanation of the results
you obtained? [4]
Question 4. [15 marks]
Simple and biased random walks.
(a) Generate 100 sample paths of the Bernoulli random walk defined by
Xi+1 = Xi + Y,
where Y is the jump variable such that P (Y = 1) = 1/2, P (Y = −1) = 1/2.
(i) Plot 10 random paths up to N = 100 time-steps (i = 0, 1, . . . , N ) using X0 = 0. [2]
(ii) Construct a histogram of the position reached at time i = 5 of 100 sample paths
and compare your result with the corresponding binomial distribution. [2]
(iii) Construct a histogram of the position reached at time i = 20 of 100 sample paths
and compare your results with the corresponding binomial distribution. [2]
(iv) Demonstrate what happens to the comparison between the binomial distribution
and the histogram when you increase the number of random walkers to 1000 and
discuss the results. [2]
(b) Generate 100 sample paths of the Bernoulli random walk defined by
Xi+1 = Xi + Y,
where Y is the jump variable such that P (Y = 1) = p, P (Y = −1) = 1− p.
(i) Plot 100 random paths up to N = 100 time-steps (i = 0, 1, . . . , N ) using X0 = 0
and p = 0.2. [3]
(ii) Construct a histogram of the position reached at time i = 5 of 100 sample paths
and compare your result with the corresponding binomial distribution with bias p.
Analyse your results. What happens if you change p? [4]
3
Part B: Programming Project and Report [55 marks]
Question 5. [55 marks]
Self-avoiding walks.
A self-avoiding walk on a regular lattice is a lattice path starting at the origin that does not
visit any site of the lattice more than once. For example, on the square lattice the number of
N -step self-avoiding walks is given by
1, 4, 12, 36, 100, 284, 780, 2172, 5916, 16268, 44100, . . .
for N = 0, 1, 2, . . . [1]. For an easily accessible introduction to self-avoiding walks, see [2].
There are numerous algorithms available to estimate the number of different self-avoiding
walks with N steps [3]. One could enumerate all N -step self-avoiding walks recursively, but
their number grows so quickly in length that this is a really inefficient algorithm.
For this programming project, you will estimate the number of self-avoiding walks on the
square lattice by considering a class of algorithms based on incomplete, or perhaps better,
stochastic enumeration. This idea goes back sixty years [4]. A modern exposition is found in
[5], which you should largely follow.
(a) Write a function recursiveSAW(N) to recursively enumerate all self-avoiding walks
of length 1..N . Using this function, reproduce the numbers in the sequence given above
for lengths up to 10. [10]
(b) Write a function simplesamplingSAW(N) to use simple sampling to estimate the
number of self-avoiding walks of length 1..N . Test your function by comparing its
output against the exact values for lengths up to 10, and give an estimate of the accuracy
of your result. Can you use your function to estimate the number of 50-step
self-avoiding walks? [10]
(c) Write a function RosenbluthSamplingSAW(N) to use Rosenbluth sampling to
estimate the number of self-avoiding walks of length 1..N . Test your function by
comparing its output against the exact values for lengths up to 10, and give an estimate
of the accuracy of your result. Can you use your function to estimate the number of
50-step self-avoiding walks? [15]
(d) Using all three methods, produce a table containing your results on the number of
N -step self-avoiding walks on the square lattice for lengths up to 100 steps, including
error estimates where appropriate. (You will find that you cannot fill all entries in this
table.) [10]
(e) Write a report describing these three algorithms and the details of your implementation.
Reflect on the relative strengths and weaknesses of these algorithms. Back up your
conclusions with evidence from specific data such as statistical errors and run times. [10]
4
[1] Number of n-step self-avoiding walks on square lattice, in The Online Encyclopedia of
Integer Sequences (http://oeis.org/A001411)
[2] B. Hayes, “How to avoid yourself,” American Scientist 86 314-319 (1998)
(https://www.jstor.org/stable/27857052)
[3] E. J. Janse van Rensburg, “Monte Carlo methods for the self-avoiding walk,” Journal of
Physics A: Mathematical and Theoretical 42 323001 (2009) (https:
//iopscience.iop.org/article/10.1088/1751-8113/42/32/323001)
[4] M. N. Rosenbluth and A. W. Rosenbluth, “Monte Carlo Calculation of the
AverageExtension of Molecular Chains,” Journal of Chemical Physics 23 356 (1955)
(https://doi.org/10.1063/1.1741967)
[5] T. Prellberg, “From Rosenbluth Sampling to PERM - rare event sampling with stochastic
growth algorithms,” in R. Leidl and A. K. Hartmann (eds), Modern Computational Science
12: Lecture Notes from the 4th International Oldenburg Summer School, pages 311-334,
BIS-Verlag der Carl von Ossietzky Universita¨t Oldenburg, 2012
(http://www.maths.qmul.ac.uk/˜tp/papers/pub084pre.pdf)
End of Paper.
5