1
$\begingroup$

Evaluate the area of a circle of radius $1= \pi$ using Monte Carlo method . Hence we can generate pairs of random numbers $(x_i,y_i) \in [-1,1]$.

Thus :

$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} X 4 $$

enter image description here

Consider a unit circle inscribed in a square, each of the small circles drawn on this figure represents a random point that was generated in the square, the red and blue circles represent points inside and outside the unit circle respectively (I got the figure from google) . If we choose a point uniformaly at random within the square, then the probability that the point lies in the unit circle is

$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} $$

We know that the area of the circumscribed square is $4$, if we knew $p$, then we could compute the area of the unit circle:

Area of the unit circle = $p$ x area of the circumscribed square = $4p$

How to Estimate $\pi$ using the Monte Carlo Method in MATLAB, then what's the error a quarter of a circle ?

$\endgroup$
3

2 Answers 2

2
$\begingroup$

The following code is vectorized, and thus typically faster than using loops:

N = 1e6;                % Number of samples
r = 1;                  % Circle radius
x = 2*rand(1,N)-1;      % N samples between -1 and 1
y = 2*rand(1,N)-1;      % N samples between -1 and 1
r2 = x.^2 + y.^2;       % compute squared distance to origin
area = mean(r2<=r^2)*4; % the area is 4 (area of square) times the proportion
                        % of points in the circle. Note that the threshold
                        % should be r^2 (1^2)
error = pi - area       % error in the estimation of pi

N_all = round(10.^(1:.5:7));      % Numbers of samples
r = 1;                            % Circle radius
errors = NaN(size(N_all));        % Preallocate result
for k = 1:numel(N_all)            % Loop over size of N_all
    N = N_all(k);                 % Current N
    x = 2*rand(1,N)-1;            % N samples between -1 and 1
    y = 2*rand(1,N)-1;            % N samples between -1 and 1
    r2 = x.^2 + y.^2;             % compute squared distance to origin
    area = mean(r2<=r^2)*4;       % the area is 4 (area of square) times the
                                  % proportion of points in the circle
    error_all(k) = pi-area;       % error in the estimation of pi
end                               % End loop
semilogx(N_all, error_all, 'o-'); % Plot with logarithmic x axis
grid                              % Use grid

enter image description here

$\endgroup$
6
  • 1
    $\begingroup$ +1 I'll add that it might be a good idea to set the random number generator seed via rng. Also, r2 could be calculated as r2 = sum((2*rand(2,N)-1).^2);, though that code might be less maintainable/understandable. $\endgroup$
    – horchler
    Jan 13, 2017 at 0:27
  • $\begingroup$ @horchler Thanks for the edit! $\endgroup$
    – Luis Mendo
    Jan 13, 2017 at 1:12
  • $\begingroup$ this code run quikly, i try to run the code with $n=10, 100, 100$ until $10e6$ and it gives different error. Then how to create the code in MATLAB to discribe the value of error.. $\endgroup$
    – user326307
    Jan 17, 2017 at 15:01
  • 1
    $\begingroup$ @user326307 Please see edit $\endgroup$
    – Luis Mendo
    Jan 17, 2017 at 15:20
  • 1
    $\begingroup$ @LuisMendo thank you for your help $\endgroup$
    – user326307
    Jan 17, 2017 at 15:26
1
$\begingroup$

Some code like this may work:

clear;
N=1000; # the experiment event number
r=1; #the circle radius
n=0; # sucessful event number 
for i=1:N
  x=-r+2*r*rand();
  y=-r+2*r*rand();
  if ((x^2+y^2)<=r^2)
    n = n+1;
  end
end

pi_sim=4*n/N
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .