Calculate pi using monte-Carlo simulation with logical vector
10 views (last 30 days)
Show older comments
I want to know how to model the script for calculating pi using Monte-Carlo simulation with using logical vectors
I already know the method using for/if, but does not know the method with logical vector
Please let me know how to do it
the below script is for the method using for/if
clear
clc
inside = 0;
nmax = 10000;
for n = 1:nmax
x = rand;
y = rand;
x1=x-0.5;
y1=y-0.5;
if sqrt(x1^2+y1^2) <= 0.5
plot(x1,y1,'b.');
inside = inside + 1;
else
plot(x1,y1,'r.');
end
if ( mod(n,1000) == 0 )
pi = 4 * inside / n;
fprintf('%8.4f\n',pi);
end
hold on;
end
hold off;
pi = 4 * inside / nmax;
fprintf('%8.4f\n',pi);
0 Comments
Accepted Answer
KSSV
on 7 Mar 2017
Edited: KSSV
on 7 Mar 2017
nmax = 10000;
x = rand(nmax,1);
y = rand(nmax,1);
x1=x-0.5;
y1=y-0.5;
r = sqrt(x1.^2+y1.^2) ;
% get logicals
inside = r<=0.5 ;
outside = r>0.5 ;
% plot
plot(x1(inside),y1(inside),'b.');
hold on
plot(x1(outside),y1(outside),'r.');
% get pi value
thepi = 4*sum(inside)/nmax ;
fprintf('%8.4f\n',thepi)
0 Comments
More Answers (2)
David Contreras
on 21 Oct 2020
Edited: David Contreras
on 3 Nov 2020
This might be a little faster by making all random points X and Y on a single function call.
nmax = 10000;
%Make x and y positions for all points
points = rand(nmax,2);
r=0.5;
%Translate to the origin
P1 = points-r;
%Check which points are insisde
inside = P1(:,1).^2+P1(:,2).^2 <= r.^2;
%Separate points in two sets, inside and outisde
IN = points(inside,:);
OUT = points(not(inside),:);
%Calculate PI
PI = 4*mean(inside);
%Plot blue inside pionts and redo outside points
plot(IN(:,1),IN(:,2),'b.',OUT(:,1),OUT(:,2),'r.')
title(['\pi = ', num2str(PI)])
axis([0 1 0 1], "equal")
0 Comments
Fuat Yilmaz
on 18 Dec 2019
hi can you make a animation for this ?:)
1 Comment
Rajan Prasad
on 15 May 2020
Try this for simulation using frame capture
clc;
clear;
%---------% starts from here__________________________
a=char({'\pi'});
n=(0:10:1000)*100;
n(1)=100;
for i=1:length(n)
np=n(i);
x=rand(np,1);
y=rand(np,1);
dis=x.^2+y.^2;
%----------Points separation---------------
in_p = dis<=1;
out_p=dis>1;
%---------------Points plot------------------
plot(x(in_p),y(in_p),'.r');
hold on;
plot(x(out_p),y(out_p),'.b');
%---------------Estimating pi----------------
tpi(i)=4*sum(in_p)/np;
str1=sprintf('n=%d,%s =%.5f',np,a,tpi(i));
title(str1,'FontName','Times New Roman');
clear x y dis in_p out_p
drawnow;
frame = getframe(gcf);
% im(i)=frame;
im2{i} = frame2im(frame);
clf;
end
filename = 'montecarlo.gif'; % Specify the output file name
for i = 1:length(n)
[A1,map] = rgb2ind(im2{i},256);
if i == 1
imwrite(A1,map,filename,'gif','LoopCount',Inf,'DelayTime',0.4);
else
imwrite(A1,map,filename,'gif','WriteMode','append','DelayTime',0.4);
end
end
See Also
Categories
Find more on Data Distribution Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!