蒲丰氏投针问题的模拟过程

合集下载

蒲丰投针――MonteCarlo算法

蒲丰投针――MonteCarlo算法

蒲丰投针 ―― Monte Carlo 算法背景:蒙特卡罗方法(Monte Carlo ),也称统计模拟方法,是在二次世界大战期间随着科学技术的发展和电子计算机的发明,而被提出的一种以概率统计理论为基础的一类非常重要的数值计算方法。

蒙特卡罗方法在应用物理、原子能、固体物理、化学、生态学、社会学以及经济行为等领域中得到广泛利用。

蒙特卡罗方法的名字来源于世界著名的赌城 —— 摩纳哥的蒙特卡罗。

其历史起源可追溯到1777年法国科学家蒲丰提出的一种计算圆周的方法 —— 随机投针法,即著名的蒲丰投针问题。

问题:设在平面上有一组平行线,间距为d ,把一根长L 的针随机投上去,则这根针和平行线相交的概率是多少?(其中 L < d )分析:由于 L < d ,所以这根针至多只能与一条平行线相交。

设针的中点与最近的平行线之间的距离为 y ,针与平行线的夹角为 θ (0 ≤ θ ≤ π)。

相交情形 不相交情形易知针与平行线相交的充要条件是:sin 2Ly x θ≤=由于1[0,], [0, ]2y d θπ∈∈,且它们的取值均满足平均分布。

建立直角坐标系,则针与平行线的相交条件在坐标系下就是曲线所围成的曲边梯形区域(见右图)。

所以有几何概率可知针与平行线相交的概率是sin d 2212LL p d d πθθππ==⎰Monte Carlo 方法:随机产生满足平均分布的 y 和 θ,其中1[0,], [0, ]2y d θπ∈∈,判断 y 是否在曲边梯形内。

重复上述试验,并统计 y 在曲边梯形内的次数 m ,其与试验次数 n 的比值即为针与平行线相交的概率的近似值。

clear;n = 100000; L = 1; d = 2; m = 0;for k = 1 : ntheta = rand(1)*pi; y = rand(1)*d/2;if y < sin(theta)*L/2m = m + 1; end endfprintf('针与平行线相交的概率大约为 %f\n', m/n)计算π的近似值利用该方法可以计算 π 的近似值:sin d 22 22 1n LL m p d m d L d n πθθπππ⇒≈==≈⎰下面是一些通过蒲丰投针实验计算出来的 π 的近似值:蒲丰投针问题的重要性并非是为了求得比其它方法更精确的π值,而是在于它是第一个用几何形式表达概率问题的例子。

蒲丰投针试验讲解课件

蒲丰投针试验讲解课件

该试验不仅在理论上具有重要意义,对 于理解随机性和几何规律的本质有重要 贡献,而且在实际应用中也有广泛的应
用价值。
蒲丰投针试验可以应用于统计学、物理 学、计算机科学等多个领域,为相关领
域的研究提供了重要的启示和工具。
蒲丰投针试验的局限性
01
02
03
04
蒲丰投针试验虽然是一个经典 的试验,但是它也存在一些局
针方向与平行线垂直。
重复投掷蒲丰投针N次,记录每 次投掷的结果。
测量与计算阶段
测量投掷后蒲丰投针 与平行线之间的距离 ,记录下来。
根据公式π=2*n/N ,计算π的近似值, 其中n为相交次数, N为投掷次数。
根据记录的数据,计 算蒲丰投针与平行线 相交的次数。
CHAPTER 03
试验结果分析
蒲丰投针试验的预期结果
蒲丰投针试验是一种估算π值的方法,其预期结果是通过投掷 一根针到一张白纸上,然后统计针与白纸边缘相交的次数, 来估算π的值。
蒲丰投针试验的预期结果是根据概率论和几何学原理推导出 来的,即当投掷次数足够多时,针与白纸边缘相交的频率接 近于π/4。
实际结果与预期结果的比较
在实际进行蒲丰投针试验时,需要记录针与白纸边缘相交的次数,并计 算出相应的π值。
限性。
首先,该试验的结果受到投针 方式、试验环境等因素的影响 ,可能导致结果存在误差。
其次,蒲丰投针试验的应用范 围相对有限,主要适用于一些 特定的几何形状和随机性问题

最后,蒲丰投针试验的结论仅 适用于理想化的模型,与实际
情况可能存在差异。
未来研究方向与展望
随着科学技术的发展和研究的深入, 蒲丰投针试验在未来仍有广阔的研究 前景。
蒲丰投针试验讲解课 件

投针实验详解

投针实验详解

一、 问题的提出在人类数学文化史中,对圆周率π精确值的追求吸引了许多学者的研究兴趣。

在众多的圆周率计算方法中,最为奇妙的是法国物理学家布丰(Boffon )在1777年提出的“投针实验”。

与传统的“割圆术”等几何计算方法不同的是,“投针实验”是利用概率统计的方法计算圆周率的值,进而为圆周率计算开辟了新的研究途径,也使其成为概率论中很有影响力的一个实验。

本节我们将借助于MATLAB 仿真软件,对“投针实验”进行系统仿真,以此来研究类比的系统建模方法和离散事件系统仿真。

二、 系统建模“投针实验”的具体做法是:在一个水平面上画上一些平行线,使它们相邻两条直线之间的距离都为a ;然后把一枚长为l (0<l <a )的均匀钢针随意抛到这一平面上。

投针的结果将会有两种,一种是针与这组平行线中的一条直线相交,一种是不相交。

设n 为投针总次数,k 为相交次数,如果投针次数足够多,就会发现公式2ln ak计算出来的值就是圆周率π。

当然计算精度与投针次数有关,一般情况下投针次数要到成千上万次,才能有较好的计算精度。

有兴趣的读者可以耐心地做一下这个实验。

为了能够快速的得到实验结果,我们可以通过编写计算机程序来模拟这个实验,即进行系统仿真。

所谓的系统仿真是指以计算机为工具,对具有不确定性因素的、可模型化的系统的一种研究方法。

建立能够反映实验情况的数学模型是系统仿真的基础。

系统建模中需解决两个问题,一个是如何模拟钢针的投掷结果,另一个是如何判断钢针与平行线的位置关系。

这里,设O 为钢针中点,y 为O 点与最近平行线之间的距离,θ为钢针与平行线之间的夹角(0180θ≤< )。

首先,由于人的投掷动作是随机的,钢针落下后的具体位置也是随机的,因此可用按照均匀分布的两个随机变量y 和θ来模拟钢针投掷结果。

其次,人工实验时可以用眼睛直接判断出钢针是否与平行线相交,而计算机仿真实验则需要用数学的方法来判别。

如下图所示,如果y 、l 和θ满足关系式1sin 2y l θ≤,那么钢针就与平行线相交,否则反之,进而可以判断钢针与平行线的位置关系。

蒲丰投针实验原理

蒲丰投针实验原理

蒲丰投针实验原理
蒲丰投针实验是一种检测泥沙粒径分布的实验方法,它是利用悬浮在水中的粒度分布模拟藉由空气流抛掷及落入平板上的控制情形来模拟河流中悬浮颗粒的粒径分布,从而进行检测的。

该实验流程是:将检测的粒料悬浮于水中,利用抛掷及落入平板上的控制条件来模拟河流中悬浮颗粒的粒径分布,然后借助投针实验来观测平面上粒料的分布情况。

最后,根据获得的结果计算出每种粒径的百分率,从而可以得出泥沙粒径分布情况。

蒲丰投针最简单的代码

蒲丰投针最简单的代码

蒲丰投针最简单的代码
蒲丰投针是一种概率统计实验,可以用来求圆周率。

这里介绍一下最简单的蒲丰投针代码。

首先,需要导入Python中的random模块来生成随机数。

然后,定义需要用到的变量和常数,如针长(L)和两根地板板缝之间的距离(d)。

接下来,生成两个随机数,分别表示针的中心点距离地板板缝的距离(x)和针与竖直方向的夹角(theta)。

利用这两个随机数可以计算出针与地板板缝相交的情况。

再用一个计数器变量count来记录针与地板板缝相交的次数,重复这个实验若干次后,圆周率的近似值就可以通过下面的公式计算出来:
pi = 2 * L / (d * p)
其中,p为相交次数与总次数之比。

代码如下:
import random
L = 1 # 针长
d = 2 # 地板板缝间距
n = 10000000 # 实验次数
count = 0 # 相交次数
for i in range(n):
x = random.uniform(0, d) # 针中心距地板板缝距离
theta = random.uniform(0, 180) # 针与竖直方向的夹角
if x <= L * 0.5 * math.sin(theta / 180 * math.pi): # 判断是否相交
count += 1
p = count / n # 相交次数与总次数之比
pi = 2 * L / (d * p) # 计算圆周率
print(pi)
需要注意的是,模拟次数越多,计算出的圆周率越接近真实值。

但是过多的模拟次数会导致程序运行时间增长,因此需要根据实际情况来选择合适的实验次数。

蒲丰(Buffon)投针试验

蒲丰(Buffon)投针试验

一、利用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<stdio.h>#include<stdlib.h>#include<time.h>void main(){int n=500000;int i,a[3]={0};srand(time(NULL));for(i=0;i<n;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(lva)的针,针与平行线相交的概率p,结果发现Ji =2*l/(a*p)・二、试验方法可以采用MATLAB软件进行模拟实验,即用MATLAB编写程序来进行“蒲丰投针实验”。

1、基本原理由丁•针投到纸上的时候,有各种不同方向和位置,但是,每一次投针时,其位置和方向都可以由两个量唯一确定, 那就是针的小点和偏离水平的角度。

以x表示针的中点到最近的一条平行线的距离,B表示针与平行线的交角。

显然有0<=x<=a/2, 0<=f3<=Pi o用边长为a/2及Pi 的长方形表示样本空间。

为使针与平行线相交, 必须x<=l*sin 3 *0.5,满足这个关系的区域面积是从0到Pi 的严sin B对B的积分,可计算出这个概率值是(21)/(Pi*a)o只要随机牛成n对这样的x和P,就可以模拟n次的投针实验, 然后统计满足x<=l*sin B *0.5的x的个数,就可以认为这是相交的次数。

然后利用公式求得兀值。

2、MATLAB 编程clear ('n')clearCa')clear('x')clear(T)clear ('y‘)clear ('m')dis"本程序用来进行投针实验的演示」代衣两线间的宽度, 针的长度l=a/2, n代表实验次数a=input(f iH 输入a: *);n=input(*W输入n: *);x=unifrnd(0,a/2,[n, 1 ]);f=unifmd(O,p i,[n,l]);y=x<0.25*a*sin(f);m=sum(y);PI=vpa(a* n/(a* m))三、实验数据(部分程序截屏见后)四、实验结论从上述数据分析可知,随着模拟次数的越来越多,PI的值逐渐稳定在H值附近,即越来越趋近于匚,故蒲丰投针实验确实可以模拟出兀的值。

蒲丰投针概率推导过程

蒲丰投针概率推导过程

蒲丰投针概率推导过程蒲丰投针是一个经典的概率问题,它的推导过程非常有趣。

本文将从一个人的视角进行叙述,以增加读者的情感共鸣。

我要向大家介绍一下蒲丰投针的背景和问题。

蒲丰投针是由法国数学家蒲丰在18世纪提出的,他提出了一个问题:如果有一根长度为L的针,将其随机地投在一块有平行线的地板上,那么这根针与平行线相交的概率是多少?现在,让我们来推导一下这个概率。

首先,我们需要假设一些条件。

假设针的长度L小于等于平行线的间距D,且针的中点与平行线之间的距离为x(0<=x<=D/2)。

这样,我们可以将针的位置分为两种情况:一种是针与平行线相交,另一种是针与平行线不相交。

让我们来计算针与平行线不相交的情况。

在这种情况下,针的中点到最近的平行线的距离大于针的半径。

针的半径可以用L/2来表示,所以针的中点到最近的平行线的距离大于L/2。

我们可以将这个问题转化为一个几何问题,即计算一个长度为L/2的线段与两条平行线之间不相交的概率。

假设针的中点到最近的平行线的距离为y(L/2<=y<=D/2),那么不相交的情况可以表示为y >= L/2。

由于y的取值范围是L/2到D/2,所以不相交的概率可以表示为 (D/2 - L/2)/(D/2)。

接下来,让我们来计算针与平行线相交的情况。

在这种情况下,针的中点到最近的平行线的距离小于等于针的半径。

也就是说,针的中点到最近的平行线的距离小于等于L/2。

我们可以将这个问题转化为一个几何问题,即计算一个长度为L/2的线段与两条平行线之间相交的概率。

假设针的中点到最近的平行线的距离为y(0<=y<=L/2),那么相交的情况可以表示为y <= L/2。

由于y的取值范围是0到L/2,所以相交的概率可以表示为 (L/2)/L/2。

现在,我们可以将不相交的概率和相交的概率相加,得到针与平行线相交的概率。

即 P = (D/2 - L/2)/(D/2) + (L/2)/L/2。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

蒲丰氏投针问题的模拟过程,随机数发生器也是自编的,以供大家参考和提出建议。

谢谢。

(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)<L*sin(theta)) counter=counter+1
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)
implicit none
integer n,seed1,seed2,i,little,m
real::R1(n),R2(n)
integer::temp1(n),temp2(n)
write(*,*) 'input seed1'
read(*,*) seed1
write(*,*) 'input seed2'
read(*,*) seed2
m=2**30
m=m*2-1
temp1(1)=397204094
temp2(1)=temp1(1)
R1(1)=(1.0*temp1(1))/(1.0*m)
R2(1)=(1.0*temp2(1))/(1.0*m)
little=0
if(R1(1)<0.5) little=little+1
do i=1,n-1
temp1(i+1)=mod(seed1*temp1(i),m)
R1(i+1)=(1.0*temp1(i+1))/(1.0*m)
temp2(i+1)=mod(seed2*temp2(i),m)
R2(i+1)=(1.0*temp2(i+1))/(1.0*m)
if(R1(i+1)<0.5) little=little+1 .
end do ;
write(*,*) 'ratio of number which is little than 0.5'
write(*,*) 1.0*little/n
return
end subroutine
拟蒙特卡洛方法(Quasi-Monte Carlo)积分实例
%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分 format long; syms x
a = sym(1/2);
f = exp(-a*x^2);
ezplot(f)
disp(int(f,-1,1));
fprintf('integral result:%1.18f.\n',double(int(f,0,1)));
%disp(double(int(f,0,1)));
复制代码%使用拟蒙特卡洛方法积分
%得到拟蒙特卡洛序列,即低偏差序列,halton法
%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset 函数实现,
x=halton(10000,2,5577);
n=length(x);
mju=0;
for i=1:n
mju=mju + exp(-0.5*x(i)^2);
end
mju=mju/n;
fprintf('Quasi-Monte Carlo result:%1.18f.\n',mju);
%disp(mju);
%使用蒙特卡洛方法积分
%得到Uniform序列,
x=random('unif',0,1,10000,1);
n=length(x);
mju=0;
for i=1:n
mju=mju + exp(-0.5*x(i)^2);
end
mju=mju/n;
fprintf('Monte Carlo result:%1.18f.\n',mju);
%=============生成HALTON序列======================== function result = halton( m,base,seeder )
%生成HALTON序列
% Check inputs
if nargin < 3
seeder = 0;
if nargin < 2
error('MATLAB:Halton:NotEnoughInputs',...
'Not enough input arguments. See Halton.'); end
end
res=0;
n=length(base);
for i=1:m
for j=1:n
element=0;
temp=seeder+i;
k=1;
while temp>0
element(k)=rem(temp,base(j));
temp=fix(temp/base(j));
k=k+1;
end
res(i,j)= 0;
for k=1:length(element)
res(i,j)=res(i,j)+element(k)/(base(j)^k); end
end
end
result=res;。

相关文档
最新文档