https://github.com/cran/pracma
Tip revision: 708a2ad382a163d1eef5af0665e3ae2aad200ced authored by HwB on 21 March 2013, 00:00:00 UTC
version 1.4.5
version 1.4.5
Tip revision: 708a2ad
integral.R
##
## i n t e g r a l . R Numerical Integration
##
integral <- function(fun, xmin, xmax,
method = c("Kronrod","Richardson","Clenshaw","Simpson","Romberg"),
vectorized = TRUE, arrayValued = FALSE, waypoints = NULL,
reltol = 1e-8, abstol = 0, ...)
{
stopifnot(is.numeric(xmin) || is.complex(xmin), length(xmin) == 1,
is.numeric(xmax) || is.complex(xmax), length(xmax) == 1)
fun <- match.fun(fun)
f <- function(x) fun(x, ...)
# if (!vectorized(f)) # check with those routines called
# f <- Vectorize(f)
if (xmin == xmax) return(0)
method <- match.arg(method)
tol <- reltol # abstol not used
if (is.complex(xmin) || is.complex(xmax) || !is.null(waypoints)) {
Q <- cintegral(f, waypoints = as.complex(c(xmin, waypoints, xmax)),
method = method, tol = reltol)
} else if (arrayValued) {
if (method != "Simpson")
warning("Only method 'Simpson' available for array-valued functions.")
Q <- quadv(f, xmin, xmax, tol = tol)$Q
} else if (is.infinite(xmin) || is.infinite(xmax)) {
Q <- quadinf(f, xmin, xmax, tol = tol, method = method)
} else {
Q <- switch(method,
"Kronrod" = quadgk(f, xmin, xmax, tol = tol),
"Clenshaw" = quadcc(f, xmin, xmax, tol = tol),
"Richardson" = quadgr(f, xmin, xmax, tol = tol)$value,
"Romberg" = romberg(f, xmin, xmax, tol = tol)$value,
"Simpson" = simpadpt(f, xmin, xmax, tol = tol)
)
}
Q <- Q
return(Q)
}
cintegral <- function (fun, waypoints, method = NULL, reltol = 1e-6, ...) {
stopifnot(is.complex(points) || is.numeric(points))
if (length(points) <= 1) return(0 + 0i)
fun <- match.fun(fun)
f <- function(z) fun(z, ...)
Q <- 0 + 0i
for (i in 2:length(waypoints)) {
a <- waypoints[i-1]
b <- waypoints[i]
d <- b - a
f1 <- function(t) Re(f(a + t*d))
f2 <- function(t) Im(f(a + t*d))
Qre <- integral(f1, 0, 1, reltol = reltol, method = method)
Qim <- integral(f2, 0, 1, reltol = reltol, method = method)
Q <- Q + d * (Qre + Qim*1i)
}
return(Q)
}