Funny behavior when computing dot product of coefficients with high-order polynomials

Question

I have a similar problem to Funny behaviour when plotting a polynomial of high degree and large coefficients. However, the thing being evaluated is not just a polynomial, but a dot product of some coefficients with a list of polynomials.

(The back story is that I want to distort a sine random variable so that it has an approximate Laplace distribution. I am doing this by approximating the nonlinear mapping with a Chebyshev approximation. This has the cool property that if I use a maxN-th order Chebyshev series then there are only maxN harmonics of the sine wave. If the exact mapping is used, there are an infinite number of harmonics.)

The dot product (coefficients times Chebyshev polynomials) begins to become inaccurate (blow up, oscillate like crazy, whatever) near $\pm 1$ for maxN greater than the mid-$40$s. I've tried the Rationalize function on the coefficients--the polynomial coefficients are all integers--but to no avail.

Here is what I've tried. I have included the exact mapping for comparison to the Chebyshev approximation. (Mathematica told me that this involves the SinIntegral so I have hard-coded that into the definition.) The mapping has odd symmetry about the origin so even-order coefficients are zero.

FLaplaceInverse[x_, β_] := 
  Piecewise[{{β Log[2 x], x < 1/2}, {-β Log[2 - 2 x], 
     x >= 1/2}}];
FSine[x_, b_] := (π + 2 ArcSin[x/b])/(2 π);
CompositeSineLaplace[x_] := FLaplaceInverse[FSine[x, 1], 1];
maxN = 49;
chebCoeffList = 
  Table[If[EvenQ[n], 0, 
    N[(4 SinIntegral[(n π)/2])/(n π), 100]], {n, 0, maxN}];
rChebCoeffList = Rationalize[chebCoeffList, 0];
chebPolyList = 
  Table[If[EvenQ[n], 0, ChebyshevT[n, x]], {n, 0, maxN}];
plotBoth = 
  Plot[{CompositeSineLaplace[x], 
    rChebCoeffList.chebPolyList} , {x, -1, 1}, PlotRange -> Automatic,
    ImageSize -> Medium];
plotDiff = 
  Plot[Evaluate[
    CompositeSineLaplace[x] - chebCoeffList.chebPolyList] , {x, -1, 
    1}, PlotRange -> Automatic, ImageSize -> Medium];
chebCoeffsWithOrder = 
  Transpose[{Table[n - 1, {n, Length[chebCoeffList]}], chebCoeffList}];
plotCoeffs = 
  ListPlot[chebCoeffsWithOrder, PlotRange -> All, ImageSize -> Medium];
{plotBoth, plotDiff, plotCoeffs}

The Plot functions give this for maxN = 49.

Answer 1

To resolve this:

Compare

Plot[chebCoeffList.chebPolyList, {x, -1, 1}]

with

Plot[chebCoeffList.Table[If[EvenQ[n], 0, ChebyshevT[n, x]], {n, 0, maxN}],
     {x, -1, 1}]

The blowup seen in the first plot is, as noted, due to the fact that chebPolyList is already converted automatically to the monomial basis, which is terribly ill-conditioned for numerical evaluations of this sort. In contrast, the evaluation in the second one did not involve this conversion; that is, ChebyshevT[] was directly evaluated at the inexact numbers fed to it by Plot[]. Setting WorkingPrecision to a higher value in the first plot can stave off this ill-conditioning, but at the expense of increased effort.

---

As also noted in the comments, a more efficient way of proceeding would be to use the Clenshaw recurrence to evaluate your Chebyshev series. This has to be modified a bit to exploit the oddness of your function. Thus,

cheboddval[c_?VectorQ, x_] := 
           Module[{n = Length[c], u, v, w, y = 4 x^2 - 2}, 
                  u = c[[n - 1]] + y (v = c[[n]]);
                  Do[w = v; v = u; u = c[[k]] + y v - w, {k, n - 2, 2, -1}];
                  x (c[[1]] + (y - 1) u - v)]

Let's generate your coefficients without the zeroes taking up space:

cc = Table[N[(4 SinIntegral[n π/2])/(n π), 50], {n, 1, maxN, 2}];

Now, the plot:

Plot[{CompositeSineLaplace[x], cheboddval[cc, x]}, {x, -1, 1}]