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

葉提芳

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

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

中圖分類號：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邊形得到小數(shù)點后的5位精度值,Ludolph Van Ceulen用262正邊形得到了小數(shù)點后的35位精度值。這種方法雖然經(jīng)典，但相當耗時。20世紀，很多數(shù)學家采取級數(shù)來近似計算圓周率的方法，已經(jīng)能把圓周率近似計算精確到了上億位，可以說，我們對圓周率的近似計算研究已經(jīng)相當成熟。本文試在利用Matlab實驗和高等數(shù)學中的知識有機結(jié)合起來，分別采取冪級數(shù)展開式的方法、隨機數(shù)的方法、數(shù)值積分的方法和公式法結(jié)合Matlab程序?qū)崿F(xiàn)對圓周率的近似計算。

1 利用冪級數(shù)展開式的計算方法

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

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

當我們?nèi)?= 時，得到一數(shù)項級數(shù)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 采取隨機數(shù)的近似計算方法

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

其中。

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

它位于長方體 = {（）|-1≤≤1}內(nèi)部，該長方體體積為8，往長方體內(nèi)部隨機投點個，然后統(tǒng)計錐形體內(nèi)的隨機點數(shù)，則≈ ≈8祝因為計算結(jié)果帶有隨機性，我們用十次重復的計算作比較，下面為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

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

3 利用定積分 = 近似計算

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

程度如下：

>> 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，可見，近似程度已經(jīng)很可觀了。

4 其它方法

我們知道 = ,= ,…= （重根號），即 = ，由此我們得到韋達公式 = …，據(jù)此= 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)

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

注釋

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

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