MATLAB: Basic monte carlo question: area of circle

monte carlo

Question: Use monte carlo method to find area of circle radius 5 inside 7×7 square.
So this is the code I've put together so far to determine how many points land inside or on the circle (hits). My output for hits is 0 but I can't for the life of me figure out why, to me everything seems fine but clearly isn't.
clear;
N= 1000; % number of points generated
a = -7;
b = 7;
hits= 0;
for k = 1:N
     x = a + (b-a).*rand(N,1);
     y = a + (b-a).*rand(N,1);
     % Count it if it is in/on the circle
     radii = sqrt(x.^2+y.^2);
         if radii <=5;
            hits = hits + 1;
         end
end
disp(hits)

Best Answer

You are mixing things up. Your code is part vectorized, part scalar.
You have a loop over N, but then at each iteration you generate N points (so you're generating N^2 points overall). The variables x, y, and radii are vectors.
Note that
 [6 7 4 2 8 5] <=5
gives
 [0 0 1 1 0 1]
Therefore condition radii<=5 has a very low probability to be met ( all of the N radii should be less than 5).
As a result hits is (almost) never accumulated.
Your code is structured for a scalar radius. It would work if you substitute the lines
 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);
with
 x = a + (b-a).*rand;
 y = a + (b-a).*rand;
In this case hits will be a number between 0 and N, and hits/N is going to be proportional to the ratio of the areas of the circle and square, pi*(5/14)^2 (note that you can use simple power operator ^ instead of elementwise .^)
You can avoid the loop altogether using the vectorized variables.
 N= 1000; % number of points generated
 a = -7;
 b = 7;
 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);
 radii = sqrt(x.^2+y.^2);
 hits = sum(radii<=5);
Related Question