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

モンテカルロ法という乱数シュミレーションで円周率を導き出してみます。

半径 r の円の面積は PI * r^2 ですよね。

半径 r を 1 とした円を中心から十文字に四等分した右上の部分を使って円周率を求めます。

円の中心から辺の長さ 1 の正方形を重ねます。

正方形の中には4等分された右上の部分の扇形の円がきっちり半径 1 の状態であります。

この時 x 軸、y 軸の 1 の位置は重ねた正方形の辺の長さ 1 になりますね。

よって、正方形の面積 : 扇の面積 = 1 : PI / 4 となります。

したがって PI を求める式は次のようになります。

PI = 4 * 扇の面積/正方形の面積

プログラムは乱数を使ってこの正方形の中にたくさんの針を落とします。

針の落ちた位置が扇の中か外かで円周率を求めていくという力技の方法です。

落とす針の数が多いほど精度は上がる理屈ですが処理時間も増えてしまいます。

とりあえず試してみましょう。 

 

jp\yucchi\montecarlo_pi\Montecarlo_PI.java
package jp.yucchi.montecarlo_pi;

import java.util.Random;

public class Montecarlo_PI {

    private static final int NEEDLE = 100000000;

    public static void main(String[] args) {
        int hitCount = 0;
        double x, y, pi;
        for (int i = 0; i < NEEDLE; i++) {
            x = Math.random();
            y = Math.random();
            if (Math.pow(x, 2.0) + Math.pow(y, 2.0) < 1.0) {
                hitCount++;
            }
        }
        pi = 4.0 * hitCount / NEEDLE;
        System.out.println("Mathクラスによる円周率は " + Math.PI);
        System.out.println("モンテカルロ法による円周率は " + pi);
    }
}

実行結果は次のようになりました。

run:
Mathクラスによる円周率は 3.141592653589793
モンテカルロ法による円周率は 3.14185944
構築成功 (合計時間: 9 秒)

何回かやって見たのですが小数点以下3位くらいまでの精度しかでませんでした。

精度を上げるにはこの処理を数回(数万回?)繰り返して平均をとればいいのかなぁ・・・?

Hatena タグ: