モンテカルロでπを計算
ある点から一定の距離にある平面上の点の集合を円といいます。この円の直径と円周の割合が円周率πなのですが、3.1415926535...ずっと続く不思議な数字です。このπの値を求めようと古くはギリシャのアルキメデスや、ドイツのルドルフ、日本の関孝和などが挑戦してきました。
一例を挙げると、アルキメデスは直径が1の円の内側にピッタリはまる正6角形(内接正6角形)の外周の長さと、円の外側にピッタリはまる正6角形(外接正6角形)の外周の長さを求めて、πはその間にあるとしました。
そして次は正12角形、次は正24角形…。そして、最後は正96角形までにして計算していきました。その最終的な式が以下になります。
3 ---- < π < 3 ----
71 70
πは左辺と右辺の間にあるというわけです。rpnで計算してみます。
3.14085 3.14286
なんと、この方法でアルキメデスは3.14まで正確に求めていました。その後、これと同様な方法で16世紀のヨーロッパにおいては、πの近似値として次の分数に至っていました。
π ≒ -----
113
この分数をrpnで計算すると…
3.14159
もはや実用には何の問題もない精度ですね。
πを求めるための2つの公式
さて、ここではアルキメデスが使った幾何学的な方法ではなくて、ちょっと変わった方法でπの計算をしてみます。モンテカルロ法というやり方です。これはでたらめな数(乱数)を使って数値計算するときに使う呼び名です。カジノの都モナコのモンテカルロに因んで付けられた名前です。
でたらめな数、乱数を使うことでπが求まるのかですが、2つの公式を知っているとギャンブルのように運に任せていてもπが求まります。まず、一つ目の公式です。円の面積の公式を覚えているでしょうか。
円の面積 = πr
rは半径ですね。まず、半径が1の円を描いてみましょう。サインとコサインで合成していますが、rpnの三角関数についてはrpn入門(初級編)の三角関数を参照ください。
+ + + ^+ + +
++ + | + + +
++ | ++
++ | ++
+ | +
+ | 半径1の円 +
+ | +
+ | +
+-------------------+------------------+
+ | +
+ | +
++ | ++
++ | ++
+ + | +
+ + | + ++
+ ++ + + + + +
ちょっと歪んでいますがテキスト文字での描画なので仕方ありません。頭の中で半径(r)が1のきれいな円を想像してください。
この円の面積を求めるとき公式によると、π×r×rになりますよね。そして、今回はrを1としましたから結局、面積はπということになります。つまり、面積を求めることはπを求めることに等しいということになりますよね。
次に円の右上四分の一だけを切り取ってみます。
^ + + +
| + +
| +
| +
| +
| +
| +
|
| +
| +
| +
|
| +
| +
|
+-------------------------------------+>
この面積はいくらかというと、面積(=π)の4分の1ですから、π÷4になります。
πを求めるためのもう一つ公式
さてπを求めるにはもう一つ公式が必要でした。それは距離の公式です。下図を見てください。
^y 1 b
|
|
|
| c
|
|
| a
|
|
|
|
|
| x
|o 1
+-------------------------------------->
x軸が1、y軸が1のエリアに点a、点b、点cがプロットされています。それぞれの座標軸は、点aが(0.5,0.5)、点bは(1,1)、点cは(0.3,0.7)です。それでは、原点oからそれぞれの点a、点b、点cへの距離はいくつでしょう。
原点o(0,0)からそれぞれの点との距離は、以下の公式を使って求められます。
距離=√(x + y )
式によると、x軸の値を2乗したものとy軸の値を2乗したものを足して、平方根を計算したものが原点oからの距離になります。では、実際に点aから点cまで計算してみましょう。
0.707107
>rpn 1 . * 1 . * + r
1.41421
>rpn .3 . * .7 . * + r
0.761577
原点oからの距離は、それぞれ点aが0.71、点bが1.41、点cが0.76くらいですね。
これでπを求める準備ができました。円の面積と原点からの距離の公式が分かっているとπが求まるのです。
πを求める考え方
上で描いた図の円の半径は1でした。すると原点oから円周までの距離はどこでも1になります。するとx軸が1、y軸が1の座標にある点a、点b、点cは、円の中にあるか外にあるか判定できることになります。
npdを使って、上の二つのグラフを合成すると以下のようなイメージです。点aと点cの原点oからの距離は1より小さいので円の中です。点bは1より大きいので円の外です。
| + +
| +
| +
| c +
| +
| +
| a
| +
| +
| +
|
| +
| +x
|o 1
+-------------------------------------+>
プログラム言語風に書くとこんな感じです。
円の中
}
else {
円の外
}
さあ、ここからがミソです。もっとたくさんの点をでたらめにプロットしてみて、円の中に入った数の割合が分かれば円の面積の四分の一と等しくなりませんか。
もともと、x軸1とy軸1に囲まれたエリアの面積は1×1=1です。すると単純に円の中に入った数を全部のプロット数で割れば面積が出てきますよね。
つまり、πとの関係で以下の式が成り立つことになります。
円の四分の一の面積 = ---- = 円の中に入った割合
4
ということは、
で、πが求まります。実際に点aから点cの3つのプロットの場合を計算してみましょう。円の中に入った割合は3点中、2点です。従って、
--- = ---- ⇒ π = 2÷3×4
3 4
ですよね。rpnで計算すると次のとおりです。
2.66667
結果は2.7です。πの3.14と全然、違いますね。でも待ってください。今回は適当に3点だけプロットしたものです。もっとたくさんプロットすれば答えは変わってくるはずです。
次は…
モンテカルロ法でπを求める原理は分かりましたね。どうやらプロット数が増えればπに近づいていくようです。そこで、rpnの乱数を使って大量にプロットしてみます。どこまでπに近似できるでしょうか。
rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。
xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。