I wonder how can authors claim such great certainty in their results, after all, it is all based on number crunching. Some floating point error, and similar is bound to creep in...
The numerical methods used for this kind of calculations are engineered to compensate numerical errors (floating point and errors inherent to integrators).
They take advantage of a priori knowledge about the laws of physics, in particular the conservation of mechanical energy and angular momentum. Predictor-corrector is one family of methods. Other methods rely on convergence as timesteps change.
There is a lot of literature on numerical integration applied to celestial mechanics and new methods being released every year. These methods are tailored to this problem space.
When I looked into making a solar system simulation awhile back, I actually found it shocking how much goes into monitoring and ensuring conservation of energy/momentum in these models.
Once I was aware of this issue though, even more shocking to me was that many papers on climate models I came across do not even mention they monitored adherence to conservation laws. There is a disconnect there, I would think this adherence either is a big deal (as would be suggested by the practice in astronomy) or not (as would be suggested by climate research), not that it would change by subfield.
When just doing point gravity dynamics (no friction, collision etc) energy and momentum take care of themselves as long as the math is clean. I say this confidently because I made a solar system simulation[1] with no attention to conservation of energy at all. In practice Star systems also shed large amounts of angular momentum through solar wind[2], but I expect this is often not modelled.
Climate models are an entirely different matter and with respect, there is a disconnect with your scepticism of them for this matter.
[1] https://strainer.github.io/fancy#00
[2] http://farside.ph.utexas.edu/teaching/plasma/Plasmahtml/node...
It says in your readme[1] though that your solar system simulation simulation is off by 70k km (regarding Earth orbit) after only 1 year relative to the JPL ephemerides, for "unknown reasons". I don't see how you can use that to justify not checking for conservation of momentum/energy.
I also do not understand your final sentence.
Come on it doesnt say "unknown reasons" it says the "It is possible the discrepancy is due to the limitation of javascripts 64bit FP numbers rather than subtle algorithmic or physics error"
70k km / a year (on a billion km long orbit) was a lot closer than I had hoped for, considering its a 64 bit newtonian model. But this is beside the fact that having worked on the project a fair bit, I understand that any model adjustments to conserve energy or momentum are non-physical 'fudges' to correct deficiencies of the model. There are no laws of physics dedicated to adjusting energy, momentum, information in order to conserve them. They are conserved by a mathematical correctness of every physical law - a proper account of that eludes me, but it is certainly explained somewhere.
Right, you don't know how much error the floating point error accounts for. You also brought up loss of angular momentum to solar wind, etc. These are all just speculations. That may be fine for the purposes of your simulation, but not for predicting if asteroids will impact the earth or the influence of human activity on Earth's climate.
Also, I know you are supposed to use a symplectic integrator to ensure conservation of energy when doing these simulations, so it is not somehow hardcoded into the laws, eg: https://en.wikipedia.org/wiki/Leapfrog_integration
We could certainly calculate how much error was introduced by rounding, if we had the time and inclination. It is useful to understand that about numerical modelling, it should not be necessary to accept any indeterminable discrepancies.
Conservation of energy, momentum, information is hardcoded into all known laws of Physics. It is our integration algorithms which do not all ensure conservation of information - but that is carefully achievable with basic first order schemes, such as Verlet and leapfrog integration.
>"Conservation of energy, momentum, information is hardcoded into all known laws of Physics. It is our integration algorithms which do not all ensure conservation of information - but that is carefully achievable with basic first order schemes, such as Verlet and leapfrog integration."
I am not sure about "conservation is encoded" part, but I can believe it. Either way, in practice the simulations will encounter issues along these lines if they are not careful. For example, I left your simulation running for what amounted to ~120 yrs and saturn lost all her moons. This is the best pic I could get, sorry (hopefully it is reproducible): https://i.imgur.com/C7emd28.png
Since you do not report anything about conservation of energy/momentum, how do we know that isn't the problem?
>For example, I left your simulation running for what amounted to ~120 yrs and saturn lost all her moons.
Thanks! heh, it is rather complex there are all sorts of things... most likely the new gravity function set for that model messed up, which can apply attraction between groups of objects if they are distant and weak enough.
But honestly I am certain that it is something in particular and not something mysterious about conservation of energy/momentum which i have not come across yet, because there is no formula in that gravity simulation which should be able to destroy momentum. In the other models, there are quasi-physical friction and pressure functions which do destroy momentum and do have some rough adjustments to compensate for that somewhat, but accuracy is lost whether things are compensated for or not. Tracking energy and momentum can be a good way for tracking accuracy of a model, but if you do any adjustments to correct it, that is a kind of fudge not based in true physics. It is poignant that you expressed scepticism of climate models because of a perceived lack of such adjustments, while the opposite is actually true. Climate models have to include a lot of calibrated adjustment and quasi-physics, because the whole system is too complex to represent at all relevant scales. This should not condemn the work of computer scientists specialised in climate modelling to undue scepticism.
It's a well established fact that numerical integration methods either gain or lose energy (in particular, the Runge-Kutta family is known to lose energy over time). For celestial mechanics simulations, a special class of numerical integration methods called "symplectic integrators" are used and their purpose is to conserve energy and angular momentum.
> but if you do any adjustments to correct it, that is a kind of fudge not based in true physics.
When you are numerically integrating differential equations that model physical phenomena, you're not doing "true physics" but an approximation thereof.
And an approximation that makes Earth drift 70 km per year or Saturn's moons drift out of orbit in a few hundred years is a very bad approximation by scientific standards.
The methods used for celestial mechanics calculations need to be precise over thousands to millions of years. And the way they work is to "fudge" with the numerical methods to preserve energy and angular momentum. It's a much better approximation of "true physics" than your toy simulation.
Your assumption that the issues are due to floating point errors is incorrect. 64 bit double precision is millimeter accurate to the orbit of Neptune. That's good enough for scientific applications.
If you're interested in this, you could take a look at this scientific grade N-body simulator and the methods it uses: https://github.com/hannorein/rebound