《用蒙特卡洛方法求解圆周率》由会员分享,可在线阅读,更多相关《用蒙特卡洛方法求解圆周率(5页珍藏版)》请在金锄头文库上搜索。
1、用蒙特卡洛方法求解圆周率兰州交大数理学院 应用物理 07 祁正荣题目:圆周率的求解关键词:蒙特卡洛方法、圆周率内容:一、蒙特卡洛基本思想:计算机模拟经常采用随机模拟方法或统计试验方法,这就是蒙特卡洛方法。它是通过 不断产生随机数序列来模拟过程。自然界中有的过程本身就是随机的过程,物理现象中如粒子的衰变过程、粒子在介质中的输运过程等。当然蒙特卡洛方法也可以借助概率模型来解 决不直接具有随机性的确定性问题。对求解问题本身就具有概率和统计性 的情况,例如中子在介质中的传播,核衰变 的过程等,我们可以直接使用直接蒙特卡洛 模拟方法。该方法是按照实际问题所遵循的概率统计规律。直接蒙特卡洛方法最充分体现出
2、 蒙特卡洛方法无可比拟的特殊性和优越性,因而在物理学的各种各样的问题中得到广泛的应用。该方法也就是通常所说的“计算机实验”。蒙特卡洛方法也可以人为地构造出一个合适的概率模型,依照对该模型进行大量的统计 实验,使它的某些统计参量正好是待求问题的解。这也就是所谓的间接蒙特卡洛方法。我在 这里求解圆周率就是用的间接蒙特卡洛方法。二、圆周率求解的基本思想:根据圆面积的求解公式S=nr2可推导出,圆周率的求解公式n=s/r2,根据此公式,可 知,只要计算出圆的面积,测得圆的半径即可求得圆周率。圆的半径可直接测得,但是圆的如图,在水平光滑的桌面上画边长为2 的正方形,一个单位圆(即半径为1 的圆)与之相切
3、,直角坐标系原点取单位圆的圆心,且二 坐标轴与正方形平行。准备一张纸,上面分别写上两个标题:N圆和N正。在桌面的上方 随便地投下缝衣针,每投下一次,就观察针尖落点的位置,如果针尖落到正方形内,在 N 正下加一条线; 如果针尖落到圆内,在 N 圆下加一条线, 落人正方形外时不记录。在多 次投针后,就可以分别统计出针尖落人正方形或圆内的次数N正和N圆,用N圆除以N正 然后乘以 4 就能得出圆周率。当重复投针实验到达几十万时,圆周率的值接近314 或 315。这其实用了概率的方法,计算圆面积。用圆面积除以半径的方,就是圆周率。打点次数 越多,得到的值就越精确。这种方法也称为巴福昂投针实验。三、用蒙特
4、卡洛方法求解圆周率:上面已经知道了蒙特卡洛应用的今本思想以及圆周率求解的基本原理与方法,那么我 们用蒙特卡洛方法求解圆周率,从0 到 2 之间任意取一组数,这一组数的每个元素包含两 个变量x和y, x代表所取点的横坐标,y代表所取点的纵坐标,x和y由计算机随机抽取, 然后编辑公式筛选符合条件的点,筛选条件为,x和y都被包含在圆内。然后设定一个计数 器m,每当有一个点符合条件,则m加1如此一直循环到取够所设定的点数,然后用m 值除以设定的点数,得到一个值,然后再乘以4,这时得到的值就是圆的面积,然后用得到 的这个面积除以半径的平方(因为设定的圆的半径为1,所以具体操作时不再考虑,得到的 圆的面积
5、直接就是圆周率的值) 说明:本例中先取0到20000之间的随机数,然后再用这些点除以 10000即可得到0到2之间的随 机数,随机数精度为五位有效数。以下是用C+语言模拟的求解圆周率的程序:#include #include viomanip /setw 函数setprecision 函数#inelude vtime.h/time 函数#include vstdlib.h /srand,rand 函数using namespace std;int main()int dwCount;coutvv请输入打点的数目vve ndl;cindwCount;/自己输入所要取的点数。double j,k;
6、double x,y;double m=0;for(int i=0;ivdwCount;i+)j = (int)(20000.0*rand()/(RAND_MAX+O);取从0 到 20000 之间的随机数 k=(int)(20000.0*rand()/(RAND_MAX+O);取从0 到20000 之间的随机数 x=j/10000; /*用刚才所取得随机数除以10000就是0到2之间的随机数,5位有 效数字,代表横坐标*/y=k/10000;/*用刚才所取得随机数除以10000就是0到2之间的随机数,5位有效数字,代表纵坐标*/if(x-1)*(x-1)+(y-1)*(y-1)v=1)/判断
7、所取得点是否在圆内m+;/如果这个所取得点在圆内,则计数器加 1coutvvxvv vvyvvendl; /*输出所取得随机数,没有实际意义,只为了直 观的看见所取得随机数*/coutvvmvvendl;double s=m/dwCount;/用圆内的点数除以所取得点数,得到一个比例coutvvsvvendl;double t=4*s;/*用刚才的点数乘以4 就是圆的面积,因为这里做了特殊,所以这个值也就是圆周率*/couttendl;/输出圆周率return 0;在VC6.0中调试好环境,运行程序,得到如下结果举例,输入点数为1000000):结论:用蒙特卡洛方法求解圆周率影响精度的几个因素
8、a. 所取得点数越多,得到的精度越高。b. 所取得随机数的精度越高,得到的圆周率的精度也就越高。四、用平均打点的方法求解圆周率:还是用了上面求解圆周率的基本思想,方法也基本相同,唯一有区别的是,上面用蒙 特卡洛方法模拟的时候是用的随机取点,而在这里,取点的方法是等概率的取点,在正方形 内,每个点只出现一次,且每个点出现的概率都相等,具体实现方法,让x从0到2之间 递增,每次增0.0001(这个值由自己设定,精度越高最后求解的圆周率的精度也就越高), 然后对应每个x的值,让y也从0到2之间递增,每次增0.0001,如此,所取得点数总 数为20000*20000=400000000个点。然后用判断
9、语句进行判断所取得点是否在圆内, 如果在园内,计数器m加1,同上,最后用m除以总的点数,然后再乘以4,所得结果就 是圆的面积,然后用圆的面积除以半径平方,因为在这里面我我们做了特殊处理,圆的半径 为 1,所以刚才求得的圆的面积就是圆周率。以下是用C+语言模拟的求解圆周率的程序:#include #inClude #inClude using namespaCe std; void main()double m=0;double x=0;/计算器,保存的是落在园内的点数。设定初始值为0。/x 的初始值为0double y=0;/y 的初始值为0for (double n=0; nv20000;
10、n+) /*这是一个循环,循环的是横坐椒,从0到2之间每次加0.0001,总共循环20000 次。*/y=0;for (double d=0;dv20000;d+)/*x坐标一定时横坐标y从0到2之间每次加0.0001,做20000 次循环。*/if(x-1)*(x-1)+(y-1)*(y-1)v = 1) /*判断点X, y)是否在圆内,如果是则计数器加1*/m+;y=y+0.0001;x=x+0.0001;/y循环完以后,x加0.0001,继续下一次循环。double s=m/400000000;double t=4*s;/t 为圆的面积,这里做了特殊,也就是所要求的圆周率coutvvtvvendl;在 vc6.0 中调试好环境,运行程序,得到如下结果,参数为上面程序中设定好的S E:circ lek ire leDebugc ire k.exe-h|x-14159Presskw寸 to continue1结论:用平均打点的方法也可以较精确的求得圆周率,但是效率比蒙特卡洛方法要低。上面两 个例子,充分的体现了这一结论。在这种方法中,影响圆周率精度的重要因素是取点的精度 也就是取点的数目。