基于R語言探求圓周率的估計(jì)值

李海

【摘 要】估算圓周率一直以來都是數(shù)學(xué)研究者感興趣的話題。本研究利用R語言，運(yùn)用已知的數(shù)學(xué)關(guān)系，設(shè)計(jì)常用的計(jì)算機(jī)模擬思路，估算π值，以探討圓周率的多種可行性估算方法，以期啟發(fā)、訓(xùn)練R語言編程思維。

【關(guān)鍵詞】R語言;實(shí)驗(yàn);圓周率

【中圖分類號(hào)】G642 ?【文獻(xiàn)標(biāo)識(shí)碼】A ?【文章編號(hào)】1671-8437（2020）04-0229-02

在小學(xué)階段學(xué)習(xí)圓的相關(guān)知識(shí)時(shí)，必然會(huì)提及圓周率，隨著進(jìn)一步的學(xué)習(xí)，發(fā)現(xiàn)有關(guān)圓周率π的學(xué)習(xí)幾乎貫穿了整個(gè)受教育過程，涉及的學(xué)科和用途也是十分廣泛。π取自希臘語，表示圓的周長(zhǎng)與直徑的比值。歷史上，許多中外數(shù)學(xué)家都試圖尋求圓周率π的近似值，而中國(guó)數(shù)學(xué)家劉徽和祖沖之采用幾何法也分別對(duì)圓周率進(jìn)行過估算，得出圓周率π的估計(jì)值，并一步一步將精確位數(shù)提升。2019年3月，谷歌利用現(xiàn)有計(jì)算機(jī)技術(shù)將圓周率精確到小數(shù)點(diǎn)后31.4萬億位。當(dāng)然也有很多民間數(shù)學(xué)愛好者借助計(jì)算機(jī)軟件技術(shù)來估算圓周率值，唐永鴿[1]用office辦公軟件excel、程登彪[2]和周輝[3]則利用數(shù)學(xué)軟件mathematica分別對(duì)圓周率進(jìn)行估算。

R語言作為一個(gè)自由、免費(fèi)、源代碼開放的軟件，是一個(gè)用于科學(xué)繪圖和統(tǒng)計(jì)計(jì)算的優(yōu)秀工具[4]。該軟件將數(shù)學(xué)思維與計(jì)算機(jī)程序設(shè)計(jì)重要思想完美結(jié)合。采用R語言軟件，根據(jù)不同數(shù)學(xué)估算公式估計(jì)圓周率的值，最終以圖形的方式呈現(xiàn)。

1 ? 估算方法

1.1 ?“劃分網(wǎng)格”估算法

“劃分網(wǎng)格”估算法：此方法類似著名的蒙特卡羅法，不同的是將正方形分成n2個(gè)小正方形，統(tǒng)計(jì)落在內(nèi)接圓內(nèi)的小正方形個(gè)數(shù)占所有正方形的比例，利用四分之一個(gè)圓來研究，看落入其中的小正方形有多少個(gè)，并計(jì)算概率，再用概率反推π的值，即π=4N/n2（N表示落入單位圓內(nèi)小正方形的個(gè)數(shù)），如圖1所示。

通過計(jì)算機(jī)模擬，運(yùn)行下列程序：

mianji=function（n）{

kk=matrix（0，n，n）

o=1/n

for（x in 1：n）{

for（y in 1：n）{

zb=c（x*o，y*o）

kk[x，y]=ifelse（sqrt（sum（zb^2））<=1，1，0）

}

}

sum（apply（kk， 1， sum））/n^2*4

}

pi=rep（0，1000）

for（s in 1：1000）{

pi[s]=mianji（s）}

plot（1：1000，pi，ylim = c（2.5，4））

abline（h=3.14，col=2）

隨著n的變大，π的值越趨于穩(wěn)定，如圖2：

1.2 ?“三角形”估算法

“三角形”估算法：利用已證明的結(jié)論，任意寫兩個(gè)小于1的數(shù)（x，y），將它和1組成一個(gè)數(shù)對(duì)（x，y，1），則由x、y和1作為線段長(zhǎng)度，能構(gòu)成一個(gè)鈍角三角形的概率為P=（π-2）/4。利用計(jì)算機(jī)選取1到1000的隨機(jī)數(shù)，并計(jì)算此概率，并反推π值，即π=4P+2。通過R軟件運(yùn)行如下程序：

pi=rep（0，1000）

for（s in 1：1000）{

i=rep（0，10000）

for（j in 1：10000）{

bianchang=runif（2，0，1）

i[j]=ifelse（sum（（bianchang）^2）<1&&sum（bianchang）>1，1，0）

}

pi[s]=4*mean（i）+2

}

plot（1：1000，pi，ylim = c（2.5，4））

abline（h=3.14，col=2）

模擬的π值大致分布如圖3：

1.3 ?“互質(zhì)”估算法

“互質(zhì)”估算法：利用數(shù)論和概率知識(shí)可知，任意兩個(gè)自然數(shù)互質(zhì)的概率為P=6/π2，即π=。利用計(jì)算機(jī)選取1到1000當(dāng)中的任意兩個(gè)自然數(shù)，并計(jì)算這兩個(gè)自然數(shù)互質(zhì)的概率，從而進(jìn)一步估算π值。通過R軟件運(yùn)行如下程序：

for（j in 1：1000）{r=rep（0，1000）

for（i in 1：1000）{

r[i]=husu（sample（1：1000，2））

}

z[j]=sum（r）

pi[j]=sqrt（6*1000/z[j]）}

plot（1：1000，pi，ylim = c（2.5，4））

abline（h=3.14，col=2）

模擬的π值大致分布如圖4：

2 ? 結(jié)論

隨著社會(huì)的高速發(fā)展，計(jì)算機(jī)的普及，越來越多的數(shù)學(xué)問題都可以采用計(jì)算機(jī)解決，再加上各種軟件的產(chǎn)生，給數(shù)學(xué)愛好者帶來了很多的樂趣，更是加強(qiáng)了各類學(xué)科之間的融合，有利于學(xué)習(xí)者開展創(chuàng)造性的工作，促進(jìn)整個(gè)科學(xué)事業(yè)的進(jìn)步。利用各種計(jì)算機(jī)軟件程序?qū)A周率推導(dǎo)公式進(jìn)行驗(yàn)證，得到了直觀的結(jié)果，如以上三種估算方法中。若將n值增大，圖形呈現(xiàn)將會(huì)更加明顯，這是人工算法無法比擬的。因此，將數(shù)學(xué)和計(jì)算機(jī)結(jié)合的學(xué)習(xí)方法，將是培養(yǎng)學(xué)習(xí)者核心素養(yǎng)的一個(gè)重要途徑，能充分顯示“數(shù)形結(jié)合”的魅力，體現(xiàn)“實(shí)踐是檢驗(yàn)真理的唯一標(biāo)準(zhǔn)”的哲學(xué)思想，更能真正體會(huì)數(shù)學(xué)之美，幫助學(xué)習(xí)者更好的學(xué)習(xí)，更快樂的學(xué)習(xí)。

【參考文獻(xiàn)】

[1]唐永鴿.用計(jì)算機(jī)模擬的方法估計(jì)圓周率的值[J].文理導(dǎo)航（中旬），2011（2）.

[2]程登彪.計(jì)算機(jī)模擬蒲豐投針試驗(yàn)近似計(jì)算圓周率研究[J].福建電腦，2012（11）.

[3]周輝.關(guān)于圓周率π的幾種計(jì)算方法[J].科技致富向?qū)В?011（5）5.

[4]薛毅等.R統(tǒng)計(jì)建模與R軟件.北京：清華大學(xué)出版社，2007.