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).
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.
Details
If expon.scaled = TRUE, exp(-x) I(x;nu),
or exp(x) K(x;nu) are returned.
For nu < 0, formulae 9.1.2 and 9.6.2 from Abramowitz &
Stegun are applied (which is probably suboptimal), except for
besselK which is symmetric in nu.
The current algorithms will give warnings about accuracy loss for
large arguments. In some cases, these warnings are exaggerated, and
the precision is perfect. For large nu, say in the order of
millions, the current algorithms are rarely useful.
Value
Numeric vector with the (scaled, if expon.scaled = TRUE)
values of the corresponding Bessel function.
The length of the result is the maximum of the lengths of the
parameters. All parameters are recycled to that length.
Author(s)
Original Fortran code:
W. J. Cody, Argonne National Laboratory
Translation to C and adaption to R:
Martin Maechler maechler@stat.math.ethz.ch.
Source
The C code is a translation of Fortran routines from
http://www.netlib.org/specfun/ribesl, ../rjbesl, etc.
The four source code files for bessel[IJKY] each contain a paragraph
“Acknowledgement” and “References”, a short summary of
which is
besselI
based on (code) by David J. Sookne, see Sookne (1973)...
Modifications... An earlier version was published in Cody (1983).
besselJ
as besselI
besselK
based on (code) by J. B. Campbell (1980)... Modifications...
besselY
draws heavily on Temme's Algol program for
Y... and on Campbell's programs for Y_ν(x)
.... ... heavily modified.
References
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:
Sockne, David J. (1973)
Bessel Functions of Real Argument and Integer Order.
NBS Jour. of Res. B.77B, 125–132.
Cody, William J. (1983)
Algorithm 597: Sequence of modified Bessel functions of the first kind.
ACM Transactions on Mathematical Software9(2), 242–245.
Campbell, J.B. (1980)
On Temme's algorithm for the modified Bessel function of the third kind.
ACM Transactions on Mathematical Software6(4), 581–586.
Campbell, J.B. (1979)
Bessel functions J_nu(x) and Y_nu(x) of float order and float argument.
Comp. Phy. Comm.18, 133–142.
Temme, Nico M. (1976)
On the numerical evaluation of the ordinary Bessel function of the second kind.
J. Comput. Phys.21, 343–350.
See Also
Other special mathematical functions, such as
gamma, Γ(x), and beta,
B(x).
Examples
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(-.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, length.out = 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+0.5)))
lines(x0, besselJ(x0, nu = nu), col = nu + 2)
legend(3, 1e50, legend = paste("nu=", paste(nus, nus+0.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+0.5)))
lines(x0, besselK(x0, nu = nu), col = nu + 2)
legend(3, 1e50, legend = paste("nu=", paste(nus, nus + 0.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)
## 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)}))
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(base)
> png(filename="/home/ddbj/snapshot/RGM3/R_rel/result/base/Bessel.Rd_%03d_medium.png", width=480, height=480)
> ### Name: Bessel
> ### Title: Bessel Functions
> ### Aliases: bessel Bessel besselI besselJ besselK besselY
> ### Keywords: math
>
> ### ** Examples
>
> 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(-.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, length.out = 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+0.5)))
+ lines(x0, besselJ(x0, nu = nu), col = nu + 2)
> legend(3, 1e50, legend = paste("nu=", paste(nus, nus+0.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+0.5)))
+ lines(x0, besselK(x0, nu = nu), col = nu + 2)
> legend(3, 1e50, legend = paste("nu=", paste(nus, nus + 0.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)
>
> ## 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)}))
>
>
>
>
>
> dev.off()
null device
1
>