bessJsmlNu: Bessel J_nu(x) for Very Small Order nu

In Bessel: Computations and Approximations for Bessel Functions

Bessel J_\nu(x) for Very Small Order \nu

Description

Use a few terms Taylor approximation for Bessel's J_\nu(x) for fixed x and \nu \to 0.

Usage

bessJsmlNu(x, nu, nTrms, m.max)

Arguments

x, nu

numeric Bessel J() arguments, see e.g., besselJ

.

nTrms	integer, one of 0, 1, or 2, denoting the number of Taylor series terms to be used, nTrms = 0 simply returning BesselJ(x, 0).
m.max	positive integer specifying the Taylor series order used.

Details

A. & S. (9.1.10) p. 360 (Ascending Series) (9.1.64) p. 362 (derivative wrt order)

............ show derivation? ...........

Value

The function is explicitly using Vectorize(<fn>, c("x", "nu")), hence a numeric vector conforming to (e.g., ) x + nu.

Author(s)

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's BesselJ, base R's besselJ

Examples

(nuSml <- 2^-c(seq(30, 53, by=1/2), 75, 100, 300, 800, 1022, 1074, Inf))
## 9.3e-10 6.6e-10 ... 1.1e-16 .. 2.22e-308 4.94e-324 0.0
## bug in R (at least up to Jan. 1 2026) :
options(digits = 14) # show more digits
Jsml <- sapply(nuSml, \(nu) besselJ(pi/2, nu))
tail(cbind(nuSml, Jsml), 20) # complete "divergence" for  nu <~ 10^-{16}
Jsm.nT.1 <- bessJsmlNu(pi/2, nuSml, nTrms = 1, m.max = 30)
Jsm.nT.2 <- bessJsmlNu(pi/2, nuSml, nTrms = 2, m.max = 30)
cbind(nuSml, Jsml, Jsm.nT.1, Jsm.nT.2)
table(Jsm.nT.2 == Jsm.nT.1) # all TRUE  (really ???)
          all.equal(Jsm.nT.2, Jsm.nT.1, tolerance = 0) # show
stopifnot(all.equal(Jsm.nT.2, Jsm.nT.1, tolerance = 4e-16))

## second example;  nu  *not* very small
x. <- 2
nu2 <- 2^-seq(1, 30, by = 1/8)
Jsm2.1 <- bessJsmlNu(x., nu2, nTrms = 1, m.max = 40)
Jsm2.2 <- bessJsmlNu(x., nu2, nTrms = 2, m.max = 40)
table(Jsm2.2 == Jsm2.1) # now they *do* differ
bessJ2 <- besselJ(x., nu2)
BessJ2 <- sapply(nu2, function(nu) BesselJ(x., nu))
## J() function values: ---------------------
matplot(nu2, cbind(bessJ2, BessJ2, Jsm2.1, Jsm2.2), log = "x", type = "l")
sum(i <- nu2 < 1e-4)
stopifnot(all.equal(BessJ2[i], Jsm2.2[i], tolerance = 4e-14),
          all.equal(BessJ2[i], Jsm2.1[i], tolerance = 4e-9))

##  |error|  wrt BesselJ() ---------------
matplot(nu2, abs(cbind(bessJ2, Jsm2.1, Jsm2.2) - BessJ2), log = "xy",
        type = "l", ylim = c(5e-17, 1e-10),
        xlab = quote("nu" == nu),
	main = substitute(list(abs("<bess J>"(x, nu) - BesselJ(x, nu)), x == XX), list(XX = formatC(x.))))
legend("topleft", lty = 1:5, lwd=2, col=1:6, bty = "n",
       legend = paste0(c("besselJ(", "bessJsmlNu(*, nTrms = 1", "bessJsmlNu(*, nTrms = 2"),")"))

Bessel documentation built on Jan. 2, 2026, 3 a.m.