Euler-Mascheroni 定数 - 努力で数論

本日は級数の分野よりハーモニックシリーズとオイラーの定数について、一つご紹介です。初めに参考文献としては

です。

オイラーの定数以下の $\lim_{n \to \infty}$ は有限な値に収束し、それはオイラーの定数と呼ばれる。 $$\gamma = \lim _{n \to \infty} \left( - \log n + \sum _{k=1}^{n} \frac {1}{k} \right)$$

定理の主張としては、ハーモニックシリーズは $\log n$ のオーダーで収束しその差はオイラーの定数分しかないということになっております。Wikipediaによると、このオイラーの定数は $0.57721566490153286060651209008240243104215933593992...$ であり、この記事ではこの極限の値を数値的に計算をする方法を見ていきたいというものです。誤差分を初めに $10^{-6}$ などとして置いて、極限値がこの誤差分分より小さくなる大きな $n$ を探していくということでもいいのですが（ $\epsilon - N$ 論法的なこと）、オイラーの定数は積分計算でも得られるので、まずはその準備としてディガンマ関数というものを用意しましょう。ここではガンマ関数はWeierstrassの無限積表示をされているものと仮定しましょう（https://okenta.hatenablog.jp/entry/2019/07/31/160242）。すると、

$$\psi(z) = \frac{d}{dz} \log (\Gamma(z)) = \frac{\Gamma'(z)}{\Gamma(z)} = \lim _{n \to \infty} (\log n - \sum_{k =0}^ {n} \frac{1}{z+k})$$ $$\psi(1) = \lim _{n \to \infty} (\log n - \sum_{k = 0}^ {n-1} \frac{1}{k+1} - \frac{1}{n+1}) = - \gamma$$

のようにオイラーの定数が得られます。他方、ここでガンマ関数の微分を積分表示の微分と捉えてみると

$$\Gamma'(z) = \frac{d}{dz} \int_{0}^{\infty} t^{z-1} e^{-t} dt = \int_{0}^{\infty} \frac{d}{dz} t^{z-1} e^{-t} dt = \int_{0}^{\infty} t^{z-1} \log t \cdot e^{-t} dt.$$

であり、 $\Gamma(1) = 1$ なのでこれらからオイラーの定数は次のように表されます。

$$\gamma = - \frac{\Gamma'(1)}{\Gamma(1)} = - \int_{0}^{\infty} \log t \cdot e^{-t} dt$$

変数変換すると

$$\gamma = - \int_{0}^{1} \log \left( \log \frac{1}{n} \right) dn$$

なので、この積分を数値計算してみましょう。数値計算はガウスの求積法により精度の良い計算が出来ます。他にも有名な定積分の数値計算としては trapezoidal rule や simpson' s rule、Euler-Maclaurinの公式などがあります。ガウスの求積法はscipyに内蔵される積分計算にも使われているのでかなり良い精度を保証してくれるようです。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def Legendre(N,x):
    P = np.zeros((N+1,len(x)))
    Q = np.zeros((N+1,len(x)))
    P[0], P[1] = 1 + 0*x, x
    Q[1] = 1 + 0*x
    for n in range(1,N):
        P[n+1] = ((2*n+1)*x*P[n] - n*P[n-1])/(n+1)
        Q[n+1] =    (2*n+1)*P[n] + Q[n-1] 
    return P, Q

def quadrature(N,imax=10):
    i = np.arange(N)
    x = np.cos(np.pi*(4*i+3)/(4*N+2))
    for i in range(imax):
        P, Q = Legendre(N,x)
        x -= P[N]/Q[N]
    w = 2/((1-x**2)*(Q[N]**2))
    return x, w 

def integrate(y,w):
    return np.sum(y*w)

def scale(a,b,x,w):
    t  = x*(b-a)/2+(b+a)/2
    dt = (b-a)*w/2
    return t, dt

def f(t):
    return np.log(np.log(1/t))

def DYI_I(a,b,N): 
    x,w = quadrature(N)
    t,dt = scale(a,b,x,w)
    return integrate(f(t),dt)

-DYI_I(0,1,10000)

するとDYI_I は次の値を返してくれます。

$$0.5772156588982373$$

これは文字通りDIYですのでargumentの3個目の値を大きなものにすると精度が高くなります。僕のパソコンでは時間がかかってそれは面倒なので、scipyにお願いしましょう。

import scipy.integrate as integrate
import scipy.special as special

def integrand(n):
    return np.log(np.log(1/n))

eulerconst = - integrate.quad(integrand, 0, 1)[0]
print("gamma = ",eulerconst)

するとこれは

$$gamma = 0.5772156649132887$$

を返します。こちらの方が初めのウィキの値と近くDYI_Iよりは精度がいいですね。本日はオイラーの定数に注目してガウスの求積法を実装してみましたがいかがだったでしょうか。上でfの部分を別の関数に置き換えても定積分の数値計算は簡単にできると思うのでそちらもやってみたい方はぜひ。ガウスの求積法についての参考文献も以下に挙げておきます。

ガウス求積 - Wikipedia