ベッセル関数とノイマン関数 │ Python 数値計算入門

1. ベッセル関数(円筒関数)
- 1.1. scipy.special.jv()
- 1.2. scipy.special.yv()
2. ハンケル関数
- 2.1. scipy.special.hankel1()
- 2.2. scipy.special.hankel2()
3. 変形ベッセル関数
- 3.1. scipy.special.iv()
- 3.2. scipy.special.kv()
4. 球ベッセル関数
- 4.1. scipy.special.spherical_jn
- 4.2. scipy.special.spherical_yn

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

　第１種ベッセル関数 (第１種円柱関数) は、ベッセルの微分方程式
　
\[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()

第1種ベッセル関数 (Bessel functions of the first kind)
　
　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")

2変数ベッセル関数のグラフ

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

第2種ベッセル関数 (Bessel function of the second kind)
　
　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")

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

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

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

球ベッセル関数

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

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

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

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

コメント（1件）

あとりえこばとより:

2023年5月15日 5:58 PM

人工知能Bardがベッセル関数についてやさしく解説します

　
[Bard] ベッセル関数は、円筒関数とも呼ばれ、円筒座標系の微分方程式を満たす関数の一種です。それらは、電磁波と音波の伝播、流体の流れ、振動などの多くの物理現象を記述するために使用されます。ベッセル関数はまた、数学や工学でも多くの用途があります。ベッセル関数はドイツの数学者フリードリヒ・ヴィルヘルム・ベッセルにちなんで名付けられました。ベッセルは1824年に彼の関数を最初に研究しました。彼は、それらが太陽系の惑星の軌道を記述するために使用できることを発見しました。ベッセル関数は、1839年にスイスの数学者ダニエル・ベルヌーイによって独立して再発見されました。ベルヌーイは、それらが電磁波の伝播を記述するために使用できることを発見しました。ベッセル関数は、それ以来、多くの物理現象を記述するために使用されてきました。それらは今日でも研究と応用が行われている重要な関数です。

返信