tests/IJKY.R

In Bessel: Computations and Approximations for Bessel Functions

library(Bessel)

#### Test cases for the Bessel functions I(), J(), K(), Y()
####                    -----------------------------------

if(getRversion() < "2.15.0")
paste0 <- function(...) paste(..., sep="")

### --- For real arguments -- Comparisons with bessel[IJYK]()

x <- c(1e-6, 0.1, 1:10, 20, 100, 200)# larger x : less (relative) accuracy (??)
nus <- c(-6, -5.9, -1.1, .1, 1:4, 10, 20)
## From R 2.10.1 (or 2.10.0 patched, >= Nov.23 2009) on, this works, previously
## 	R"s  besselJ() had an inaccuracy that we have corrected here !
## FIXME ? ---- we currently need to fudge for negative nu
## note that (nu != -6 but nu ~ -6, |x| << |nu|) is still unstable,
## since sin(nu*pi) ~ 0 (but not exactly!) and besselK() is large
nS <- 10
for(F in c("I","J","K","Y")) {
    cat("\n", F," ---  nu : ")
    zF <- get(paste0("Bessel", F))
    FF <- get(paste0("bessel", F))
    stopifnot(is.function(zF), is.function(FF))
    for(nu in nus) {
        x. <- x
        spec <- FALSE ## nu < 0 && F %in% c("I","J")
        cat(format(nu), if(spec)"* " else " ", sep="")
        zr <- zF(x., nu, nSeq = nS)
        rr <- outer(x., nu+ sign(nu)*(seq_len(nS)-1), FF)
        stopifnot(all.equal(zr, rr, tol = 500e-16))
    }; cat("\n")
}

zr[,1]

### *Large* arguments ,
### However, base::bessel*(): only I() and K() have 'expon.scaled'

x <- 1000*(1:20)
str(rI <- BesselI(x, 10, nSeq = 5, expon.scaled=TRUE))
rI
str(rK <- BesselK(x, 10, nSeq = 5, expon.scaled=TRUE))

if(getRversion() >= "2.8.2") { ## besselI(), besselJ() with larger working range
    r2 <- outer(x, 10+seq_len(5)-1,
                besselI, expon.scaled=TRUE)
    stopifnot(all.equal(rI, r2, tol = 8e-16))
} ## R >= 2.8.2

r2 <- outer(x, 10+seq_len(5)-1,
            besselK, expon.scaled=TRUE)
stopifnot(all.equal(rK, r2, tol = 8e-16))

###--------------------- Complex  Arguments ------------------------------

## I would expect that   besselIasym()  work too {should according to Abramowitz&Stegun!}
## but they seem complete "nonsense":

## e.g.
z <- 10000*(1 + 1i)
BesselI       (z, 10, expon.scaled=TRUE) ##  -0.003335734 + 0.000266159i
besselIasym   (z, 10, expon.scaled=TRUE) ## very different
besselI.nuAsym(z, 10, expon.scaled=TRUE, k.max=3) ## (same)

##  J(i * z, nu) =  c(nu) * I(z, nu)

## for *integer* nu, it's simple
N <- 100
set.seed(1)
for(nu in 0:20) {
    cat(nu, "")
    z <- complex(re = rnorm(N),
                 im = rnorm(N))
    r <- BesselJ(z * 1i, nu) / BesselI(z, nu)
    stopifnot(all.equal(r, rep.int(exp(nu/2*pi*1i), N)))
}; cat("\n")

nus <- round(sort(rlnorm(20)), 2)

if(FALSE) { ## Bug ??
## For fractional nu, there's a problem (?) :
for(nu in nus) {
    cat("nu=", formatC(nu, wid=6),":")
    z <- complex(re = rnorm(N),
                 im = rnorm(N))
    r <- BesselJ(z * 1i, nu) / BesselI(z, nu)
    r.Theory <- exp(nu/2*pi*1i)
    cat("correct:"); print(table(Ok <- abs(1 - r /r.Theory) < 1e-7))
    cat("log(<wrong>) / (i*pi/2) :",
        format(unique(log(signif(r[!Ok], 6)))/(pi/2 * 1i)),
        "\n")
}

}# not for now


### Replicate some testing "ideas" from  zqcbi.f (TOMS 644 test program)

##   zqcbi is a quick check routine for the complex i bessel function
##   generated by subroutine zbesi.
##
##   zqcbk generates sequences of i and k bessel functions from
##   zbesi and cbesk and checks the wronskian evaluation
##
##         I(nu,z)*K(nu+1,z) + I(nu+1,z)*K(nu,z) = 1/z
##
##   in the right half plane and a modified form
##
##        I(nu+1,z)*K(nu,zr) - I(nu,z)*K(nu+1,zr) = c/z
##
##   in the left half plane where zr: = -z and c := exp(nu*sgn*pi*i) with
##   sgn=+1 for Im(z) >= 0 and sgn=-1 for Im(z) < 0.          ^^^
##                                                          ( ||| corrected, MM)
N <- 100
nS <- 20
set.seed(257)

nus. <- unique(sort(c(nus,10*nus)))
## For exploration nus. <- (1:80)/4
for(i in seq_along(nus.)) {
    nu <- nus.[i]
    cat(nu, "")
    z <- complex(re = rnorm(N),
                 im = rnorm(N))
    P <- Re(z) >= 0
    rI  <- BesselI( z,     nu, nSeq = 1+nS) # -> for (nu, nu+1, ...,nu+nS)
    rKp <- BesselK( z[ P], nu, nSeq = 1+nS)
    rKm <- BesselK(-z[!P], nu, nSeq = 1+nS)
    ##
    sgn <- ifelse(Im(z) >= 0, +1, -1)
    Izp <- 1 / z [ P]
    for(j in 1:nS) {
        nu.. <- nu + j-1
        allEQ <- function(x,y) all.equal(x,y, tol= max(1,nu..)*nS*2e-15)
        c. <- exp(pi*nu..*sgn*1i)
        Izm <- (c./z)[!P]
        stopifnot(allEQ(rI[ P,j  ]*rKp[,j+1] + rI[ P,j+1]*rKp[,j  ], Izp),
                  allEQ(rI[!P,j+1]*rKm[,j  ] - rI[!P,j  ]*rKm[,j+1], Izm) )
    }
}; cat("\n")


### Replicate some testing "ideas" from  zqcbk.f (TOMS 644 test program)

##  part 1) in the right half plane ----> see above (I & K)
##  part 2)
##      the analytic continuation formula
##     for H(nu,2,z) in terms of the K function
##
##           K(nu,z) = c3*H(nu,2,zr) + c4*H(nu,1,zr)    Im(z) >= 0
##                   =  conjg(K(nu,conjg(z)))           Im(z) < 0
##
##     in the left half plane where c3=c1*conjg(c2)*c5, c4 = c2*c5
##     c1=2*cos(pi*nu),   c2=exp(pi*nu*i/2),   c5 =-pi*i/2   and
##     zr = z*exp(-3*pi*i/2) = z*i


### Replicate some testing "ideas" from  zqcbj.f (TOMS 644 test program)

##     zqcbj is a quick check routine for the complex J bessel function
##     generated by subroutine zbesj.
##
##     zqcbj generates sequences of J and H bessel functions from zbesj
##     and zbesh and checks the wronskians
##
##     J(nu,z)*H(nu+1,1,z) - J(nu+1,z)*H(nu,1,z) =  2/(pi*i*z)   y >= 0
##     J(nu,z)*H(nu+1,2,z) - J(nu+1,z)*H(nu,2,z) = -2/(pi*i*z)   y < 0
##
##     in their respective half planes,  y = Im(z)

N <- 100
nS <- 20
set.seed(47)

for(i in seq_along(nus.)) {
    nu <- nus.[i]
    cat(nu, "")
    z <- complex(re = rnorm(N),
                 im = rnorm(N))
    P <- Im(z) >= 0
    rJ  <- BesselJ(  z,     nu, nSeq = 1+nS) # -> for (nu, nu+1, ...,nu+nS)
    rH1 <- BesselH(1,z[ P], nu, nSeq = 1+nS)
    rH2 <- BesselH(2,z[!P], nu, nSeq = 1+nS)
    ##
    sgn <- ifelse(Im(z) >= 0, +1, -1)
    Iz <- 2/(pi*1i*z)
    for(j in 1:nS) {
        nu.. <- nu + j-1
        allEQ <- function(x,y) all.equal(x,y, tol= max(1,nu..)*nS*1e-15)
        stopifnot(allEQ(rJ[ P,j]*rH1[,j+1] - rJ[ P,j+1]*rH1[,j],  Iz[ P]),
                  allEQ(rJ[!P,j]*rH2[,j+1] - rJ[!P,j+1]*rH2[,j], -Iz[!P]) )
    }
}; cat("\n")


### TODO __FIXME__

### Replicate some testing "ideas" from  zqcby.f (TOMS 644 test program)

##     zqcby generates sequences of y bessel functions from zbesy and
##     zbesyh and compares them for a variety of values in the (z,nu)
##     space. zbesyh is an old version of zbesy which computes the y
##     function from the h functions of kinds 1 and 2.

###--------> zbesyh() in ../src/zbsubs.c  is completely unneeded (otherwise) !


##  "limit  z -> 0  does not exist (there are many complex "Inf"s),
##		    but for z = real, z >=0 is -Inf
stopifnot(BesselY(0,1) == -Inf,# == besselY(0,1),
	  is.nan(BesselY(0+0i, 1)))


all.eq <- function(x,y, tol=1e-15,...) all.equal(x,y, tol=tol, ...)

## Subject: bug in Bessel package
## From: Hiroyuki Kawakatsu <...@gmail.com>, 18 Mar 2015

z <- c(0.23+1i, 1.21-1i)
nu <- -1/2
stopifnot(length(Bz.s <- BesselI(z, nu, expon.scaled=TRUE)) == length(z))
##
## Check that the exp() scaling is correct:
Bzu <- BesselI(z, nu)
stopifnot(abs(Im(sc <- Bz.s / Bzu)) < 1e-15,
	  all.eq(Re(sc), exp(-abs(Re(z)))))
## Using  nSeq > 1  -- and checking it:
options(warn=2)# warning = error {had warning of incompatible length}:
(Bz.s3 <- BesselI(z, nu, expon.scaled=TRUE, nSeq=3))
stopifnot(length(dim(Bz.s3)) == 2,
	  dim(Bz.s3) == c(length(z), 3),
	  all.eq(Bz.s3[,1], Bz.s),
	  all.eq(Bz.s3[,2], BesselI(z, nu-1, expon.scaled=TRUE)),
	  all.eq(Bz.s3[,3], BesselI(z, nu-2, expon.scaled=TRUE)))



cat('Time elapsed: ', proc.time(),'\n') # for ''statistical reasons''