damian-666
commented
2 years ago because this in an inherent issue with floating points numbers , each machine and or language implements it differently, for orbits it will barely ever matter esp if you can wrap the c engline and pass params from python script via s standard serializer, but not execute the integrator in phyton, , but with chaos, black holes, flybys, collisions, fluids- like motion, any sort of jerky interaction, will happen right away regardless.. not deterministic, and not reversible at all. also the timestep is crucial in explicit integrators, and no way to predict the maximum in advance.. one might see.. https://www.josstam.com/reversible where its immediately demonstable... his method is not unconditionally stable, it does require a CFL condition, ( max time step) and it does require using of fixed point., but can be used with existing integrators that use floats, and doesnt need to take snapshots of state , just start reversing with -dt..
this might help fix the issues wihtout overhauling the whole thing... https://arxiv.org/abs/2207.07695 its just a preprint but from Jos Stam... here is an excerpt: ( hope it might help, i his contributioins are legendary, im just a game maker using robotics)
Enter the Fixed Point numbers. This format was popular before floats were implemented in hardware. They work well if your data is bounded. Since our integrators are bounded for reasonable time step sizes h we can assume without loss of generality that our reals lie in the interval [1,−1]. We can then approximate a real number R by the following integer I: I = bR × BIG_INTc, (10) where BIG_INT is a large integer and b·c is the integer part of a real number. The only operation needed in our integration algorithm is addition as shown below in the implementation. To get a real number from a fixed number, we simply “invert” Eq. 10: R = I / BIG_INT. Actually, our implementation is hybrid and uses floating points to compute forces and fixed points when integrating. 4 Implementation In the following C++ code only the integration steps use the fixed format. Both the coordinates and the velocities are represented in fixed format. The time step can be any arbitrary float. The simulation is completely controlled by the force function which uses only floats. Existing simulations can therefore easily be implemented in this framework. We represent our fixed points as 64 bit integers “int64_t”, a type standardized since C++11. t y p e d e f int64_t f i x e d ; c o n s t f i x e d c_max_fixed = 0 x1000000000000000 ; s t a t i c f l o a t x 2 f ( c o n s t f i x e d i_f ){ r e t u r n i_f / ( f l o a t ) c_max_fixed ; 6 } s t a t i c f i x e d f 2 x ( c o n s t f l o a t i_f ) { r e t u r n ( f i x e d ) ( i_f c_max_fixed ) ; } s t a t i c f i x e d x , v ; s t a t i c i n t s i z e ; // d e f i n e your f a v o r i t e f o r c e f i e l d h e r e . s t a t i c f l o a t f o r c e ( c o n s t i n t i , c o n s t f l o a t i_x ) ; v o i d i n t e g r a t e ( c o n s t f l o a t h ) { c o n s t f l o a t h2 = 0 . 5 f h ; f o r ( i n t i = 0 ; i < s i z e ; i++) { x [ i ] += f 2 x ( h2 x 2 f ( v [ i ] ) ) ; } f o r ( i n t i = 0 ; i < s i z e ; i++) { f l o a t f = f o r c e ( i , x 2 f ( x [ i ] ) ) ; v [ i ] += f 2 x ( h f ) ; } f o r ( i n t i = 0 ; i < s i z e ; i++) { x [ i ] += f 2 x ( h2 x 2 f ( v [ i ] ) ) ; } } We also implemented the integrator in JavaScript (JS) to create a web based demo that is easily shared. Ironically, we were faced with the problem that JS does not differentiate between integers and float numbers. It only has the monolithic Number type. We found a solution by treating fixed point numbers as strings. We lost some efficiency of course. The implementation is as follows: v a r XN = 1000000000; f u n c t i o n X2F( x ) { r e t u r n (+x ) / XN; } f u n c t i o n F2X( f ) { ( ( f XN) | 0 ) . t o S t r i n g ( ) ; } f u n c t i o n XADD( x1 , x2 ) { (+x1 + +x2 ) . t o S t r i n g ( ) ; } The code is a bit obfuscated but JS programmers will recognize the hack to get an integer from a Number by taking the bitwise “or” with zero. The “+” operator transforms a string into a Number. See the accompanying HTML file from the link given below for the entire implementation. We would like to hear from readers who have a better solution. Note that this is just a quirk of the JS language. For most typed languages, like in our C++ implementation above,
dmhernan
hannorein
DorianSchuylerAbbot
Hi, I've been doing some tests, and there's some strange behavior. I was trying to test WHFast's time symmetry. With the C interface, I got the expected behavior. I set up some initial conditions, and integrated for some timestep dt and time t. I reversed the sign of the timestep dt -> -dt, and again integrated for t. The error is close to machine precision. WHFast is time-symmetric, as expected. But then, the behavior was different with the python version. After reversing the timestep, the energy error magnitude didn't decrease to machine precision. I did confirm the time goes back to 0, and used (for C and python) sim.steps() rather than sim.integrate(), to ensure we arrive exactly at the origin. I'll update the thread if I find something new, and happy to post more specific code if that helps.
Also, a question/suggestion: I think an index of all Rebound commands and functions would be helpful. This is probably available somewhere, but I couldn't find it.