Calculate pi using monte-Carlo simulation with logical vector
9 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
KSSV
le 7 Mar 2017
Modifié(e) : KSSV
le 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 commentaires
Plus de réponses (2)
David Contreras
le 21 Oct 2020
Modifié(e) : David Contreras
le 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 commentaires
Fuat Yilmaz
le 18 Déc 2019
hi can you make a animation for this ?:)
1 commentaire
Rajan Prasad
le 15 Mai 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
Voir également
Catégories
En savoir plus sur Data Distribution Plots dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!