概要

template <class T1, class T2>
calculated-result-type cyl_bessel_i(T1 v, T2 x);

template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_i(T1 v, T2 x, const Policy&);

template <class T1, class T2>
calculated-result-type cyl_bessel_k(T1 v, T2 x);

template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_k(T1 v, T2 x, const Policy&);

說明

函數cyl_bessel_i 和 函數cyl_bessel_k 分別返回第一類修正貝賽爾函數和第二類修正貝賽爾函數:

cyl_bessel_i(v, x) = I_v(x)

cyl_bessel_k(v, x) = K_v(x)

其中:

當T1和T2是不同類型的時候，函數的返回類型使用返回值類型推導法則 來確定. 這些函數也針對T1是整數的情況進行了優化.

最後一個策略 參數是可選 的並且可以用來控制函數的行為: 如何處理錯誤, 使用哪種層次的精度等等. 參見策略文檔來瞭解更多信息。

當結果是未定義或是複數的時候函數返回定義域錯誤. 對於函數cyl_bessel_j 這會在x < 0 且 v不是一個整數的情況下發生, 或者當x == 0 且 v != 0時發生. 對於函數cyl_neumann 這會在x <= 0時發生。

下面的圖像顯示了 I_v 的指數特性：

下面的圖像顯示了 K_v.的指數下降特性：

測試

有兩種測試數據: 使用functions.wolfram.com計算出的抽樣測試數據以及使用這種實現的簡化版本來計算出的數據量更大的測試數據 (去掉所有的對特殊情況的處理).

精確性

下面的表顯示了這些函數在不同平台上的準確性的變化， 以及與 GSL-1.9 和Cephes 庫的比較. 注意: 這些函數的週期性意味著它們具有無限個數的無理根: 一般來說當參數足夠接近根的時候這些函數具在任意大的相對誤差. 當然在這種情況下絕對誤差永遠都很小. 注意：在那些將寬浮點類型給定為窄浮點類型的系統上，函數的返回值具有有效的零誤差. 所有結果的相對誤差都是10的-5次方數量級。

實現

下面的情況首先作為特殊情況進行處理:

當 x < 0時計算I_v , ν 必須是一個整數 否則會產生一個定義域錯誤. 如果ν 是一個整數, 那麼當 ν 是奇數的時候函數是奇函數，當ν 是偶數的時候函數是偶函數, 這裡要求符合x > 0.

對於 I_v ，當 v等於 0, 1 或 0.5 時，作為特殊情況處理。

v = 0 和 v = 1的情況使用在有限和無限區間上的有理數極值逼近的方法. 係數來自於:

而 v = 0.5 的情況是一個簡單的三角函數:

I_0.5(x) = sqrt(2 / πx) * sinh(x)

對於 K_v ，當 v 是一個整數時, 結果使用遞歸關係來計算:

從 K₀ 和 K₁ 開始使用上面的有理數逼近來計算. 這些有理數逼近的精度為 19 個數字, 因此僅當 T 的精度為不超過64個有效數字時才使用這種方法.

通常情況下, 我們首先使用反射方程將ν正規化到區間 [0, [inf]) :

設μ = ν - floor(ν + 1/2), 那麼μ 就是 ν 的小數部分且滿足 |μ| <= 1/2 (我們需要這一點用於後面的收斂性). 這個思路是計算K_μ(x) 和 K_μ+1(x), 然後用它們來計算I_ν(x) 和 K_ν(x).

這個算法求Temme 在 N.M. Temme, On the numerical evaluation of the modified bessel function of the third kind, Journal of Computational Physics, vol 19, 324 (1975)中發表, 需要兩個連分數以及Wronskian:

連分數使用修正的 Lentz's 算法來獲得 (W.J. Lentz, Generating Bessel functions in Mie scattering calculations using continued fractions, Applied Optics, vol 15, 668 (1976)). 收斂半徑依賴於x,因此對於大數x 和小數x. 我們需要使用不同的策略

x > v, CF1 需要 O(x)步來收斂, CF2 快速收斂

x <= v, CF1 快速收斂, 當 x -> 0 時，CF2不能收斂

當x 較大時(x > 2), 連續分數收斂 (對於 大數 x ,CF1可能收斂很慢). K_μ 和 K_μ+1 可以通過下面的方程計算

其中

S是一個與 CF2同時求和的級數, 參見 I.J. Thompson and A.R. Barnett, Modified Bessel functions I_v and K_v of real order and complex argument to selected accuracy, Computer Physics Communications, vol 47, 245 (1987).

當x 較小 (x <= 2), CF2 可能無法收斂 (但 CF1 可以收斂). 解為 Temme 級數:

其中

f_k 和 h_k 也使用遞歸來計算 (包括 γ 函數), 但是方程稍有些複雜, 讀者可以參見 N.M. Temme, On the numerical evaluation of the modified Bessel function of the third kind, Journal of Computational Physics, vol 19, 324 (1975). 注意 Temme 級數 僅對 |μ| <= 1/2收斂.

K_ν(x) 使用前向遞歸來計算, K_ν+1(x)也是這樣的.使用這兩個值和 f_ν, Wronskian 直接生成 I_ν(x) .