besselJs: Bessel J() function Simple Series Representation

In Bessel: Computations and Approximations for Bessel Functions

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 of 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

See Also

This package BesselJ, base besselJ, etc

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

Example output

Warning messages:
1: In FUN(X[[i]], ...) :  'nterm=800' may be too small for x=4
2: In FUN(X[[i]], ...) :  'nterm=800' may be too small for x=5
3: In FUN(X[[i]], ...) :  'nterm=800' may be too small for x=6
4: In FUN(X[[i]], ...) :  'nterm=800' may be too small for x=7
 [1]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
 [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[11]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[16]  0.000000e+00  0.000000e+00  0.000000e+00 3.378127e-305 1.673577e-300
[21] 4.760010e-296
There were 19 warnings (use warnings() to see them)
Loading required package: Rmpfr
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following object is masked from 'package:gmp':

    outer

The following objects are masked from 'package:stats':

    dbinom, dnorm, dpois, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

Bessel documentation built on May 2, 2019, 4:38 p.m.