| Bessel {base} | R Documentation |
Bessel Functions of integer and fractional order, of first and second kind, J(nu) and Y(nu), and Modified Bessel functions (of first and third kind), I(nu) and K(nu).
gammaCody is the (Γ) function as from the Specfun
package and originally used in the Bessel code.
besselI(x, nu, expon.scaled = FALSE) besselK(x, nu, expon.scaled = FALSE) besselJ(x, nu) besselY(x, nu) gammaCody(x)
x |
numeric, >= 0. |
nu |
numeric; The order (maybe fractional!) of the corresponding Bessel function. |
expon.scaled |
logical; if TRUE, the results are
exponentially scaled in order to avoid overflow
(I(nu)) or underflow (K(nu)),
respectively. |
The underlying C code stems from Netlib (http://www.netlib.org/specfun/r[ijky]besl).
If expon.scaled = TRUE, exp(-x) I(x;nu),
or exp(x) K(x;nu) are returned.
gammaCody may be somewhat faster but less precise and/or robust
than R's standard gamma. It is here for experimental
purpose mainly, and may be defunct very soon.
For nu < 0, formulae 9.1.2 and 9.6.2 from the
reference below are applied (which is probably suboptimal), unless for
besselK which is symmetric in nu.
Numeric vector of the same length of x with the (scaled, if
expon.scale=TRUE) values of the corresponding Bessel function.
Original Fortran code:
W. J. Cody, Argonne National Laboratory
Translation to C and adaption to R:
Martin Maechler maechler@stat.math.ethz.ch.
Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. Dover, New York; Chapter 9: Bessel Functions of Integer Order.
Other special mathematical functions, such as
gamma, Γ(x), and beta,
B(x).
nus <- c(0:5, 10, 20)
x <- seq(0, 4, len = 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, len=801); yl <- c(-.8, .8)
plot(x, x, ylim = yl, ylab = "", type = "n",
main = "Bessel Functions J_nu(x)")
for(nu in nus) lines(x, besselJ(x, nu=nu), col = nu+2)
legend(32,-.18, legend = paste("nu=", nus), col = nus+2, lwd = 1)
## Negative nu's :
xx <- 2:7
nu <- seq(-10, 9, len = 2001)
op <- par(lab = c(16, 5, 7))
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=seq(xx))
par(op)
x0 <- 2^(-20:10)
plot(x0, x0^-8, log="xy", ylab="",type="n",
main = "Bessel Functions J_nu(x) near 0\n log - log scale")
for(nu in sort(c(nus, nus+.5)))
lines(x0, besselJ(x0, nu=nu), col = nu+2)
legend(3, 1e50, legend = paste("nu=", paste(nus, nus+.5, sep=",")),
col = nus + 2, lwd = 1)
plot(x0, x0^-8, log="xy", ylab="", type="n",
main = "Bessel Functions K_nu(x) near 0\n log - log scale")
for(nu in sort(c(nus, nus+.5)))
lines(x0, besselK(x0, nu=nu), col = nu+2)
legend(3, 1e50, legend = paste("nu=", paste(nus, nus+.5, sep=",")),
col = nus + 2, lwd = 1)
x <- x[x > 0]
plot(x, x, ylim=c(1e-18, 1e11), log = "y", ylab = "", type = "n",
main = "Bessel Functions K_nu(x)")
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)