ベッセル関数とノイマン関数

ベッセル関数(円筒関数)

第１種ベッセル関数 (第１種円柱関数) は、ベッセルの微分方程式
 \[x^2y”+xy’+(x^2-v^2)y=0\tag{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}\]
をもう１つの解として採用し、これを 第２種ベッセル関数 (第２種円柱関数) 、または ノイマン関数 とよびます。すなわちベッセル微分方程式の一般解は $J_v(x)$ と $Y_v(x)$ の線形結合
 \[y=AJ_v(x)+BY_v(x)\tag{6}\]
で表されることになります。

scipy.special.jv()

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

# PYTHON_BESSEL_01

# In[1]

# ライブラリをインポート
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()

3 次元関数 $J_0(x^2+y^2)$ のグラフを描くコードを載せておきます。

# In[2]

from mpl_toolkits.mplot3d import Axes3D

# FigureとAxes
fig = plt.figure(figsize = (10, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("x", size = 15)
ax.set_ylabel("y", size = 15)
ax.set_zlabel("j0(x**2+y**2)", size = 14)

# 分割数
n = 513

# 連続データと格子点を作成
x = np.linspace(-12, 12, n)
y = np.linspace(-12, 12, n)
X, Y = np.meshgrid(x, y)

# 2変数ベッセル関数
Z = jv(0, np.sqrt(X**2 + Y**2))

# 曲面を描画(カラーマップは'plasma_r')
ax.plot_surface(X, Y, Z, cmap="afmhot")

scipy.special.yv()

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

# PYTHON_BESSEL_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")

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

ハンケル関数

$J_v(x)$ と $Y_v(x)$ を以下のように結びつけた関数 $H_v^{(1)}(x)$ と $H_v^{(2)}(x)$ を、それぞれ第１種ハンケル関数、第２種ハンケル関数とよびます。
 \[\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) は第１種ハンケル関数 $H_v^{(1)}(z)$ を計算します。

# PYTHON_BESSEL_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) は第２種ハンケル関数 $H_v^{(2)}(z)$ を計算します。

# PYTHON_BESSEL_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)

変形ベッセル関数

ベッセルの微分方程式
 \[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$ の係数の符号が反転します。この変形されたベッセル微分方程式の特殊解の１つ
 \[I_v(x)=i^{-v}J_v(ix)\tag{10}\]
を第１種変形ベッセル関数とよびます。次数 $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}\]
を用意します。この関数を第２種変形ベッセル関数とよびます。$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) は第１種変形ベッセル関数 $I_v(z)$ の値を返します。

# PYTHON_BESSEL_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")

scipy.special.kv()

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

# PYTHON_BESSEL_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")

球ベッセル関数

球ベッセル関数と球ノイマン関数 (第２種ベッセル関数) は半奇数次のベッセル関数を使って
 \[\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) は球ベッセル関数の値を返します。

# PYTHON_BESSEL_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")

scipy.special.spherical_yn

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

# PYTHON_BESSEL_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")