scott wrote:

> Firstly, I wouldn't use integration like this, it's called Euler 
> integration and is basically rubbish and well known for being unstable.  
> A far better method, and the de-facto standard in real-time 
> applications, is Runge-Kutta 4, or RK4 for short.  It is far more stable 
> and much more accurate per CPU time than the Euler method.

True indeed. Now, if I could just figure out how the hell to implement 
RK4...

[Sure, I have a dictionary of mathematics, and there's Wikipedia and 
there's MathWorld, but all involve pretty abstract mathematics that is 
difficult to figure out how to apply to a really complicated 
differential euqation.]

> There are numerous tutorials and explanations online of how to use RK4 
> in exactly the sort of simulation you are doing, but one of my 
> favourites is here, it's very clear.  It's got small code snippets in 
> C++ but they are very simple and should be easy enough to convert to any 
> language.
> 
> http://www.gaffer.org/game-physics/integration-basics

Mmm, that's a nice article. I'm still having trouble following it 
though. Perhaps a Haskell example?

   integrate t dt a (v,p) = (v + dt*(a t v p), p + dt*v)

OK, well that seems trivial enough. So how about RK4?

   evaluate (v,p) t dt (dv,dp) =
     let
       v' = v + dt*dv
       p' = p + dt*dp
     in  (v', a (t+dt) v' p')

   integrate t dt a (v,p) =
     let
       a = evaluate (v,p) t dt -- er... where's the rest?
       b = evaluate (v,p) t (dt/2) a
       c = evaluate (v,p) t (dt/2) b
       d = evaluate (v,p) t dt     c
       dv = (fst a + 2 * (fst b + fst c) + fst d) / 6
       dp = (snd a + 2 * (snd b + snd c) + snd d) / 6
     in (v + dt*dv, p + dt*dp)

Er... is that actually correct? And what's with the formula for A? I'm 
not really understanding how this is supposed to work...

> Your physics time-step is too big.

Yep, I've seen this happen lots before.

> There are several ways to fix this (I would do all of them if I were you):
> 
> 1) Use RK4 instead of Euler (this will let you get away with much bigger 
> time-steps without explosion)
> 2) Use smaller time-steps, eg update the physics 10 times for every 
> graphics update
> 3) Limit the maximum force artificially in the code to avoid explosions

Indeed.

For my chaos pendulum simulations, I added a constant term to the 
denominator in the "inverse square law" force calculations to prevent 
division by zero and absurdly large forces.

I intended (but never got round to) implementing a system where the 
computer integrates a time step, then integrates it again with twice the 
resolution, and if the results differ too much, it subdivides further...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message