> OTOH, almost all physical systems are of this type.

Indeed, and the exact solution of that is a damped sine wave, which is why 
so much stuff vibrates and resonates in the real world.

>> 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.
>
> Right. I don't really understand why it has to be these exact points...
>
>> The RK4 algorithm then tells us to take a weighted average of these 
>> velocities:
>
> ...or why this specific weighting is the correct one. But hey, I guess 
> somebody had a reason for that.

It matches the Taylor Series approximation of the function:

http://grb.physics.unlv.edu/~zbb/files/RungeKuttaProof.pdf

> Clearly I'm going to have to meditate on exactly why that's the case. I 
> would have thought that taking 4 seperate steps would be more accurate 
> than computing 4 steps and then blurring them all into one step, but hey.

RK4 is not just taking 4 blind steps though, it's using the result of the 
previous ones to get the next ones, which is how it gets to be more accurate 
than a simple linear estimate.

>> 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.
>
> Hmm... How is momentum and velocity different?
>
> Oh, wait - the mass could be non-unital. I forgot about that...

Yes that's rather key :-) I also find that keeping momentum as a state 
variable rather than velocity (which would also work) avoids some of the 
confusion with the differential of position in the algorithm.

> Thinking about Euler integration, you end up doing something like
>
>   x' = x' + dt x''
>   x  = x  + dt x'
>
> In other words, you compute the new velocity from the old one, and then 
> assume that velocity is constant while you compute the new position from 
> the old one. But if you assume that acceleration is constant, you can 
> actually *exactly solve* for position. I can't figure out the equation off 
> the top of my head,

x = x + dt x' + 1/2 x'' dt^2

That looks a bit like a Taylor Series expansion, doesn't it :-)  You can 
extend that to:

x = x + dt x' + 1/2 x'' dt^2 + 1/6 x''' dt^3

and also write it for x' instead of x:

x' = x' + dt x'' + 1/2 x''' dt^2 + 1/6 x'''' dt^3

And that is exactly what your RK4 algorithm is estimating :-)

Post a reply to this message