# Estimate $\pi$ using the Monte Carlo Method in MATLAB

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$$

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 ?

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


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


• +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. – horchler Jan 13 '17 at 0:27
• @horchler Thanks for the edit! – Luis Mendo Jan 13 '17 at 1:12
• 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.. – user326307 Jan 17 '17 at 15:01
• @user326307 Please see edit – Luis Mendo Jan 17 '17 at 15:20
• @LuisMendo thank you for your help – user326307 Jan 17 '17 at 15:26