Description Usage Arguments Details Value Note Note Author(s) References See Also Examples
View source: R/integrateRomberg.R
Numerical integration of onedimensional functions in pure R,
with care so it also works for "mpfr"
numbers.
Currently, only classical Romberg integration of order ord
is
available.
1 2 3  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.2e4, 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.5e28)
if abs.tol <= 0
.
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 

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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37  ## See more from ?integrate
## this is in the region where integrate() can get problems:
integrateR(dnorm,0,2000)
integrateR(dnorm,0,2000, rel.tol=1e15)
(Id < integrateR(dnorm,0,2000, rel.tol=1e15, 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=1e10), digits=15); .N
### Using highprecision functions 
## Polynomials are very nice:
integrateR(function(x) (x2)^4  3*(x3)^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.10543e15
## 20 with absolute error < 7.1e15
## Now, using higher accuracy:
I < integrateR(function(x) (x2)^4  3*(x3)^2, 0, mpfr(5,128),
rel.tol = 1e20, verbose=TRUE)
I ; I$value ## all fine
## with floats:
integrateR(exp, 0 , 1, rel.tol=1e15, verbose=TRUE)
## with "mpfr":
(I < integrateR(exp, mpfr(0,200), 1, rel.tol=1e25, verbose=TRUE))
(I.true < exp(mpfr(1, 200))  1)
## true absolute error:
stopifnot(print(as.numeric(I.true  I$value)) < 4e25)
## Want absolute tolerance check only (=> set 'rel.tol' very high, e.g. 1):
(Ia < integrateR(exp, mpfr(0,200), 1, abs.tol = 1e6, 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))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.