モンテカルロ法による計算
シミュレーションプログラムの基礎として、モンテカルロ法を用いた円周率の計算を紹介します。
モンテカルロ法による円周率の計算
円周率 π の値は、解析的に求めることができません。 そこで、乱数を用いたモンテカルロシミュレーションで円周率を求めてみます。
半径１の円の 1/4 とそれに隣接する四角形を考えます。
四角形内にランダムに点を打ちます。
円の中心と打った点との距離から、打った点が円内にある(距離が１未満)か、円外にある(距離が１以上)かが分かります。
四角形と円の面積比は、全ての打った点の数と円内に打った点の数の比に一致します(完全にランダムで打った点の数が無限の場合)。
四角形の面積は１、円の面積はπ/4、打った点の数を n 、円内の点の数を m とすると、円周率は、4m/n になります。
ブラウザで動作させる JavaScript では、計算時間が長いと警告が発生します。 そこで、インターバルタイマーを使い、1000ポイントづつ経過を表示しながら、1000ステップ計算します。
var canvas = document.getElementById('jsccanvas');
var context = canvas.getContext('2d');
context.clearRect(0, 0, 640, 400);
function point(x, y, c) {
if(c) {
context.strokeStyle='rgb(255, 255, 0)';
} else {
context.strokeStyle='rgb(0, 255, 255)';
}
x = 400*x + 200;
y = 400 - 400*y;
context.beginPath();
context.moveTo(x, y);
context.lineTo(x+1, y+1);
context.stroke();
}
var n = 0;
var m = 0;
var step = 0;
function calcpi() {
var i;
for(i = 0; i < 1000; ++ i) {
var x = Math.random();
var y = Math.random();
var r = x*x + y*y;
if(r < 1.0) {
++ m;
point(x, y, true);
} else {
point(x, y, false);
}
++ n;
}
console.log(4.0*m/n);
}
jscsetintervalid = setInterval(function() {
if(step >= 1000) {
clearInterval(jscsetintervalid);
jscsetintervalid = null;
}
++ step;
calcpi();
}, 10);