乱数を生成して円周率を求めてみよう！っていうのを過去にどこかでやった記憶があったので、思い出しながらPythonでやってみました。
乱数を使って近似値を求めることをモンテカルロ法っていうらしいです。
モンテカルロ法による円周率の計算
モンテカルロ法とは...
シミュレーションや数値計算を乱数を用いて行う手法の総称
モンテカルロ法 - Wikipedia より引用
こちら↓のサイトにてモンテカルロシミュレーションで円周率を求める方法について詳しく説明されています。
今回プログラムを組むのに使った方法を紹介しておきます。
-1≦ｘ≦1
-1≦ y ≦1
となるように乱数x,yを浮動小数で生成します。
点(x,y)は、下の画像の赤または黄色で塗られた範囲に落ちます。
x^2+y^2≦1 の点は黄色で塗られた範囲
x^2+y^2>1 の点は赤色で塗られた範囲
c個の異なる点(x,y)を用意し、そのうち黄色の範囲に落ちた点を数え、その数をn個とします。
n/c が、c個用意した点のうち黄色の範囲に落ちた割合となります。
(黄色の範囲の面積)＋(赤色の範囲の面積)=4
であることから(一辺2の正方形の面積)、
cの数が大きいほど
4×(n/c) が黄色の範囲の面積に近づきます。...(1)
円の面積が πr^2 (π:円周率、r:半径) と表わされ、
r=1であるので、黄色の範囲の面積はπです。
以上より
(1)で求まった値が、円周率πの近似値となります。
c=10000000 とし、この方法で円周率をpythonに計算させてみましょう。
私の環境では、この計算に10秒ちょっとかかりました。
cの大きさはお使いの環境に合わせて適当に設定してください。
乱数生成のために、randomモジュールをimportして使います。
.pyファイルの中身
import random
n=0
c=10000000
for z in range(c):
x=random.uniform(-1,1)
y=random.uniform(-1,1)
if x**2+y**2<=1: #円の内部に落ちたものをカウント
n=n+1
print(4*n/c)
出力例
3.1419892
最後に
小数第三位まではあってますね。
(実際の円周率の値は、π=3.141592...)
あくまで近似計算なので、ピッタリとはいきません。
誤差の検討などは冒頭に紹介した「高校数学の美しい物語」にて数学的に説明されているので、気になる人はどうぞ。
それではまた。
