> http://www.pacificwebguy.com/downloads/TestingRK4.zip

BTW, I get a file not found runtime error with that...

Post a reply to this message

scott wrote:
>> http://www.pacificwebguy.com/downloads/TestingRK4.zip
> 
> BTW, I get a file not found runtime error with that...
> 

Dagnabit, I wish the packaging system would work properly, otherwise I 
don't know what dependencies these executables have...

-- 
...Ben Chambers
www.pacificwebguy.com

Post a reply to this message

Chambers wrote:

> Dagnabit, I wish the packaging system would work properly, otherwise I 
> don't know what dependencies these executables have...

Try distibuting executables dynamically linked against GTK+ ;-)

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

Post a reply to this message

"scott" <sco### [at] laptopcom> wrote:
> If I am guessing correctly, "particles_stepped" includes velocity
> information?  Well you need to keep 4 separate copies of the velocity
> information just like you did with the force information.
>
> >  for (int i=0; i<size; i++) {
> >   forces_final[i] =
> >    (forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])/6.0f;
>
> Then you need to calculate the weighted sum of the velocities in this manner
> too here, and use the weighted velocity in the final step call.

Sorry, but I have trouble believing that without seeing a proof.  Intuitively,
if I am tracking the velocity information, then why do I need the force
information separate?  But if I integrate the forces, why should the velocity
information be separate?  What I want is a final force value that I can use to
modify my velocity value.

Anyway, looking for other sources on RK4, the Wikipedia article at
http://en.wikipedia.org/wiki/RK4 mentions something that was left out here...
that the weighted average of the forces should be multiplied by the timestep.

That would change the above line to
(forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])*dt/6.0f;

This would account for the strange scaling effect I saw.  When I get home
tonight, I'll try it out.

......Chambers
www.pacificwebguy.com

Post a reply to this message

Chambers wrote:
> Here's a test app for Windows:
> 
> http://www.pacificwebguy.com/downloads/TestingRK4.zip
> (Note that, although the app itself is named "Threaded_Physics.exe", I 
> haven't gotten around to implementing the "threaded" part yet.  Right 
> now, I'm just trying to make it work *correctly*).
> 
> Click anywhere on the screen to shoot particles in that direction.
> Press F1 to switch between Newtonian integration and Rk4 integration.
> 
> The main thing to notice, other than the fact that RK4 is *much* slower 
> than Newtonian, is that using RK4 seems to amplify the forces somehow. 
> The result is that everything seems rather... well, crazy.  In contrast, 
> I've reached a point where the Newtonian doesn't look quite so bad 
> (although I'd still like to work on the speed factor a bit).
> 
> Here's the inside of the RK4 integrator, as I implemented it:
> 
> {
>  collide(particles_now, forces_now, gravity);
>  step(particles_now, particles_stepped, forces_now, dt/2);
>  collide(particles_stepped, forces_A, gravity);
>  step(particles_now, particles_stepped, forces_A, dt/2);
>  collide(particles_stepped, forces_B, gravity);
>  step(particles_now, particles_stepped, forces_B, dt);
>  collide(particles_stepped, forces_C, gravity);
>  for (int i=0; i<size; i++) {
>   forces_final[i] =
>    (forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])/6.0f;
>  }
>  step(particles_now, particles_now, forces_final, dt);
> }
> 
> And stubbed headers for the functions:
> 
> collide(data, result, gravity)
> This calculates the collisions between particles, and the resultant 
> forces (gravity is used as the initial, or base, force to start with)
> 
> step(data, result, forces, dt)
> Applies forces to data with a timestep of dt
> 
> PS I just got a collection of all of Chopin's piano music (something 
> like 15 CDs), and it's great "productivity" music (ie, have on in the 
> background while I'm working on something else) :)

The last time I did a physics simulator (springs and weights), I found 
that calculating the change in acceleration, and figuring that into the 
new location and velocity data, greatly reduced the error that accumulates.

Regards,
John

Post a reply to this message

John VanSickle wrote:
> The last time I did a physics simulator (springs and weights), I found 
> that calculating the change in acceleration, and figuring that into the 
> new location and velocity data, greatly reduced the error that accumulates.

That's why I'm calculating forces :)

-- 
...Ben Chambers
www.pacificwebguy.com

Post a reply to this message

>> If I am guessing correctly, "particles_stepped" includes velocity
>> information?  Well you need to keep 4 separate copies of the velocity
>> information just like you did with the force information.
>>
>> >  for (int i=0; i<size; i++) {
>> >   forces_final[i] =
>> >    (forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])/6.0f;
>>
>> Then you need to calculate the weighted sum of the velocities in this 
>> manner
>> too here, and use the weighted velocity in the final step call.
>
> Sorry, but I have trouble believing that without seeing a proof. 
> Intuitively,
> if I am tracking the velocity information, then why do I need the force
> information separate?

Because you are performing two integrations in parallel here.  You are 
integrating the accelerations to get velocities, and you are also 
integrating velocities to get positions.  The way RK4 works is to calculate 
the gradients at certain points and then take a weighted sum to work out the 
most accurate gradient.  You must do this for the gradients of velocity (ie 
acceleration or force) *AND* position (ie velocity).  If you only do it for 
one you are not doing it properly.

> But if I integrate the forces, why should the velocity
> information be separate?

Because you are also integrating the velocities to get position.

> What I want is a final force value that I can use to
> modify my velocity value.

But you also need a final velocity value to modify your position.  This also 
has to be calculated using the RK4 method.  So you need to work out the 
velocity at each step and remember it, just like you did with the forces.

> Anyway, looking for other sources on RK4, the Wikipedia article at
> http://en.wikipedia.org/wiki/RK4 mentions something that was left out 
> here...
> that the weighted average of the forces should be multiplied by the 
> timestep.
>
> That would change the above line to
> (forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])*dt/6.0f;
>
> This would account for the strange scaling effect I saw.  When I get home
> tonight, I'll try it out.

No that is wrong, I assume your "step" function already multiplies the 
forces by dt when it does the Euler integration?  If not then it should, you 
certainly don't need dt to calculate the weighted force.  The weighted force 
is the result of the RK4 method, and simply says "this is the best guess of 
constant force to be applied over the next dt to get nearest to the true 
value of velocity dt time later".

Your code is missing the bit that says "this is the best guess of constant 
velocity to be applied over the next dt to get nearest to the true value of 
position dt time later".

That is why I said it is best to keep your state information separate from 
the derivative information, it is confusing because velocity appears in 
both, and if you don't have them separate it can be hard to know which one 
to use where.

Have something like:

State {position , velocity }
Derivative {velocity , acceleration }

You can put anything else in your state block, but you must put the 
corresponding derivative in the derivative block.  For example, I could 
have:

State {position , angle , temperature}
Derivative {velocity , angluar velocity , temperature gradient }

And that would work fine with RK4.

Your algorithm for RK4 is then:

Derivative A = EvaluateDerivative( currentState , dt )

tempState    = Step( A , currentState, dt*0.5 )
Derivative B = EvaluateDerivative( tempState , dt*0.5 )

tempState    = Step( B , currentState, dt*0.5 )
Derivative C = EvaluateDerivative( tempState , dt*0.5 )

tempState    = Step( C , currentState, dt )
Derivative D = EvaluateDerivative( tempState , dt )

BestGuessDerivative = (A + 2*(B+C) + D)/6
newState = Step( BestGuessDerivative , currentState , dt )

And then, Step would simply do an Euler integration step, of *all* members 
of your state block, using the corresponding members of the derivative 
block.

Step( derivative , state , dt )
{
  newState = state;
  newState.position += derivative.velocity     * dt;
  newState.velocity += derivative.acceleration * dt;
  newState.angle    += derivative.angularVelocity * dt;
  newState.temp     += derivative.tempGradient * dt;
  ...
  return newState;
}

and EvaluateDerivative will work out the derivative block at that point

EvaluateDerivative( state , dt )
{
  derivative deriv;
  deriv.velocity        = state.velocity;
  deriv.acceleration    = CalculateAccelerations( state );
  deriv.angluarVelocity = CalculateAngularVelocity( state );
  deriv.tempGradient    = 50;  // example
  ...
  return Out;
}

The line that says "deriv.velocity = state.velocity" is the key to linking 
the two integrations of acceleration and velocity.  It makes the output from 
the acceleration integration be used as the input for the velocity 
integration.  But it does this at each step inside RK4, not just at the end.

Post a reply to this message

Thanks, I'll give it another go tomorrow (too late tonight, unfortunately).

-- 
...Ben Chambers
www.pacificwebguy.com

Post a reply to this message

Chambers wrote:
> John VanSickle wrote:
>> The last time I did a physics simulator (springs and weights), I found 
>> that calculating the change in acceleration, and figuring that into 
>> the new location and velocity data, greatly reduced the error that 
>> accumulates.
> 
> That's why I'm calculating forces :)

Well, yeah.  But whereas force goes to the second derivative of location 
over time, it's also helpful to calculate the third derivative over 
time, and figure that into the updates for both position and velocity.

The differential calculus behind it isn't that hard, and provides a 
better solution than dampening to the wild oscillations produced by sims 
that only go to the second derivative over time.

Regards,
John

Post a reply to this message

OK, so I did all of that, and I *still* have the problem of the forces 
in RK4 appearing to be a different scale from the forces using the 
Newtonian method.

Shouldn't I be able to just swap out the two, choosing Newtonian for 
speed, or RK4 for accuracy, and see relatively equal* results?


*Obviously not exactly equal, else there's no point.  And also 
obviously, the errors inherent in Newtonian will compound over time, 
meaning the two sets of results will diverge.

-- 
...Ben Chambers
www.pacificwebguy.com

Post a reply to this message