# basic monte carlo question: area of circle

86 views (last 30 days)

Show older comments

Question: Use monte carlo method to find area of circle radius 5 inside 7x7 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)

##### 0 Comments

### Accepted Answer

pfb
on 1 May 2015

Edited: pfb
on 1 May 2015

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);

##### 0 Comments

### More Answers (3)

Eric
on 2 May 2015

##### 5 Comments

pfb
on 2 May 2015

Edited: pfb
on 2 May 2015

ok, then you can write the function sum with a loop. After creating the vector i (which you need for the plotting)

hits = 0;

for j = 1:N

hits=hits+i(j);

end

misses= N-hits;

If you're not allowed to use logical indexing in the plot you can pretend you do not know it and do this

% vectors for the hits

xh = zeros(1,hits);

yh = xh;

% vector for the misses

xm = zeros(1,misses);

ym = xm;

% counters

mc=1;

hc=1;

for j=1:N

if(i(j))

% put the point in the hits vector

xh(hc)=x(j);

yh(hc)=y(j);

% increment the hit counter

hc=hc+1;

else

% put the point in the misses vector

xm(mc)=x(j);

ym(mc)=y(j);

% increment the miss counter

mc=mc+1;

end

end

plot(xh,yh,'g.');

hold

plot(xm,ym,'r.');

etcetera

munesh pal
on 27 Feb 2018

Calculate the area of the circle using Monte Carlo Simulation

clc

clear all

close all

%%input circle radius and coordinate

r=2;

c_x=7;

c_y=7;

%%position of the outer rectangle

p_x=c_x-r;

p_y=c_y-r;

pos=[p_x,p_y,r^2,r^2]

%%random number generation within the sqaure

N=1000;

a=p_x;

b=p_x+(2*r);

x = a + (b-a).*rand(N,1);

y = a + (b-a).*rand(N,1);

radii = sqrt((x-c_x).^2+(y-c_y).^2);

i = radii <= r;

hits = sum(i);

misses = N-hits;

plot(x(i),y(i),'.g');

hold;

plot(x(~i),y(~i),'.r');

rectangle('Position',pos,'Curvature',[1 1])

rectangle('Position',pos,'EdgeColor','r')

actual_a=22/7*r^2;

ttl = sprintf('Estimate Actual Area of circle: %1.3f, Area of the circle: %1.3f',actual_a,hits/N*(2*r)^2);

title(ttl);

axis equal

##### 1 Comment

hosein bashi
on 30 Jul 2018

Edited: hosein bashi
on 30 Jul 2018

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!