Calculating integrals for a function approximated by Chebyshev polynomials

Question

Setup (complete, but all very standard):

My problem is how to best calculate the cumulative integral of a function which comes out of Spectral Collocation with a chebyshev basis. Take some function $f : [0, \bar{z}] \to \mathbb{R}$ approximated by $$ f(z) \approx \sum_{n=0}^{N-1}d_n T_n(z),\, z \in [0, \bar{z}] $$ where $T_n(z)$ is a basis of Chebyshev polynomials adapted to the $[0, \bar{z}]$ domain. Denote the vectors of coefficients as $d\in \mathbb{R}^N$. Calculate the Chebyshev polynomial roots (adapted to the $[0, \bar{z}]$ domain) and define them as, $$ \vec{z}_{\mathrm{int}} \equiv \{z_1,\ldots z_{N}\}\in \mathbb{R}^{N} $$

And the complete set of nodes including boundary values as (with $z_0 \equiv 0$ and $z_{N+1} \equiv \bar{z}$ as $\vec{z}\equiv \{0, z_1,\ldots z_N, \bar{z}\}\in \mathbb{R}^{N+2}$ Now, define the basis matrices as \begin{align} B &\equiv \begin{bmatrix} T_0(z_0)& \ldots & T_{N-1}(z_0)\\ \ldots & & \ldots\\ T_0(z_{N+1}) & \ldots & T_{N-1}(z_{N+1}) \end{bmatrix} \in \mathbb{R}^{(N+2)\times N}\\ B' &\equiv \begin{bmatrix} T_0'(z_0)& \ldots & T'_{N-1}(z_0)\\ \ldots & & \ldots\\ T_0'(z_{N+1}) & \ldots & T_{N-1}'(z_{N+1}) \end{bmatrix} \in \mathbb{R}^{(N+2)\times N} \end{align} Then, given the coefficient matrix $d$, you can find $f$ or the derivative $f'$ at every point in $\vec{z}$ with \begin{align} \vec{f} &\equiv \{f(z_n)\}_{n=0}^{N+1} = B \cdot d\in \mathbb{R}^{N+2}\\ \vec{f}' &\equiv \{f'(z_n)\}_{n=0}^{N+1} = B' \cdot d\in \mathbb{R}^{N+2}\\ \end{align} with the function at the interior nodes as $\vec{f}_{\mathrm{int}} \equiv \vec{f}(1:N) \in \mathbb{R}^N$.

Finally, we can find the weighting vector for Chebyshev-Gauss quadrature (https://en.wikipedia.org/wiki/Gaussian_quadrature) on the exact same roots $\vec{z}_{\mathrm{int}}$, and call it $\omega \in \mathbb{R}^N$. With this, we can approximate integrals for some $g(z)$ with $g(\vec{z}_{\mathrm{int}}) \equiv \{g(z) | z \in \vec{z}_{\mathrm{int}}\}$ $$ \int_0^{\bar{z}} g(z) d z \approx \omega \cdot g(\vec{z}_{\mathrm{int}}) $$ (note that this quadrature scheme does not use the endpoints).

My Problem: Define the cumulative function, $$ F(z) \equiv \int_0^z f(\tilde{z}) d\tilde{z} $$ If all I cared about was $F(\bar{z})$, then I have a nice approximation, $$ F(\bar{z}) \approx \omega \cdot \vec{f}_{\mathrm{int}} $$ But what an approximation for the cumulative integral, $F(z)$ for all $z \in \vec{z}$ given only the above?

Current Solution: Note that as I am not able to evaluate $f(z)$ at other points, I cannot naively use Gauss-Chebyshev quadrature as the quadrature roots depend on the domain of integration. My, grossly imprecise, method is to use trapezoidal integration at the unevenly spaced chebyshev roots. Or, $$ \Delta \vec{z} = \{\vec{z}(n) - \vec{z}(n-1)\}_{n=1}^{N+1} $$ Then, to calculate the integral to one of the nodes, $$ \int_0^{z_{n-1}}f(\hat{z}) d \hat{z} \approx \frac{1}{2}\Delta \vec{z}(0:n) \cdot (\vec{f}(0:n-1) + \vec{f}(1:n)) $$ I can even get fancier and calculate this all in one step,

\begin{align} \Omega \equiv \frac{1}{2}\begin{bmatrix}0 &0 & 0 &0 & \ldots & 0 \\ \Delta \vec{z}(0:1) & & 0 & 0 & \ldots & 0 \\ \Delta \vec{z}(0:2) & & & 0 & \ldots & 0 \\ \ldots & & & & & \ldots \\ \Delta \vec{z}(0:N)& & & & & \end{bmatrix}\in \mathbb{R}^{(N+2)\times (N+1)} \end{align}

Which gives the complete set of integrals as, $$ \vec{F} =\Omega \cdot (\vec{f}(0:N) + \vec{f}(1:N+1))\in\mathbb{R}^{N+2} $$ This is especially useful for me because I am solving a spectral collocation method, so having auto-differentiation within the calculation of the residual makes the problem solvable.

Are there better approaches? I am hoping that there is a more precise way of calculating these partial integrals. For example, I found https://math.stackexchange.com/questions/344073/integrating-non-uniform-grid-data-from-an-accelerometer which gives some ideas for a non-uniform Simpson's Rule, but it seems intractable to get a quadratic-form anything like my $\Omega$ above. Has anyone done that work, or is there a method I am forgetting about?

Answer 1

If you know the Chebyshev expansion for $f(z)$, why don't you formally integrate the polynomials using the recurrence relation for Chebyshev polynomials ? The Clenshaw-Curtis method is based on this approach (combined with an intelligent use of the FFT).

Answer 2

There seems to be a problem with your setup. Chebyshev polynomials are defined on the interval [-1,1], not on a general domain as you have specified.

So you need to scale your domain to [-1,1] to be able to use collocation.

The other problem of doing a general quadrature could simple be broken down to doing the quadrature over the domain [0,z] scaled to [-1,1].

On the other hand, if you want to fix the number of points, then you have to basically use some kind of further approximation. One way would be to break down the domain into smaller sub-domains, scale each one to [-1, 1], do the quadrature to the nearest point to z, and then do an approximate quadrature for the remaining part.

Answer 3

Based on the accepted answer from GertVdE I wrote the following Python routine that returns the quadrature weights for the Chebyshev nodes of second kind in the default range [-1,1]:

import numpy as np
from math import cos, pi

def cq(n,x):
    """Cumulative quadrature weights for a Chebyshev series

    Returns the weights w, such that the dot product w^T a is
    equal to the integral int_-1^x sum_{i=0}^n a_i T_i(x) dx
    """

    w = np.empty(n+1)
    w[0] = x + 1
    w[1] = (x**2 - 1)/2
    Tnm = 1
    Tn = x
    for i in range(2,n+1):
        Tn, Tnm = 2*x*Tn - Tnm, Tn
        w[i] = -(i*Tnm + (x - i*x)*Tn + cos(i*pi))/(i**2 - 1)
    return w

def cqb(n):
    """Cumulative quadrature matrix for a Chebyshev series

    Returns the (n+1) by (n+1) matrix B, which when multipled with
    the vector of coefficients of a Chebyshev series, returns a vector of 
    cumulative integrals evaluated at the Chebyshev nodes of second 
    kind: x_i = cos((n-i) pi/n), i = 0,...,n
    """

    if (n < 1):
        raise ValueError("The polynomial degree n should be >= 1")

    # we skip the knot at -1
    x = [cos((n - i)*pi/n) for i in range(1,n+1)]

    # at -1, the weights are zero
    B = np.zeros((n+1,n+1))
    for i, xi in enumerate(x):
        B[i+1,:] = cq(n,xi)
    return B

if __name__ == '__main__':
    
    B = cqb(3)
    print(B)