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?
Calculating integrals for a function approximated by Chebyshev polynomials
quadrature;collocation
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). 