pi的計(jì)算 實(shí)驗(yàn)報(bào)告

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

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

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

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

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

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

7、大數(shù)的乘除運(yùn)算時(shí)間由O(n2)縮短為O(nlog(n)。2、 Ramanujan公式1914年，印度數(shù)學(xué)家Srinivasa Ramanujan在他的論文里發(fā)表了一系列共14條圓周率的計(jì)算公式，這是其中之一。這個(gè)公式每計(jì)算一項(xiàng)可以得到8位的十進(jìn)制精度。1985年Gosper用這個(gè)公式計(jì)算到了圓周率的17,500,000位。3、AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式：初值： 重復(fù)計(jì)算： 最后計(jì)算：這個(gè)公式每迭代一次將得到雙倍的十進(jìn)制精度，比如要計(jì)算100萬位，迭代20次就夠了。1999年9月Takahashi和Kanada用這個(gè)算法計(jì)算

8、到了圓周率的206,158,430,000位，創(chuàng)出新的世界紀(jì)錄。4、Borwein四次迭代式：初值： 重復(fù)計(jì)算： 最后計(jì)算：這個(gè)公式由Jonathan Borwein和Peter Borwein于1985年發(fā)表，它四次收斂于圓周率。5、Bailey-Borwein-Plouffe算法這個(gè)公式簡稱BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995年共同發(fā)表。它打破了傳統(tǒng)的圓周率的算法，可以計(jì)算圓周率的任意第n位，而不用計(jì)算前面的n-1位。這為圓周率的分布式計(jì)算提供了可行性。1997年，F(xiàn)abrice Bellard找到了一個(gè)比BBP快40的

9、公式：第三部分：對于的幾種計(jì)算的研究和討論：1、數(shù)值積分法（I）利用積分公式計(jì)算n=10 ans =3.1045183262;n=20 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的圓稱為單位圓，它的面積等于。只要計(jì)算出單位圓的面積，就算出了。在坐標(biāo)軸上畫出以圓點(diǎn)為圓心，以1為半徑的單位圓（如下圖），則這個(gè)單

10、位圓在第一象限的部分是一個(gè)扇形，而且面積是單位圓的1/4，于是，我們只要算出此扇形的面積，便可以計(jì)算出。但計(jì)算量較大。2、數(shù)值積分法（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.141592

11、67442313;設(shè)分點(diǎn)x1,x2,xn-1將積分區(qū)間0,1分成n等分。所有的曲邊梯形的寬度都是h=1/n。記yi=f(xi).則第i個(gè)曲邊梯形的面積A近似地等于梯形面積，即：A=(y(i-1)+yi)h/2。將所有這些梯形面積加起來就得到：A2/n2(y1+y2+yn-1)+y0+yn3、利用復(fù)化梯形算法求Pi的近似值.=1/2 (Tn+h )4、泰勒級數(shù)法計(jì)算：利用反正切函數(shù)的泰勒級數(shù)來計(jì)算(I)x=1時(shí) =4* n = 10 ans= 3.04183961892940; n = 20 ans = 3.09162380666784; n = 50 ans = 3.1215946525910

12、1; n = 100 ans = 3.13159290355855; n = 200 ans = 3.13659268483882; n = 500 ans = 3.13959265558979; n = 1000 ans = 3.14059265383979; n = 2000 ans = 3.14109265362104;x=1時(shí)得到的的展開式收斂太慢,逼近速度太慢，運(yùn)算龐大，對速度造成了很大影響。5、泰勒級數(shù)法計(jì)算： 利用反正切函數(shù)的泰勒級數(shù)來計(jì)算(II)進(jìn)一步精細(xì)化 n = 10 ans= 3.1415925796063512110; n = 20 ans = 3.1415926535

13、897574098; 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;當(dāng)x=1時(shí)得到的的展開式收斂太慢。要使泰勒級數(shù)收斂得快，容易想到，應(yīng)當(dāng)使x的絕對值小于1，最好是遠(yuǎn)比1小。例如，因?yàn)?，所以我們可以?jì)算出的

14、值，從而得到的值。這樣，就使得收斂速度加快，逼近的速度大大增加.6、利用麥琴給出，推出=4() n = 10 ans= 3.1415926535897916969; 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.141592653589793238

15、5; n = 2000 ans =3.1415926535897932385;對泰勒級數(shù)，隨著x的減小，級數(shù)的收斂速度明顯加快，這啟示我們另外構(gòu)造相關(guān)級數(shù)來逼近。7、蒙特卡羅法計(jì)算單位圓的1/4是一個(gè)扇形，它是邊長為1的單位正方形的一部分，單位正方形的面積。只要能夠求出扇形的面積在正方形的面積中所占的比例，就能立即得到，從而得到的值。下面的問題歸結(jié)為如何求的值，這就用到了一種利用隨機(jī)數(shù)來解決此種問題的蒙特卡羅法，其原理就是在正方形中隨機(jī)的投入很多點(diǎn)，是所投的每個(gè)點(diǎn)落在正方形中每一個(gè)位置的機(jī)會(huì)均等，看其中有多少個(gè)點(diǎn)落在扇形內(nèi)。降落在扇形內(nèi)的點(diǎn)的個(gè)數(shù)與所投店的總數(shù)的比可以近似的作為的近似值。 n

16、= 100 ans= 3.200000000000000; n = 200 ans = 3.040000000000000; n = 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ù)模擬算法收斂的速度很慢，從運(yùn)行結(jié)果來看，蒙特卡

17、羅法的計(jì)算結(jié)果為3.108，雖然精確度不太高，但運(yùn)行時(shí)間短，在很多場合下，特別是在對精確度要求不高的情況下很有用的。除以上幾種方法之，還有下列一些的求其他的方法：*利用高斯公式=48+3220*、=2+1/3*(2+2/5*(2+3/7*(2+ （2+k/(2k+1)*(2+.)）).).(當(dāng) k=2799時(shí)可精確到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+.

18、(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ù)點(diǎn)后第10000位數(shù)字的幾種方法：1、反正切函數(shù)的泰勒級數(shù)得到；推出：=4 *Mathematica程序小數(shù)點(diǎn)后10000位數(shù)字為82、沃里斯(Wallis)方法Mathematica程序小數(shù)點(diǎn)后10000位數(shù)字為83、基于的級數(shù)麥琴給出：；推出=4()MATLAB程序小數(shù)點(diǎn)后10000位數(shù)字為8求取的方法很多，過程類似，不再累述。驗(yàn)證結(jié)論：的小數(shù)點(diǎn)后10000位數(shù)字為8第五部分：心得體會(huì)對于的值我們早已知曉，但是對于這個(gè)數(shù)值如何得到卻不清楚。通過這次學(xué)習(xí)，我們知道原來計(jì)算的值有如此多的方法，并且知道了級數(shù)在一些計(jì)算中的作用。通過本次實(shí)驗(yàn)，我們知道了運(yùn)用各種方法計(jì)算的近似值，知道了matlab和mathematica中關(guān)于級數(shù)的運(yùn)算，學(xué)習(xí)是一個(gè)枯燥的過程，有些東西是必須懂的。也就是一些基礎(chǔ)知識(shí)吧，包括基礎(chǔ)的操作和相關(guān)的理論知識(shí)。因此，有一本教材看