エルミート多項式
エルミート多項式 (Hermite polynomial) は $g(t,x)=\exp(2tx-t^2)$ を指数型母関数として生成される直交多項式です。すなわち、
\[g(t,x)=\sum_{n=0}^{\infty}H_n(x)\frac{t^n}{n!}\]
と展開したときの係数をエルミート多項式 $H_n(x)$ と定義します。エルミート多項式は以下の漸化式を満たします。
\[\begin{align*}&H_{n+1}(x)=2xH_n(x)-2nH_{n-1}(x)\\[6pt]&H'_n(x)=2nH_{n-1}(x)\\[6pt]&H'_n(x)=2xH_n(x)-H_{n+1}(x)\end{align*}\]
この漸化式を使って $H_n(x)$ の具体的な表式を求めると
\[\begin{align*}H_0(x)&=1\\[6pt]H_1(x)&=2x\\[6pt]H_2(x)&=4x^2-2\\[6pt]H_3(x)&=8x^3-12x\\[6pt]H_4(x)&=16x^4-48x^2+12\end{align*}\]
のようになります。
scipy.special.eval_hermite()
scipy.special.eval_hermite() は任意の点におけるエルミート多項式の値を返します。
scipy.special.eval_hermite(n, x, out=None)
n には多項式の次数、x には配列もしくは配列に相当するオブジェクトを渡します。
# SCIPY_EVAL_HERMITE
# In[1]
import numpy as np
from scipy.special import eval_hermite
# 3次のエルミート多項式のx=0, 1, 2, 3, 4における値を計算
x = np.arange(5)
h = eval_hermite(3, x)
print(h)[ -0. -4. 40. 180. 464.]
$n=0,\ 1,\ 2,\ 3$ のエルミート多項式のグラフをプロットしてみます。
# In[2]
import matplotlib.pyplot as plt
# x座標データ
x = np.linspace(-4, 4, 129)
# グラフ描画領域を設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_title("Hermite polynomials", size = 15)
ax.grid()
ax.set_xlim(-4, 4)
ax.set_ylim(-10,10)
ax.set_xlabel("x", size = 14, labelpad = 8)
ax.set_ylabel("Hn(x)", size = 14, labelpad = 8)
# 色のリスト
c = ["red", "blue", "green", "darkorange"]
# H0,H1,H2,H3のグラフを描画
for i in range(4):
ax.plot(x, eval_hermite(i, x),
label = "n = {}".format(i), color = c[i-1])
# 凡例を表示
ax.legend()
plt.show()
scipy.special.hermite()
scipy.special.hermite() は、次数 n を指定して直交多項式オブジェクトを生成する関数です。
4 次のエルミート多項式を作ってみます。
# SCIPY_HERMITE
# In[1]
import numpy as np
from scipy.special import hermite
# 4次のエルミート多項式
h4 = hermite(4)
print(h4)4 2 16 x - 48 x - 8.882e-16 x + 12
生成されたインスタンスの正式なクラス名を確認しておきます。
# In[2]
# h4のクラスを表示
print(type(h4))<class 'scipy.special.orthogonal.orthopoly1d'>
生成したオブジェクトから、任意の点における値を取得することができます。
# In[3]
# x=[0 1 2 3 4 5]
x = np.arange(0, 6)
# xの各要素におけるh4の値を表示
print(h4(x))[12. -20. 76. 876. 3340. 8812.]
各項の係数はインスタンス変数 c に格納されています。
# In[4]
# h4の各項の係数を表示
print(h4.c)[1.6000000e+01 0.0000000e+00 -4.8000000e+01 -8.8817842e-16 1.2000000e+01]
多項式の根 (方程式 h4=0 の解) を得ることもできます。
# In[5]
# h4の根を表示
print(h4.r)[-1.65068012 1.65068012 -0.52464762 0.52464762]
正規化されたエルミート多項式
エルミート多項式に重み関数 $e^{-x^2/2}$ を乗じて、$\sqrt{\sqrt{\pi}2^nn!}$ で割った関数
\[u_n(x)=\frac{1}{\sqrt{\sqrt{\pi}2^nn!}}H_n(x)\exp\left(-\frac{x^2}{2}\right)\]
は正規直交系をなします。すなわち $u_n$ を全区間で積分すると
\[\int_{-\infty}^{\infty}u_m(x)u_n(x)=\delta_{mn}\]
となります ($\delta_{mn}$ は、$m=n$ のときに $1$、それ以外は $0$ の値をとるクロネッカーのデルタ記号です)。$u_n(x)$ は量子力学に登場する１次元調和振動子の固有関数です。
scipy.special.hermitenorm()
scipy.special.hermitenorm() は正規化されたエルミート多項式 $u_n$ を返します。
scipy.special.hermitenorm(n, monic=0)
5 次の正規化されたエルミート多項式を生成してみます。
# SCIPY_HERMITE_NORM
# In[1]
import numpy as np
from scipy.special import hermitenorm
# 正規化された5次のエルミート多項式
u5 = hermitenorm(5)
print(u5)5 3 1 x - 10 x + 15 x
u5 の次数、係数、根を調べてみます。
# In[2]
# u5の次数
print("u5の次数 {}".format(u5.order))
# u5の各項の係数
print("\nu5の係数\n{}".format(u5.c))
# u5=0の根
print("\nu5の根\n{}".format(u5.r))u5の次数 5 u5の係数 [1. 0. -10. 0. 15. 0.] u5の根 [-2.85697001 2.85697001 -1.35562618 1.35562618 0.]
正規化されたエルミート多項式のグラフを描いてみます。
# In[3]
import matplotlib.pyplot as plt
# x座標データ
x = np.linspace(-4, 4, 129)
# グラフ描画領域を設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_title("Normalized Hermite Polynomial", size = 14)
ax.grid()
ax.set_xlim(-4, 4)
ax.set_ylim(-10,10)
ax.set_xlabel("x", size = 14, labelpad = 8)
ax.set_ylabel("un(x)", size = 14, labelpad = 8)
# 色のリスト
c = ["red", "blue", "green", "darkorange"]
# u0,u1,u2,u3のグラフを描画
for i in range(4):
u = hermitenorm(i)
ax.plot(x, u(x),
label = "n = {}".format(i), color = c[i-1])
# 凡例を表示
ax.legend()
plt.show()