はじめに
円周率の近似値を求める方法は、たくさんある。そこで、できるだけマニアックな方法で円周率の近似値を出すことができたら喜ばしい限りである。そこで、前半戦では、今回はガウス関数の面積は円周率で表されることを利用して、モンテカルロ法で面積を近似し、円周率を推定する。後半戦では、三角関数のsin関数をマクローリン展開で多項式近似し、二分法を用いることで円周率を推定する。
ガウス関数による推定
ガウス関数の広義積分は以下のように表される。
\int_{-\infty}^{\infty}e^{-x^2} = \sqrt{\pi}
したがって、上記のガウス関数$y=f(x)=e^{-x^2}$と$y=0$で囲まれる面積をモンテカルロ法で求める。ここで、$0<f(x)<1$であることから、$-m<x<m,0<y<1$という範囲内における$y<f(x)$の領域の面積を求めればよい。ただし$L$は十分大きな正の実数と仮定する。
そこで、以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
import math
import random
n=100
m=30
N=3
num=10**N
x= np.linspace(-m, m, 100)
y=np.exp(-x**2)
x_S=[]
y_S=[]
x_S_not=[]
y_S_not=[]
for i in range(num):
x_seed=random.random()*2*m-m
y_base=np.exp(-x_seed**2)
y_seed=random.random()
if y_seed<y_base:
x_S.append(x_seed)
y_S.append(y_seed)
else:
x_S_not.append(x_seed)
y_S_not.append(y_seed)
S_gauss=(len(x_S)/num )*1*2*m
#print(S_gauss**2)
pai_gauss=S_gauss**2
plt.title("ガウス分布と円周率の近似値%f"% (pai_gauss))
plt.plot(x,y)
plt.scatter(x_S,y_S,color="red")
plt.scatter(x_S_not,y_S_not,color="blue")
plt.savefig("gaussian_S.png")
plt.show()
これを実行すると以下のようになる。領域内に入った、サンプル点を赤で、外れてしまったものを青で示してある。
もちろん、サンプル点の数を十分に大きく($1.0\times 10^7$個程度)にすれば高い近似精度を示すが、計算時間が比例的に大きくなる。
二分法による推定
そこで、次は三角関数を用いた推定について考える。
g(x)=sin(x)-0.5=0
という方程式を考える。ただし、初期解は$0.5<x<0.9$とし二分法を用いて解くものとする。
この方程式の解は上記の範囲ならば$x=\frac{\pi}{6}$であるは明らかである。
そこで、$sin(x)$はマクローリン展開を用いて展開すると以下のようになる。
sin(x)\approx x-\frac{x^3}{3!}+\frac{x^5}{5!}+\cdot \cdot \cdot +(-1)^n \frac{x^{2n+1}}{(2n+1)!}
ただし、$x$は小さく$n$については十分に大きい値とする。
(近似の目安について厳密な数値を用いた説明はここでは省く)
そこで、まずマクローリン展開の精度を確かめるために以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
a=3
b=4
n=1
p=10
num_ary=[]
val_ary=[]
x_ary=[0]
def sin_near(x):
num_ary=[]
val_ary=[]
# # valは(2n+1)の階乗である
for m in range(2*n):
val=1
num=2*m+1
for i in range(1,num+1):
val=val*i
num_ary.append(num)
val_ary.append(val)
B=0
for i in range(len(num_ary)):
A=(-1)**i*(x**num_ary[i]/val_ary[i])
B=B+A
return B
theta=np.linspace(0,math.pi/2,100)
y_real=np.sin(theta)
y=sin_near(theta)
plt.plot(theta,y_real,label="sin(theta)")
plt.plot(theta,y,label="sin_near(theta)")
plt.xlabel("theta[rad]")
plt.ylabel("sin(theta)と近似値")
plt.legend()
plt.title("sin(theta)の近似値と真値の比較")
plt.savefig("sin(theta)の近似値と真値の比較.png")
plt.show()
このプログラムを実行すると以下のようになる。
このように、$0<x<\frac{\pi}{6}$程度では、$sin(x)$の近似式による計算と真値はほぼ等しくなることが分かる。
したがって、これをもとにして以下のような方程式を考える。
g(x)=sin(x)-0.5=0
ここで$y=g(x)$は、$x=\frac{\pi}{6}$付近において単調な関数であり、$y=0$と交わることから二分法を用いることができる!!
ゆえに、以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
a=0.5
b=0.9
n=3
i_ary=[]
err=1.0*10**-5
err_ary=[]
err_pi_ary=[]
def sin_near(x):
num_ary=[]
val_ary=[]
# # valは(2n+1)の階乗である
for m in range(2*n):
val=1
num=2*m+1
for i in range(1,num+1):
val=val*i
num_ary.append(num)
val_ary.append(val)
B=0
for i in range(len(num_ary)):
A=(-1)**i*(x**num_ary[i]/val_ary[i])
B=B+A
return B
i=0
while True:
i=i+1
# 二分法の解の範囲指定
c=(a+b)/2
if (sin_near(a)-0.5)*(sin_near(c)-0.5)<0:
b=c
else:
a=c
# 終了条件
if abs(sin_near(c)-0.5)<err:
break
i_ary.append(i)
err_ary.append(abs(sin_near(c)-0.5))
err_pi =abs(c-math.pi/6)/(math.pi/6)*100
err_pi_ary.append(err_pi)
plt.title("円周率の近似値:{}".format(c*6))
plt.plot(i_ary,err_pi_ary)
plt.xlabel("試行回数")
plt.ylabel("誤差率[%]")
plt.savefig("二分法とマクローリン展開による円周率の近似.png")
plt.show()
これを実行すると以下のようなグラフが出力される。
このように、かなり良い精度で円周率の近似値を推定できていることが分かる。二分法は単調増加かつ求めたい領域において解が一つであるときにしか求めることができない。しかし、今回のようにそれらを満たせば単純なアルゴリズムでも高い精度で解を推定できることが分かった。
まとめ
今回は、2種類の方法をもとにして円周率の近似値を推定することを試みた。結果、三角関数をマクローリン展開で近似し、二分法で方程式を解くアルゴリズムを用いれば、比較的計算量がいらずに高い精度(小数点下3桁レベル)で円周率の近似値が推定できた。
参考文献
ガウス関数の広義積分
三角関数のマクローリン展開