Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

第一類貝賽爾函數和第二類貝賽爾函數

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

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

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

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

函數cyl_bessel_j 和函數cyl_neumann 分別返回第一類貝賽爾函數和第二類貝賽爾函數的結果:

cyl_bessel_j(v, x) = Jv(x)

cyl_neumann(v, x) = Yv(x) = Nv(x)

其中:

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

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

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

下面的圖像顯示了函數 Jv的週期性:

下面的圖像顯示了函數Yv的情況: 這同樣也是大數 x, 以及趨近於 -∞ 的小數的週期 :

測試

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

準確性

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

表36.函數 cyl_bessel_j 的出錯率(Errors Rates)

有效數字位數

平台和編譯器

J0 和 J1

Jv

Jv (大數值 x > 1000)

53

Win32 / Visual C++ 8.0

峰值=2.5 均值=1.1

GSL 峰值=6.6

Cephes 峰值=2.5 均值=1.1

峰值=11 均值=2.2

GSL 峰值=11

Cephes 峰值=17 均值=2.5

峰值=413 均值=110

GSL 峰值=6x1011

Cephes 峰值=2x105

64

Red Hat Linux IA64 / G++ 3.4

峰值=7 均值=3

峰值=117 均值=10

峰值=2x104 均值=6x103

64

SUSE Linux AMD64 / G++ 4.1

峰值=7 均值=3

峰值=400 均值=40

峰值=2x104 均值=1x104

113

HP-UX / HP aCC 6

峰值=14 均值=6

峰值=29 均值=3

峰值=2700 均值=450


表7.函數 cyl_neumann 的出錯率

有效數字位數

平台和編譯器

J0 和 J1

Jn (整數次數)

Jv (小數次數)

53

Win32 / Visual C++ 8.0

峰值=330

均值=54

GSL 峰值=34 均值=9

Cephes 峰值=330 均值=54

峰值=923 均值=83

GSL 峰值=500 均值=54

Cephes 峰值=923 均值=83

峰值=561 均值=36

GSL 峰值=1.4x106 均值=7x104

Cephes 峰值=+INF

64

Red Hat Linux IA64 / G++ 3.4

峰值=470 均值=56

峰值=843 均值=51

峰值=741 均值=51

64

SUSE Linux AMD64 / G++ 4.1

峰值=1300 均值=424

峰值=2x104 均值=8x103

峰值=1x105 均值=6x103

113

HP-UX / HP aCC 6

峰值=180 均值=63

峰值=340 均值=150

峰值=2x104 均值=1200


注意 對於大數x 這些函數的精度極大的依賴於函數 std::sin 和函數std::cos.

與 GSL 庫和 Cephes 庫的比較非常有趣: Cephes 庫和這個庫都對整數次數的情況做了優化 - 產生很明顯的結果 - 對大部分處理使用一般的處理情況就可以獲得稍微好的精確性, 就像由GSL庫在整數參數的情況下獲得的更好的準確性所顯示的那樣. 這個庫的實現趨向於在參數變大的時候獲得更好的表現, Cephes 庫在一些測試數據中產生了一些明顯的不精確的結果 (沒有有效數字校正),甚至GSL庫對於函數 Jv的一些輸入的表現也很差。 注意:通過雙重檢查這些結果的方式,Cephes 庫和 GSL 庫的最差運行情況使用functions.wolfram.com進行了重新計算 並且這些結果與我們的測試數據進行了重新比較: 測試數據中沒有發現誤差.

實現

庫的實現大部分是在處理各種特殊情況:

x 負數的時候, 次數v 必須是一個整數否則結果是一個定義域錯誤. 如果次數是一個整數,那麼如果次數是奇數則函數是一個奇函數,如果次數是一個偶數則函數是一個偶函數, 這裡要求滿足x > 0.

當次數 v 是負數的時候,那麼反射方程可以用來調整到v > 0:

注意如果次數是一個整數, 那麼這些方程簡化為:

J-n = (-1)nJn

Y-n = (-1)nYn

然而, 一般來說, 一個負的次數暗示我們必須計算 J 和 Y.

x比次數 v 大,那麼將會使用 M. Abramowitz 和 I.A. Stegun, 數學函數手冊 (Handbook of Mathematical Functions 9.2.19)中的大數x的漸近展開(這比 A&S 9.2.5中的更值得依賴).

當係數 v 是一個整數 這種方法首先依據哪一個更穩定來使用前向循環或後向循環將結果與 J0, J1, Y0 and Y1關聯起來(密勒算法). J0, J1, Y0 和 Y1使用對於小數|x| 的有理數極值逼近以及對於大數|x|的漸近展開來計算係數:

W.J. Cody, ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of Special Function Routines and Test Drivers, ACM Transactions on Mathematical Software, vol 19, 22 (1993).

以及

J.F. Hart et al, Computer Approximations, John Wiley & Sons, New York, 1968.

這種逼近精確於 19 位十進制數字: 因此當T類型的二進制位超過64們的時候就不會使用這種方法.

當x較小的時候, Jx使用下面的級數來計算:

通常情況下我們同時計算 Jv 和 Yv.

為了獲得初始值, 設 μ = ν - floor(ν + 1/2), 那 μ 就是 ν的小數部分滿足 |μ| <= 1/2 (我們需要這一點在稍後滿足收斂性). 方法是計算 Jμ(x), Jμ+1(x), Yμ(x), Yμ+1(x) 然後使用它們來計算 Jν(x), Yν(x).

算法稱作Steed 方法, 需要兩個連續的小數 和 Wronskian:

參見: F.S. Acton, Numerical Methods that Work, The Mathematical Association of America, Washington, 1997.

連分數使用修正的 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可能收斂很慢). Jμ, Jμ+1, Yμ, Yμ+1 可以通過下面的方程計算

其中

Jν 和 Yμ 然後分別使用前向和後向循環來計算 (Miller 算法) .

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

其中:

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

像前面的情況一樣, Yν使用前向遞歸來計算, 因此是 Yν+1. 使用這兩個值以及 fν, Wronskian 不需要使用後向遞歸面直接計算 Jν(x) .


PrevUpHomeNext