>> I wonder how fast it would be in C++ :-)
>
> Go ahead, knock yourself out. ;-)
>
>> Are the three attractors just applying a 1/r^2 force on each particle?
>
> Nope. Currently just 1/r for the "magnets", plus an r force for "gravity".
>
>> What resolution are you running this at?  "a few thousand" sounds like 
>> 64x64 or something.
>
> It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)

OK so I quickly modified my existing RK4 code to do integration steps on 
1728*972 particles in 2D with 4x 1/r attractors.  So it takes about 2.5 
seconds per integration step, does this compare with your results?

Post a reply to this message

scott wrote:
>>> What resolution are you running this at?  "a few thousand" sounds 
>>> like 64x64 or something.
>>
>> It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)
> 
> Ah ok, so "a few" means "over 1600" :-)

1,679,616 pixels, baby! :-D

(Or what you DSLR people would call "1.7 megapixels". :-P )

>> I *said* YouTube doesn't do it justice. ;-)
> 
> I can well imagine.

When I get home I'll put up some still frames in full resolution.

> Also I can imagine, this is the sort of stuff I like playing about with 
> too. I wonder if the plane seat will have a power socket for my laptop 
> ... hmmmm.

Muhuhuhu! >:-D

Post a reply to this message

scott wrote:

> OK so I quickly modified my existing RK4 code to do integration steps on 
> 1728*972 particles in 2D with 4x 1/r attractors.  So it takes about 2.5 
> seconds per integration step, does this compare with your results?

It's roughly comparable, certainly.

(What hardware are you running it on?)

Post a reply to this message

>> OK so I quickly modified my existing RK4 code to do integration steps on 
>> 1728*972 particles in 2D with 4x 1/r attractors.  So it takes about 2.5 
>> seconds per integration step, does this compare with your results?
>
> It's roughly comparable, certainly.
>
> (What hardware are you running it on?)

That was using 1 core of an Intel E8500 @ 3.16 GHz.

Of course someone has already done particle system simulation on the GPU 
before:

http://2ld.de/gdc2004/

Running option 2 of that demo on my GPU gives 15fps for a million particles, 
so for 1728x972 I make that about 9fps.  If you restrict it to 2D rather 
than 3D it should go a little faster.

Post a reply to this message

scott wrote:
>>> OK so I quickly modified my existing RK4 code to do integration steps 
>>> on 1728*972 particles in 2D with 4x 1/r attractors.  So it takes 
>>> about 2.5 seconds per integration step, does this compare with your 
>>> results?
>>
>> It's roughly comparable, certainly.
>>
>> (What hardware are you running it on?)
> 
> That was using 1 core of an Intel E8500 @ 3.16 GHz.

OK. My hardware is considerably less powerful.

Is that RK4 with a single step, or with adaptive step subdivision?

> Of course someone has already done particle system simulation on the GPU 
> before:
> 
> http://2ld.de/gdc2004/
> 
> Running option 2 of that demo on my GPU gives 15fps for a million 
> particles, so for 1728x972 I make that about 9fps.  If you restrict it 
> to 2D rather than 3D it should go a little faster.

Heh, nice... Coming to a gaming station near you, eh?

Post a reply to this message

>> That was using 1 core of an Intel E8500 @ 3.16 GHz.
>
> OK. My hardware is considerably less powerful.
>
> Is that RK4 with a single step, or with adaptive step subdivision?

Just one complete integration step of plain RK4.  I didn't do any checking 
to see if sub-steps are needed or not.

Post a reply to this message

scott wrote:
>>> That was using 1 core of an Intel E8500 @ 3.16 GHz.
>>
>> OK. My hardware is considerably less powerful.
>>
>> Is that RK4 with a single step, or with adaptive step subdivision?
> 
> Just one complete integration step of plain RK4.  I didn't do any 
> checking to see if sub-steps are needed or not.

That's probably why my version is slower then.

It computes one step forward. It then recomputes that step as two 
half-steps. Finally, it computes the difference between the two answers. 
If the difference is too large, the first result is thrown away, and 
further subdivision is done to check the second result is OK. (There is 
a maxmimum number of subdivisions though - or rather, a minimum step size.)

I can try turning that off to see what speed I get - but I'm still 
running on different hardware to you, so it's probably not very meaningful.

Post a reply to this message

> It computes one step forward. It then recomputes that step as two 
> half-steps.

That sounds like a waste of time.  Within the RK4 algorithm you already have 
4 estimates of slope, just check the variance of these to determine whether 
you need to use a smaller time step or not.

Post a reply to this message

scott wrote:
>> It computes one step forward. It then recomputes that step as two 
>> half-steps.
> 
> That sounds like a waste of time.  Within the RK4 algorithm you already 
> have 4 estimates of slope, just check the variance of these to determine 
> whether you need to use a smaller time step or not.

Well, the stability depends on how sharply the force acting on the 
particle changes as a function of position...

I have to say, I don't really comprehend how RK4 works. I copied the 
formulas, and it seems to produce the right answers, but I have no idea 
*why* it works or what it's doing.

Still, I could try what you suggest and see what the accuracy / 
performance is like. (I already had to decrease the error tollerance 
once to get crisp images...)

Post a reply to this message

> I have to say, I don't really comprehend how RK4 works. I copied the 
> formulas, and it seems to produce the right answers, but I have no idea 
> *why* it works or what it's doing.

It's easier to understand if you consider a single function to integrate 
first, rather than doing acceleration/velocity and velocity/position in one 
go.  There is a clever trick that lets you do both of these together without 
much extra work, but it's not too obvious and usually confuses things.

Let's suppose you have a function that calculates the instantaneous velocity 
of an object based on its position (x) and a time (t):

v = F(x,t)

Now using Euler integration you can easily track the position of the object 
like this:

> t = 0
> x = 0
> repeat
>  average_velocity = F(x,t)
>  x = x + average_velocity * dt
>  t = t + dt

Where dt is some time step.

As you are aware, the problem is that with large time-steps big errors 
appear because you are assuming the average velocity during the time step is 
equal to the instantaneous velocity at the beginning of the time step - this 
is only approximately correct!

All RK4 does is to find a much better estimate of the average_velocity 
during the time step, then you can use that to advance x with much higher 
accuracy:

> t = 0
> x = 0
> repeat
>  average_velocity = RK4_Algorithm( x , t , dt )
>  x = x + average_velocity * dt
>  t = t + dt

The RK4 algorithm looks like this:

> Function RK4_Algorithm( x , t , dt )
> AVE1 = F(x,t)

First step there is to compute an Average Velocity Estimate by simply using 
F(x,t) like we did in Euler.

> AVE2 = F( x + AVE1 * dt/2 , t + dt/2)

This is the first clever bit, we are going to calculate another estimate for 
the velocity *half* a time step forward.  We use the first estimate of 
velocity (AVE1) to work out what the position will be half a time step 
later, then we then use this new position in our velocity finding function 
to find the velocity there.  This gives us our estimate of velocity at the 
middle of the time step, call it AVE2.

> AVE3 = F( x + AVE2 * dt/2 , t + dt/2)

This does the same as before, but instead of using AVE1 as the assumption 
for average velocity, we use AVE2.

> AVE4 = F( x + AVE3        , t + dt  )

This is the 4th estimate, it uses the 3rd one to do a *full* time step 
forward, and calculate the velocity there (at the *end* of the time step).

Now we have four estimates for the velocity, one at the beginning of the 
time step (standard Euler), two in the middle, and one at the end of the 
time step.

The RK4 algorithm then tells us to take a weighted average of these 
velocities:

> AVE = (AVE1 + 2*AVE2 + 2*AVE3 + AVE4) / 6

That is it!  That is the estimate of average velocity over the entire time 
step.  You will have noticed we called F four times, but doing this is MORE 
accurate than simply doing four dt/4 Euler time steps.

> return AVE

You should be able to simply compare the values of AVE1,AVE2,AVE3,AVE4 to 
determine if you need to go back and do a smaller time step.

Now to make things more general, you can replace the position "x" by the 
generic "state" of the system, and the velocities by the "derivative of 
state" of the system.  This way, the functions F and RK4_Algorithm will 
return objects of type "derivative of state".

If you are dealing with forces acting on particles, then the "state" of each 
particle is its position and momentum (4 variables in 2D), and the 
"derivative of state" is its velocity and force (another 4 variables).

Now, your function F must return the velocity and force.  The force it 
calculates based on your magnet formula or whatever, and the velocity it can 
easily calculate by dividing the supplied momentum value by the mass.  That 
last part is the key that links the two differential equations together, and 
why it is sometimes complicated how RK4 can do both at the same time for 
acceleration and velocity.

Post a reply to this message