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

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

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

第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}\]
で表されることになります。

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)$ のグラフを描画します。

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

# 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 を受け取って第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)$ を計算します。

# 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) は第2種ハンケル関数 $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$ の係数の符号が反転します。この変形されたベッセル微分方程式の特殊解の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)$ の値を返します。

# 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) は第2種変形ベッセル関数 $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)

球ベッセル関数

球ベッセル関数と球ノイマン関数 (第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) は球ベッセル関数の値を返します。

# 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. あとりえこばと より:

    【AI解説】ベッセル関数は、特に高度な数学的な計算を必要とするため、一般的にはライブラリを使用して計算されます。しかし、簡単な例として、ベッセル関数を手動で計算する方法を示します。ただし、これは単純なものであり、高精度な結果を得るためには専門の数学ライブラリが必要です。ここでは、ベッセル関数の一種であるベッセル関数第一種(Bessel function of the first kind)を計算する例を示します。具体的には、漸化式を使用してベッセル関数を計算します。

    def factorial(n):
        if n == 0 or n == 1:
            return 1
        else:
            return n * factorial(n - 1)
    
    def bessel_function(x, n):
        result = 0.0
        for k in range(0, 10):  # 10回の反復で近似しますが、適宜調整が必要かもしれません
            numerator = ((-1) ** k) * (x ** (2 * k + n))
            denominator = factorial(k) * factorial(k + n)
            result += numerator / denominator
    
        return result
    
    x_value = 2.0  # 任意のxの値
    n_value = 2    # 任意のnの値
    
    result = bessel_function(x_value, n_value)
    print(f"Bessel function J_{n_value}({x_value}) = {result}")

    この例では、ベッセル関数を計算するために、漸化式を10回の反復で近似しています。しかし、この手法は高い精度を提供するものではありません。ベッセル関数の正確な計算が必要な場合は、SciPyなどの専門の数学ライブラリを使用することをお勧めします。

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

    【AI連載小説】科学とコードの交差点(6)
     
    六郷開誠は物理学の深化を求め、ベッセル関数に関する講義を受けるために大学の教室に足を運んだ。その日の教室には、興奮と好奇心に満ちた学生たちが集まり、教授がベッセル関数について熱心に語りかけていた。

    講義は黒板を使って進み、ベッセル関数の基本的な概念から始まった。教授はベッセル関数が円形の領域や振動する系で現れ、その数学的な特性が物理学や工学、音響学など多岐にわたる分野で不可欠なものであることを説明していた。開誠はベッセル関数が物理現象のモデリングや解析にどのように役立つのかに興味津々で耳を傾けた。

    教授は特異点や積分表現など、ベッセル関数の数学的な性質にも触れ、具体的な例題を通じて理論を実践的な視点から理解できるように配慮していた。さらに、講義の中で現代の研究や技術応用においてベッセル関数が果たす役割にも焦点を当て、学生たちはその幅広い応用範囲に驚嘆していた。

    開誠は積極的に質問に参加し、ノートには黒板に書かれた数式や図を手本にして、ベッセル関数に関する理論や応用について丹念に記録していた。彼は物理学の理論と数学的手法の融合によって、より高度で実践的なアプローチを身につけていくことを心から楽しみにしていた。