"Orchid XP v2" <voi### [at] devnull> schrieb im Newsbeitrag 
news:43079838$1@news.povray.org...
> There comes a time in every POVer's lifetime when they have to know...
>
> I've Googled, I've visited Mathworld and I've checked out Wikipedia. (That 
> last one actually made sense, impressively.) So I now know what RK4 
> *is*... and roughly how it works... but... Integral calculus isn't my 
> strong point. (I'm much better at differential calculus.) The notation 
> confuses me. Plus, the equations of motion usually aren't first-order.
>
> So, I'm *sure* this is in a VFAQ _somewhere_... anybody have an 
> explanation for how you do physics simulations with this thing?

i'll try to explain it:
i'll take a simple harmonic oscillator as an example.
you start with your differential equations in a vectorised form
so that the timederivative is on the left side and the right side
is a function of all your variables. In our example you start with
F = m*a = -k*x  =>  a = -k/m*x
Acceleration is the second derivative of the distance with respect to

but we need the differential equation in a form that contains only
first order time derivatives. Therefore we just introduce a helper
variable, namely the velocity. Then we can formulate our second order
diff. eq. with two first order diff. eqs.:
instead of

we write
 dx/dt = v
 dv/dt = -k/m*x
written more generally as matrix:
    d /x\   / 0    1\ /x\
   -- | | = |       |.| |
   dt \v/   \-k/m  0/ \v/

now to the runge-kutta part:
The right side of the last equation can be thought of as a function
of the coordinates in phasespace (in our case x,v) which tells us
the first timederivative. In general the time itself 't' could also
influence the rightside (for example if the springconstant of our
oscillator changes with time). Now we have all what we need to
simulate our differential equation since our function gives us the
timederivative of all systemcoordinates for every possible state
of the mechanical (or whatever) system.
We call this function f(x,t) (where x is an n-dimensional vector)
and it returns an n-dimensional vector with the timederivatives.

Now we can formulate the various integration algorithms in an easy
way with our function f(x,t):
Let the stepsize in time be 'h'.
Then an ordinary Euler integration step can be seen as:

 x(t+h)=x(t) + h*f(x,t)

we just take the old coordinates and add our stepsize times the
timederivative which we get from our function f
In this way the Runge-Kutta 4th order step looks like this:

k1=h*f(x(t)     ,t    )
k2=h*f(x(t)+k1/2,t+h/2)
k3=h*f(x(t)+k2/2,t+h/2)
k4=h*f(x(t)+k3  ,t+h  )
x(t+h)=x(t) + 1/6 * (k1+2*k2+2*k3+k4)

the k1,...,k4 are multiple guesses for the time-derivative
which build upon each other. Then they are combined in
a weighted average to give the final estimate for the slope.

so the essential idea is to formulate your physical problem
as function f(x,t), then it's very easy to implement.
If the programming language doesn't support vector-valued
returnvalues you have to make a function f for every vector-
component, so in the example we would have fx and fv and then

kx1=h*fx(x(t)      ,v(t)      ,t    )
kv1=h*fv(x(t)      ,v(t)      ,t    )
kx2=h*fx(x(t)+kx1/2,v(t)+kv1/2,t+h/2)
kv2=h*fv(x(t)+kx1/2,v(t)+kv1/2,t+h/2)
kx3=h*fx(x(t)+kx2/2,v(t)+kv2/2,t+h/2)
kv3=h*fv(x(t)+kx2/2,v(t)+kv2/2,t+h/2)
kx4=h*fx(x(t)+kx3  ,v(t)+kv3  ,t+h  )
kv4=h*fv(x(t)+kx3  ,v(t)+kv3  ,t+h  )

x(t+h)=x(t) + 1/6 * (kx1+2*kx2+2*kx3+kx4)
v(t+h)=v(t) + 1/6 * (vx1+2*vx2+2*vx3+vx4)

if you have a 3d-problem, you'll have
6 formulas like that for px,py,pz,vx,vy,vz
it's just more complicated to write, but not really anything new,
you get the idea...

hope that get's you started :)

Thies

Post a reply to this message