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

多种方法计算圆周率并比较精确度【摘要】本文介绍了多种方法求圆周率的近似值并对各种方法进行精确度的比较得出具体情况选择的方法，且通过mathematica编程模拟实验过程，得出各种方法的特点。【关键字】圆周率数值积分法泰勒级数法蒙特卡罗法拉马努金公式法0.引言平面上圆的周长与直径之比是一个常数，称为圆周率，记作。在很长的一段时期，计算的值是数学上的一件重要的事情。有数学家甚至说：“历史上一个国家所得的圆周率的准确程度，可以作为衡量一个国家当时数学发展的一面旗帜。”足以见圆周率扮演的是角色是如此举足轻重。作为经常使用的数学常数，它的计算已经持续了2500多年了，到今天都依然在进行着，中间涌现出许多的计算方法，它们都各有千秋，在此，我们选择几种较典型的方法，包括数值积分法，泰勒级数法，蒙特卡罗法，韦达公式法，拉马努金公式法以及迭代法来和大家一起体验的计算历程，同时利用mathematica通过对各种方法作精确度的比较得出选择的优先顺序，为相关的理论研究提供一定的参考价值。1.数值积分法以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G是一个扇形，由曲线y=211x(x[0,1])及两条坐标轴围成，它的面积S=4。算出了S的近似值，它的4倍就是的近似值。（411102dxx）用一组平行于y轴的直线x=ix(1in-1，a=0x1x2x...nx=b)将曲边梯形T分成n个小曲边梯形，总面积S分成这些小曲边梯形的面积之和。如果取n很大，使每个小曲边梯形的宽度很小，可以将它上方的边界f(x)(1ixxix)近似地看作直线，将每个小曲边梯形近似当作梯形来求面积，就得到梯形公式，如果更准确一些将每个小曲边梯形的上边界近似地看作抛物线，就得到辛普森公式。梯形公式：S]2...[10121yyyyynabn辛普森公式：S)]...(4)...(2)[(62123211210nnnyyyyyyyynabMathematica程序如下：n=1000;y[x_]:=4/(1+x^2);s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/(6*n);Print[{N[s1,20],N[s2,20],N[Pi,30]}]注：以上s1,s2分别是用梯形共识和辛普森公式计算出的。最后一句中的N[s1,20]表示s1的前20为准确有效数字组成的近似值。N[Pi,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.泰勒级数法利用反正切函数的泰勒级数：Arctanx=x-...12)1(...5312153kxxxkkX=1代入上式得到莱布尼兹级数：1112)1(4kkk令a=arctan21,b=4-a,则tanb=tan(4-a)=aatan4tan1tan4tan=2111211=31因此b=arctan31,即4-arctan21=atctan31，从而得到4=arctan21+arctan31即=4(arctan21+arctan31）令a=arctan51，tana=51，则容易算出tan2a=125,tan4a=119120,Tan(4a-4)=aa4tan4tan14tan4tan=11912011119120=23914a-4=arctan23914=4a-arctan2391=4arctan51-arctan2391即=16arctan51-4arctan2391利用mathematica编写程序如下：T[x_,n_]:=Sum[(-1)^(k-1)x^(2k-1)/(2k-1),{k,1,n}]4N[T[1,1000],20]N[4(T[1/2,1000]+T[1/3,1000]),50]N[16T[1/5,1000]-4T[1/239,1000],50]N[Pi,50]注：4N[T[1,1000],20]为用的方法求解近似值，N[4(T[1/2,1000]+T[1/3,1000]),50]表示用的方法求解近似值，N[16T[1/5,1000]-4T[1/239,1000],50]表示用的方法求解近似值。输出结果为：3.14059265383979292603.14159265358979323853.14159265358979323853.1415926535897932385观察可见，的方法精确度不高，当n=5000,10000,时用方法所得结果为{3.1413926535917932384,3.1414926535900432385}，即当n的值增大较大幅度时近似值的精确度提高不大，有效数依然只有4位。而由上面表格中的结果可知当n取1000时，两种方法所得的前20位有效数字完全与的前20位数相同，对两种方法取150个有效数字时所得结果为：3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408133.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408133.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813观察可知，用两种方法所求近似值前150位有效数字均和依然相同，可见用这两种方法求解的近似值的精确度非常高！当计算要求精确度比较高时，可以选取用两种方法均可。3.蒙特卡罗法单位圆的41是一个扇形G,它是边长为1的单位正方形G1的一部分，如图，0.00.20.40.60.81.00.00.20.40.60.81.0单位正方形G1的面积S1=1。只要能够求出扇形G的面积S在正方形G1的面积S1中所占的比例k，即可以求得S，乘与4即为的值。此处利用正方形中随机地投入很多点，使所投的每个点落在正方形中每个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投点的总数n的比nm可以作为k的近似值。在mathematica中用Random[]函数产生[0,1]之间的随机数。Mathematica编写程序如下：Module[{n=1000,p={},m,x,y,k},Do[m=0;Do[x=Random[];y=Random[];If[x^2+y^21,m++],{k,1,n}];AppendTo[p,N[4m/n]],{t,1,10}];Sum[p[[t]],{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条圆周率的计算公式。nnnnn404396)26391103()!()!4(980221这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。Mathematica编写程序如下：a[n_]:=Sum[(4k)!*(1103+26390k)/((k!)^4*369^(4k)),{k,0,n}]N[1/(2Sqrt[2]/9801*a[1000]),15]N[Pi,15]输出结果为：3.14159262864498，3.14159265358979可见用该公式可以精确到8位有效数字。1989年，大卫·丘德诺夫斯基和格雷高里·丘德诺夫斯基兄弟将拉马努金公式改良，这个公式被称为丘德诺夫斯基公式，每计算一项可以得到15位的十进制精度。1994年丘德诺夫斯基兄弟利用这个公式计算到了4,044,000,000位。以下是丘德诺夫斯基公式的一个便于计算机计算的版本。Mathematica编写程序如下：b[n_]:=Sum[((6k)!(13591409+545140134k))/((3k)!(k!)^3(-640320)^(3k)),{k,0,n}]N[426880*Sqrt[10005]/b[1000],150]N[Pi,150]输出结果：近似值3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813准确值3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813由表格可知，当n取1000时，该公式计算所得的近似值前150位与圆周率的精确值前150位完全相同，可见该公式计算所得的精确度已经非常高，同时这个公式也非常适合计算机编程，是目前计算机使用较快的一个公式。5.算法比较由以上