pi的计算试验报告

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

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

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

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

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

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

7、的乘除运算，要用FFT(Fast Fourier Transform由法。FFT可以将两个大数的乘除运算时间由 O(n2)缩短为O(nlog(n)。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最后

8、计算：1，这个公式每迭代一次将得到双倍的十进制精度，比如要计算 100万位，迭代20次就够了, 1999年9月Takahashi和Kanada用这个算法计算到了圆周率的 206,158,430,000位，创出 新的世界纪录。Borwein四次迭代式:卜瑜镯细K鲍评舞嘲我脾乂 7祗蹄|+北做最后计算：国入初化总0-4调为审重复计算：卧1”枷”这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。5、Bailey-Borwein-Plouffe 算法16n 8n 1 8n 4 8n 5 8n 6这个公式简称 BBP公式，由David Baile

9、y, Peter Borwein和Simon Plouffe于1995年共同发 表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到了一个比BBP快40%的公式：61召1024'（ 4门+14h+3 + 10n+l IOti+3 10n+d10n+7 + 10f?+95 n = U第三部分:对于兀的几种计算的研究和讨论:1、数值积分法（I）利用积分公式4 x/TTdx计算Command window» n=lCOOO：inspace(Oj lj ir+1),y=

10、sq.rt (1 一%"2);h=l/n;ari5=vpsr aps (y) ,11)ans -3. 1415914770n=10ansn=20ansn=50ansn=100ansn=200ansn=500ansn=1000ansn=2000ans半径为1的圆称为单位圆，它的面积等于。只要计算出单位圆的面积，就算出了在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第 一象限的部分是一个扇形，而且面积是单位圆的 1/4,于是，我们只要算出此扇形的面 积，便可以计算出 。但计算量较大。2、i 1数值积分法（II）利用公式42dx01 xCommand Window

11、» n=50 ：i=D; 1/n； 1;s=0;fou k=i; lenetliti)-1 s=s+(l/(l-K(i(i)+i(k-hl)V2) 2)*1/;号ndans =3.141625班89230。n=10ans =;n=20ans =;n=50ans =;n=100ans =;n=200ans =;n=500ans =;n=1000ans =;n=2000ans =;设分点Xi,X2Xn-1将积分区间0,1分成n等分。所有的曲边梯形的宽度都是h=1/n0记yi=f(xi).则第i个曲边梯形的面积A近似地等于梯形 面积，即：A=(y(i-1)+yi)h/2将所有这些梯形面积加

12、起来就得到：2 Zn2(y+y2+ "-1)+丫0+y13、利用复化梯形算法求Pi的近似值.1 - Xf a i h hl01 x2 =1/2 (Tn+h i 12 )Command Windowyy clear；a-0;b=l f=inliae4/( l-hc*x')；1= (b-aj /2* (£ (a) +f (bD ;cr=l;n=l ;uhile er>l. Oe-6 h(b-a)/n.5=0； for i=l;ns=s-H(a+i*h-h/2> ;endt2= (t l+h*£)/?； eu=abs (t2-t 1);fprin+f

13、 ( of j r=9t. 6£n t 12n er); n=2*n; tl=t2;endendt-3. lOCOOO, r=0.100000t=5. 131170J r=C( 031176t=3. 13S9S8, r=0. 007912：i=3.140942, r=0. 001955t=3. 141430,匚0. COC409t=3. 141552, r=0. 00C122t=3.141B82, r=0.000031t=3. 141590, r=C. 000003t=3.141592, r=0.000002t=3. 141592, r=0. 000000»4、泰勒级数法计

14、算利用反正切函数的泰勒级数arctan x2k 1（1）k11h来计算(I)x=1 时11143 5n=4*171 k 12 k 1C3irnnand Window» n=lQOO;>> 5=0；for k=1：ne=s44« (-1) * (k+1) / (2*k' 1)；endypa(Ej 20)ans -J. 1405926538397941350n = 10ans=;n =20ans =n =50ans =;n =100ans =;n =200ans =;n =500ans =;n =1000ans =;n =2000ans =;Gemm4nd

15、window» eyms n;f2=C-l) * (n-l) *(1/3) V (2tn-1);ans 1 - symsujii (f 1; 1, 79),anssynsujnCf2, % I, 79);ari£=vp a (4* (51+52), 20)ans -3. H15926535897932385n = 10ans=;n :=20ansn =二 50ansn :=100ansn :=200ansn :=500ansn :=1000ansn =二 2000ans当x=1时得到的arctanl的展开式收敛太慢。要使泰勒级数收敛得快，容易想到，1 1应当使x的绝对值小于

16、1,取好是远比1小。例如，因为arctan1 arctan- arctan-,2 3所以我们可以计算出arctan工,arctan1的值，从而得到arctan1的值。这样，就使得收敛 23速度加快，逼近的速度大大增加.6、利用攵今给出一4 arctan11r r r ， .11.1 、arctan,推出九-4(4arctan arctan )452395239Command WindowI » syiris n；T1二(-1)" Cn.l)*(l/5) (2*ir-l)/ (2»n-1);£2二(-1)、Cn-D* (1/233) * (2*n-l)/C

17、2*n-l)；ansl=s7nsujn (f lf n, 1,9);ans2=s!fflk5uni(f2,n, 79);an£=vpa(4* (4*ansl-ans2)j 2Q)ans =3.141592663B8Q7932386n = 10ans=;n = 20ans =;n = 50ans =;n = 100ans =;n = 200ans =;n = 500ans =;n = 1000ans =;n = 2000ans =;对泰勒级数，随着1 X 1的减小，级数的收敛速度明显加快，这启示我们另 外构造相关级数来逼近兀。7、蒙特卡罗法计算单位圆的1/4是一个扇形，它是边长为1的单

18、位正方形的一部分，单位正方形的面积S 1。只要能够求出扇形的面积S在正方形的面积中所占的比例k S/S,就能立即得到S,从而得到 的值。下面的问题归结为如何求k的值，这就用到了一种利用随机 数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投 的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落 在扇形内的点的个数m与所投店的总数n的比可以近似的作为k的近似值。I Command Window» k=2000,» mrD;for n= 1 ikif randU) " 2+-r and (1) '1m=m+1;er

19、ulendvpa(4*VlJ)ans =3. 1080000000000000000000000000000n = 100ans=;n = 200ans =;n = 500ans =n = 1000ans =;n = 2000ans =;n = 5000ans =;n = 10000ans =;n = 20000ans =;这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为, 虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情 况下很有用的。除以上几种方法之，还有下列一些的求其他的方法:*利用高斯公式1”11=48arctan +32arctan 2

20、0arctan* 、兀=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)八6+(1/3)八6+(1/4)八6+.(1/rf)6/945 兀* 、1+(1/2)八8+(1/3)八8+(1/

21、4)八8+.(1/rf)8/S450 兀* 、1+(1/2)八10+(1/3)八10+(1/4)八10+-.(1/0r110793555第四部分：求 的小数点后第位数字的几种方法:1、反正切 函 3Xarctan x x 3数5 X5的泰勒级(1)k12k 1X2k 1Mathematics 程序上打得到一14推出：1 =4 *1)k(1)kk 02k 12k 1i=i ；回3.1415926535S"93Z3S426433S32"502S341971fi9339J7510552097444592307514052&6205962603402&343小数点后

22、10000位数字为82、沃里斯(Wallis)方法c 2 24 46 621 33 55 7c2k2k2k1 2k 1 2k 1Mathematica 程序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 344S923C73146622 620 2 93250 3 4 3Z5S42-G未令名1* *35y?y5U545U7ir

23、awr2UgTTU53B7UlJT745tUUI7T742BZr二三二三：”;肥匕 4 斐="5 二 4三1 二二二 4 2 3Emm：9S922 5«9S9Eaai59205600101fi5525fi 37567«小数点后10000位数字为83、基于arctanx的级数麦琴给出：414 arctan51,推出兀=4( 4 arctan arctan5,1 . arctan;239 _1_ 239)MATLAB程序小数点后10000位数字为8求取的方法很多，过程类似，不再累述验证lr- < BPir 10 002j141592«535B9793Z

24、3e462m3a3Z79028841971693933-2g>49445铝配HLGM及萌20的9找之酊34醛53屯2；110 67 92214£0£5132E23(j647Q938460955C5:223172-qq耳4代4汽11 i,idAn?*di仆。7什i 口才白彳寸力与:白片工口书*?。；*日 回 於35鹫/颐970704号乳工口乳4。;二。1*41口口162=02109276457 州 1M57922955249犯7275后架 比126躺3 的黑3f?225c95965£150560010ie55256 3757e6结论：兀的小数点后10000位数字为8第五部分：心得体会对于兀的值我们早已知晓，但是对于这个数值如何得到却不清楚。通过这次学习，我们知道原来计算兀的值有如此多的方法，并且知道了级数在一些计 算中的作用。