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