pi的计算 实验报告

1、班级：9131138502 姓名：张海洋 学号:913113850232圆周率（选做）要求：先查资料看看古人是怎样计算的，再对的各种计算方法进行研究和讨论（收敛速度等），并给出不同算法算出的小数点后第10000位的数字是什么，你觉得该数字应该是多少？第一部分：圆周率简介 圆周率是指平面上圆的周长与直径之比 （ratio of the circumference of a circle to the diameter） 。用符号（读音：pài）表示。中国古代有圆率、圆率、周等名称。它是一个常数（约等于3.141592654）。它是一个无理数，即无限不循环小数。在日常生活中，通常都用3.

2、14代表圆周率去进行近似计算。而用十位小数3.141592654便足以应付一般计算。即使是工程师或物理学家要进行较精密的计算，充其量也只需取值至小数点后几百个位。计算圆周率的方法 “历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。” 历史上最马拉松式的计算，其一是德国的Ludolph Van Ceulen，他几乎耗尽了一生的时间，计算到圆的内接正262边形，于1609年得到了圆周率的35位精度值，以至于圆周率在德国被称为Ludolph数；其二是英国的William Shanks，他耗费了15年的光阴，在1874年算出了圆周率的小数点后707位。可惜，后人发现

3、，他从第528位开始就算错了。把圆周率的数值算得这么精确，实际意义并不大。现代科技领域使用的圆周率值，有十几位已经足够了。如果用Ludolph Van Ceulen算出的35位精度的圆周率值，来计算一个能把太阳系包起来的一个圆的周长，误差还不到质子直径的百万分之一。以前的人计算圆周率，是要探究圆周率是否循环小数。自从1761年Lambert证明了圆周率是无理数，1882年Lindemann证明了圆周率是超越数后，圆周率的神秘面纱就被揭开了。 在中国，公元263年前后，刘徽提出著名的 “割圆术”求出了比较精确的圆周率。他发现：当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而多

4、边形的面积也会越来越逼近圆面积。于是，刘徽利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加倍：正十二边形、正二十四边形，正四十八边形，一直到正三七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第四位。在刘徽研究的基础上，祖冲之进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正24576边形，而得到一个结论： 3.1415926 3.1415927 同时得到 的两个近似分数：约率为227；密率为355113。他算出的 的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。 现在的人计算圆周率, 多

5、数是为了验证计算机的计算能力，还有，就是为了兴趣。第二部分：古人计算圆周率方法 古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了1.Machin公式 这个公式由英国天

6、文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Fourier Transform

7、)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nlog(n)。2、 Ramanujan公式1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。3、AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式：初值：        重复计算：    

8、;    最后计算：    这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。4、Borwein四次迭代式：初值：        重复计算：         最后计算：    这个公式

9、由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。5、Bailey-Borwein-Plouffe算法这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到了一个比BBP快40的公式：第三部分：对于的几种计算的研究和讨论：1、数值积分法（I）利用积分公式计算n=10 ans =3.1045183262;n=2

10、0 ans =3.1284648798;n=50 ans =3.1382685111;n=100 ans =3.1404170318;n=200 ans =3.1411769448;n=500 ans =3.1414874770;n=1000 ans =3.1415554669;n=2000 ans =3.1415795059;半径为1的圆称为单位圆，它的面积等于。只要计算出单位圆的面积，就算出了。在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4，于是，我们只要算出此扇形的面积，便可以计算出。但计算量较大。2、数值积分法

11、（II）利用公式 n=10 ans = 3.14242598500110; n=20 ans = 3.14180098689309; n=50 ans =3.14162598692300; n=100 ans =3.14160098692312; n=200 ans = 3.14159473692312; n=500 ans = 3.14159298692312; n=1000 ans =3.14159273692313; n=2000 ans =3.14159267442313;设分点x1,x2,xn-1将积分区间0,1分成n等分。所有的曲边梯形的宽度都是h=1/n。记yi=f(xi).则第

12、i个曲边梯形的面积A近似地等于梯形面积，即：A=(y(i-1)+yi)h/2。将所有这些梯形面积加起来就得到：A2/n2(y1+y2+yn-1)+y0+yn3、利用复化梯形算法求Pi的近似值. =1/2 (Tn+h )4、泰勒级数法计算：利用反正切函数的泰勒级数来计算(I)x=1时 =4* n = 10 ans= 3.04183961892940; n = 20 ans = 3.09162380666784; n = 50 ans = 3.12159465259101; n = 100 ans = 3.13159290355855; n = 200 ans = 3.13659268483882

13、; n = 500 ans = 3.13959265558979; n = 1000 ans = 3.14059265383979; n = 2000 ans = 3.14109265362104;x=1时得到的的展开式收敛太慢,逼近速度太慢，运算庞大，对速度造成了很大影响。5、泰勒级数法计算： 利用反正切函数的泰勒级数 来计算(II)进一步精细化 n = 10 ans= 3.1415925796063512110; n = 20 ans = 3.1415926535897574098; n = 50 ans =3.1415926535897932385; n = 100 ans = 3.14

14、15926535897932385; n = 200 ans = 3.1415926535897932385; n = 500 ans = 3.1415926535897932385; n = 1000 ans = 3.1415926535897932385; n = 2000 ans =3.1415926535897932385;当x=1时得到的的展开式收敛太慢。要使泰勒级数收敛得快，容易想到，应当使x的绝对值小于1，最好是远比1小。例如，因为，所以我们可以计算出的值，从而得到的值。这样，就使得收敛速度加快，逼近的速度大大增加.6、利用麦琴给出，推出=4() n = 10 ans= 3.14

15、15926535897916969; n = 20 ans = 3.1415926535897932385; n = 50 ans =3.1415926535897932385; n = 100 ans = 3.1415926535897932385; n = 200 ans = 3.1415926535897932385; n = 500 ans = 3.1415926535897932385; n = 1000 ans = 3.1415926535897932385; n = 2000 ans =3.1415926535897932385;对泰勒级数，随着x的减小，级数的收敛速度明显加快，

16、这启示我们另外构造相关级数来逼近。7、蒙特卡罗法计算单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积。只要能够求出扇形的面积在正方形的面积中所占的比例，就能立即得到，从而得到的值。下面的问题归结为如何求的值，这就用到了一种利用随机数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落在扇形内的点的个数与所投店的总数的比可以近似的作为的近似值。 n = 100 ans= 3.200000000000000; n = 200 ans = 3.040000000000000; n

17、 = 500 ans =3.152000000000000000 n = 1000 ans = 3.0800000000000000; n = 2000 ans = 3.1560000000000000; n = 5000 ans =3.14400000000000000; n = 10000 ans = 3.12040000000000000; n = 20000 ans =3.1620000000000000;这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为3.108，虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情况下很有用的。除以上几种方

18、法之，还有下列一些的求其他的方法：*利用高斯公式 =48+3220*、=2+1/3*(2+2/5*(2+3/7*(2+ （2+k/(2k+1)*(2+.)）).).(当 k=2799时可精确到800位)*、/6=1/2+1/2*1/(3*23)+(1*3)/(2*4)*(1/(5*25)+*、e(*i)+1=0 (欧拉公式，也称世界上最杰出的公式)*、1+(1/2)2+(1/3)2+(1/4)2+.(1/n)2=2/6*、1+(1/2)4+(1/3)4+(1/4)4+.(1/n)4=4/90*、1+(1/2)6+(1/3)6+(1/4)6+.(1/n)6=6/945*、1+(1/2)8+(1/3)8+(1/4)8+.(1/n)8=8/9450*、1+(1/2)10+(1/3)10+(1/4)10+.(1/n)10=10/93555第四部分：求的小数点后第10000位数字的几种方法：1、反正切函数的泰勒级数得到；推出：=4 *Mathematica程序小数点后10000位数字为82、沃里斯(Wallis)方法Math