Matlab仿真两种方法求圆周率π

本文用到的matlab源码下载地址

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变大蒙特卡洛演示结果如下：