[SciPy] ベッセル関数とノイマン関数

[SciPy] ベッセル関数とノイマン関数

ベッセル関数

 第1種ベッセル関数 (第1種円柱関数) は、ベッセルの微分方程式
 
\[x^2y''+xy'+(x^2-v^2)y=0\tag{1}\]
の特殊解の1つであり、具体的には
 
\[Jv(x)=\sum_{s=0}^{\infty}\frac{(-1)^2}{s!\Gamma(s+v+1)}
\left(\frac{x}{2}\right)^{v+2s}\tag{2}\]
という級数で表されます。$v$ が非整数のときには、$J_{-v}$ もまた $J_v$ とは独立な解です。したがって、ベッセルの微分方程式の一般解は
 
\[y=AJ_v(x)+BJ_{-v}(x)\tag{3}\]
と表せます。しかし $v$ が整数 $n$ であるときには、
 
\[J_{-n}(x)=(-1)^nJ_n(x)\tag{4}\]
という関係式が成り立つので、$J_n$ と $J_{-n}$ は独立な解ではありません。そこで $v$ が整数であっても $J_v$ と独立となる
 
\[Y_v(x)=\frac{J_v(x)\cos\pi v-J_{-v}(x)}{\sin \pi v}\tag{5}\]
をもう1つの解として採用し、これを 第2種ベッセル関数 (第2種円柱関数) 、または ノイマン関数 とよびます。すなわちベッセル微分方程式の一般解は $J_v(x)$ と $Y_v(x)$ の線型結合
 
\[y=AJ_v(x)+BY_v(x)\tag{6}\]
で表されることになります。ベッセル関数とノイマン関数については、『物理のための応用数学 (小野寺嘉孝)』に詳しい説明があります。

物理のための応用数学

中古価格
¥1,393から
(2019/10/24 10:26時点)

scipy.special.jv()

 scipy.special.jv(v, z) は第1種ベッセル関数の値を返します。
 v には 次数、z には変数 (float または complex) を渡します。
 たとえば、scipy.special.jv(3, 0) は $J_3(0)$ の値を計算します。
 以下のコードでは、Matplotlib を使って、$J_0(x),\ J_1(x),\ J_2(x)$ のグラフを描画します。

# リストBSL_01

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv

# グラフ描画領域を設定
fig = plt.figure(figsize = (8, 5))
ax = fig.add_subplot(111)
ax.set_title("Bessel functions of the first kind", size = 15)
ax.grid()
ax.set_xlim(-10, 10)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("Jv(x)", size = 15, labelpad = 8)

x = np.linspace(-10, 10, 65)

# 第1種ベッセル関数をプロット
for v in range(3):
    ax.plot(x, jv(v, x), label = "v={}".format(v))

# 凡例を表示
ax.legend()

第1種ベッセル関数 (Bessel functions of the first kind)
 

scipy.special.yv()

 scipy.special.yv(v, z) は第2種ベッセル関数の値を返します。
 v には 次数、z には変数 (float または complex) を渡します。たとえば、scipy.special.yv(1, 5) は $Y_1(5)$ の値を計算します。以下のコードでは、Matplotlib を使って第2種ベッセル関数のグラフをプロットします。

# リストBSL_02

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import yv

# グラフ描画領域を設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("Bessel functions of the second kind", size = 15)
ax.grid()
ax.set_xlim(0, 12)
ax.set_ylim(-3, 1)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("Yv(x)", size = 15, labelpad = 8)

# 区間[0,12]を513分割
x = np.linspace(0, 12, 513)

# 第2種ベッセル関数をプロット
for v in range(5):
    ax.plot(x, yv(v, x), label = "v={}".format(v))

# 凡例を表示
ax.legend(loc = "lower right")

第2種ベッセル関数 (Bessel function of the second kind)

 v に integer, z に float を受け取って第2種ベッセル関数を計算する scipy.special.yn(v, z) という関数も用意されています。
 

ハンケル関数

 $J_v(x)$ と $Y_v(x)$ を以下のように結びつけた関数 $H_v^{(1)}(x)$ と $H_v^{(2)}(x)$ を、それぞれ第1種ハンケル関数、第2種ハンケル関数とよびます。
\[\begin{align*}
H_v^{(1)}(x)&=J_v(x)+iY_v(x)\tag{7}\\[6pt]
H_v^{(2)}(x)&=J_v(x)-iY_v(x)\tag{8}\end{align*}\]

scipy.special.hankel1()

 scipy.special.hankel1(v, z) は第1種ハンケル関数 $H_v^{(1)}(z)$ を計算します。

# リストBSL_03

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hankel1

# 第1種ハンケル関数のv=1,x=2における値を計算
hf = hankel1(1, 2)

# hfを丸める
hf = np.round(hf, 3)

print(hf)
(0.577-0.107j)

 

scipy.special.hankel2()

 scipy.special.hankel2(v, z) は第2種ハンケル関数 $H_v^{(2)}(z)$ を計算します。

# リストBSL_04

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hankel1, hankel2

# 第2種ハンケル関数のv=3,x=-1における値を計算
hs = hankel2(3, -1)

# hsを丸める
hs = np.round(hs, 3)

print(hs)
(-0.059-5.822j)

 

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~

新品価格
¥3,218から
(2019/8/6 11:41時点)

変形ベッセル関数

 ベッセルの微分方程式
 
\[x^2y''+xy'+(x^2-v^2)y=0\tag{1}\]
において、$x$ を $ix$ に置き換えると、
 
\[x^2y''+xy'-(x^2+v^2)y=0\tag{9}\]
のように、$x^2$ の係数の符号が反転します。この変形されたベッセル微分方程式の特殊解の1つ
 
\[I_v(x)=i^{-v}J_v(ix)\tag{10}\]
を第1種変形ベッセル関数とよびます。次数 $v$ が非整数のときには、$I_{-v}(x)$ は $I_v(x)$ と独立した解であり、変形ベッセル微分方程式の一般解は
 
\[y=AI_v(x)+BI_{-v}(x)\tag{11}\]
と表せますが、$v$ が整数のときは $I_{-n}(x)=I_n(x)$ という関係が成り立つので、$I_n(x)$ と独立した別の解
 
\[K_v(x)=\frac{\pi}{2}\frac{I_{-v}(x)-I_v(x)}{\sin\pi v}
=\frac{\pi i}{2}i^{v+1}H_v^{(1)}(ix)\tag{12}\]
を用意します。この関数を第2種変形ベッセル関数とよびます。$K_v(x)$ については、$v$ が整数・非整数に関わらず $K_v(x)=K_{-v}(x)$ が成立します。変形ベッセル微分方程式の一般解は
 
\[y=AI_v(x)+BK_v(x)\tag{13}\]
で表されることになります。
 

scipy.special.iv()

 scipy.special.iv(v, z) は第1種変形ベッセル関数 $I_v(z)$ の値を返します。

# リストBSL_05

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import iv

# グラフ描画領域を設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("Modified Bessel functions of the first kind", size = 15)
ax.grid()
ax.set_xlim(0, 8)
ax.set_ylim(0, 400)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("i_v(x)", size = 15, labelpad = 8)

x = np.linspace(0, 8, 65)

# 第1種修正ベッセル関数をプロット
for v in range(5):
    ax.plot(x, iv(v, x), label = "v={}".format(v))

# 凡例を表示
ax.legend(loc = "upper left")

第1種ベッセル関数 (Modified Bessel function of the first kind)

scipy.special.kv()

 scipy.special.kv(v, z) は第2種変形ベッセル関数 $K_v(z)$ の値を返します。

# リストBSL_06

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import kv

# グラフ描画領域を設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("Modified Bessel functions of the second kind", size = 15)
ax.grid()
ax.set_xlim(0, 8)
ax.set_ylim(0, 4)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("k_v(x)", size = 15, labelpad = 8)

x = np.linspace(0, 10, 513)

# 第2種修正ベッセル関数をプロット
for v in range(5):
    ax.plot(x, kv(v, x), label = "v={}".format(v))

# 凡例を表示
ax.legend(loc = "upper right")

変形ベッセル関数 (Modified Bessel function of the second kind)

球ベッセル関数

 球ベッセル関数と球ノイマン関数 (第2種ベッセル関数) は半奇数次のベッセル関数を使って
 
\[\begin{align*}
j_n(x)&=\sqrt{\frac{\pi}{2x}}J_{n+1/2}(x)\tag{14}\\[6pt]
y_n(x)&=(-1)^n\sqrt{\frac{\pi}{2x}}Y_{n+1/2}(x)\tag{15}
\end{align*}\]
と定義されます。量子力学の自由粒子のシュレーディンガー方程式を球座標で解いたときに、動径方向の解としてこれらの関数が現れます。
 

scipy.special.spherical_jn

 scipy.special.spherical_jn(n, z) は球ベッセル関数の値を返します。

# リストBSL_07

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import spherical_jn

# グラフ描画領域を設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("Spherical Bessel functions of the first kind", size = 15)
ax.grid()
ax.set_xlim(0, 12)
ax.set_ylim(-0.4, 1)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("jn(x)", size = 15, labelpad = 8)

x = np.linspace(0, 12, 65)

# 球ベッセル関数をプロット
for n in range(3):
    ax.plot(x, spherical_jn(n, x), label = "n={}".format(n))

# 凡例を表示
ax.legend(loc = "upper right")

球ベッセル関数 (spherical Bessel function)

scipy.special.spherical_yn

 scipy.special.spherical_yn(n, z) は球ノイマン関数の値を返します。

# リストBSL_08

# ライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import spherical_yn

# グラフ描画領域を設定
fig = plt.figure(figsize = (7, 5))
ax = fig.add_subplot(111)
ax.set_title("Spherical Bessel functions of the second kind", size = 15)
ax.grid()
ax.set_xlim(0, 12)
ax.set_ylim(-0.8, 0.4)
ax.set_xlabel("x", size = 15, labelpad = 10)
ax.set_ylabel("yn(x)", size = 15, labelpad = 8)

# 球ノイマン関数をプロット
for n in range(3):
    x = np.linspace(0, 12, 65)
    ax.plot(x, spherical_yn(n, x), label = "n={}".format(n))

# 凡例を表示
ax.legend(loc = "lower right")

球ノイマン関数 (spherical Neumann function)