Numerical stability in continuous physics simulation [closed]

Pretty much all of the game development I have been involved with runs afoul of simulating a physical world in discrete time steps. This is of course very simple, but hardly elegant (not to mention mathematically inaccurate). It also has severe disadvantages when large values are involved (either very large speeds, or very large time intervals).

I'm trying to make a continuous physics simulation, just for learning, which goes like this:

time = get_time()
while true do
new_time = get_time()
update_world(new_time - time)
render()
time = new_time
end

And update_world() is a continuous physical simulation. Meaning that for example, for an accelerated object, instead of doing

object.x = object.x + object.vx * timestep
object.vx = object.vx + object.ax * timestep -- timestep is fixed

I'm doing something like

object.x = object.x + object.vx * deltatime + object.ax * ((deltatime ^ 2) / 2)
object.vx = object.vx + object.ax * deltatime

However, I'm having a hard time with the numerical stability of my solutions, especially for very large time intervals (think of simulating a physical world for hundreds of thousands of virtual years). Depending on the framerate, I get wildly different solutions.

How can I improve the numerical stability of my continuous physical simulations?

-

closed as too localized by Tetrad♦Dec 24 '12 at 20:18

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

You're assuming that acceleration is constant across a whole timeslice, which is why changing the size of the timeslices gives you different results. – Trevor Powell Dec 6 '12 at 3:46
@TrevorPowell: In this example, acceleration is constant, so that's not the issue here. What I'm saying is that if you run this the continuous program 1 million times for deltatime==1 you will get a slightly different answer than if you run it once with deltatime==1000000 because of numerical stability problems. The actual answer should be exactly the same. – Panda Pajama Dec 6 '12 at 3:49
@TrevorPowell I'm not comparing the discrete solution to the continuous solution. I'm comparing different continuous solutions for different framerates. That's where numerical stability becomes a problem. – Panda Pajama Dec 6 '12 at 3:56
Floating point precision varies with the size of the number you are representing. For (way) more information see: Float Precision–From Zero to 100+ Digits. (The whole series can be found here on Bruce Dawson's blog.) – Eric Dec 6 '12 at 8:12
Never mind. I think I made a mistake writing the question. However, I found my exact question over here: stackoverflow.com/questions/11265081/… so feel free to close this question. – Panda Pajama Dec 6 '12 at 9:49

Plain C code:

#include <stdio.h>
int main(int argc, char **argv)
{
int i;
double manyX = 0.0;
double oneX = 0.0;
double vx = 0.0;
double ax = 0.5;

double oneTimeSlice = 100000.0;
double manyTimeSlice = 1.0;

// calculate one big 1,000,000 timeslice
oneX += oneTimeSlice * vx + (0.5 * ax * oneTimeSlice * oneTimeSlice);

// calculate 1,000,000 1-unit-long timeslices
for ( i = 0; i < 1000000; i++ )
{
manyX += manyTimeSlice * vx + (0.5 * ax * manyTimeSlice * manyTimeSlice);
vx += manyTimeSlice * ax;
}

printf("One Time Slice: %f\n", oneX);
printf("Many Time Slices: %f\n", manyX);
return 0;
}

output:

One Time Slice: 2500000000.000000
Many Time Slices: 250000000000.000000

I'm not seeing the difference you're claiming.

I think you're more likely to be having simple floating point precision problems than numerical stability problems, under the constraints you've specified.

-
Okay, that was not a good example. My program is much larger, but I'm definitely having floating point precision problems (I thought it was the same as "numerical stability"?), and I would like to know if there is a way to mitigate them, or have them somewhat cancel each other. – Panda Pajama Dec 6 '12 at 7:55