李洪毅，賀 陵，周盛康

(吉首大學數(shù)學與統(tǒng)計學院，湖南 吉首 416000)

當前，在工程、系統(tǒng)科學、高新技術(shù)等諸多領(lǐng)域的數(shù)值計算都離不開統(tǒng)計軟件，常用的統(tǒng)計軟件有R、Matlab、SAS等。R軟件是一個開源、免費的統(tǒng)計軟件，具有強大的統(tǒng)計分析和數(shù)值計算功能[1-3]。以圓周率的近似計算為例，詳細介紹了R軟件在數(shù)值計算中的廣泛應(yīng)用。

1 基于Buffon投針實驗計算圓周率

圖1 Buffon投針實驗的幾何概型

基于R軟件在計算機上實現(xiàn)Buffon投針實驗并近似計算圓周率π的步驟為：

第一，產(chǎn)生隨機數(shù)。首先產(chǎn)生n個相互獨立的隨機變量θ，x的抽樣序列θi，xi，i=1，2，…，n，其中θi～U(0,θ),xi～U(0,a/2)。

基于R軟件將上述步驟編寫模擬程序Buffon.r如下：

Buffon<-function(n, l=0.8, a=1){

k<-0; i<-1; pai=rep(0,n); set.seed(666)

while(i<=n){

theta_i<-runif(1, 0, pi); x_i<-runif(1, 0, a/2)

if(x_i<=l/2*sin(theta_i))

k<-k+1

pai[i]=2*l*i/(k*a); i=i+1}

return(pai)}

圖2 圓周率π的估計值隨實驗次數(shù)n(≤10 000)變化的動態(tài)圖

圖3 基于buffon.needle函數(shù)Buffon投針實驗的動態(tài)模擬結(jié)果(n=50)

2 基于概率分析法計算圓周率

圖4 基于概率分析計算圓周率的幾何概型

基于上述結(jié)果計算圓周率π的近似值，具體步驟為：

第一，產(chǎn)生隨機數(shù)。首先產(chǎn)生n個相互獨立的隨機變量X與Y，X與Y的抽樣序列xi，yi，i=1, 2,...,n，其中xi～U(0,1),yi～U(0,1)。

基于R軟件將上述步驟編寫模擬程序PA.r如下：

PA<-function(n){

k<-0; i<-1; pai=rep(0,n); set.seed(666)

while(i<=n){

x_i<-runif(1);y_i<-runif(1)

if(x_i^2+y_i^2<1)

k<-k+1

pai[i]=4*k/i; i=i+1}

return(pai)}

圖5 圓周率π的估計值隨實驗次數(shù)n(≤10 000)變化的動態(tài)圖

3 基于Monte Carlo方法計算圓周率

基于上述結(jié)果計算圓周率π的近似值，具體步驟為：

第一，產(chǎn)生隨機數(shù)。首先產(chǎn)生n個相互獨立的隨機變量X，X的抽樣序列xi，i=1,2,…,n，其中xi～U(0,1)。

基于R軟件將上述步驟編寫模擬程序MC.r如下：

MC <-function(n){

i <-1; pai=rep(0,n); set.seed(66)

x <-runif(n)

for(i in 1:n)

pai[i]=4*sum(sqrt(1-x[1:i]^2))/i

return(pai)}

圖6 圓周率π的估計值隨實驗次數(shù)n(≤10 000)變化的動態(tài)圖

4 結(jié)語

以圓周率的近似計算為例，詳細介紹了R軟件在數(shù)值計算和數(shù)值模擬中的具體應(yīng)用，由文中實例可以發(fā)現(xiàn)：R軟件能夠非常高效、便捷地解決數(shù)值計算和數(shù)值模擬中的近似計算問題。教學實踐和研究實踐證明，R軟件可以為統(tǒng)計學、數(shù)學專業(yè)課程教學和科學研究提供有力支撐，一方面可以加深學生對基本概念和算法的理解，更好地掌握數(shù)學理論方法，另一方面還可以通過R軟件優(yōu)秀的圖形功能加強計算結(jié)果的展示，使學生加深對數(shù)學理論方法的理解，更好地處理實際問題，激發(fā)學生學習興趣和動力，為學生更好地應(yīng)用專業(yè)知識處理實際問題奠定堅實基礎(chǔ)，對其今后的工作和學習產(chǎn)生積極作用。