Pythonで円周率の計算を可視化する求め方[モンテカルロ法]

ちょっと前書いた手の絵をお絵かきしたら割と上手にできました、おはようございます、まがりかどです。

そんなことより今日は円周率をPythonで求めてそれを可視化していきましょうの回でございます。

円周率が何かはまあ大丈夫なはずなので、どうやって求めていくかまず書いていきますじゃ。

モンテカルロ法

知っている人多いのか少ないのかなんとも微妙な立ち位置の単語だなあと思いますけど今回はそんなモンテカルロ法で円周率を求めていきたいと思います。

モンテカルロ法というのはですね簡単に言うと乱数を使って数値計算していきましょうというやつですね。

もちろん乱数を使うのでモンテカルロ法では近似値しか求まりません。
というかどの方法でも円周率は近似値しか求められないですけど。

具体的には今回はこうします

ここに二次元の座標があります。

そして一片の長さが2の正方形とその中に半径が1の円がありますね。

正方形の中に点を任意の個数（これをloopと置く）打っていく　多いほうが精度が高い
そのうち円の中に入った点の数（これをcnt を置く）を数える
$\frac{cnt}{loop}\times4$で円周率のできあがり

なぜでしょう。

めちゃくちゃ点いっぱい打ったら大体まんべんなくなりますよね。

なので点を打てば打つほど正方形と円の面積の比が、loop とcntが近づいていきますよね。

なので正方形の面積が4、円の面積が$\pi$より

無限回点を打った時、

$$\frac{cnt}{loop} = \frac{\pi}{4}$$

\[\pi = \frac{cnt}{loop}\times4\]

となります。

方針が決まったということでコードを書くぞ

モンテカルロ法で円周率のコード

import random
import numpy as np
import matplotlib.pyplot as plt
loop = int(input("How many times?"))
cnt = 0#円の中に入った点の数
ran = 100000
fig = plt.figure()
ax = fig.add_subplot()
ax.set_aspect('equal')#表示するやつを正方形にする
plt.xlabel("x")
plt.ylabel("y")
plt.xlim([-ran, ran])
plt.ylim([-ran, ran])
for i in range(loop):
    x = random.randint(-ran, ran)
    y = random.randint(-ran, ran)
    dis = x * x + y * y
    if dis <= ran * ran:
        cnt += 1
        plt.scatter(x, y, c = 'b', s = 10)
    else: plt.scatter(x, y, c = 'r', s = 10)
ans = (cnt / loop) * 4
print(ans)
plt.show()

はい、こんな感じです。

ran = 100000
fig = plt.figure()
ax = fig.add_subplot()
ax.set_aspect('equal')#表示するやつを正方形にする
plt.xlabel("x")
plt.ylabel("y")
plt.xlim([-ran, ran])
plt.ylim([-ran, ran])

この辺でグラフの設定をしていますね。

上の図では-1から１までの範囲でしたがそうすると少数が出てきて嫌なので-100000から100000にしておきました。

これの二行目から四行目は縦横のアスペクト比を1:1にするために書いています。

そしてつぎ、

for i in range(loop):
    x = random.randint(-ran, ran)
    y = random.randint(-ran, ran)
    dis = x * x + y * y
    if dis <= ran * ran:
        cnt += 1
        plt.scatter(x, y, c = 'b', s = 10)
    else: plt.scatter(x, y, c = 'r', s = 10)

上で何回ループするかをloopにインプットしているのでその回数ループします。

まず打つ点の座標を-100000から100000の間でx,yそれぞれ適当に決めます。

その点が円の中に入っているかどうかは原点からの距離が円の半径以下かどうかで判定します。

半径以下なら入っています。

式にすると三平方の定理より、

$$x^2 + y^2 <= 100000^2$$

のとき円の中にあります。

普通に平方根で距離求めて100000以下か判定してもいいですが少数が出てきて嫌なのとちょっと遅くなりそうなのでやめときます。

円の中に入っていたら青で、そうじゃなかったら赤で点を打ちます。

点のサイズは10くらいがちょうどよかったです。

ということで実行すると、