# Why do we use Chebyshev polynomials to predict positions of solar system bodies?

I've seen information for how to use the Chebyshev coefficients that JPL publishes, but nothing that explains why. Why isn't a simpler mathematical model used?

I'm also surprised that JPL uses a numerical model rather than doing a n-body physics simulation of the Solar System.

• Welcome to Space! search results for @DavidHammen Chebyshev But I'm confused by your last sentence, because JPL uses a numerical model to do n-body physics simulation of the Solar System. There's no "rather than" here.
– uhoh
Jan 20 at 1:38
• Thank you! Re: my last sentence, I was under the impression that JPL generated the ephemeris by examining historical observed data and extrapolating, rather than using the known masses of the bodies and applying Newton's laws. Thanks for correcting me. Jan 20 at 4:17
• See en.wikipedia.org/wiki/… for a good summary of how JPL produce their ephemerides. Briefly, they integrate the equations of motion, with relativistic corrections. But that requires good data for the body masses and initial locations & velocities. So the generated ephemerides are verified against ground- & space- based observational data. The Chebyshev coefficients are simply the method used to store the generated ephemeride data so that it can be precisely interpolated as necessary. Jan 20 at 7:03
• @prideout Oh, I see what you mean. What's done is a fitting procedure; many many parameters (masses and higher order gravitational terms, and initial conditions) are allowed to vary in order to find the best fit between the numerical model and hundreds of years of observations. See the documentation for the latest DE at The JPL Planetary and Lunar Ephemerides DE440 and DE441 (found here)
– uhoh
Jan 20 at 10:39
• Note that Chebyshev series are commonly used to construct fast, precise, numerically stable approximations to functions that are otherwise difficult and expensive to calculate. Ephemerides are merely one application of many. Jan 21 at 18:19

I'm also surprised that JPL uses a numerical model rather than doing a n-body physics simulation of the Solar System.

JPL does use an n-body physics simulation of the solar system. It's a very complex and very time consuming simulation. Then they run it over and over and over, each time adjusting various elements of the simulation (e.g., planet masses, position and velocity at epoch time) so as to decrease the mismatch between computed versus recorded observations. It is only after running the simulation many times over that they generate the Chebyshev coefficients that are released as a new version of the Development Ephemerides.

Why do we use Chebyshev polynomials to predict positions of solar system bodies?

Suppose instead of releasing sets of interpolating coefficients that JPL released all of those magic numbers: masses, epoch time, and states (positions and velocities) at the epoch time. To re-create the states as propagated in the JPL simulation, one would need to run the exact simulation JPL uses (which JPL does not release) using

• The same type of computer,
• The same operating system,
• The same compiler / linker / system library, and
• The same compiler / linker options

that JPL uses for running its simulation. Any deviation from that sameness will result in a simulation run that diverges non-linearly from the JPL results. To make matters worse, running that simulation will be computationally expensive. Releasing an ephemeris in terms of sets of interpolating coefficients avoids this entire mess.

JPL uses Chebyshev polynomial coefficients as opposed to other interpolation schemes in its released ephemerides for many reasons. One reason is that calculating position and velocity from those Chebyshev coefficient sets is rather fast. The underlying software needs to determine the time block that contains the relevant coefficients, load the coefficients, and apply them.

Another reason for using Chebyshev polynomials is that a Chebyshev polynomial fit comes rather close to a magical polynomial fit called the minimax polynomial. The minimax polynomial minimizes the maximum absolute deviation from the true curve. This is exactly what a consumer of a fit wants. I, for one, care less about the goodness of a fit in a least squares sense than I care about the worst case performance. The minimax polynomial minimizes worst case performance.

I don't know whether JPL uses it, but there's a complicated and computationally expensive algorithm, the Remez exchange algorithm, that tweaks the Chebyshev coefficients even closer to the magical minimax polynomial.

• FWIW, here's a plot comparing the errors in approximating $\frac{3x^2}{5-2\sqrt[3] x}$ over $[0, 1]$ with 4th degree Chebyshev & Remez minimax polynomials, see here for the poly coefficients. Jan 20 at 13:24
• That's a nice plot, @PM2Ring. It shows several key characteristics of the Chebyshev and minimax polynomials. (1) They're fairly close. A little bit of Remez exchange and voila: there's the minimax polynomial. (2) The minimax polynomial has six error extrema, each of which is exactly the same distance from the zero error x axis. (3) The Chebyshev polynomial is "better" over about 80% of the domain. (4) Murphy's law dictates that my randomly chosen points will inevitably fall in the 20% portion of the domain where the minimax polynomial is better. Jan 20 at 14:17
• This is a great answer, thanks so much! Maybe part of the reason for why I was finding it unintuitive is that position is a 3D quantity but many introductory explanations (e.g. the plot from @PM2Ring) are 1D. Obviously these approximation techniques extend naturally to multiple dimensions, so it's just a mental block. Jan 20 at 17:03
• @prideout The Chebyshev coefficients in the JPL Development Ephemerides are for one dimension as a function of time. Each planet has three sets of coefficients in a JPL Development Ephemerides, one for x position as a function of time, another for y position, and a third for z position. Jan 20 at 19:12
• @uhoh There's no particular need for velocity coefficients. The Chebyshev polynomials are functions of (scaled) time. It's a trivial matter to use the same coefficients used for the x, y, & z position elements to compute the time derivatives of those elements. That said, the SPICE SPK system does provide the ability to have separate coefficients for position and velocity. A type 2 kernel contains position only coefficients while a type 3 kernel contains position and velocity coefficients. The Development Ephemerides contain what are essentially type 2 kernels. Jan 22 at 3:29