EVOLUTION-MANAGER
Edit File: Bessel.html
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml"><head><title>R: Bessel Functions</title> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> <link rel="stylesheet" type="text/css" href="R.css" /> </head><body> <table width="100%" summary="page for Bessel {base}"><tr><td>Bessel {base}</td><td style="text-align: right;">R Documentation</td></tr></table> <h2>Bessel Functions</h2> <h3>Description</h3> <p>Bessel Functions of integer and fractional order, of first and second kind, <i>J(nu)</i> and <i>Y(nu)</i>, and Modified Bessel functions (of first and third kind), <i>I(nu)</i> and <i>K(nu)</i>. </p> <h3>Usage</h3> <pre> besselI(x, nu, expon.scaled = FALSE) besselK(x, nu, expon.scaled = FALSE) besselJ(x, nu) besselY(x, nu) </pre> <h3>Arguments</h3> <table summary="R argblock"> <tr valign="top"><td><code>x</code></td> <td> <p>numeric, <i>≥ 0</i>.</p> </td></tr> <tr valign="top"><td><code>nu</code></td> <td> <p>numeric; The <em>order</em> (maybe fractional and negative) of the corresponding Bessel function.</p> </td></tr> <tr valign="top"><td><code>expon.scaled</code></td> <td> <p>logical; if <code>TRUE</code>, the results are exponentially scaled in order to avoid overflow (<i>I(nu)</i>) or underflow (<i>K(nu)</i>), respectively.</p> </td></tr> </table> <h3>Details</h3> <p>If <code>expon.scaled = TRUE</code>, <i>exp(-x) I(x;nu)</i>, or <i>exp(x) K(x;nu)</i> are returned. </p> <p>For <i>nu < 0</i>, formulae 9.1.2 and 9.6.2 from Abramowitz & Stegun are applied (which is probably suboptimal), except for <code>besselK</code> which is symmetric in <code>nu</code>. </p> <p>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 <code>nu</code>, say in the order of millions, the current algorithms are rarely useful. </p> <h3>Value</h3> <p>Numeric vector with the (scaled, if <code>expon.scaled = TRUE</code>) values of the corresponding Bessel function. </p> <p>The length of the result is the maximum of the lengths of the parameters. All parameters are recycled to that length. </p> <h3>Author(s)</h3> <p>Original Fortran code: W. J. Cody, Argonne National Laboratory <br /> Translation to C and adaptation to <span style="font-family: Courier New, Courier; color: #666666;"><b>R</b></span>: Martin Maechler <a href="mailto:maechler@stat.math.ethz.ch">maechler@stat.math.ethz.ch</a>. </p> <h3>Source</h3> <p>The C code is a translation of Fortran routines from <a href="http://www.netlib.org/specfun/ribesl">http://www.netlib.org/specfun/ribesl</a>, <span class="samp">../rjbesl</span>, etc. The four source code files for bessel[IJKY] each contain a paragraph “Acknowledgement” and “References”, a short summary of which is </p> <dl> <dt>besselI</dt><dd><p>based on (code) by David J. Sookne, see Sookne (1973)... Modifications... An earlier version was published in Cody (1983).</p> </dd> <dt>besselJ</dt><dd><p>as <code>besselI</code></p> </dd> <dt>besselK</dt><dd><p>based on (code) by J. B. Campbell (1980)... Modifications...</p> </dd> <dt>besselY</dt><dd><p>draws heavily on Temme's Algol program for <i>Y</i>... and on Campbell's programs for <i>Y_ν(x)</i> .... ... heavily modified.</p> </dd> </dl> <h3>References</h3> <p>Abramowitz, M. and Stegun, I. A. (1972). <em>Handbook of Mathematical Functions</em>. Dover, New York; Chapter 9: Bessel Functions of Integer Order. </p> <p>In order of “Source” citation above: </p> <p>Sockne, David J. (1973). Bessel Functions of Real Argument and Integer Order. <em>Journal of Research of the National Bureau of Standards</em>, <b>77B</b>, 125–132. </p> <p>Cody, William J. (1983). Algorithm 597: Sequence of modified Bessel functions of the first kind. <em>ACM Transactions on Mathematical Software</em>, <b>9</b>(2), 242–245. doi: <a href="https://doi.org/10.1145/357456.357462">10.1145/357456.357462</a>. </p> <p>Campbell, J.B. (1980). On Temme's algorithm for the modified Bessel function of the third kind. <em>ACM Transactions on Mathematical Software</em>, <b>6</b>(4), 581–586. doi: <a href="https://doi.org/10.1145/355921.355928">10.1145/355921.355928</a>. </p> <p>Campbell, J.B. (1979). Bessel functions J_nu(x) and Y_nu(x) of float order and float argument. <em>Computer Physics Communications</em>, <b>18</b>, 133–142. doi: <a href="https://doi.org/10.1016/0010-4655(79)90030-4">10.1016/0010-4655(79)90030-4</a>. </p> <p>Temme, Nico M. (1976). On the numerical evaluation of the ordinary Bessel function of the second kind. <em>Journal of Computational Physics</em>, <b>21</b>, 343–350. doi: <a href="https://doi.org/10.1016/0021-9991(76)90032-2">10.1016/0021-9991(76)90032-2</a>. </p> <h3>See Also</h3> <p>Other special mathematical functions, such as <code><a href="Special.html">gamma</a></code>, <i>Γ(x)</i>, and <code><a href="Special.html">beta</a></code>, <i>B(x)</i>. </p> <h3>Examples</h3> <pre> 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)})) </pre> <hr /><div style="text-align: center;">[Package <em>base</em> version 3.6.0 <a href="00Index.html">Index</a>]</div> </body></html>