Calculate pi using monte-Carlo simulation with logical vector

10 visualizaciones (últimos 30 días)
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);

Respuesta aceptada

KSSV
KSSV el 7 de Mzo. de 2017
Editada: KSSV el 7 de Mzo. de 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)

Más respuestas (2)

David Contreras
David Contreras el 21 de Oct. de 2020
Editada: David Contreras el 3 de Nov. de 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")

Fuat Yilmaz
Fuat Yilmaz el 18 de Dic. de 2019
hi can you make a animation for this ?:)
  1 comentario
Rajan Prasad
Rajan Prasad el 15 de Mayo de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Data Distribution Plots en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by