R: Bessel J() function Simple Series Representation

besselJs {Bessel}

R Documentation

Bessel J() function Simple Series Representation

Description

Computes the modified Bessel J function, using one of its basic definitions as an infinite series, e.g. A. & S., p.360, (9.1.10). The implementation is pure R, working for numeric, complex, but also e.g., for objects of class "mpfr" from package Rmpfr.

Usage

besselJs(x, nu, nterm = 800, log = FALSE,
         Ceps = if (isNum) 8e-16 else 2^(-x@.Data[[1]]@prec))

Arguments

`x`	numeric or complex vector, or of another `class` for which arithmetic methods are defined, notably objects of class `mpfr`.
`nu`	non-negative numeric (scalar).
`nterm`	integer indicating the number of terms to be used. Should be in the order of `abs(x)`, but can be smaller for large x. A warning is given, when `nterm` was possibly too small. (Currently, many of these warnings are wrong, as
`log`	logical indicating if the logarithm `log J.()` is required.
`Ceps`	a relative error tolerance for checking if `nterm` has been sufficient. The default is “correct” for double precision and also for multiprecision objects.

Value

a “numeric” (or complex or "mpfr") vector of the same class and length as x.

Author(s)

Martin Maechler

References

Abramowitz, M., and Stegun, I. A. (1955, etc). Handbook of mathematical functions (NBS AMS series 55, U.S. Dept. of Commerce). http://people.math.sfu.ca/~cbm/aands/page_360.htm

Examples

stopifnot(all.equal(besselJs(1:10, 1), # our R code
                    besselJ (1:10, 1)))# internal C code w/ different algorithm

## Large 'nu' ...
x <- (0:20)/4
(bx <- besselJ(x, nu=200))# base R's -- gives (mostly wrong) warnings
if(require("Rmpfr")) { ## Use high precision, notably large exponent range, numbers:
  Bx <- besselJs(mpfr(x, 64), nu=200)
  all.equal(Bx, bx, tol = 1e-15)# TRUE -- warnings were mostly wrong; specifically:
  cbind(bx, Bx)
  signif(asNumeric(1 - (bx/Bx)[19:21]), 4) # only [19] had lost accuracy

  ## With*out* mpfr numbers -- using log -- is accurate (here)
  lbx <- besselJs(     x,      nu=200, log=TRUE)
  lBx <- besselJs(mpfr(x, 64), nu=200, log=TRUE)
  cbind(x, lbx, lBx)
  stopifnot(all.equal(asNumeric(log(Bx)), lbx, tol=1e-15),
	    all.equal(lBx, lbx, tol=4e-16))
} # Rmpfr

[Package Bessel version 0.6-0 Index]