rk4 algorithm sent to you chuck

[Errol_Summerlin]

Just wanted to make sure chuck saw he has a private message. I sent him the files via pm.

For those that are unaware, chuck asked me to write a 4th order runge-kutta integrator for the back end physics engine. It should eliminate orbital drift and allow gravity to be put back in. I also included two other algorithms that are less accurate including a simple Euler's method and a midpoint method that is slightly more accurate. The less accurate methods might be used for unimportant objects in order to save on computing power (rk4 is roughly 5 times more expensive than Euler)...but many more times accurate.

Hopefully, I understood what he wanted from me and what is there will be useful.

[chuck_starchaser]

Yes, I got a pop-up as soon as I logged in. I think the first time you messaged me phpBB got confused because you deleted some messages.
Looks nice and tidy. Had a first quick look and messaged you back with questions and comments.

[Errol_Summerlin]

just an update, linear position, velocity, and acceleration components of the rk4 algorithm are fully implemented and tested. I am still testing rotation and orientation integration. Don't listen to that gaffer on games guy. He doesn't understand quaternions at all!. I wasted 3 days trying to follow his non-sensical routines and make sense of them. So far, my implementation is working correctly, but I haven't finished testing yet

[chuck_starchaser]

Got your files on svn. Sent you a PM with the details.
For public viewing, username is username; password is password:
https://svn.wcjunction.com/code_projects/

[chuck_starchaser]

Templetized version of Vector.h committed as vec3.h and vec3.cpp.

[chuck_starchaser]

Progress continues, btw. We're at r12.
Templetizations of vector and quaternion classes done, plus a cleanup of class matrix; though these
updated classes are in new files with new names and not yet being used.
Got started in the cleanup and modernization of the physics and RK4 parts that Errol wrote (the
other files were from other sources); which mostly involved a major cleanup of Mathematics.h, a
totally EVIL file; --though it will not be used; but for now just to make sure the physics part is
compiling and running without errors.
Next step (tonight) will be to replace the old vector, quaternion and matrix classes with the new.
Next step after that will be to refactor the physics, to separate concerns.
And somewhere along the lines we need to figure out whether we use angle and axis instead of
quaternions, or in addition to... Note: The problem with quats is that they only represent rotations
in a sphere, but can't represent multiple turns, which can be a problem for a fast spinning entity
that we simulate seldom but may need to interpolate often for graphics frames.
And somewhere along the line I need to start working on the new scheduler.

[klauss]

And somewhere along the lines we need to figure out whether we use angle and axis instead of
quaternions, or in addition to... Note: The problem with quats is that they only represent rotations
in a sphere, but can't represent multiple turns, which can be a problem for a fast spinning entity
that we simulate seldom but may need to interpolate often for graphics frames.

Just a wild idea, but:

Perhaps a better solution would be a representation with quaternions as rotational steps and a scalar specifying a stepping speed.

So if I have a quaternion that rotates along the (x,y,z) axis 820°/s, I could use the quaternion <cos(102.5°) + xi + yj + zk> with a stepping speed of 4 steps/s. I'll note it <cos(102.5°) + xi + yj + zk>@4

If I compose two rotations, I first find the "GCD form" of that representation:

<cos(a) + xi + yj + zk>@s * <cos(a') + x'i + y'j + z'k>@s' =
(<cos(a/s') + xi + yj + zk> * <cos(a'/s) + x'i + y'j + z'k>)@(s*s')

This would mean that quaternions+stepping speed should always interpolate nicely (ie: it would acelerate/decelerate if all that changes is the stepping speed), and the representation would still be compact, weighting only at 5 floats (yep, the stepping speed could be real).

Interpolation would be accomplished by interpolating the step quaternions to their "stepping power" (which is the elapsed time times the stepping speed). Quaternion exponentiation is also rather cheap.

[chuck_starchaser]

Yeah; this may be the ticket; and I think Errol was thinking along the same lines; not 100% sure.
But, speaking of the devil, why the devil do we need to interpolate anything at all?
If we have physics frames A and B, 0.07 seconds apart, and a graphics frame 0.04 seconds after A,
and want to interpolate x, it boils down to,

Code: Select all

float delta_x = B.x - A.x;
float delta_t = B.t - A.t;
float dx_dt = delta_x / delta_t;
float our_x = A.x + 0.04 * dx_dt;

The first three lines of code are superfluos. In fact, they involve a division, so about 90% of the
cpu time there is superfluous. Why? Because we already know dx_dt, from the physics; if we only
cared to read its value...

[Errol_Summerlin]

@chuck the delta_x/delta_t that you get from the code you gave is NOT going to be the same as either v_0 or v_1 (the state values of v). It will, in most cases be somewhere in between, but not necessarily the midpoint value. Despite having v, and a at the endpoints of the region, you need to compute the average dx/dt over the interval that gets you from state 0 to state 1. The code you quoted, is a very simple linear interpolation.

See, this is kind of what I was thinking.

The way I set up the rk4 algorithm, the state carries with it the "time", so the future state knows the time elapsed since the last state and it knows the new angular_velocity and angular_acceleration values, but it doesn't know the values in between.

However, a standard interpolation routine does the following for linear motion.

x1-x0=v*t+.5*a*t^2 +j*(t^3)/6
v1-v0=a*t+j*(t^2)/2
a1-a0=j*t

j=da/dt, sometimes called the "jerk"

since you know x0,x1,v0,v1,a0,a1, and t=elapsed time for the whole time step, you can solve for v,a, and j easily. Once you have v,a, and j, the equations can be written as

x(t')=x0+v*t'+.5*a*t'^2 +j*(t'^3)/6
v(t')=v0+a*t'+j*(t'^2)/2
a(t')=a0+j*t'

where t' is any value of t during the time step. And, x(t)=x1 is guaranteed because of how we solved for v,a, and j.

For angular motion, this works exactly the same way except that in quaternion representation for orientation, the number of rotations is lost, so
2*pi*n+theta_1-theta_0=v*t+.5*a*t^2 +j*(t^3)/6

where v and a are now angular velocity and angular acceleration values. We have introduced another unknown into the equations, n. But we still only have 3 equations. However, j=(a1-a0)/t is unchanged and the solution for a is also unchanged.

That just leaves 2*pi*n+theta_1-theta_0=v*t+computed_value

This gives a bunch of different values for v(remember v is the average value of angular velocity used for interpolation): one for each potential value of n. But, we know that the average value of v should be somewhere between v0 and v1. This will greatly restrict the potential values of n. In fact, in most cases, it will leave only one possible value of n. In the rare cases where there are multiple values of n that result in v0<v<v1 (this requires v1/v0>(n+1)/n, where n is the actual number of revolutions)we can even compute the average value of v(t')=v0+a*t'+j*(t'^2)/2 over the interval 0 to t to compute the average value of v and determine the value of n that results in v closest to the average value.

That solves the problem, but it is just a TON of computing. It seems for more simple have the physics keep track of n in the state. The physics.h file know exactly how many rotation it took, this info is just lost because the 1st element of the rotation quaternion is cos(delta_theta/2). If instead, I just return an (angle,axis) object that represents the orientation_delta rather than(or in addition to) the end result orientation quaternion, the entire problem disappears. Keeping track of delta_theta%(2*Pi) also works.

@klauss, the way I see it, that works, but you have to always have a check in to make sure the delta_theta per step is < 2*Pi. If you start getting more and more delta_theta, you will need to keep changing your stepping speed. The solution of just carrying n around in the state such that future_state.n = number of complete revolutions between currenstate.orientation and futurestate.orientation.

[chuck_starchaser]

Errol, you're da man. Your formulas go far over my head, though I wasn't entirely unaware that my
dx_dt proposition could run into trouble.
Trying to speak my mortal language, what vegastrike does now is interpolate positions linearly.
About a month or two ago, Klauss suggested we could interpolate velocities instead, to get smoother results.
You're probably suggesting interpolation of accelerations, by the mesmerizing looks of the math involved.
Is this it?
I think interpolating speeds would be orders of magnitude better than what we got now, and a lot cheaper to
compute. No?
After all, interpolations are for graphics only; no cumulative physical consequences. No cumulative errors, period.

EDIT:
No, I'm starting to understand what you're talking about, basically about which dxdt I use: old frame or new frame?
Well, how about we keep RK4's mid-point data stored somewhere, and use that?
So, to interpolate speed between frames A and B, we don't even look at A and B; we just consider the position,
velocity and acceleration at AB_midpoint, and assume acceleration constant.
I think the problems you foresee could be nullified by making it a rule of the software that wherever forces are
changing fast, that unit will get faster simulation steps.
In other words, the simulation speed priority should be proportional to force deltas and inversely proportional to
distance to the camera; or some such.

EDIT2:
This is a bit OT, but I've been meaning to say it for a long time; gotta get it out of my system:
We should NOT assume that all motions are necessarily modeled as mesh motions; if we do.
Case in point: Fast rotating objects whose geometry is a perfect circular extrusion don't need to be rotated; and
we can rotate their textures, instead. Suppose we were coding a car-racing game: Car wheels could be frozen
solid in place, and if they are cylindrically projected onto the UV plane and spanning the entire x or entire y and
the texture used tilingly, then all we need to do is add offsets to the uv coords in the shader, and the wheel will
seem to rotate. This would have the fringe benefit that we could implement motion blur on the texture simply
by fetching it anisotropically.
It would have yet another fringe benefit: We could rotate the diffuse, specular and normalmap textures, but not
rotate the AO or PRT's, so the shadows of other car parts onto the wheels stay where they are.
I think this car-wheel solution might transpose pretty well to a rotating station. In fact, if there are geometric
detail that do need to rotate, they could rotate even as the less detailed parts of the station have their texture
spun around, instead. Not that it really matters all that much, either way; but just to keep in mind.

[chuck_starchaser]

Errol: EJphysics is now working with the new files/classes (vec3.h/cpp, tmat.h/.cpp, quat.h/.cpp, and mat4x4.h/.cpp).
Output:

Code: Select all

1, 0, 0, energy =-0.5
0.995004, 0.0998333, 0, energy =-0.5
0.909505, -0.415626, 0, energy =-0.500014
0.579116, -0.815176, 0, energy =-0.500028
0.097231, -0.995176, 0, energy =-0.500042
-0.406545, -0.913506, 0, energy =-0.500056
-0.801896, -0.597225, 0, energy =-0.50007
-0.99108, -0.131981, 0, energy =-0.500084
-0.931716, 0.362643, 0, energy =-0.500098
-0.644249, 0.764521, 0, energy =-0.500111
-0.203546, 0.978808, 0, energy =-0.500125

Added a "build_all" script that compiles AND links to an EJphysics binary.

There is a division by zero tho, vec3<double> / 0.0 ... not sure where. Probably an uninitialized moment of inertia.

Code: Select all

1, 0, 0, energy =-0.5
EJphysics: vec3.h:347: vec3<double> operator/(const vec3<double>&, double): Assertion `not_zero(s) ||! "trying to divide a vector by zero"' failed.
Aborted

I simply commented out the assertion on line 366 of vec3.h and the test program worked; but you might want to
un-comment it back in and look into the matter.
r13

EDIT:
Ran the uncrustify code-formatter on all the current files; and deleted the old ones no longer used.
r15
EDIT2:
Removed some unused code and did a bit of manual cosmetic fixing here and there.
r16

[chuck_starchaser]

Started to refactor the physics code.
Errol, I've got rid of the distinction between "state variables" and "state derivatives", as I thought the split should be 3-fold
if done at all, adding "state second-derivatives". Anyways, having all state-variables in one struct seems clearer to me, and
as a fringe benefit now we got a power-of-two size structure (if you remember that I added a fourth var to vec3 for
padding). This is what it now looks like:

Code: Select all

class RigidBodyPhysics
{
public:
    struct ObjectState
    {
        vec3< double > acceleration;
        vec3< double > angular_acceleration;
        vec3< double > velocity;
        vec3< double > angular_velocity;
        vec3< double > position;
        quat< double > orientation;
        vec3< double > tensor_of_inertia;
        double mass;
        double inverse_mass;
        double mdot; //rate of mass gain (<0 for loss)
        double time; //current time
    };
    ...................

32 variables exactly, when added up.
And as you see, I changed scalar_moment_of_inertia to tensor_of_inertia; initialized to (1,1,1) for now.
Everything is double precision, for now; I figured might as well worry about optimizations later; let's get
it to work for now (and I'm not even sure we CAN avoid double precision anywhere; long story).
But I made all auto variables in rk4Step() static, to avoid huge stack allocations at each call.
I also began to change uses of operator + into += forms; --i.e. where,

Code: Select all

c = b + a;

I've started substituting,

Code: Select all

c  = b;
c += a;

Which a) avoids calls of constructors, b) makes the code (in this case) a bit more readable, as many of
the lines are rather long, and c) maps more directly to SSE intrinsics, --which I want to use once it's
all in the engine and working.
r17

EDIT:
Added -O2 to the compile flags in the build_all script.
r18

[Errol_Summerlin]

woah...slow down bro.

I gotta go see what you did so I can find the source of the divide by zero and make sure it still works as intended. I had the derivatives and the variables separate for a reason. It has to do with the rk4 look-ahead. When rk4 does its look ahead partial steps, the only thing that gets returned back are the derivatives at that point. So, returning the position and orientation is superfluous. I separated them so I wouldn't have to be constantly passing around variables that are not getting used.

Also, you have to be very careful using += with rk4 because you need to preserve the initialState through each of the rk4 steps.

Using += will increment the initialState and produce an error in the 4th rk4 step which will gradually build up causing that subtle bug I told you about in private message.

*= may also be hairy with quaternion rotations as the order of the multiplication matters. *= is usually such that q*=w -> q=q*w, but that is actually wrong. It has to be q=w*q if you want to change the value of q through a rotation represented by w.

I had reasons for doing the things I did the way I did them. I wish you would have asked me why before changing the implementation so substantially.


P.S. I am having some trouble with svn up. It keeps saying that it is "skipping" the directory i want it to update with no explanation as to why. The only command I can get to work is co. Presumably that gets the latest version. I will keep trying to figure it out.

[chuck_starchaser]

woah...slow down bro.

I gotta go see what you did so I can find the source of the divide by zero and make sure it still works as intended.

I think it's probably something not getting initialized, like after each step you copy current to previous, and next to current; but, what happens the very first time? I think you need to initialize current AND previous, for the simulation to start correctly, or something like that.

I had the derivatives and the variables separate for a reason. It has to do with the rk4 look-ahead. When rk4 does its look ahead partial steps, the only thing that gets returned back are the derivatives at that point. So, returning the position and orientation is superfluous. I separated them so I wouldn't have to be constantly passing around variables that are not getting used.

I realized that's why you'd separated them; but considered it irrelevant, because these routines shouldn't be passing entire states by value, anyways. That's a lot of time spent copying data to and from the stack.

Also, you have to be very careful using += with rk4 because you need to preserve the initialState through each of the rk4 steps.

Using += will increment the initialState and produce an error in the 4th rk4 step which will gradually build up causing that subtle bug I told you about in private message.

No, man, I wouldn't change the semantics in any way. As per the example I gave you, where,

Code: Select all

c = a+b;

substitute,

Code: Select all

c = a;
c += b;

*= may also be hairy with quaternion rotations as the order of the multiplication matters. *= is usually such that q*=w -> q=q*w, but that is actually wrong. It has to be q=w*q if you want to change the value of q through a rotation represented by w.

Ditto. I wouldn't change the order of the operations even if they were commutative; --just a matter of family honor.

I had reasons for doing the things I did the way I did them. I wish you would have asked me why before changing the implementation so substantially.

I can change back anything that needs to be; don't get jumpy; but wait till I'm done with it, first; I think you'll see the light.

P.S. I am having some trouble with svn up. It keeps saying that it is "skipping" the directory i want it to update with no explanation as to why. The only command I can get to work is co. Presumably that gets the latest version. I will keep trying to figure it out.

Hmmm...
Did you try svn cleanup?
What's the exact message?
SVN is not terribly helpful with its error messages; that's for sure.

[charlieg]

Errol_Summerlin: since SVN is a source control management tool, you can always see what changes were made, in what order, and compare different revisions - thus allowing you to work at a revision you still understand and (attempt to) forward port anything you change.

You just need a decent SVN client. Eclipse's SVN plugin is pretty decent for this kind of thing. (But that would require you to install Eclipse.)

[chuck_starchaser]

Almost done refactoring the physics.
r21

• File test.cpp added, which has the testing code (main())
• Big functions formerly inlined in EJphysics.h moved to EJphysics.cpp
• RK4helper used to do nothing more than call EulerHelper; so now there's a single helper() function
• NO structs whatsoever are passed by value to functions or returned from them, anymore
• A number of function variables were eliminated by writing directly to class states
• ALL remaining function auto variables are now static

The last step I'm about to take is the most difficult:
Right now there are three states:

Code: Select all

    ObjectState futurestate, paststate, currentstate;

The RK4 code works from paststate and currentstate, and writes to futurestate.
Then, a sweep of copying happens...

Code: Select all

    void AdvanceState()
    {
        paststate    = currentstate;
        currentstate = futurestate;
    }

I want to --just try-- to avoid the copying, and reduce memory footprint, by having only two
states:

Code: Select all

    ObjectState  past_and_future_state, current_state;

That is, when the RK4Step code writes the futurestate, it will be overwriting the paststate. The hard
part is making sure no part of the paststate is overwritten before it is used.
Then, instead of sweeping through all the objects copying, we sweep quickly toggling a boolean, or
swapping two pointers; I'm not sure yet.
But so, in reality there will be,

Code: Select all

    ObjectState  A, B;
    ObjectState * pPastAndFutureState, const * pCurrentState;

That is, two "real" states, and two pointers (which we can then swap)...

Code: Select all

    void init()
    {
        pCurrentState       = &B;
        pPastAndFutureState = &A;
    }
    void AdvanceState()
    {
        static ObjectState * pTemp = pCurrentState;
        pCurrentState = pPastAndFutureState;
        pPastAndFutureState = pTemp;
    }

We've made 3 copy operations 8 bytes each (assuming 64-bit pointers), or 24 bytes; versus, let me see...
32 x 8 x 3 = 768 bytes. This is per-object simulated, so it matters.

But having these pointers to states would cause all accesses to be derefenced; so I'm not sure if having
swappable pointers is the right way, or having two identical sets of RK4Step() and helper functions, and
choose between one set and the other based on a boolean at the entry point.

[chuck_starchaser]

Well, after a first study of the code, I cannot find a place where previous state is used in computing next state,
at all. Looks like the purpose of previous state is merely for benefit of interpolation by the clients; so this will
be easier than I thought.
Thing is, though, the previous state will no longer be available to clients; so clients MUST use current state
only. Interpolation could be linear (based on position at "current" time (as per time-stamp) and velocity; or
we could get a bit more sophisticated and use acceleration.
But so, just to clarify my intent, when rk4 is running along a container of objects doing physics updates, it
uses current data, and overwrites "previous" with "next", so the previous/next state is not available to clients;
only current state is. Once all physics updates are done for a given frame, prvious/next becomes current, and
current becomes previous/next.
This swapping of roles has to happen when rk4 is finished computing next AND clients are done using current.
In other words, the sweep through units toggling booleans has to be semaphored. What I'm talking about is
multithreading, of course. When rk4 is reading from current and updating next state, it changes nothing of
current, and therefore current can be used concurrently by clients; but if clients are still using current data
when rk4 has already finished, or viceversa, the state swap has to wait until all parties are done with current.
Once this happens, the physics has to mark the container as off-line, quickly sweep through it toggling the
boolean flags that determine which state is current, which is prev/next; then make the container available
to clients, again, and can now run rk4 on it, again, when it wishes to.
Since the toggling of the boolean flags happens immediately after the rk4 run, and another rk4 run may not
happen for a while, chances are all users of current data will be done once the next rk4 is done.
Alternatively, the clients could check for a "wish to swap states" flag after each unit they process, so as to
let the physics thread do the swap as soon as it wants to, without delay.
And you might ask, if all units in a given container are being simulated synchronously, wouldn't it make more
sense to have a single, per-container flag?
Good question!
The problem is that units will occasionally migrate from one container to another; so although they are being
simulated synchronously, it is not guaranteed that their current-prev/next states are the same.
But perhaps a smarter idea would be to have a per-container toggle flag, and then a per-unit mask flag that
is xor-ed with the container flag when brought into it. Then, rk4 uses the container flag xor the unit flag to
decide which state is current. This would save us a whole sweep through the container; AND, toggling a single
flag is an atomic operation that doesn't need to be mutexed.

[Errol_Summerlin]

On how to interpolate. I would NOT recommend using past states for purposes of interpolation. This can be problematic for esoteric mathematical reasons related to the discontinuity of various derivatives due to user input (thrust on, thrust off, thrust on, thrust off).

You suggested that we could "carry around" the deltas from the rk4 step. While this is possible, it is not recommended either because it is not self consistent.

Consider an object initially stationary that experiences acceleration = 10 for a time step of 1 second. The "delta" for x is 5. The "delta" for v is 10. The delta for a is zero.

So, if we compute x(.5), we get 2.5 with a simple linear interpolation. v(.5)=5 and a(.5)=10. So, if the thing traveled 2.5 is .5 seconds, its average velocity during that .5 second interval must have been 5, but v has only just reached 5, so how could the average value of v have been 5. It was actually 2.5. The value of x is changing faster at the start of the interval than it should be and it will change slower than it should at the end of the interval. The reason this happens is because the interpolation is not self-consistent. You need to treat all the data points you have together.

Consider the same situation, but now we interpolate self-consistently
We start with j=(a1-a0)/t. This is a simple linear interpolation of the acceleration. We have no higher-order derivatives to use, so this is the best we can do.

j=0 because a1=a0
Then, we have v1=v0+a*t+j*t^2/2. j=0, so a=(v1-v0)/t=10. We could kind of figure this out since a0=10 and a1=10, but this is just a simple example to show the process.
Now, we write x1=x0+v*t+.5*a*t^2+j*(t^3)/6=x0+v*t+5*t^2. We know from our rk4 that x1-x0=5 for t=1. So, 5=v+5. This gives v=0.

So, to summarize, v=0,a=10,j=0 are the values we go for our coefficients. Now we plug them into the equations of motion

a(t)=a0+j*t=a0
v(t)=v0+a*t+j*(t^2)/2=10*t
x(t)=x0+v*t+a*(t^2)/2+j*(t^3)/6=5*t^2

So, now, at t=.5, we have a=10,v=5,x=1.25 as we should. at t=.8, a=10,v=8, x=3.2 and at t=1, a=10,v=10,x=5 just as rk4 told us we would.

It is essentially doing a linear interpolation on the acceleration, but carrying that down to the lower order derivatives. If you don't do this, you will get "lurch" in between time steps as the initial changes in x will over account for the acceleration with later times during the physics time step under accounting for the acceleration. You still end up in the same place either way, but the latter method gets you there more smoothly.

Hopefully that makes more sense than the first time around. I usually find stepping through the algorithm helps.

As for the state swap dealio. Are you sure it is worth going through all this trouble to eliminate one state variable from floating around in memory? I don't understand most of programming lingo you are talking about, but it seems like you are putting a lot of restrictions on how users can interface with the class and that this change will make various other parts of the simulation more complicated. Anyway, its your call. I don't see any reason to need both the future and past states at any point in time, but it just seems like it may not be worth it.

[chuck_starchaser]

On how to interpolate. I would NOT recommend using past states for purposes of interpolation. This can be problematic for esoteric mathematical reasons related to the discontinuity of various derivatives due to user input (thrust on, thrust off, thrust on, thrust off).

You suggested that we could "carry around" the deltas from the rk4 step. While this is possible, it is not recommended either because it is not self consistent.

I don't remember if I said "deltas", but if I did, what I meant was "derivatives".
If, for interpolation purposes, I look at a unit's "current" state (ONLY), and its time-stamp is 777.7, and position is X, and velocity is V, and acceleration is A, and the current, interpolation time is 777.8, then:

• dt = 777.8 - 777.7 = 0.1
• Vdelta = dt * A
• newV = V+Vdelta
• averageV = 0.5 * ( V + newV )
• deltaX = dt * averageV
• newX = X + deltaX

That should be good enough for interpolation. Excessive, I'd say.
Consider this: You have a

Code: Select all

 double dt

as input to rk4Step(), but the truth you probably don't know is that the finest time resolution we can get out of system timers is 1 millisecond. So, the time-graininess of interpolations at 60 FPS is 1 millisecond in 16.67, or roughly 6%. Does it make sense to worry about 5th derivatives, anymore?

Consider an object initially stationary that experiences acceleration = 10 for a time step of 1 second. The "delta" for x is 5. The "delta" for v is 10. The delta for a is zero. Interpolation is destined to be for the birds, no matter.

So, if we compute x(.5), we get 2.5 with a simple linear interpolation. v(.5)=5 and a(.5)=10. So, if the thing traveled 2.5 is .5 seconds, its average velocity during that .5 second interval must have been 5, but v has only just reached 5, so how could the average value of v have been 5. It was actually 2.5. The value of x is changing faster at the start of the interval than it should be and it will change slower than it should at the end of the interval. The reason this happens is because the interpolation is not self-consistent. You need to treat all the data points you have together.

Consider the same situation, but now we interpolate self-consistently
We start with j=(a1-a0)/t. This is a simple linear interpolation of the acceleration. We have no higher-order derivatives to use, so this is the best we can do.

j=0 because a1=a0
Then, we have v1=v0+a*t+j*t^2/2. j=0, so a=(v1-v0)/t=10. We could kind of figure this out since a0=10 and a1=10, but this is just a simple example to show the process.
Now, we write x1=x0+v*t+.5*a*t^2+j*(t^3)/6=x0+v*t+5*t^2. We know from our rk4 that x1-x0=5 for t=1. So, 5=v+5. This gives v=0.

So, to summarize, v=0,a=10,j=0 are the values we go for our coefficients. Now we plug them into the equations of motion

a(t)=a0+j*t=a0
v(t)=v0+a*t+j*(t^2)/2=10*t
x(t)=x0+v*t+a*(t^2)/2+j*(t^3)/6=5*t^2

So, now, at t=.5, we have a=10,v=5,x=1.25 as we should. at t=.8, a=10,v=8, x=3.2 and at t=1, a=10,v=10,x=5 just as rk4 told us we would.

It is essentially doing a linear interpolation on the acceleration, but carrying that down to the lower order derivatives. If you don't do this, you will get "lurch" in between time steps as the initial changes in x will over account for the acceleration with later times during the physics time step under accounting for the acceleration. You still end up in the same place either way, but the latter method gets you there more smoothly.

Hopefully that makes more sense than the first time around. I usually find stepping through the algorithm helps.

I'll try to go through it, later; just going to say for now that I think you're talking about simulation, rather than interpolation. Interpolation is not something terribly important; it's just for the graphics, which often will run at a faster rate than the physics. If we're going to do something as complex as rk4 for interpolation, we might as well run the physics inside the graphics loop; but that's what we're trying to avoid by interpolating.

As for the state swap dealio. Are you sure it is worth going through all this trouble to eliminate one state variable from floating around in memory?

Completely sure.
It's not to reduce memory footprint, mainly; --though it IS a concern--, but to avoid COPYING states.

I don't understand most of programming lingo you are talking about, but it seems like you are putting a lot of restrictions on how users can interface with the class

None whatsoever.

and that this change will make various other parts of the simulation more complicated.

Nope. Totally transparent.

Anyway, its your call. I don't see any reason to need both the future and past states at any point in time, but it just seems like it may not be worth it.

Done.
r23

Not sure if you've figured out svn or not, yet. Just in case, here are the files gzipped:
http://wcjunction.com/temp_images/RK4.tar.gz

All yours; I'm going to work on the scheduler, now.

[chuck_starchaser]

By the way, here's a list of things that I think need yet to be done with the physics code:

• As I was saying to you in my last PM, I think the mass loss thing could be simplified. The only thing we REALLY need RK4 for is gravity --so that docking at stations works. Other forces are much less relevant because the user doesn't know how hard the engines are pushing at any given time. We don't know that they produce a perfectly even thrust. So the fact that other forces benefit from RK4 is icing on the cake, but not really important. And it so happens that mass affects accelerations from all forces EXCEPT gravity. So, your getForceTorqueandMassLoss() function could be changed to getForceTorqueandMass(), and that would be by far sufficient. Right now we're applying RK4 to mass, which seems to me ridiculously excessive, considering that mass changes ever so gradually (except when you launch missiles, or drop cargo)...

  Code: Select all

  futurestate.mass+=(1.0/6.0)*dt*( d1.dm_dt+2.0*( d2.dm_dt+d3.dm_dt )+d4.dm_dt);

  !!!
• Another thing we'll need to do with this physics code is switch from simulating speeds to simulating moments; then make speed just equal to moment/mass. This is important to me because I have another plan for mass: Mass Hiding as an FTL technology. Hiding 99.999 of inertial mass should multiply your speed by 10,000 without spending any fuel or propellant, and that's the ONLY way, I believe, we can explain such ridiculous speeds AND accelerations as are needed to travel to another planet in a matter of minutes. So, for this I would have an inertial mass dividing factor that would look like the SPEC multiplier in the dashboard. So, if we simulate momentum, then we divide it by this divided inertial mass, resulting in higher speeds. The really neat thing about this solution is that if we leave gravitational mass unaffected, the ship's trajectory will bend around planets as if this was some kind of "time compression".
• Another thing that needs to be done; but this I'm assigning to myself, is to have a "physics_object" that can be as simple as a rigid body, or as complex as a graph of weights and springs and motors and whatnots. Not that we're going to implement this any time soon, but just so that the class system is organized properly from the start. The work I just started on the scheduler assumes there are physics objects in a container (array allocator), but doesn't assume they are rigid bodies.
• The division by zero problem is still there. The commented out asserts are in vec3.h, lines 383, 407, 408 and 409. I'm sure it's something not initialized.
• We need code that moves planets around stars using simple Kepler/Newton rules, but that as you enter a system smoothly changes to simulating all the planets and stations; then goes back to the simpler formula after you leave for another system.
• This is for myself: Harmonizing physics code with scheduler (reminder to self: freq_bin toggler xor-ing).

[chuck_starchaser]

Something that just occurred to me:
I think I would move gravitation into the physics code; --i.e.: separate it from other forces.

Why? Because,

• a) It would be faster for us to do it, because we have very few branches, smaller code, and better data management than the Unit class.
• b) In my evolving physics scheduler, I have an unused (9th) frequency slot I'm thinking of using for planets, and if I do this, we have all our planets in one place and can compute gravitational forces rather quickly, rather than the way the Unit class would probably do it (scanning across all units for units of type "planet")
• c) Most importantly, because gravity needs to be computed at double precision, whereas all other forces can be computed at single (float) precision.
• d) Gravitation, unlike most other forces, is not conditional on anything, so it is easily separable from Unit.

[Errol_Summerlin]

I don't remember if I said "deltas", but if I did, what I meant was "derivatives".
If, for interpolation purposes, I look at a unit's "current" state (ONLY), and its time-stamp is 777.7, and position is X, and velocity is V, and acceleration is A, and the current, interpolation time is 777.8, then:
dt = 777.8 - 777.7 = 0.1
Vdelta = dt * A
newV = V+Vdelta
averageV = 0.5 * ( V + newV )
deltaX = dt * averageV
newX = X + deltaX
That should be good enough for interpolation. Excessive, I'd say.

If you ONLY look at the current state, you will not end up where you need to at the end. In fact, THIS is actually simulation. Its a simple Euler simulation, and the position your object is at a time dt later using this interpolation will not match what was calculated via the rk4 physics calculations.

IF, as I believe you meant to say, you use the rk4 derivatives, (x1-x0)/deltaT, for example, then you will end up in the right position at the end, but your velocities and accelerations will be discontinuous. You will have "lurch" where the velocity suddenly changes at each physics step but is constant in between.

Remember, x(t) DETERMINES all of the subsequent observed derivatives. And if you make x(t) and simple linear interpolation such that x(t)=((x1-x0)/deltaT)*t, then your x will be continuous, but at each point,(x0,x1,x2, etc.), the velocity will instantaneously change from (x1-x0)/deltaT to (x2-x1)/deltaT. The values you compute for the velocity are meaningless in terms of the motion of the object. The observed motion is determined solely by x(t). The actual motion the ship experiences will be constant velocity motion with sudden changes in the velocity at each physics step. Yet, of course, the velocity you calculate will be changing over the interval. This is the inconsistency I am talking about.

I'll try to go through it, later; just going to say for now that I think you're talking about simulation, rather than interpolation. Interpolation is not something terribly important; it's just for the graphics, which often will run at a faster rate than the physics. If we're going to do something as complex as rk4 for interpolation, we might as well run the physics inside the graphics loop; but that's what we're trying to avoid by interpolating.

No, in the algoirthm I presented, I already know x1,v1, and a1. These were calculated via the rk4 algorithm. But since the frame rate is faster than the physics simulation frequency, we need to figure out what path the object took to get from x0 to x1, unless acceleration was constant, and 0 over the entire time step, a simple linear interpolation is not what the motion actually was, and it isn't the best we can do with the information we have.

If all we had was x0 and x1, then a linear interpolation would be the best we can do, but we also have the derivatives at the end points.

Consider an object in a cicular orbit with time step .25*the period.

The rk4 algorithm calculates positions (1,0),(0,1),(-1,0),(0,-1)...all perfectly correct values.

With a linear interpolation of position, the observed position of the object in game will be along a diamond with vertexes at the points that the physics engine actually calculated. The "observed" velocity is then obviously a vector pointing along the direction of motion. At each vertex of the diamond, the velocity that is observed in game changes discontinuously.

Furthermore, the linear interpolated velocity, which you calculate via v=v0+average_acc*t is not consistent with the direction of motion which the object is experiencing. The calculated value varies linearly in time, but the observed value, based on the linear interpolation of x is a constant.

For example, at (1,0), the unit vector is (0,1). The velocity you are caculating then changes from pointing in the y-direction to pointing in the x-direction over the interval. However, the objects actual trajectory is just a straight line. This is another demonstration of the inconsistency of linear interpolation.

So, hopefully I have demonstrated that we ONLY need an interpolation for x(t) (the velocity and acceleration observed by players is purely related to the positions observed at different times), but a linear interpolation created discontinuities in the derivatives and it doesn't use all the information we have. It ONLY uses x1 and x0. With two points, the best you can do is a linear interpolation. With 6 data points the best you can do is a 5th order polynomial.

Think of it like this. For a linear interpolation what you are actually doing is saying x1=x0+someconstant*t and you are solving for that constant. 1 equation, 1 unknown..ez pz. But, if you compute the derivative of this equation of motion, you get v=someconstant. But v should be changing over the interval. The derivative should equal v0 at t=0 and V1 at the end of the interval, but it doesn't. It is someconstant over the whole interval. The acceleration should be a0 at t=0 and a1 at the end of the interval. It is 0 over the entire interval. In other words, a linear interpolation of position does not match the values of any of the derivatives at the end points.

I should also note that there is no point to interpolating velocity or acceleration values. Those are only important for the physics. The actual velocity that a player observes will be based on (x(t1)-x(t2))(t1-t2) where t1-t2 if the time between frames. You could, for example, display a linear interpolated velocity on his hud, but it wouldn't match the actual time derivative of his motion because that would be a constant over the interval.

So lets say i want motion in which the velocity does not change discontinuously.

then my constraints are
x(0)=x0
x'(0)=v0
x(deltaT)=x1
x'(deltaT)=v1

So, I write my equation of motion

x(t)=x0+v0*t+a*(t^2)/2
x(0)=x0 as I required
and x'(0)=v0 as I required

Now, I need to make sure that the end points are correct.

x(deltaT)=x1=x0+v0*deltaT+a*(deltaT^2)/2
x'(deltaT)=v1=v0+a*deltaT

but now we have two equations and only one unknown. It may be, and almost always is, impossible to satisfy both equations simultaneously for a constant a. So, we add another term, the jerk, to the equation of motion and solve for a and j. The results is a= (-2*(2*deltaT*v0 + deltaT*v1 + 3*x0 - 3*x1))/Power(deltaT,2) and j= (6*(deltaT*v0 + deltaT*v1 + 2*x0 - 2*x1))/Power(deltaT,3).

This satisfies all 4 of our constraints. However, notice now that x''(0)=a=(-2*(2*deltaT*v0 + deltaT*v1 + 3*x0 - 3*x1))/Power(deltaT,2). This is almost never going to be a0. In other words, now, instead of the velocity being discontinuous, the acceleration is discontinuous. While it does change in time, it does not change to match the end points. This can be fixed adding two more derivatives position the second derivative of acceleration (sometimes called snap) and the third derivative of acceleration (sometimes called crackle...I am totally serious), but I am not advocating that. It would be over kill.

For reference, this particular solution is called a hermite polynomial interpolation. It "looks" much better than a simple linear interpolation and is not terribly computationally expensive. The values a and j only need to be calculated once per physics step and they can then apply to any value in between. Furthermore, their formulas are given above. The resulting position equation is

x(dt)=x0+v0*dt+a*(dt^2)/2+j*(dt^3)/6 for dt=t-t0 with t<t1 and a=(-2*(2*deltaT*v0 + deltaT*v1 + 3*x0 - 3*x1))/Power(deltaT,2) and j= (6*(deltaT*v0 + deltaT*v1 + 2*x0 - 2*x1))/Power(deltaT,3) where deltaT=t1-t0

[chuck_starchaser]

Errol, I can't argue with the math.
What I was trying to say is that the purpose of "interpolation" is to get approximate value for finer subdivision.
It should be optimized for speed, though hopefully produce pleasing results.
I'm glad you say that the Hermite solution would be overkill. No kidding.
But you mention me as having suggested a linear solution... I did NOT. Look at my algorithm again:

Code: Select all

dt = 777.8 - 777.7 = 0.1
Vdelta = dt * A
newV = V+Vdelta
averageV = 0.5 * ( V + newV )
deltaX = dt * averageV
newX = X + deltaX

I'm using position, speed AND acceleration there.
You say that I'm suggesting a constant speed between simulation frames.
No; my code above assumes constant, non-zero acceleration.

You also mention that that is simulation, rather than interpolation. I agree. The whole point was, and is,
that there won't be more than one frame to interpolate. If the right name is "simulation", so be it.
What I was trying to say is that we have enough info to plot smooth points connecting our RK4 frames.
That's all we need to do.

In your example of 4 points around an attractor, this would produce 4 parabolas stitched together.
Not a circle; but a lot closer to a circle than a square.

And in any case, the example is a gross exaggeration of the actual situation. Our phyisics RK4 simulation
speeds will range from 16 frames per second to one frame every 8 seconds for things like planets far away.
That's far from a "4 dots around a circle" situation. Very far from it.

Anyways, at the beginning of this work you thought that RK4 wouldn't be necessary; and now you seem
to think it is necessary for everything, even for fuel consumption related mass changes, and want to
compute I'm not sure how many terms for the (should be quick) interpolations. Could we get back to
some happy middle ground? We needed RK4 for one main reason: So that we could bring back gravity
without the drift problems that caused it to be removed. We don't need accurate simulation of mass,
since gravitation and inertial (orbital kinematics) are mass agnostic. And interpolations ... err, simulations
using second term (derivative) will be far better than what we have now, and will rock.

What needs to be done now is to optimize the code a bit, add new scheduling, and integrate in the engine.

[Errol_Summerlin]

What you have posted there is simply not correct.

What are you using for A? How is it determined? Is it a1? a0? (a1+a0)/deltaT?. You said you were ONLY using the current state.

So, I assumed that means A=currentstate.acceleration
So, proceeding from there, we get down to
averageV=V+A*dt/2
and deltaX=V*dt+A*(dt^2)/2

First, this is a midpoint step. It is an integrator. It is not an interpolator. Second, when dt=physics time step, you won't get futurestate.position because futurestate.position was calculated with an rk4 integrator, and yours was calculated with a midpoint method integrator.

So, what will happen with this "interpolation" routine is that at every physics time step, the object will suddenly shift its position instantaneously to match the value computed by the rk4step.

Consider the circle problem, A=(-1,0), so deltaX=(0,1)*Pi/2+(-1,0)*Pi^2/8=(-1.23,1.57). It should be (-1,1) in order to end up at the place the rk4 algorithm told you you should be. So, after this time step is done "interpolating", the object suddenly shifts to (0,1) and starts the same errant trajectory again.

The motion would be smoother using a simple linear interpolator that just says x=t*(x1-x0)/physics_time_step.

In interpolation, there are two thing to consider about the values at the nodes. 1) are they continuous and 2) is it smooth (is the first derivative continuous). The algorithm you gave is neither. If it is working in the game, it is merely because your physics time steps are small enough that the sudden shift at the end of each time step is not perceived because it is small.

My objection is not to the simplicity or inaccuracy of the algorithm. It is that the algorithm is just plain wrong. And it is wasting cpu cycles doing something that was already done by the rk4 algorithm..integrating

If A is not currentstate.acceleration and is, in fact, (a0+a1)/2, then A=(-.5,-.5). deltaX=(-0.61,0.95). This is still NOT what the physics step gave you.

Interpolation is NOT integration. Interpolation is taking two known values of the function, and approximating what the values were in between subject to the constraint that the function has the right values at the nodes (points where the value is known).

A linear interpolation is continuous, but not smooth as the first derivative, velocity is not continuous. Using that may be sufficient if the physics time step is fast enough that the velocity lurches are small, but it may not be. I was suggesting, as Klauss did, that a linear interpolation of velocities would be smoother. Where v(t)=v0+t*(v1-v0)/physics_time_step. Then x(t)=x0+v0*t+(v1-v0)*(t^2)/(2*physics_time_step)=(1,0)+(0,1)*t+(-1,-1)*(t^2)/Pi. Evaluating this at t=Pi/2(the physics_time_step), we get (0.215,.785) instead of the expected (0,1). So, now the velocity is continuous, but we messed up the position interpolation. To have both satisfy the values at the nodes, the only options is a 3rd order polynomial for x(t). This allows us to match both the positions, x1, and x0 and the velocities, v1 and v0. This gives a smooth and continuous function for position.

Anyways, at the beginning of this work you thought that RK4 wouldn't be necessary; and now you seem
to think it is necessary for everything, even for fuel consumption related mass changes, and want to
compute I'm not sure how many terms for the (should be quick) interpolations. Could we get back to
some happy middle ground? We needed RK4 for one main reason: So that we could bring back gravity
without the drift problems that caused it to be removed. We don't need accurate simulation of mass,
since gravitation and inertial (orbital kinematics) are mass agnostic. And interpolations ... err, simulations
using second term (derivative) will be far better than what we have now, and will rock.

I don't think rk4 is necessary for every object, but once an object is using rk4, you can't just pick some variables out and use them as euler steps and others as rk4 steps without being VERY careful how you do it. Assuming mass is constant over the rk4 step is fine and then doing the euler step on the mass AFTER the rk4 step is done CAN be done, but when I suggested this, you seemed to say that we could use an euler step during the rk4 step..but that is what rk4 already DOES. It uses an euler step 4 times to calculate derivatives at future times. rk4 is a series of Euler steps followed up by an average of the results so when you say you want to compute the mass with an euler step, yet not have it remain constant during the rk4 step, that doesn't make any sense. I am opened to the idea of cutting out unnecessary computational steps, but I don't understand how you intend to do it, and I am not convinced it isn't just plain wrong. The same goes for the interpolation. If a linear interpolation works well enough that nobody notices, then just do a linear interpolation, but the next step up is a 3rd order polynomial fit(a hermite polynomial interpolation), and I was just trying to explain how one would do it because you seemed to be looking for a better interpolator. However, it seems that you don't need a better interpolator so much as you need an interpolator as opposed to an integrator, but maybe that is a communication problem. But from the information I have read, that appears to be the case.

I apologize belatedly if my tone sounds severe. It is 2:30 am and I need sleep. It is hard to relay complicated concepts like this over the internet.

[Errol_Summerlin]

One last thing, though we don't have any use for this, it may provide a reference point.

When multiple nodes are strung together via a series of interpolations between each set of data points such that the function is smooth and continuous, it is called a cubic spline. We only have two nodes, but we know the values of the derivatives at the end points and can require smoothness. This is the subtle difference between a hermite polynomial interpolation and a cubic spline.