Examples including integration

[ivan-pi]

For completeness, we need some examples that include

• integration of the PDE variables over the entire domain or individual elements/material domains, for $m = 0, 1, 2$
• calculation of cumulative integrals, at the knots of the Chebyshev basis (see this question: https://scicomp.stackexchange.com/q/23398/37438)

The coefficients in the Chebyshev expansion can be obtained from the work array wkres as follows (TOMS Algorithm 660, pg. 195):

C     
      NPTL = NPOLY + 1
      DO 100 I = 1,NEL
C      ITH ELEMENT
C      IU IS THE COMPONENT OF SOLUTION VECTOR AT LHS
C      OF ELEMENT I
       IU = (I - 1)*(NPOLY)*NPDE + 1
       DO 80 IS = 1,NPTL
        DO 80 JK = 1,NPDE
         COEFF(JK,IS,I) = 0.0D0
   80  CONTINUE
       DO 90 IS = 1,NPTL
        DO 90 JS = 1,NPTL
         DO 90 JK = 2,NPDE
         COEFF(JK,IS,I) = COEFF(JK,IS,I) +
     1                    WKRES(IS + (JS-1)*NPTL)*Y(IU + JS + JK - 2)
   90  CONTINUE
  100 CONTINUE
C

Or using free-form Fortran:

nptl = npoly + 1
do i = 1, nel
   ! ith element
   ! iu is the component of solution vector at LHS of element i
   iu = (i - 1)*npoly*npde + 1
   do is = 1, nptl
      do jk = 1, npde
         coeff(jk,is,i) = 0.0d0
      end do
   end do
   do is = 1, nptl
      do js = 1, nptl
         do jk = 2, npde
            coeff(jk,is,i) = coeff(jk,is,i) + &
               wkres(is + (js - 1)*nptl) * y(iu + js + jk - 2)
         end do
      end do
   end do
end do

The coefficient $a_{j,n}(t)$ of the $k$-th PDE is stored in COEFF(K,J,N), where $j$ runs along elements in the mesh, and $n$ counts the Chebyshev basis polynomials $T_n(x)$.

[ivan-pi]

Here is a draft code in Python for calculating the cumulative quadrature weights,

import scipy as sp
import numpy as np

# Based on the answer from the Computational Science Stack Exchange:
# https://scicomp.stackexchange.com/questions/23398/calculating-integrals-for-a-function-approximated-by-chebyshev-polynomials

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

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

    Arguments:
        n ... order of the Chebyshev polynomial
        x ... upper bound of the integral

    Returns:
        d ... a (n+1)-vector of coefficients
    """

    d = np.empty(n+1)
    
    d[0] = x + 1
    d[1] = x**2/2 - 1.0/2

    Tnm = 1.0
    Tn  = x

    Rnm, Rn = 1.0, -1.0

    for i in range(2,n+1):

        Tn, Tnm = 2*x*Tn - Tnm, Tn
        du = x*Tn/(i+1) - i*Tnm/(i**2 - 1)

        Rn, Rnm = Rnm, Rn
        dl = (-1.0)*Rn/(i+1) - i*Rnm/(i**2 - 1)

        d[i] = du - dl

    return d

def cquadmat(n):
    """Cumulative quadrature matrix at Chebyshev points

    Returns a matrix, which when multiplied with the coefficients of 
    a Chebyshev polynomials, returns a vector of integrals at the
    Chebyshev points (Chebyshev extrema or roots of the Chebyshev polynomial
    of second kind)
    """

    assert n >= 1

    B = np.zeros((n+1,n+1))
    for i in range(0, n+1):
        xi = np.cos((n - i)*np.pi/n)
        B[i,:] = cqw(n,xi)
    return B    

if __name__ == '__main__':
    
    B = cquadmat(3)

    print(B)