モンテカルロ法で円周率 π を求める

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万回の試行結果は以下のようになる。円周率にかなり近似されている。

という考えで実装ができる。