モンテカルロでπを計算 | 数学アラカルト [実践] | 逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓rpnの実践ウェブ   
rpn hacks! アールピーエヌ・ハックスサイトマップ

rpn | 実戦 | 数学アラカルト | 数学に関する話題をrpnで探求!数字と遊んでみよう。

HOME > 実践 > 数学アラカルト > モンテカルロでπを計算

hatena twitter facebook rss ソーシャルブックマーク

モンテカルロでπを計算

 ある点から一定の距離にある平面上の点の集合を円といいます。この円の直径と円周の割合が円周率πなのですが、3.1415926535...ずっと続く不思議な数字です。このπの値を求めようと古くはギリシャのアルキメデスや、ドイツのルドルフ、日本の関孝和などが挑戦してきました。

一例を挙げると、アルキメデスは直径が1の円の内側にピッタリはまる正6角形(内接正6角形)の外周の長さと、円の外側にピッタリはまる正6角形(外接正6角形)の外周の長さを求めて、πはその間にあるとしました。

  内接正6角形の外周 < π < 外接正6角形の外周


そして次は正12角形、次は正24角形…。そして、最後は正96角形までにして計算していきました。その最終的な式が以下になります。

     10              10
  3 ---- < π < 3 ----
     71              70


πは左辺と右辺の間にあるというわけです。rpnで計算してみます。

  >rpn 3 10 71 / + 3 10 70 / +
  3.14085 3.14286

情報rpnの簡単な説明はrpn基礎に、詳しい説明はrpn入門にありますので、ご参照ください。

なんと、この方法でアルキメデスは3.14まで正確に求めていました。その後、これと同様な方法で16世紀のヨーロッパにおいては、πの近似値として次の分数に至っていました。

         355
  π ≒ -----
         113


この分数をrpnで計算すると…

  >rpn 355 113 /
  3.14159


もはや実用には何の問題もない精度ですね。

πを求めるための2つの公式

 さて、ここではアルキメデスが使った幾何学的な方法ではなくて、ちょっと変わった方法でπの計算をしてみます。モンテカルロ法というやり方です。これはでたらめな数(乱数)を使って数値計算するときに使う呼び名です。カジノの都モナコのモンテカルロに因んで付けられた名前です。

でたらめな数、乱数を使うことでπが求まるのかですが、2つの公式を知っているとギャンブルのように運に任せていてもπが求まります。まず、一つ目の公式です。円の面積の公式を覚えているでしょうか。

                 2
  円の面積 = πr


rは半径ですね。まず、半径が1の円を描いてみましょう。サインとコサインで合成していますが、rpnの三角関数についてはrpn入門(初級編)三角関数を参照ください。

  >rpn @t .1 + #t @t s @t c -r 62 | xyp -d+ | npd
                + + + ^+ + +
           ++ +       |      + + +
        ++            |           ++
     ++               |              ++
    +                 |                +
   +                  |   半径1の円     +
  +                   |                  +
  +                   |                  +
  +-------------------+------------------+
  +                   |                 +
   +                  |                 +
    ++                |               ++
      ++              |             ++
        + +           |            +
           + +        |       + ++
               + ++ + + + + +


ちょっと歪んでいますがテキスト文字での描画なので仕方ありません。頭の中で半径(r)が1のきれいな円を想像してください。

この円の面積を求めるとき公式によると、π×r×rになりますよね。そして、今回はrを1としましたから結局、面積はπということになります。つまり、面積を求めることはπを求めることに等しいということになりますよね。

次に円の右上四分の一だけを切り取ってみます。

  >rpn @t .1 + #t @t s @t c -r 62 | xyp -x,1 -y,1 -d+
  ^  +   +   +
  |              +  +
  |                     +
  |                        +
  |                          +
  |                             +
  |                               +
  |
  |                                 +
  |                                   +
  |                                    +
  |
  |                                     +
  |                                     +
  |
  +-------------------------------------+>


この面積はいくらかというと、面積(=π)の4分の1ですから、π÷4になります。

πを求めるためのもう一つ公式

 さてπを求めるにはもう一つ公式が必要でした。それは距離の公式です。下図を見てください。

  >rpn .5 .5 1 1 .3 .7 | xyp -x,1 -y,1 -k3 -da -m
  ^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)からそれぞれの点との距離は、以下の公式を使って求められます。

            2    2
  距離=√(x  + y )


式によると、x軸の値を2乗したものとy軸の値を2乗したものを足して、平方根を計算したものが原点oからの距離になります。では、実際に点aから点cまで計算してみましょう。

  >rpn .5 . * .5 . * + r
  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より大きいので円の外です。

  ^y +   +   +                           b
  |              +  +
  |                     +
  |                        +
  |          c               +
  |                             +
  |                               +
  |                  a
  |                                 +
  |                                   +
  |                                    +
  |
  |                                     +
  |                                     +x
  |o                                     1
  +-------------------------------------+>


プログラム言語風に書くとこんな感じです。

  if (√(x^2+y^2) <= 1) {
     円の中
  }
  else {
     円の外
  }


さあ、ここからがミソです。もっとたくさんの点をでたらめにプロットしてみて、円の中に入った数の割合が分かれば円の面積の四分の一と等しくなりませんか。

もともと、x軸1とy軸1に囲まれたエリアの面積は1×1=1です。すると単純に円の中に入った数を全部のプロット数で割れば面積が出てきますよね。

つまり、πとの関係で以下の式が成り立つことになります。

                         π
  円の四分の一の面積 = ---- = 円の中に入った割合
                         4


ということは、

  π = 円の中に入った割合 × 4


で、πが求まります。実際に点aから点cの3つのプロットの場合を計算してみましょう。円の中に入った割合は3点中、2点です。従って、

   2      π
  --- = ----    ⇒   π = 2÷3×4
   3      4


ですよね。rpnで計算すると次のとおりです。

  >rpn 2 3 / 4 *
  2.66667


結果は2.7です。πの3.14と全然、違いますね。でも待ってください。今回は適当に3点だけプロットしたものです。もっとたくさんプロットすれば答えは変わってくるはずです。

次は…

 モンテカルロ法でπを求める原理は分かりましたね。どうやらプロット数が増えればπに近づいていくようです。そこで、rpnの乱数を使って大量にプロットしてみます。どこまでπに近似できるでしょうか。

続き(part2)はこちらLinkIcon

警告rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。

情報xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。

数学アラカルト

数字と遊んでみよう

※実践コーナーのTOP

紹介 rpnの利用シーンはこちら…

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

part2

数字と遊んでみよう

※実践コーナーのTOP

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

書籍紹介

記事に関連した書籍

本ウェブサイトで扱った話題に関連した書物で、スタッフが実際に読了したものを紹介。

書籍数学の書籍
数の世界は思ったよりもエキサイティング

書籍投資の書籍
失敗しない投資には広範囲で実践的な知識が必要

警告バックスラッシュはエンマークに読み替えて下さい(IE)。
バックスラッシュとエンマーク

警告文字で作られた図表や式が崩れることがあります。ブラウザによっては固定幅フォントをMSゴシックにするときれいに表示されます。それでも崩れる場合は図表や式をメモ帳にコピー後、閲覧下さい。

警告rpn試用版と標準版(2kリビジョン)はダブルクォートで囲って下さい。

rpn 1 2 + ⇒ rpn "1 2 +"
rpn 1 -c foo ⇒ rpn "1" -c "foo"

ダブルクォートは省略できることが多いのですが、慣れない間は囲んだほうが無難です。なお、本ウェブサイトの記事ではrpn標準版(98リビジョン)を使用しているため囲っていません。詳しくは技術サポートの「rpn TIPS参照ください。

注意rpnの障害情報はこちら

警告rpn試用版の場合、複雑なプログラムや処理時間のかかるプログラムの一部には動作しないものがあるかもしれません。あくまで無料提供であることを勘案・了承ください。rpn標準版は、すべてのプログラムが動作します。