【動手做做看】蒙地卡羅演算法-計算圓周率π

蒙地卡羅演算法是一種以機率統計理論為基礎的數值計算方法,主要是將求解的問題轉化為某種隨機分布的特徵數或問題本身具有隨機性,以程式模擬隨機事件出現的頻率並估計其機率,該方式通常用來解多維積分問題。雖然蒙地卡羅方法存在著精確度的問題,但該方法的求解中心思想是值得我們學習的。

1求圓周率π值

首先我們先複習一下底下兩個圖形的面積如何計算,假如有一個圓型半徑為1,那四分之一個圓面積就是半徑(r)乘以半徑(r)乘以π再除以4。

四分之一圓面積
= r * r * π / 4
= π / 4

在這個四分之一圓所在的正方形面積,就是1。

圓形面積
正方形面積

接下來我們來思考一個問題,我們已經知道四分之一圓面積及四分之一圓所在的正方形面積如何計算了,但假設π值我們並不知道實際數值是多少,這樣有什麼方法可以知道π值呢?

想像一個簡單遊戲,如果我們在上述正方形中隨機產生非常多的座標點,那有多少座標點會落在四分之一圓內呢?

假設總共產生了N個座標點,且有C個點落在四分之一圓內,則我們可以得到如下公式:
π / 4 : 1 = C : N
π = 4 * C / N

按照上述的公式推算,如果產生的座標可以平均的分散在每個地方,那麼只要知道N跟C的數值,我們就能知道π到數直為多少了。

2程式模擬

要知道N跟C的數值,除了我們可以實際玩個幾萬遍這個遊戲之外,還可以寫程式模擬來計算,其實程式碼就只有短短的9行,參考如下:

        
            var N = 30000;
            var C = 0;
            for (var i = 1; i < N; i++) {
                var point = {x: Math.random(), y: Math.random()}
                if (Math.pow(point.x, 2) + Math.pow(point.y, 2) < 1) {
                    C++;
                }
            }
            console.log(4 * C / N);
        
    

行5:這裡是判斷座標點是否落在四分之一圓內。

馬克有寫一個簡單的模擬程式,大家也可以透過以下連結玩看看。

3討論

相信您將此程式多跑個幾次應該會發現其實還蠻不精準的,這個在維基百科上面有一段解釋如下:

用蒙地卡羅方法近似計算圓周率的先天不足是:第一,電腦產生的亂數是受到儲存格式的限制的,是離散的,並不能產生連續的任意實數;上述做法將平面分割成一個個網格,在空間也不是連續的,由此計算出來的面積當然與圓或多或少有差距。

另外維基百科還提到即使是取10的9次方個隨機座標點,其結果也僅在前4位與圓周率吻合,不過馬克是沒有真的跑過10的9次方,如果大家有興趣可以自己試試看。

4附註

參考資料

Facebook 留言