目录:
1.1 蒲丰投针法求π
1.1.1 蒲丰投针算法介绍
1.1.2 混合同余算法介绍
1.1.3 程序及注释
混合同余法产生随机数——Mixrand()函数文件Mixrand.m
function [q] = Mixrand(n)
%A是取余函数(Ax+n)的参变量
a=1;A=4*a+1;
%循环实验
x=1;
%当满足
%x=rem(A*x+n,2^k)
%A=4*a+1
%y=x/2^k
%T>2^m(优化后的混合同举法)
for i = 1:n
x=rem(A*x+i,2^20);
d=x/(2^20);
X(i)= d;
A=4*a+1;%该值舍去(提高精度)
%进行两次该算法的原因:因为发现y的离散性更好,故猜测多次迭代可提高精度
x=rem(A*x+i,2^20);
d=x/(2^20);
X(i)= d;
a=a+1;
x=rem(A*x+i,2^20);
d=x/2^20;
Y(i)= d;
a=a+1;
end
%第二次混合同余法,取得更精准的数据求平均
for i = 1:n
x=rem(A*x+i,2^20);
d=x/(2^20);
X(i)= d;
A=4*a+1;%该值舍去(提高精度)
%进行两次该算法的原因:因为发现y的离散性更好,故猜测多次迭代可提高精度
x=rem(A*x+i,2^20);
d=x/(2^20);
X(i)= d;
a=a+1;
x=rem(A*x+i,2^20);
d=x/2^20;
Y(i)= d;
a=a+1;
end
q=Y;
end
模拟投针过程文件.test6.m
clear all
d=1;%两条边的距离
l=0.6;%针的长度
%总的实验次数
n = input('请输入n:');
N=n;
m=0;
mpi=0;
x=Mixrand(3*N);
for i=1:n
x1(i)=x(i);
y1(i)=x(i+N);
Mixpi(i)= pi*x(i+2*N);%mpi为接收到的(0-2pi/)的随机数
if x1(i)<l*sin(Mixpi(i))/2
m=m+1;
end
end
rectangle('Position',[-0.6 -0.6 2.2 2.2],'EdgeColor','r','LineWidth',3);%画矩形
hold on;
rectangle('Position',[0 0 1 1],'EdgeColor','m','LineWidth',7);
axis([-0.7 1.7 -0.7 1.7]);
for i=1:n
P(i)=cos(2*Mixpi(i));
Q(i)=sin(2*Mixpi(i));
x11(i)=x1(i)+l*P(i);
y11(i)=y1(i)+l*Q(i);
hold on;
end
plot(x1,y1,'X');
title('模拟蒲丰投针');
xlabel('X轴');
ylabel('Y轴');
X=[x1;x11];
Y=[y1;y11]
line(X,Y);
q=l*n/(d*m);
主函数文件tz.m
d=1;%两条边的距离
%总的实验次数
N=n;
m=0;
mpi=0;
x=Mixrand(3*n);
for i=1:n
x1(i)=x(i);
y1(i)=x(i+N);
Mixpi(i)= pi*x(i+2*N);%mpi为接收到的(0-2pi/)的随机数
if x1(i)<l*sin(Mixpi(i))/2
m=m+1;
end
end
q=l*n/(d*m);
变量分析文件anlgnse.m
clear all;
t=10^4:10^4:1*10^6;
M=length(t);
for i=1:length(t) %采集Pi的个数
m(i)=tz(t(i),0.6); %1000-10000中选择
end
plot(m,'X');
hold on;
for i=1:M
x(i)=i;
y(i)=3.1415926;
end
plot(x,y,'g');
title('Pi随着N变化的变化图');
xlabel('N点');
ylabel('Pi值');
分析l/a(针长/边长)对实验结果的影响文件anlgnse1.m
clear all;
%t=10^4:10^4:1*10^6;
I=0.01:0.01:1
M=length(I);
for i=1:length(I) %采集Pi的个数
m(i)=tz(100000,I(i)); %1000-10000中选择
end
for i=1:M
x(i)=i;
y(i)=pi;
end
plot(I,m,'X');
hold on;
axis([0 1 3 3.3]);
title('Pi随着l/a变化的变化图');
xlabel('l/a');
ylabel('Pi值');
1.1.4 仿真结果
1.2 蒙特卡洛法求π
1.2.1 蒙特卡洛算法介绍
1.2.2 关键程序
算法2——蒙特卡洛方法:
%每次重置所有变量
clear all;
m = 0;
n=input('请输入n:');
%使用的圆的半径
r = 1;
%A是取余函数(Ax+n)的参变量
a=1;A=4*a+1;
%循环实验
x=1;
d1=Mixrand(2*n);
N=n;
for i = 1:n
X(i)= d1(i);
Y(i)= d1(i+N);
if (X(i)^2 + Y(i)^2 <= r^2)
X1(i)=X(i);
Y1(i)=Y(i);
m = m + 1;
else
X2(i)=X(i);
Y2(i)=Y(i);
end
end
p1=4 * m / n;
%显示结果
X0=0;Y0=0;
t=0:pi/40:2*pi;
r=1;
XX=abs(X0+r*cos(t));
YY=Y0+r*sin(t);
plot(X1,Y1,'.',X2,Y2,'.',XX,YY,'-');
fprintf('当总实验次数n = %d时计算出来的圆周率:Pi = %d\n',n, p1);
anlgnse函数
clear all;
t=5*10^4:10^4:1*10^5;
M=length(t);
for i=1:length(t) %采集Pi的个数
m(i)=mt(t(i)); %1000-10000中选择
end
plot(m,'X');
hold on;
for i=1:M
x(i)=i;
y(i)=3.1415926;
end
plot(x,y,'g');
title('Pi随着N变化的变化图');
xlabel('N点');
ylabel('Pi值');
1.2.3 仿真结果
当投针N变大蒙特卡洛演示结果如下: