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.