用蒙特卡洛方法求解圆周率

-1-用蒙特卡洛方法求解圆周率兰州交大数理学院应用物理07祁正荣题目：圆周率的求解关键词：蒙特卡洛方法、圆周率内容：一、蒙特卡洛基本思想：计算机模拟经常采用随机模拟方法或统计试验方法，这就是蒙特卡洛方法。它是通过不断产生随机数序列来模拟过程。自然界中有的过程本身就是随机的过程，物理现象中如粒子的衰变过程、粒子在介质中的输运过程…等。当然蒙特卡洛方法也可以借助概率模型来解决不直接具有随机性的确定性问题。对求解问题本身就具有概率和统计性的情况，例如中子在介质中的传播，核衰变的过程等，我们可以直接使用直接蒙特卡洛模拟方法。该方法是按照实际问题所遵循的概率统计规律。直接蒙特卡洛方法最充分体现出蒙特卡洛方法无可比拟的特殊性和优越性，因而在物理学的各种各样的问题中得到广泛的应用。该方法也就是通常所说的“计算机实验”。蒙特卡洛方法也可以人为地构造出一个合适的概率模型，依照对该模型进行大量的统计实验，使它的某些统计参量正好是待求问题的解。这也就是所谓的间接蒙特卡洛方法。我在这里求解圆周率就是用的间接蒙特卡洛方法。二、圆周率求解的基本思想：根据圆面积的求解公式s=πr2可推导出，圆周率的求解公式π=s/r2，根据此公式，可知，只要计算出圆的面积，测得圆的半径即可求得圆周率。圆的半径可直接测得，但是圆的面积并不好求，传统的做法有微圆分割法，如右图，将圆周无限分割成三角形这些三角形都是等边的，可以通过求解三角形的面积间接求得圆的面积，但是这样的圆周率并不精确。在本题中采用计算机模拟来间接求得圆的面积，具体做法如下：如图，在水平光滑的桌面上画边长为2的-2-正方形，一个单位圆(即半径为1的圆)与之相切，直角坐标系原点取单位圆的圆心，且二坐标轴与正方形平行。准备一张纸，上面分别写上两个标题：N圆和N正。在桌面的上方随便地投下缝衣针，每投下一次，就观察针尖落点的位置，如果针尖落到正方形内，在N正下加一条线；如果针尖落到圆内，在N圆下加一条线，落人正方形外时不记录。在多次投针后，就可以分别统计出针尖落人正方形或圆内的次数N正和N圆，用N圆除以N正然后乘以4就能得出圆周率。当重复投针实验到达几十万时，圆周率的值接近3．14或3．15。这其实用了概率的方法，计算圆面积。用圆面积除以半径的方，就是圆周率。打点次数越多，得到的值就越精确。这种方法也称为巴福昂投针实验。三、用蒙特卡洛方法求解圆周率：上面已经知道了蒙特卡洛应用的今本思想以及圆周率求解的基本原理与方法，那么我们用蒙特卡洛方法求解圆周率，从0到2之间任意取一组数，这一组数的每个元素包含两个变量x和y，x代表所取点的横坐标，y代表所取点的纵坐标，x和y由计算机随机抽取，然后编辑公式筛选符合条件的点，筛选条件为，x和y都被包含在圆内。然后设定一个计数器m，每当有一个点符合条件，则m加1，如此一直循环到取够所设定的点数，然后用m值除以设定的点数，得到一个值，然后再乘以4，这时得到的值就是圆的面积，然后用得到的这个面积除以半径的平方（因为设定的圆的半径为1，所以具体操作时不再考虑，得到的圆的面积直接就是圆周率的值）说明：本例中先取0到20000之间的随机数，然后再用这些点除以10000即可得到0到2之间的随机数，随机数精度为五位有效数。以下是用c++语言模拟的求解圆周率的程序：#includeiostream#includeiomanip//setw函数,setprecision函数#includetime.h//time函数#includestdlib.h//srand,rand函数usingnamespacestd;intmain(){intdwCount;-3-cout请输入打点的数目endl;cindwCount;//自己输入所要取的点数。doublej,k;doublex,y;doublem=0;for(inti=0;idwCount;i++){j=(int)(20000.0*rand()/(RAND_MAX+0));//取从0到20000之间的随机数k=(int)(20000.0*rand()/(RAND_MAX+0));//取从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))=1)//判断所取得点是否在圆内{m++;//如果这个所取得点在圆内，则计数器加1}coutxyendl;/*输出所取得随机数，没有实际意义，只为了直观的看见所取得随机数*/}coutmendl;doubles=m/dwCount;//用圆内的点数除以所取得点数,得到一个比例coutsendl;doublet=4*s;/*用刚才的点数乘以4就是圆的面积，因为这里做了特殊，所以这个值也就是圆周率*/couttendl;//输出圆周率return0;}在vc6.0中调试好环境，运行程序，得到如下结果-4-（举例，输入点数为1000000）：结论：用蒙特卡洛方法求解圆周率影响精度的几个因素a.所取得点数越多，得到的精度越高。b.所取得随机数的精度越高，得到的圆周率的精度也就越高。四、用平均打点的方法求解圆周率：还是用了上面求解圆周率的基本思想，方法也基本相同，唯一有区别的是，上面用蒙特卡洛方法模拟的时候是用的随机取点，而在这里，取点的方法是等概率的取点，在正方形内，每个点只出现一次，且每个点出现的概率都相等，具体实现方法，让x从0到2之间递增，每次增0.0001（这个值由自己设定，精度越高最后求解的圆周率的精度也就越高），然后对应每个x的值，让y也从0到2之间递增，每次增0.0001，如此，所取得点数总数为20000*20000=400000000个点。然后用判断语句进行判断所取得点是否在圆内，如果在园内，计数器m加1，同上，最后用m除以总的点数，然后再乘以4，所得结果就是圆的面积，然后用圆的面积除以半径平方，因为在这里面我我们做了特殊处理，圆的半径为1，所以刚才求得的圆的面积就是圆周率。以下是用c++语言模拟的求解圆周率的程序：#includeiostream#includecmath#includewindows.husingnamespacestd;voidmain(){doublem=0;//计算器，保存的是落在园内的点数。设定初始值为0。doublex=0;//x的初始值为0-5-doubley=0;//y的初始值为0for(doublen=0;n20000;n++)/*这是一个循环，循环的是横坐标x，从0到2之间每次加0.0001，总共循环20000次。*/{y=0;for(doubled=0;d20000;d++)/*x坐标一定时横坐标y从0到2之间每次加0.0001，做20000次循环。*/{if(((x-1)*(x-1)+(y-1)*(y-1))=1)/*判断点（x，y）是否在圆内，如果是则计数器加1*/{m++;}y=y+0.0001;}x=x+0.0001;//y循环完以后，x加0.0001，继续下一次循环。}doubles=m/400000000;doublet=4*s;//t为圆的面积，这里做了特殊，也就是所要求的圆周率couttendl;}在vc6.0中调试好环境，运行程序，得到如下结果，参数为上面程序中设定好的结论：用平均打点的方法也可以较精确的求得圆周率，但是效率比蒙特卡洛方法要低。上面两个例子，充分的体现了这一结论。在这种方法中，影响圆周率精度的重要因素是取点的精度，也就是取点的数目。