pythonでバーゼル問題を考える
次のような数列の合計を考えます。
$${\displaystyle\sum_{k=1}^{n} \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{1}{k^2}}$$
これも、簡単化と思いきや、とんでもなく奥が深いことがわかりました。まずは積分の知識を使い、この数列の合計の範囲を調べます。
$${y=\frac{1}{k^2}\quad(x\gt 0 )}$$は単調減少なので
$${\displaystyle\int_{k}^{k+1}\frac{1}{x^2}dx \lt \frac{1}{k^2} \lt \int_{k-1}^{k}\frac{1}{x^2}dx}$$
から、
$${\displaystyle 1-\dfrac{1}{n}+\dfrac{1}{n^{2}} \lt \sum_{k=1}^{n} \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{1}{k^2} \lt 2-\frac{1}{n}}$$
となります。それでは早速これらの関係をグラフにしてみます。
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-poster')
n = 100
x = np.arange(1, n+1)
#1/nの合計
sigma=0
func=[]
for i in x:
sigma+=(1/i**2)
func.append(sigma)
func_np=np.array(func)
#グラフの作成
formula=[2-1/x,func_np,1-1/x+1/(x**2)]
labels = [r'$2-\frac{1}{n}$',r'$\sum \frac{1}{x^2}$',r'$1-\frac{1}{n}+\frac{1}{n^2}$']
colors=['k','r','y','b']
plt.figure(figsize = (10,8))
for fx, label, color in zip(formula, labels,colors):
plt.plot(x,fx, color, label = label)
plt.grid()
plt.title(r'$\sum \frac{1}{x^2}$とnの関係',fontname="MS Gothic")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
極めて大雑把ですが、$${\displaystyle 2-\dfrac{1}{n}}$$は2に収束するので、これより小さい$${\displaystyle \sum_{k=1}^{n} \dfrac{1}{k^{2}}}$$も収束します。そしてグラフを見ると1.6より少し大きな値になります。
そして、この1.6より少し大きな値はいくつになるかと調べてみると、驚くべきことに、nを無限大にすると
$${\displaystyle \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}}$$
という結果になります。実際に計算してみます。
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.6449340668482264n=1000で計算すると小数点第2位まで正しく計算することができます。正の整数2乗して1を割った値を順次足しこんでいくと円周率が現れるというのは驚きです。なぜこのようなことが起このか探ってみたいと思います。まず、唐突ですが三角関数のsinカーブを考えます。
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-poster')
n = 150
theta=np.arange(-3*np.pi,3*np.pi,0.1)
formula=[np.sin(theta)]
labels = [r'$sin(x)$']
colors=['b']
plt.figure(figsize = (10,8))
for fx, label, color in zip(formula, labels,colors):
plt.plot(theta,fx, color, label = label)
plt.grid()
plt.title(r'sinカーブ',fontname="MS Gothic")
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='lower right')
plt.xticks([-np.pi*3,-np.pi*2,-np.pi,0,np.pi,np.pi*2,np.pi*3], ['$-3\pi$','$-2\pi$','$-\pi$','0','$\pi$','$2\pi$','$3\pi$'])
plt.show()
見てのとおり、$${x=\cdots-3\pi,-2\pi,-\pi,0,\pi,2\pi,3\pi\cdots}$$のとき$${\sin x}$$は0になります。ところで、例えば$${y=4x^2-4}$$のような関数は、y=0にするためにはx=1,-1とする必要があるので因数分解すると$${4(x+1)(x-1)}$$のようにと表現することができます。ということは、$${\sin x}$$を次の通り書き換えることを考えます。
$${\displaystyle \sin x =\alpha({x+n\pi})\cdots({x+2\pi})({x+\pi})x({x-\pi})({x-2\pi})\cdots({x-n\pi})}$$のような多項式と置きます。ここで少し雑ですがある理由で$${\alpha=1}$$と考えることができ、この式は、$${\displaystyle \sin x =\left({1+\frac{x}{n\pi}}\right)\cdots\left({1+\frac{x}{2\pi}}\right)\left({1+\frac{x}{\pi}}\right)x \left({1-\frac{x}{\pi}}\right) \left({1-\frac{x}{2\pi}}\right) \cdots\left({1-\frac{x}{n\pi}}\right)}$$となります。ところで$${\displaystyle \left({1+\frac{x}{n\pi}}\right)\left({1-\frac{x}{n\pi}}\right)=1-\frac{x^2}{n^2\pi^2}}$$なので、上記の式はさらに次の通り書き換えることができます。
$${\displaystyle 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)}$$
$${\displaystyle=x\prod_{n=1}^{\infty}\left(1-\frac{x^2}{n^2\pi^2}\right)}$$
このような展開をワイエルシュトラスの因数分解定理（無限乗積展開）と呼ばれています。それでは、早速グラフにしてみます。
import matplotlib.pyplot as plt
import numpy as np
def sin_w(x,n):
mul = x
for r in range(1,n):
mul*=1-x**2/((r**2)*(np.pi**2))
return mul
plt.style.use('seaborn-poster')
n = 150
x=np.arange(-6*np.pi,6*np.pi,0.1)
sin_v=np.vectorize(sin_w)
#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'ワイエルシュトラスの因数分解定理',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()
xの値が0から遠くなると、yの絶対値が1の付近で差が発生しますが、xが$${-2\pi}$$から$${2\pi}$$の間ではキレイに一致することがわかります。それでは、なぜここから
$${\displaystyle \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}}$$が言えるのかは、また次回。