円周率を近似計算してみた with Python
今まで全く気にしてこなかったのですが、円周率ってどうやって導き出すんだろうと思ったので、Pythonでプログラムしてみました。
色々な求め方があるんですね。。。
その中でなんとなくPythonで組めそうな感じがした以下の２つを試してみました。
グレゴリー級数
- 【参考】Python で円周率を求める
end = 1000
phi = 0
isTurnOver = True
for i in range(1, end, 2):
if isTurnOver:
g = 4 * 1 / i
else:
g = -1 * 4 * 1 / i
phi += g
isTurnOver = not isTurnOver
# 最後の方だけ「パイ」を表示してみる
if i > end*0.99:
print('π=', i, phi)
π= 991 3.1395765266063047
π= 993 3.1436047239879765
π= 995 3.139584623485464
π= 997 3.143596659593789
π= 999 3.139592655589785
結果は・・・π= 3.139592655589785
モンテカルロ法
- 【参考】モンテカルロ法で円周率を求る
bR = .5 # 半径
sS = 0 # 四角のプロット数
cS = 0 # 円のプロット数
count = 500
np.random.seed(0) # 乱数を固定
ps = np.random.uniform(-1*bR, bR, (count, 2))
x = ps[:, 0]
y = ps[:, 1]
r0 = np.sqrt(x**2 + y**2) # ランダムにプロットした点
for r in r0:
if bR > r:
cS += 1; # 円の中
sS += 1; # 四角形の中
print('π=', 4. * cS / sS)
結果は・・・π= 3.064
乱数を固定せずに何度かやれば、下3, 4桁位までは精度が出ることがあったので、誤差のせいか今回の500というプロット数は圧倒的に少ないということでしょうか・・・
意外と精度が出なかった
グレゴリー級数で1/1000まで、モンテカルロ法で500点で実験してみましたが意外と3.141592…とはならないんだなって思いました。バグがなければですが・・・
割り算した後足したり引いたりの所で、ごく僅かですが桁が落ちてるのが精度が出ない原因なのでしょうか・・・？
また時間があったら調べてみようかと思います。
実は
ファインマンポイントという円周率の下762桁目から「9」が6回続くらしいのですが、この「999999」が出るまでやってみたかったのですが、グレゴリー級数の方を10億までやってみたのですがまだ駄目でした。いつかチャレンジしてみます。
