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?


3 Answers 3


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).

  • $\begingroup$ Thank you so much. So are you saying that I am overthinking this, and all I need to do is find a new basis matrix above for the integra like my $B$ and $B'$l, say $B^{int}$, where $B^{int}[i,j] \equiv \int_0^{z_i}T_j(z) dz$? Then my cumulative integrals is just $\vec{F} = B^{int} \cdot d$? $\endgroup$
    – jlperla
    Commented Mar 20, 2016 at 14:06
  • $\begingroup$ I was able to build a basis matrix, as described, given the standard recurrence formula $\int T_n(x) dx = \frac{1}{2}\left( \frac{T_{n+1}(x)}{n+1} - \frac{T_{n-1}(x)}{n-1} \right)$ with special cases for $n=0,1$ and with the affine transformation from my domain. Thanks for the suggestion. $\endgroup$
    – jlperla
    Commented Mar 25, 2016 at 6:26

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.

  • $\begingroup$ Sorry, I meant to say that the nodes and the polynomials were rescaled to my domain. I added that comment above. As for breaking down the quadrature for a particular $[0,z]$ other than $\bar{z}$, the problem is that the quadrature nodes would no longer line up. I can only calculate my function at the $[0, \bar{z}]$ nodes. Thanks for your comments $\endgroup$
    – jlperla
    Commented Mar 20, 2016 at 14:01

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.