R Bessel 貝塞爾函數

說明

第一類和第二類整數階和小數階貝塞爾函數 J_{\nu} 和 Y_{\nu} 以及修正貝塞爾函數(第一類和第三類)I_{\nu} 和 K_{\nu} 。

用法

besselI(x, nu, expon.scaled = FALSE)
besselK(x, nu, expon.scaled = FALSE)
besselJ(x, nu)
besselY(x, nu)

參數

x	數字，\ge 0 。
nu	數字;相應貝塞爾函數的階數(可能是分數和負數)。
expon.scaled	邏輯性；如果是 TRUE ，則結果按指數縮放，以避免溢出 ( I_{\nu} ) 或下溢 ( K_{\nu} )。

細節

如果返回 expon.scaled = TRUE 、 e^{-x} I_{\nu}(x) 或 e^{x} K_{\nu}(x) 。

對於 \nu < 0 ，應用 Abramowitz & Stegun 的公式 9.1.2 和 9.6.2(這可能不是最佳的)，但 besselK 除外，它在 nu 中是對稱的。

當前的算法將針對大型參數的準確性損失發出警告。在某些情況下，這些警告被誇大了，而精確度卻是完美的。對於大型 nu ，比如數百萬個，當前的算法很少有用。

值

具有相應貝塞爾函數的(縮放，如果 expon.scaled = TRUE )值的數值向量。

結果的長度是參數長度中的最大值。所有參數都會循環到該長度。

例子

require(graphics)

nus <- c(0:5, 10, 20)

x <- seq(0, 4, length.out = 501)
plot(x, x, ylim = c(0, 6), ylab = "", type = "n",
     main = "Bessel Functions  I_nu(x)")
for(nu in nus) lines(x, besselI(x, nu = nu), col = nu + 2)
legend(0, 6, legend = paste("nu=", nus), col = nus + 2, lwd = 1)

x <- seq(0, 40, length.out = 801); yl <- c(-.5, 1)
plot(x, x, ylim = yl, ylab = "", type = "n",
     main = "Bessel Functions  J_nu(x)")
abline(h=0, v=0, lty=3)
for(nu in nus) lines(x, besselJ(x, nu = nu), col = nu + 2)
legend("topright", legend = paste("nu=", nus), col = nus + 2, lwd = 1, bty="n")

## Negative nu's --------------------------------------------------
xx <- 2:7
nu <- seq(-10, 9, length.out = 2001)
## --- I() --- --- --- ---
matplot(nu, t(outer(xx, nu, besselI)), type = "l", ylim = c(-50, 200),
        main = expression(paste("Bessel ", I[nu](x), " for fixed ", x,
                                ",  as ", f(nu))),
        xlab = expression(nu))
abline(v = 0, col = "light gray", lty = 3)
legend(5, 200, legend = paste("x=", xx), col=seq(xx), lty=1:5)

## --- J() --- --- --- ---
bJ <- t(outer(xx, nu, besselJ))
matplot(nu, bJ, type = "l", ylim = c(-500, 200),
        xlab = quote(nu), ylab = quote(J[nu](x)),
        main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
abline(v = 0, col = "light gray", lty = 3)
legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)

## ZOOM into right part:
matplot(nu[nu > -2], bJ[nu > -2,], type = "l",
        xlab = quote(nu), ylab = quote(J[nu](x)),
        main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
abline(h=0, v = 0, col = "gray60", lty = 3)
legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)


##---------------  x --> 0  -----------------------------
x0 <- 2^seq(-16, 5, length.out=256)
plot(range(x0), c(1e-40, 1), log = "xy", xlab = "x", ylab = "", type = "n",
     main = "Bessel Functions  J_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
for(nu in sort(c(nus, nus+0.5)))
    lines(x0, besselJ(x0, nu = nu), col = nu + 2, lty= 1+ (nu%%1 > 0))
legend("right", legend = paste("nu=", paste(nus, nus+0.5, sep=", ")),
       col = nus + 2, lwd = 1, bty="n")

x0 <- 2^seq(-10, 8, length.out=256)
plot(range(x0), 10^c(-100, 80), log = "xy", xlab = "x", ylab = "", type = "n",
     main = "Bessel Functions  K_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
for(nu in sort(c(nus, nus+0.5)))
    lines(x0, besselK(x0, nu = nu), col = nu + 2, lty= 1+ (nu%%1 > 0))
legend("topright", legend = paste("nu=", paste(nus, nus + 0.5, sep = ", ")),
       col = nus + 2, lwd = 1, bty="n")

x <- x[x > 0]
plot(x, x, ylim = c(1e-18, 1e11), log = "y", ylab = "", type = "n",
     main = "Bessel Functions  K_nu(x)"); axis(2, at=1)
for(nu in nus) lines(x, besselK(x, nu = nu), col = nu + 2)
legend(0, 1e-5, legend=paste("nu=", nus), col = nus + 2, lwd = 1)

yl <- c(-1.6, .6)
plot(x, x, ylim = yl, ylab = "", type = "n",
     main = "Bessel Functions  Y_nu(x)")
for(nu in nus){
    xx <- x[x > .6*nu]
    lines(xx, besselY(xx, nu=nu), col = nu+2)
}
legend(25, -.5, legend = paste("nu=", nus), col = nus+2, lwd = 1)

## negative nu in bessel_Y -- was bogus for a long time
curve(besselY(x, -0.1), 0, 10, ylim = c(-3,1), ylab = "")
for(nu in c(seq(-0.2, -2, by = -0.1)))
  curve(besselY(x, nu), add = TRUE)
title(expression(besselY(x, nu) * "   " *
                 {nu == list(-0.1, -0.2, ..., -2)}))

作者

Original Fortran code: W. J. Cody, Argonne National Laboratory
Translation to C and adaptation to R: Martin Maechler maechler@stat.math.ethz.ch.

來源

C 代碼是 Fortran 例程的翻譯https://netlib.org/specfun/ribesl, '⁠../rjbesl⁠bessel[IJKY] 的四個源代碼文件每個都包含一段 “Acknowledgement” 和 “References”，其簡短摘要為

貝塞爾

基於 David J. Sookne 的(代碼)，請參閱 Sookne (1973)... 修改... 早期版本發表於 Cody (1983)。

貝塞爾

如besselI

貝塞爾K

基於 J. B. Campbell (1980) 的(代碼)...修改...

貝塞爾

大量借鑒了 Temme 的 Algol 程序 Y ... 以及 Campbell 的 Y_\nu(x) 程序 .... ... 進行了大量修改。

參考

Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions. Dover, New York; Chapter 9: Bessel Functions of Integer Order.

In order of “Source” citation above:

Sookne, David J. (1973). Bessel Functions of Real Argument and Integer Order. Journal of Research of the National Bureau of Standards, 77B, 125-132. doi:10.6028/jres.077B.012.

Cody, William J. (1983). Algorithm 597: Sequence of modified Bessel functions of the first kind. ACM Transactions on Mathematical Software, 9(2), 242-245. doi:10.1145/357456.357462.

Campbell, J.B. (1980). On Temme's algorithm for the modified Bessel function of the third kind. ACM Transactions on Mathematical Software, 6(4), 581-586. doi:10.1145/355921.355928.

Campbell, J.B. (1979). Bessel functions J_nu(x) and Y_nu(x) of float order and float argument. Computer Physics Communications, 18, 133-142. doi:10.1016/0010-4655(79)90030-4.

Temme, Nico M. (1976). On the numerical evaluation of the ordinary Bessel function of the second kind. Journal of Computational Physics, 21, 343-350. doi:10.1016/0021-9991(76)90032-2.

也可以看看

其他特殊數學函數，例如 gamma 、 \Gamma(x) 和 beta 、 B(x) 。