蒙地卡羅法(Monte Carlo Method)求圓周率的原理示意圖如下。正方形邊長為1單位長,面積為1平方單位;黃色扇形面積等於半徑為1單位長的1/4圓,面積為pi/4。在正方形內均勻隨機丟石頭,落在扇型內的機率 = 扇型面積÷正方形面積=pi/4。所以只要隨機產生N個座標(x,y),看看座標(x,y)落在扇形中(x2+y2<1)的次數有幾次。落在扇形中的次數除以N再乘上4的數值理論上就會接近圓周率PI。

蒙地卡羅法求圓周率

  程式中我們要不斷產生0<=x,y<1的座標點,使用frand()的方式來產生即可。完整程式碼如下。

C_code.GIF
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double throwPI(int N){
    int i, count;
    double x, y;
    
    for( count = 0, i = 0; i < N; ++i ){
        x = rand()/((double)RAND_MAX+1);
        y = rand()/((double)RAND_MAX+1);
        if( x*x + y*y < 1 ) ++count;
    }
    return 4.0 * count / N;
}

int main(void){
    int i;

    srand( time(NULL) );
    for( i = 10; i <= 10000000; i *= 10 )
        printf("%10d : %10.6lf\n", i, throwPI(i) );
}

  下面是執行結果。理論上每次執行都會有點不同,但趨勢應該是相同的,也就是N愈大,得到的結果越接近PI。從數學上可以推估答案的收斂誤差為1/(2*Sqrt(N))。

CMD.GIF
        10 :   3.600000
       100 :   3.240000
      1000 :   3.176000
     10000 :   3.131200
    100000 :   3.138960
   1000000 :   3.140680
  10000000 :   3.141832
 100000000 :   3.141683

第一次使用亂數就上手-目錄for C語言:

  第一次使用亂數就上手 - 1. 即用篇for C
  第一次使用亂數就上手 - 2. 基本篇for C
    基本篇範例一:蒙地卡羅法(Monte Carlo Method)求圓周率
    基本篇範例二:產生不重複樂透號碼
  第一次使用亂數就上手 - 3. 推廣篇for C
    產生常態分布亂數
    RAND_MAX的試煉:原理x限制x修正
  第一次使用亂數就上手 - 4. 原理篇for C
  第一次使用亂數就上手 - 5. 基本篇for C


創作者介紹
創作者 latinboy 的頭像
latinboy

阿賢的部落格

latinboy 發表在 痞客邦 留言(9) 人氣()


留言列表 (9)

發表留言
  • 阿旭
  • WOW!
    最近在研究一些亂數公式。
    你的這系列文章讓我受益非淺,萬分感謝!
    期待之後的幾篇!
  • 亂數很有趣 :)

    latinboy 於 2010/03/12 23:07 回覆

  • 小黃
  • 謝謝你讓我完成作業..泣...
    話說...我寫的圓周率隨時都會改變....
    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    #define N 32767
    int main(int argc, char *argv[])
    {srand(time(NULL));
    int sum = 0;
    int i;
    for(i = 1; i < N; i++) {
    double x = (double) rand() / RAND_MAX;
    double y = (double) rand() / RAND_MAX;
    if((x * x + y * y) < 1)
    sum++; }
    printf("PI = %f\n", (double) 4 * sum / N);
    system("PAUSE");
    return 0;
    }
  • 每次執行都不一樣才是正常的...

    latinboy 於 2010/07/24 16:36 回覆

  • 訪客
  • 請問一下,計算座標(x,y)落在扇形中(x2+y2<1)的次數有幾次的目的是什麼呢?這一句看不太懂。
  • 計算與原點的距離,小於1就是落在圓內

    latinboy 於 2011/04/09 09:24 回覆

  • GRIFFIN
  • 請問有VISUAL STUDIO C#2010的語法寫的嗎?
  • 訪客
  • 用vba寫了一個macro如下,可參考看看,順便看看電腦行不行,如要停止按
    control break<^C>
    '=====================================
    Sub throwPI()

    Dim i, count, n, n1, x, y As Double
    Dim AA, BB As String

    ActiveWorkbook.Sheets("Sheet1").Activate
    Range("A11").Select
    n1 = Cells(Rows.count, 1).End(xlUp).Row '找A欄最後列號
    AA = "A11"
    BB = "B" + CStr(n1)
    If n1 > 11 Then
    Range(AA & ":" & BB) = "" '清空
    End If

    'Range("A11:B65536").Select '設定不鎖定時可用

    Randomize
    count = 0
    n = 0
    x = 0
    y = 0


    For i = 1 To 500000000
    x = Rnd()
    y = Rnd()
    If (x * x + y * y) <= 1 Then '在第一象限四分圓內時
    count = count + 1
    If (count Mod 10000) = 0 Then '每壹萬次登記一次資料
    n = n + 1
    Cells(10 + n, 1) = i
    Cells(10 + n, 2) = (count / i) * 4 'pi=圓內/總次數*4
    End If

    End If

    If Range("A1") = "END" Then Exit For

    Next i

    End Sub
  • robert
  • 剛貼上文,順道回應樓上問題
    squart(x^2+y^2)=squart(r^2)
    由於r=1,並取掉squart
    是圓的數學式(x^2+y^2)=1
    也是畢氏定理,在圓上做任一直角三角形可知
  • 訪客
  • ((double)RAND_MAX+1)

    這一行請檢查,我認為不需要加一,直接 RAND_MAX 即可。
  • cctsao
  • 很棒的範例, 讚!!!