蒙地卡羅法(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()的方式來產生即可。完整程式碼如下。
#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))。
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
文章標籤
全站熱搜
WOW! 最近在研究一些亂數公式。 你的這系列文章讓我受益非淺，萬分感謝！ 期待之後的幾篇！
亂數很有趣 :)
謝謝你讓我完成作業..泣... 話說...我寫的圓周率隨時都會改變.... #include
#include
#include
#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;
}
每次執行都不一樣才是正常的...
我的問題是我COMPILER之後都不知哪裡出錯，可以幫我看看嘛~~真得很不好意思，麻煩你><謝謝你感激不盡 http://tw.knowledge.yahoo.com/question/question?qid=1610121303266
請問一下，計算座標(x,y)落在扇形中(x2+y2<1)的次數有幾次的目的是什麼呢？這一句看不太懂。
計算與原點的距離，小於1就是落在圓內
請問有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
剛貼上文,順道回應樓上問題 squart(x^2+y^2)=squart(r^2) 由於r=1,並取掉squart 是圓的數學式(x^2+y^2)=1 也是畢氏定理,在圓上做任一直角三角形可知
((double)RAND_MAX+1) 這一行請檢查，我認為不需要加一，直接 RAND_MAX 即可。
很棒的範例, 讚!!!