多种方法计算Pi并且精确度比较.doc

多种方法计算圆周率并比较精确度 【摘要】本文介绍了多种方法求圆周率的近似值并对各种方法进行精确度的比较得出具体情况选择的方法，且通过mathematica编程模拟实验过程，得出各种方法的特点。【关键字】圆周率 数值积分法 泰勒级数法 蒙特卡罗法 拉马努金公式法0.引言 平面上圆的周长与直径之比是一个常数，称为圆周率，记作。在很长的一段时期，计算的值是数学上的一件重要的事情。有数学家甚至说：“历史上一个国家所得的圆周率的准确程度，可以作为衡量一个国家当时数学发展的一面旗帜。”足以见圆周率扮演的是角色是如此举足轻重。 作为经常使用的数学常数，它的计算已经持续了2500多年了，到今天都依然在进行着，中间涌现出许多的计算方法，它们都各有千秋，在此，我们选择几种较典型的方法，包括数值积分法，泰勒级数法，蒙特卡罗法，韦达公式法，拉马努金公式法以及迭代法来和大家一起体验的计算历程，同时利用mathematica通过对各种方法作精确度的比较得出选择的优先顺序，为相关的理论研究提供一定的参考价值。1. 数值积分法以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G是一个扇形，由曲线y=(x0,1)及两条坐标轴围成，它的面积S=。算出了S的近似值，它的4倍就是的近似值。（）用一组平行于y轴的直线x=(1in-1，a=.=b)将曲边梯形T分成n个小曲边梯形，总面积S分成这些小曲边梯形的面积之和。如果取n很大，使每个小曲边梯形的宽度很小，可以将它上方的边界f(x)()近似地看作直线，将每个小曲边梯形近似当作梯形来求面积，就得到梯形公式，如果更准确一些将每个小曲边梯形的上边界近似地看作抛物线，就得到辛普森公式。梯形公式：S辛普森公式：SMathematica程序如下：n=1000;yx_:=4/(1+x2);s1=(Sumyk/n,k,1,n-1+(y0+y1)/2)/n;s2=(y0+y1+2*Sumyk/n,k,1,n-1+4*Sumy(k-1/2)/n,k,1,n)/(6*n);PrintNs1,20,Ns2,20,NPi,30注：以上s1,s2分别是用梯形共识和辛普森公式计算出的。最后一句中的Ns1,20表示s1的前20为准确有效数字组成的近似值。NPi,30是的前30位有效数字组成的近似值。Print语句表示将方括号内的数显示出来。再取n=5000,10000时分别算出的近似，记录入表格如下：梯形公式辛普森公式圆周率的准确值n=10003.14159248692312657183.14159265358979323853.14159265358979323846264338328n=50003.14159264692312657183.14159265358979323853.14159265358979323846264338328n=100003.14159265192312657183.14159265358979323853.14159265358979323846264338328注：下划线标记的数字为用对应方法所求得的精确到的位数。观察发现，当n的取值越大，即当a,b之间分的越细时，计算所得的的精确度越高。n 从1000到10000，梯形公式所求得的近似值精确度增加两位，而辛普森公式所求的的近似值有效位数没有变化，可见n的变化对梯形公式影响较大。但总的而言利用辛普森公式所求得的近似值精度比用梯形公式所求得的近似值很多。当要求精确度较高时，两种方法中可优先选择辛普森公式法。2. 泰勒级数法利用反正切函数的泰勒级数：Arctan x=x-X=1代入上式得到莱布尼兹级数：令a=arctan,b=-a,则tan b=tan(-a)=因此b=arctan,即-arctan=atctan，从而得到=arctan+arctan 即=4(arctan+arctan）令a=arctan，tan a=，则容易算出tan2a=,tan4a=,Tan(4a-)=4a-=arctan=4a-arctan=4arctan-arctan 即=16arctan-4arctan利用mathematica编写程序如下：Tx_,n_:=Sum(-1)(k-1)x(2k-1)/(2 k-1),k,1,n4NT1,1000,20N4(T1/2,1000+T1/3,1000),50N16T1/5,1000-4 T1/239,1000,50NPi,50注：4NT1,1000,20为用的方法求解近似值，N4(T1/2,1000+T1/3,1000),50表示用的方法求解近似值，N16T1/5,1000-4 T1/239,1000,50表示用的方法求解近似值。输出结果为：3.14059265383979292603.14159265358979323853.14159265358979323853.1415926535897932385观察可见，的方法精确度不高，当n=5000,10000,时用方法所得结果为3.1413926535917932384,3.1414926535900432385，即当n的值增大较大幅度时近似值的精确度提高不大，有效数依然只有4位。而由上面表格中的结果可知当n取1000时，两种方法所得的前20位有效数字完全与的前20位数相同，对两种方法取150个有效数字时所得结果为：3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408133.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408133.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813观察可知，用两种方法所求近似值前150位有效数字均和依然相同，可见用这两种方法求解的近似值的精确度非常高！当计算要求精确度比较高时，可以选取用两种方法均可。3. 蒙特卡罗法单位圆的是一个扇形G,它是边长为1的单位正方形G1的一部分，如图，单位正方形G1的面积S1=1。只要能够求出扇形G的面积S在正方形G1的面积S1中所占的比例k，即可以求得S，乘与4即为的值。此处利用正方形中随机地投入很多点，使所投的每个点落在正方形中每个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投点的总数n的比可以作为k的近似值。在mathematica中用Random函数产生0,1之间的随机数。Mathematica编写程序如下：Modulen=1000,p=,m,x,y,k,Dom=0;Dox=Random;y=Random;Ifx2+y21,m+,k,1,n;AppendTop,N4m/n,t,1,10; Sumpt,t,1,10/10由于程序中采用随机数思想，故而所得结果不唯一，运行三次程序所得结果为：3.158，3.1348，3.1504。观察可知，用该方法所得的近似值精确度不是很高，只是得出了两位有效数字。但思想方法比较简单，比较直观。进一步投点5000,10000。投点5000，三次运行结果为：3.1516，3.15056，3.14168投点10000，三次运行结果为：3.15032，3.1412，3.13656观察发现，随着投点数的增加，会增加结果的精确度，但是效果非常不明显。对比数值积分法和泰勒级数法，蒙特卡罗法的精确度确实比较低，但是当所求不是一个扇形面积而是其它比较特殊的图形时（如100个已知圆的公共部分G的面积）时，利用蒙特卡罗法利用计算机实现是非常简单的事情。所以，蒙特卡罗法在很多场合，特别是在对精确度要求不太高的情况下是大有用武之地的。4.拉马努金公式法1914年，印度天才数学家拉马努金在他的论文里发表了一系列共14条圆周率的计算公式。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。Mathematica编写程序如下：an_:=Sum(4k)!*(1103+26390k)/(k!)4*369(4k),k,0,nN1/(2Sqrt2/9801*a1000),15NPi,15输出结果为：3.14159262864498 ，3.14159265358979可见用该公式可以精确到8位有效数字。1989年，大卫丘德诺夫斯基和格雷高里丘德诺夫斯基兄弟将拉马努金公式改良，这个公式被称为丘德诺夫斯基公式，每计算一项可以得到15位的十进制精度。1994年丘德诺夫斯基兄弟利用这个公式计算到了4,044,000,000位。以下是丘德诺夫斯基公式的一个便于计算机计算的版本。Mathematica编写程序如下：bn_:=Sum(6k)!(13591409+545140134k)/(3k)!(k!)3(-640320)(3k),k,0,nN426880*Sqrt10005/b1000,150NPi,150输出结果：近似值3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813准确值3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813由表格可知，当n取1000时，该公式计算所得的近似值前150位与圆周率的精确值前150位完全相同，可见该公式计算所得的精确度已经非常高，同时这个公式也非常适合计算机编程，是目前计算机使用较快的一个公式。5.算法比较由以上计算圆周率的方法和计算结果，我们可以得知，当要求计算的精确度非常高时，我们可以选择泰勒级数法中的公式=4(arctan+arctan）以及=16arctan-4arctan和丘德诺夫斯基公式;当精确度要求较高时我们可以选择数值积分法中的梯形公式，辛普森公式以及拉马努金公式;当精确度要求不高时，我们可以选择泰勒级数法中的莱布尼兹级数以及蒙特卡罗法，特别是当所求不是一个扇形面积而是其它比较特殊的图形时（如100个已知圆的公共部分G的面积）时，可利用蒙特卡罗法快速简单的求出近似值。6.结束语随着计算机技术的飞速发展和计算方法的突破与创新, 计算的世界纪录正在迅速地被刷新。1949 年计算机 ENIAC 用 70 分钟的时间把 P计算到 2037位以来, 到 1999 年东京大学的 Kana-da 和Takahashi 用两台超级计算机花费了 73 小时把 P计算到 68, 719, 470, 000 位, 这期间由于计算机运算速度的加快, 计算位数越来越多, 相对使用时间越来越短。虽然这样高的精确度已经没有太多的实际意义,但这反映了现代数学科学的日新月异,反映了人类智慧向极限的挑战。就像东京大学的 Kanada所说的，“我们进行运算的主要原因是检验我们计算机的可靠性和对于运算、程序和算法的正确性,成为世界记录的保持者只是它的一个副产品而己”。参考文献1张晓贵.圆周率计算的四个时期J.辽宁教育学院学报，2009，9(17)，32 33.2何光.用蒙特卡罗方法计算圆周率的近似值.内江师范学院学报,第23卷第4期.3陶珊珊，孟凤娟.计算的几种方法. 科技信息，2010年第35期4强春晨,刘兴祥,岳育英.圆周率计算方法发展史.延安大学学报(自然科学版), 2012年6月第31卷第2期5王广超. 圆周率的数学实验算法设计. 博士专家论坛6李尚志，陈发来，张韵华，吴耀华著.数学实验 第二版7王维平.实验数学导引:数学原来可以这样做（2004）Several ways to calculate pi and compare the accuracy LUO Yan(School of Mathematics and Applied Mathematics, Sun Yat-sen university,GuangZhou,Guangdong province)Abstract:This article describes several ways to seek an approximation of pi and tells us whe