pythonでバーゼル問題を考える2

前回、$${\displaystyle\sum_{n=1}^{\infty} \dfrac{1}{n^{2}}=1+\frac{1}{2^2}+\frac{1}{3^2}+\frac{1}{4^2}+\frac{1}{5^2}+\frac{1}{6^2}\cdots=\frac{\pi^2}{6}}$$を言うために、
$${\displaystyle \sin x =\alpha({x+n\pi})\cdots({x+2\pi})({x+\pi})x({x-\pi})({x-2\pi})\cdots({x-n\pi})}$$

という多項式を想定し、$${\alpha}$$を適当に定めることで、

 $${\displaystyle \sin x =x\left(1-\frac{x^2}{\pi^2}\right)\left(1-\frac{x^2}{2^2\pi^2}\right)\cdots\left(1-\frac{x^2}{n^2\pi^2}\right)}$$

$${=x\prod_{n=1}^{\infty}\left(1-\frac{x^2}{n^2\pi^2}\right) \cdots(1)}$$

とあらわすことがきることを見てきました。ここからは、その続きです。

$\sin x$をTayler展開すると次のようになります。

$${\displaystyle \sin x = \sum _{n=0}^{\infty }{\frac {(-1)^{n}}{(2n+1)!}}x^{2n+1} =x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+\cdots(2)}$$

import matplotlib.pyplot as plt
import numpy as np
def sin_t(x, n):
   sin_approx = 0
   fact=1
   for i in range(1,n,2):
       fact*=i*(i-1) if i !=1 else 1
       sin_approx += (-1)**((i-1)/2) *(x**i/fact)
   return sin_approx

plt.style.use('seaborn-poster')
n = 50
x=np.arange(-6*np.pi,6*np.pi,0.1)
sin_v=np.vectorize(sin_t)

#1/nの合計
#グラフの作成
formula=[sin_v(x,n),np.sin(x)]
labels = [r'fx',r'$\sin x$']
colors=['b','r']
alphas=[1,0.2]
widths=[1,7]
plt.figure(figsize = (10,8))
for fx, label, color ,alpha, w in zip(formula, labels,colors,alphas,widths):
   plt.plot(x,fx, color, label = label ,alpha=alpha,linewidth=w)
plt.grid()
plt.title(r'Tayler展開でのsinの計算',fontname="MS Gothic")
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='lower center',ncol=2)
plt.xticks([-np.pi*6,-np.pi*5,-np.pi*4,-np.pi*3,-np.pi*2,-np.pi,0,np.pi,np.pi*2,np.pi*3,np.pi*4,np.pi*5,np.pi*6], 
          ['$-6\pi$','$-5\pi$','$-4\pi$','$-3\pi$','$-2\pi$','$-\pi$','0','$\pi$','$2\pi$','$3\pi$','$4\pi$','$5\pi$','$6\pi$'])
plt.show()

無限乗積展開とは異なり、$${5\pi,6\pi}$$と0から離れても正しく計算できることがわかります。
ここで(1)式と(2)式のxの係数を比較すると、両方とも1になるので(1)式の$\alpha$の設定は上手くいっているようです。そのうえで、$${x^3}$$の係数を比較します。(1)式は一番初めに$${x}$$があるので、掛け算して残り$x^2$となるのは各括弧の中身で1つだけ$${\displaystyle\frac{1}{n^2 \pi^2}}$$で、残りは1を掛けた項の和になります。一方、(2)の$${x^3}$$の係数は$${3!=6}$$となります。これらの係数は等しいと考えられるので、ここだけに注目すると次の関係になります。

$${\displaystyle -\frac{1}{\pi^2}-\frac{1}{2^2\pi^2}-\frac{1}{3^2\pi^2}-\frac{1}{4^2\pi^2}=-\frac{1}{3!}=-\frac{1}{6}}$$
両辺に$${-\pi^2}$$を掛けると、

$${\displaystyle\sum_{k=1}^{\infty} \dfrac{1}{k^{2}}=1+\frac{1}{2^2}+\frac{1}{3^2}+\frac{1}{4^2}+\frac{1}{5^2}+\frac{1}{6^2}\cdots=\frac{\pi^2}{6}}$$

実際に計算すると、大体等しくなるようです。

import numpy as np
sigma=0
for i in range(1,1000):
   sigma+=(1/i**2)
print(f'{sigma},{np.pi**2/6}')

#1.6439335666815615,1.644934066848226

それではこの方法で、どのくらいの精度で円周率を計算することができるか試してみます。
まず、上の式は
$${\displaystyle\pi=\left(\sum_{n=1}^{\infty}\dfrac{1}{n^{2}}\times6\right)^{1/2}}$$変換することができるので、10万回までループさせてみます。

d=5
sigma_2=0
kubun = 10
for i in range(1,10**d+1):
   sigma_2 +=1/(i** 2)
   if i == kubun:
       print(f'{i:>6}回,pi={(sigma_2*6)**(1/2)}')
       kubun*=10
       
#output       
   10回,pi=3.04936163598207
  100回,pi=3.1320765318091053
 1000回,pi=3.1406380562059946
10000回,pi=3.1414971639472147
100000回,pi=3.141583104326456

円周率は3.141592653589なので、1000回ループして3.14で小数点2桁、100000回ループして3.1415小数点4桁までということで、話としては面白くもありますが、効率は余りよくはないようです。