モンテカルロ法を用いた円周率の近似値算出

標準

今日は先日暇つぶしでやってみた、モンテカルロ法を用いた円周率の近似値を求めるプログラムについて書いていきます。

モンテカルロ法とは

モンテカルロ法とは何かしらの数値計算を乱数を用いて近似値を得る手法の総称です。「モンテカルロ」という名前はカジノで有名な都市の名前からとられたそうです。このモンテカルロ法は円周率を求める以外にも「巡回セールスマン問題」などの数値計算も可能にする、使いによってはとても優れた計算手法です。ではまずはこれを使ってどう円周率の近似値を得るか考えていきます。

どうやって円周率を求めるか

まずは半径r=1、中心角θ=90°の扇形を描きます。この時の扇形の面積Sは、それをちょうど覆う正方形の面積 1×1=1 に対してS=pi/4になります。

上記の図の正方形の左下の頂点を座標平面の原点と定めて、点(x,y) (但し、0≦x≦1, 0≦y≦1をいずれも満たす)を乱数を用いて定めます。そうすると、試行回数Nに対して定めた点が扇形に含まれる回数Xの割合は、正方形に対する扇形の面積の割合に近似して、これを用いるとpiは以下のように近似させることが出来ます。

ちなみに、定めた点が扇形に含まれているかどうかは原点からの距離dが1以下であるか否かで判別することが出来ます。

実際に求める

今回は自分が一番得意なPHPを使って円周率の近似値を得ようと思います。

<?
//変数初期化
$count=0;
$start=1;
$max=100000000;
for($i=$start;$i<=$max;$i++) {
    //点を定める
    $x=(double)(mt_rand(0,1000000000)/1000000000);
    $y=(double)(mt_rand(0,1000000000)/1000000000);

    //距離を三平方の定理を用いて算出
    $d=$x*$x + $y*$y;

    //扇形に含まれる場合
    if($d <= 1.00000000) {
        $count++;
    }

    //結果出力
    print(sprintf("( %.9f , %.9f ) d = %.9f ( %d / %d ) pi = %.9f\n",$x,$y,$d,$count,$i,4.0*($count/$i)));
}
//最終結果出力
$pi=4.0*($count/$max);
print(sprintf("%.9f",$pi));
?>

PHPでは乱数(正確には「疑似乱数」)を生成する関数として、rand()mt_rand()random_int()などありますが、今回はより高速に、より良い(偏らない)乱数を生成出来るという特性をもつメルセンヌ・ツイスターを利用したmt_randを用いてみました。

早速これで1億回試行してみると以下のような値を取得することが出来ました。

pi = 3.141529200

(たまたま)小数点以下4桁まで正確に算出することが出来ました。この算出結果は、1億回の試行でも大きく外れてしまうこともありますし、精度良く求めることもあるという揺らめきの中での算出とはなりますが、いかにして精度良く算出できるか、極める価値がありそうですよね。

今回はこれで以上です。

投稿者プロフィール

yoshipc
コンピューター関連を得意としています。PHPが専門です(尚、技量はお察し)。このブログとMastodonのインスタンスを運営・管理しています。よろしくお願いいたします。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です