https://github.com/cran/pracma
Tip revision: 85eeff77b2e878cc9dbd7659e4acbc035be93c28 authored by Hans W. Borchers on 22 September 2022, 13:50:02 UTC
version 2.4.2
version 2.4.2
Tip revision: 85eeff7
clenshaw_curtis.R
##
## q u a d c c . R Clenshaw-Curtis Quadrature
##
clenshaw_curtis <- function(f, a = -1, b = 1, n = 1024, ...) {
fun <- match.fun(f)
f <- function(x) fun(x, ...)
if (a == b) return(0)
eps <- .Machine$double.eps # assume a < b
if (!is.finite(f(a))) a <- a + eps * sign(b-a)
if (!is.finite(f(b))) b <- b - eps * sign(b-a)
# Evaluate f at Chebyshev points
x <- cos(pi*(0:n)/n)
fx <- f(0.5*((b-a)*x + (b+a)))/(2*n)
# Fast Fourier transform
g <- Re(fft(fx[c(1:(n+1), n:2)]))
# Chebyshev coefficients and weight vector
d <- c(g[1], g[2:n] + g[(2*n):(n+2)], g[n+1])
w <- 0 * d
w[seq(1, n+1, by=2)] <- 2/(1-(seq(0, n, by=2))^2)
# Return the integral
Q <- sum(w * d) * (b-a)/2
return(Q)
}