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 http://math.stackexchange.com/questions/344073/integrating-non-uniform-grid-data-from-an-accelerometerhttps://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?