OK, so over here in the UK the clocks just changed over from British Summer
Time to Grenwich Mean Time. Hence, although it was 10PM just a while ago, it
is now 9PM [again]. All of which means that I will be awake for another
hour. And what better reason to ask a random question? (OK, maybe having a
*point* would help, but hey...)

As some of you may or may not have noticed, I am currently at work on
another motion simulation. It began ages and ages ago with the "chaos
pendulumn" simulation, but current work is on a grid of particles held
together by elastic connections and allowed to deform under gravity.

Anyway, someone pointed out that my method is "Euler's method", the simplest
but least accurate method. (How exactly someone can tell just from the
pictures what algorithm I'm using I don't know... lol) So I was wondering...
what other methods are there? Have I missed something??

All I know is this: POV-Ray wants to know the positions of all the objects.
I know the initial positions - I *chose* them ;-) But an object's position
can change over time - something we usually refer to as "velocity".
Additionally, an object's velocity can change over time also - but only if a
"force" is applied to it. This force produces an "acceleration" which is
equal to the quotient of the force and the object's mass...

From the beginning, velocity is the rate at which an object's position
changes over time. Taking the movement over 1 dimension only, if an object's
velocity is 3 meters per second, that would mean that in 1 second it will be
3 meters away from where ever it is now. SO... between each frame, I can
take the object's previous position, and add on the velocity multiplied by
the number of seconds that have elapsed since the last frame. (Or, more
likely, the fraction of a second that has elapsed!)

In a similar way, an acceleration of 3 meters per second per second would
mean that in 1 seconds' time, the object will be moving at 3 meters per
second *faster* than however fast it is moving now...

In short, (working in 1 dimension) once I calculate the force acting on an
object (which depends on what I'm simulating), I can do this:
#declare Acceleration = Force / Mass;
#declare Velocity = Velocity + Acceleration*TimeStep;
#declare Position = Position + Velocity*TimeStep;
Now just add some file I/O to load/save the particle positions and
velocities between frames. (As far as I can see, this is the only
information you need from the previous frame.)

Even better, if I'm not mistaken, if I write each of the forces acting on
the object as a vector, then the resultant force is equal to the sum of the
vectors. (Correct?) Also, if I let Acceleration, Velocity and Position be
vectors, the above formulas still work. (Unless I missed something...)

I have found that for my current simulation, if I set certain parameters
certain ways, the simulation seems to go into "infinite acceleration" mode.
(I've actually had POV-Ray CRASH a few times due to number overflow
problems - I say chaps, should it really do that??) Anyway, following some
experiments I did with Excell in my lunch break at work, I tried changing
the simulation so that it does *several* updates each frame - in other
words, the TimeStep is now smaller tham 1/25 seconds.

I tried one set of parameters, and it went totally nutz. I changed the
simulation so that it does 10 updates per frame - 1/250 seconds interval -
and suddenly it's totally stable and works perfectly. Well, *almost*
perfectly... there are still little instabilities in it that seem to take
forever to dissapate. I'm not sure if using a finer grid would solve this,
or whether I need an even smaller TimeStep... or whether I need to use
something other than Euler's method (if, indeed, that is what I'm currently
using).

Basically, I can't really see what I could do any differently here... Any
suggestions folks?

Thanks.
Andrew.

PS. Yes, there *are* dissapative forces. In fact, last time round, I did it
by multiplying each object's velocity by (1 - exp(-8)) (because it's easier
to adjust the exponent). But that means if you change the frame rate, the
dissapation changes. So now it's done by creating a [small] force in the
opposite direction to the object's current velocity. Because it's a force,
it goes through that formula above, which depends on TimeStep. So if I
change the frame rate (or indeed the number of updates per frame) it
shouldn't affect the simulation in that way.

Post a reply to this message

The time-step you use will determine the accuracy of the system you are trying
to simulate. If it is too long, you will be missing some of the 'detail' in the
simulation, this can result in instability in the simulation. On the other hand,
if it is too short you will waste processing time unnecessarily.

The instability can be caused by a long time-step introducing an error in your
calculations, which builds up over succesive time-steps, becoming ever larger.
Ultimately you will have to find a time step that will produce an acceptably low
error (and hence reduce or eliminate instability) over the full length of the
simulation.

You should also note that some simulations produce stable oscillations in the
system which don't ever dissipate. This will be because the system you are
simulating would in real life oscillate in this way. If increasing the time-step
does not remove all oscillation then you can assume they are there to stay.

The rate of change of the forces on the particles will determine the time-step
you need. If the forces change rapidly then your system will require a very
short time-step to correctly simulate the motion of the particle under those
conditions. Otherwise a longer time-step will not properly account for those
changes in force.

For example:
If you were to apply an oscillatory acceleration of the form:
acc=A*cos(2*pi*t)
If you sample the acceleration every 1 second, your simulation will see a
constant acceleration of magnitude A for the entire duration of the simulation
(resulting in particles that start moving and just get faster and faster). This
is obviously not the kind of behaviour that an oscillatory acceleration should
produce, and is a result of the sample rate being too long.

Post a reply to this message

> The time-step you use will determine the accuracy of the system you are
trying
> to simulate. If it is too long, you will be missing some of the 'detail'
in the
> simulation, this can result in instability in the simulation. On the other
hand,
> if it is too short you will waste processing time unnecessarily.

Yup, sure... the classic accuracy/speed tradeoff found in all computer
graphics algorithms. (Strange, that... *wink*)

> You should also note that some simulations produce stable oscillations in
the
> system which don't ever dissipate. This will be because the system you are
> simulating would in real life oscillate in this way. If increasing the
time-step
> does not remove all oscillation then you can assume they are there to
stay.

Yes... particularly my simulation of elastic springs... you would certainly
expect elastic to oscilate to some degree or other.

> The rate of change of the forces on the particles will determine the
time-step
> you need. If the forces change rapidly then your system will require a
very
> short time-step to correctly simulate the motion of the particle under
those
> conditions. Otherwise a longer time-step will not properly account for
those
> changes in force.

Yes, that's clear.

> For example:
> If you were to apply an oscillatory acceleration of the form:
> acc=A*cos(2*pi*t)
> If you sample the acceleration every 1 second, your simulation will see a
> constant acceleration of magnitude A for the entire duration of the
simulation
> (resulting in particles that start moving and just get faster and faster).
This
> is obviously not the kind of behaviour that an oscillatory acceleration
should
> produce, and is a result of the sample rate being too long.

Mmm... that's just the Nyquist limit comming in to play... sadly, for other
types of simulation, it's not to easy to calculate an adiquate time step
value... I am currently working with 1/250 seconds, and it seems to work
well.

Further to other messages I've posted, I tried a much finer grid (with
smaller masses so the total mass is still unchanged). It seems a lot more
stable...

I was really reacting to what Chris Johnson said - "Euler Integration is the
simplest, least accurate method. [...] Most people here use the 4th-order
Runge-Kutta method". I have no idea what that means! lol

Thanks.
Andrew.

Post a reply to this message

> Hmm, having written all that out, I'm not entirely sure how you go about
> applying it to your system. This could be tricky due to the fact you're
> calculating both position and velocity from your acceleration values, which
> ifself is position/velocity dependent. Hmm, I'm going to have to think about
> this for a while...

Right, worked out what was confusing me. There's something not quite right about
your equations. Between each time step you are assuming that the acceleration is
constant. This is fine, and is due to your use of Euler's integration. However
your equations aren't reflecting this. Typically the equations of motion of a
particle with constant accleration are of the form:

v = v0 + a * t
p = p0 + v0 * t + 0.5 * a * t^2

For your system the following substitutions should be made:
t = ts (= time step)
v = v[n+1] (= velocity at next time step)
p = p[n+1] (= position at next time step)
v0 = v[n] (= current velocity)
p0 = p[n] (= current position)
a = a[n] (= acceleration calculated from resultant of forces at current time)

Therefore the equations of motion should be of the form:

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

However your equations were of the form:

v[n+1] = v[n] + a[n] * ts
p[n+1] = p[n] + v[n+1] * ts

which is incorrect. Your equations were assuming that the position would change
according to the future velocity, whereas with a constant acceleration it would
change according to the average of the current and future velocity.

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...

Ryan

Post a reply to this message

> 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