What is the Chebyshev method for `InterpolationFunction`?

Question

With InterpolationOrder -> All, some methods of NDSolve returns interpolations with Method -> "Chebyshev". E.g.

NDSolve[{x'[t] == x[t], x[0] == 1}, x[t], {t, 0, 1},
  InterpolationOrder -> All, Method -> "ExplicitRungeKutta"]

What is the Chebyshev method? I couldn't find it explained in the documentation nor on this site. Does Method -> "Chebyshev" refer to some widely used method with a more specific name to google?

Answer 1

I guess my comment should be an answer:

Basically, my understanding is this, based more on textbook descriptions of DE solving (and extensive playing with Interpolation) than on Mathematica documentation: InterpolationOrder -> All causes NDSolve to construct a polynomial of degree equal to the order of the method for each step. For the integration methods "ExplicitRungeKutta", "ImplicitRungeKutta", "Extrapolation" and "DoubleStep", this polynomial is stored as a Chebyshev series in the InterpolatingFunction. (For other integration methods, it is stored as a local Taylor Series. When the order of the method is small, cubic Hermite interpolation is used. See The only usage for the option InterpolationOrder in NDSolve is to be set to All?) I tested the evaluation of the InterpolatingFunction, and it agrees with the Clenshaw algorithm (at least when using MachinePrecision). I presume that is what is used under the hood.

I used the interpolation method "Chebyshev" for chebInterpolation. Unfortunately, currently one has to construct the InterpolatingFunction[] directly. There is no option setting such as Method -> "Chebyshev" supported by Interpolation, which I think is a pity. I think users should have access to all the standard interpolation methods. The problem with using the undocumented construction via InterpolatingFunction is that if you make a mistake, you find out in one of two ways: The InterpolatingFunction does not format properly or the kernel crashes. Maybe sometimes you can get an error message, but I can't seem to recall getting one.

The basic idea is what one would expect. For each pair of successive grid points $x_i, x_{i+1}$, the Chebyshev series coefficients $\{c_0,c_1,\dots,c_k\}$ are stored in the InterpolatingFunction[]; then, scaling $z$ to run from $-1$ to $1$ as $x$ runs from $x_i$ to $x_{i+1}$, the polynomial $p(x) = c_0\,T_0(z) + c_1\,T_1(z) + \cdots + c_k\,T_k(z)$, where $T_j(z)$ is the Chebyshev polynomial of order $j$ (ChebyshevT[j, z]), is used to interpolate over the interval between the grid points. It appears that the order of the series $k$ has to be same for each interval in a given InterpolatingFunction[].