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