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 would like to expand a function $f(r)$ in the domain $[0,R]$, around the points $r =0$, and $r = R$ in the following manner

$f(r = 0) = \Sigma_{i=0,i = even}^{imax} f_i (r/R)^i$

and

$f(r = R) = \Sigma_{k=0}^{kmax} f_k (1 - r/R)^k$

I would like to get the coefficients $f_i$ and $f_k$ by first decomposing $f(r)$ into Chebyshev's and then turn that into the power series. Here is what I have so far for a random function :

g[x_] = 0.17768*x^5 + 0.115594*x^4 - 0.049490*x^3 - 0.659085*x^2 - 0.2254209;

The function looks like

enter image description here

Now to expand $g[x]$ in terms of Chebyshev I do this

In[8]:= GS = 6; (*gridsize*)

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

In[10]:= A = Array[a, GS + 1, 0]

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

In[11]:= For[i = 1, i <= GS + 1, i++, a[i] = 2/\[Pi]*   
 (Integrate[g[y]*ChebyshevT[i, y] 1/Sqrt[1 - y^2], {y, -1, 1}])]

In[12]:= a[0] = 1/\[Pi]* (Integrate[g[y]*ChebyshevT[0, y] 1/Sqrt[1 - y^2],    
 {y, -1, 1}])

Out[12]= -0.511616

In[13]:= A

Out[13]= {-0.511616, 0.0739325, -0.271746, 0.0431525, 0.0144492, 0.011105, -1.00221*10^-15}

In[14]:= t1 = Table[a[i]*ChebyshevT[i, y]  , {i, 0, GS}]

Out[14]= {-0.511616, 0.0739325 y, -0.271746 (-1 + 2 y^2), 0.0431525 
(-3 y + 4 y^3), 0.0144492 (1 - 8 y^2 + 8 y^4), 0.011105 (5 y - 20 y^3 + 16    
y^5), -1.00221*10^-15 (-1 + 18 y^2 - 48 y^4 + 32 y^6)}

In[15]:= g1[y_] = Sum[t1[[i]], {i, 1, GS + 1}]  

Out[15]= -0.511616 + 0.0739325 y - 0.271746 (-1 + 2 y^2) 
+ 0.0431525 (-3 y + 4 y^3) + 0.0144492 (1 - 8 y^2 + 8 y^4) 
+ 0.011105 (5 y - 20 y^3 + 16 y^5) -1.00221*10^-15 (-1 + 18 y^2 - 48 y^4 +      
32 y^6)

Now I plot the function and the Chebyshev approximation $g1[x]$.

enter image description here

and you can not really see the difference between the two plots, so this makes me think I am doing the Chebyshev part correctly.

After this I attempt to get the power series approximation using (Numerical Recipes sectinon 5.10, Clenshaw recurrence)

In[17]:= Clear[D1, d]

In[18]:= D1  = Array[d, GS]

Out[18]= {d[1], d[2], d[3], d[4], d[5], d[6]}

In[19]:= d[GS + 1] = 0

Out[19]= 0

In[20]:= d[GS] = 0

Out[20]= 0

In[23]:= For[i = GS - 1, i >= 1, i--, d[i] = 2*x*d[i + 1] - d[i + 2] + a[i]]

In[24]:= D1

Out[24]= {0.041885 - 2 (0.0144492 + 0.02221 x) x 

+ 2 x (-0.286195 - 0.02221 x + 2 x (0.0320475 + 2 (0.0144492 + 0.02221 x) x)), -0.286195 - 0.02221 x + 2 x (0.0320475 + 2 (0.0144492 + 0.02221 x) x), 0.0320475 + 2 (0.0144492 + 0.02221 x) x, 0.0144492 + 0.02221 x, 0.011105, 0}

In[27]:= d[0] = x*d[1] - d[2] + a[0]

Out[27]= -0.225421 + 0.02221 x - 2 x (0.0320475 + 2 (0.0144492 + 0.02221 x) x) + x (0.041885 - 2 (0.0144492 + 0.02221 x) x + 2 x (-0.286195 - 0.02221 x + 

2 x (0.0320475 + 2 (0.0144492 + 0.02221 x) x)))

In[28]:= f[x_] = Simplify[d[0]]

Out[28]= -0.225421 + 5.48173*10^-16 x - 0.659085 x^2 - 0.04949 x^3 + 0.115594 x^4 + 0.17768 x^5

Now I plot all there versions

enter image description here

and again I can see no difference between different versions.

I take this to mean that I have used the Clenshaw recurrence correctly.

Question 1: Did I really use it correctly or did I just get lucky?

Question 2: How would I change my code to deal with the function $f(r)$ around the origin which has to be normalized in its domain?

Question 3: Where do I even begin with $f(r)$ around $r =R$ since that does not have the normal power series form?

share|improve this question
    
What do you mean in Q2 that $f(r)$ has to be normalized in its domain? –  rcollyer Sep 27 '12 at 18:04
    
@rcollyer. The question I am trying to ask is: knowing how to take a function $g[x]$ in this case and expanding it in terms of a power series. How would I do the same thing, but instead of expanding in $x$, expand in $x/X$ where $X$ is the maximum value $x$ can have? So how would I have to change the code above to allow for expansions over the variable $(x/X)$? –  tau1777 Sep 27 '12 at 22:12
    
In this answer, I gave some Mathematica code for doing a Chebyshev expansion, starting from power series coefficients. On the other hand, it would seem to me that what you really want is a two-point Padé approximant, as opposed to a Chebyshev approximant. In that approximation, things can be set up such that the Taylor expansion of the approximant agrees with the Taylor expansion of the original function at two different points. –  J. M. Sep 27 '12 at 23:18

1 Answer 1

I am going to go out on a limb here as my mathematics is not as solid as I would like in this area. But, the primary difference between a Taylor series and expansion in terms of Chebyshev polynomials is the Chebyshev expansion is global while the Taylor series is not. Hence, the phrase "expanding around a point" is not applicable in the Chebyshev case. Instead, you are changing the basis your function is expressed in. Quite literally, the integral

$$f_i = \int^1_{-1} \frac{T_i(x) f(x)}{\sqrt{1 - x^2}} dx$$

is the projection of $f(x)$ onto the polynomial $T_i(x)$. Fourier series operates the same way.

As to your first question of whether you got lucky: no, my own experimentation with this several years ago reveals the same behavior. For instance, here's a comparison of the absolute errors ($\operatorname{abs}(f(x) - \sum_i f_i(x))$) for a 14 term Chebyshev vs. a 14 term Taylor expansion about 0 for $\sin x$:

enter image description here

The Taylor expansion does quite well for a broad range, but the Chebyshev is consistently at the $10^{-18}$ level which is on par with machine precision. In fact, to plot it correctly, I had to increase the WorkingPrecision to see anything at all!


To turn this into a power series, simply substitute in the Chebyshev polynomials into the expansion

$$f(x) \approx \sum^N_i f_i T_i(x)$$

and rearrange. For the first three terms of $\sin(x)$, the Chebyshev expansion is

0.999979 x - 0.166497 x^3 + 0.00799225 x^5

which the Taylor series closely resembles

x - 0.166667 x^3 + 0.00833333 x^5

Obviously, they diverge more extensively as the number of terms increases.

share|improve this answer
    
Great, so I did use the Clenshaw recurrence formula correctly. Now the confusion that remains for me is, how to get power series expansions if I want them. I have learned that Chebyshevs are better because they are global fcns, etc. But if I still wanted a power series expansion, for lets say $f(r=R)$ how can I get that from the Clenshaw recurrence? –  tau1777 Sep 27 '12 at 22:14
    
The plot on the left demonstrates the equi-ripple error property of a Chebyshev approximation. –  J. M. Sep 27 '12 at 23:15
    
@tau1777 I updated my answer. It is really straightforward. –  rcollyer Sep 28 '12 at 1:25
    
@J.M. right, but at that level the ripple doesn't matter. –  rcollyer Sep 28 '12 at 1:26

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.