# R/integrate-Romberg.R In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

#### Documented in integrateRprint.integrateR

```####  Romberg integration in pure R
####  ===================    ====== so it can be used with  Rmpfr

## TODO: Lauren K  would like to get return()ed all the intermediate sums as well
## ----  I agree: but only if  'all.sums = TRUE'

integrateR <- function(f, lower, upper, ..., ord = NULL,
rel.tol = .Machine\$double.eps^0.25, abs.tol = rel.tol,
max.ord = 19, verbose = FALSE)
{
stopifnot(length(lower) == 1, length(upper) == 1,
is.finite(lower), is.finite(upper))
f <- match.fun(f)
ff <-
## if(verbose) function(x) { cat("f(x), x="); str(x) ; f(x, ...) } else
function(x) f(x, ...)
null.ord <- is.null(ord) ## ord := Romberg order
has.tol <- !missing(rel.tol) || !missing(abs.tol)# one of them specified
if(chkConv <- (null.ord || has.tol)) {
## will use rel.tol and abs.tol
if (abs.tol <= 0 && rel.tol < max(50 * .Machine\$double.eps, 5e-29))
stop("invalid tolerance values")
## but need (maximal) order for Bauer's algorithm:
}
if(null.ord && !has.tol) {
## This is "approximate" (and too large, typically); but if it's
## too small, t[.] will be extended automatically:
ord <- ## == max(3, min(25, ceiling(-log2(rel.tol))))  with default rel.tol
13
if(verbose) cat("ord =", ord, "(as no '*.tol' specified)\n")
}
useOrd <- !null.ord || !has.tol
if(useOrd) {
stopifnot(ord >= 0)
if(verbose) cat(sprintf(
" ord = %d; ==> evaluating integrand at %s 2^(ord+1)-2 = %d locations\n",
ord, if(chkConv) "up to" else "", 2^(ord+1)-2))
}

### Bauer(1961)  "Algorithm 60 -- Romberg Integration" Comm.ACM 4(6), p.255
m <- le <- upper - lower # 'l'
##  a "hack", but really improves the result:
if(!is.numeric(m)) {
if(is(m, "mpfr")) { # should get *same* precision
if(is.numeric(lower)) lower <- 0*m+ lower
if(is.numeric(upper)) upper <- 0*m+ upper
} else { ## other high-precision...
if(is.numeric(lower)) lower <- as(lower, class(m))
if(is.numeric(upper)) upper <- as(upper, class(m))
}
}
t1 <- (ff(lower) + ff(upper))/2
t <- rep(t1, if(useOrd) ord+1 else 10)## <- must work for "mpfr" numbers
one <- 1 + 0*t1 # for "mpfr"

r. <- t[1]*le
n <- 1 # 'n'(Bauer) = 2^n (Romberg Algo)
if(verbose) { ## In "mpfr" case, cannot use sprintf("%g");
## ==> rather use format(.) with higher number of digits
prDigs <- max(10, min(50, 2 + ceiling(-log10(rel.tol))))
FORM <- paste0("n=%2d, 2^n=%9.0f | I = %",(5+prDigs), "s, abs.err =%14s\n")
}
h <- 1L
repeat {
if(verbose >= 2) {
cat("range(le*t[1:h]):\n\t");
print(format(range(le*t[1:h]), digits=15), quote=FALSE)
}
u <- 0
m <- m/2 # == le/(2*n)
## here, we require f(.) to be vectorized:
u <- sum(ff(lower+ seq(1, 2*n-1, by=2)*m))
t[h+1L] <- (u/n + t[h])/2
f. <- one
for(j in h:1) {
f. <- 4*f.
t[j] <- t[j+1] + (t[j+1] - t[j])/ (f. - 1)
}
r <- t[1]*le
aErr <- abs(r - r.)
if(verbose)
cat(sprintf(FORM, h, 2*n, format(r, digits = prDigs),
format(aErr, digits = max(7, getOption("digits")))))
if(chkConv) { ## check if we converged: |r - r.| < min(.,.): *both* tolerances must be satisfied
if(converged <- (aErr < min(abs(r)*rel.tol, abs.tol)))
break
}
if((useOrd && h >= ord) || h >= max.ord) break
r. <- r
n <- 2*n # == 2^h
h <- h+1L
}
if(chkConv && !converged) {
relE <- format(aErr/abs(r), digits=getOption("digits"))
msg <- paste0("no convergence up to order ", ord,
"; last relative change = ", relE, "\n",
"Consider setting 'ord = <n>' (e.g. = ", ord+1,").")
warning(msg)
} else msg <- "OK"
r <- list(value = r, abs.error = aErr, subdivisions = 2*n+1,
"message" = msg, call = match.call())
class(r) <- c("integrateR", "integrate")
r
}

## This is  such that  print.integrate() calls our format(<mpfr>) method
## (and do *not* hide it via S3method() in NAMESPACE):
## print.integrate <- getS3method("print","integrate")# from 'stats' possibly not exported
## environment(print.integrate) <- environment()
## setMethod(show, "integrate", function(object) print.integrate(object))

print.integrateR <- function (x, digits = max(3, getOption("digits")-2), ...)
{
if(x[["message"]] != "OK")
cat("Non-convergence message ", sQuote(x\$message), "\n", sep = "")
## The "Ok" message:
cat(format(x\$value, digits = digits),
" with absolute error < ", format(x\$abs.error, digits=digits),
"\n", sep = "")
invisible(x)
}
setMethod(show, "integrateR", function(object) print.integrateR(object))
```

## Try the Rmpfr package in your browser

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

Rmpfr documentation built on Jan. 13, 2018, 1:11 a.m.