Rmpfr (version 0.5-3)

integrateR: One-Dimensional Numerical Integration - in pure R

Description

Numerical integration of one-dimensional functions in pure R, with care so it also works for "mpfr"-numbers.

Currently, only classical Romberg integration of order ord is available.

Usage

integrateR(f, lower, upper, ..., ord = NULL,
           rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol,
           verbose = FALSE)

Arguments

f

an Rfunction taking a numeric or "mpfr" first argument and returning a numeric (or "mpfr") vector of the same length. Returning a non-finite element will generate an error.

lower, upper

the limits of integration. Currently must be finite.

...

additional arguments to be passed to f.

ord

integer, the order of Romberg integration to be used. If this is NULL, as per default, the order is increased until convergence, see rel.tol and abs.tol.

rel.tol

relative accuracy requested. The default is 1.2e-4.

abs.tol

absolute accuracy requested.

verbose

logical or integer, indicating if and how much information should be printed during computation.

Value

• A list of class "integrate" with components
• value
  the final estimate of the integral.
• abs.error
  estimate of the modulus of the absolute error.
• subdivisions
  for Romberg, the number of function evaluations.
• message
  "OK" or a character string giving the error message.
• call
  the matched call.

Details

Note that arguments after ... must be matched exactly.

rel.tol cannot be less than max(50*.Machine$double.eps, 0.5e-28) if abs.tol <= 0<="" code="">.

References

Bauer, F.L. (1961) Algorithm 60 -- Romberg Integration, Communications of the ACM 4(6), p.255.

See Also

R's standard, integrate, is much more adaptive, also allowing infinite integration boundaries, and typically considerably faster for a given accuracy.

We use the same (actually a copy) print S3 method, print.integrate(), as provided by R.

Examples

## See more  from  ?integrate
## this is in the region where  integrate() can get problems:
integrateR(dnorm,0,2000)
integrateR(dnorm,0,2000, rel.tol=1e-15)
integrateR(dnorm,0,2000, rel.tol=1e-15, verbose=TRUE)

## Demonstrating that 'subdivisions' is correct:
Exp <- function(x) { .N <<- .N+ length(x); exp(x) }
.N <- 0; str(integrateR(Exp, 0,1, rel.tol=1e-10), digits=15); .N

### Using high-precision functions -----

## Polynomials are very nice:
integrateR(function(x) (x-2)^4 - 3*(x-3)^2, 0, 5, verbose=TRUE)
# n= 1, 2^n=        2 | I =            46.04, abs.err =      98.9583
# n= 2, 2^n=        4 | I =               20, abs.err =      26.0417
# n= 3, 2^n=        8 | I =               20, abs.err =  7.10543e-15
## 20 with absolute error < 7.1e-15
I <- integrateR(function(x) (x-2)^4 - 3*(x-3)^2, 0, mpfr(5,128),
                rel.tol = 1e-20, verbose=TRUE)
I ; I$value  ## all fine

## with floats:
integrateR(exp,      0     , 1, rel.tol=1e-15, verbose=TRUE)
## with "mpfr":
(I <- integrateR(exp, mpfr(0,200), 1, rel.tol=1e-25, verbose=TRUE))
(I.true <- exp(mpfr(1, 200)) - 1)
## true absolute error:
stopifnot(print(as.numeric(I.true - I$value)) < 4e-25)