> v[n+1] = v[n] + a[n] * ts
> p[n+1] = p[n] + v[n] * ts + 0.5 * a[n] * ts^2
>
> Factoring in the equation for velocity into the equation for position you get:
>
> p[n+1] = p[n] + 0.5 * ( v[n] + v[n+1] ) * ts

> This should make it easier to apply Runge-Kutta, but not by much. I'm still
not
> too clear on how myself, but I'm working on it...

Right I've come up with a solution using the above equations. Whether it's the
correct solution remains to be seen. I'm not entirely happy with the integration
from velocity to position as I have this feeling it's not a Runge-Kutta
integration, but Euler instead. That said, I ran quick test of my own, and after
800 steps, the error was 7e-11% for the velocity and 7e-4% for the position
(compared to Euler: 0.4% for velocity and 0.5% for position). I just hope the
way I've applied it to your system actually works. Here goes...

The differential function in your system is the acceleration, which is
calculated from the forces acting on the particles, which in turn is calculated
from the position of all the particles. We can therefore say that the
differential function is a function of position: A(p)

The velocity after one time-step is calculated from the current velocity plus
the integral of the acceleration for one time-step:

v[n+1] = v[n] + V( p[n] )
where
V(p) is the integral of A(p) for one time-step at the current position

Using Runge-Kutta to integrate A(x), the intermediate values k0 to k3 are:

k0 = ts * A( p[n] )
k1 = ts * A( p[n] + k0 / 2 )
k2 = ts * A( p[n] + k1 / 2 )
k3 = ts * A( p[n] + k2 )

And the equation of velocity becomes:

v[n+1] = v[n] + ( k0 + 2 * ( k1 + k2 ) + k3 ) / 6

Note that in the calculation of k0 to k3, you need to recalculate the
acceleration using the new estimates of all the positions of all the particles.
A( p[n] ) is just the acceleration calculated from current position of all the
particles. To calculate k1 you must calculate the accleration of the all
particles in the position: p[n] + k0 / 2 i.e. the current position plus the
first estimate. And the same for k2 and k3. Basically, you have to recalculate
the position of all the particles using the current set of k values in order to
calculate the acceleration of all the particles so that you can calculate the
next k value for each particle.

Once you have all the values of k, you just plug them into the equation of
velocity. The equation of position for now is simply:

p[n+1] = p[n] + 0.5 * ( v[n] + v[n+1] ) * ts

I can't find a better way to calculate it. This is equivalent to:

p[n+1] = p[n] + v[n] * ts + 0.5 * ts * ( k0 + 2 * ( k1 + k2 ) + k3 ) / 6

but the former is quicker and easier to calculate. I'd like to find a way of
applying Runge-Kutta to the integration of velocity into position, but things
get too complicated when I try that, and I don't know if this is necessary or
not. For now, this is still very accurate, and far more so that using Euler.
Also, due to the increased accuracy of Runge-Kutta, you'll be able to reduce the
step size.

Hope that all helps. I really don't feel like replying to myself any more...

Ryan

Post a reply to this message