パソコンで数値計算

モンテカルロ法(積分計算)

モンテカルロ法

数学で学ぶ積分計算を用いて2次元上に描かれた閉曲線に囲まれた面積を求まる事ができますが、この場合閉曲線が具体的な関数で書ける事が必要です。もし関数で表せない複雑な閉曲線に囲まれた面積を求めるときは、モンテカルロ法を用いる方法があります。

モンテカルロ法で面積を求める方法のイメージを説明する前に、実際に身の回りにある物を使った面積を求める方法を紹介します。四角い容器と、閉曲線で囲まれた形をした容器を用意します。後者の容器を四角い容器の中に入れ、さらに雨が降っている外に設置します。しばらくたって容器に水がたまったら、たまった全体の水の重さMと、閉曲線の形をした容器にたまった水の重さmを量ります。四角い容器の面積をSとすると、閉曲線で囲まれた面積は、Sm/Mとなることがわかると思います。実際この方法で面積を求めるには単に水を注げばいいのですが、雨を使った考え方を利用したのがモンテカルロ法です。

モンテカルロ法はパソコンを使うので、さすがに雨は使えません。そこで雨の変わりに乱数を使います。ある求めたい積分領域を覆うことができる範囲に一様乱数を発生させます。求めたい積分領域に乱数が発生した場合カウントをとり、このカウント数と乱数を発生した回数との比をとります。この比の値と乱数を発生させた範囲の面積の値を掛ければ、求めたい積分領域の面積の値が出せます。またこの方法は3次元に拡張することも可能です。

モンテカルロ法を用いて積分計算をする場合、精度を上げるには非常に多くの質の高い乱数を発生させる必要があります。無限回発生すれば限りなく精度が高くなりますがコンピューターで乱数を発生する回数は有限の値となります。精度を上げるには非常に多くの乱数発生回数が必要ですが、積分領域の範囲が決める事が可能ならは比較的簡単に積分計算が行えます。

プログラムソース

ここではモンテカルロ法を用いて半径1の円の面積、すなわち円周率を求めてみます。考え方ですが、下図のように0<=x,y<=1の範囲の正方形内に乱数(下図では赤い点)を発生させます。そして中心原点半径1の 扇形の中に乱数が発生したらカウントします。そのカウント数を乱数の発生回数で割り、4培すれば円周率が求まります。



#include <stdio.h>
#include <stdlib.h>

int main()
{
        double i,imax,n;
        double x,y,pi;

        n=0.0;
        imax=10000000.0;        //      乱数の発生回数

        for(i=0;i<=imax;i++) {
/*0以上1未満の乱数を生成*/
                x=rand()/(RAND_MAX+1.0);
                y=rand()/(RAND_MAX+1.0);        

/*0<=x,y<=1の範囲にある中心原点半径1の
扇形の中に乱数による点が入ったらカウントする*/
                if((x*x+y*y)<1.0) {                     
                        n+=1.0;                                 
                }
        }

        pi=n/imax*4.0;          //半径1の円の面積すなわち円周率

        printf("%f\n",pi);

        return 0;
}


数値計算実行結果

上記のプログラムで求めた結果は

3.141131


となりました。乱数の発生回数を増やすほど、正確な円周率の値に近づくことがわかります。

2009/04/19 更新