Calculate pi using monte-Carlo simulation with logical vector

조회 수: 10 (최근 30일)
Seungman Kim
Seungman Kim 2017년 3월 7일
편집: David Contreras 2020년 11월 3일
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);

채택된 답변

KSSV
KSSV 2017년 3월 7일
편집: KSSV 2017년 3월 7일
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)

추가 답변 (2개)

David Contreras
David Contreras 2020년 10월 21일
편집: David Contreras 2020년 11월 3일
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 2019년 12월 18일
hi can you make a animation for this ?:)
  댓글 수: 1
Rajan Prasad
Rajan Prasad 2020년 5월 15일
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

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Data Distribution Plots에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by