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 − 6×2 ) 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) = 5×4, 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) = 4×3, f2(x) = 3×2 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
