# Chebyshev Approximation

Is there functionality in Mathematica to expand a function into a series with Chebyshev polynomials?

The Series function only approximates with Taylor series.

• Do you mean by using the formula's in this paper, specifically these? If so, ChebyshevT would be the place to start – Jason B. Jan 5 '16 at 16:52
• (1) Post a specific example. (2) Have a look at ChebyshevT documentation, in particular "Generalizations and Extensions", and "Applications" – Daniel Lichtblau Jan 5 '16 at 16:54
• This question applied to WhitakkerW may be found here: mathematica.stackexchange.com/a/98350 – Michael E2 Jan 5 '16 at 18:14
• And I got the technique from here: mathematica.stackexchange.com/a/89692 – Michael E2 Jan 5 '16 at 18:17
• Here is another related question with very good answers. – Oleksandr R. Jan 5 '16 at 22:13

You can just take Bob Hanlon's answer from 2006 directly, and modify the plot just a bit to update it.

ChebyshevApprox[n_Integer?Positive, f_Function, x_] :=
Module[{c, xk}, xk = Pi (Range[n] - 1/2)/n;
c[j_] = 2*Total[Cos[j*xk]*(f /@ Cos[xk])]/n;
Total[Table[c[k]*ChebyshevT[k, x], {k, 0, n - 1}]] - c[0]/2];

f = 3*#^2*Exp[-2*#]*Sin[2 #*Pi] &;

ChebyshevApprox[3, f, x] // Simplify

((-(3/4))*((-E^(2*Sqrt[3]))*(Sqrt[3] - 2*x) - 2*x - Sqrt[3])*x*
Sin[Sqrt[3]*Pi])/E^Sqrt[3]

GraphicsGrid[
Partition[
Table[Plot[{f[x], ChebyshevApprox[n, f, x]}, {x, -1, 1},
Frame -> True, Axes -> False, PlotStyle -> {Blue, Red},
PlotRange -> {-2, 10},
Epilog -> Text["n = " <> ToString[n], {0.25, 5}]], {n, 9}], 3],
ImageSize -> 500]


• I am having trouble using this with non-pure functions. For example, I can not expand Sin with ChebyshevApprox[3, Sin, x]. It also does not plot anything. – Sogartar Jan 6 '16 at 14:11
• You can make a pure function out of Sin, via ChebyshevApprox[3, Sin[#] &, x]. And the result is a polynomial in x, so you would plot it in the normal way: Plot[{Sin[x], ChebyshevApprox[3, Sin[#] &, x]}, {x, -1, 1}] - but for Sin only one term survives so you have a linear function. Going up to n=4 you get a cubic function and it is indistinguishable on that range. – Jason B. Jan 6 '16 at 14:27

Here's a way to leverage the Clenshaw-Curtis rule of NIntegrate and Anton Antonov's answer, Determining which rule NIntegrate selects automatically, to construct a piecewise Chebyshev series for a function. It also turns out that InterpolatingFunction implements a Chebyshev series approximation as one of its interpolating units (undocumented). With IntegrationMonitor, you can save the sampling on the subintervals and use FourierDCT[] to convert the function values to Chebyshev coefficients. The error estimate for the integral ($E \approx |I_{n/2}-I_{n}|$) is not exactly the same as an approximation norm, but most numerical procedures have pitfalls.

fn = 3*#^2*Exp[-2*#]*Sin[2 #*Pi] &[x];
{ifn, {{series}}} = Reap[
chebApprox[fn, {x, 0, 2}],             (* see below for code *)
"ChebyshevSeries"]; // AbsoluteTiming
(*  {0.003718, Null}  *)

Plot[{fn, ifn[x]}, {x, 0, 2},
PlotStyle -> {AbsoluteThickness[5], Automatic}],
LogPlot[(fn - ifn[x])/fn // Abs, {x, 0, 2}]


The reaped series are in a list of the form {{{x1, x2}, {coefficients}},...}.

Short[series, 12]


In this example there are four series over the intervals

series[[All, 1]]
(*  {{0, 0.5}, {0.5, 1.}, {1., 1.5}, {1.5, 2}}  *)


Code dump

ClearAll[chebInterpolation];
(* Constructs a piecewise InterpolatingFunction,
* whose interpolating units are Chebyshev series *)
(* data0 = {{x0,x1},c1},..} *)
chebInterpolation[data0 : {{{_, _}, _List} ..}] :=
Module[{data = Sort@data0, domain1, coeffs1, domain, grid, ngrid,
coeffs, order, y0, yp0},
domain1 = data[[1, 1]];
coeffs1 = data[[1, 2]];
domain = List @@ Interval @@ data[[All, 1]];
y0 = chebFunc[coeffs1, domain1, First@domain1];
yp0 = chebFunc[dCheb[coeffs1, domain1], domain1, First@domain1];
grid = Union @@ data[[All, 1]];
ngrid = Length@grid;
coeffs = data[[All, 2]];
order = Length[coeffs[[1]]] - 1;
InterpolatingFunction[
domain,
{5, 1, order, {ngrid}, {4; order}, 0, 0, 0, 0, Automatic, {}, {}, False},
{grid},
{{y0, yp0}} ~Join~ coeffs,
{{{{1}}~Join~Partition[Range@ngrid, 2, 1] ~Join~ {{ngrid - 1, ngrid}},
{Automatic } ~Join~ ConstantArray[ChebyshevT, ngrid]}}] /;
Length[domain] == 1 && ArrayQ@coeffs
];

Clear[chebApprox];
(* Uses NIntegrate's Clenshaw-Curtis Rule
* to construct Chebyshev series approximations to a function
* over the subintervals created by NIntegrate *)
Options[chebApprox] = {"Points" -> 17} ~Join~ Options[NIntegrate];
chebApprox[f_, {x_, a_, b_}, opts : OptionsPattern[]] :=
Module[{t, samp, sampling},
With[{pg = OptionValue[PrecisionGoal] /. Automatic -> 14,
ag = OptionValue[AccuracyGoal] /. Automatic -> 15},
t = Reap[
NIntegrate[f, {x, a, b},
Method -> {"ClenshawCurtisRule",
"Points" -> OptionValue["Points"]},
PrecisionGoal -> pg, AccuracyGoal -> ag,
IntegrationMonitor :> (Sow[
Map[{First[#1@"Boundaries"], #1@"GetValues"} &, #1],
samp] &),
Evaluate@FilterRules[{opts}, Options[NIntegrate]]],
samp];
sampling = With[{steps = t[[2, 1]]},
Flatten[
Table[If[MemberQ[steps[[n + 1]], {{s[[1, 1]], _}, __}],
Nothing, s], {n, Length@steps - 1}, {s, steps[[n]]}], 1] ~Join~
DeleteCases[Last@steps, {{-Infinity, Infinity}, __}]
];
sampling = Sort@MapAt[chebSeries@*Reverse, sampling, {All, 2}];
Sow[sampling, "ChebyshevSeries"];
chebInterpolation[sampling]
]];

chebSeries[y_] := Module[{cc},
cc = Sqrt[2/(Length@y - 1)] FourierDCT[y, 1]; (* get coeffs from values *)
cc[[{1, -1}]] /= 2;  (*
adjust first & last coeffs *)
cc
]

(* Differentiate a Chebyshev series *)
(* Recurrence: $2 r c_r = c'_{r-1} - c'_{r+1}$ *)
Clear[dCheb];
dCheb::usage =
"dCheb[c, {a,b}] differentiates the Chebyshev series c scaled over \
the interval {a,b}";
dCheb[c_] := dCheb[c, {-1, 1}];
dCheb[c_, {a_, b_}] := Module[{c1 = 0, c2 = 0, c3},
2/(b - a) MapAt[#/2 &,
Reverse@ Table[
c3 = c2;
c2 = c1;
c1 = 2 (n + 1)*c[[n + 2]] + c3,
{n, Length[c] - 2, 0, -1}],
1]
];


Notes and references:

One can get a single Chebyshev series by setting MaxRecursion -> 0.

{{domain}, values} =
Reap[NIntegrate[fn, {x, 0, 2}, PrecisionGoal -> 14,
AccuracyGoal -> 15,
Method -> {"ClenshawCurtisRule",
"Points" -> 1 + 2^5},     (* adjust "Points" to achieve desired accuracy *)
MaxRecursion -> 0,
WorkingPrecision -> 40,
IntegrationMonitor :> (Sow[
Map[{#1@"Boundaries", #1@"GetValues"} &, #1]] &)]
][[2, 1, 1, 1]];
cs = Module[{n = 1, max, sum = 0, ser, len},
ser = chebSeries[Reverse@values];
max = Max[ser];
len = LengthWhile[Reverse[ser], (sum += Abs@#) < 10^-22*max &];
Drop[N@ser, -len]
];

approx[x_?NumericQ] :=
cs.Table[ChebyshevT[n - 1, Rescale[x, domain, {-1, 1}]], {n, Length@cs}]

Plot[{fn, approx[x]}, {x, 0, 2},
PlotStyle -> {AbsoluteThickness[5], Automatic}],
LogPlot[(fn - approx[x])/fn // Abs, {x, 0, 2}]


• Now, you've previously told me how you discovered the hidden support for piecewise Chebyshev interpolants, but I'm still amazed that we're just about close to having chebfun being natively supported. Thanks for this! – J. M. is in limbo May 6 '16 at 0:33

One slick way to derive the analytic Chebyshev series of a function is to use the relationship between the Chebyshev polynomials and the cosine, and then use the built-in FourierCosSeries[]. As an example:

f[x_] := Exp[x];
n = 5; (* degree of approximation *)
approx[x_] = FourierCosSeries[f[Cos[t]], t, n] /. Cos[k_. t] :> ChebyshevT[k, x]


(Note that the result of that evaluation contains modified Bessel functions of the first kind, which arise as the coefficients.)

{Plot[{f[x], approx[x]}, {x, -1, 1}],
Plot[approx[x] - f[x], {x, -1, 1}, PlotStyle -> ColorData[97, 3]]} // GraphicsRow


See how good the approximation is? Note the equiripple behavior of the error at the right.

Of course, not every function will admit a closed form Chebyshev series representation, since the Fourier integrals involved won't necessarily have a closed form known to Mathematica. In that case, you can of course use NIntegrate[] instead. In fact, Mathematica does provide a package for numerically evaluating those integrals. Thus,

Needs["FourierSeries"]

f[x_] := 3 x^2 Exp[2 x] Sin[2 π x];

n = 12;
cof = Table[If[k == 0, 1/2, 1] NFourierCosCoefficient[f[Cos[t]], t, k, Method -> "LevinRule"],
{k, 0, n}];

(* Clenshaw recurrence for a Chebyshev series *)
chebval[c_?VectorQ, x_] := Module[{n = Length[c], u, v, w},
u = c[[n - 1]] + 2 x (v = c[[n]]);
Do[w = v; v = u; u = c[[k]] + 2 x v - w, {k, n - 2, 2, -1}];
c[[1]] + x u - v]

approx[x_] = chebval[cof, x];
{Plot[{f[x], approx[x]}, {x, -1, 1}],
Plot[approx[x] - f[x], {x, -1, 1}, PlotStyle -> ColorData[97, 3]]} // GraphicsRow


though the approximation in this case is not too good.

• The coefficients can also be computed using FourierDCT` as shown here – Stelios Feb 22 '16 at 13:21
• @Stelios, yes, using the discrete cosine transform corresponds to using the trapezoidal rule to discretize the underlying Fourier integrals. The answers linked by Michael in the comments use that strategy as well. – J. M. is in limbo Feb 22 '16 at 13:33