0%

蒙地卡羅方法求圓周率

蒙地卡羅方法簡介

隨機點10個點, 這些點可能在原湖內,也有可能在圓弧外,例如有10個點,會有7個點在圓內,機率是0.7
pi-1
如果將每個點看作有一定大小的面積的話,則點在圓弧內的概率可以接近途中圓弧的面積
所以只要將概率x4就可以得出近似圖形的面積.所以將概率x4就可以得出圖形面積
但是如何判斷一個點在圓弧內還是外?根據畢氏定理,x^2+y^2小於1在圓內,大於1在圓外
所以我們只要一直弄出亂數放置隨機點去統計,透過佔總點比例和生成範圍就可以求出pi了

上code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <iostream>
#include <stdlib.h>
#include <math.h>
using namespace std;

int main()
{
int a = 0;
int n = 30000;
double x, y, pi;
for(int i = 0; i < n; i++)
{
x = (double)rand()/RAND_MAX; //隨機生成x座標
y = (double)rand()/RAND_MAX; //隨機生成y座標
if (pow(x,2)+pow(y,2) <= 1) //如果x^2+y^2<=1計數器增加一次,共運行n次
{
a++;//如果在扇形之內增加一次
}
pi = 4* (double)a / (double)i; //計算出pi
cout << pi << endl; //輸出
}
return 0;

延伸思考

那我們要不要來看看蒙地卡羅方法求出的pi和真實的pi誤差值的變化?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <stdlib.h>
#include <math.h>
using namespace std;

int main()
{
int a = 0;
int n = 10000;
double x, y, pi,err;
for(int i = 0; i < n; i++)
{
x = (double)rand()/RAND_MAX; //隨機生成x座標
y = (double)rand()/RAND_MAX; //隨機生成y座標
if (pow(x,2)+pow(y,2) <= 1) //如果x^2+y^2<=1計數器增加一次,共運行n次
{
a++;//如果在扇形之內增加一次
}
pi = 4* (double)a / (double)i; //計算出pi
err = fabs(pi - M_PI);//計算誤差值
cout << err << endl; //輸出誤差值
}
return 0;
}

使用gnuplot繪製圖表

1
2
3
4
5
gnuplot> set xlabel "i" //設定x軸名稱為i
gnuplot> set ylabel "error" //設定y軸名稱為error
gnuplot> set logscale xy //設定雙對數座標軸
gnuplot> set yrange [1:0.00003] //設定y軸範圍1到0.00003
gnuplot> plot "pi.dat" with lines[1] //使用連線描圖

gnuplot