clc;clear all;
M=1000;n=0;
for j=1:M
N=10000000;a=8;L=3; %必须 a>=L>0
B=2*pi*rand(N,1); %针与X轴正方向所成的角度。
y0=-a+2*a*rand(N,1); %针上某定点在y轴的坐标。
r=2*L*rand; %针上某定点到针一端之长。
R=2*L-r; %针上某定点到针另一端之长。
y1=y0+r*sin(B); %针一端在y轴的坐标。
y2=y0-R*sin(B); %针另一端在y轴的坐标。
for i=1:N
if (y1(i)-a)*(y2(i)-a)<=0|(y1(i)+a)*(y2(i)+a)<=0 %针与两平行线相交的条件。
n=n+1;
end
end
end
disp('圆周率=');
p=(2*M*N*L)/(n*a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%若取r=2*L,(即y0是针最下端在y轴的坐标)程序简化如下:
clc;clear all;
M=1000;n=0;
for j=1:M
a=8;L=3;N=1000000;
B=pi*rand(N,1); %此处用0.5*pi亦可。
y0=-a+2*a*rand(N,1);
for i=1:N
if y0(i)+2*L*sin(B(i))>=a %针与直线y=a相交的条件。
n=n+1;
end
end
end
disp('圆周率=');
p=(2*M*N*L)/(n*a)
文章评价:
优秀
0
好
0
一般
0
差
0
喜欢
0
收藏
0