0. 概要
やりたいこと：
MonteCalro methodを用いて、円の円周率を求める。
動機：
大学で履修している計算物理の講義で登場してきたから、ここでその考え方についてまとめたいと思ったから。
1. 理論的なところ
半径 r の円を準備する
x^2 + y^2 = r^2
と記述できる円(中心が原点、半径 $r$ )を用意する。
この円を4分割したうちの1つを準備しよう。
切り取ったものを正方形にぶち込む
切り取った円を、以下のように一辺 r の正方形に収納することを考える。
以下、簡単のために r = 1 とする。
乱数を落として考える
ここで、
0 \le x \le 1 \\
0 \le y \le 1
で範囲を設定した乱数を発生させる。
(半径の任意の $r$ とするのであれば、乱数の範囲を $0\le x\le r$ までにしよう。)
この範囲内に乱数を $n$ 個落としたとする。そのうち、円の範囲内に $a$ 個落下したとする。すると、
正方形の面積：1/4円の面積 = $n : a$
r^2 : \frac{1}{4}\pi r^2 = n : a \\
\pi = \left(\frac{a}{n}\right)\times 4
というふうに、円周率を求めることができる。これを実装していこう。
2. 実装
Mntclo.cpp
#include <iostream>
#include <cstdllib>
#include <cmath>
using namespace std;
//乱数の定義
#define RAND()(float(rand())/(RAND_MAX + 1.0))
main() {
int i;
int ncout = 0;
int ntry;
float x , y , π;
cout << "ntry=?";
cin >> ntry;
for(i = 1; i <= ntry; ++i){
x = RAND();
y = RAND();
if (x*x + y*y <= 1) {
ncout = ncout + 1;
}
π = (float(ncout)/i)*4;
cout << "i=" << i << "π=" << π << "\n";
}
}
3.結果
1万回の試行結果は以下のようになる。円周率にかなり近似されている。
という考えで実装ができる。