Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to get approximate the derivative of a function from its Chebyshev expansion.

I start out with the following random function

In[4]:= f[x_] = -4.1261 + .99213*x - 0.562172*x^2 - 3.7120 x^3 + 3.06508*x^4

Which looks like graph of random function

Now I try the Chebyshev approximations as follows:

In[105]:= cgs = 5

In[107]:= Clear[A, a]

In[108]:= A = Array[a, cgs + 1, {0}]

Out[108]= {a[0], a[1], a[2], a[3], a[4], a[5]}

In[109]:=   a[0] = (1/Pi)* Integrate[(f[x]* ChebyshevT[0, x]/Sqrt[1-   
x^2]), {x, -1, 1}]

Out[109]= -3.25778

In[110]:= For[i = 1, i <= cgs, i++,  a[i] = (2/Pi)*Integrate[(f[x]* 
ChebyshevT[i, x]/Sqrt[1 - x^2]), {x, -1, 1}]]

In[111]:= A

Out[111]= {-3.25778, -1.79187, 1.25145, -0.928, 0.383135, -1.97901*10^-15}

In[113]:= cf = Table[a[i]*ChebyshevT[i, x], {i, 0, cgs}]

Out[113]= {-3.25778, -1.79187 x, 1.25145 (-1 + 2 x^2), -0.928 (-3 x + 4     
x^3), 0.383135 (1 - 8 x^2 + 8 x^4), -1.97901*10^-15 (5 x - 20 x^3 + 16 x^5)}

In[115]:= CF[x_] = Sum[cf[[i]], {i, 1, cgs + 1}]

Out[115]= -3.25778 - 1.79187 x + 1.25145 (-1 + 2 x^2) - 0.928 (-3 x + 4 x^3)     
+ 0.383135 (1 - 8 x^2 + 8 x^4) -  1.97901*10^-15 (5 x - 20 x^3 + 16 x^5)

Plot[{CF[x], f[x]}, {x, -1, 1}]

gives

graph of function and Chebyshev expansion

So I think I have done everything correctly up until now, because I cannot see a difference between the two graphs.

As for the derivative approximation, I use the discussion in Numerical Recipes section 5.9.

First I take the derivative of f[x]

In[200]:= df[x_] = D[f[x], x]

Out[200]= 0.99213 - 1.12434 x - 11.136 x^2 + 12.2603 x^3

In[201]:= A

Out[201]= {-3.25778, -1.79187, 1.25145, -0.928, 0.383135, -1.97901*10^-15}

In[236]:= Clear[dA, da]

In[237]:= dA = Array[da, cgs + 1, {0}]

Out[237]= {da[0], da[1], da[2], da[3], da[4], da[5]}

Now I use equation (5.9.2)

In[273]:= da[cgs] = da[cgs - 1] = 0;

In[275]:= For[i = 4, i >= 1, i--, da[i - 1] = 2*(i)*a[i] + da[i + 1]]

to get

In[276]:= dA

Out[276]= {-9.15174, 8.0709, -5.568, 3.06508, 0, 0}

In[277]:= dcf = Table[da[i]*ChebyshevT[i, x], {i, 0, cgs}]

Out[277]= {-9.15174, 8.0709 x, -5.568 (-1 + 2 x^2), 3.06508 (-3 x + 4 x^3),       
0, 0}

In[278]:= dCF[x_] = Sum[dcf[[i]], {i, 1, cgs + 1}]

Out[278]= -9.15174 + 8.0709 x - 5.568 (-1 + 2 x^2) + 3.06508 (-3 x + 4 x^3)

Plot[{dCF[x], df[x]}, {x, -1, 1}]

which gives

plot of wrong derivative

which is clearly not the same graph. I was able to figure out that the problem is in the constant term. Namely, it is twice the value it should be. To test this idea I created the test function

In[279]:=  test[x_] = -9.151740000000007`/2 + 8.070896000000015` x - 
5.568000000000007` (-1 + 2 x^2) + 3.065080000000011` (-3 x + 4 x^3);

and plotted Plot[{test[x], df[x]}, {x, -1, 1}] to get

plot of derivative

which shows that test is the same as f[x].

If anyone could help spot my mistake, that would be greatly appreciated.

share|improve this question
    
Do not use For[] function in Mathematica, use Table[] instead. In this case it will make Your code more readable. –  mmal Oct 25 '12 at 23:47
add comment

1 Answer 1

Yes, I believe the NR authors neglected to halve the last coefficient in their formulae (the halving is due to the fact that $T_1(x)=x$, not $2x$). As a check on the correctness of your coefficients, I'll use (a special case of) Salzer's algorithm to produce Chebyshev coefficients from the monomial coefficients:

c = CoefficientList[-4.1261 + .99213 x - 0.562172 x^2 - 3.7120 x^3 + 3.06508 x^4, x];
n = Length[c] - 1; Remove[a];
a[0, 2] = c[[n - 1]] + c[[n + 1]]/2;
a[1, 2] = c[[n]]; a[2, 2] = c[[n + 1]]/2;
Do[
   a[0, k + 1] = c[[n - k]] + a[1, k]/2;
   a[1, k + 1] = a[0, k] + a[2, k]/2;
   Do[
      a[m, k + 1] = (a[m + 1, k] + a[m - 1, k])/2,
      {m, 2, k - 1}];
   a[k, k + 1] = a[k - 1, k]/2;
   a[k + 1, k + 1] = a[k, k]/2,
   {k, 2, n - 1}];

ccof = Table[a[m, n], {m, 0, n}]
   {-3.25778, -1.79187, 1.25145, -0.928, 0.383135}

To (symbolically or numerically) evaluate a polynomial expressed in the Chebyshev basis, we can use Clenshaw's algorithm like so:

chebeval[cof_?VectorQ, x_] := Module[{m = Length[cof], j, u, v, w},
  w = v = 0;
  Do[
     u = 2 x v - w + cof[[j]];
     w = v; v = u,
     {j, m, 2, -1}];
  Expand[x v - w + First[cof]]]

chebeval[ccof, u]
   -4.1261 + 0.99213 u - 0.562172 u^2 - 3.712 u^3 + 3.06508 u^4

and we see that the coefficients are neatly reproduced.

If one wants to evaluate the polynomial and its derivative simultaneously, rather than separately generating the Chebyshev coefficients of the derivative, one can instead symbolically differentiate Clenshaw's algorithm, like so:

chebdeval[cof_List, x_] := 
 Module[{m = Length[cof], d, e, f, j, u, v, w},
  w = v = 0;
  f = e = 0;
  Do[
   d = 2 (x e + v) - f;
   f = e; e = d;
   u = 2 x v - w + cof[[j]];
   w = v; v = u;,
   {j, m, 2, -1}];
  Expand[{x v - w + First[cof], x e + v - f}]]

chebdeval[ccof, x]
   {-4.1261 + 0.99213 x - 0.562172 x^2 - 3.712 x^3 + 3.06508 x^4,
    0.99213 - 1.12434 x - 11.136 x^2 + 12.2603 x^3}

(This approach is also readily extended if one wants the second, third... derivatives.)

share|improve this answer
    
Thanks you very much J.M. I'm not that familiar with the Clenshaw algorithm, but I will check the references you put up, and work through your code. This is so very helpful. –  tau1777 Oct 26 '12 at 18:07
    
Hello again. So I should have asked this earlier, but why does numerical recipes approximate functions as (Sum c_i T_i(y) - 1/2*c_0 where i=0, i=m-1)? When I approximated the function above I did not need to subtract off half the value of the c_0 constant. Sorry if this is a really basic question. Is it simply because I am doing integrals and they are actually doing the sum over a grid of values? Is that the significance of the m-1 in the sum? –  tau1777 Nov 13 '12 at 18:36
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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