とつとつとしてろうとせず

ひまつぶしにどうぞ。

次数の異なる2つの変形Bessel関数の比の近似

はじめに

変形Bessel関数はBessel関数と同様、円柱座標系や平面極座標系における偏微分方程式を変数分離法によって解いたときの半径方向の分布を表す関数として知られている。たとえば、平面極座標系$(r,\theta)$における変形Helmholtz方程式
\[ \paren{\nabla^2 - \alpha^2}\psi = 0 \]の一般解は、次式のような関数の重ね合わせによって表される($\alpha$は定数)。
\[ \psi_{n}(r,\theta) = \paren{ A_{n} I_{n}(\alpha r) + B_{n} K_{n}(\alpha r) } e^{i n \theta} \]ここで、$I_{n}(r)$と$K_{n}(r)$はそれぞれ$n$次の第一種および第二種の変形Bessel関数であり、$A_n$と$B_n$は境界条件から定まる未定係数である。$n$は分離定数であり、周期境界条件を考慮する場合は整数である。もし解析領域が単位円盤であり、$r = 1$でNeumann境界条件が与えられたとすると、未定係数は次のように次数の隣り合う2つの変形Bessel関数の比を含んだ形になる(変形Bessel関数の微分は異なる次数の変形Bessel関数との漸化式によって与えられるため)。
\[ \frac{I_{n \pm 1}(\alpha) }{I_{n}(\alpha) } \ ,\ \frac{K_{n \pm 1}(\alpha) }{K_{n}(\alpha) } \]一般に、異なる次数$\nu \neq \mu$を持つ同種の変形Bessel関数の比は、次数を固定して変数$x$を大きくしていくと、$x \rightarrow \infty$の極限で$1$に収束する。
\[ \lim_{x \rightarrow \infty} \frac{I_{\mu}(x) }{I_{\nu}(x) } = 1 \ ,\ \lim_{x \rightarrow \infty} \frac{K_{\mu}(x) }{K_{\nu}(x) } = 1 \]例として、次数が$0$と$1$の組み合わせについてグラフをプロットすると、下図のようになる。

一方、変数を固定して次数を大きくしていった場合、$\nu \rightarrow \infty$における比の極限値は次数の組み合わせによって$0$または$\infty$となる。例えば、$\mu = \nu + n$、$n > 0$とすると、次式のようになる。
\begin{align*}
\lim_{\nu \rightarrow \infty} \frac{I_{\nu + n}(x) }{I_{\nu}(x) } = \lim_{\nu \rightarrow \infty} \frac{K_{\nu}(x) }{K_{\nu + n}(x) } = 0 \ ,\ \lim_{\nu \rightarrow \infty} \frac{I_{\nu}(x) }{I_{\nu + n}(x) } = \lim_{\nu \rightarrow \infty} \frac{K_{\nu + n}(x) }{K_{\nu}(x) } = \infty
\end{align*}例として、次数の組み合わせが$\nu$と$\nu+1$の場合のグラフを示す。

どちらの場合であっても、変形Bessel関数は原点もしくは無限遠で発散する関数であり、数値的な取り扱いが難しい。また、特殊関数の厳密値の計算にはライブラリやソフトウェアの利用が必要であり、時間や手間がかかる。そこで、漸近べき級数を用いて2つの変形Bessel関数の比の近似値を求めることを考える。なお、ここでは主に変数を固定して次数を大きくする場合を想定している。

変形Bessel関数の漸近級数展開

$x$を実変数、$\nu$を実数の次数としたとき、$x \rightarrow \infty$における$I_{\nu}(x)$および$K_{\nu}(x)$の漸近級数展開は次式のように表される。
\begin{align}
& I_{\nu}(x) \sim \frac{e^{x}}{\sqrt{2 \pi x}} \sum_{k=0}^{\infty} \paren{-1}^{k} \frac{a_{k}(\nu)}{x^{k}} \qquad x \rightarrow \infty \\
& K_{\nu}(x) \sim \sqrt{\frac{\pi}{2 x}} e^{-x} \sum_{k=0}^{\infty} \frac{a_{k}(\nu)}{x^{k}} \qquad x \rightarrow \infty
\end{align}展開係数$a_{k}(\nu)$は次式で表される。
\begin{align}
a_{0}(\nu) = 1 \ ,\ a_{k}(\nu) = \frac{1}{k!8^{k}} \prod_{i=1}^{k} \paren{4\nu^{2}-\paren{2i-1}^2}
\end{align}最初の数項を具体的に書き表すと、次のようになる。
\begin{align}
\label{eq:asympt I}
& I_{\nu}(x) \sim \frac{e^{x}}{\sqrt{2 \pi x}} \paren{ 1 - \frac{4\nu^{2}-1}{8x} + \frac{\paren{4\nu^{2}-1}\paren{4\nu^{2}-9}}{2! \paren{8x}^{2}} - \frac{\paren{4\nu^{2}-1}\paren{4\nu^{2}-9}\paren{4\nu^{2}-25}}{3! \paren{8x}^{3}} +- \cdots } \qquad x \rightarrow \infty \\
\label{eq:asympt K}
& K_{\nu}(x) \sim \sqrt{\frac{\pi}{2 x}} e^{-x} \paren{ 1 + \frac{4\nu^{2}-1}{8x} + \frac{\paren{4\nu^{2}-1}\paren{4\nu^{2}-9}}{2! \paren{8x}^{2}} + \frac{\paren{4\nu^{2}-1}\paren{4\nu^{2}-9}\paren{4\nu^{2}-25}}{3! \paren{8x}^{3}} + \cdots } \qquad x \rightarrow \infty
\end{align}これより、変形Bessel関数は$x \rightarrow \infty$のとき、主要項$e^{x} / \sqrt{2 \pi x}$、$e^{-x} / \sqrt{2 \pi^{-1} x}$と$1/x$のべき級数の積で表される。また、次数の異なる2つの変形Bessel関数の比の漸近展開は、次数に無関係な主要項の部分が打ち消され、$1/x$のべき級数のみで近似できることがわかる。
ちなみに$\nu \rightarrow \infty$に対する漸近級数展開の主要項は次式で与えられる($x \neq 0$)。
\begin{align}
& I_{\nu}(x) \sim \frac{1}{\sqrt{2 \pi \nu}} \paren{ \frac{e x}{2 \nu} }^{\nu} \\
& K_{\nu}(x) \sim \sqrt{\frac{\pi}{2 \nu}} \paren{ \frac{e x}{2 \nu} }^{-\nu}
\end{align}

漸近べき級数の四則演算

一般に、ある関数$y(x)$の漸近級数
\[ y(x) \sim \sum_{n} a_{n} \varphi_{n}(x) \]のように表され、Taylor展開のように常にべき級数で表されるわけではない。しかし、もしも2つの関数$f(x)$、$g(x)$が
\[ f(x) \sim \sum_{n=0}^{\infty} a_{n} x^{n} \ ,\ g(x) \sim \sum_{n=0}^{\infty} b_{n} x^{n} \]ような漸近べき級数を持つとすると、それらの間の四則演算は通常の級数のように行うことができる。
\begin{align*}
& \alpha f(x) + \beta g(x) \sim \sum_{n=0}^{\infty} \paren{ \alpha a_{n} + \beta b_{n} } x^{n} \\
& f(x) g(x) \sim \sum_{n=0}^{\infty} c_{n} x^{n} \\
& \frac{f(x)}{g(x)} \sim \sum_{n=0}^{\infty} d_{n} x^{n}
\end{align*}ここで、$\alpha$、$\beta$は定数であり、積と商の係数$c_{n}$および$d_{n}$は次式で与えられる。
\begin{align}
c_{n} &= \sum_{m=0}^{n} a_{n} b_{n-m} \\
\label{eq:coef d}
d_{n} &= \frac{1}{b_{0}} \paren{ a_{n} - \sum_{m=0}^{n-1} d_{m} b_{n-m} }
\end{align}ただし、$b_{0} \neq 0$、$d_{0} = a_{0} / b_{0}$である。これらの四則演算は$1/x$のべき級数に対しても成り立つ。したがって、
\[ I_{\nu}(x) \sim \frac{e^{x}}{\sqrt{2 \pi x}} \sum_{k=0}^{\infty} (-1)^{k} a_{k}(\nu) x^{-k} \ ,\ K_{\nu}(x) \sim \sqrt{\frac{\pi}{2 x}} e^{-x} \sum_{k=0}^{\infty} a_{k}(\nu) x^{-k} \]とおくと、2つの変形Bessel関数の比の漸近べき級数は同じ形式で次式のように表される。
\begin{align}
\frac{Z_{\mu}(x) }{Z_{\nu}(x) } &\sim \sum_{k=0}^{\infty} \delta_{k}(\nu,\mu) x^{-k} = \sum_{k=0}^{\infty} \frac{1}{\alpha_{0}(\nu)} \paren{ \alpha_{k}(\mu) - \sum_{m=0}^{k - 1} \delta_{m}(\nu,\mu) \alpha_{k-m}(\nu) } x^{-k}
\end{align}ここで、$Z_{\nu}(x)$は$I_{\nu}(x)$または$K_{\nu}(x)$であり、$I_{\nu}(x)$のとき$\alpha_{k}(\nu) = (-1)^{k} a_{k}(\nu)$、$K_{\nu}(x)$のとき$\alpha_{k}(\nu) = a_{k}(\nu)$である。また、$\delta_{k}(\nu,\mu)$は$\alpha_{k}(\nu)$と$\alpha_{k}(\mu)$を含む再帰的に求められる係数である。
$I_{\mu}(x) / I_{\nu}(x)$に対応する具体的な係数$\delta_{k}(\nu,\mu)$を4項目($k = 0 , 1, 2, 3$)まで表すと、次のようになる。
\begin{align*}
&\delta_{0}(\nu,\mu) = 1 \\
& \delta_{1}(\nu,\mu) = - \paren{ a_{1}(\mu) - a_{1}(\nu) } \\
& \delta_{2}(\nu,\mu) = a_{2}(\mu) - \paren{ a_{2}(\nu) + \delta_{1}(\nu,\mu) } \\
& \delta_{3}(\nu,\mu) = - a_{3}(\mu) - \paren{ - a_{3}(\nu) + \delta_{1}(\nu,\mu)a_{2}(\nu) - \delta_{2}(\nu,\mu)a_{1}(\nu) }
\end{align*}次数が隣り合う場合、$\mu = \nu \pm 1$として計算する。
例として、$I_{\nu - 1}(x) / I_{\nu}(x)$および$K_{\nu - 1}(x) / K_{\nu}(x)$に対し、$x = 100$としたときの元の関数と近似関数のグラフを示す。近似関数は4次の項まででそれぞれ次式のようになる。
\begin{align}
& \frac{I_{\nu - 1}(x) }{I_{\nu}(x) } \approx 1 + \frac{\nu - 1/ 2}{x} + \frac{4 \nu^{2} - 1}{8 x^{2}} + \frac{4 \nu^{2} - 1}{8 x^{3}} - \frac{16 \nu^{4} - 104 \nu^{2} + 25}{128 x^{4}} \\
& \frac{K_{\nu - 1}(x) }{K_{\nu}(x) } \approx 1 - \frac{\nu - 1/ 2}{x} + \frac{4 \nu^{2} - 1}{8 x^{2}} - \frac{4 \nu^{2} - 1}{8 x^{3}} - \frac{16 \nu^{4} - 104 \nu^{2} + 25}{128 x^{4}}\end{align}


漸近展開の条件より、この近似は$x$が$\nu$に対して十分大きな領域でしか成り立たない。このため、$\nu = x$付近になると誤差が大きくなっていることがわかる。4次の項まで取った場合は広い範囲でよい近似となっているが、$\nu < x / 2$の領域に限れば、2次の項まででも十分精度のよい近似になっている。