One-Dimensional Numerical Integration - in pure R
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.
integrateR(f, lower, upper, ..., ord = NULL,
rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol,
max.ord = 19, verbose = FALSE)f |
an R function taking a numeric or |
lower, upper |
the limits of integration. Currently must
be finite. Do use |
... |
additional arguments to be passed to |
ord |
integer, the order of Romberg integration to be used. If
this is |
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 |
verbose |
logical or integer, indicating if and how much information should be printed during computation. |
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.
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 |
|
call |
the matched call. |
f must accept a vector of inputs and produce a vector of function
evaluations at those points. The Vectorize function
may be helpful to convert f to this form.
If you want to use higher accuracy, you must set lower or
upper to "mpfr" numbers (and typically lower the
relative tolerance, rel.tol), see also the examples.
Note that the default tolerances (rel.tol, abs.tol) are
not very accurate, but the same as for integrate, which
however often returns considerably more accurate results than
requested. This is typically not the case for
integrateR().
We use practically the same print S3 method as
print.integrate, provided by R,
with a difference when the message component is not "Ok".
Martin Maechler
Bauer, F.L. (1961) Algorithm 60 – Romberg Integration, Communications of the ACM 4(6), p.255.
R's standard, integrate, is much more adaptive,
also allowing infinite integration boundaries, and typically
considerably faster for a given accuracy.
## 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))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.