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,
           max.ord = 19, verbose = FALSE)

Arguments

an R function 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. Do use "mpfr"-numbers to get higher than double precision, see the examples.

…

additional arguments to be passed to f.

ord

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

rel.tol

relative accuracy requested. The default is 1.2e-4, about 4 digits only, see the Note.

abs.tol

absolute accuracy requested.

max.ord

only used, when neither ord or one of rel.tol, abs.tol are specified: Stop Romberg iterations after the order reaches max.ord; may prevent infinite loops or memory explosion.

verbose

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

Value

A list of class "integrateR" (as from standard R's integrate()) with a print method and 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.

For convergence, both relative and absolute changes must be smaller than rel.tol and abs.tol, respectively.

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

References

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

Examples

Run this code

# NOT RUN {
## 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)
(Id <- integrateR(dnorm,0,2000, rel.tol=1e-15, verbose=TRUE))
Id$value == 0.5 # exactly

## 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
## Now, using higher accuracy:
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)

## Want absolute tolerance check only (=> set 'rel.tol' very high, e.g. 1):
(Ia <- integrateR(exp, mpfr(0,200), 1, abs.tol = 1e-6, rel.tol=1, verbose=TRUE))

## Set 'ord' (but no  '*.tol') --> Using 'ord'; no convergence checking
(I <- integrateR(exp, mpfr(0,200), 1,  ord = 13, verbose=TRUE))

# }

Run the code above in your browser using DataLab

Get 50% off unlimited learning