# How is the Chebyshev method used by JPL?

In the DE packages, NASA gives us a series of coefficients for the Chebyshev approximation. As I understand it, those approximations are given by a series of polynomials $T_n$ for the interval $[-1,1]$. Then we consider $f(t)$ to be $\sum_i a_i T_i(\tau)$ (where $\tau$ is the normalised time). What is $a_i$ and what are its values?

Then we consider an approximation using a Taylor series defined as $f(t)=\sum_i b_i\tau^i$. What is $b_i$? Is it a Bessel function?

The last two relations represent the Clenshaw algorithm. Here I don't exactly understand how it works. Also does the number of 13 coefficients have something to do with Runge's Phenomenon?

Below are Eqs 3.52 through 3.59, from section 3.3.3 (Chebychev Approximation) in the book Satellite Orbits: Models, Methods and Applications* by De Oliver Montenbruck and Eberhard Gill Springer 2000, ISBN 978-3-642-58351-3 which is viewable in books.google.com    • It would be great if you could add a link to the source that you are reading. – uhoh Sep 9 '18 at 3:29
• books.google.com/… – Alexandru Lapusneanu Sep 9 '18 at 10:05
• That's a fantastic book by the way! – uhoh Sep 9 '18 at 10:46

The $$a_i$$ are the coefficients for the Chebyshev approximation. As you say, NASA gives us those. That is what you find in the DE files, e.g. de430.bsp. (Don't click on that unless you want to download a >100MB file.) NASA/JPL needed a way to provide high resolution and high accuracy functions of time for the positions of the planets, and the most compact way was coefficients of Chebyshev polynomials. The functions are broken up into small intervals, within which the coefficients and resulting functions are valid.

No, the $$b_i$$ are not Bessel functions. They are the coefficients of a Taylor series that estimates $$f$$ over the $$\tau$$ for that interval. That equation was taken out of context, where that context was discussing why Chebyshev polynomials, $$T_n\!\left(\tau\right)$$ are used for the approximation instead of a Taylor series, which uses the polynomials $$\tau^n$$. The reason is that a Chebyshev approximation will require fewer terms for the same accuracy.

The Clenshaw algorithm is simply a way to both generate the Chebyshev polynomials and to multiply by the $$a_i$$ and add them up, all at the same time to minimize the number of operations required. It is well described in the Wikipedia page linked, and can be easily derived from the recurrence relation for the Chebyshev polynomials near the top of your image.

• The Planetary and Lunar Ephemerides DE430 and DE431 – uhoh Sep 9 '18 at 5:20
• The error in a Chebyshev approximation is not flat. It oscillates and is close to that of the minimax polynomial of the same order. The difference between a well-behaved function $f(t)$ and the minimax polynomial of order $n$ that approximates this function over that interval exhibits $n+2$ extrema, where $n$ is the order of the polynomial, with the absolute value of the difference at each extremum being the same. There's no guarantee of how the error function behaves for a Taylor series, or even for a least square error approximation. ... – David Hammen Jan 23 '19 at 22:48
• A user of an approximation typically cares about the worst thing the approximation can do to them rather least square behavior. That Chebyshev polynomials are close to the ideal minimax polynomial is why they're used so widely. – David Hammen Jan 23 '19 at 22:50

The article Format of the JPL Ephemeris Files has a detailed breakdown of how to use the Development Ephemeris including an example walkthrough and example source code.

I think you have a lot of superfluous information, all you really need is:

$$\sum_{i=1}^{n}a_i T_i(x)$$

Where $$a_i$$ are the coefficients, $$x$$ is the time variable normalized to the interval $$[-1,1]$$, and $$T_i(x)$$ is defined as follows:

$$T_0(x)=0$$

$$T_1(x)=x$$

$$T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)$$

Code in Javascript to perform this computation:

function computePolynomial(x,coefficients){
let T=new Array();

T=1;
T=x;
for(let n=2;n<coefficients.length;n++)  {
T[n]=2*x*T[n-1] - T[n-2];
}

let v=0;
for(let i=coefficients.length-1;i>=0;i--){
v+=T[i]*coefficients[i];
}
return v;
}


Note the summation of all of the variables is done in reverse order, from smallest to largest to avoid floating point rounding.

The Github repository gmiller123456/jpl-development-ephemeris has non-optimized source code in several languages which implement the whole process.