matplotlibによるガンマ関数の描画

参考：Wikipedia:Gamma Function

ガンマ関数って？

これ。

$\Gamma(n) = (n-1)!$

ただの階乗なんだけど実数や複素数を考慮すると面白いものが見えてくる。

何が見えるの？

一般式
$\Gamma(z)=\int_0^\infty \exp(t)\exp(z-1) dt$

まず実数範囲のグラフ。
$\Gamma(Re)$

普通にラインを引くと発散している箇所がつながってしまうので点描のグラフに差し替えた。

ソースコード

あまりプログラミングらしいことをしていません。
ソースにほとんど色がつかないのはNumPyとMatplotlibの記述力が高いおかげです。

from scipy.special import gamma
from matplotlib import pyplot as plt
import numpy as np

C = 500
x = np.linspace( -5., 5., C)
y = gamma(x)
plt.gca().set_autoscale_on(False)
plt.axis([-5., 5., -5., 5.])
plt.plot(x, y, "bo")
plt.show()

続いて複素数範囲。
$\Gamma(Re + Im)$

複素数平面で色付けされたグラフ。きれいですね！わたしが描いたわけじゃないけど！これはmpmathのサンプルコードなんだけど短すぎてよく分からないんだよね…

ソースコード

from mpmath import *
cplot(gamma, points=10000)

複素数の絶対値
$\|\Gamma(Re+Im)\|$

おもしろいですね！これが描きたかっただけです！おわり！

ソースコード

複素数の扱いが面倒になって変数を余計に定義することにした。ひどいコード＞＜

from scipy.special import gamma
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

C = 500
x = np.linspace( -5., 5., C)
y = np.linspace( -5, 5, C)
_y = np.linspace( -5.j, 5.j, C)
X, Y = np.meshgrid(x, y)
_X, _Y = np.meshgrid(x, _y)
z = gamma(X+_Y)
plt.gca().set_autoscale_on(False)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, Y, np.absolute(z), cmap="autumn")
ax.set_zlim3d(0, 10.)
ax.autoscale(False)
plt.axis([-6., 6., -6., 6.])
plt.show()

トリビア

ガンマ関数の数値計算はランチョスの近似法というのを用いるのが標準だそうです。
Wikipedia:Lanczos approximation