11
$\begingroup$

I have been looking at how to interpolate a function using Chebyshev polynomials. There are several good posts such as this one by Michael E2 and this one by J. M. There is also an example in the help for FourierDCT. By studying these posts I have assembled a module that does interpolation of a function.

ClearAll[ChebInterpolation];
ChebInterpolation::usage = 
  "ChebInterpolation[f,{xmin,xmax},prec,n] developes an interpolation \
of the function f in the interval [xmin, xmax] based on n Chebyshev \
functions calculated with a precision prec. f must be a pure \
function. ";
ChebInterpolation[f_, {xmin_, xmax_}, prec_, n_] := 
 Module[{nodes, coeff, c},
  nodes = 
   Rescale[N[Cos[Pi Range[0, n]/n], prec], {-1, 1}, {xmin, xmax}];
  coeff = Sqrt[2/n] FourierDCT[f /@ nodes, 1];
  coeff[[{1, -1}]] /= 2;
  c = Compile[{x}, 
    coeff.Table[Cos[(m - 1) ArcCos[x]], {m, Length@coeff}]];
  With[{w = c}, w[Rescale[#, {xmin, xmax}, {-1, 1}]] &]
  ]

With this module I have to supply a precision (prec) and the number of Chebyshev polynomials (n) to use. My question is how to automate finding these two parameters.

I use the module like this. First define a function and the interval within which it will be approximated.

f = BesselJ[1, #^(3/2)] Sin[#] &;
{xmin, xmax} = {25, 35};
Plot[f[x], {x, xmin, xmax}]

Mathematica graphics

Now we run the module to calculate the interpolation function and plot the result.

fint = ChebInterpolation[f, {xmin, xmax}, 25, 85];
Plot[fint[x], {x, xmin, xmax}]

Mathematica graphics

The plot looks very similar and to check the accuracy we can look at the error

Plot[fint[x] - f[x], {x, xmin, xmax}, WorkingPrecision -> 100, 
 PlotRange -> All, PlotPoints -> 100]

Mathematica graphics

As can be seen the error is small. The parameters in the interpolation were a precision of 25 and 85 Chebyshev polynomials. I used these values by looking at various alternatives. The following calculation determines the log of the RMS error between the exact function and the approximate function for a range of precisions and number of polynomials.

 errs = Table[fint = ChebInterpolation[f, {xmin, xmax}, prec, n];
   e = Table[fint[x] - f[x], {x, xmin, xmax, (xmax - xmin)/100}];
   {prec, n, Log[10, RootMeanSquare[e]]},
   {n, 65, 100}, {prec, 16, 32}];
ListPlot3D[Flatten[errs, 1], PlotRange -> All, 
 AxesLabel -> {"Precision", "No of Chebs.", 
   "\!\(\*SubscriptBox[\(Log\), \(10\)]\) RMS Error"}]

Mathematica graphics

A floor is reached when the precision is about 25 and the number of polynomials is 85.

Question 1: Is my module a good approach or can it be improved?

Question 2: Is looking at the errors in the above plot a robust way of finding the required precision and number of polynomials? If yes, what would be a good way of automating this?

$\endgroup$
3
  • 1
    $\begingroup$ If you are interested in approximation using Chebyshev polynomials, you may find the documentation (if not the software) at chebfun.org very interesting. $\endgroup$
    – mikado
    Jul 16, 2018 at 21:34
  • $\begingroup$ @mikado Thanks for the reference which I am aware of. However it does not give algorithms. At least I have not seen them on that site. $\endgroup$
    – Hugh
    Jul 16, 2018 at 21:44
  • $\begingroup$ It is all based on published academic work, so I think that you will find the essentials described in detail. E.g. people.maths.ox.ac.uk/trefethen/publication/PDF/2004_107.pdf $\endgroup$
    – mikado
    Jul 16, 2018 at 22:13

1 Answer 1

4
$\begingroup$

Since we're compiling at the end, doing everything in machine precision seemed okay. We need some cushion for the stopping criterion, since ordinary roundoff error might make achieving a relative error of $MachineEpsilon impossible. I chose a factor of 10, which might not be high enough in every case.

If we think of the order-$n$ Chebyshev points as the $x$ coordinates of $n+1$ points equispaced around the upper half of the unit circle beginning and ending at the endpoints, then to get the order-$2n$ points, we need only bisect each arc. We can add the function values at the new points to the function values at the old points with Riffle[]. Hence we can recursively bisect the half circle until the desired precision is reached.

Since the Chebyshev polynomials satisfy $-1 \le T_j(x) \le 1$ for $-1 \le x \le 1$, the error of truncating the Chebyshev series at order $n$ is bounded by $$\sum_{j=n+1}^\infty |c_j|\,$$ where $c_j$ are the Chebyshev coefficients. For an analytic function, the coefficients decrease geometrically, and the last coefficient $|c_n|$ can be expected to be bound the error, except that any particular coefficient can potentially have any value (just add $a\,T_j(x)$ to your function and you increase the $j$-th coefficient by an arbitrary amount $a$). But this is rare. In practice there's one common such snafu, namely an even or an odd function. If you take the maximium of the last few coefficients (at least two), you can avoid all but the most exceptional problems.

ClearAll[ChebInterpolation];
ChebInterpolation::usage = 
  "ChebInterpolation[f,{xmin,xmax},prec,n] developes an interpolation \
of the function f in the interval [xmin, xmax] based on n Chebyshev \
functions calculated with a precision prec. f must be a pure \
function. ";
ChebInterpolation[f_, {xmin_, xmax_}, prec_] := Module[{nodes, coeff, c
   , n = 8, y, ynew, error = Infinity, tol, len, sum, abscc, maxcc},
  
  (* Need cushion for ordinary roundoff error *)
  tol = Max[10.^-prec, 10*$MachineEpsilon];
  nodes = Rescale[N[Cos[Pi Range[0, n]/n]], {-1, 1}, {xmin, xmax}];
  y = f /@ nodes;
  While[
   error > tol && n < 2^24,
   n = 2*n;
   nodes = Rescale[N[Cos[Pi Range[1, n, 2]/n]], {-1, 1}, {xmin, xmax}];
   ynew = f /@ nodes;
   y = Riffle[y, ynew]; (* riffle new func vals *)
   coeff = Sqrt[2/n] FourierDCT[y, 1];
   coeff[[{1, -1}]] /= 2;
   (* error estimate, save abs/max cc for after *)
   abscc = Abs@coeff;
   maxcc = Max@abscc;
   error = Max[abscc[[-3 ;;]]]/maxcc (* use last 3 coeff to est. 
   error *)
   ];
  (* trim coeffs *)
  sum = 0;
  len = LengthWhile[Reverse@abscc, (sum += #) <= tol*maxcc &];
  coeff = If[TrueQ[len > 1], Drop[coeff, -len], coeff];
  (*foo=coeff;*)
  c = Compile[{x},
    Block[{cc = coeff},
     cc . Table[Cos[(m - 1) ArcCos[x]], {m, Length@cc}]
     ], CompilationOptions -> {"InlineExternalDefinitions" -> True}];
  With[{w = c}, w[Rescale[#, {xmin, xmax}, {-1, 1}]] &]]

Example:

f = BesselJ[1, #^(3/2)] Sin[#] &;
{xmin, xmax} = {25, 35};
fint = ChebInterpolation[f, {xmin, xmax}, 25];
Plot[fint[x] - f[x], {x, xmin, xmax}, WorkingPrecision -> 100, 
 PlotRange -> All, PlotPoints -> 100]

Another example (precision goal = 8):

ff = ChebInterpolation[Sin, {-Pi, Pi}, 8];
Plot[ff[x] - Sin[x], {x, -Pi, Pi}]

Looks pretty close to the minimax approximation!

$\endgroup$
1
  • $\begingroup$ A very useful development. Thank you. $\endgroup$
    – Hugh
    Apr 2, 2022 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.