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
quadcc.R
##
## q u a d c c . R Adaptive Clenshaw-Curtis Quadrature
##
.quadcc <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {
stopifnot(is.numeric(a), is.numeric(b),
length(a) == 1, length(b) == 1)
fun <- match.fun(f)
f <- function(x) fun(x, ...)
eps <- .Machine$double.eps
if (a == b) {
return(list(value = 0, error = 0))
} else if (a > b) {
Q <- quadcc(f, b, a, tol = tol)
return(list(value = -Q$value, error = Q$error))
}
# Break integration interval into four regions
N <- 4
p <- (b+a)/2
q <- (b-a)/2
# Write out [ cos((0:4)/4*pi) ] explicitly to avoid cos(pi/2) roundoff
x <- p - q * c(1, 0.7071067811865475, 0, -0.7071067811865475, -1)
F <- f(x)
if (!is.finite(F[1])) F[1] <- f(x[1] + eps)
if (!is.finite(F[5])) F[5] <- f(x[5] - eps)
# Modify the first and last term of F
F[1] <- F[1]/2
F[5] <- F[5]/2
y <- 0; yinit <- Inf
while (abs(y - yinit) > tol && N < 2^13) {
yinit <- y
# Double the size of the approximation
N <- 2*N
Nvec <- 0:N
x2 <- numeric(N+1)
s1 <- seq(1, N+1, by = 2)
s2 <- seq(2, N+1, by = 2)
x2[s1] <- x
x2[s2] <- p - q * cos(seq(1,N,by=2)/N * pi)
G <- numeric(N+1)
G[s1] <- F
G[s2] <- f(x2[s2])
if (any(is.infinite(G))) {
warning("Integrand has a pole between integration limits.")
y <- NA
break
}
# Evaluate the Chebyshev polynomial coefficients
M <- cos(as.matrix(Nvec[seq(1, length(Nvec)-1, by=2)]) %*% Nvec/N*pi)
M[abs(M) < eps] <- 0
C <- 4/N * M %*% G
# Compute the quadrature coefficients
coeffs <- -1/Nvec[s2]/(Nvec[s2]-2)
coeffs[1] <- 0.5
# Compute the integral estimate
y <- q * sum(coeffs * C)
# Reset the variables
x <- x2
F <- G
}
return(list(value = y, error = abs(y - yinit)))
}