# Integral approximation using a matrix operator

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.00416667 f +
0.444945 f[-(1/Sqrt)] + 0.621722 f[1/Sqrt],
0.134763 f[-1] + 0.824264 f - 0.115237 f +
0.431658 f[-(1/Sqrt)] + 0.431658 f[1/Sqrt],
0. + 0.0958333 f[-1] + 0.4 f - 0.0291667 f +
0.531832 f[-(1/Sqrt)] + 0.00150162 f[1/Sqrt],
0.181904 f[-1] - 0.0242641 f - 0.0680964 f +
0.101675 f[-(1/Sqrt)] + 0.101675 f[1/Sqrt],
0. + 0.0625 f[-1] - 0.0625 f - 0.0883883 f[-(1/Sqrt)] +
0.0883883 f[1/Sqrt]}


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.

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.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[]...

• Amazing! You never disappoint! – MathX Jan 13 '17 at 16:29