Bessel (version 0.7-1)

besselI.nuAsym: Asymptotic Expansion of Bessel I(x,nu) and K(x,nu) for Large nu (and x)

Description

Compute Bessel functions \(I_{\nu}(x)\) and \(K_{\nu}(x)\) for large \(\nu\) and possibly large \(x\), using asymptotic expansions in Debye polynomials.

Usage

besselI.nuAsym(x, nu, k.max, expon.scaled = FALSE, log = FALSE)
besselK.nuAsym(x, nu, k.max, expon.scaled = FALSE, log = FALSE)

Value

a numeric vector of the same length as the long of x and

nu. (usual argument recycling is applied implicitly.)

Arguments

x

numeric or complex, with real part \(\ge 0\).

nu

numeric; The order (maybe fractional!) of the corresponding Bessel function.

k.max

integer number of terms in the expansion. Must be in 0:5, currently.

expon.scaled

logical; if TRUE, the results are exponentially scaled, the same as in the corresponding BesselI() and BesselK() functions in order to avoid overflow (\(I_{\nu}\)) or underflow (\(K_{\nu}\)), respectively.

log

logical; if TRUE, \(\log(f(.))\) is returned instead of \(f\).

Author

Martin Maechler

Details

Abramowitz & Stegun , page 378, has formula 9.7.7 and 9.7.8 for the asymptotic expansions of \(I_{\nu}(x)\) and \(K_{\nu}(x)\), respectively, also saying When \(\nu \to +\infty\), these expansions (of \(I_{\nu}(\nu z)\) and \(K_{\nu}(\nu z)\)) hold uniformly with respect to \(z\) in the sector \(|arg z| \le \frac{1}{2} \pi - \epsilon\), where \(\epsilon\) is an arbitrary positive number. and for this reason, we require \(\Re(x) \ge 0\).

The Debye polynomials \(u_k(x)\) are defined in 9.3.9 and 9.3.10 (page 366).

References

Abramowitz, M., and Stegun, I. A. (1964, etc). Handbook of mathematical functions, pp. 366, 378, e.g.,

See Also

From this package Bessel: BesselI(); further, besselIasym() for the case when \(x\) is large and \(\nu\) is small or moderate (using A.&S. (9.7.1) and (9.7.2)).

Further, from base: besselI, etc.

Examples

x <- c(1,2, 10, 20, 50, 100, 1000, 1e4)
nu <- c(outer(c(1,2,5), 10^(0:9)))
ks <- 0:4; ks <- setNames(ks, paste0("k=",ks))

if(requireNamespace("sfsmisc")) { # Mostly for plotting
    formatN <- if(packageVersion("sfsmisc") >= "1.1.26")
                   sfsmisc::frmtNum else sfsmisc::formatN
    eaxis    <- sfsmisc::eaxis
    mult.fig <- sfsmisc::mult.fig
} else { ## no {sfsmisc}
    formatN <- function(x) formatC(x, width=1)
    eaxis <- axis
    mult.fig <- function(nr.plots, main=NULL, ...) { # cheap but working substitute
        mfrow <- n2mfrow(nr.plots)
        tit.wid <- if(is.null(main)) 0 else 1 + 1.5 * par("cex.main")
        old.par <- par(mfrow = mfrow, oma = c(0, 0, tit.wid, 0),
                       mar = 0.1 + c(4, 4, 2, 1), mgp = c(1.5, 0.6, 0))
        if (!is.null(main)) {
            plot.new(); mtext(main, side = 3, outer = TRUE, cex=par("cex.main")); par(new = TRUE)
        }
        invisible(list(old.par = old.par))
      }
}

nus <- setNames(nu, paste0("nu=",formatN(nu)))

mI <- sapply(simplify = "array", ks, function(k.)
            sapply(nus, function(n.)
                   besselI.nuAsym(x, nu=n., k.max = k., log = TRUE)))

mK <- sapply(simplify = "array", ks, function(k.)
            sapply(nus, function(n.)
                   besselK.nuAsym(x, nu=n., k.max = k., log = TRUE)))
str(mK) # 7 x 30 x 5 array

##------------ Plotting  besselI :  I_[nu](x = x0) vs nu -----------------------------------

str(op <- mult.fig(length(x))$old.par); par("mfrow")
for(ix in seq_along(x)) {
    matplot(nu, mI[ix,,], type = "l", log = "x",
            xaxt = "n", xlab = quote(nu), ylab = quote(I[nu](x)),
            main = substitute(I[nu](x==X, nu, k.max == 0:4), list(X = (xi <- x[ix]))))
    eaxis(1)
}
par(op)

##------------ Plotting  besselK :  K_[nu](x = x0) vs nu -----------------------------------

str(mult.fig(length(x), main = expression(log(K[nu](x)) ~ vs ~ nu ~~ "x in log-scale"))$old.par)
par("mfrow")
for(ix in seq_along(x)) {
    y <- mK[ix,,]
    matplot(nu, y, type = "l", log = "x", ylim = c(min(y), min(max(y), 3*abs(min(y)))),
            xaxt = "n", xlab = quote(nu), ylab = quote(K[nu](x)),
            main = substitute(K[nu](x==X, nu, k.max == 0:4), list(X = (xi <- x[ix]))))
    abline(h = 0, col = "gray20", lty = 3)
    eaxis(1)
}
par(op)

str(mult.fig(length(x), main = expression(log(K[nu](x)) ~ vs ~ nu ~~ "in log-log"))$old.par)
for(ix in seq_along(x)) {
    matplot(nu, mK[ix,,], type = "l", log = "xy", # --> 14 warnings: each plot h
            xaxt = "n", yaxt = "n", xlab = quote(nu), ylab = quote(K[nu](x)),
            main = substitute(K[nu](x==X, nu, k.max == 0:4), list(X = (xi <- x[ix]))))
    eaxis(1); eaxis(2)
}
par(op)

## Compare with  besselK() -- which I can use expon.scaled !
dimnames(mI)[[1]] <- dimnames(mK)[[1]] <- paste0("x=", formatN(x))
str(mK)

rbind(t(mK[,"nu=200",]), log.besselK = log(besselK(x, nu=200, expon.scaled=TRUE)) - x)
##--> apart from k=0  the asymptotic approximation seems good already for k=1 here

rbind(t(mK[,"nu=500",]), log.besselK = log(besselK(x, nu=500, expon.scaled=TRUE)) - x)
##--> apart from k=0  the asymptotic approximation seems good already for k=1 here
##---> need other, smaller x[] range -  to get non-Inf besselK()