sdetorus (version 0.1.10)

logBesselI0Scaled: Efficient computation of Bessel related functions

Description

Computation of \(\log(I_0(x))-x\) and the inverse of \(A_1(k)=\frac{I_0(k)}{I_1(k)}\).

Usage

logBesselI0Scaled(x, splineApprox = TRUE)
a1Inv(x, splineApprox = TRUE)

Value

A vector of the same length as x.

Arguments

x

evaluation vector. For logBesselI0Scaled, x must contain non-negative values. For a1Inv, x must be in \([0, 1]\).

splineApprox

whether to use a pre-computed spline approximation (faster) or not.

Details

Both functions may rely on pre-computed spline interpolations (logBesselI0ScaledSpline and a1InvSpline). Otherwise, a call to besselI is done for \(\log(I_0(x))-x\) and \(A_1(k)=x\) is solved numerically. The data in which the interpolation is based is given in the examples.

For x larger than 5e4, the asymptotic expansion of besselIasym is employed.

Examples

# \donttest{
# Data employed for log besselI0 scaled
x1 <- c(seq(0, 1, by = 1e-4), seq(1 + 1e-2, 10, by = 1e-3),
        seq(10 + 1e-1, 100, by = 1e-2), seq(100 + 1e0, 1e3, by = 1e0),
        seq(1000 + 1e1, 5e4, by = 2e1))
logBesselI0ScaledEvalGrid <- log(besselI(x = x1, nu = 0,
                                         expon.scaled = TRUE))
# save(list = "logBesselI0ScaledEvalGrid",
#      file = "logBesselI0ScaledEvalGrid.rda", compress = TRUE)

# Data employed for A1 inverse
x2 <- rev(c(seq(1e-04, 0.9 - 1e-4, by = 1e-4),
            seq(0.9, 1 - 1e-05, by = 1e-5)))
a1InvEvalGrid <- sapply(x2, function(k) {
  uniroot(f = function(x) k - besselI(x, nu = 1, expon.scaled = TRUE) /
          besselI(x, nu = 0, expon.scaled = TRUE),
          lower = 1e-06, upper = 1e+05, tol = 1e-15)$root
})
# save(list = "a1InvEvalGrid", file = "a1InvEvalGrid.rda", compress = TRUE)
# }
# Accuracy logBesselI0Scaled
x <- seq(0, 1e3, l = 1e3)
summary(logBesselI0Scaled(x = x, splineApprox = TRUE) -
        logBesselI0Scaled(x = x, splineApprox = FALSE))

# Accuracy a1Inv
y <- seq(0, 1 - 1e-4, l = 1e3)
summary(a1Inv(x = y, splineApprox = TRUE) -
        a1Inv(x = y, splineApprox = FALSE))