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 have a function:

f[x_] := x + 31 x^3 + 5 x^25

Which I want to find an expansion for:

ex[x_, a_, b_, c_, d_] := a + b x + c x^3 + d x^25

I do this by solving to find the coefficients a,b,c,d by setting the expansion equal to the function at the appropriate number of points:

sol = Solve[{ex[1, a, b, c, d] == f[1], ex[3, a, b, c, d] == f[3], 
ex[2, a, b, c, d] == f[2], ex[4, a, b, c, d] == f[4]}]

This works fine, but its ugly. I've only used this to illustrate what I want to do... but actually I want to find an expansion involving 'N' expansion coefficients, and perform solve at N different points to find what they are. However, I don't understand how to do that without having functions with N+1 arguments and so on.

If I haven't expressed myself clearly let me know and I'll add more.

EDIT: Perhaps I'll add that really the expansion I'm looking for is something like (though how I put the coefficients into the function argument I don't know):

ex[x_,a[i]_,b[i]_,N_,M_]:= Sum[a[n] ChebyshevT[n,x],{n,0,N}] + Sum[b[n] * y[n,x],{n,1,M}]

Where the functions y[n,x] are some other functions which is not orthogonal to the Chebyshev polynomials.

share|improve this question
    
Is your function always a polynomial ? –  b.gatessucks Feb 12 '13 at 15:38
    
The function isn't but the expansion basis functions will be Chebyshev Polynomials, I just used that to illustrate my point. –  user5866 Feb 12 '13 at 15:40
1  
If your basis isn't orthogonal, your expansion won't (probably) be unique ... –  belisarius Feb 12 '13 at 16:11
    
As in more than one set of a[n] and b[n] will give f[x] to arbitrary accuracy? (that wouldn't necessarily be a problem but I hadn't realised that) –  user5866 Feb 12 '13 at 16:16
    
Something like Solve[Table[ex[j, a, b, c, d] == f[j],{j,nn+mm}],Join[Table[a[n],{n,nn}],Table[b[n],{n,mm}]]]? –  Daniel Lichtblau Feb 12 '13 at 16:21
show 3 more comments

2 Answers

(Edited; original at the end)

Gram-Schmidt orthogonalization provides an answer. Let's use a running example to illustrate. It begins even before the Chebyshev polynomials, with their domain--the interval $[-1,1]$--and the kernel for which they are orthogonal:

limits = {-1, 1};
k[x_] := 2/Sqrt[1 - x^2]/Pi;

The Chebyshev polynomials are obtained by means of the Gram-Schmidt orthogonalization process from the monomial functions $(1, x, x^2, \ldots, x^n, \ldots)$. For instance, here are the first four of them as obtained from the powers $0$ through $3$:

 Orthogonalize[{1, x, x^2, x^3}, Integrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}] &]  
  // Expand

$\left\{\frac{1}{\sqrt{2}},x,2 x^2-1,4 x^3-3 x\right\}$

(If you compare this to ChebyshevT[#, x] & /@ Range[0, 3] // TraditionalForm you will see that the first is rescaled but all the rest are the same; this is because Mathematica for some reason does not use an orthonormal basis for ChebyshevT.)

Thus, in order to obtain an expansion in terms of a finite set of Chebyshev polynomials and other functions, first orthogonalize them, normalize those if necessary, and obtain the coefficients in terms of the inner products.

There's a problem: unless the collection of functions is very nice, Mathematica will be unable to obtain closed-form expressions for the orthogonal basis. Instead, we can numerically integrate.

For example, let's start with a few power functions--they will generate Chebyshev polynomials automatically--together with some stranger transcendental ones (plotted below):

vectors = Function /@ Table[#^k, {k, 0, 6}];
y = {Sqrt[1 - Abs[#]] &, Exp, # Abs[Log[#]] &};
Plot[Evaluate@Through[y[x]], {x, limits[[1]], limits[[2]]}]

Interesting functions

In this fashion we have recreated the situation presented in the problem: we seek to expand arbitrary functions in terms of Chebyshev polynomials--which are mutually orthogonal--and some other functions y which are not necessarily orthogonal to the Chebyshev polynomials.

Obtaining an orthogonal basis takes a few seconds and some fiddling to compute the integrals with sufficient accuracy:

t = Orthogonalize[Through[(vectors~Join~y)[x]], 
 NIntegrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}, 
   MaxRecursion -> 12] &] // Simplify // Chop // Rationalize

Let's check that this basis is really orthonormal (if not, we would have to normalize its elements):

Outer[NIntegrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}, MaxRecursion -> 12] &, t, t]

I applied Chop[%, 10^-6] to remove some small near-zero values and obtained

$$\left( \begin{array}{cccccccccc} 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1. \end{array} \right)$$

The off-diagonal zeros confirm orthogonality. The lengths of the vectors are on the diagonal: already normalized, as intended. Notice the zero on the diagonal! To within the limits of numerical error, it appears one of these functions is a linear combination of the others. Nevertheless, we can still proceed without having to take any special steps. Let's pause first to look at these basis functions:

GraphicsGrid[{Plot[Evaluate@#, {x, limits[[1]], limits[[2]]}] & /@ Partition[t, 5]}]

Basis functions

Definitely some of the later ones (on the right) aren't pure Chebyshev polynomials!

Let's expand a function in terms of this basis; some function that clearly is not within the linear span of the original set of functions:

f[x_] := x + 31 x^3 + 5 x^25
a = NIntegrate[f[x] # k[x], {x, limits[[1]], limits[[2]]}] & /@ t // Chop;
a . t // Simplify

$2.29258 \sqrt{1-|x|}-6.43222 x |\log (x)|-3.82443 x^6+21.3062 x^5+10.7913 x^4+6.02688 x^3-14.6504 x^2+18.7429 x-2.6231$

(Numerical error is not as critical in this calculation, so we don't have to be so fiddly.) The vector a is the desired set of expansion coefficients. Its inner product with t expresses the expanded value: it is the orthogonal projection of $f$ onto the space of functions spanned by the functions we defined earlier in vectors and y. ("Orthogonal" means relative to the kernel $k$.)

How does the original function compare to its expansion?

Plot[{f[x], a.t}, {x, limits[[1]], limits[[2]]}]

Plots of the function and its expansion

This looks pretty good. However, because f gets large, we can't see small differences. Let's look at the residuals of the expansion instead. There are two meaningful ways to do this: (a) by subtracting f from its expansion and (b) by dividing those differences by the values of the kernel: these "standardized" differences ought to be fairly uniform across the domain. Here are both plots:

GraphicsGrid[{{Plot[a.Through[nCn[x]] - f[x], {x, limits[[1]], limits[[2]]}, PlotStyle -> Thick], 
   Plot[ (a.Through[nCn[x]] - f[x]) / k[x], {x, limits[[1]], limits[[2]]}]}}]

Residual plots

These errors are pretty small compared to the original range of $f$ and, as expected, the standardized errors are reasonably uniform.


Original Reply

Chebyshev polynomials $\psi_n(x)$ are orthogonal on $[-1,1]$ with respect to $\frac{dx}{\sqrt{1-x^2}}$. The normalization used by Mathematica makes all but the zeroth have norm $\pi/2$; the zeroth has norm $\pi$; e.g.,

Integrate[ChebyshevT[#, x]^2 / Sqrt[1 - x^2] , {x, -1, 1}] & /@ Range[0, 5]

$\left\{\pi ,\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2}\right\}$

Brute Force Method

Therefore the coefficient of the $n$th Chebyshev polynomial is found by integrating $f$ against ChebyshevT[n, x] / Sqrt[1 - x^2] and dividing by $\pi/2$ (by $\pi$ when $n=0$). (For an odd function $f$ only the odd coefficients will be nonzero: knowing this can halve the execution time.) Only the Chebyshev polynomials up to the degree of $f$ have to be evaluated:

f[x_] := x + 31 x^3 + 5 x^25;
degree = Length[CoefficientList[p[x], x]] - 1;
a = Integrate[(2/Pi) p[x] ChebyshevT[#, x] / Sqrt[1 - x^2] , {x, -1, 1}] & /@ Range[0, degree];
a = a Prepend[ConstantArray[1, degree], 1/2]

$\left\{0,\frac{108212247}{4194304},0,\frac{19038803}{2097152},0,\frac{2042975}{2097152}, \ldots, \frac{125}{16777216},0,\frac{5}{16777216}\right\}$

As a check, compute the expansion to verify it equals the original polynomial:

a . (ChebyshevT[#, x] & /@ Range[0, degree]) // Expand

$5 x^{25}+31 x^3+x$

Fast Method

Using the preceding, it's straightforward to check and prove that the coefficient of the $n$th Chebyshev polynomial in the expansion of $x^k$ is related to a Binomial coefficient:

Clear[c];
c[k_, n_] /; n <= k && EvenQ[n - k] := c[k,n] 
  = 2^(1 - k) Binomial[k, (k - n)/2] If[n == 0, 1/2, 1];
c[k_, n_] := 0

Whence, to expand any polynomial f, just replace the powers of its argument by their expansions:

f[x] /. {Times[a_, Power[x, e_]] :> 
 a Sum[c[e, k] Subscript[\[Psi], k][x], {k, Mod[e, 2], e, 2}], 
 x -> Subscript[\[Psi], 1][x]} // Expand 

$\frac{108212247 \psi _1(x)}{4194304}+\frac{19038803 \psi _3(x)}{2097152}+\frac{2042975 \psi _5(x)}{2097152}+\ldots+\frac{125 \psi _{23}(x)}{16777216}+\frac{5 \psi _{25}(x)}{16777216}$


A similar approach will work with other systems of orthogonal polynomials.

share|improve this answer
    
I posted a similar answer earlier and deleted it because the OP commented Sorry, that is completely true, but it won't work for the actual expansion I want to do which is a mixture of Chebyshev and other singular functions. I'm sorry, i'll add more information, I was just trying to keep it simple for readability –  belisarius Feb 12 '13 at 22:06
    
@belisarius Thank you for the clarification. Because the techniques in this answer might still have a bearing on the reformulated question, I will leave it up in the meantime. –  whuber Feb 12 '13 at 22:57
    
But .. why don't you post a question about that and then answer it yourself? The OP's problem is different! –  belisarius Feb 12 '13 at 23:29
    
@belisarius That's a good suggestion, but I still think that what I originally wrote contributes to the question. Instead, I have--I believe--now fully answered the question in an edit. I have retained the original because for some readers it might help clarify or streamline the application of the first part. –  whuber Feb 13 '13 at 19:25
add comment

You could do something like:

b[i_, x_] := ChebyshevT[i - 1, x];(*Tch basis as an example*)
app[f_, pSet_]:= 
     Solve[Table[Sum[a@i b[i, x], {i,Length@#}]==f@x, {x, #}], Array[a,Length@#]] &@  pSet

Where pSet is the list of points (x-values) you use for equaling the functions.

Use it like this:

app[# + 31 #^3 + 5 #^25 &, Range@5]

Edit

Note that for example:

app[# + 31 #^3 + 5 #^25 &, Range@3]

Expands to:

Solve[{ a[1] +   a[2] +    a[3] == 37, 
        a[1] + 2 a[2] + 7  a[3] == 167772410, 
        a[1] + 3 a[2] + 17 a[3] == 4236443048055}, 
     {a[1], a[2], a[3]}]

So you could also try some linear method directly

share|improve this answer
    
Thanks a lot for your response What is the role of pSet? Is it just a list of points where we set the expansion equal to the function? If so can I just use pSet=Table[n,{n,1,N}] to give a list of N points? also should the 1st and 2nd arguments of your Solve function be enclosed in { } or not? –  user5866 Feb 12 '13 at 16:56
    
never mind I think its a pointer to a set, and you've used Range in the place of my suggested table. so much still to learn –  user5866 Feb 12 '13 at 17:09
    
@user5866 pSet is a ... points set :) –  belisarius Feb 12 '13 at 17:35
    
@user5866 There are no pointers natively in Mathematica –  belisarius Feb 12 '13 at 18:25
    
This works as an answer, but I lack the intelligence to expand it to the case I actually need, which isn't just a 1,2,3 range but rather N set of points. What is it that you are passing with the Range@3 argument? I had it in my head that this was the address of a range of 3 numbers that you were pointing to with pSet.... :S –  user5866 Feb 12 '13 at 19:49
show 2 more comments

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.