checkGravity.py : reduce error tolerance from 1e-3 → 1e-10
@bchareyre I remember that you said once that in our leapfrog implementation the positions are at the Δt/2 while velocities are at Δt. I did not delve (it's on my todo list, see below ;) right now into that code, but by applying your random comment to the starting positions of spheres in this test I was able to reduce error tolerance from 1e-3 to 1e-10. Can you comment on this :)
One day we could think about some other time integration algorithms, like for example PEFRL, Eq.22 in [1]. It has both position and velocity at Δt (not at Δt/2), its Δt would be 4 times larger than Δt in our leapfrog to have similar computational effort (*) but would have 340 times smaller error (see page 14 in [2]) and it has error order of Δt⁴. I'm just mentioning this if someone is interested.
That's an interesting research direction, now that we have a very nicely working high precision implementation it would make sense to eliminate time integration errors :) Otherwise using float128 with leapfrog algorithm that has an error in the order of Δt² necessitates an extremely small Δt ;) Though Δt⁴ error could also be still too large for high precision calculations. More research needed here.
The [1] uses a very nice approach of employing exponential propagators, the same approach which is used in quantum mechanics time integrators. And this makes me wonder because among QM time integrators there is a Kosloff-Chebychev [3][4] method which has error always guaranteed to be smaller than numerical precision (regardless of how small the ULP is [5]). If Kosloff-Chebychev method, which works with exponential propagator, could be converted via approach demonstrated in [1] to classical mechanics then we would also have a classical time integrator with error smaller than numerical precision [5]. That's certainly an interesting thing to do. I am now finishing the QM branch with this algorithm.
(*) PEFRL evaluates force four times in one Δt iteration increment. If we wanted to allow changes in contact detection (and drastic force changes) in these mid-time-steps then some extra mid-steps and knowledge that only "every fourth" iteration one can read position and velocity of particles. Huh. But contact detection with these mid-step positions might require extra work e.g. to determine to which time these positions correspond exactly. So with PEFRL I guess that at first it will be simpler to detect contact every four force evaluations.