I want to convert a polynomial in "standard form" to Chebyshev form. Here's one way to do it:
(* My polynomial *)
pa = Sum[a[i]*t^i, {i, 0, 5}]
(* "standard form" coefficients of Chebyshev polynomial of same degree *)
pb = CoefficientList[Sum[b[i]*ChebyshevT[i, t], {i, 0, 5}], t]
(* helper variable *)
bs = Table[b[i], {i, 0, 5}]
Solve[Table[a[i - 1] == pb[[i]], {i, 1, 6}], bs]
This works, but is there a better (more efficient, more Mathematica-ish) way to do it?
Answer
Let's say you have a polynomial:
f[x_] := 1 - x + 2 x^2 - 5 x^4
n = 5; (* Polynomial order + 1 *)
Now you have to solve a system of linear equations:
a = Transpose @ PadRight @ Table[CoefficientList[ChebyshevT[k - 1, x], x], {k, n}];
c = LinearSolve[a, CoefficientList[f[x], x]]
(* {1/8, -1, -3/2, 0, -5/8} *)
Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand
(* 1 - x + 2 x^2 - 5 x^4 *)
If getting the numerical values of coefficients is sufficient, you can use the discrete Fourier transform, which is an extremely fast method.
Tn(x)=cos(narccosx)
c = 2 MapAt[#/2 &, #, 1]/Sqrt[n] & @ FourierDCT @ f @ Cos[π Range[.5, n]/n]
(* {0.125, -1., -1.5, -9.93014*10^-17, -0.625} *)
Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand // Chop
(* 1. - 1. x + 2. x^2 - 5. x^4 *)
P.S. The last method can be used for the Chebyshev expansion of an arbitrary function (there can be a problem with Gibbs oscillations, but that is another story).
Comments
Post a Comment