Adaptive Numerical Integration

Combines several approaches to adaptive numerical integration of functions of one variable.

integral(fun, xmin, xmax,
         method = c("Kronrod", "Clenshaw","Simpson"),
         no_intervals = 8, random = FALSE,
         reltol = 1e-8, abstol = 0, ...)

Arguments

fun: integrand, univariate (vectorized) function.
xmin,xmax: endpoints of the integration interval.
method: integration procedure, see below.
no_intervals: number of subdivisions at at start.
random: logical; shall the length of subdivisions be random.
reltol: relative tolerance.
abstol: absolute tolerance; not used.
...: additional parameters to be passed to the function.

Details

integral combines the following methods for adaptive numerical integration (also available as separate functions):

Kronrod (Gauss-Kronrod)
Clenshaw (Clenshaw-Curtis; not yet made adaptive)
Simpson (adaptive Simpson)

Recommended default method is Gauss-Kronrod. Also try Clenshaw-Curtis that may be faster at times.

Most methods require that function f is vectorized. This will be checked and the function vectorized if necessary.

By default, the integration domain is subdivided into no_intervals subdomains to avoid 0 results if the support of the integrand function is small compared to the whole domain. If random is true, nodes will be picked randomly, otherwise forming a regular division.

If the interval is infinite, quadinf will be called that accepts the same methods as well. [If the function is array-valued, quadv is called that applies an adaptive Simpson procedure, other methods are ignored – not true anymore.]

Value

Returns the integral, no error terms given.

References

Davis, Ph. J., and Ph. Rabinowitz (1984). Methods of Numerical Integration. Dover Publications, New York.

Note

integral does not provide `new' functionality, everything is already contained in the functions called here. Other interesting alternatives are Gauss-Richardson (quadgr) and Romberg (romberg) integration.

Examples

##  Very smooth function
fun <- function(x) 1/(x^4+x^2+0.9)
val <- 1.582232963729353
for (m in c("Kron", "Clen", "Simp")) {
    Q <- integral(fun, -1, 1, reltol = 1e-12, method = m)
    cat(m, Q, abs(Q-val), "\n")}
#> Kron 1.582233 3.197442e-13 
#> Clen 1.582233 3.199663e-13 
#> Simp 1.582233 3.241851e-13 
# Kron 1.582233 3.197442e-13 
# Rich 1.582233 3.197442e-13  # use quadgr()
# Clen 1.582233 3.199663e-13 
# Simp 1.582233 3.241851e-13 
# Romb 1.582233 2.555733e-13  # use romberg()

##  Highly oscillating function
fun <- function(x) sin(100*pi*x)/(pi*x)
val <- 0.4989868086930458
for (m in c("Kron", "Clen", "Simp")) {
    Q <- integral(fun, 0, 1, reltol = 1e-12, method = m)
    cat(m, Q, abs(Q-val), "\n")}
#> Kron 0.4989868 2.220446e-16 
#> Clen 0.4989868 2.259304e-14 
#> Simp 0.4989868 6.27276e-15 
# Kron 0.4989868 2.775558e-16 
# Rich 0.4989868 4.440892e-16  # use quadgr()
# Clen 0.4989868 2.231548e-14
# Simp 0.4989868 6.328271e-15 
# Romb 0.4989868 1.508793e-13  # use romberg()

## Evaluate improper integral
fun <- function(x) log(x)^2 * exp(-x^2)
val <- 1.9475221803007815976
Q <- integral(fun, 0, Inf, reltol = 1e-12)
#> For infinite domains Gauss integration is applied!
# For infinite domains Gauss integration is applied!
cat(m, Q, abs(Q-val), "\n")
#> Simp 1.947522 2.015876e-11 
# Kron 1.94752218028062 2.01587635473288e-11 

## Example with small function support
fun <- function(x)
            ifelse (x <= 0 | x >= pi, 0, sin(x))
integral(fun, -100, 100, no_intervals = 1)      # 0
#> [1] 0
integral(fun, -100, 100, no_intervals = 10)     # 1.99999999723
#> [1] 2
integral(fun, -100, 100, random=FALSE)          # 2
#> [1] 2
integral(fun, -100, 100, random=TRUE)           # 2 (sometimes 0 !)
#> [1] 0
integral(fun, -1000, 10000, random=FALSE)       # 0
#> [1] 0
integral(fun, -1000, 10000, random=TRUE)        # 0 (sometimes 2 !)
#> [1] 0