はじめに
ガンマ関数は、$\Gamma[p]$と表されるが、$p>0$のときは、広義積分で定義される。一方で、$p<0$の場合は、$p>0$の場合の$\Gamma[p]=\frac{\Gamma[p+1]}{p}$を用いて帰納的に定義される。
そこで、今回は長方形分割による積分値の算出を数値解析を行うことで、ガンマ関数のグラフを描写することができるかを試みる。
以下結果を示す。
1.散布図で表した場合$(-4<p<4)$のとき
上手く無限大を表現することができなかったので、あまりいいグラフは描写することができなかった。
2.散布図で表した場合$(-2<p<2)$のとき
それっぽいグラフは描写できているが、境界付近の値の誤差が大きすぎる。
以下、サンプル点の分布を等間隔ではなく、対数間隔にした場合のグラフを複数示す。
このように、局部だけで見れば、かなり高い精度でガンマ関数を再現できていることが分かる。
ガンマ関数について
p>0のとき
ガンマ関数は、広義積分を用いて以下のように定義される。
\Gamma[p]=\int_{0}^{\infty}x^{p-1}e^{-x}dx
p<0のとき
ガンマ関数の性質として、以下を認めることにする。
\Gamma[p+1]=p\Gamma[p]
この場合、
\Gamma[p]=\frac{\Gamma[p+1]}{p}
となる。したがって、$p>0$のときのガンマ関数の値を上式に帰納的に代入していくと、$p<0$のときのガンマ関数の値を求めることができる。
区分求積法について
区間、$a<x<b$における以下の積分は、無限に分割される短冊状の長方形面積の和であらわすことができる。
\int_a^b f(x)dx = \lim_{n\to \infty} \frac{1}{n}\sum^{n} _{k=1}f(a+\frac{b-a}{n}k)
ただし、数値計算では、無限大を上手く扱えないので、有限の大きな数値で代用する。
プログラム
上記の考えを用いて以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
num=100
L=4
p_ary=np.linspace(-L,L,num)
#p_ary=np.logspace(np.log10(1e-6),np.log10(1),num)-1
S_ary=[]
#区分求積法による数値積分
##積分範囲
a=0
#b→∞が理想だが、無限大を表現できないので、大きな値を代入
b=10000
##分割数
n=10000
# 短冊の幅Δx
dx = (b-a)/n
def gamma_s(p):
# 面積の総和
if p > 0:
s = 0
#print(p)
for i in range(1,n):
# x,yの値
#print(s)
x1 = a + dx*i
#print(x1)
f1 = x1**(p-1)*np.exp(-x1)
# 面積
s += dx*f1
return s
#print(s)
else:
return gamma_s(p+1)/p
# 積分する関数の定義
# def f(x):
# return x**(p-1)*np.exp(-x)
for i in range(num):
S_ary.append(gamma_s(p_ary[i]))
#plt.plot(p_ary,S_ary)
plt.scatter(p_ary,S_ary)
plt.xlabel("p")
plt.ylabel("Γ(p)")
plt.savefig("gamma_sum.png")
plt.show()
これを実行する。（定義域やプロットの分布を変えることで、冒頭に示したようなグラフを得ることができる。）
まとめ
今回は、ガンマ関数の値を数値積分で算出することができるかを調査した。具体的には、Pythonを用いて、正の引数で定義されるガンマ関数の値をまず計算した。次に、帰納的に推定される負の引数のガンマ関数の値も求めることで、ガンマ関数のグラフ化を試みた。しかし、上手く境界付近や無限小,無限大を扱うことができなかったため、誤差が大きくなってしまった。