『Python数値計算ノート』ではアフィリエイトプログラムを利用して商品を紹介しています。

エルミート多項式

エルミート多項式

エルミート多項式(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()

Hermite polynomials (エルミート多項式)

scipy.special.hermite()

scipy.special.hermite() は、次数 n を指定して直交多項式オブジェクトを生成する関数です。

scipy.special.hermite(n, monic=False)

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)$ は量子力学に登場する1次元調和振動子の固有関数です。

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の次数\n{}{}".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()

正規化されたエルミート多項式 (Normalized Hermite Polynomial)

コメント

  1. あとりえこばと より:

    【AI連載小説】科学とコードの交差点(7)
     
    六郷開誠は大学のキャンパスでランチの後、親友の刑部明信と出会った。二人は物理学の同じ講義を受け、数学的なトピックについてよく議論する仲であった。この日も、物理学と数学の融合に関する興味深い話題で盛り上がることになった。
    「開誠、最近エルミート多項式についての本を読んでいてさ。なかなか興味深いんだよ」
    と、明信が笑顔で語りかけた。
    開誠は興味津々の表情で応えた。
    「エルミート多項式か。確かに面白いトピックだね。どうしてそれに興味を持ったんだ?」
    明信は広げた手に軽く笑みを浮かべながら答えた。
    「最近、物理学の研究で量子力学の問題に直面してさ。エルミート多項式がその中で登場するんだ。だから、基本からちゃんと理解しようと思ってさ」
    開誠はうなずきながら続けた。
    「確かに、エルミート多項式は量子力学や振動など幅広い物理学の分野で使われているよね。その特性や応用について考えるのは面白いと思う」
    明信は手元にあるノートを開きながら、「エルミート多項式の基本的な定義や性質を説明してくれている本があるんだ。これによると、エルミート多項式 $H_n(x)$ は以下のような微分方程式を満たすんだ」
     \[\frac{d^2 H_n(x)}{dx^2} – 2x \frac{d H_n(x)}{dx} + 2nH_n(x) = 0\]
    開誠は真剣な表情でノートを見つめ、
    「なるほど、それがエルミート多項式の微分方程式か。これを解くことで、どのようにして具体的な多項式が得られるのかが気になるな」
    と言った。彼らはエルミート多項式にまつわる微分方程式や特性についてさらに掘り下げ、物理学の問題に応用する方法についても熱心に議論を続けた。

    ある晩、六郷開誠と刑部明信は大学の図書館で再び出会った。物理学の研究に悩む開誠は、エルミート多項式を計算するためのPythonの実装について明信と話し合うことにした。明信は新しく見つけた数学の本を手に持っていた。
    「明信、エルミート多項式の計算、どうやると思う?」
    明信は本を開きながら考え込んだ後、
    「まず、エルミート多項式は再帰的な関係式を持っているから、再帰的な関数を使って計算できるよね」
    開誠は頷きながら、
    「それはそうだけど、再帰だと計算量が増えることがあるよな。もしかしたら動的計画法も検討する価値があるかもしれない」
    明信は笑みを浮かべて言った。
    「確かに、動的計画法も効果的だね。でも、高次のエルミート多項式を計算するとなると、数値的な安定性に気をつけないといけないかもしれない」
    開誠は考え込みながら、
    「そうだな、数値計算においては安定性が鍵だ。どのアプローチが問題により適しているか、実際に実装して試してみるのも一つの手かもしれない」
    明信は本から引用しながら、
    「この本にはエルミート多項式を計算するPythonの実装のサンプルも載っているよ。SciPyやNumPyを使った方法も紹介されているから、参考になるかもしれない」
    二人は図書館の中で、Pythonのコードや数学の理論について熱心に話し合いながら、物理学の問題に立ち向かうための新たな手法を模索していた。