# integrateR: One-Dimensional Numerical Integration - in pure R In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

## 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

 ```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) ```

## Arguments

 `f` 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.

## 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`.

## 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.

## Note

`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()`.

## Note

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

## References

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=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)) ```