利用Matlab近似計算圓周率的若干方法

葉提芳

摘 要 本文分別采取冪級數展開式的方法、隨機數的方法、數值積分的方法和公式法結合Matlab程序實現對圓周率的近似計算，分析實驗結果，比較每種方法的近似程度的高低，實現了Matlab實驗和數學理論的很好結合。

關鍵詞 Matlab實驗 圓周率 冪級數 隨機數 數值積分

中圖分類號：TP312 文獻標識碼：A

Some Methods of Approximatively Calculating Using Matlab Experiments

YE Tifang

(Industrial and Commercial College, Wuhan Polytechnic University, Wuhan, Hubei 430065)

Abstract In this paper, we used expanding power series, random number, numerical integration and formula methods combining matlab experiments to achieve the approximate value of. Then, we analyzed the experiments results, and compared the degree approximation of every method. It achieved the satisfying results of combination withMatlab experiment and mathematics theory.

Key words Matlab experiments; ; power series; random number; numerical integration

古今中外，歷史上有許多人積極致力于圓周率的研究與計算。我國的劉徽用正3072邊形得到小數點后的5位精度值,Ludolph Van Ceulen用262正邊形得到了小數點后的35位精度值。這種方法雖然經典，但相當耗時。20世紀，很多數學家采取級數來近似計算圓周率的方法，已經能把圓周率近似計算精確到了上億位，可以說，我們對圓周率的近似計算研究已經相當成熟。本文試在利用Matlab實驗和高等數學中的知識有機結合起來，分別采取冪級數展開式的方法、隨機數的方法、數值積分的方法和公式法結合Matlab程序實現對圓周率的近似計算。

1 利用冪級數展開式的計算方法

設是以2 為周期的周期函數，在[- , )上的表達式如下：

顯然，為奇函數，利用我們在高等數學冪級數一章知識①，可以將其展開為正弦級數為：= 4

當我們取 = 時，得到一數項級數4，因為，則有4 =，求其部分和為 = 4，顯然有 =，下面我們分別給出 = 1000,10000,50000時，利用Matlab程序求得的 的近似值。程序如下：

>> s=0;

>> n=50000;

>> digits(22)

>> for k=1:n

s=s+(-1)^(k-1)/(2*k-1);

end

>> s=vpa(4*s,20)

s =

3.1415726535897814387

= 1000,10000,50000時， 的近似值如下：

2 采取隨機數的近似計算方法

設一制作均勻的冰激凌可以看做是由圓錐面 = 和球面 ++ （）2= 1圍成②。我們利用積分知識求得它的體積為：

其中。

我們還可以采取隨機數的方法，由于所求錐形體可表示為：= {（）| + ≤ （）,+ ≤}，

它位于長方體 = {（）|-1≤≤1}內部，該長方體體積為8，往長方體內部隨機投點個，然后統計錐形體內的隨機點數，則≈ ≈8祝因為計算結果帶有隨機性，我們用十次重復的計算作比較，下面為N取100000時，用Matlab近似 的程序：

>> for k=1:10

r=rand(100000,3);

x=2*r(:,1)-1;

y=2*r(:,2)-1;

z=2*r(:,3);

fl=x.^2+y.^2;

p(k)=8*sum(fl-z.*z<=0&fl-z.*(2-z)<=0)/100000;

end

>> p

運行程序，得到 的近似值如下：

3.11363.12003.10003.16963.15603.12483.14643.09203.20003.1752

當取5000000時，運行程序，得到的 的近似值如下，可見近似程度是較好的：

3.14293.14293.14373.14053.14183.13973.13983.14093.14323.1430

從計算結果看：這種方法雖然簡單可行，但收斂的速度慢，距離真實值誤差較大。

3 利用定積分 = 近似計算

根據定積分的定義，（積分的結果和區間的分法及的取法無關，現在采取特殊分割和特殊取法不影響結果）將區間[0,1]分成等份，在每個小區間上，選取中點為，有

程度如下：

>> n=1000;

>> i=0:1/n:1;

>> s=0;

>> for k=1:length(i)-1

s=s+(1/(1+((i(k)+i(k+1))/2)^2))*1/n;

end

>> vpa(4*s,20)

運行程序，得到的 的近似值為3.1415927369231306798，可見，近似程度已經很可觀了。

4 其它方法

我們知道 = ,= ,…= （重根號），即 = ，由此我們得到韋達公式 = …，據此= 2/（…），

我們編寫Matlab程序如下（下面為計算式分母中取前10項乘積的程序）:

a=sqrt(2);

>> s=1;

>> for i=1:10

s=(s*a)/2;

a1=sqrt(2+a);

a=a1;

end

>> vpa(2/s,20)

運行程序，得到結果：3.1415914215111997443。當計算式分母中取前100項乘積時，運行程序，得到結果：3.1415926535897932385，可見，這種方法收斂速度快，近似精度高。

注釋

① 高純一,周勇.高等數學[M].上海:復旦大學出版社,2006:249-259.

② 李繼成.數學實驗[M].北京:高等教育出版社,2006:92.