MatLab homework - Monte Carlo

Simulation on the Nanoscale EGNM03

Lab. 4

In this lab we will develop a suite of Monte Carlo routines to enable us to evaluate

integrals with complex domains. These methods can easily be extended to N

dimensions. First we will use an intuitive graphical/numerical method to determine an

estimate of π that utilises the main components of a Monte Carlo routine.

1. Estimation of π with rejection sampling

A circle with unit radius, r = 1, is inscribed within a square of side, l = 2. The area of

the circle is AC =π r 2 =π 1

2 =π , and the area of the square is A S=( 2 r )

2 =2

2 =4 .

The ratio of the area of the circle to the area of the square is

AC AS

= π r 2

( 2 r )2 =

π 4

≈0 .7853981633974483 .. . [1]

By computing this ratio, we can easily obtain an estimate for π by multiplying the

ratio by 4. One particularly simple way do this is to randomly pick points, (x, y), in

the square and count how many of them lie inside the circle, NC. Dividing this number

by the total number of trials, NT, we obtain an estimate of ρ, where

ρ= N C N T

≈π / 4 [2]

Thus, we can estimate the value of π:

AC≈ ρ∗ AS =ρ∗4≈π [3]

The criterion and ML code for determining if a random point is within the unit circle

is simply:

x=rand(1)*2-1; y=rand(1)*2-1; R=sqrt(x^2+y^2); if R<=1 count=count+1; end

In ML, the function rand generates a random number between 0 and 1. We need to

generate random values between -1 and 1. To do this we multiply rand by 2 to get a

random number within the interval [0 2]. If we simply subtract 1 from this interval we

will obtain random numbers in the interval [-1 1] as desired.

(a) Using the pseudo-code above, develop a Monte Carlo routine to estimate the value

of π for 10, 100, 1000, 10000 iterations.

(b) Estimate the accuracy of the above computation with the following:

dA≈4∗r 2∗√ ρ− ρ 2

n . [4]

2. Evaluation of Integrals

As discussed in the lectures, Monte Carlo methods can be used to approximate

integralsthe area, A, under a curve y=f ( x ) for a≤x≤b . First we must define a

rectangular box R containing A as follows:

R= {( x,y ) : a≤x≤b and 0≤ y≤d } [4]

where

d= max a≤ x ≤b f ( x ) [5]

Equation 4 simply states the domain of R and the allowed values of x and y in those

domains and equation 5 states that the limit d is simply the maximum value of f(x)

within the allow range of x.

Secondly, we must randomly pick points within R, i.e. {( xi ,yi )}i= 1 n

, xi and yi are

chosen from independent uniformly distributed random variables over [a b] and [0 d]

respectively.

Third calculate the ratio ρ as follows:

ρ= m n

[6]

Where m is the number of points that lie in A. The area is then computed using the

approximation:

A≈ ρ∗R=ρ∗( b−a)∗( d−0) =ρ∗( b−a )∗d . [7]

An ‘estimate’ for the accuracy of the above computation is

dA≈( b−a )∗d∗√ ρ− ρ 2

n . [8]

(a) Write a subroutine to calculate the area below a function, f(x), with the aid of the

following pseudo code:

function [rho,ABC,dABC]=MC_Int_Eval_2(f,a,b,c,n) % Input: f function to be evaluated % a,b limits of x domain % c lower limit of y domain % n number of samples % Output: rho equation 6 % ABC equation 7 % dABC equation 8

% Initialise Counters, m in equation 6 which counts the number of % points below the function f(x) & k which counts the number of % points above the function

% Determine the unknown limit, d, by calculating the maximum value of % y, i.e. equation 5, using:

x=linspace(a,b,100)'; y=feval(f,x);

% and the built-in function max; max(x) determines the maximum value % in the array x

% Determine fraction above/below function, i.e. % Loop from 1:n % calculate random variables x & y using

x=a+(b-a)*rand(1); y=d*rand(1);

% Check if y is less than or equal to or greater than f(x) % Use an if else statement to achieve this % Update the appropriate counter, i.e. m = m + 1; or k = k + 1;

% Calculate the sample ratio, equation 6

% Calculate Domain area, R, i.e.

BoxArea=(b-a)*(d-c);

% Calculate the Area Below the Curve, i.e. ABC, equation 7

% Calculate the accuracy in ABC, i.e. dABC equation 8

(b) Define the following functions:

(i) f MC 1( x )=√ 4−x 2

(ii) f MC 2( x )= 4

1+x 2

(iii) f MC 3( x )= 6

√ 4− x2

(c) Create a script to evaluate the integral, I=∫ a

b

f ( x ) dx for the functions defined in

part (b) using the following respective integration limits:

(i) a = 0, b = 2, c = 0;

(ii) a = 0, b = 1, c = 0;

(iii) a = 0, b = 1, c = 0;

For each function determine an estimate of π and an estimate of the accuracy of the

result for the following sample values n = 100, 1000, 10000.