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) :)

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

Post a reply to this message

>  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);

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.

I find it easier if you keep two distinct data structures, one for the state 
information (position, velocity), and one for the derivatives 
(momentum,force or velocity,acceleration it doesn't matter for point 
masses).

Then it is clear that "collide" must return a "derivative of state" data 
structure which must be then used in the following step function.  At the 
end, you do the weighted average of the "derivative of state" data structure 
as a whole.

This allows you to add more stuff easily in future (eg orientation and 
angular momentum) without having to worry about the actual RK4 core code.

Post a reply to this message

> 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