3
$\begingroup$

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:

correct matrix

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

$\endgroup$

1 Answer 1

5
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ Amazing! You never disappoint! $\endgroup$
    – MathX
    Jan 13, 2017 at 16:29

Your Answer

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

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