円周率を近似計算してみた with Python

今まで全く気にしてこなかったのですが、円周率ってどうやって導き出すんだろうと思ったので、Pythonでプログラムしてみました。

色々な求め方があるんですね。。。
その中でなんとなくPythonで組めそうな感じがした以下の2つを試してみました。

グレゴリー級数

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億までやってみたのですがまだ駄目でした。いつかチャレンジしてみます。

こういうブログもありました。

人気のあるブログ:

コメントを残す

メールアドレスが公開されることはありません。

CAPTCHA


このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください