1% File src/library/base/man/Bessel.Rd
2% Part of the R package, https://www.R-project.org
3% Copyright 1995-2018 R Core Team
4% Distributed under GPL 2 or later
5
6\name{Bessel}
7\title{Bessel Functions}
8\alias{bessel}
9\alias{Bessel}
10\alias{besselI}
11\alias{besselJ}
12\alias{besselK}
13\alias{besselY}
14\usage{
15besselI(x, nu, expon.scaled = FALSE)
16besselK(x, nu, expon.scaled = FALSE)
17besselJ(x, nu)
18besselY(x, nu)
19}
20\description{
21  Bessel Functions of integer and fractional order, of first
22  and second kind, \eqn{J_{\nu}}{J(nu)} and \eqn{Y_{\nu}}{Y(nu)}, and
23  Modified Bessel functions (of first and third kind),
24  \eqn{I_{\nu}}{I(nu)} and \eqn{K_{\nu}}{K(nu)}.
25}
26\arguments{
27  \item{x}{numeric, \eqn{\ge 0}.}
28
29  \item{nu}{numeric; The \emph{order} (maybe fractional and negative) of
30    the corresponding Bessel function.}
31
32  \item{expon.scaled}{logical; if \code{TRUE}, the results are
33    exponentially scaled in order to avoid overflow
34    (\eqn{I_{\nu}}{I(nu)}) or underflow (\eqn{K_{\nu}}{K(nu)}),
35    respectively.}
36}
37\value{
38  Numeric vector with the (scaled, if \code{expon.scaled = TRUE})
39  values of the corresponding Bessel function.
40
41  The length of the result is the maximum of the lengths of the
42  parameters.  All parameters are recycled to that length.
43}
44\details{
45  If \code{expon.scaled = TRUE}, \eqn{e^{-x} I_{\nu}(x)}{exp(-x) I(x;nu)},
46  or \eqn{e^{x} K_{\nu}(x)}{exp(x) K(x;nu)} are returned.
47
48  For \eqn{\nu < 0}{nu < 0}, formulae 9.1.2 and 9.6.2 from Abramowitz &
49  Stegun  are applied (which is probably suboptimal), except for
50  \code{besselK} which is symmetric in \code{nu}.
51
52  The current algorithms will give warnings about accuracy loss for
53  large arguments.  In some cases, these warnings are exaggerated, and
54  the precision is perfect.  For large \code{nu}, say in the order of
55  millions, the current algorithms are rarely useful.
56}
57\source{
58  The C code is a translation of Fortran routines from
59  \url{https://www.netlib.org/specfun/ribesl}, \samp{../rjbesl}, etc.
60  The four source code files for bessel[IJKY] each contain a paragraph
61  \dQuote{Acknowledgement} and \dQuote{References}, a short summary of
62  which is
63  \describe{
64    \item{besselI}{based on (code) by David J. Sookne, see Sookne (1973)\dots
65      Modifications\dots An earlier version was published in Cody (1983).}
66    \item{besselJ}{as \code{besselI}}
67    \item{besselK}{based on (code) by J. B. Campbell (1980)\dots Modifications\dots}
68    \item{besselY}{draws heavily on Temme's Algol program for
69      \eqn{Y}\dots and on Campbell's programs for \eqn{Y_\nu(x)}
70      \dots. \dots heavily modified.}
71  }
72}
73\references{
74  Abramowitz, M. and Stegun, I. A. (1972).
75  \emph{Handbook of Mathematical Functions}.
76  Dover, New York;
77  Chapter 9: Bessel Functions of Integer Order.
78
79  In order of \dQuote{Source} citation above:
80
81  Sookne, David J. (1973).
82  Bessel Functions of Real Argument and Integer Order.
83  \emph{Journal of Research of the National Bureau of Standards},
84  \bold{77B}, 125--132.
85  \doi{10.6028/jres.077B.012}.
86
87  Cody, William J. (1983).
88  Algorithm 597: Sequence of modified Bessel functions of the first kind.
89  \emph{ACM Transactions on Mathematical Software}, \bold{9}(2), 242--245.
90  \doi{10.1145/357456.357462}.
91
92  Campbell, J.B. (1980).
93  On Temme's algorithm for the modified Bessel function of the third kind.
94  \emph{ACM Transactions on Mathematical Software}, \bold{6}(4), 581--586.
95  \doi{10.1145/355921.355928}.
96
97  Campbell, J.B. (1979).
98  Bessel functions J_nu(x) and Y_nu(x) of float order and float argument.
99  \emph{Computer Physics Communications}, \bold{18}, 133--142.
100  \doi{10.1016/0010-4655(79)90030-4}.
101
102  Temme, Nico M. (1976).
103  On the numerical evaluation of the ordinary Bessel function of the
104  second kind.
105  \emph{Journal of Computational Physics}, \bold{21}, 343--350.
106  \doi{10.1016/0021-9991(76)90032-2}.
107}
108\seealso{
109  Other special mathematical functions, such as
110  \code{\link{gamma}}, \eqn{\Gamma(x)}, and \code{\link{beta}},
111  \eqn{B(x)}.
112}
113\author{
114  Original Fortran code:
115  W. J. Cody, Argonne National Laboratory \cr
116  Translation to C and adaptation to \R:
117  Martin Maechler \email{maechler@stat.math.ethz.ch}.
118}
119\examples{
120require(graphics)
121
122nus <- c(0:5, 10, 20)
123
124x <- seq(0, 4, length.out = 501)
125plot(x, x, ylim = c(0, 6), ylab = "", type = "n",
126     main = "Bessel Functions  I_nu(x)")
127for(nu in nus) lines(x, besselI(x, nu = nu), col = nu + 2)
128legend(0, 6, legend = paste("nu=", nus), col = nus + 2, lwd = 1)
129
130x <- seq(0, 40, length.out = 801); yl <- c(-.5, 1)
131plot(x, x, ylim = yl, ylab = "", type = "n",
132     main = "Bessel Functions  J_nu(x)")
133abline(h=0, v=0, lty=3)
134for(nu in nus) lines(x, besselJ(x, nu = nu), col = nu + 2)
135legend("topright", legend = paste("nu=", nus), col = nus + 2, lwd = 1, bty="n")
136
137## Negative nu's --------------------------------------------------
138xx <- 2:7
139nu <- seq(-10, 9, length.out = 2001)
140## --- I() --- --- --- ---
141matplot(nu, t(outer(xx, nu, besselI)), type = "l", ylim = c(-50, 200),
142        main = expression(paste("Bessel ", I[nu](x), " for fixed ", x,
143                                ",  as ", f(nu))),
144        xlab = expression(nu))
145abline(v = 0, col = "light gray", lty = 3)
146legend(5, 200, legend = paste("x=", xx), col=seq(xx), lty=1:5)
147
148## --- J() --- --- --- ---
149bJ <- t(outer(xx, nu, besselJ))
150matplot(nu, bJ, type = "l", ylim = c(-500, 200),
151        xlab = quote(nu), ylab = quote(J[nu](x)),
152        main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
153abline(v = 0, col = "light gray", lty = 3)
154legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)
155
156## ZOOM into right part:
157matplot(nu[nu > -2], bJ[nu > -2,], type = "l",
158        xlab = quote(nu), ylab = quote(J[nu](x)),
159        main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
160abline(h=0, v = 0, col = "gray60", lty = 3)
161legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)
162
163
164##---------------  x --> 0  -----------------------------
165x0 <- 2^seq(-16, 5, length.out=256)
166plot(range(x0), c(1e-40, 1), log = "xy", xlab = "x", ylab = "", type = "n",
167     main = "Bessel Functions  J_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
168for(nu in sort(c(nus, nus+0.5)))
169    lines(x0, besselJ(x0, nu = nu), col = nu + 2, lty= 1+ (nu\%\%1 > 0))
170legend("right", legend = paste("nu=", paste(nus, nus+0.5, sep=", ")),
171       col = nus + 2, lwd = 1, bty="n")
172
173x0 <- 2^seq(-10, 8, length.out=256)
174plot(range(x0), 10^c(-100, 80), log = "xy", xlab = "x", ylab = "", type = "n",
175     main = "Bessel Functions  K_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
176for(nu in sort(c(nus, nus+0.5)))
177    lines(x0, besselK(x0, nu = nu), col = nu + 2, lty= 1+ (nu\%\%1 > 0))
178legend("topright", legend = paste("nu=", paste(nus, nus + 0.5, sep = ", ")),
179       col = nus + 2, lwd = 1, bty="n")
180
181x <- x[x > 0]
182plot(x, x, ylim = c(1e-18, 1e11), log = "y", ylab = "", type = "n",
183     main = "Bessel Functions  K_nu(x)"); axis(2, at=1)
184for(nu in nus) lines(x, besselK(x, nu = nu), col = nu + 2)
185legend(0, 1e-5, legend=paste("nu=", nus), col = nus + 2, lwd = 1)
186
187yl <- c(-1.6, .6)
188plot(x, x, ylim = yl, ylab = "", type = "n",
189     main = "Bessel Functions  Y_nu(x)")
190for(nu in nus){
191    xx <- x[x > .6*nu]
192    lines(xx, besselY(xx, nu=nu), col = nu+2)
193}
194legend(25, -.5, legend = paste("nu=", nus), col = nus+2, lwd = 1)
195
196## negative nu in bessel_Y -- was bogus for a long time
197curve(besselY(x, -0.1), 0, 10, ylim = c(-3,1), ylab = "")
198for(nu in c(seq(-0.2, -2, by = -0.1)))
199  curve(besselY(x, nu), add = TRUE)
200title(expression(besselY(x, nu) * "   " *
201                 {nu == list(-0.1, -0.2, ..., -2)}))
202}
203\keyword{math}
204
205