Integral approximation using a matrix operator

Question

In the paper "Chebyshev solution of differential, integral and integro-differential equations" (it is freely accessible and can be downloaded from the link), El-gendi uses a method to approximate the integral

$$\int_{-1}^x f(t) \,dt$$

by a Chebyshev approximation at points

$$x_j = \cos\left(\frac{j \pi}{N}\right)$$

by a matrix $B_{(N+1)\times (N+1)}$:

$$\left[\int_{-1}^x f(t) \,dt\right]=B.[f]$$

where $B$ is a function of $N$ only and that's what I am trying to find ultimately. In other words, I want to be able to produce this matrix for different values of $N$ . But I can't reproduce the example $B$ that the paper has reported for $N=4$ in Mathematica.

Here is what I have so far:

xp[ip_, np_] := Cos[ip Pi/np]

c01[rr_, np_] := If[rr == 0 || rr == np, 0.5, 1]

a[r_, np_] := 
 2/np Sum[c01[k, np] f[xp[k, np]] ChebyshevT[r, xp[k, np]], {k, 0, np}]

c[kk_, nn_] := 
 Which[kk == 0, 
  0.5 a[0, nn] + 
   Sum[c01[jj, nn] ((-1^(jj + 1) a[jj, nn])/(jj^2 - 1)), {jj, 2, 
     nn}] - 0.25 a[1, nn], kk >= 1 && kk <= nn - 2, (
  a[kk - 1, nn] - a[kk + 1, nn])/(2 kk), kk == nn - 1, (
  a[nn - 2, nn] - 0.5 a[nn, nn])/(2 (nn - 1)), kk == nn, 
  a[nn - 1, nn]/(2 (nn)), kk == nn + 1, (0.5 a[nn, nn])/(2 (nn + 1))]

Now when I produce the approximation to the integral for $N=4$ by:

ct4 = Sum[c[r, 4] ChebyshevT[r, x], {r, 0, 5}] // Simplify

and then by replacing $x$ values:

ct4 /. x -> {xp[0, 4], xp[1, 4], xp[2, 4], xp[3, 4], xp[4, 4]}

I get:

{0.129167 f[-1] + 0.8 f[0] + 0.00416667 f[1] + 
  0.444945 f[-(1/Sqrt[2])] + 0.621722 f[1/Sqrt[2]], 
 0.134763 f[-1] + 0.824264 f[0] - 0.115237 f[1] + 
  0.431658 f[-(1/Sqrt[2])] + 0.431658 f[1/Sqrt[2]], 
 0. + 0.0958333 f[-1] + 0.4 f[0] - 0.0291667 f[1] + 
  0.531832 f[-(1/Sqrt[2])] + 0.00150162 f[1/Sqrt[2]], 
 0.181904 f[-1] - 0.0242641 f[0] - 0.0680964 f[1] + 
  0.101675 f[-(1/Sqrt[2])] + 0.101675 f[1/Sqrt[2]], 
 0. + 0.0625 f[-1] - 0.0625 f[1] - 0.0883883 f[-(1/Sqrt[2])] + 
  0.0883883 f[1/Sqrt[2]]}

which is different from the results in the paper:

Any ideas what I'm doing wrong? I'd appreciate it if someone has a better way to derive this table too.

Answer 1

Here is the integration analog of chebmat[] from this answer. It is not very efficiently written, but at least the matrix is assembled from recognizable pieces:

chebint[n_Integer?Positive] := 
 Module[{cm = Cos[π Array[Times, {n + 1, n + 2}, {0, 0}]/n], v = 1/(2 Range[n])},
        Reverse[cm].
        SparseArray[{{1, j_} :> Switch[j - 1, 0, 1/2, 1, -1/4,
                                       n, (-1)^(n + 1)/(2 (n^2 - 1)),
                                       _, (-1)^j/(j (j - 2))], 
                     Band[{2, 1}] -> Append[v, 1/(4 (n + 1))], 
                     Band[{2, 3}] -> -Append[Drop[v, -2], 1/(4 (n - 1))]},
                    {n + 2, n + 1}].
        Drop[Reverse[cm, 2], None, 1].
        DiagonalMatrix[ArrayPad[ConstantArray[2/n, n - 1], 1, 1/n]]]

Test:

n = 4;
cnodes = N[-Cos[π Range[0, n]/n]];
fc = Exp[cnodes];

chebint[4].fc
   {0., 0.12511070068305466, 0.6324668678294361, 1.660103967519322, 2.350375376931479}

Function[x, Exp[x] - 1/E] /@ cnodes
   {0., 0.1251892502237975, 0.6321205588285577, 1.6602355404760298, 2.3504023872876028}

A more efficient implementation would likely involve the use of FourierDCT[]...