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

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

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

Yes, if your time-step is small enough for Newton to work roughly correctly. 
Do you still get the same difference if you reduce the time step by a factor 
of 100?  If you do, then it's a coding error, if you don't then it's simply 
the fact the Newton is not doing very well at your previous time step 
amount.

Did you try a simple system first, eg constant force or very soft spring 
from one point?  With a constant force you should be able to calculate the 
exact velocity and position by hand and compare with what your algorithm 
gives out.

Post a reply to this message

scott wrote:
>> 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?
> 
> Yes, if your time-step is small enough for Newton to work roughly 
> correctly. Do you still get the same difference if you reduce the time 
> step by a factor of 100?  If you do, then it's a coding error, if you 
> don't then it's simply the fact the Newton is not doing very well at 
> your previous time step amount.
> 
> Did you try a simple system first, eg constant force or very soft spring 
> from one point?  With a constant force you should be able to calculate 
> the exact velocity and position by hand and compare with what your 
> algorithm gives out.

OK, I "fixed" the timestep at 1/200th of a second, and removed all 
particle collisions - now the only force is that of gravity.

Using the two integrators, the particles follow roughly the same path, 
but using RK4 they appear to have an increased velocity - they follow 
the same path much more quickly.

Which means I've got an error somewhere in my code...

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

Post a reply to this message

> OK, I "fixed" the timestep at 1/200th of a second, and removed all 
> particle collisions - now the only force is that of gravity.
>
> Using the two integrators, the particles follow roughly the same path, but 
> using RK4 they appear to have an increased velocity - they follow the same 
> path much more quickly.
>
> Which means I've got an error somewhere in my code...

It should be very easy to debug, simply print out (to a file?) the "state" 
and "derivative of state" information at each step of RK4 for one particle 
(that starts at zero velocity would make things simpler).  eg Just print out 
position and velocity in the gravity direciton.

Attached is a spreadsheet I made that matches up with the code I posted 
before, it shows you the values at each step for a constant acceleration, 
and compares the results with the Euler method and the exact solution. 
Hopefully you can spot the mistake in your code...

Post a reply to this message

Attachments:
Download 'rk4.xls.dat' (44 KB)

scott wrote:
> It should be very easy to debug, simply print out (to a file?) the 
> "state" and "derivative of state" information at each step of RK4 for 
> one particle (that starts at zero velocity would make things simpler).  
> eg Just print out position and velocity in the gravity direciton.

Thanks, I got it sorted out.  I still don't know what the error was, but 
switching from preallocating the arrays to allocating new ones every 
step & throwing them away when done solved the problem.

So there must have been some carryover somewhere that I missed.

Anyway, I'll post a binary soon.  I've been going through the forums for 
XNA studio, and apparently there's no easy way to distribute games for 
windows yet (the whole thing is a work in progress, and they're about to 
release yet another version of it).

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

Post a reply to this message

Chambers wrote:

> Thanks, I got it sorted out.  I still don't know what the error was, but 
> switching from preallocating the arrays to allocating new ones every 
> step & throwing them away when done solved the problem.
> 
> So there must have been some carryover somewhere that I missed.

Ah, the evils of mutable state...


...er, wait, ignore me. [Oh, you already are...]

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

Post a reply to this message

Orchid XP v8 wrote:
> Chambers wrote:
> 
>> Thanks, I got it sorted out.  I still don't know what the error was, 
>> but switching from preallocating the arrays to allocating new ones 
>> every step & throwing them away when done solved the problem.
>>
>> So there must have been some carryover somewhere that I missed.
> 
> Ah, the evils of mutable state...
> ....er, wait, ignore me. [Oh, you already are...]

You know, it's not all that hard to accidentally carry over incorrect 
state in a language without mutable state, too. :-)

-- 
   Darren New / San Diego, CA, USA (PST)
     "That's pretty. Where's that?"
          "It's the Age of Channelwood."
     "We should go there on vacation some time."

Post a reply to this message

Orchid XP v8 wrote:
> Chambers wrote:
> 
>> Thanks, I got it sorted out.  I still don't know what the error was, 
>> but switching from preallocating the arrays to allocating new ones 
>> every step & throwing them away when done solved the problem.
>>
>> So there must have been some carryover somewhere that I missed.
> 
> Ah, the evils of mutable state...
> 
> 
> ...er, wait, ignore me. [Oh, you already are...]
> 

I'm not ignoring you.  It was from reading all your posts on Haskell 
that I got interested in immutable values to begin with.  While I'm 
still not using them very much, at least now they're in the back of my 
mind :)

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

Post a reply to this message

> Thanks, I got it sorted out.  

Cool :-)

Post a reply to this message

Darren New wrote:

> You know, it's not all that hard to accidentally carry over incorrect 
> state in a language without mutable state, too. :-)

Heh. Do you have *any idea* how long it took me to finish debugging this 
JavaScript snippet?

function Bindings_Lookup(key)
{
   if (this.Key  == key) return this.Value;
   this.Next.Lookup(key);
}

It took an absurd amount of time to figure this out - especially given 
the strange way Firefox responded to it...

Sure, they'll never make a programming language where mistakes are 
impossible. But mutable state just offers *so many* ways to shoot 
yourself in the foot. :-/

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

Post a reply to this message

Invisible wrote:
> Sure, they'll never make a programming language where mistakes are 
> impossible. But mutable state just offers *so many* ways to shoot 
> yourself in the foot. :-/

I don't understand what that has to do with mutable state, since you're 
not mutating anything there. It looks like you have the bug of too much 
recursion and the bug of not checking for the end of the list, assuming 
I'm reading it right.

And yes, firefox seems to crash my machine at least once a day, which is 
rather annoying. I've taking to running anything on a remote machine 
that takes more than 5 minutes as a nohup job, just in case my machine 
locks up.

-- 
   Darren New / San Diego, CA, USA (PST)
     "That's pretty. Where's that?"
          "It's the Age of Channelwood."
     "We should go there on vacation some time."

Post a reply to this message