R/I-fn.R

Defines functions besselK.nuAsym besselI.nuAsym besselI.ftrms besselIasym bI besselIs

Documented in besselIasym besselI.ftrms besselI.nuAsym besselIs besselK.nuAsym bI

###--- This used to be part /u/maechler/R/MM/NUMERICS/bessel-fn.R

#### "Modified" Bessel Function I_nu(x)
####  ----------------------------------------------
## besselI() - Definition as infinite sum  -- working for "mpfr" numbers, too!
## ---------
besselIs <-
    function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
	     Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
    ## Purpose: besselI() primitively
    ## ----------------------------------------------------------------------
    ## Arguments: (x,nu) as besselI;  nterm: number of terms for "infinite" sum
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 21 Oct 2002, 21:59
    if(length(nu) > 1)
        stop(" 'nu' must be scalar (length 1)!")
    j <- (nterm-1):0 # sum smallest first!
    n <- length(x)
    if(n == 0) return(x)
    if(is(nu, "mpfr")) x <- mpfr(x, precBits = max(64, .getPrec(nu)))
    l.s.j <- outer(j, (x/2), function(X,Y) X*2*log(Y))##-> {nterm x n} matrix
    ##
    isNum <- is.numeric(x) || is.complex(x)
    ## improve accuracy for lgamma(j+1)  for "mpfr" numbers
    ## -- this is very important [evidence: e.g. bI(10000, 1)]
    if(is(l.s.j, "mpfr"))
	j <- mpfr(j, precBits = max(sapply(l.s.j@.Data, slot, "prec")))
    else if(!isNum) j <- as(j, class(x))

    ## underflow (64-bit AMD) for x > 745.1332
    ## for large x, this already overflows to Inf :
    ## s.j <- outer(j, (x^2/4), function(X,Y) Y^X) ##-> {nterm x n} matrix
    ## s.j <- s.j / (gamma(j+1) * gamma(nu+1 + j))  but without overflow :
    log.s.j <-
	if(expon.scaled)
	    l.s.j - rep(x, each=nterm) - lgamma(j+1) - lgamma(nu+1 + j)
	else
	    l.s.j		       - lgamma(j+1) - lgamma(nu+1 + j)
    s.j <-
        if(log) # NB: lsum() works on whole matrix
            lsum(log.s.j) # == log(sum_{j} exp(log.s.j) )
        else ## log I_nu(x) -- trying to avoid overflow/underflow for large x OR large nu
            ## log(s.j) ; e..x <- exp(-x) # subnormal for x > 1024*log(2) ~ 710;
            exp(log.s.j)
    if(log) {
	if(any(i0 <- x == 0)) s.j[i0] <- 0
	if(any(lrgS <- log.s.j[1,] > log(Ceps) + s.j))
	    lapply(x[lrgS], function(x)
		warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
			call.=FALSE))
	nu*log(x/2) + s.j
    } else { ## !log
	s <- colSums(s.j)
	if(any(i0 <- x == 0)) s[i0] <- 0
	if(!all(iFin <- is.finite(s)))
	    stop(sprintf("infinite s for x=%g", x[!iFin][1]))
	if(any(lrgS <- s.j[1,] > Ceps * s))
	    lapply(x[lrgS], function(x)
		warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
			call.=FALSE))
	(x/2)^nu * s
    }
}

## old name [back compatibility]:
bI <- function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
               Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
    .Deprecated("besselIs")
    isNum <- is.numeric(x) || is.complex(x)
    besselIs(x, nu, nterm=nterm, expon.scaled=expon.scaled, log=log, Ceps=Ceps)
}

###--------------- besselI() for large x or also large nu ---
###
###  see --> ~/R/MM/NUMERICS/bessel-large-x.R  for code usage
###          ================================

besselIasym <- function(x, nu, k.max = 10, expon.scaled=FALSE, log=FALSE)
{
    ## Purpose: Asymptotic expansion of Bessel I_nu(x) function   x -> oo
    ##        by Abramowitz & Stegun (9.7.1), p.377 :
    ##
    ## I_a(z) = exp(z) / sqrt(2*pi*z) * f(z,..)  where
    ##   f(z,..) = 1 - (mu-1)/ (8*z) + (mu-1)(mu-9)/(2! (8z)^2) - ...
    ##           = 1- (mu-1)/(8z)*(1- (mu-9)/(2(8z))*(1-(mu-25)/(3(8z))*..))
    ## where  mu = 4*a^2  *and*  |arg(z)| < pi/2
    ## ----------------------------------------------------------------------
    ## Arguments: x, nu, expon.scaled:  as besselI()
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008

    ## Note that for each x the series eventually *DIVERGES*
    ## it should be "stopped" as soon as it *converged* (:-)

    stopifnot(k.max == round(k.max))

    ## I_\nu(z) = exp(z) / sqrt(2*pi*z) * f(z, nu)

    ## First compute  f(x, nu) to order k.max
    ## f <- 1
    d <- 0 ## d = 1 - f  <==>  f = 1 - d
    if(k.max >= 1) {
	## m. <- 4*nu^2
	x8 <- 8*x
        for(k in k.max:1) {
            ## mu <- 4*nu^2; d <- (1 - d)*(mu - (2*k-1)^2)/(k*x8)
            ## mu - (2k-1)^2 = (2nu - (2k-1)) (2nu + (2k-1)) =
            ##               = (2(nu-k)+1)(2(nu+k)-1)
            d <- (1 - d)*((2*(nu-k)+1)*(2*(nu+k)-1))/(k*x8)
        }
    }
    if(log) {
        ## f = 1 - d  ==> log(f) = log1p(-d) :
        (if(expon.scaled) log1p(-d) else x + log1p(-d)) - log(2*pi*x) / 2
    } else {
        (if(expon.scaled) (1-d) else exp(x)*(1-d)) / sqrt(2*pi*x)
    }
}

besselI.ftrms <- function(x, nu, K = 20)
{
    ## Purpose: all *Terms* in besselIasym()
    ##        by Abramowitz & Stegun (9.7.1), p.377 :
    ##   f(x,..) = 1 - (mu-1)/ (8*x) + (mu-1)(mu-9)/(2! (8x)^2) - ...
    ##           = 1- (mu-1)/(8x)*(1- (mu-9)/(2(8x))*(1-(mu-25)/(3(8x))*..))
    ## where  mu = 4*nu^2
    ## ----------------------------------------------------------------------
    ## Arguments: x, nu, expon.scaled:  as besselI()
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008

    stopifnot(length(x) == 1, length(nu) == 1, length(K) == 1,
              x > 0, K == round(K), K >= 0)

    ## I_\nu(x) = exp(x) / sqrt(2*pi*x) * f(x, nu)

    kk <- seq_len(K)
    mu <- 4*nu^2
    x8 <- 8*x
    multf <- - (mu - (2*kk-1)^2)/(kk*x8)
    ser <- cumprod(multf)

    ## now, for k-term approx. f_k  of  f(x,nu)   we have
    ##  f_k = 1 + sum(ser[1:k]), i.e.    ##  (f_k)_k = 1 + cumsum(c(0, ser))[k+1]  for k= 0,1,...,K
    ser
}


###----------------------------------------------------------------------------

### When  BOTH  nu and x  are large :

besselI.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
    ## Purpose: Asymptotic expansion of Bessel I_nu(x) function
    ##        when BOTH  nu and x  are large

    ## Abramowitz & Stegun , p.378, __ 9.7.7. __
    ##
    ## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
    ##                * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }

    ## where   t = 1 / sqrt(1 + z^2),
    ##       eta = sqrt(1 + z^2) + log(z / (1 + sqrt(1+z^2)))
    ##
    ##
    ## and u_k(t)  from  p.366  __ 9.3.9 __

    ## u0(t) = 1
    ## u1(t) = (3*t - 5*t^3)/24
    ## u2(t) = (81*t^2 - 462*t^4 + 385*t^6)/1152
    ## ...

    ## with recursion  9.3.10    for  k = 0, 1, .... :
    ##
    ## u_{k+1}(t) = t^2/2 * (1 - t^2) * u'_k(t) +
    ##            1/8  \int_0^t (1 - 5*s^2)* u_k(s) ds

    ## ----------------------------------------------------------------------
    ## Arguments: x, nu, expon.scaled:  as besselI()
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 22 Nov 2008, 15:22

    stopifnot(k.max == round(k.max),
              0 <= k.max, k.max <= 4)

    z <- x/nu # -> computing I_nu(z * nu)
    sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
    t <- 1/sz
    eta <- (if(expon.scaled) ## <==> * exp(-x)  ==  exp(- z * nu)  <==> "eta - z"
            ## sz - z = sqrt(1 + z^2) - z  =  1/(sqrt(..) + z)
            1/(z + sz) else sz) +
                log(z / (1 + sz))

    ## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
    ##                * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }
    if(k.max == 0)
        d <- 0
    else { ## k.max >= 1 --- Find the  Debye  polynomials u_j(t) in
           ##                /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
        t2 <- t^2
        u1.t <- (t*(3 - 5*t2))/24
        d <-
            if(k.max == 1) {
                u1.t/nu
            } else { ## k.max >= 2
                u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
                    t2*(81 + t2*(-462 + t2*385)) / 1152
                if(k.max == 2)
                    (u1.t + u2.t/nu)/nu
                else { ## k.max >= 3
                    u3.t <- t*t2*(30375 +
                                  t2*(-369603 +
                                      t2*(765765 - t2*425425)))/ 414720
                    if(k.max == 3)
                        (u1.t + (u2.t + u3.t/nu)/nu)/nu
                    else { ## k.max >= 4
			u4.t <- t2*t2*(4465125 +
				       t2*(-94121676 +
					   t2*(349922430 +
					       t2*(-446185740 + t2*185910725))))/39813120
			if(k.max == 4)
			    (u1.t + (u2.t + (u3.t + u4.t/nu)/nu)/nu)/nu
                    }
                }
            }
    }

    if(log) {
        log1p(d) + nu*eta - (log(sz) + log(2*pi*nu))/2
    } else {
        (1+d) * exp(nu*eta)/sqrt(2*pi*nu*sz)
    }
} ## {besselI.nuAsym}

besselK.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
    ## Purpose: Asymptotic expansion of Bessel K_nu(x) function
    ##        when BOTH  nu and x  are large

    ## Abramowitz & Stegun , p.378, __ 9.7.8. __
    ##
    ## K_nu(nu * z) ~ sqrt(pi/(2*nu)) * exp(-nu*eta)/(1+z^2)^(1/4) *
    ##                * {1 - u_1(t)/nu + u_2(t)/nu^2 - ... }

    ## { see besselI.nuAsym() above, for t, eta, u_k ...}
    ## ----------------------------------------------------------------------
    ## Arguments: x, nu, expon.scaled:  as besselK()
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 23 Nov 2009, 13:39

    stopifnot(k.max == round(k.max),
              0 <= k.max, k.max <= 4)

    z <- x/nu # -> computing K_nu(z * nu)
    sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
    t <- 1/sz
    eta <- (if(expon.scaled) ## <==> * exp(-x)  ==  exp(- z * nu)  <==> "eta - z"
            ## sz - z = sqrt(1 + z^2) - z  =  1/(sqrt(..) + z)
            1/(z + sz) else sz) +
                log(z / (1 + sz))

    ## K_nu(nu * z) ~  [see above]
    if(k.max == 0)
        d <- 0
    else { ## k.max >= 1 --- Find the  Debye  polynomials u_j(t) in
           ##                /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
        t2 <- t^2
        u1.t <- (t*(3 - 5*t2))/24
        d <-
            if(k.max == 1) {
                - u1.t/nu
            } else { ## k.max >= 2
                u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
                    t2*(81 + t2*(-462 + t2*385)) / 1152
                if(k.max == 2)
                    (- u1.t + u2.t/nu)/nu
                else { ## k.max >= 3
                    u3.t <- t*t2*(30375 +
                                  t2*(-369603 +
                                      t2*(765765 - t2*425425)))/ 414720
                    if(k.max == 3)
                        (- u1.t + (u2.t - u3.t/nu)/nu)/nu
                    else { ## k.max >= 4
                        u4.t <- t2*t2*(4465125 +
                                           t2*(-94121676 +
                                                   t2*(349922430 +
                                                           t2*(-446185740 + t2*185910725))))/39813120
                        if(k.max == 4)
                            (- u1.t + (u2.t + (-u3.t + u4.t/nu)/nu)/nu)/nu
                        else stop("k.max >= 5 not yet implemented")
                    }
                }
            }
    }

    if(log) {
        log1p(d) - nu*eta - (log(sz) - log(pi/(2*nu)))/2
    } else {
        (1+d) * exp(-nu*eta)*sqrt(pi/(2*nu * sz))
    }
} ## {besselK.nuAsym}

Try the Bessel package in your browser

Any scripts or data that you put into this service are public.

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