實驗二怎樣計算Pi

數(shù)學實驗實驗報告學院：數(shù)學與統(tǒng)計學院班級：數(shù)學與應用數(shù)學3班學號：4:康萍時間：2021.04.05實驗二怎樣計算丸一、實驗目的分別用以下三種方法計算丸的近似值，并比擬三種方法的準確度：數(shù)值積分法：通過使用Mathematica7.0編寫梯形公式和辛普森公式的程序語言計算丸。泰勒級數(shù)法：利用反正切函數(shù)泰勒級數(shù)計算丸。蒙特卡羅〔MonteCarlo〕法：通過使用Mathematica7.0編寫蒙特卡羅公式的程序語言來計算丸。二、實驗環(huán)境基于Windows環(huán)境下的Mathematica7.0軟件。三、實驗的根本理論和方法1、數(shù)值積分法以單位圓的圓心為原點建立直角坐標系，那么單位圓在第一象限的局部G是一個扇形，由曲線j=、1-x2(%e[0,1])及兩條坐標軸圍成，它的面積S=4。算出了S的近似值，它的4倍就是丸的近似值。而扇形面積S實際上就是定積分j1*1-%2d%=—。04與丸有關的定積分有很多，比方的定積分j11d%=-就比*m的+%201+%24定積分更容易計算，更適合于用來計算丸。一般地，要計算定積分jbfGM，也就是計算曲線j=f(%)與直線aj=0,%=a,%=b所圍成的曲邊梯形G的面積S。為此，用一組平行于y軸的直線%=%(1<i<n-1,a=x<%<%<<%<%=b)將曲邊梯形T分成n個小曲邊i012n-1n梯形，總面積S分成這些小曲邊梯形的面積之和。如果取「很大，使每個小曲邊

梯形的寬度都很小，可以將它上方的邊界f(x)X<x<x近似的看作直線段，將每個小曲邊梯形近似的看作梯形來求面積，就得到梯形公式。如果更準確些，將每個小曲邊梯形的上邊界近似的看作拋物線段，就得到辛普森公式。具體公式如下：梯形公式設分點X1,…,x1將積分區(qū)間［%b］分成n等份，即x.=a+iG-a)/n,0<i<n。所有的曲邊梯形的寬度都是h-G-a)/n。記y=f(x)那么第i個曲邊梯形的面積s近似的等于梯形面積上(y+y)h。將所有這些梯i2i-1i形面積加起來就得到y(tǒng)+y+…+y+yo+y”TOC\o"1-5"\h\z12n-12這就是梯形公式。辛普森公式仍用分點x.=a+i_aG<i<n—1)將區(qū)間偵b】分成n等份,直線x=xG<i<n-1)將曲邊梯形分成n個小曲邊梯形。再做每個小區(qū)間L,x】的ii-1i中點x.1i-2<x<x.)近似的看作經過三點(x,f(x)x=x中點x.1i-2<x<x.)近似的看作經過三點(x,f(x)x=xi-1,x,x.1ii-27的拋物線段，那么可求得r)y,-i+4yi+y,ki-2，于是得到其中y1=fxi-2s牝土6，于是得到s牝土6n(y+y)+2(y+y+…y0n12,)+4n-12、泰勒級數(shù)法

利用反正切函數(shù)的泰勒級數(shù)X3X5X2k-1TOC\o"1-5"\h\z\o"CurrentDocument"arctanx=x一一+—nV—1六-1H—52k-1當x的絕對值小于1,最好是遠小于1,這樣，隨著指數(shù)的增加，x的冪快速接近于0，泰勒級數(shù)就會很快收斂，比方，取X=1得到的arctan1就收斂的快，在\o"CurrentDocument"2++(—1)-12n-1++(—1)-12n-1中取2n-1=63得到的arctan1的近似值的誤差就小于上，準確度度已經非常高22651—,、.了。我們并不知道arctan是—的多少倍，但是卻能計異出a=arctan與一相差24那么—多少。記a=arctan,P=—一以，4tanp=tan(—-q]"4J丸tan—-tana41+tan—tan那么tanp=tan(—-q]"4J丸tan—-tana41+tan—tana41-121+1x132因此&.1=arctan—，31=arctan上，從而得到23——=arctan—+arctan-423arctan比arctan收斂得更快。利用泰勒級數(shù)計算出arctan與arctan的近似3423〔1〕值再相加，然后再乘以4，就得到—的近似值。還可以考慮用a=arctan5來計算—，它收斂的更快。由a=5易算出tan2a=Mtan4a120tanf4a-—"4Jtan4a-tan—也4_1191+tan4atan—1+12041239119—1-—=arctan239從而得到—=4a-arctan上4239=4arctan!-arctan51239TOC\o"1-5"\h\z即兀=16arctan-4arctan〔2〕239稱為Maqin公式，利用arctan工的泰勒展開式求出arctan^,ractan口一的近似值，5—=4a-arctan上4239=4arctan!-arctan51239由于丸是通過計算arctan工工=—,上」,口一等算出來的，只要計算這些"235239)arctanx的近似值到足夠的準確度，就能保證所得到的丸的近似值到足夠的準確度，我們是通過計算氣G)=x-1++(-*打來得到arctanx的近似值的。當0<x<1時，這個近似值的誤差為e=|arctanx-T(x)<；"+1]3、蒙特卡羅〔MonteCarlo〕法在數(shù)值積分中，我們利用求單位圓面積的1來得到，，從而得到丸，單位圓44的4是一個扇形G，它是邊長為1的單位正方形§的一局部，單位正方形§的面積S廣1。只要能夠求出扇形G的面積在正方形§的面積§中所占的比例k=—，就能立即得到S，從而得到丸的值。S1為了求出扇形面積在正方形面積中所占的比例k，一個方法是在正方形中隨機的投入很多點，使所投的每個點落在正方形每一個位置的時機均等，看其中有多少個點落在扇形。將落在扇形的點的個數(shù)m與所投點的總數(shù)n的比m可以看n作k的近似值。而任何一種計算機語言都有這樣的語言可以實現(xiàn)這樣的隨機點，能夠產生在區(qū)間k1]均勻分布的隨機數(shù)。在Mathematica中，產生區(qū)間t),」均勻分布的隨機數(shù)的語句是Random]產生兩個這樣的隨機數(shù)x,y，那么以〔x,y〕為坐標的點就是單位正方形的一點P,它落在正方形每個點的時機均等。P落在扇形的充分必要條件是x2+y2<1。這樣利用隨機數(shù)來解決問題的數(shù)學方法稱為蒙特卡羅法。四、實驗容與步驟及得到的結果分析實驗1數(shù)值積分計算兀1、實驗容分別取n=5000,n=6000,n=1000,用梯形公式和辛普森公式計算兀=j14dx01+x2的近似值，取20位有效數(shù)字，將所得的結果與兀進展比擬。2、實驗步驟在Mathematica中輸入語句如下：In[13]：=n=50a0：：=4/(1+x*jr)：si=(Suni[y[k/n]T{k；lrn-1)]+(y[0]+y[l])/2)/n;s2=<￥[?]+y[l]*{虹Ln-1)]+4wSum[y[(k-l/2)/n],{k,1.n}])/(6沖;Frint[(ff[sl^20Kff[s2,30],H[Pi,30]}]ln[17]:=n=6000;y[x_]:=4/'(1+xirjc)；sl=(Sum[r[k/n],(kr1.n-l>]-h(y[0]-ny[1])/2)/ji;s2=(y[0]+y[l]+2^Sun[y[k/n]f{k.1,-1}]+4*Sun[y[<k-1/2)/n],1』n}])/(6*n)：Print[(NEsl；20]；H[sZT30]TN[Pi；30]}]In[21]：=n=10BO：y[x_]:=4/(1+：si=(Sum[y[k/n],(k,Ln-1}]+C/[O]+y[U)/2)/n;s2=<￥[?]+y[l]42*Sum[y[k/3in{虬Ln-U3+4wSum[y[(k-l/2)/n],(k,1,n}])r(6wn);Print[{IT[S1.TM].IT[展.30]fM[Pi,30]}]3、實驗結果p.L4159264592312S571S,3.14159255350979323546264334360,.3.14159265358979323846264338328){3.1415926489601636088；3.14159265358979323846264336999；3.14159255358979323S4626433S328}{3.141592486923126571Sf3.14159265358979323846202334360f3.1415526535S979323S4626433S328}4、結果分析三次實驗結果所得的第一個數(shù)字是利用梯形公式計算出的丸，結果保存了20位有效數(shù)字，實驗結果所得的第二個數(shù)字是利用辛普森公式計算出的丸，結果保存了30位有效數(shù)字，第三個結果是丸的前30位有效數(shù)字組成的近似值。而且隨著「取值的增大，兀的準確度也隨之不斷增高。實驗2泰勒級數(shù)法計算丸1、實驗容利用反正切函數(shù)的泰勒級數(shù)：X3X5X2k-1TOC\o"1-5"\h\zarctanx=x一一+—nV—1六-1H—52k-1計算丸，分別將X=L,X=i,X=i,X=—代入上面的級數(shù)，并分別對「取值為235239n=5000,n=20000,n=40000計算丸的值，觀察所得的結果和所花的時間。2、實驗步驟在Mathematica中輸入語句如下In[55]：=T[吃,nJ:=Swik[(-1)1)/(2k+l)f{k,0,n>];ff[4*T[1,20000],20]//TimingT[x_,也|:=Sum[(-1}"KwxWK+1)/(2K+I),g0；n>];Print[K[4*(T[l/2.260]-nT[l./3,170]),150]];Print[R[16*CT[1/5.110]-4*T[1/239.30])^15H]]：Print[H[Pi,150]]郵1]：=T[x_.nJ:=Sinri[(-l)^k^x^(2k+1)/(2k+1).{kf",時];H[4*T[lr5000].20]//TimingT[x_,nJ:=Sum[(-1)^kwx^(2k+1)/(2k+1).(k,afn}];Print[N[4ir(T[l/2,260]+T[l/3,170]),150]];Print[II[16^(丁[1/5.110]-4frT[l/239,30])f150]];Print[II[Pi,250]]In[67]:=T[x_r《]:=SiDiL[(-l)-kirffA(2K+1}/(ZK+1),{K,0,n}];H[4*![!,40000]^M]//TimLngT[x_r心:=SiniL[(-1)(2k-F1)/(2k+1)f{k.0,n>]:Print[M[4n-CT[1/2,260]+T[1/3,LU]),150]]：Print[H[16t(T[l/5；110]-4irT[l/233；30])』150]];Print[H[Pi,250]]3、實驗結果Out[56]=(0.485,3.1416426510898659669}3.141592653589793238462643383279502884197169399375105820974944592307816■?.40628620899862S03482534211706?9821480865L32823066470938446095SOS8223L7'-.2535940SL32.35C543D334653033C655035C43572237571573423-5^30^171037713-53767519Q769Q5-.7O56QOS354O22O1CI764765247425OCI603621Q46526-59O4SL26394634O922191S7O32S1-.30031673143.141592653-5B979323B4626433S3279502884197169399375105820974944592307BL6■?.4O62B62O899B62SO34B25342117O67902148O065L32B23O6647O938446O955O5B223L7-.253.5940SL3Out[62]=(0.063^3.141792613595792B384}3.141592653509793230462643303279502884197169399375105020974944-592307016-.4Ofi2062O099B62BO34B2534211r7O67982148OS65132823O664,7a93044fiO95-5O582231'7-.2S35940S13890546093465309806590350485722375719734285480917183771353767619076909-.70-5S0083540220107847652474250068362L04652659048L2833463409221918703281'-.900316781414L592653-5S979323S4626433S32795O2S041971693993751Q532O974944-5923O7316-.4062862089986280348253421170679321480865132823066470938446095^05822317-.2535940812848111745028410270193852110555964462294895493036196442881097-.5665933446128475648233786783165271201909Out[68]=U-672,3.1416176523646049571}3.141592653-58979323S4626433S3279502SS41971633993751053203749445923Q7SL6\4062S62089986280348253421170679821460665132823066470936446095505622317■.2^359405132.S9054SO93465309805.590350485T22375719,7342S5480517LSS771353757819C76909-.7058008354022010784765247425006836210465265904612839463409221918703231■.90031678143.1415926SS-BBg'jgsa334626433S32^95023841971653993^5105320974944592307016-.4062862089986280348253421170679821480865133823066470938446095505622317■.2535940812S481LL74S0284L0270193852110555964462294895493038156442S81037\566593344612S47564S2337867831652712019094、結果分析在第一個實驗結果中，0.469指的是所用的時間是0.469s,后面的數(shù)值指的是取20位近似值所得出的兀的近似值。后面的三個數(shù)值，第一個數(shù)值是將X=1和21X=3代入所得的結果，結果保存了150位有效數(shù)字；第二個數(shù)值是將11X=尺和X=-—代入所得的結果，結果保存了150位有效數(shù)字；第三個數(shù)值是兀的前150位有效數(shù)字組成的近似值。第二個實驗比第一個實驗所用時間少，為0.047s，第三個實驗比前兩個實驗所用時間多，為1.672s。可見，隨著「取值的增多，所用時間會隨著增多，同時兀的準確度也會隨之增高。而且，通過比照數(shù)值積分法與泰勒級數(shù)法的計算結果，我們發(fā)現(xiàn)，當n的取值一樣時，相對于泰勒級數(shù)法而言，數(shù)值積分法的準確度較高。實驗3蒙特卡羅法計算兀1、實驗容取n=10000，n=5000,n=20000利用隨機投點的方法來計算兀的值。2、實驗步驟在Mathematica中輸入語句如下：In[?]：=hi=10000：p=4}:Bo[m=0：Do[x=RaiuLain.[];y=Rajuk?ni[];工瓦天以+丫以§I.m++]「{跖1,n}]：j^peiulTo[p,H[4m/n]K(t,1,10}]:Print[p];s?m[P[[t]].{七，LIn[47]:-n-5000?I>-EJfDo[in-0;Do[k-RcukcLoiTt-[]fy—TiRUiclom[]：Ifs■舊++],lr辿jgulT口[f,t<[4m/nJ]{七j比10}];Print[p];sumLpEtt]J,{t-,:L.1WJJIn[51]：=n.=20000：p={}：Do[n.=D;Do[x=RandiHiif];y=Raiukm[];If[xrt2+yA2v1,m+T,(k,1,nJ]:碩EiulTomH[4in/n]]；(t；lr10}];Print[p];Sun[r[[t]],(t,1,10}]/103、實驗結果{3.114r3.1292,3.1228,3.1804,3.1492.3.1284.3.1584.3.1476,3.1472,3.168}Out[42]=3-14452{3.1312,3.1288,3.1032,3.1752,3.11S4,3.1136,3.172^3.1152,3.1576,3.112}Out[5D]=3-13272{3.1462^3.141Z,3.1156,3.1432,3.1348r3.16r3.1374,3.1474,3.1534,3.126}Out[54]-3.140723、結果分析每次投10000個點得出兀的一個近似值存放在數(shù)組P中，一共做10次得到10個近似值，最后通過print[p]語句將其全部輸出得到:3.1264,3.156,3.1376,3.1616,3.1312,3.134,3.1624,3.1616,3.144,3.1448.最后求出這10個近似值的平均值，相當于隨機投點10000次得到的近似值，即3.14596。第二個積第三個實驗結果原理與第一個一樣，通過比照，我們發(fā)現(xiàn)，隨著n取值的增大，兀的準確度也越來越高。四三種計算方法的比擬通過上面的實驗，我們發(fā)現(xiàn)：當n=10000時精度很低，取更大的「，準確度會更高一些。但總的來說，蒙特卡羅法的準確度比數(shù)值積分和泰勒級數(shù)法低，基于此，我們依然用