pi的計算實驗報告

1、實驗報告班級：02姓名：張海洋學(xué)號：9圓周率（選做）要求：先查資料看看古人是怎樣計算兀的，再對兀的各種計算方法進行研究和討論（收斂速度等），弁給出不同算法算出.的小數(shù)點后第10000位的數(shù)字是什么，你覺得該數(shù)字應(yīng)該是多少第一部分:圓周率簡介圓周率是指平面上圓的周長與直徑之比（ratioofthecircumferenceofacircletothediameter）。用符號冗（讀音：pa）i表示。中國古代有圓率、圓率、周等名稱。它是一個常數(shù)（約等于）。它是一個無理數(shù)，即無限不循環(huán)小數(shù)。在日常生活中，通常都用代表圓周率去進行近似計算。而用十位小數(shù)便足以應(yīng)付一般計算。即使是工程師或物理學(xué)家要進行較

2、精密的計算，充其量也只需取值至小數(shù)點后幾百個位。計算圓周率的方法“歷史上一個國家所算得的圓周率的準確程度，可以作為衡量這個國家當時數(shù)學(xué)發(fā)展水平的指標?！睔v史上最馬拉松式的計算，其一是彳惠國的LudolphVanCeulein他幾乎耗盡了一生的時間，計算到圓的內(nèi)接正262邊形，于1609年得到了圓周率的35位精度值，以至于圓周率在彳惠國被稱為Ludolph數(shù)；其二是英國的WilliamShanks,他耗費了15年的光陰，在1874年算出了圓周率的小數(shù)點后707位?？上В笕税l(fā)現(xiàn)，他從第528位開始就算錯了。把圓周率的數(shù)值算得這么精確，實際意義并不大?，F(xiàn)代科技領(lǐng)域使用的圓周率值，有十幾位已經(jīng)足夠了

3、。如果用LudolphVanCeulenB出的35位精度的圓周率值，來計算一個能把太陽系包起來的一個圓的周長，誤差還不到質(zhì)子直徑的百萬分之一。以前的人計算圓周率，是要探究圓周率是否循環(huán)小數(shù)。自從1761年Lambert證明了圓周率是無理數(shù)，1882年Lindemann證明了圓周率是超越數(shù)后，圓周率的神秘面紗就被揭開了。在中國，公元263年前后，劉徽提出著名的“割圓術(shù)”求出了比較精確的圓周率。他發(fā)現(xiàn)：當圓內(nèi)接正多邊形的邊數(shù)不斷增加后，多邊形的周長會越來越逼近圓周長，而多邊形的面積也會越來越逼近圓面積。于是，劉徽利用正多邊形面積和圓面積之間的關(guān)系，從正六邊形開始，逐步把邊數(shù)加倍：正十二邊形、正二十

4、四邊形，正四十八邊形，一直到正三。七二邊形，算出圓周率等于三點一四一六，將圓周率的精度提高到小數(shù)點后第四位。在劉徽研究的基礎(chǔ)上，祖沖之進一步地發(fā)展，經(jīng)過既漫長又煩瑣的計算，一直算到圓內(nèi)接正24576邊形，而得到一個結(jié)論：<兀<同時得到冗的兩個近似分數(shù)：約率為22/7;密率為355/113。他算出的冗的8位可靠數(shù)字，不但在當時是最精密的圓周率，而且保持世界記錄九百多年。以致于有數(shù)學(xué)史家提議將這一結(jié)果命名為“祖率”。現(xiàn)在的人計算圓周率，多數(shù)是為了驗證計算機的計算能力，還有，就是為了興趣。第二部分：古人計算圓周率方法古人計算圓周率，一般是用割圓法。即用圓的內(nèi)接或外切正多邊形來逼近圓的周長

5、。Archimedes用正96邊形得到圓周率小數(shù)點后3位的精度；劉徽用正3072邊形得到5位精度；LudolphVanCeuleM正262邊形得到了35位精度。這種基于幾何的算法計算量大，速度慢，吃力不討好。隨著數(shù)學(xué)的發(fā)展，數(shù)學(xué)家們在進行數(shù)學(xué)研究時有意無意地發(fā)現(xiàn)了許多計算圓周率的公式。下面挑選一些經(jīng)典的常用公式加以介紹。除了這些經(jīng)典公式外，還有很多其它公式和由這些經(jīng)典公式衍生出來的公式，就不一一列舉了1.Machin公式16arctg154argtg品3572n1xxxn1xarctg(x)xL(1)3572n1這個公式由英國天文學(xué)教授JohnMachin于1706年發(fā)現(xiàn)。他利用這個公式計算到

6、了100位的圓周率。Machin公式每計算一項可以得到位的十進制精度。因為它的計算過程中被乘數(shù)和被除數(shù)都不大于長整數(shù)，所以可以很容易地在計算機上編程實現(xiàn)。還有很多類似于Machin公式的反正切公式。在所有這些公式中，Machin公式似乎是最快的了。雖然如此，如果要計算更多的位數(shù)，比如幾千萬位，Machin公式就力不從心了。下面介紹的算法，在PC機上計算大約一天時間，就可以得到圓周率的過億位的精度。這些算法用程序?qū)崿F(xiàn)起來比較復(fù)雜。因為計算過程中涉及兩個大數(shù)的乘除運算，要用FFT(FastFourierTransform由法。FFT可以將兩個大數(shù)的乘除運算時間由O(n2)縮短為O(nlog(n)。

7、4、Ramanujan 公式980142 kw /3、1914年，印度數(shù)學(xué)家Srinivasa Ramanujan在他的論文里發(fā)表了一系列共14條圓周率的計算公式，這是率的 17,500,000 位。其中之一。這個公式每計算一項可以得到8位的十進制精度。1985年Gosper用這個公式計算到了圓周AGM(Arithmetic-Geometric Mean)算法Gauss-Legendre 公式：WEQ 二 1 二1重復(fù)計算：STWHHE HFW最后計算：1，這個公式每迭代一次將得到雙倍的十進制精度，比如要計算 100萬位，迭代20次就夠了, 1999年9月Takahashi和Kanada用這個

8、算法計算到了圓周率的 206,158,430,000位，創(chuàng)出 新的世界紀錄。Borwein四次迭代式:卜瑜鐲細K鮑評舞嘲我脾乂 7祗蹄|+北做最后計算：國入初化總0-4調(diào)為審重復(fù)計算：臥1”）嚴這個公式由JonathanBorwein和PeterBorwein于1985年發(fā)表，它四次收斂于圓周率。5、Bailey-Borwein-Plouffe算法16n8n18n48n58n6這個公式簡稱BBP公式，由DavidBailey,PeterBorwein和SimonPlouffe于1995年共同發(fā)表。它打破了傳統(tǒng)的圓周率的算法，可以計算圓周率的任意第n位，而不用計算前面的n-1位。這為圓周率的分布

9、式計算提供了可行性。1997年，F(xiàn)abriceBellard找到了一個比BBP快40%的公式：61召1024'（4門+14h+3+10n+lIOti+310n+d10n+7+10f?+95n=U第三部分:對于兀的幾種計算的研究和討論:1、數(shù)值積分法（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數(shù)值積分法（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=;設(shè)分點Xi,X2Xn-l將積分區(qū)間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、利用復(fù)化梯形算法求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、泰勒級數(shù)法計算利用反正切函數(shù)的泰勒級數(shù)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、泰勒級數(shù)法計算：利用反正切函數(shù)的泰勒級數(shù)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的展開式收斂太慢。要使泰勒級數(shù)收斂得快，容易想到，11應(yīng)當使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=;對泰勒級數(shù)，隨著1X1的減小，級數(shù)的收斂速度明顯加快，這啟示我們另外構(gòu)造相關(guān)級數(shù)來逼近兀。7、蒙特卡羅法計算單位圓的1/4是一個扇形，它是邊長為1的單位正方形的一部分，單位正方形的面積S1。只要能夠求出扇形的面積S在正方形的面積中所占的比例kS/S,就能立即得到S,從而得到的值。下面的問題歸結(jié)為如何求k的值，這就用到了一種利用隨機數(shù)來解決此種問題的蒙特卡羅法，其原理就是在正方形中隨機的投入很多點，是所投的每個點落在正方形中每一個位置的機會均等，看其中有多少個點落在扇形內(nèi)。降落在扇形內(nèi)的點的個數(shù)m與所投店的總數(shù)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=;這種數(shù)據(jù)模擬算法收斂的速度很慢，從運行結(jié)果來看，蒙特卡羅法的計算結(jié)果為,雖然精確度不太高，但運行時間短，在很多場合下，特別是在對精確度要求不高的情況下很有

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第四部分：求的小數(shù)點后第位數(shù)字的幾種方法:1、反正切函3 X arctan x x 3數(shù)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小數(shù)點后10000位數(shù)字為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«小數(shù)點后10000位數(shù)字為83、基于arctanx的級數(shù)麥琴給出：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芮&小數(shù)點后10000位數(shù)字為8求取的方法很多，過程類似，不再累述驗證lr-<BPir10002j141592«535B9793Z3e462m3a3Z79028841