Bessel (version 0.7-0)

besselJs: 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,
         use.log = log || any(nu * log(x/2) > .Machine$double.max.exp * log(2)),
         Ceps = if(isNum) 8e-16 else 2^(- x.@.Data[[1]]@prec))

Value

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

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.

use.log

logical indicating if the logarithm should work in logarithmic scale notably even when log = FALSE. Its smart default, makes besselJs(4, 1e10) to work, correctly giving 0, where it gave NaN previously.

Ceps

a relative error tolerance for checking if nterm has been sufficient. The default is “correct” for double precision and also for multiprecision objects.

Author

Martin Maechler

References

Abramowitz, M., and Stegun, I. A. (1964--1972). Handbook of mathematical functions (NBS AMS series 55, U.S. Dept. of Commerce). https://personal.math.ubc.ca/~cbm/aands/page_360.htm

See Also

This package BesselJ(), base besselJ(), etc

Examples

stopifnot(all.equal(besselJs(1:10, 1), # our R code --> 4 warnings, for x = 4:7
                    besselJ (1:10, 1)))# internal C code w/ different algorithm

## Large 'nu' ...
x <- (0:20)/4
if(interactive()) op <- options(nwarnings = 999)
(bx <- besselJ(x, nu=200))# base R's -- gives 19 (mostly wrong) warnings about precision lost
## Visualize:
bj <- curve(besselJ(1, x), .001, 2^10, log="xy", n=1001,
            main=quote(J[nu](1)), xlab = quote(nu), xaxt="n", yaxt="n") # 50+ warnings
eaxis <- if(!requireNamespace("sfsmisc")) axis else sfsmisc::eaxis
eaxis(1, sub10 = c(-2,3)); eaxis(2)
abline(h = (J01 <- besselJ(1, 0)), col=adjustcolor(2, 1/2))
axis(4, at=J01, quote(J[0](1)), las=2, col=2, col.axis=2, line = -1)
bj6 <- curve(besselJ(6, x), add=TRUE, n=1001, col=adjustcolor(2, 1/2), lwd=2)
plot(y~x, as.data.frame(bj6), log="x", type="l", col=2, lwd=2,
     main = quote(J[nu](6)), xlab = quote(nu), xaxt="n")
eaxis(1, sub10=3); abline(h=0, lty=3)

## large nu (non-log case) -- gave NaN wrongly (typically when it should have given 0):
stopifnot(exprs = {
    besselJs(c(.1, 1:10),  999) == 0
    besselJs(c(.1, 1:10), 9999) == 0
    besselJs(c(.1, 1:1000), 1e10) == 0
 ## besselJs(c(.1, 1:10), 1e20) == 0 -- another *BUG*: Error ... lssum found non-positive sums
})


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
if(interactive()) options(op) # reset 'nwarnings'