pi的計算 實驗報告

1、班級：9131138502 姓名：張海洋 學號:913113850232圓周率（選做）要求：先查資料看看古人是怎樣計算的，再對的各種計算方法進行研究和討論（收斂速度等），并給出不同算法算出的小數(shù)點后第10000位的數(shù)字是什么，你覺得該數(shù)字應該是多少？第一部分：圓周率簡介 圓周率是指平面上圓的周長與直徑之比 （ratio of the circumference of a circle to the diameter） 。用符號（讀音：pài）表示。中國古代有圓率、圓率、周等名稱。它是一個常數(shù)（約等于3.141592654）。它是一個無理數(shù)，即無限不循環(huán)小數(shù)。在日常生活中，通常都用3.

2、14代表圓周率去進行近似計算。而用十位小數(shù)3.141592654便足以應付一般計算。即使是工程師或物理學家要進行較精密的計算，充其量也只需取值至小數(shù)點后幾百個位。計算圓周率的方法 “歷史上一個國家所算得的圓周率的準確程度，可以作為衡量這個國家當時數(shù)學發(fā)展水平的指標?！?歷史上最馬拉松式的計算，其一是德國的Ludolph Van Ceulen，他幾乎耗盡了一生的時間，計算到圓的內(nèi)接正262邊形，于1609年得到了圓周率的35位精度值，以至于圓周率在德國被稱為Ludolph數(shù)；其二是英國的William Shanks，他耗費了15年的光陰，在1874年算出了圓周率的小數(shù)點后707位?？上?，后人發(fā)現(xiàn)

3、，他從第528位開始就算錯了。把圓周率的數(shù)值算得這么精確，實際意義并不大?，F(xiàn)代科技領域使用的圓周率值，有十幾位已經(jīng)足夠了。如果用Ludolph Van Ceulen算出的35位精度的圓周率值，來計算一個能把太陽系包起來的一個圓的周長，誤差還不到質(zhì)子直徑的百萬分之一。以前的人計算圓周率，是要探究圓周率是否循環(huán)小數(shù)。自從1761年Lambert證明了圓周率是無理數(shù)，1882年Lindemann證明了圓周率是超越數(shù)后，圓周率的神秘面紗就被揭開了。 在中國，公元263年前后，劉徽提出著名的 “割圓術”求出了比較精確的圓周率。他發(fā)現(xiàn)：當圓內(nèi)接正多邊形的邊數(shù)不斷增加后，多邊形的周長會越來越逼近圓周長，而多

4、邊形的面積也會越來越逼近圓面積。于是，劉徽利用正多邊形面積和圓面積之間的關系，從正六邊形開始，逐步把邊數(shù)加倍：正十二邊形、正二十四邊形，正四十八邊形，一直到正三七二邊形，算出圓周率等于三點一四一六，將圓周率的精度提高到小數(shù)點后第四位。在劉徽研究的基礎上，祖沖之進一步地發(fā)展，經(jīng)過既漫長又煩瑣的計算，一直算到圓內(nèi)接正24576邊形，而得到一個結(jié)論： 3.1415926 3.1415927 同時得到 的兩個近似分數(shù)：約率為227；密率為355113。他算出的 的8位可靠數(shù)字，不但在當時是最精密的圓周率，而且保持世界記錄九百多年。以致于有數(shù)學史家提議將這一結(jié)果命名為“祖率”。 現(xiàn)在的人計算圓周率, 多

5、數(shù)是為了驗證計算機的計算能力，還有，就是為了興趣。第二部分：古人計算圓周率方法 古人計算圓周率，一般是用割圓法。即用圓的內(nèi)接或外切正多邊形來逼近圓的周長。Archimedes用正96邊形得到圓周率小數(shù)點后3位的精度；劉徽用正3072邊形得到5位精度；Ludolph Van Ceulen用正262邊形得到了35位精度。這種基于幾何的算法計算量大，速度慢，吃力不討好。隨著數(shù)學的發(fā)展，數(shù)學家們在進行數(shù)學研究時有意無意地發(fā)現(xiàn)了許多計算圓周率的公式。下面挑選一些經(jīng)典的常用公式加以介紹。除了這些經(jīng)典公式外，還有很多其它公式和由這些經(jīng)典公式衍生出來的公式，就不一一列舉了1.Machin公式 這個公式由英國天

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

7、)算法。FFT可以將兩個大數(shù)的乘除運算時間由O(n2)縮短為O(nlog(n)。2、 Ramanujan公式1914年，印度數(shù)學家Srinivasa Ramanujan在他的論文里發(fā)表了一系列共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位，創(chuàng)出新的世界紀錄。4、Borwein四次迭代式：初值：        重復計算：         最后計算：    這個公式

9、由Jonathan Borwein和Peter Borwein于1985年發(fā)表，它四次收斂于圓周率。5、Bailey-Borwein-Plouffe算法這個公式簡稱BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995年共同發(fā)表。它打破了傳統(tǒng)的圓周率的算法，可以計算圓周率的任意第n位，而不用計算前面的n-1位。這為圓周率的分布式計算提供了可行性。1997年，F(xiàn)abrice Bellard找到了一個比BBP快40的公式：第三部分：對于的幾種計算的研究和討論：1、數(shù)值積分法（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、數(shù)值積分法

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將積分區(qū)間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、泰勒級數(shù)法計算：利用反正切函數(shù)的泰勒級數(shù)來計算(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、泰勒級數(shù)法計算： 利用反正切函數(shù)的泰勒級數(shù) 來計算(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時得到的的展開式收斂太慢。要使泰勒級數(shù)收斂得快，容易想到，應當使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;對泰勒級數(shù)，隨著x的減小，級數(shù)的收斂速度明顯加快，

16、這啟示我們另外構(gòu)造相關級數(shù)來逼近。7、蒙特卡羅法計算單位圓的1/4是一個扇形，它是邊長為1的單位正方形的一部分，單位正方形的面積。只要能夠求出扇形的面積在正方形的面積中所占的比例，就能立即得到，從而得到的值。下面的問題歸結(jié)為如何求的值，這就用到了一種利用隨機數(shù)來解決此種問題的蒙特卡羅法，其原理就是在正方形中隨機的投入很多點，是所投的每個點落在正方形中每一個位置的機會均等，看其中有多少個點落在扇形內(nèi)。降落在扇形內(nèi)的點的個數(shù)與所投店的總數(shù)的比可以近似的作為的近似值。 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;這種數(shù)據(jù)模擬算法收斂的速度很慢，從運行結(jié)果來看，蒙特卡羅法的計算結(jié)果為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第四部分：求的小數(shù)點后第10000位數(shù)字的幾種方法：1、反正切函數(shù)的泰勒級數(shù)得到；推出：=4 *Mathematica程序小數(shù)點后10000位數(shù)字為82、沃里斯(Wallis)方法Math