蒲丰投针实验模拟

蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。谢谢。（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project->settings->link->

output，将stack allocations reserve:设为1000000000） 󰀀

program getpi

implicit none

real,parameter::a=5,L=4,pi=3.14159

integer::n1,i,counter=0

real,allocatable::R1(:),R2(:)

real::theta,x,pi1

write(*,*) 'input the size of the array:'

read(*,*) n1

allocate(R1(n1))

allocate(R2(n1))

call random(n1,R1,R2)

do i=1,n1

x=a*(2*R1(i)-1)

theta=pi*R2(i)

if(abs(x)

end do

pi1=(1.0*n1)/(counter*1.0)*(2.0*L/a)

write(*,*) 'this is PI:'

write(*,*) pi

write(*,"('this is ratio of counter to total number',F10.6)")

1.0

&*counter/n1

stop

end program

subroutine random(n,R1,R2)

学生所在学院：理学院 专业：应用数学 班级：091

姓 名 镇凡迪 学 号 0907010278 实验组 五

实验时间 2012.5.9 指导教师 龚敏庆 成

绩

实验项目名称 随机模拟Buffon投针试验

实验目的及要求：

掌握Buffon随机模拟的原理及算法

实验（或算法）原理：

理论依据：客观世界的某些现象之间存在着某种相似性，因而可以从一种现象出发研究另一种现象。比如在分析一个系统时，可先构造一个与该系统相似的模型，通过在模型上进行实验来研究原模型，这就是模拟。随机系统可以用概率模型来描述并进行实验，称为随机模拟方法。

步骤分为（1）建立恰当模型 （2）设计实验方法 （3）从一个或者多个概率分布中重复生成随机数

（4）分析模拟结果

实验硬件及软件平台： 计算机 MATLAB VC 网络

实验步骤：

复习第四章有关知识，熟悉Buffon投针试验的原理

运用MATLAB，编写相关程序，输入代码

观察实验结果，进行讨论

撰写实验报告

实验内容（包括实验具体内容、算法分析、源代码等等）：

（法一）Buffon投针试验求圆周率：

function y=Buffon(a,l,N) //a为平行线间的距离；l为针长；N为投针次数

M=0;

i=0;

while i<=N

xi=rand()*a/2;

yi=rand()*pi;

if xi<=l/2*sin(yi)

M=M+1; end

i=i+1;

end

phi=2*l*N/(a*M);

y=phi;

end

法二：

function [p0,pm,pival]=BuffonMonteCarlo(d,h,N) %Buffon投针试验

if h>=d %判断针的长度h与相邻平行线间的距离d是否满足指定条件

一、利用Matlab计算机语言验证蒲丰(Buffon)投针试验问题

给定a=10，b=5时，模拟100万次投针实验的Matlab程序如下：

a=10;

b=5;

n=1000000;

p=10; % a为平行线间距，b为针的长度，n为投掷次数，p为有效数字位数

x=unifrnd(0,a/2,[n,1]);

phi=unifrnd(0,pi,[n,1]); % 产生均匀分布的随机数，分别模拟针的中点与最近平行线的距离和针的倾斜角

y=x<0.5*b*sin(phi); m=sum(y); % 计数针与平行线相交的次数

PI=vpa(2*b*n/(a*m),p)

运行结果

PI =

3.138919145

二、利用C++计算机语言编程通过大量重复实验验证以下结论：三个阄， 其

中一个阄内写着“有”字，两个阄内不写字， 三人依次抓取，各人抓到“有”字阄的概率均为1/3。

程序如下：

#include

#include

#include

void main()

{

int n=500000;

int i,a[3]={0};

srand(time(NULL));

for(i=0;i

a[rand()%3]++;

printf(

"共测试%d次，其中有字事件有%d次, 占%.2f%%\n"

"抓到无字事件1有%d次，占%.2f%%\n"

"抓到无字事件2有%d次，占%.2f%%\n"

"抓到无字事件共%d次，占%.2f%%",

n,a[0],a[0]*100.0/n,

a[1],a[1]*100.0/n,

a[2],a[2]*100.0/n,

a[1]+a[2],(a[1]+a[2])*100.0/n

);

return 0;

}

一、蒲丰投针问题

在平面上画有等距离的一些平行线，平行线间的距离为

a(a>0) ，向平面上随机投一长为 l(l

交的概率 p，结果发现π =2*l/(a*p).

二、试验方法

能够采纳 MATLAB软件进行模拟实验，即用 MATLAB编写程

序来进行“蒲丰投针实验” 。

1、 基来源理

因为针投到纸上的时候，有各样不一样方向和地点，但

是，每一次投针时， 其地点和方向都能够由两个量独一确立，

那就是针的中点和偏离水平的角度。

以 x 表示针的中点到近来的一条平行线的距离，β表

示针与平行线的交角。明显有 0<=x<=a/2 ，0<=β <=Pi 。用边

长为 a/2 及 Pi 的长方形表示样本空间。为使针与平行线相

交，一定

x<=l*sin

β * ，知足这个关系的地区面积是从

0 到

Pi

的

l*sin

β 对 β 的 积 分 ， 可 计 算 出 这 个 概 率 值 是

(2l)/(Pi*a)

。只需随机生成

n 对这样的

x 和β，就能够模

拟 n 次的投针实验， 而后统计知足 x<=l*sin β * 的 x 的个数，就能够以为这是订交的次数。而后利用公式求得π值。

2、 MATLAB编程

clear ('n')

clear('a')

clear('x')

clear('f')

clear ('y')

clear ('m')

disp(' 本程序用来进行投针实验的演示， a 代表两线间的宽度，针的长度 l=a/2 ，n 代表实验次数 '); a=input(' 请输入 a：');

n=input(' 请输入 n：');

x=unifrnd(0,a/2,[n,1]);

f=unifrnd(0,pi,[n,1]);

y=x<*a*sin(f);

m=sum(y);

PI=vpa(a*n/(a*m))