Skip to main content

Browse the archive

https://github.com/cran/pracma
21 November 2023, 19:34:12 UTC
Revision 26e049d70b4a1c237987e260cba68f6a9413736c authored by Hans W. Borchers on 09 April 2019, 04:10:07 UTC, committed by cran-robot on 09 April 2019, 04:10:07 UTC 
version 2.2.5
1 parent bf07673
Raw File
Tip revision: 26e049d70b4a1c237987e260cba68f6a9413736c authored by Hans W. Borchers on 09 April 2019, 04:10:07 UTC
version 2.2.5
Tip revision: 26e049d
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)
}
back to top