> Well, I don't know about that, but physics is mostly about forces, and 
> forces generate acceleration, and force is usually a function of position 
> (or perhaps velocity). Hence differential equations seem to be kind of 
> ubiquitous in physics.
>
> Interestingly, a trivial differential equation has exponential growth as 
> its solution, but you never see this in mechanical systems.

Diesel engines would destroy themselves without a control system to keep the 
speed down, that's the best I can come up with now of a mechanical system 
with exponential growth.  I think the key point is that most things have a 
force that acts in the opposite direction of the movement, which usually 
keeps things stable and also produces the oscillations.

> I thought a Taylor series is just a way of constructing a polynomial 
> approximation to a function that you don't have any better way to compute.

Indeed, so if we can estimate the derivatives by some numerical method, we 
should be able to estimate how the actual function goes forwards, at least 
for a little while.

> And taking 4 linear steps uses the result of the previous steps for the 
> next one.

But it doesn't use the results to go back and work on the *same* part of the 
function to improve the accuracy, it just blindly moves onto the next time 
step ignoring whatever errors might have been present.

> I don't get why taking those steps, throwing the result away and just 
> keeping the intermediate results, and then mixing that all together gives 
> a better result. (Altough it clearly *does*.)

Another way to look at it I guess is that RK4 is effectively estimating 
"deeper" derivatives, like x'' and x''' to better estimate how the function 
goes.

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

It's also crucial when you come do doing rotations as well as linear 
movements.  For linear movements you can simply add on the accerlation times 
dt to the velocity, then the same for position.  But you cannot do this for 
angular movements, you need to add on the torque to angluar momentum, use 
that to calculate a thing called spin, then add that onto the angular 
orientation.

> I don't follow. (I.e., I don't see how the final formula is remotely 
> related to RK4.)

Try writing a function to estimate x' x'' and x''' of a function at a 
certain time and position.  Then use it to predict what the value of x will 
be at the next time step based on the Taylor Series expansion.  Then, 
assuming you are estimating the derivatives in the "correct" way, you can 
simplify your solution and you'll end up with the RK4 algorithm :-)  I have 
no idea what mathematical jiggery pokery you'll need to do that, but that is 
essentially what RK4 is doing.

Euler integration estimates the next value by assuming:

x = x + dt x'

ie it will perfectly match a function like x = 4 + 5t

RK4 estimates the next value by assuming:

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

ie it will perfectly match a function like x = 4 + 5t - 3t^2 + 45t^4

Post a reply to this message

> I think the key point is that most 
> things have a force that acts in the opposite direction of the movement, 
> which usually keeps things stable and also produces the oscillations.

Yeah.

Of course, the fun part happens when objects hit each other. Then you 
have a sudden and drastic change in the forces applied. Good luck 
integrating that...

>> I thought a Taylor series is just a way of constructing a polynomial 
>> approximation to a function that you don't have any better way to 
>> compute.
> 
> Indeed, so if we can estimate the derivatives by some numerical method, 
> we should be able to estimate how the actual function goes forwards, at 
> least for a little while.

So... well, we don't need to "estimate" the acceleration at all. We can 
compute it directly. (Assuming we have an exact position...) But you're 
saying you compute the acceleration several different times in order to 
estimate its integral or something?

> Another way to look at it I guess is that RK4 is effectively estimating 
> "deeper" derivatives, like x'' and x''' to better estimate how the 
> function goes.

...so you compute acceleration several times to try to estimate the 
derivative of acceleration? (I think the term is "jerk".) That makes 
slightly more sense.

> It's also crucial when you come do doing rotations as well as linear 
> movements.  For linear movements you can simply add on the accerlation 
> times dt to the velocity, then the same for position.  But you cannot do 
> this for angular movements, you need to add on the torque to angluar 
> momentum, use that to calculate a thing called spin, then add that onto 
> the angular orientation.

o_O

Let us pray I never need to handle anything other than point masses...

>> I don't follow. (I.e., I don't see how the final formula is remotely 
>> related to RK4.)
> 
> Try writing a function to estimate x' x'' and x''' of a function at a 
> certain time and position.

Well, if you have a way to compute f(x), you could try using

   f'(x) ~= f(x + dx) - f(x)

By analogy you have

   f''(x) ~= f'(x + dx) - f'(x)

and substituting the previous you eventually arrive at

   f''(x) ~= f(x + 2 dx) - 2 f(x + dx) + f(x)

By a similar process you have

   f'''(x) ~= f(x + 3 dx) - 3 f(x + 2 dx) + 3 f(x + dx) - f(x)

Of course, none of this does you any good if you can't compute f in the 
first place.

> Then use it to predict what the value of x 
> will be at the next time step based on the Taylor Series expansion.

Ouch. That sounds like fun...

> Then, assuming you are estimating the derivatives in the "correct" way, 
> you can simplify your solution and you'll end up with the RK4 algorithm 
> :-)  I have no idea what mathematical jiggery pokery you'll need to do 
> that, but that is essentially what RK4 is doing.

Mmm, right.

Of course, now I'm trying to comprehend how this *explicit* method 
differs from the more accurate *implicit* method I keep hearing about...

Post a reply to this message

> Of course, the fun part happens when objects hit each other. Then you have 
> a sudden and drastic change in the forces applied. Good luck integrating 
> that...

Usually in real-time simulations you can just apply an impulse to the body 
and bypass the integrator altogether.  Your collision response algorithm 
will calculate exactly what impulse to apply to each body so that they 
"bounce" realistically, you then directly modify the momentum values, and in 
the next step RK4 doesn't know what happened :-)

A more accurate way is to actually simulate the contact deformation with 
forces, but usually due to the big forces involved you need very small time 
steps which is beyond real time simulation.

> So... well, we don't need to "estimate" the acceleration at all. We can 
> compute it directly. (Assuming we have an exact position...) But you're 
> saying you compute the acceleration several different times in order to 
> estimate its integral or something?

By computing x' several times (which you can exactly) you are in effect 
estimating x'' and x''', which gives you a better estimate for x than just 
using x' alone.  Remember that here "x" can be position or velocity as far 
as RK4 is concerned.

> ...so you compute acceleration several times to try to estimate the 
> derivative of acceleration? (I think the term is "jerk".) That makes 
> slightly more sense.

Yes.

> Let us pray I never need to handle anything other than point masses...

If you have set up your RK4 code nice and object orientatedly, it is really 
very simple to add in.

> Well, if you have a way to compute f(x), you could try using

You don't though, you are trying to find f(x+dt) given that you can easily 
calculate f'(x).

f''(x)  ~= (f'(x+dt) - f'(x)) / dt
f'''(x) ~= (f''(x+dt) - f''(x)) / dt
f''''(x) ~= (f'''(x+dt) - f'''(x)) / dt

If you substitute those expressions into the Taylor series expansion and 
write it in terms of just f'(x) then you will probably end up with something 
similar to RK4.  Maybe :-)

f(x) = f(x) + f'(x)*dt + 1/2*f''(x)*dt^2 + ...

> Of course, now I'm trying to comprehend how this *explicit* method differs 
> from the more accurate *implicit* method I keep hearing about...

You can usually rewrite your equations of motion for your whole simulation 
as one big matrix equation, something like Ax + Bx' + Cx'' = D where A-D are 
matrices and x is a vector of all your state variables.  C would typically 
contain masses and inertias, A stiffnesses, B damping and D external forces. 
An implicit integrator works on this entire formula in one stab, solving 
usually by matrix inversion, iterative methods or whatever.  For big 
interacting systems it is much better as the constraints (eg two masses 
connected by a spring) are forced throughout the calculation, and not just 
adjusted at the end of each time step.  This is how mechanical finite 
element simulation software works, essentially solving a massive matrix 
equation for each time step.

Post a reply to this message

Invisible wrote:

> PS. I still need to fix the rendering program. I wrote it while at work, 
> and it worked OK, but was rather inflexible and brittle. So when I got 
> home, I wrote it again. But for some reason, the version I wrote at home 
> is about 10x slower - and I can't figure out why! o_O

I just spent about 2 hours peering over the source code for the two 
programs, back to back, side by side. Considering that they were written 
on different ways on different computers, it's remarkable how many of 
the files are identical character-by-character. (Most of them differ by 
variable names, function names, order of arguments, etc.)

I even tried replacing source files from the defective version with 
files from the optimal version, but to no avail. The program is big, but 
not *that* damned big! There's only so many places to look for bugs...

I tried everything I could think of to fix the problem. *Finally*, I 
discovered that replacing one block of code with an apparently 
equivilent block causes a 10x speedup. Whiskey Tango Foxtrot?

At least, I know the bug is inside one tiny function. But where?

Well, I just nailed it. It was in the adaptive subsampling code. When 
computing the error in the current step, I wrote

   error = mag (position1 - position3) + mag (velocity1 - position3)

Looks fine, right?

...until you put it on consecutive lines:

   errorP = mag (position1 - position3)
   errorV = mag (velocity1 - position3)

PWN3D!

Change the second line to

   errorV = mag (velocity1 - velocity3)

and suddenly the program speeds up massively.

D'OH!! >_<

Post a reply to this message

And now for some speed tuning.

A simple way to measure the program's speed is to count how many frames 
per minute it can render. So, set it to 320x480 so it's moderately fast 
[on this ancient single-core 32-bit 1.5 GHz CPU] and count.

After fixing the stupid bug: ~24 frames/minute.

Changed array indexing from tuples to integers: ~24 frames/minute.

Merged the four arrays into one interleaved array: ~26 frames/minute.

[I'm thinking this might improve cache behaviour slightly... but 
apparently not much. Or just the cache isn't the bottleneck.]

Removed array bounds checks: ~27 frames/minute.

Removed adaptive time step: ~80 frames/minute.

[Well there's a surprise. Pitty the frames look awful...]

Added adaptive time step within RK4 as per Scott's suggestion: ~1 
frame/minute.

[Ah, but now the "error" is being calculated in a totally different way, 
so the error tolerance required to get a good image may have changed. 
What happens if I use a less strict tolerance?]

Turned down error bound: ~68 frames/minute.

[I made the maximum permissible error 10,000x larger, and the images 
still come out looking flawless. But they render vastly faster...!]

Post a reply to this message

scott wrote:
> I think the key point is that most 
> things have a force that acts in the opposite direction of the movement, 

I've always found this to be quite fascinating. Why should that be?

Why is back EMF always in the direction that reduces the EMF? Why does 
friction work in a direction that reduces the effects of friction?

Almost everything one can think of works this way. I only remember 
encountering two or three effects in the last few decades (other than 
economics and some odd quantum stuff, IIRC) that don't damp themselves. 
Unless you get something like thermal runaway, where you have two separate 
effects reinforcing each other that would individually damp each other out.

Funkily enough, if you look at wikipedia under "Negative_mass", there's some 
stuff in there about what happens if mass is negative, with fun observations 
like "negative masses would produce an attractive force on one another, but 
would be repelled because of their negative inertial masses". I.e., the two 
minus signs on the mass cancel out, so two negative masses are attracted to 
each other, *but* they'll move apart, because negative mass moves towards 
the direction you push it from. :-)

</ramble>

-- 
   Darren New, San Diego CA, USA (PST)
   There's no CD like OCD, there's no CD I knoooow!

Post a reply to this message

Invisible wrote:
> Change the second line to
>   errorV = mag (velocity1 - velocity3)
> and suddenly the program speeds up massively.

Shouldn't velocity and position be incompatible types?

-- 
   Darren New, San Diego CA, USA (PST)
   There's no CD like OCD, there's no CD I knoooow!

Post a reply to this message

Darren New wrote:
> Invisible wrote:
>> Change the second line to
>>   errorV = mag (velocity1 - velocity3)
>> and suddenly the program speeds up massively.
> 
> Shouldn't velocity and position be incompatible types?

Well, if you wanted to go that far you could. I just had them both as 2D 
vectors...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message

Orchid XP v8 wrote:
> Well, if you wanted to go that far you could. 

Well, now you know why people do, when it's important to be correct. :-)

Isn't it pretty trivial to do that in Haskell, tho?

-- 
   Darren New, San Diego CA, USA (PST)
   There's no CD like OCD, there's no CD I knoooow!

Post a reply to this message

Invisible wrote:

> And now for some speed tuning.
> 
> A simple way to measure the program's speed is to count how many frames
> per minute it can render. So, set it to 320x480 so it's moderately fast
> [on this ancient single-core 32-bit 1.5 GHz CPU] and count.
> 
> After fixing the stupid bug: ~24 frames/minute.
> 
> Changed array indexing from tuples to integers: ~24 frames/minute.
> 
> Merged the four arrays into one interleaved array: ~26 frames/minute.
> 
> [I'm thinking this might improve cache behaviour slightly... but
> apparently not much. Or just the cache isn't the bottleneck.]
> 
> Removed array bounds checks: ~27 frames/minute.
> 
> Removed adaptive time step: ~80 frames/minute.
> 
> [Well there's a surprise. Pitty the frames look awful...]
> 
> Added adaptive time step within RK4 as per Scott's suggestion: ~1
> frame/minute.
> 
> [Ah, but now the "error" is being calculated in a totally different way,
> so the error tolerance required to get a good image may have changed.
> What happens if I use a less strict tolerance?]
> 
> Turned down error bound: ~68 frames/minute.
> 
> [I made the maximum permissible error 10,000x larger, and the images
> still come out looking flawless. But they render vastly faster...!]

A clearer (for a human) way to see the speed is to show the inverse of that:
seconds per frame.

Post a reply to this message