matlab模拟蒲丰掷针实验计算pi的值

蒲丰投针――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 πθθπππ⇒≈==≈⎰下面是一些通过蒲丰投针实验计算出来的 π 的近似值：蒲丰投针问题的重要性并非是为了求得比其它方法更精确的π值，而是在于它是第一个用几何形式表达概率问题的例子。

投针实验详解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

蒲丰投针试验

平行线的距离
针长
投掷次数
相交次数
π的近似值
1
0.7
1000
462
3.0303
1
0.9
2000
1146
3.1414
2
1.7
3000
1649
3.0928
2
1.9
4000
2401
3.1653
3
2.7
5000
2932
3.0696
3
2.9
6000
3696
3.1385
4
3.8
7000
4238
3.1383
4
3.6
考核结果
教师签名：年月日
模拟法(也称为Monte-Carlo法).
实验的目的与作用
(1)理解频率具有客观稳定性；
(2)理解概率是频率的稳定值；
(3)知道我们常用频率作为计算概率的近似值；
(4)掌握通过设计一个随机实验，使一个事件的概率与某一未知数有关，然后通过重复实验，以频率近似概率，即可求得未知数的方法。
实验方法
运用计算机模拟蒲丰投针实验，通过重复实验，以频率近似概率，并通过公式推导，求的圆周率π的近似值。
实验报告
课程名称：___概率论与数理统计___
学院名称：____数学与统计学院____
班级：__________
姓名：
学号：
2013-2014______学年第 _____2____学期
数学与统计学院制
实验地点ห้องสมุดไป่ตู้
课程类别
①公共课□ ②专业课□
实验日期
2014.4.28
实验编组

matlab function实现pi运算

matlab function实现pi运算如何用Matlab实现pi的计算在科学计算中，π（pi）是一个非常重要的数学常数。

它代表了一个圆的周长与其直径之间的比例关系。

不同的领域，如数学、物理、工程等，经常需要使用到π。

Matlab是一款功能强大的科学计算软件，我们可以利用其编程能力来实现π的计算。

本文将一步一步地介绍如何使用Matlab来编写一个计算π的函数。

步骤一：了解π的计算方法计算π的方法有很多种，其中一种著名的方法是蒙特卡洛方法。

该方法基于随机抽样，利用圆的特性进行估算。

我们首先在一个正方形内随机生成大量的点，然后计算这些点落在圆内的比例。

通过该比例乘以正方形的面积，我们就可以估算出圆的面积，从而得到π的近似值。

步骤二：编写Matlab函数首先，我们需要定义一个函数来实现π的计算。

打开Matlab编辑器，并创建一个新的函数文件pi_calculation.m。

接下来，我们将在函数内部编写代码。

步骤三：生成随机点我们可以使用Matlab的rand函数来生成随机的x和y坐标。

由于我们只需要在正方形内生成点，所以我们可以将随机数限制在[0,1]之间。

为了方便起见，我们可以一次性生成大量的点，而不是逐个生成。

matlabfunction pi_value = pi_calculation(num_points)% Generate random pointspoints = rand(num_points, 2);end步骤四：计算落在圆内的点接下来，我们需要判断这些随机生成的点是否落在圆内。

根据圆的特性，我们可以通过判断点到原点的距离是否小于半径来判断点是否在圆内。

matlab% Count points inside the circlecount_inside = sum(points(:, 1).^2 + points(:, 2).^2 < 1);步骤五：计算π的近似值通过计算落在圆内的点的个数，我们可以得到圆的面积的近似值。

计算pi

一、实验目的探索精确计算π值的方法，并且比较不同方法之间的不同之处和优缺点。

掌握数值积分的辛普森公式。

二、问题描述1. 任务11) 用反正切函数的幂级数展开式结合有关公式求π，若要精确到40位、50位数字，试比较简单公式和Machin 公式所用的项数。

2) 验证公式111=arctan arctan arctan 4258π++ 试试此公式右端做幂级数展开完成任务1所需要的步数。

2. 任务2用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson 法精确到15位数字。

3. 任务3用Monte Carlo 法计算π，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得5位精确数字？设计方案用计算机模拟Buffon 实验4. 任务4利用积分20(1)!!sin !!2n n xdx n ππ-=⎰ ，n 为奇数推导公式224422213352121n n n n π=-+ ……… 用此公式计算π的近似值，效果如何？5. 任务5利用学过的知识（或查阅资料），提出其他计算π的方法（先用你学过的知识证明），然后实践这种方法。

对你在实验中应用的计算π的方法进行比较分析。

6. 任务6e 是一个重要的超越数1e lim 1)n n n→∞=+（ 1111...2!!(1)!e e n n θ=++++++ 试用上述公式或其他方法近似计算e 。

三、问题解法1. 任务11) 根据幂级数展开的相关知识，易知：24122211(1)1n n x x x x--=-+-+-++……… 因为21(arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 35211arctan (1)3521n n x x x x x n --=-+-+-+-……… 当x=1时，-11111--(-1)4352-1n n π=+⋯⋯++⋯ 当叠加了十万次以后得到结果π=3.141582654…只有五位有效数字，可见其精度与效率极低。

Π的计算

π的计算一、实验指导书解读及目的在本次试验中，我们将追溯关于圆周率的计算历程。

通过古典方法，对数值积分法，级数加速法、蒙特卡洛法等计算方法的介绍和计算体验，感受数学思想和数学方法的发展过程，提高对极限和级数收敛性及收敛速度的综合认识，同时使我们看到数学家对科学真理的永无止境的追求。

实验目的：1.用多种方法计算圆周率错误！未找到引用源。

的值；2.通过实验来说明各种方法的优劣；3.尝试提出新的计算方法。

π的简介圆周率是一个常数（约等于 3.1415926），是代表圆周长和直径的比例。

它是一个无理数，即是一个无限不循环小数。

但在日常生活中，通常都用3.14来代表圆周率去进行计算，即使是工程师或物理学家要进行较精密的计算，也只取值至小数点后约20位。

π(pai)是第十六个希腊字母，本来它是和圆周率没有关系的，但大数学家欧拉在一七三六年开始，在书信和论文中都用π来代表圆周率。

既然他是大数学家，所以人们也有样学样地用π来表示圆周率了。

但π除了表示圆周率外，也可以用来表示其他事物，在统计学中也能看到它的出现。

二、实验计划、过程与结果1.古典方法：用圆内接正多边形和圆外切正多边形来逼近以阿基米德的圆内接96边形和圆外切96边形逼近为例已知：sin错误！未找到引用源。

<错误！未找到引用源。

<tan错误！未找到引用源。

，另错误！未找到引用源。

=错误！未找到引用源。

/96推出：96 sin<错误！未找到引用源。

< 96tan错误！未找到引用源。

编写matlab程序format longx=sin(pi/96)y=96*x得：96sin=3.141031950890509Format longx=tan(pi/96)y=96*x得：96tan错误！未找到引用源。

=3.1427145996453683.141031950890509<<3.1427145996453682、"投针实验"求圆周率的方法1777年，法国数学家蒲丰取一根针，量出它的长度，然后在纸上画上一组间距相等的平行线，这根针的长度是这些平行线的距离是的一半。

如何计算π的值(MATLAB)

如何计算π的值1、蒙特卡罗(Monte Carlo)法思想:取一正方形A,以A 的一个顶点为圆心,A 的边长为半径画圆,取四分之一圆(正方形内的四分之一圆)为扇形B 。

已知A 的面积,只要求出B 的面积与A 的面积之比B AS k S =,就能得出B S ,再由B 的面积为圆面积的四分之一,利用公式2=S R π圆即可求出π的值。

因此,我们的目的就就是要找出k 的值。

可以把A 与B 瞧成就是由无限多个点组成,而B 内的所有点都在A 内。

随机产生n 个点,若落在B 内的有m 个点(假定A 的边长为1,以扇形圆心为坐标系原点。

则只要使随机产生横纵坐标x 、y 满足221x y +≤的点,就就是落在B 内的点),则可近似得出k 的值,即m k n=,由此就可以求出π的值。

程序(1):i=1;m=0;n=1000;for i=1:na=rand(1,2);if a(1)^2+a(2)^2<=1m=m+1;endendp=vpa(4*m/n,30)程序运行结果:p =3、14000002、泰勒级数法思想:反正切函数arctan x 的泰勒级数展开式为:35211arctan (1)3521k k x x x x x k --=-+-⋅⋅⋅+-+⋅⋅⋅- 将1x =代入上式有1111arctan11(1)43521n n π-==-+-⋅⋅⋅+--、利用这个式子就可以求出π的值了。

程序(2):i=1;n=1000;s=0;for i=1:ns=s+(-1)^(i-1)/(2*i-1);endp=vpa(4*s,30)程序运行结果:p =3、1494126当取n 的值为10000时,就会花费很长时间,而且精度也不就是很高。

原因就是1x =时,arctan1的展开式收敛太慢。

因此就需要找出一个x 使得arctan x 收敛更快。

若取12x =,则我们只有找出α与4π的关系,才能求出π的值。

令1arctan 2α=,4πβα=-, 根据公式tan tan tan()1tan tan αβαβαβ++=-有1tan 3β=,则有11arctan arctan 423π=+。

matlab 圆周率

nl n 2 , 特别取针的长度 l d / 2 时， md m
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
x表示针的中点与最近的一条平行线间的距离, 为针与平行线的夹角.则0 x d 2 ,0 , 确定x 平面内一个矩形 .如图：
平面内一个矩形确定的夹角为针与平行线一条平行线间的距离表示针的中点与最近的针与平行线相交本实验涉及的方法3蒙特卡罗法蒲丰buffon掷针法针与平行线相交时针在平行线矩形区域内所组成的区域面积与矩形区域面积比即为针与直线相交的概率
实验三
圆周率的求法
问
题
• 设半径为1的圆的周长与圆的直径比为A; • 设半径为2的圆的周长与圆的直径比为B. • 请问A与B哪个大?
计算历程的几个阶段
3.分析方法时期 17世纪出现了数学分析，这锐利的工具使得许多初等数学束手无策的问题迎刃而解。 π 的计算历史也随之进入了一个新的阶段。
计算历程的几个阶段
1593年，韦达给出
1671年，苏格兰数学家詹姆斯.格雷戈里公开了他发现的公式： 3 5 7
arctan x x x3 x5 x7 (1 x 1)

故
2 nl md
本实验涉及的方法
3、蒙特卡罗法——蒲丰(Buffon)掷针法
• 随机生成一个针的中点与最近一条平行直线间的距离x,满足0≤x≤d/2; • 随机生成针与平行线的夹角a,满足0≤a≤ ； • 判断x≤L/2*sin[a],若成立，则针与直线相交，否则不相交。
本实验涉及的方法
计算历程的几个阶段
圆周长大于内接正四边形而小于外切正四边形
计算历程的几个阶段
在我国，首先是由数学家刘徽得出较精确的圆周率。他采用割圆术。

Matlab学习系列18-蒙特卡罗方法

Matlab学习系列18-蒙特卡罗⽅法18. 蒙特卡罗⽅法（⼀）概述⼀、原理蒙特卡罗（Monte Carlo ）⽅法，是⼀种基于“随机数”的计算机随机模拟⽅法，通过⼤量随机试验，利⽤概率统计理论解决问题的⼀种数值⽅法。

其理论依据是：⼤数定律、中⼼极限定理（估计误差）。

常⽤来解决如下问题：1. 求某个事件的概率，基于“频率的极限是概率”；2. 可以描述为“随机变量的函数的数学期望”的问题，⽤随机变量若⼲个具体观察值的函数的算术平均值代替。

⼀般形式为求积分：[]()()()d ba E f X f x x x ?=? X 为⾃变量（随机变量），定义域为[a,b], f (x )为被积函数，()x ?为概率密度函数。

蒙特卡罗法步骤为：(1) 依据概率分布()x ?不断⽣成N 个随机数x , 依次记为x 1, …, x N , 并计算f (x i );(2) ⽤f (x i )的算术平均值1()Nii f x N =∑近似估计上述积分类⽐：“积分”同“求和”，“d x ”同“1/N ”，“()x ?”同“服从()x ?分布的随机数”；(3) 停⽌条件：⾄⾜够⼤的N 停⽌；或者误差⼩于某值停⽌。

3. 利⽤随机数模拟各种分布的随机现象，进⽽解决实际问题。

⼆、优缺点优点：能够⽐较逼真地描述具有随机性质的事物的特点及物理实验过程；受⼏何条件限制⼩；收敛速度与问题的维数⽆关；误差容易确定。

缺点：收敛速度慢；误差具有概率性；进⾏模拟的前提是各输⼊变量是相互独⽴的。

三、应⽤随机模拟实验，随机最优化问题，含有⼤量不确定因素的复杂决策系统进⾏风险模拟分析（⾦融产品定价、期权）。

（⼆）⽤蒙特卡罗法求事件概率⼀、著名的“三门问题”源⾃博弈论的⼀个数学游戏：参赛者⾯前有三扇关闭的门，其中⼀扇门的后⾯藏有⼀辆汽车，⽽两扇门的后⾯各藏有⼀只⼭⽺。

参赛者从三扇门中随机选取⼀扇，若选中后⾯有车的那扇门就可以赢得该汽车。

当参赛者选定了⼀扇门，但尚未开启它的时候，节⽬主持⼈会从剩下的两扇门中打开⼀扇藏有⼭⽺的门，然后问参赛者要不要更换⾃⼰的选择，选取另⼀扇仍然关着的门。

数学建模 π的计算

注：泰勒级数法中的
可将算到 200 位。
二、泰勒级数法考虑反正切函数的泰勒级数
取 x =1，n =1e3(3,4,5,6,7)
程序： clear tic k=1e3(3,4,5,6,7); x=1; n=1:k; pai=sum((-1).^(n-1).*x.^(2.*n-1)./(2.*n-1)) toc 结果显示：花费时间长！准确度差！！原因：x =1 时得到的 arctan1 的展开式收敛太慢！ Maqin 公式法提高收敛速度：
1、梯形公式法将扇形 G 分为 n 个宽度相等的部分 Gk（1 ≤ k ≤ n） Gk：曲边梯形左右边界平行，上方边界为曲线 n→∞ Gk→梯
G G
形 Gk 梯形面积：
等价于
π
取 n =10k （k=3, 4, …,7），按上述方法，通过计算扇形面积计算 π 的近似值。程序
Clear tic n=1e3(3, 4, …,7); x=linspace(0,1,n+1); y=(1-x.^2).^0.5; pai=(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)))*4/n toc
rand(n) 产生 n×n 的随机方阵 rand(m,n) 产生 m×n 的随机矩阵需要产生两个[0, 1]上均匀分布的随机数 x，y，表示一个点的坐标，这个点落在单位正方形内每个位置的机会均等，当 x2 + y2 ≤ 1 时，该点落在扇形内。程序： clear tic n=1e7; x=rand(1,n); y=rand(1,n); c=x.^2+y.^2; m=sum(c<=1); k=m/n; pai=4*k toc 四、蒲丰（Buffon）掷针试验 1777 年法国数学家蒲丰提出了一种用于计算 π 的随机掷针试验，步骤如下：（1）取一张白纸，在上面画许多间距为 d 的等距平行线；（2）取一根长度为 l ( l < d ) 的均匀直针，随机地向画有平行线的纸

实验9 随机模拟

（5）二项分布随机数 1) binornd(n, p)：产生一个二项分布随机数 2) binornd(n,p,m,n)产生m行n列的二项分布随机数例4、产生B(10, 0.8)上的一个随机数，15个随机数， 3行6列的随机数。命令 (1) y1=binornd(10,0.8) (2) y2=binornd(10,0.8,1,15) (3) y3=binornd(10,0.8,3,6)
生成Y = min{ X 1 , X 2 , , X n }的随机数
（2）离散分布的直接抽样法设分布律为P(X=xi)= pi , i=1, 2, ... ① 产生均匀随机数u，即u~U(0,1) ②
2,3, xi 若p1 + ... + pi −1 < u ≤ p1 + ... + pi , (i = X = 若u ≤ p1 x1
练习3: 掷一枚骰子两次,比较掷出的点数之和为9 和为10这两个事件何者更容易发生.,
二中心极限定理中心极限定理讨论的是充分大的n，互相独立的随机变量X1,X2,……..,Xn的和的分布问题。即：
∑X
i =1
n
i
近似服从正态分布。
{ X i , i = 1, 2,3}
例7、设{Xi ,i=1,2,3…} 是一些独立同分布的随机变量且它们都服从泊松分布P(λ)，则部分和
（1）均匀分布U(a,b) 1)unifrnd (a,b)产生一个[a,b] 均匀分布的随机数 2)unifrnd (a,b,m, n)产生m行n列的均匀分布随机数矩阵当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它。

21——matlab——用蒙特卡罗方法计算pi的值

用蒙特卡罗( Monte Carlo ) 投点法计算的值程序：n=10000;m=0; a=2;for i=1:nx=rand(1)*a/2;y=rand(1)*a/2;if x^2+y^2<=(a/2)^2m=m+1;endendfprintf('蒙特卡罗投点法计算的pi为：%f\n',4*m/n)结果：蒙特卡罗投点法计算的pi为：3.117200蒙特卡罗投点法计算的pi为：3.132400蒙特卡罗投点法计算的pi为：3.165600蒙特卡罗投点法计算的pi为：3.135200蒙特卡罗投点法计算的pi为：3.146800蒙特卡罗投点法计算的pi为：3.140400蒙特卡罗投点法计算的pi为：3.114800蒙特卡罗投点法计算的pi为：3.117200 蒙特卡罗投点法计算的pi为：3.154800 蒙特卡罗投点法计算的pi为：3.140400 蒙特卡罗投点法计算的pi为：3.121200 蒙特卡罗投点法计算的pi为：3.134400 蒙特卡罗投点法计算的pi为：3.122800 蒙特卡罗投点法计算的pi为：3.127600 蒙特卡罗投点法计算的pi为：3.147200 蒙特卡罗投点法计算的pi为：3.145200 蒙特卡罗投点法计算的pi为：3.158400 蒙特卡罗投点法计算的pi为：3.147600 蒙特卡罗投点法计算的pi为：3.147600 蒙特卡罗投点法计算的pi为：3.146400 蒙特卡罗投点法计算的pi为：3.112000 蒙特卡罗投点法计算的pi为：3.180000 蒙特卡罗投点法计算的pi为：3.120000 蒙特卡罗投点法计算的pi为：3.153600 蒙特卡罗投点法计算的pi为：3.144000 蒙特卡罗投点法计算的pi为：3.148000 >>。

圆周率的计算-MATLAB

n0
1 16n
4 8n
1
2 8n
4
1 8n
5
1 8n
6
.
该公式的最大优点在于：经后来人将该公式变形后打破了传统的计算方法，可以直接计算圆周率的任意第n位数，而不是先计算前面的n-1位数．
1997年，Fabrice Bellard发表了一个比BBP算法更快的公式
1
1 n
6边形
12边形
24边形
圆... 3.142857..., (n 96)
71
7
刘徽魏晋数学家
所谓“割圆术”，是用圆内接正多边形的面积去无限逼近圆面积以求取圆周率的方法。在刘徽看来， “割之弥细，所失弥少，割之又割，以至于不可割，则与圆台体而无所失矣。”
公式推导：
y
b f ( x)dx
a
s1 s2
sn
n
(
i 1
f ( xi1 ) 2
f ( xi )( xi
xi1 )) o
s4 s3 a x1s2 s1 x2 x3
bx
当区间划分为n等分时
——trapz(x,y)
b a
f (x)dx [
f
(a) 2
f (b)
n i 1
ba f (xi )] n , xi
慢
！能否算得更快一点，更精确一点
arctan1 arctan1
4
2
3
4arctan1 arctan 1
4
5
239
(Machin公式)
三、数值积分法
基本思想
1 0
1
1 x2
dx
4
4
11 0 1 x2

综合实验三蒲丰投针问题实验

综合实验三蒲丰投针问题实验一、实验目的1. 掌握几何概型、熟悉Monte Carlo 方法的基本思想；3．会用MATLAB 实现简单的计算机模拟二、实验内容在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述．由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难．有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用．在这种情况下，可以考虑采用Monte Carlo 方法。

下面通过例子简单介绍Monte Carlo 方法的基本思想．Monte Carlo 方法是计算机模拟的基础，它的名字来源于世界著名的赌城——摩纳哥的蒙特卡洛，其历史起源于1777年法国科学家蒲丰提出的一种计算圆周π的方法——随机投针法，即著名的蒲丰投针问题。

这一方法的步骤是：1) 取一张白纸，在上面画上许多条间距为d 的平行线，见图8.1(1)2) 取一根长度为()l l d <的针，随机地向画有平行直线的纸上掷n 次，观察针与直线相交的次数，记为m3）计算针与直线相交的概率．由分析知针与平行线相交的充要条件是 ϕs i n 21≤x 其中πϕ≤≤≤≤0,20d x 建立直角坐标系),(x ϕ,上述条件在坐标系下将是曲线所围成的曲边梯形区域，见图 8．l （2）．由几何概率知(*)22s i n 210d l d d G g p ππϕϕπ===⎰的面积的面积 4）经统计实验估计出概率,n m P ≈由(*)式即?2=⇒=ππd l n m Monte Carlo 方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量．然后通过模拟一统计试验，即多次随机抽样试验（确定m 和n ），统计出某事件发生的百分比．只要试验次数很大，该百分比便近似于事件发生的概率．这实际上就是概率的统计定义．利用建立的概率模型，求出要估计的参数．蒙特卡洛方法属于试验数学的一个分支．问题：(1) 经过n次试验后圆周率估计与的圆周之间的差的绝对值的规律是？其中n分别取100，1000，2000，5000，10000，20000，50000(2) 参数l，d的不同选择，会对圆周率的估计有什么影响？可以选择d为l.5倍，2倍，3倍，4倍，5倍，8倍，10倍，20倍，50倍三、实验要求写出实验步骤、结果显示及分析四、实验分析以x 表示针的中点与最近一条平行线的距离，以j表示针与此线间的交角．显然0≤x≤a/20≤j≤p针与平行线相交的充要条件是x≤lsin(j)/2因(x,j)在图(4)中下面的矩形中等可能地取点，可见针与平行线相交的概率p 为图(4)正弦曲线线段与横轴围成的面积同图(4)中矩形面积的比．经计算得p= 另一方面得到如大量得投针实验，利用大数定理知：随着实验次数的增加，针与平行线相交的频率依概率收敛到概率p．那么在上式中以频率代替相应的概率p，则可以获得圆周率p的近似值．下面的程序是用matlab语言编写的计算机模拟投针以计算p 的近似值的程序．五、实验步骤1.编写MATLAB程序cleard=2l=0.5counter=0n=100x=unifrnd(0,d/2,1,n)fi=unifrnd(0,pi,1,n)for i=1:nif x(i)<1*sin(fi(i))/2counter=counter+1endendfren=counter/npihat=2*1/(d*fren)sqrt((pihat-pi)^2)结果显示:fren = 0.3300pihat =3.0303ans =0.1113以此类推：将n=1000，2000，5000，10000，20000，50000分别代入，可得：当n=1000时，fren =0.3240pihat =3.0864ans =0.0552当n=2000时，fren =0.3230pihat =3.0960ans =0.0456当n=5000时，fren =0.3204pihat =3.1211ans =0.0205当n=10000时，fren =0.3190pihat =3.1348ans =0.0068当n=20000时，fren =0.3172pihat =3.1521ans =0.0105当n=50000时，fren =0.3177pihat =3.1478ans =0.00622.改变d的取值，分别为1.5，2 ，3 ，4，5，8，10，20，50倍仍用1中的程序：cleard=3l=0.5counter=0n=100x=unifrnd(0,d/2,1,n)fi=unifrnd(0,pi,1,n)for i=1:nif x(i)<1*sin(fi(i))/2counter=counter+1endendfren=counter/npihat=2*1/(d*fren)sqrt((pihat-pi)^2)结果显示：d为1.5倍时fren =0.2300pihat =2.8986ans =0.2430d为2倍时fren =0.1700pihat =2.9412ans =0.2004d为3倍时fren =0.1100pihat =3.0303ans =0.1113d为4倍时fren =0.0800pihat =3.1250ans =0.0166d为5倍时fren =0.0600pihat =3.3333ans =0.1872d为8倍时fren =0.0400pihat =3.1250ans =0.0211d为10倍时fren =0.0300pihat =3.3333ans =0.1872d为20倍时fren =0.0100pihat =5ans =1.8539d为50倍时fren =0pihat =Infans =Inf六、结果分析1．经过n次试验后圆周率估计与的圆周π之间的差的绝对值的规律是：n的次数取值越多，圆周率估计与的圆周π之间的差的绝对值越小：圆周率越接近真值。

蒲丰氏问题

x ≤ l ⋅ sin θ
针在平行线间的位置
蒲丰氏问题的算法
如何产生任意的（x,θ）？x 在［0，a］上任意取值，表示，］上任意取值，表示x 在［0，a］上是均匀分布的，］上是均匀分布的，其分布密度函数为：类似地，θ的分布的分布密度函数的分布为：因此，产生任意的（x,θ）产生任意的（）产生任意的的过程就变成了由f1(x)抽样及的过程就变成了由抽样x及抽样抽样θ的过程了由f2(θ)抽样的过程了抽样的过程了。由此得到：
当 −1 ≤ x ≤ 1 其它
令φ=2θ，则θ在［0,π］上均匀分布在］上均匀分布，作变换 x = ρ cos θ y = ρ sin θ 其中0≤ρ≤1，0≤ρ≤π，则 x cos θ = x2 + y2 y sin θ = x2 + y2 (x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单位圆内均匀分布，则θ在［0,π］上均匀分布，由于 x2 − y2 2 2 cos ϕ = cos 2θ = cos θ − sin θ = 2 x + y2 2 xy sin ϕ = sin 2θ = 2 sin θ cos θ = 2 x + y2
2l 2l ≈ π= aP as N
请输入平行线相间的半距离（默认值请输入平行线相间的半距离（默认值a=4.0） a=4 ）请输入针的半长度（默认值请输入针的半长度（默认值l=3.0） l=3 ）请输入模拟试验次数（默认值请输入模拟试验次数（默认值N=100000） N=100000 ）
s=0; for n=1:1:N; x2=mod(aa*x1,MM); x1=x2; randomx=x2*1.0/MM; x=a*randomx; x2=mod(aa*x1,MM); x1=x2; randomx=x2*1.0/MM; phi=pi*randomx; if (x<=l*sin(phi)) s=s+1; end end

蒙特卡罗方法详解与MATLAB实现

通过某种试验, 得到Ｎ个观察值r1, r2, …, rN（用概率语言来说, 从分布密度函数f(r)中抽取Ｎ个子样r1, r2, …, rN, ）, 将相应的Ｎ个随机变量的值g(r1), g(r2), …, g(rN) 的算术平均值
1 N
g N N i1 g(ri )
作为积分的估计值（近似值）。
0.2 命中8环
用计算机作随机试验（射击） 0.5 命中9环
的方法为, 选取一个随机数ξ, 按右边所列方法判断得到成绩。
命中10环
这样, 就进行了一次随机试验（射击）, 得到了一次成绩ｇ(r), 作Ｎ次试验后, 得到该运动员射击成绩的近似值
g N
1 N
N
g(ri )
i 1
2. 蒙特卡罗方法的收敛性, 误差
显然, 当给定置信度α后, 误差ε由σ和N决定。要减小ε, 或者是增大N, 或者是减小方差σ2。在σ固定的情况下, 要把精度提高一个数量级, 试验次数N需增加两个数量级。因此, 单纯增大N不是一个有效的办法。
另一方面, 如能减小估计的均方差σ, 比如降低一半, 那误差就减小一半, 这相当于N增大四倍的效果。因此降低方差的各种技巧, 引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。
4) 具有同时计算多个方案与多个未知量的能力
对于那些需要计算多个方案的问题, 使用蒙特卡罗方法有时不需要像常规方法那样逐个计算, 而可以同时计算所有的方案, 其全部计算量几乎与计算一个方案的计算量相当。例如, 对于屏蔽层为均匀介质的平板几何, 要计算若干种厚度的穿透概率时, 只需计算最厚的一种情况, 其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。
➢ 计算机模拟试验过程