This is an implementation of the Clenshaw-Curtis quadrature for solving definite integrals. It is a direct translation from the MATLAB code presented in (Trefethen, L. N. (2008). Is Gauss quadrature better than Clenshaw–Curtis?. SIAM review, 50(1), 67-87).

clenshaw_curtis.py

import numpy as np
from scipy.fft import fft


def clenshaw_curtis(f, a, b, n):
    # Obtain Chebyshev points
    x = np.cos(np.pi * np.arange(n + 1) / n)
    # Rescale the points to the interval
    scaled_x = ((a + b) + ((b - a) * x)) * 0.5
    # Evalute the integrand
    fx = f(scaled_x) / (2.0 * n)
    # We need symmetric values for the FFT
    indices = np.concatenate((np.arange(n), np.arange(n, 0, -1)))
    g = np.real(fft(fx[indices]))
    # Now, compute the weights
    aug_g = np.concatenate(
        (
            np.array(g[0]).reshape(-1),
            g[1:n] + g[2 * (n + 1) : n : -1],
            np.array(g[n]).reshape(-1),
        )
    )
    w = np.zeros_like(aug_g)
    w[::2] = 2.0 / (1.0 - np.arange(0, n + 1, 2) ** 2)
    # Finally, rescale the weights to the integration interval
    w *= (b - a) * 0.5

    return w @ aug_g


def main():
    # The result should be around 44.3365294622606835310
    print(clenshaw_curtis(lambda x: x**2 + np.sin(x), -2, 5, 11))


if __name__ == "__main__":
    main()