球面調和関数を用いた3次元波動方程式の解【音響理論基礎】

今回は3次元波動方程式の解について考えていきます。

導出の流れをざっくり追っていきますので常微分方程式の解き方など細かいところは記事の最後に挙げる本などを参考にしてください。

3次元波動方程式の導出については以下の記事をご参照ください。

【音響理論基礎】波動方程式の導出

導出後，イメージを掴むためにpythonを使って登場する特殊関数を見てみましょう。

では早速本題です。

極座標における3次元波動方程式解の導出

次の極座標における波動方程式

$$\Delta p(r, \theta, \phi, t) – \frac{1}{c^{2}} \frac{\partial^{2} p(r, \theta, \phi, t)}{\partial t^{2}} = 0 \tag{1}$$

の解を変数分離法を用いて導出していきます。
ここで極座標の$\Delta$は

$$\Delta = \frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right) + \frac{1}{r^{2} \sin \theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial}{\partial \theta}\right) + \frac{1}{r^{2} \sin^{2} \theta}\frac{\partial^{2}}{\partial \phi^{2}}$$

であることに注意してください。※幾つか表現があります

周波数領域で考えてHelmholtz方程式の解を導出する流れになっている本も多いですが，ここでは時間領域で考えましょう。

まず解の形を次の変数分離形で仮定します。

$$p(r, \theta, \phi, t) = R(r) \Theta(\theta) \Phi(\phi) T(t)$$

これを式(1)に代入して$p(r, \theta, \phi, t)$で割ると

$$\frac{1}{r^{2}R}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial R}{\partial r}\right) + \frac{1}{r^{2} \sin \theta \Theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + \frac{1}{r^{2} \sin^{2} \theta \Phi}\frac{\partial^{2} \Phi}{\partial \phi^{2}} = \frac{1}{c^{2}T}\frac{\partial^{2} T}{\partial t^{2}}$$

となります。これが恒等的に成り立つので(左辺)=(右辺)=定数です。すなわち定数を$-k^{2}$とすると

\begin{align}
\frac{1}{r^{2}R}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial R}{\partial r}\right) + \frac{1}{r^{2} \sin \theta \Theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + \frac{1}{r^{2} \sin^{2} \theta \Phi}\frac{\partial^{2} \Phi}{\partial \phi^{2}} = -k^{2} \tag{2} \\
\frac{\partial^{2} T}{\partial t^{2}} = -c^{2} k^{2} T \tag{3}
\end{align}

式(3)の解は$\omega := ck$とすれば

$$T(t) = T_{1} e^{i\omega t} + T_{2} e^{-i\omega t} \tag{4}$$

となります。

今回時間依存項は2番目を採用する事にすれば$T_{1} = 0$です。
これで時刻$t$の関数$T(t)$が決まりましたね。

次は式(2)の両辺に$r^{2}$をかけて

$$\frac{1}{R}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial R}{\partial r}\right) + k^{2} r^{2} = – \left\{ \frac{1}{\sin \theta \Theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + \frac{1}{\sin^{2} \theta \Phi}\frac{\partial^{2} \Phi}{\partial \phi^{2}} \right\}$$

先ほどと同じ議論から両辺は定数となるので定数を$n(n+1)$として整理すると

\begin{align}
\frac{\partial}{\partial r}\left(r^{2} \frac{\partial R}{\partial r}\right) + k^{2} r^{2}R = n(n+1)R \tag{5} \\
\frac{\sin \theta}{\Theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + n(n+1) \sin^{2} \theta = – \frac{1}{\Phi}\frac{\partial^{2} \Phi}{\partial \phi^{2}} \tag{6}
\end{align}

とできます。式(5)が$r$だけの式になっているので解を求めていきます。$r^{2}$で割ると

$$\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial R}{\partial r}\right) + k^{2}R – \frac{n(n+1)}{r^{2}} R$$

この形の微分方程式の解は第1種及び第2種の球ベッセル関数を用いて

$$R(r) = R_{1} j_{n}(kr) + R_{2} y_{n}(kr) \tag{7}$$

と書けます。若しくは第1種及び第2種の球ハンケル関数を用いて

$$R(r) = R_{3} h_{n}^{(1)}(kr) + R_{4} h_{n}^{(2)}(kr) \tag{8}$$

と書けます。

最後に角度方向の関数形を求めていきましょう。
式(6)も両辺ともに定数となりますので，その定数を$m^{2}$として

\begin{align}
\frac{\sin \theta}{\Theta}\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + n(n+1) \sin^{2} \theta = m^{2} \tag{9} \\
\frac{\partial^{2} \Phi}{\partial \phi^{2}} = -m^{2} \Phi \tag{10}
\end{align}

と整理できます。

式(10)は先ほどの調和振動子の方程式と同じ形なのでその解は

$$\Phi(\phi) = \Phi_{1} e^{imt} + \Phi_{2} e^{-imt} \tag{11}$$

ですね。
ここで$\Phi(\phi)$は当然周期$2\pi$の関数のはずなので$m \in \mathbb{Z}$である事に注意しましょう。

$\Theta(\theta)$について考えていきましょう。
式(9)を整理して

$$\sin \theta\frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \Theta}{\partial \theta}\right) + \left\{n(n+1) \sin^{2} \theta – m^{2}\right\} \Theta = 0$$

ここで$\eta := \cos \theta$として何やかんや整理すると

$$\frac{d}{d\eta} \left\{(1-\eta^{2}) \frac{d\Theta}{d\eta}\right\} + \left\{n(n+1) – \frac{m^{2}}{1-\eta^{2}}\right\} \Theta = 0 \tag{12}$$

とできます。
この形の微分方程式の解は第1種及び第2種ルジャンドル関数より与えられて

$$\Theta(\theta) = \Theta_{1} P_{n}^{m}(\cos \theta) + \Theta_{2} Q_{n}^{m}(\cos \theta) \tag{13}$$

となります。

ちなみに$Q_{n}^{m}(\eta)$は$\eta=\pm 1$で発散する関数で，$P_{n}^{m}(\eta)$も$n\notin\mathbb{Z}$であれば$\eta= 1$で発散します。

また，$m > n$で$P_{n}^{m}(\eta)=0$です。

角度関数$\Theta(\theta)\Phi(\phi)$を纏めて

$$Y_{n}^{m}(\theta, \phi) = \sqrt{\frac{2n+1}{4\pi} \frac{(n-|m|)!}{n+|m|)!}} P_{n}^{|m|}(\cos \theta) e^{im\phi} \tag{14}$$

と書いてこれを球面調和関数と呼びます。

式(4)，(7)，(8)，(14)を纏めて最終的に3次元波動方程式の解は

$$p(r, \theta, \phi, k) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \left\{ A_{n}^{m} j_{n}(kr) + B_{n}^{m} y_{n}(kr) \right\} Y_{n}^{m}(\theta, \phi) \tag{15}$$

とか

$$p(r, \theta, \phi, k) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \left\{ C_{n}^{m} h_{n}^{(1)} + D_{n}^{m} h_{n}^{(2)}(kr) \right\} Y_{n}^{m}(\theta, \phi) \tag{16}$$

と書くことが出来ます。

ちなみに式(15)は物理的に定在波の解を，式(16)は進行波の解を表しています。
これでざっくりとした導出完了です！

登場する特殊関数の可視化

この記事をご覧になっている方で特殊関数にとっても馴染みがあるという方は少ないのではないでしょうか？

僕もこの辺りの分野を勉強するまでは全然知らなくて凄くとっつき難かった覚えがあります。

こういったよくわからない関数を理解するための近道は目で見てイメージをつけることだと思いますので，ここでは本記事で登場した特殊関数をpythonを使って可視化してみましょう。

pythonのscipyには結構な数の特殊関数が入っていますので簡単に実装できます。

以降プログラムコード→結果という形で進めていきます。

球ベッセル関数

まずは動径方向の関数である球ベッセル関数です。

ここでは低次の幾つかだけをお示ししますが，同じようにして更に高次のものも可視化出来ますのでお試しください。
今回書いたコードは以下です。

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import spherical_jn

N = 3 # order
x_lim = 14
step = 0.1
x = np.arange(0,x_lim,step,dtype=np.float64)
order_list = range(N+1)

fig, ax = plt.subplots()
ax.set_title(r"Spherical Bessel functions $\{j_{n}(x)\}$")
ax.set_xlabel("x")
ax.set_ylabel("$\{j_{n}(x)\}$")
ax.grid()
ax.plot([0,x_lim], [0,0], color="black")
for n in order_list:
    ax.plot(x, spherical_jn(n,x), label="$j_{%d}(x)$" %(n))
    ax.legend()

出力は以下のようなものになります。

これを見ると変数が大きくなるにつれて値が小さくなっていくのがわかります。
次数が上がったときにゼロ点が交互に出てきているのも面白いですね。

球ノイマン関数

次は球ノイマン関数(第2種球ベッセル関数)です。
同様にコードと出力結果をお示しします。

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import spherical_yn

N = 3 # order
x_lim = 14
step = 0.1
x = np.arange(0,x_lim,step,dtype=np.float64)
order_list = range(N+1)
#order_list = [0,2,4,6]

fig, ax = plt.subplots()
ax.set_title(r"Spherical Neumann functions $\{y_{n}(x)\}$")
ax.set_xlabel("x")
ax.set_ylabel("$\{y_{n}(x)\}$")
ax.grid()
ax.set_ylim([-0.8,0.5])
ax.plot([0,x_lim], [0,0], color="black")
for n in order_list:
    ax.plot(x, spherical_yn(n,x), label="$y_{%d}(x)$" %(n))
    ax.legend()

出力結果は以下です。

ここで重要なのは$x=0$で値が発散しているところです。

動径関数は他に球ハンケル関数がありますが複素領域にあって分かりづらいので本記事では省略します。

球面調和関数

最後に球面調和関数です。

この球面調和関数は$L^{2}(S^{2})$の正規直交基底でもあります。

複素関数ですので描画する場合には少し工夫が必要です。 ※実数値球面調和関数というのもあります

今回は各$n$に対して$m$の方が負の場合については虚部を，

$m$が正の場合については実部を取る事で可視化する事にします。

from tqdm import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.special import sph_harm as _sph_harm

def sph_harm(m, n, theta, phi):
    return _sph_harm(m, n, phi, theta)

N = 3
row = N+1
colmn = 2*N+1
theta_num = 50
phi_num = 100
theta, phi = np.mgrid[0:np.pi:1j*theta_num, 0:2*np.pi:1j*phi_num]
fig_scale = 1.5
color_list = ["royalblue", "tomato"]
colormap = mpl.colors.ListedColormap(color_list)

fig = plt.figure(dpi=100, figsize=(colmn*fig_scale,row*fig_scale))
for n in tqdm(range(N+1)):
    for m in range(-n, n+1):
        if m &amp;lt; 0:
            r = np.imag(sph_harm(m, n, theta, phi))
        else:
            r = np.real(sph_harm(m, n, theta, phi))
        x = np.abs(r) * np.sin(theta) * np.cos(phi)
        y = np.abs(r) * np.sin(theta) * np.sin(phi)
        z = np.abs(r) * np.cos(theta)
        fcolors = np.sign(r)
        ax = fig.add_subplot(row, colmn, n*colmn+(m+N+1), projection="3d")
        ax.plot_surface(x, y, z, facecolors=colormap(fcolors), shade=True, rstride=1, cstride=1)
        if m &amp;lt; 0:
            ax.set_title(r"$Im[\it{Y_{%d}^{%d}}(\theta, \phi)]$" %(n,m))
        elif m &amp;gt; 0:
            ax.set_title(r"$Re[\it{Y_{%d}^{%d}}(\theta, \phi)]$" %(n,m))
        else:
            ax.set_title(r"$\it{Y_{%d}^{%d}}(\theta, \phi)$" %(n,m))
        # limit = np.max([x.max(),y.max(),z.max()])
        limit = np.sqrt((2 *n+1)/(8*np.pi)) 
        if m==0:
            limit *= np.sqrt(2)
        ax.set_xlim(-limit,limit)
        ax.set_ylim(-limit,limit)
        ax.set_zlim(-limit,limit)
        ax.set_axis_off()
fig.tight_layout()
fig.savefig("./spherical_harmonics.png")

以下が出力結果です。