pracma (version 1.4.5)

integral: Adaptive Numerical Integration

Description

Combines several approaches to adaptive numerical integration of functions of one variable, and provides complex line integrals.

Usage

integral(fun, xmin, xmax,
            method = c("Kronrod","Richardson","Clenshaw","Simpson","Romberg"),
            vectorized = TRUE, arrayValued = FALSE, waypoints = NULL,
            reltol = 1e-8, abstol = 0, ...)
cintegral(fun, waypoints, method = NULL, reltol = 1e-6, ...)

Arguments

fun

integrand, univariate (vectorized) function.

xmin,xmax

endpoints of the integration interval.

method

integration procedure, see below.

vectorized

logical; is the integrand a vectorized function; not used.

arrayValued

logical; is the integrand array-valued.

waypoints

complex integration: points on the integration curve.

reltol

relative tolerance.

abstol

absolute tolerance; not used.

...

additional parameters to be passed to the function.

Value

• Returns the integral, no error terms given.

Details

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

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

Most methods require that function f is vectorized.

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.

If one of the endpoints is complex, or waypoints is not NULL, function cintegral for complex line integration is called which accepts the same set of methods.

cintegral realizes complex line integration, in this case straight lines between the waypoints. By passing discrete points densely along the curve, arbitrary line integrals can be approximated.

Recommended default method is Gauss-Kronrod.

References

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

See Also

quadgk, quadgr, quadcc, simpadpt, romberg, quadv, quadinf

Examples

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

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

## Evaluate improper integral
fun <- function(x) log(x)^2 * exp(-x^2)
val <- 1.9475221803007815976
for (m in c("Kron", "Rich", "Clen", "Simp", "Romb")) {
    Q <- integral(fun, 0, Inf, reltol = 1e-12, method = m)
    cat(m, Q, abs(Q-val), "")}
# Kron 1.947522 1.101341e-13 
# Rich 1.947522 2.928655e-10 
# Clen 1.948016 1.960654e-13 
# Simp 1.947522 9.416912e-12 
# Romb 1.952683 0.00516102

## Array-valued function
log1 <- function(x) log(1+x)
fun <- function(x) c(exp(x), log1(x))
Q <- integral(fun, 0, 1, reltol = 1e-12, arrayValued = TRUE, method = "Simpson")
abs(Q - c(exp(1)-1, log(4)-1))          # 2.220446e-16 4.607426e-15

##  Complex integration examples
points <- c(0, 1+1i, 1-1i, 0)           # direction mathematically negative
f <- function(z) 1 / (2*z -1)
I <- cintegral(f, points)
abs(I - (0-pi*1i))                      # 3.788081e-13 ; residuum 2 pi 1i * 1/2

f <- function(z) 1/z
points <- c(-1i, 1, 1i, -1, -1i)
I <- cintegral(f, points)               # along a rectangle around 0+0i
abs(I - 2*pi*1i)                        #=> 0 ; residuum: 2 pi i * 1

N <- 100
x <- linspace(0, 2*pi, N)
y <- cos(x) + sin(x)*1i
J <- cintegral(f, waypoints = y)        # along a circle around 0+0i
abs(I - J)                              #=> 2.230056e-17; same residuum