ODE Solver / Adaptive stepsize ... how?

Started by
5 comments, last by orbano 19 years, 2 months ago
Hi there! I have finally implemented a 4th order Runge-Kutta approximation in my physics simulation. My step function looks something like this: Step(ObjectArray,ObjectCount,Stepsize) { For each object { Collect forces acting on the actual object (from global variables); Calculate displacement with RK4 (or something similar); Set new object state; } } ---------- The probles is, i really have no idea, how adaptive stepsize comes in. Now in my simulation i use no framerate limitig, physical rendering is called each time graphics rendering is called. I want my simulation to run deltaT = 1.0 in each second. If i use framerate limiting, its quite easy to achieve this, but how should i do it with adaptive stepsizing? Could anyone help? some of my questions: -should i use different stepsize for each object to reach de desired final deltaT (and set the desired deltaT according to the best approximation error)? -should i compute the worst approximation error, and adapt the whole simulation to that? please help :)
"Knowledge is no more expensive than ignorance, and at least as satisfying." -Barrin
Advertisement
1. If your force is fixed during a step you don't need RK4. What do you mean by collect forces? Values or functions?
2. When to do collision detection?

You need some kind of error metric which you calculate every ode-step and from which you decide to adapt your stepsize. You store the stepsize between gameloop-steps.

If you adapt your stepsize per object:
instead of Calculate displacement with RK4 (or something similar);
you do something like:
while(rk4.stepsize()<dt) {
dt-=rk4.step(); // returns actual stepsize, updates predicted maximum stepsize for next step
}
rk4.step(dt);
You need an adaptive stepsize to get a solution with an error under some threshold while reducing calculations.

Perform this algorithm to find the best stepsize:

1)x0 = ?y0 = ?H = ?2)calculate: e(x0 + H, H)e(x0 + H, H/2)h3)if H/h is much greater than 2, goto 4)else goto 5)4)H = 2*hgoto 2)5)x0 = x0 + Hy0 = e(x0 + H, H/2)H = 2*hfinished this step


The function e(x0 + H, H) is an estimation for y, where the second parameter is the current stepsize. So e(x0 + H, H) is one big step, while e(x0 + H, H/2) are two half steps. It is used to approximate the error. You can use RK4 for those functions.

While H is the currently selected stepsize, h is used to decrease it until the error goes below some threashold:

H / h = ( 2^p/(2^p - 1)  *  abs( e(x0 + H, H) - e(x0 + H, H/2) ) / epsilon )^(1/(p+1))




Here epsilon should not be smaller than K*eps, where eps is the machine accuracy and K is:
K = max { abs(y(x)) | x0 <= x <= x0+h }


Good luck, its not as complicated as it looks.



Edit:
You may also want to look for Runge-Kutta-Fehlberg.
nmi! thank you for your description, but i know, how to approximate the error, the only thing i dont know is how to deal with it :)
Im still trying to understan what Trap is saying :) My english is not so good...
"Knowledge is no more expensive than ignorance, and at least as satisfying." -Barrin
1. yes, i need RK4 even with fixed forces/step. its much more accurate. ive tested it.
2. collision detection is after the Step. it doesnt interfere with ODE solver (yet).

I can calculate the error metric you said, but dont know what to do with it :)

Bit i dont really understand this:

while(rk4.stepsize()<dt) {
dt-=rk4.step();
}
rk4.step(dt);

step(); - steps with its own estimated required stepsize to maintain accuracy?
step(dt); - what does this do? it steps back if stepsize exceeds dt?

BTW: how does collision detection come into play? now i use a very simple collision detection:

for each object(a)
{
for each object(b)
{
if a intersects b and a and b are closing, then do the collision.
}
}

//of course im not doing O(n^2) check, im using AABB trees...//
"Knowledge is no more expensive than ignorance, and at least as satisfying." -Barrin
If forces are fixed during a step you can explicitly solve the ode: x_new=x_old+v_old*dt+0.5*a*dt*dt, v_new=v_old+a*dt
Way better than rk4 :)

If your adaptive stepsize is bigger than dt, you just do a step of dt length. If it is smaller than dt you step forward with the adaptive stepsize till the remaining step is smaller than your adaptive stepsize and do a final step of dt length.
oh yes got it :)) just didnt look at the loop...

i think your equation is worse (its something like RK2). but tonight i will sit down and calculate the equations on paper... but the tests showed, that RK4 need ten times larger stepsize than Midpoint method, to achieve same accuracy... we only have to sit down and calculate the errors...
"Knowledge is no more expensive than ignorance, and at least as satisfying." -Barrin

This topic is closed to new replies.

Advertisement