Python（SymPy）で学ぶガンマ関数とベータ関数

数学の具体的な計算にPython（SymPy）を使って、数学もPython（SymPy）も同時に学んでしまいましょう。今回はPython（SymPy）を使ってガンマ関数の性質をいくつか調べます。ガンマ関数の仲間であるベータ関数についても簡単に紹介します。

以下全体を通して

import sympy as sp

としたうえで、SymPyの組み込み関数の呼び出しを sp と省略している点を注意しておきます。

ガンマ関数（Gamma function）
ガンマ関数の性質
ガウスの乗法公式
倍数公式
ベータ関数（Euler Beta function）
ウォリス積分

ガンマ関数（Gamma function）

ガンマ関数（第２種オイラー積分ともよばれます）は次の積分によって定義されます。

$\qquad \displaystyle{\Gamma(x)=\int_{0}^{\infty} t^{x-1} e^{-t} \mathrm{~d} t\qquad\left( {\rm Re}\,x>0\right)}$

ガンマ関数のSymPyでの呼び出しコマンドは

sp.var('x')
sp.gamma(x)

です。あるいは定義式の右辺の積分を計算すると

sp.var('x, t', positive=True )
sp.integrate(t**(x-1)*sp.exp(-t), (t, 0, sp.oo))

から、積分結果としてガンマ関数を出力します。

ガンマ関数の性質

ガンマ関数の特徴的な性質として

$\qquad \displaystyle{\Gamma(x+1)=x\,\Gamma(x)}$

という恒等式をみたします。この関係式は次のように部分積分することで得られます：

$\qquad\displaystyle{\begin{aligned}\Gamma(x+1)&=\int_0^\infty t^{x}e^{-t}dt\\ &=\left[t^{x}(-e^{-t})\right]_{t=0}^\infty-\int_0^\infty (xt^{x-1})(-e^{-t})dt\\ &=x\,\Gamma(x)\end{aligned}}$

SymPyでたしかめてみると

sp.simplify(sp.gamma(x+1)/sp.gamma(x))

です。とくに $z=n$ （自然数）とするとき、関係式を逐次的に用いて

$\qquad \displaystyle{\Gamma(n+1)=n\,\Gamma(n)=n(n-1)\Gamma(n-1)=\dots=n!\,\Gamma(1)}$

がわかります。ここで

sp.gamma(1)

から

$\qquad\displaystyle{\Gamma(1)=1}$

より

$\qquad\displaystyle{\Gamma(n+1)=n!}$

となります。よってGamma関数は階乗を非整数に一般化したものだと思えます。また、半整数値においては

$\qquad \displaystyle{\begin{aligned}\Gamma\left(n+\frac{1}{2}\right)&=\left(n-\frac{1}{2}\right)\Gamma\left(n-\frac{1}{2}\right)\\&=\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right)\Gamma\left(n-\frac{3}{2}\right)\\&=\dots=\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right)\dots\frac32\frac12\Gamma\left(\frac{1}{2}\right)\\&=\frac{(2n-1)!!}{2^{n}}\Gamma\left(\frac{1}{2}\right)\end{aligned}}$

ここで

sp.gamma(sp.Rational(1,2))

から

$\qquad\displaystyle{\Gamma(1/2)=\sqrt{\pi}}$

より

$\qquad\displaystyle{\Gamma\left(n+\frac{1}{2}\right)=\frac{(2n-1)!!\sqrt{\pi}}{2^{n}}}$

となります。

ガウスの乗法公式

ガンマ関数はまた次の極限によっても得られます。

$\qquad \displaystyle{\Gamma(x)=\displaystyle{\lim_{n\to\infty}\frac{n!\,n^{x}}{x(x+1)\dots(x+n)}}}$

この右辺の極限値がガンマ関数であることをSymPyは知っています。

sp.var('k, n')
sp.limit(sp.factorial(n)*n**x/sp.product(x+k, (k, 0, n) ), n, sp.oo)

この公式の詳しい導出については以下の過去記事で解説しています。

pianofisica.hatenablog.com

倍数公式

ガンマ関数で引数を２倍にしたとき、以下の公式が成り立ちます。

$\qquad \displaystyle{\Gamma(2x)=\frac{2^{2x}}{2\sqrt{\pi}}\Gamma(x)\Gamma\left(x+\frac12\right)}$

導出の詳細は以下の過去記事を見てください。

pianofisica.hatenablog.com

ここでは、公式の成立をSymPyで確かめてみましょう。

sp.simplify(sp.gamma(2*x)/(sp.gamma(x)*sp.gamma(x+sp.Rational(1,2))))

ベータ関数（Euler Beta function）

ベータ関数は次の積分で定義される２つの引数をもつ関数です。

$\qquad \displaystyle{\mathrm{B}(x, y)=\int_{0}^{\infty} \frac{t^{y-1}}{(1+t)^{x+y}} \mathrm{~d} t}$

ベータ関数は第１種オイラー積分ともよばれ、第２種オイラー積分（ガンマ関数）と密接に関係しています。実際、ベータ関数はガンマ関数を用いて書き表せます。ベータ関数とガンマ関数を結び付ける関係式は、定積分の計算などに力を発揮します。

変数変換 $s=t/(1+t)$ を施すと上記の積分表示は

$\qquad \displaystyle{\mathrm{B}(x, y)=\int_{0}^{1} \ {s^{y-1}}{(1-s)^{x-1}} \mathrm{~d} s}$

とも表せます。実際

$\qquad \displaystyle{1-s=1-\frac{t}{1+t}=\frac{1}{1+t}}$

から

$\qquad \displaystyle{\frac{ds}{dt}=\frac{1}{(1+t)^2}}$

すなわち

$\qquad \displaystyle{\frac{dt}{ds}=(1+t)^2}$

に注意して

$\qquad\displaystyle{\begin{aligned}B(p,q)&=\int_{0}^{\infty} \frac{t^{y-1}}{(1+t)^{x+y}} \mathrm{~d} t\\&=\int_0^1 ds\,\frac{dt}{ds}\,\frac{s^{y-1}}{(1+t)^{\,x+1}}\\&=\int_0^1 ds\,(1-s)^{-2}\,s^{y-1}(1-s)^{x+1}\\&=\int_0^{1} ds\,s^{y-1}(1-s)^{x-1}\end{aligned}}$

として示せます。さらに $s=\sin^2\theta$ と変数変換して

$\qquad\displaystyle{B\left(\frac{p+1}{2},\frac{q+1}{2}\right)=2\int_0^{\frac{\pi}{2}} d\theta\,\sin^{q}\theta\cos^{p}\theta}$

という第３の積分表示が得られます。

ベータ関数のSymPyでの呼び出しコマンドは

sp.beta(x,y)

です。あるいは定義式の右辺の積分を計算すると

sp.var('x, y, t', positive=True )
sp.integrate(t**(y-1)*(1+t)**(-x-y), (t, 0, sp.oo))

から、積分結果としてガンマ関数の積を出力します。

$\qquad \displaystyle{\int_{0}^{\infty} \frac{t^{y-1}}{(1+t)^{x+y}} \mathrm{~d} t=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}}$

実はベータ関数はガンマ関数によって

$\qquad \displaystyle{\mathrm{B}(x, y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}}$

と表せます。詳細については以下の過去記事を参照してください。

pianofisica.hatenablog.com

ここでは上記の関係式が成り立っていることをSymPyで確かめましょう。

sp.simplify(sp.beta(x,y)/(sp.gamma(x)*sp.gamma(y)/sp.gamma(x+y)))

ウォリス積分

ここまでで示した関係式から、次のウォリス（Wallis）積分の公式を得ます。

$\qquad \displaystyle{\int_0^{\frac{\pi}{2}} \sin^p\theta \,d\theta = \begin{cases}\frac{(2n-1)!!}{(2n)!!}\frac{\pi}{2}\quad (p=2n) \\ \frac{(2n)!!}{(2n+1)!!} \quad (p=2n+1) \end{cases}}$

公式の証明の詳細は上で紹介した過去記事内でみることができます。ここでは $p=5,10$ の２つの場合で公式が成り立っていることを確認してみます。

sp.var('theta')
sp.integrate((sp.sin(theta))**5, (theta, 0, sp.pi/2 ))-sp.Rational(sp.factorial2(2*2),sp.factorial2(2*2+1))

sp.integrate((sp.sin(theta))**10, (theta, 0, sp.pi/2 ))-sp.Rational(sp.factorial2(2*5-1),2*sp.factorial2(2*5))*sp.pi

キーワード：Python、SymPy、特殊関数、ガンマ関数、ガウスの乗法公式、ベータ関数