特殊関数: ベッセル・Airy・直交多項式

概要

このページは、ベッセル型関数 (円柱関数)・Airy 関数・古典直交多項式という、互いに密接に関連する 3 つの関数族の数値評価アルゴリズムを解説する。 いずれも 2 階線形常微分方程式の解として現れ、評価には級数展開・漸近展開・三項漸化式という共通の道具立てが使われる。

• ベッセル型関数: ベッセル方程式 $x^2 y'' + x y' + (x^2 - \nu^2) y = 0$ の解。 円柱座標での波動・拡散・ポテンシャル問題に現れる。第 1/2 種 $J_\nu, Y_\nu$、変形 $I_\nu, K_\nu$、球ベッセル、Hankel 関数を含む。
• Airy 関数: $y'' = x y$ の解 $\mathrm{Ai}, \mathrm{Bi}$。 転回点近傍の漸近 (WKB)・光学の焦線・量子力学の一様近似に現れる。次数 $\pm 1/3$ のベッセル関数で表せる。
• 直交多項式: 各々の重み関数に対して直交する多項式列。 Legendre は球面調和の角度部分、Hermite は量子調和振動子、Laguerre は水素原子の動径波動関数に現れる。

関連 API: SpecialFunctions。 実数引数の標準的なベッセル・球ベッセルは C++17 標準ライブラリ (std::cyl_bessel_j 等) に委譲し、 複素引数・任意精度 (Complex<double>, Complex<Float>) は sangi が独自に級数評価する。

ベッセル関数 $J_\nu, Y_\nu, I_\nu, K_\nu$

4 種類の円柱関数を扱う:

• 第 1 種 $J_\nu(x)$ ・第 2 種 (Neumann) $Y_\nu(x)$ ― 通常のベッセル方程式の 2 つの独立解
• 変形第 1 種 $I_\nu(x)$ ・変形第 2 種 $K_\nu(x)$ ― 変形ベッセル方程式 $x^2 y'' + x y' - (x^2 + \nu^2) y = 0$ の解

小引数: 冪級数

原点近傍では Taylor 級数が高速に収束する。第 1 種は

$$J_\nu(z) = \left(\frac{z}{2}\right)^{\nu} \sum_{m=0}^{\infty} \frac{(-1)^m}{m!\,\Gamma(m+\nu+1)} \left(\frac{z}{2}\right)^{2m}$$

変形第 1 種は符号交代がなく

$$I_\nu(z) = \left(\frac{z}{2}\right)^{\nu} \sum_{m=0}^{\infty} \frac{1}{m!\,\Gamma(m+\nu+1)} \left(\frac{z}{2}\right)^{2m}$$

sangi の複素整数次の実装はこの級数を漸化的に評価する。各項を前項から

$$\text{term}_m = \text{term}_{m-1} \cdot \frac{\mp (z/2)^2}{m\,(m+n)}$$

で更新するので、階乗やガンマ関数を毎回計算する必要がない ($J$ は分子に $-(z/2)^2$、$I$ は $+(z/2)^2$)。 相対許容 $|\text{term}_m| < \varepsilon\,|\text{sum}|$ で打ち切る。

整数次の $Y_n, K_n$: 対数項を含む級数

整数次の第 2 種関数は原点で対数発散するため、$\ln(z/2)$ を含む級数 (A&S 9.1.11 / 9.6.11) を用いる。$Y_n$ の構造は

$$Y_n(z) = \frac{2}{\pi} J_n(z)\left(\gamma + \ln\frac{z}{2}\right) - \frac{1}{\pi}\left(\frac{z}{2}\right)^{-n}\!\!\sum_{k=0}^{n-1}\frac{(n-k-1)!}{k!}\left(\frac{z^2}{4}\right)^{k} - \frac{1}{\pi}\left(\frac{z}{2}\right)^{n}\!\!\sum_{k=0}^{\infty}\frac{(-1)^k (H_k + H_{k+n})}{k!\,(k+n)!}\left(\frac{z^2}{4}\right)^{k}$$

ここで $\gamma$ は Euler-Mascheroni 定数、$H_k = \sum_{j=1}^{k} 1/j$ は調和数。 $K_n$ も同様に $\ln(z/2)$ と $(H_k + H_{k+n} - 2\gamma)$ を含む形になる (符号と係数が異なる)。 sangi では調和数とガンマ項を漸化的に積み上げ、上の $J_n$ / $I_n$ を主要項として再利用する。

大引数: 漸近展開

引数が大きい領域では振動・減衰の漸近形が使われる。第 1 種は

$$J_\nu(x) \sim \sqrt{\frac{2}{\pi x}}\,\cos\!\left(x - \frac{\nu\pi}{2} - \frac{\pi}{4}\right), \qquad x \to \infty$$

これは発散級数なので、項の絶対値が最小になる手前で打ち切る (最適打ち切り)。 sangi の実数次ベッセルは標準ライブラリに委譲しており、これらの領域分割は実装内部で扱われる。

後退漸化 (Miller のアルゴリズム)

整数次まで一括して $J_n$ を得たいとき、三項漸化式

$$J_{\nu+1}(x) = \frac{2\nu}{x} J_\nu(x) - J_{\nu-1}(x)$$

を次数の増える向き (前進) に回すと不安定になる。$J_n$ と $Y_n$ は同じ漸化式を満たし、前進方向では $Y_n$ 成分 (これは次数とともに指数増大する) が丸め誤差から励起されて解を覆い隠すからである。

安定なのは後退方向。Miller のアルゴリズムは十分高い次数 $N \gg n$ から $J_{N+1} = 0,\ J_N = \varepsilon$ と仮定して

$$J_{\nu-1}(x) = \frac{2\nu}{x} J_\nu(x) - J_{\nu+1}(x)$$

を $\nu = N, N-1, \ldots$ と下げながら回し、最後に総和規則 $J_0(x) + 2\sum_{k\geq 1} J_{2k}(x) = 1$ で全体を正規化する。 誤差成分は後退方向では減衰するため、初期値が任意でも比 $J_n / J_0$ は正しく得られる。

Wronskian による正規化と導関数

第 1 種・第 2 種の Wronskian は

$$W\{J_\nu, Y_\nu\}(x) = J_\nu(x) Y_\nu'(x) - J_\nu'(x) Y_\nu(x) = \frac{2}{\pi x}$$

で、これは独立解の規格化と相互整合性の検証に使える。導関数は漸化式 (DLMF 10.6.1, 10.29.1 等) で隣接次数から求める:

$$J_\nu'(x) = \tfrac{1}{2}\big(J_{\nu-1}(x) - J_{\nu+1}(x)\big), \qquad I_\nu'(x) = \tfrac{1}{2}\big(I_{\nu-1}(x) + I_{\nu+1}(x)\big)$$

$Y_\nu', K_\nu'$ も同様 ($K$ は $K_\nu'(x) = -\tfrac12(K_{\nu-1}+K_{\nu+1})$)。 sangi の besselJPrime, besselIPrime 等はこの漸化式をそのまま実装する。

球ベッセルと Hankel 関数

球ベッセル関数は半整数次のベッセルで、初等関数で閉じる:

$$j_n(x) = \sqrt{\frac{\pi}{2x}}\, J_{n+1/2}(x), \qquad y_n(x) = \sqrt{\frac{\pi}{2x}}\, Y_{n+1/2}(x)$$

$j_0(x) = \sin x / x$、$y_0(x) = -\cos x / x$ から出発し、漸化式 (DLMF 10.51.1)

$$j_{n+1}(x) = \frac{2n+1}{x} j_n(x) - j_{n-1}(x)$$

を前進方向に適用する (小次数・適度な引数では実用上安定)。 Hankel 関数 (第 3 種) は進行波を表す複素結合で、第 1/2 種を組み合わせて作る:

$$H_n^{(1)}(z) = J_n(z) + i\,Y_n(z), \qquad H_n^{(2)}(z) = J_n(z) - i\,Y_n(z)$$

球 Hankel 関数 $h_n^{(1,2)} = j_n \pm i\,y_n$ も同様に作られる。 sangi では hankelH1, hankelH2, sphericalHankelH1 等として提供する。

Airy 関数 $\mathrm{Ai}, \mathrm{Bi}$

Airy 関数は方程式 $y'' = x y$ の 2 つの独立解。$\mathrm{Ai}(x)$ は $x \to +\infty$ で減衰、$\mathrm{Bi}(x)$ は発散する。

原点付近の Maclaurin 級数

Airy 方程式の 2 つの整級数解を取る:

$$f(x) = 1 + \frac{x^3}{2\cdot 3} + \frac{x^6}{2\cdot 3\cdot 5\cdot 6} + \cdots, \qquad g(x) = x + \frac{x^4}{3\cdot 4} + \frac{x^7}{3\cdot 4\cdot 6\cdot 7} + \cdots$$

これらは超幾何関数 $_0F_1$ で表せる ($f = {}_0F_1(;\tfrac23;\tfrac{x^3}{9})$ など) 整関数で、漸化式

$$f_{\text{term}_{k+1}} = f_{\text{term}_{k}} \cdot \frac{x^3}{(3k+2)(3k+3)}, \qquad g_{\text{term}_{k+1}} = g_{\text{term}_{k}} \cdot \frac{x^3}{(3k+3)(3k+4)}$$

で各項を更新する。Airy 関数はこの 2 解の線形結合:

$$\mathrm{Ai}(x) = \mathrm{Ai}(0)\, f(x) + \mathrm{Ai}'(0)\, g(x), \qquad \mathrm{Bi}(x) = \mathrm{Bi}(0)\, f(x) + \mathrm{Bi}'(0)\, g(x)$$

結合係数はガンマ関数で与えられる:

$$\mathrm{Ai}(0) = \frac{1}{3^{2/3}\,\Gamma(2/3)}, \quad \mathrm{Ai}'(0) = -\frac{1}{3^{1/3}\,\Gamma(1/3)}, \quad \mathrm{Bi}(0) = \frac{1}{3^{1/6}\,\Gamma(2/3)}, \quad \mathrm{Bi}'(0) = \frac{3^{1/6}}{\Gamma(1/3)}$$

$f, g$ は整関数なので級数は全平面で収束し、複素引数にもそのまま適用できる。 導関数 $\mathrm{Ai}', \mathrm{Bi}'$ も $f', g'$ を同じ枠組みで級数評価する。 任意精度版では $\Gamma(1/3), \Gamma(2/3)$ と $3^{1/3}, 3^{1/6}$ を作業精度で計算してから結合係数を作る。

ベッセル関数による表現

Airy 関数は次数 $\pm 1/3$ のベッセル関数で表せ、両者の関係を明示する:

$$\mathrm{Ai}(x) = \frac{1}{\pi}\sqrt{\frac{x}{3}}\, K_{1/3}\!\left(\tfrac{2}{3}x^{3/2}\right) \quad (x>0)$$

$$\mathrm{Ai}(-x) = \frac{\sqrt{x}}{3}\left[J_{1/3}\!\left(\tfrac{2}{3}x^{3/2}\right) + J_{-1/3}\!\left(\tfrac{2}{3}x^{3/2}\right)\right] \quad (x>0)$$

これは Airy 関数とベッセル関数族が同じ円柱関数の枠組みに属することを示すと同時に、 大引数での評価に変形・通常ベッセルの漸近を流用する道を与える。

大 $|x|$ の漸近

$\zeta = \tfrac{2}{3}|x|^{3/2}$ として、原点から離れた領域では漸近形が使われる:

$$\mathrm{Ai}(x) \sim \frac{e^{-\zeta}}{2\sqrt{\pi}\, x^{1/4}} \quad (x\to +\infty), \qquad \mathrm{Ai}(-x) \sim \frac{1}{\sqrt{\pi}\, x^{1/4}}\sin\!\left(\zeta + \tfrac{\pi}{4}\right) \quad (x\to +\infty)$$

正の側では指数減衰、負の側では振動するという Airy 関数特有の挙動 (転回点を境にした遷移) がここに表れる。

Kelvin・Struve・Anger-Weber

ベッセル関数から派生する、あるいは隣接する関数族をまとめて扱う。

Kelvin 関数 ber, bei, ker, kei

Kelvin 関数は複素引数のベッセル関数の実部・虚部として定義される (次数 0):

$$\mathrm{ber}(x) + i\,\mathrm{bei}(x) = J_0\!\left(x\,e^{3i\pi/4}\right), \qquad \mathrm{ker}(x) + i\,\mathrm{kei}(x) = K_0\!\left(x\,e^{i\pi/4}\right)$$

sangi の実装は、複素引数の $J_0$ / $K_0$ をそのまま呼び出して実部・虚部を取り出す。 $\mathrm{ber}, \mathrm{bei}$ は偶関数で全実数で定義され、$\mathrm{ker}, \mathrm{kei}$ は $x>0$ でのみ定義される ($\mathrm{ker}$ は $x\to 0^+$ で対数発散、$\mathrm{kei}(0) = -\pi/4$)。 電気工学の表皮効果やソレノイドコイル中の磁性コアの磁場計算に現れる。

Struve 関数 $H_\nu, L_\nu$

Struve 関数は非斉次ベッセル方程式の特殊解で、小引数では級数 (DLMF 11.2.1):

$$H_\nu(z) = \left(\frac{z}{2}\right)^{\nu+1}\sum_{k=0}^{\infty} \frac{(-1)^k}{\Gamma(k+\tfrac32)\,\Gamma(\nu+k+\tfrac32)}\left(\frac{z}{2}\right)^{2k}$$

変形 Struve $L_\nu$ は符号交代を除いた同形。大引数では漸近展開 (DLMF 11.6.1)

$$H_\nu(z) \sim Y_\nu(z) + \frac{1}{\pi}\sum_{k=0}^{N} \frac{\Gamma(k+\tfrac12)}{\Gamma(\nu+\tfrac12-k)}\left(\frac{z}{2}\right)^{\nu-2k-1}$$

を使い ($L_\nu$ は $I_\nu$ を主要項に持ち符号が異なる)、漸近級数は項の絶対値が増え始めた時点で打ち切る。 sangi では $|z| \lesssim 16$ を境に級数と漸近を切り替える。現状は実数引数・$z \geq 0$ のみ対応。

Anger 関数 $J_\nu$ ・ Weber 関数 $E_\nu$

Anger・Weber 関数は積分で定義される (DLMF 11.10):

$$\mathbf{J}_\nu(z) = \frac{1}{\pi}\int_0^{\pi}\cos(\nu t - z\sin t)\,dt, \qquad \mathbf{E}_\nu(z) = \frac{1}{\pi}\int_0^{\pi}\sin(\nu t - z\sin t)\,dt$$

次数 $\nu$ が整数 $n$ のとき Anger 関数はベッセル関数 $J_n$ に一致するが、非整数次では独立な別関数族になる。 被積分関数は滑らかな解析関数なので、sangi は定義積分を合成 Simpson 則で直接数値積分する (級数表現は整数次極限の扱いが煩雑なため採用しない)。

直交多項式

古典直交多項式は、各々の重み関数 $w(x)$ に対して区間上で直交する多項式列。 いずれも三項漸化式で安定に評価できる (これは Favard の定理の帰結)。

Legendre 多項式 $P_n(x)$

重み $w(x) = 1$、区間 $[-1, 1]$。Bonnet の漸化式で評価する:

$$(n+1) P_{n+1}(x) = (2n+1)\, x\, P_n(x) - n\, P_{n-1}(x), \qquad P_0 = 1,\ P_1 = x$$

Chebyshev 多項式 $T_n(x), U_n(x)$

第 1 種 $T_n$ (重み $1/\sqrt{1-x^2}$) と第 2 種 $U_n$ (重み $\sqrt{1-x^2}$) は同じ漸化式を共有し、初期値だけが異なる:

$$T_{n+1}(x) = 2x\, T_n(x) - T_{n-1}(x), \quad T_0=1,\ T_1=x$$

$$U_{n+1}(x) = 2x\, U_n(x) - U_{n-1}(x), \quad U_0=1,\ U_1=2x$$

三角関数表現 $T_n(\cos\theta) = \cos(n\theta)$、$U_n(\cos\theta) = \sin((n+1)\theta)/\sin\theta$ も成り立つ。 $|x| \leq 1$ で安定、$|x| > 1$ でも漸化式はそのまま使える。

Hermite 多項式 (物理流/確率流)

物理流 $H_n$ (重み $e^{-x^2}$、量子調和振動子の固有関数) と確率流 $He_n$ (重み $e^{-x^2/2}$) で漸化式が異なる:

$$H_{n+1}(x) = 2x\, H_n(x) - 2n\, H_{n-1}(x), \quad H_0=1,\ H_1=2x$$

$$He_{n+1}(x) = x\, He_n(x) - n\, He_{n-1}(x), \quad He_0=1,\ He_1=x$$

両者は $H_n(x) = 2^{n/2} He_n(x\sqrt{2})$ で結ばれる。

Laguerre 多項式 $L_n(x), L_n^{\alpha}(x)$

重み $e^{-x}$ (一般化版は $x^\alpha e^{-x}$)、区間 $[0, \infty)$。水素原子の動径波動関数に現れる。

$$(k+1) L_{k+1}^{\alpha}(x) = (2k+1+\alpha-x)\, L_k^{\alpha}(x) - (k+\alpha)\, L_{k-1}^{\alpha}(x), \quad L_0^{\alpha}=1,\ L_1^{\alpha}=1+\alpha-x$$

$\alpha = 0$ で通常の Laguerre $L_n(x)$ に帰着する。

Jacobi 多項式 $P_n^{(\alpha,\beta)}(x)$ と Gegenbauer $C_n^{\lambda}(x)$

Jacobi 多項式は重み $(1-x)^{\alpha}(1+x)^{\beta}$ ($\alpha, \beta > -1$) に対する直交多項式で、上記の多くを統一的に含む:

$$P_n^{(0,0)} = P_n \ (\text{Legendre}), \quad P_n^{(-1/2,-1/2)} \propto T_n, \quad P_n^{(1/2,1/2)} \propto U_n, \quad P_n^{(\alpha,\alpha)} \propto C_n^{\alpha+1/2}$$

評価は DLMF 18.9.2 の三項漸化式 (係数が $\alpha, \beta, n$ に依存する) を回す。 Gegenbauer (超球) 多項式 $C_n^{\lambda}$ は対称な Jacobi に比例し、専用の漸化式

$$n\, C_n^{\lambda}(x) = 2(n+\lambda-1)\, x\, C_{n-1}^{\lambda}(x) - (n+2\lambda-2)\, C_{n-2}^{\lambda}(x)$$

を持つ。sangi では Chebyshev・Hermite・Laguerre を専用実装 (Jacobi 経由より低オーバーヘッド) として提供し、一般の $\alpha, \beta$ には jacobiP を用いる。

Clenshaw アルゴリズム

直交多項式の級数 $S(x) = \sum_{k=0}^{N} c_k\, p_k(x)$ を一括評価するとき、各 $p_k(x)$ を作って掛け足すよりClenshaw アルゴリズムが効率的かつ安定。 漸化式 $p_{k+1} = (\alpha_k x + \beta_k) p_k - \gamma_k p_{k-1}$ に対し、補助列 $b_k$ を後退方向に

$$b_k = c_k + (\alpha_k x + \beta_k)\, b_{k+1} - \gamma_{k+1}\, b_{k+2}, \qquad b_{N+1} = b_{N+2} = 0$$

と計算すると、$S(x)$ は $b_0, b_1$ と低次の $p_0, p_1$ だけで表せる。 Chebyshev 級数や球面調和展開の評価で標準的に使われる手法。

随伴 Legendre と球面調和

随伴 Legendre 関数 $P_n^m(x)$ は球面調和の角度部分に現れ、対角項から漸化式で立ち上げる:

$$P_m^m(x) = (2m-1)!!\,(1-x^2)^{m/2}, \quad P_{m+1}^m(x) = x(2m+1) P_m^m(x)$$

$$(n-m+1) P_{n+1}^m(x) = (2n+1)\, x\, P_n^m(x) - (n+m)\, P_{n-1}^m(x)$$

正規化球面調和 (実部) は

$$Y_n^m(\theta) = \sqrt{\frac{2n+1}{4\pi}\cdot\frac{(n-m)!}{(n+m)!}}\; P_n^m(\cos\theta)$$

で与えられる。sangi の sphLegendre は階乗比 $(n-m)!/(n+m)!$ をオーバーフロー回避のため対数空間で計算してから正規化係数を作る。 この実装は Condon-Shortley 位相を含まない (C++17 標準ライブラリ std::assoc_legendre と整合)。

比較表

共通する設計指針は、級数は漸化的に各項を更新して階乗・ガンマの再計算を避けること、 漸化式は安定な向きに回すこと (ベッセルは後退、直交多項式は前進)、 そして領域で手法を切り替える (小引数=級数、大引数=漸近) ことである。

参考文献

• Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards. (§9 ベッセル関数, §10 Airy 関数, §12 Struve・Anger-Weber, §22 直交多項式)
• Olver, F. W. J., Lozier, D. W., Boisvert, R. F. & Clark, C. W. (eds.) (2010). NIST Handbook of Mathematical Functions. Cambridge University Press. (DLMF, dlmf.nist.gov §10, §11, §18)
• Olver, F. W. J. (1974). Asymptotics and Special Functions. Academic Press.
• Watson, G. N. (1944). A Treatise on the Theory of Bessel Functions. 2nd ed. Cambridge University Press.
• Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing. 3rd ed. Cambridge University Press. (§6 特殊関数, Miller のアルゴリズム, Clenshaw 漸化)
• Gil, A., Segura, J. & Temme, N. M. (2007). Numerical Methods for Special Functions. SIAM.