Fixing rounding errors

[Neo-13]

I am trying to model planetary orbits, and have hit a snag. I get the orbits almost perfect, but every loop, it gets slightly out, which, I have been told, is due to rounding errors, which would make sense.

I am further told that the rounding errors lie where I work out the distance between the bodies. I am using Pythagoras' theorem:

sqrt(dx^2 + dy^2), where dx and dy are delta x and delta y respectively.

Apparently this is likely to give me rounding errors. Bare in mind that this is done 30 times every seconds, for an indefinite period of time, so even the smallest rounding errors become significant.

Thanks in advance for any advice.

[dELtaluca]

It is most likely not due to rounding errors, but due simply to the fact that numerical integration is not perfect, because you cannot use an infinitesemally small time step in your simulation, you are always going to get innacurate results, different methods of numerical integration will give you differeing results, some are more stable, some will conserve energy in physical simulations better, but none the less they will never be perfect.

[dELtaluca]

To add, rounding errors will, over a long enough period of time manifest themselves in anomalous results, but the fact you said it happens every time it loops to a notable degree tells me it is much more likely to be due to the numerical integration being too inaccurate

[Neo-13]

At 9/21/09 01:53 PM, dELtaluca wrote: To add, rounding errors will, over a long enough period of time manifest themselves in anomalous results, but the fact you said it happens every time it loops to a notable degree tells me it is much more likely to be due to the numerical integration being too inaccurate

Hmm, OK. Well, there is surely a way to resolve that, as I have seen an example of a simulation that does perfect orbits, even when you fiddle with it.

Have you got any tips on how to deal with inaccuracies?

Thanks.

[dELtaluca]

As i said, use a better numerical integrator and or make your time steps smaller.

A simple example, whilst you might be tempted to treat acceleration as a constant over the time step: aka:

s += vt + 0.5at^2
v += at

whilst this may give you perfect accuracy when modelling freefall, since acceleration is a constant, using simpletic euler integration is much better when modelling varying forces as it will conserve energy a lot better.

v += at
s += vt

i have even found that simpletic euler can give better results than RK4 (Runge kutta integration) when modelling vastly varying accelerations, such as a body orbiting a sun in a highly elliptical orbit

[dELtaluca]

you might also consider using a variable time step, so that as accelerations become larger, you decrease the time step to increase accuracy, doing so means that as accelerations become larger you do smaller time steps so that it remains accurate, whilst using larger time steps when accelerations are smaller to avoid unnecesary accuracy :P

[Neo-13]

At 9/21/09 02:08 PM, dELtaluca wrote: As i said, use a better numerical integrator and or make your time steps smaller.

A simple example, whilst you might be tempted to treat acceleration as a constant over the time step: aka:

s += vt + 0.5at^2
v += at

whilst this may give you perfect accuracy when modelling freefall, since acceleration is a constant, using simpletic euler integration is much better when modelling varying forces as it will conserve energy a lot better.

v += at
s += vt

i have even found that simpletic euler can give better results than RK4 (Runge kutta integration) when modelling vastly varying accelerations, such as a body orbiting a sun in a highly elliptical orbit

Have you got any useful links on Euler integration? I haven't met it before.

Also, at the moment, I work out the velocity like this:

vx += a*cos(theta);
vy += a*sin(theta);

planet.x += vx;
planet.y += vy;

I leave out t because I'm looping this code (I assume that's fine).

Theta is the angle to the horizontal, i.e. I am resolving a into its cartesian components.

[dELtaluca]

At 9/21/09 02:29 PM, Neo-13 wrote: vx += a*cos(theta);
vy += a*sin(theta);

... that is not what your acceleration calculations should look like.

given constant G, mass of body acting on the current body m, the acceleration should be calculated something like this:

var dx:Number = p1.x - p0.x;
var dy:Number = p1.y - p0.y;
var sc:Number = G*m*Math.pow(dx*dx+dy*dy,-1.5);

var ax:Number = dx*sc;
var ay:Number = dy*sc;

taking the squared magnitude of the vector to -1.5, will return the value of 1/r^3, r being the length of the vector, which is a composition of 1/r to normalise the vector, and 1/r^2 which is the scaling value from the gravitational acceleration formula a = -GM/r^2

from then, you should have a time step 'dt' over which your integration takes place whicih for simpletic euler integration would be:

p0.x += (v0.x += ax*dt)*dt;
p0.y += (v0.y += ay*dt)*dt;

[Neo-13]

At 9/21/09 02:35 PM, dELtaluca wrote:

At 9/21/09 02:29 PM, Neo-13 wrote: vx += a*cos(theta);
vy += a*sin(theta);

... that is not what your acceleration calculations should look like.

given constant G, mass of body acting on the current body m, the acceleration should be calculated something like this:

var dx:Number = p1.x - p0.x;
var dy:Number = p1.y - p0.y;
var sc:Number = G*m*Math.pow(dx*dx+dy*dy,-1.5);

var ax:Number = dx*sc;
var ay:Number = dy*sc;

taking the squared magnitude of the vector to -1.5, will return the value of 1/r^3, r being the length of the vector, which is a composition of 1/r to normalise the vector, and 1/r^2 which is the scaling value from the gravitational acceleration formula a = -GM/r^2

from then, you should have a time step 'dt' over which your integration takes place whicih for simpletic euler integration would be:

p0.x += (v0.x += ax*dt)*dt;
p0.y += (v0.y += ay*dt)*dt;

I am using F = Gm1m2/r^2 and then a = F/m, which is the same as your equation, plus the force - is that not missing from your one? Forgive my ignorance.
The way I see it is if you combine the two equations you are left with a = FGm/r^2.
Also, the time step, why do I need it since the calculation is looped at a constant rate?

Thanks for all the help.

[dELtaluca]

The point i was making is that rather than using the radial vector as the basis of your acceleration vector, you are calculating the angle of the radial vector, then using sine/cosine to build the vector which is totally uneccesary, and slower and less accurate.

the point of having an explicit time step, is that you can easily change the speed of the simulation, and as per my advice, have a variable time step, rather than having it bound to a single value: 1

[henke37]

There is also the potentional issue of coordinates of displayobjects being exact to exactly a twentieth of a pixel, a twip.

[LeechmasterB]

Sorry to differ delta!

Anyhow, not sure how it is with flash, but i have done physical simulations of orbiting bodies ect. and euler proved to be far off compared to RK4. So i highly suggest trying at least both methods and comparing the error! :)

I got a nasty java example to prove the difference of result. ^^

[LeechmasterB]

And also chose the number to model the physics wisely, very large or very small numbers give you bad results aswell, but it doesn't hurt tuning the numbers to your needs to prevent numerical innacuracy. Just keep em in a nice range...

[dELtaluca]

At 9/21/09 04:58 PM, LeechmasterB wrote: Sorry to differ delta!

You say Euler, i say Simpletic Euler.

They're not the same.

[LeechmasterB]

At 9/21/09 05:12 PM, dELtaluca wrote:

At 9/21/09 04:58 PM, LeechmasterB wrote: Sorry to differ delta!

You say Euler, i say Simpletic Euler.

They're not the same.

Whats the difference?

Anyhow from my experience (and also AAA game developers point of view) its good to have a fixed timestep (preferably 1 as the time multiplication can be left out that way in flash, thats always good...) and using RK4.

Why?

The faster your object travels the more inaccurate euler integration comes, because the gravitation force will change its angle and strength within the interval of the velocity, thats why RK 4 is trying to get more close to the real deal. So basically using 'normal euler integration' will feel like the physics are frozen for the time of the object passing from x to x' which its actually not because of the force working on it all the time throughout.

That explanation is simply to make clear whats the advantage of RK4 to euler, yet i can't see a difference in your formula above and regular euler integration (except for your suggestion to vary dt, which would mean you d also have to change the actual velocity according to dt ect. ect. and that doesn't fix the other problem it even can make it much worse).

To maintain a constant physical simulation framerate there are niftier tricks btw. I posted some code ages ago about an alternative way to fix the timestep compared to 8bitrockets aproach (mine s based on his aswell, but trying to give it a thread like aproach from separating graphics and physics ect.).

[LeechmasterB]

OOH you mean Semi-Implicit-Euler.. i will read up on it, but calling it by its name would have helped ;D

[LeechmasterB]

And now i see why searching didn't show much results on this:

semi-implicit Euler method, also called symplectic Euler :P

[LeechmasterB]

Neo13, you could also have a look at Velocity Verlet integration. VV and SE are claimed to be energy preserving whereas RK4 is not. So giving all 3 a shot and checking which works best might be the way to go..

[dELtaluca]

At 9/21/09 05:33 PM, LeechmasterB wrote: And now i see why searching didn't show much results on this:

semi-implicit Euler method, also called symplectic Euler :P

hehehe, and yeh; i did another test, and RK4 this time did come out slightly ontop in terms of stability.

However, one thing that is interesting to note, is that considering the following situation:

a body orbiting a static sun in a highly elliptical orbit:

using RK4 integration, the body will over time enter a shorter and shorter orbit, as RK4 loses energy over time.

using sympletic integration, the body will over time, have it's perhilion advance, and gain energy at a slower rate than the RK4 approach will gain energy.

http://spamtheweb.com/ul/upload/210909/8 2284_compare_eueler_rungekutta_solar.php

red is RK4, black is sympletic euler, fixed time step, in this situation at whch both approaches fail, i would consider sympletic euler to give the best results, not only because it doesn't explode, but because in reality orbits do have advancing perihelion (the orbit itself rotates) (look at mercury) which although this is not modelling it, it's an indirect result :P

However i do retract my previous statement, RK4 is still in general more stable (if only by a small margin for the massively increased level of computational complexity) than sympletic euler.

[dELtaluca]

although i guess by that same argument you could say RK4 is better, since in reality orbits do decay due to solar wind friction and gravitational radiation etc, but the effects are magnitudes less noticebal than perihillion advancement

[LeechmasterB]

Yeh, to quote my insight:

At 9/21/09 05:47 PM, LeechmasterB wrote: Neo13, you could also have a look at Velocity Verlet integration. VV and SE are claimed to be energy preserving whereas RK4 is not. So giving all 3 a shot and checking which works best might be the way to go..

Leaves Velocity Verlet to be tried aswell, yet.

Still your earlier post was quiet confusing delta, since:

At 9/21/09 02:35 PM, dELtaluca wrote: from then, you should have a time step 'dt' over which your integration takes place whicih for simpletic euler integration would be:

p0.x += (v0.x += ax*dt)*dt;
p0.y += (v0.y += ay*dt)*dt;

To me thats still regular euler integration since simplectic euler would use v n +1 istead of v n whereas v0 is vn... <.<

[LeechmasterB]

At 9/21/09 05:58 PM, dELtaluca wrote: although i guess by that same argument you could say RK4 is better, since in reality orbits do decay due to solar wind friction and gravitational radiation etc, but the effects are magnitudes less noticebal than perihillion advancement

Thats why people add thrusters to satelites/space stations and whatever and it takes fuel to stay up in orbit.

[dELtaluca]

At 9/21/09 06:17 PM, LeechmasterB wrote: Thats why people add thrusters to satelites/space stations and whatever and it takes fuel to stay up in orbit.

Satelites/space stations are usually still skimming the earth's atmosphere so that the atmospheric drag is quite significant, wheras drag from solar wind is practicly negligable unless you're talking about keeping a satellite in orbit about the sun for 10 million years :P in which case the thrusters are more likely being used to counteract the influence of over attractors.

And if you looked more closely you'd see my code snippet IS sympletic, since:

p0.x += (v0.x += ax*dt)*dt;

is equivalent to:

v0.x += ax*dt;
p0.x += v0.x*dt;

which is sympletic, normal euler is:

p0.x += v0.x*dt;
v0.x += a0.x*dt;

p0,v0 here are the labeling for which body i'm referring to in the code.

[Neo-13]

Woah. Now I guess I'll spend the next few days trying to understand all that! :P

[LeechmasterB]

At 9/21/09 06:25 PM, dELtaluca wrote: Satelites/space stations...

^^

And if you looked more closely you'd see my code snippet IS sympletic, since:

Okay i guess it looks okay on paper, we ll just have to trust your sample. Somehow I am too lazy to try it out lol

@neo-13 have fun and good luck! ^^

[Neo-13]

At 9/21/09 06:33 PM, LeechmasterB wrote:

At 9/21/09 06:25 PM, dELtaluca wrote: Satelites/space stations...

^^

And if you looked more closely you'd see my code snippet IS sympletic, since:

Okay i guess it looks okay on paper, we ll just have to trust your sample. Somehow I am too lazy to try it out lol

@neo-13 have fun and good luck! ^^

Thanks guys, even though I don't understand most of it.

*Rolls up Sleeves...*