pythonで円周率を求める
pythonで円周率を求めるアルゴリズム
○ライプニッツの方法
収束が遅い。
100万回計算して、6桁くらい。
pi4 = 0 for i in range(1000000): pi4 += (1 / (i * 4 + 1) - 1 / (i * 4 + 3)) print(pi4 * 4) # 3.141592153589902
○区分求積を使う。
収束が速い。
25回の計算で15桁くらいまで出ます。
pi4 = 0 for i in range(1, 50, 2): pi4 += (-1)**(i//2) * (1/2**i + 1/3**i) / i print(pi4 * 4) # 3.1415926535897922
○ガウス・ルジャンドルの方法
収束が非常に速く、4回も計算すれば15桁まで出ます。
a, b, t, p = 1, 1/2**0.5, 1/4, 1 for i in range(4): a, b, t, p = (a+b)/2, (a*b)**0.5, t-p*((b-a)/2)**2, p*2 print((a + b)**2 / (4 * t)) # 3.1405792505221686 # 3.141592646213543 # 3.141592653589794 # 3.141592653589794
○モンテカルロ法
乱数を使って円周率を求める方法。
原点中心、半径1の円とその円に外接する正方形の面積比がπ/4であることを利用し、1万個の座標を乱数で生成して円周率を求めてみた。
サンプルコード(1) randomモジュールで。
import random inner = 0 for i in range(10000): x, y = random.random(), random.random() if x**2 + y**2 < 1: inner += 1 print(inner *4 / 10000) # 3.1446
サンプルコード(2) NumPyのrandomで。
np.random.seed(1) # 乱数の初期化。好みの数で。 arr = np.random.rand(2, 10000) # 2行10000列の乱数ndarray distance = (arr**2).sum(axis=0) # 全要素を2乗し縦に合計すると原点からの距離 count = (distance < 1).sum() # 距離が1未満か判断しTrueの個数を調べる print(4 * count / 10000) # 3.1348
1万個の散布図を描くとキャンバスが埋まってしまうので、1000個で描画
import numpy as np import matplotlib.pyplot as plt c1 = plt.Circle((0, 0), radius=1, fc="None", ec="r", linewidth=3) ax = plt.gca() ax.add_patch(c1) plt.axis("scaled") plt.xlim(0, 1) plt.ylim(0, 1) x = np.random.rand(1000) y = np.random.rand(1000) plt.scatter(x, y, marker="x") plt.show()
○逆三角関数を使う方法
cos(π)=-1なので
π=arccos(-1)
import math print(math.acos(-1)) # 3.141592653589793