pi的计算实验报告

1、实验报告班级：02姓名：张海洋学号：9圆周率（选做）要求：先查资料看看古人是怎样计算兀的，再对兀的各种计算方法进行研究和讨论（收敛速度等），弁给出不同算法算出.的小数点后第10000位的数字是什么，你觉得该数字应该是多少第一部分:圆周率简介圆周率是指平面上圆的周长与直径之比（ratioofthecircumferenceofacircletothediameter）。用符号冗（读音：pa）i表示。中国古代有圆率、圆率、周等名称。它是一个常数（约等于）。它是一个无理数，即无限不循环小数。在日常生活中，通常都用代表圆周率去进行近似计算。而用十位小数便足以应付一般计算。即使是工程师或物理学家要进行较

2、精密的计算，充其量也只需取值至小数点后几百个位。计算圆周率的方法“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。”历史上最马拉松式的计算，其一是彳惠国的LudolphVanCeulein他几乎耗尽了一生的时间，计算到圆的内接正262边形，于1609年得到了圆周率的35位精度值，以至于圆周率在彳惠国被称为Ludolph数；其二是英国的WilliamShanks,他耗费了15年的光阴，在1874年算出了圆周率的小数点后707位。可惜，后人发现，他从第528位开始就算错了。把圆周率的数值算得这么精确，实际意义并不大。现代科技领域使用的圆周率值，有十几位已经足够了

3、。如果用LudolphVanCeulenB出的35位精度的圆周率值，来计算一个能把太阳系包起来的一个圆的周长，误差还不到质子直径的百万分之一。以前的人计算圆周率，是要探究圆周率是否循环小数。自从1761年Lambert证明了圆周率是无理数，1882年Lindemann证明了圆周率是超越数后，圆周率的神秘面纱就被揭开了。在中国，公元263年前后，刘徽提出著名的“割圆术”求出了比较精确的圆周率。他发现：当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而多边形的面积也会越来越逼近圆面积。于是，刘徽利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加倍：正十二边形、正二十

4、四边形，正四十八边形，一直到正三。七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第四位。在刘徽研究的基础上，祖冲之进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正24576边形，而得到一个结论：<兀<同时得到冗的两个近似分数：约率为22/7;密率为355/113。他算出的冗的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。现在的人计算圆周率，多数是为了验证计算机的计算能力，还有，就是为了兴趣。第二部分：古人计算圆周率方法古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长

5、。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；LudolphVanCeuleM正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了1.Machin公式16arctg154argtg品3572n1xxxn1xarctg(x)xL(1)3572n1这个公式由英国天文学教授JohnMachin于1706年发现。他利用这个公式计算到

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

7、4、Ramanujan 公式980142 kw /3、1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是率的 17,500,000 位。其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周AGM(Arithmetic-Geometric Mean)算法Gauss-Legendre 公式：WEQ 二 1 二1重复计算：STWHHE HFW最后计算：1，这个公式每迭代一次将得到双倍的十进制精度，比如要计算 100万位，迭代20次就够了, 1999年9月Takahashi和Kanada用这个

8、算法计算到了圆周率的 206,158,430,000位，创出 新的世界纪录。Borwein四次迭代式:卜瑜镯细K鲍评舞嘲我脾乂 7祗蹄|+北做最后计算：国入初化总0-4调为审重复计算：卧1”）严这个公式由JonathanBorwein和PeterBorwein于1985年发表，它四次收敛于圆周率。5、Bailey-Borwein-Plouffe算法16n8n18n48n58n6这个公式简称BBP公式，由DavidBailey,PeterBorwein和SimonPlouffe于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布

9、式计算提供了可行性。1997年，FabriceBellard找到了一个比BBP快40%的公式：61召1024'（4门+14h+3+10n+lIOti+310n+d10n+7+10f?+95n=U第三部分:对于兀的几种计算的研究和讨论:1、数值积分法（I）利用积分公式4x/TTdx计算Commandwindow»n=lCOOO：inspace(Ojljir+1),y=sq.rt(1一%"2);h=l/n;ari5=vpsraps(y),11)ans-3.1415914770n=10ansn=20ansn=50ansn=100ansn=200ansn=500ansn=1

10、000ansn=2000ans半径为1的圆称为单位圆，它的面积等于。只要计算出单位圆的面积，就算出了在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4,于是，我们只要算出此扇形的面积，便可以计算出。但计算量较大。2、11数值积分法（II）利用公式42dx01xCommandWindow»n=50：i=D;1/n；1;s=0;fouk=i;lenetliti)-1s=s+(l/(l-K(i(i)+i(k-hl)V2)2)*1/;号ndans=3.141625班89230。n=10ans=;n=20ans=;n=50

11、ans=;n=100ans=;n=200ans=;n=500ans=;n=1000ans=;n=2000ans=;设分点Xi,X2Xn-l将积分区间0,1分成n等分。所有的曲边梯形的宽度都是h=1/n0记yi=f(xi).则第i个曲边梯形的面积A近似地等于梯形面积，即：A=(y(i-1)+yi)h/2将所有这些梯形面积加起来就得到：22n2(yi+y2+-yn-i)+yo+yn3、利用复化梯形算法求Pi的近似值.1 -Xfaihhl01x2=1/2(Tn+hi12)CommandWindowyyclear；a-0;b=lf=inliae4/(l-hc*x')；1=(b-aj/2*(&#

12、163;(a)+f(bD;cr=l;n=l;uhileer>l.Oe-6h(b-a)/n.5=0；fori=l;ns=s-H(a+i*h-h/2>endt2=(tl+h*£)/?；eu=abs(t2-t1);fprin+f(ofjr=9t.6£nt12ner);n=2*n;tl=t2;endendt-3.lOCOOO,r=0.100000t=5.131170Jr=C(031176t=3.13S9S8,r=0.007912：i=3.140942,r=0.001955t=3.141430,匚0.COC409t=3.141552,r=0.00C122t=3.141B8

13、2,r=0.000031t=3.141590,r=C.000003t=3.141592,r=0.000002t=3.141592,r=0.000000»4、泰勒级数法计算利用反正切函数的泰勒级数arctan x2k 1（1）k11h来计算(I)x=1 时1 -43 5n=4*171 k 12 k 1C3irnnandWindow»n=lQOO;>>5=0；fork=1：ne=s44«(-1)*(k+1)/(2*k'1)；endypa(Ej20)ans-J.1405926538397941350n=10ans=;n=20ans=n=50ans=;

14、n=100ans=;n=200ans=;n=500ans=;n=1000ans=;n=2000ans=;x=1时得到的arctanl的展开式收敛太慢，逼近速度太慢，运算庞大，对速度造成了很大影响。5、泰勒级数法计算：利用反正切函数的泰勒级数352k1xxk1xarctanxx(1)采计舁352k11 1(II)进一步精细化arctan1arctan-arctan-2 3Gemm4ndwindow»eymsn;f2=C-l)*(n-l)*(1/3)V(2tn-1);ans1-symsujii(f1;1,79),anssynsujnCf2,%I,79);ari£=vpa(4*(

15、51+52),20)ans-3 .H15926535897932385n=10ans=;n:=20ansn=二50ansn:=100ansn:=200ansn:=500ansn:=1000ansn=二2000ans当x=1时得到的arctanl的展开式收敛太慢。要使泰勒级数收敛得快，容易想到，11应当使x的绝对值小于1,取好是远比1小。例如，因为arctanlarctan-arctan-,所以我们可以计算出arctan工,arctan1的值，从而得到arctan1的值。这样，就使得收敛23速度加快，逼近的速度大大增加.6、利用攵今给出一4arctan114*Grrr，.11.1、arctan,

16、推出九=4(4arctanarctan)452395239CommandWindowI»syirisn；T1二(-1)"Cn.l)*(l/5)(2*ir-l)/(2»n-1);£2二(-1)、Cn-D*(1/233)*(2*n-l)/C2*n-l)；ansl=s7nsujn(flfn,1,9);ans2=s!fflk5uni(f2,n,79);an£=vpa(4*(4*ansl-ans2)j2Q)ans=3.141592663B8Q7932386n=10ans=;n=20ans=;n=50ans=;n=100ans=;n=200ans=;n=5

17、00ans=;n=1000ans=;n=2000ans=;对泰勒级数，随着1X1的减小，级数的收敛速度明显加快，这启示我们另外构造相关级数来逼近兀。7、蒙特卡罗法计算单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积S1。只要能够求出扇形的面积S在正方形的面积中所占的比例kS/S,就能立即得到S,从而得到的值。下面的问题归结为如何求k的值，这就用到了一种利用随机数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落在扇形内的点的个数m与所投店的总数n的比可以近似的作为k的近似

18、值。ICommandWindow»k=2000,»mrD;forn=1ikifrandU)"2+-rand(1)'1m=m+1;erulendvpa(4*VlJ)ans=3.1080000000000000000000000000000n=100ans=;n=200ans=;n=500ans=n=1000ans=;n=2000ans=;n=5000ans=;n=10000ans=;n=20000ans=;这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为,虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情况下很有

19、用的。除以上几种方法之，还有下列一些的求其他的方法:*利用高斯公式1”11=48arctan+32arctan20arctan* 、兀=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*2八3)+(1*3)/(2*4)*(1/(5*2八5)+* 、a(裙i)+1=0(欧拉公式，也称世界上最杰出的公式)* 、1+(1/2)八2+(1/3)八2+(1/4)八2+.(1/rf)/6=兀* 、1+(1/2)八4+(1/3)八4+(1/4)八4+.(1/rf)440兀* 、1+(1/2)

20、八6+(1/3)八6+(1/4)八6+.(1/rf)6/945兀* 、1+(1/2)八8+(1/3)八8+(1/4)八8+.(1/rf)8/S450兀* 、1+(1/2)八10+(1/3)八10+(1/4)八10+-.(1/0r110793555第四部分：求的小数点后第位数字的几种方法:1、反正切函3 X arctan x x 3数5X5的泰勒级(1)k12k 1X2k 1Mathematica 程序出：=n = In finity;上打得到一14推出：1 =4 *1)k(1)kk 02k 12k 1i=i ；回3.1415926535S"93Z3S426433S32"50

21、2S341971fi9339J7510552097444592307514052&6205962603402&343小数点后10000位数字为82、沃里斯(Wallis)方法c2244662133557c2k2k2k12k12k1Mathematica 程序r-rnS3ii丁.二 n s Infinity;5 = 2 P (2 t / C2上一1) t 2 1c / (2 J： + 1)；Ns, 10 0021二r :二 3,141552«53559"532324626433532795020641971«9399375L05E2C 974 344

22、S923C73146622 620 2 93250 3 4 3Z5S42-G未令名1*35y?y5U545U7irawr2UgTTU53B7UlJT745tUUI7T742BZr二三二三：”;肥匕 4 斐="5 二 4三1 二二二 4 2 3Emm：9S922 5«9S9Eaai59205600101fi5525fi 37567«小数点后10000位数字为83、基于arctanx的级数麦琴给出：414 arctan5.山 ，1推出兀=4( 4 arctan arctan5+1 .arctan;239 _1_ 239)MATLAB程序yyeynsii；f1=(-1)xtn-l)*(l/b)A(2*n-l)/(2*n-l);arx51=sjTnsm(fljns1irf),ms2=syikEun(f2,m1inf);ans=vpa(4*(4*anpl-ais2)；10002)STlE3.14159265355979323840264338327050285419?1693991592056。101比52563756芮&小数点后10000位数字为8求取的方法很多，过程类似，不再累述验证lr-<BPir10002j141592«535B9793Z3e462m3a3Z79028841