> 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