R軟件在數值計算中的應用
——以圓周率的近似計算為例

李洪毅，賀 陵，周盛康

(吉首大學數學與統計學院，湖南 吉首 416000)

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

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

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

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

第一，產生隨機數。首先產生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 圓周率π的估計值隨實驗次數n(≤10 000)變化的動態圖

圖3 基于buffon.needle函數Buffon投針實驗的動態模擬結果(n=50)

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

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

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

第一，產生隨機數。首先產生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 圓周率π的估計值隨實驗次數n(≤10 000)變化的動態圖

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

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

第一，產生隨機數。首先產生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 圓周率π的估計值隨實驗次數n(≤10 000)變化的動態圖

4 結語

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