rk4 algorithm sent to you chuck

Development directions, tasks, and features being actively implemented or pursued by the development team.
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

Errol_Summerlin wrote: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,
Nope; right now we only have one. If we really need two, then we'll have to bring back the copying of states (oldstate=currentstate); which I'm hoping we can avoid.
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.
I see.
What I was trying to do in my little code piece was to NOT depend on two ends; but just use one frame's data, including position, velocity and acceleration.
I KNOW it wouldn't be "correct".
I'm just hoping it is "good enough" (or can be made good enough), given the smallness of our physics steps.
But if you think it would not be "good enough", then we'll have to go back to copying states, have two
points to interpolate, and do the Hermite thing, I guess.
But we should quantify the problem; not just push for correctness to the point of doing rk4 yet again.

Furthermore, there's another reason why I wanted to avoid having 2 states and interpolating: In most cases, we'd be extrapolating, anyways. If the last physical frame for some capital ship was computed half a second ago, and the previous one 2.5 seconds ago, then you don't have "interpolation"; you have "extrapolation", and would extrapolation be any better than the little code I wrote using a single frame with velocity and acceleration? I doubt it... With extrapolation, the older of the two data-frames is too old to be too relevant.
Well, I just recently learned that our graphics are slightly delayed relative to the physics and collisions, so probably we do have "interpolation" in the case of objects that are simulated at the highest speed (16 simulation steps per second); but most objects would be simulated at lower speeds, so in most cases we'd be doing "extrapolation".
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

I've committed a first batch of code files for the scheduler. Note that I haven't tried to compile anything; this is just
text editor work, with probably gazillions of typos and major bugs; just for the curious to browse. The svn requires
username and password. If you don't have an account, you can login as username / password, literally.
https://svn.wcjunction.com/code_projects/scheduler
The only file there that I've worked on is bit_array, so far. The rest are just ideas, comments, etc. Ironically, I've not
much in the way of comments in bit array. Basically, there's functions for setting and clearing bits; and then the big
thing is the iterator. The iterator traverses the bit array, hopefully optimally, and pushes the offsets of bits it finds
set into a small queue.
The queue has two pull functions: One for pulling a "next" bit offset, and one for pulling a bit offset two places ahead
of the next, for prefetching purposes. The prefetch pull occasionally precipitates a search for a whole bunch of set
bits to refill the queue, when it's running out of bit offsets. The bit offset should be identical to the index of the
physical unit to simulate.
So the code is a bit complex but should be ridiculously fast.
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

@errol: so much I won't read (seems like a conversation between you and chuck anyways)...

About the <axis, angle> representation, how would you interpolate between two different such vectors?

Because <axis, angle> has no speed implied, and so interpolating naively would produce artifacts IMHO. Suppose a unit that has been spinning a long time in a certain way, so you have angle = <huge>*360° + <small>. Then you have the new rotation vector computed out of a collision that has angle = <small> - naive interpolation would "roll back" the huge rotation, causing a severe interpolation artifact.

What we care most is the rotational speed, rather... so we can interpolate speeds, with no discernible artifact. I was hoping (though I haven't proven it) that my 5-tuple would achieve that precisely.
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

I think errol already addressed the 5-tuple idea. It's not bad, if I'm understanding it correctly; but it seems to me a bit hard to manage. Basically, you'd settle for a quaternion and an "angle multiplier" --essentially, a quaternion power; right? But how do you decide how big an angle to allow --i.e. where to switch from 4-tuple to 5-tuple? I foresee if statements starting to crop up everywhere...
Right now, Errol's physics code is 100% if-free, and I'd like to keep it that way.
Unless we scale all our quaternios to one radian, and then use fractional 5th floats; but then we're going to be exponentiating most quaternions for nothing.
Then again, my suggestion of using angle and axis requires us to convert between quats and angle/axis; but that's probably less computationally involved than computing quaternion powers, all other things being equal.
In any case, the current question has metamorphosed during your absence, to whether to interpolate, or as Errol calls it, "simulate".
I kind of made up my mind about a week ago that interpolation was for the birds and decided to get rid of the old frame and just use one frame's data for "interpolation" ("simulation"), basically, if you want to "interpolate" x at time 777.8, and the last physics frame is time-stamped 777.7

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;
But Errol thinks it is too incorrect.
I'm aware of its incorrectness; but as I mentioned in my last post, I think interpolation, even Hermite, would be just as incorrect, since what we really do is extrapolate, in most cases, anyways.
If we do this, then we won't have two quaternions to interpolate between; so, problem solved :mrgreen:
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

chuck_starchaser wrote:...But Errol thinks it is too incorrect.
I'm aware of its incorrectness; but as I mentioned in my last post, I think interpolation, even Hermite, would be just as incorrect, since what we really do is extrapolate, in most cases, anyways.
If we do this, then we won't have two quaternions to interpolate between; so, problem solved :mrgreen:
Ehm... I thought the idea was to simulate one frame ahead in physics, and show stuff with a small lag (one physics frame lag) so we're always interpolating rather than extrapolating.

In any case, true second-order inter/extrapolation (interpolate speed and apply linearly variable speed equations for inter/extrapolation) isn't terribly incorrect IMHO, and error (of inter/extrapolation vs simulation) can be interpolated (linearly) just as well, to hide artifacts because of discrepance between second-order integration in the interpolator and rk4 integration.

What I mean about error interpolation is that whatever the graphics engine comes up with second-order linearly variable velocity equations get blended in (interpolated according to timestamp) with rk4 physical keyframes. Extrapolation can be made to work like this just as well: whenever the graphics engine has to extrapolate past a physics keyframe, it records the extrapolated state so when it goes back to interpolating it can gracefully transition from the extrapolated position back to the simulated state, hiding integration errors over time.

In fact, what you're proposing is what I've been saying all along, only you fail to realize that whenever the physical state's timestamp changes, discrepancy between rk4 simulation and your second-order integration will create discontinuities across physical frames, which is bad - not sure how perceptible it will be, since it depends how rk4 error compares to said second-order integration (integration that is imprecise, btw, you can use quadratic polynomials with little computational overhead: newV = V + dt * A, newP = P + V * dt + 0.5 * A * A - easily vectorizable and minimal data dependencies)
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

klauss wrote:@errol: so much I won't read (seems like a conversation between you and chuck anyways)...
Please jump in. We are having a fundamental misunderstanding about how the physics and interpolator work together.

About the <axis, angle> representation, how would you interpolate between two different such vectors?

Because <axis, angle> has no speed implied, and so interpolating naively would produce artifacts IMHO. Suppose a unit that has been spinning a long time in a certain way, so you have angle = <huge>*360° + <small>. Then you have the new rotation vector computed out of a collision that has angle = <small> - naive interpolation would "roll back" the huge rotation, causing a severe interpolation artifact.

What we care most is the rotational speed, rather... so we can interpolate speeds, with no discernible artifact. I was hoping (though I haven't proven it) that my 5-tuple would achieve that precisely.
Ok, perhaps I should clarify what "I" mean when i say an <angle,axis> object. Quaternions can represent two things. They can represent an orientation or they can represent a rotation. If you think about it, an orientation is just an object rotated in some way with respect to a default orientation. The default orientatation is (1,0,0,0). Any quaternion q that acts on (1,0,0,0) always returns itself. q=q*(1,0,0,0). So the orientation quaternion q is just the rotation quaternion q acting on the default orientation. A rotation quaternion is when an orientation quaternion acts on something other than the default orienation. It is a rotation with respect to something other than the default orientation. For example q2=r*q1 are all quaternions. Now a quaternion is constructed such that the first term is cos(angle/2). The next 3 terms are the unit vector of the axis of rotation times sin(angle/2) However, for angles >360 degrees, the quaternion does not tell you how many rotation were made. Thus, it is useful to simply give an (angle,axis) object from which you can construct a quaternion when you need to. But, the angle,axis object can be either a rotation or an orientation just like a quaternion.


So, this is how "I" would do it. 1st option: linear interpolation. rk4 needs to give you the rotation <angle,axis>that gets you from the current state to the future state. If you just have the current state <angle,axis> and future state <angle,axis>, you could, in principle extract the rotation that got you there, but it is computationally expensive and complicated. It is more efficient to just carry it around in the state. So, the rotation <angle,axis> object is given. The value of angle is just the average angular velocity (av) times the time step(ts). So, I write interp_angle=(angle/ts)*t where t is the elapsed time since current state. You can see that when t=ts, we are at the future state and the value of angle is what it should be. Then, I construct the interp_rot_quaternion=quaternion(cos(interp_angle/2),axis*sin(interp_angle/x)). I then say interp_quaternion=interp_rot_quaternion*current_state_orientation_quaternion. This gives you the orientation at a time in between the current state and future state. And, you are guaranteed to end up where you are supposed to at t=ts. I see no reason for artifacts to occur. I don't understand the scenario you are proposing. I think you are imagining that the <angle,axis> object will just keep building up the value of angle as it spins, but the object's orientation is NOT being stored in an <angle,axis> object. It is being stored in a quaternion. The only use for the <angle, axis> object is to avoid losing information about the number of revolutions in a single time step for purposes of interpolation. If an object is spinning at 3.2 revolutions per physics step, then when we go to interpolate for frames, if we just have the orientation before, and the orientation after and a time step, it will appear as if the thing is spinning at only .2 revolution per physics step. So, we make sure to give the total angle that the thing spun during that physics step so we can properly calculate the interpolated orientations. For example, you may have 10 interpolated steps in one physics step. If we remember the angle of the rotation during this physics step was 3.2 revs, then we can easily see that each of our 10 interpolated steps should involve a rotation of .32 revolutions about the axis and we can then construct the right rotation quaternion that will get us to the new orientation. If we didn't keep track of the angle, and just kept track of the orientations, each of the interpolated steps would only be .02 revolutions. Both scenarios get you to the right orientation in the end, but one misrepresents the angular velocity of the object. Once the physics step is done, all we have is an orientation that has no memory of how many times it has been spun. However, even if it did, it wouldn't cause any problems. So you have angle =360*n+10 degrees. The orientation is 10 degrees off of the default orientation. Then it experiences a -10 degree rotation. It is now at the default orientation. Where is the artifact?

2nd option is a cubic interpolation. It is the same principle, but I also require that av be continuous. It is 3 times as computationally expensive, but it is about 30 times more complex to code. I am having enough trouble conveying cubic interpolation for linear motion, and I would have to sit down and work out how to do it for rotations, so I won't even try to explain it here until after I have done it myself. However, I know it is doable. It is just a matter of dealing with the quaternions.
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

Ehm... I thought the idea was to simulate one frame ahead in physics, and show stuff with a small lag (one physics frame lag) so we're always interpolating rather than extrapolating.
YES! This is what I am trying to say.
to hide artifacts because of discrepance between second-order integration in the interpolator and rk4 integration.
There would not be any discrepancies if he stops integrating in his interpolator. That is the job of the rk4 algorithm. The only job of the interpolator is to go from where it currently is to where the rk4 algorithm told it to go in a continuous smooth way. The current algorithm fails at both. A simple linear interpolation succeeds at being continuous and a cubic interpolation is additionally smooth.
whenever the graphics engine has to extrapolate past a physics keyframe, it records the extrapolated state so when it goes back to interpolating it can gracefully transition from the extrapolated position back to the simulated state, hiding integration errors over time.
It shouldn't ever have to interpolate past a physics keyframe. There should always be a computed future state towards which the object is traveling. How it gets there is determined by interpolation
In fact, what you're proposing is what I've been saying all along, only you fail to realize that whenever the physical state's timestamp changes, discrepancy between rk4 simulation and your second-order integration will create discontinuities across physical frames, which is bad - not sure how perceptible it will be, since it depends how rk4 error compares to said second-order integration (integration that is imprecise, btw, you can use quadratic polynomials with little computational overhead: newV = V + dt * A, newP = P + V * dt + 0.5 * A * A - easily vectorizable and minimal data dependencies)
YES! ^ is the problem. It is not the inaccuracy of your extrapolation that is the problem. The problem is that it doesn't match what is being calculated by the physics.
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

@errol: oh... now I get it. I indeed was thinking you were going to use <axis,angle> for state as well, hence the misunderstanding.

It all seems to add up rather nicely the way you describe it. :D
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

klauss wrote:
chuck_starchaser wrote:...But Errol thinks it is too incorrect.
I'm aware of its incorrectness; but as I mentioned in my last post, I think interpolation, even Hermite, would be just as incorrect, since what we really do is extrapolate, in most cases, anyways.
If we do this, then we won't have two quaternions to interpolate between; so, problem solved :mrgreen:
Ehm... I thought the idea was to simulate one frame ahead in physics, and show stuff with a small lag (one physics frame lag) so we're always interpolating rather than extrapolating.
But you're only thinking about the fastest simulation bin; or are you saying that a unit simulated every 8 seconds would compute frames 8 seconds ahead?
(This might be okay, as long as collisions are computed as often and as interpolated as necessary... just asking.)
In any case, true second-order inter/extrapolation (interpolate speed and apply linearly variable speed equations for inter/extrapolation) isn't terribly incorrect IMHO, and error (of inter/extrapolation vs simulation) can be interpolated (linearly) just as well, to hide artifacts because of discrepance between second-order integration in the interpolator and rk4 integration.

What I mean about error interpolation is that whatever the graphics engine comes up with second-order linearly variable velocity equations get blended in (interpolated according to timestamp) with rk4 physical keyframes. Extrapolation can be made to work like this just as well: whenever the graphics engine has to extrapolate past a physics keyframe, it records the extrapolated state so when it goes back to interpolating it can gracefully transition from the extrapolated position back to the simulated state, hiding integration errors over time.

In fact, what you're proposing is what I've been saying all along, only you fail to realize that whenever the physical state's timestamp changes, discrepancy between rk4 simulation and your second-order integration will create discontinuities across physical frames, which is bad - not sure how perceptible it will be, since it depends how rk4 error compares to said second-order integration (integration that is imprecise,
I think the discrepancies will be negligible. What we all seem to be forgetting is that the reason we want the extra precision rk4 affords us is NOT (immediately) perceptual, but rather for cumulative error reasons. The step to step differences between RK4 and a much simpler simulation should be far beyond our ability to perceive. And if it's not; our stepping is wrong, in the first place.
IOW, I agree about the incorrectnes; but I doubt its import.
Let us not forget that presently the way simulation steps are assigned is horrendously buggy. Sometimes ships in plain sight are obviously using very long steps, even when there's hardly any ships around. I'm assuming we will have a functional way of prioritizing simulation frequencies. And if we do so, I doubt even the inclusion of acceleration will have a very noticeable effect.
But if what we end up doing is always compute a frame ahead (1 second ahead for the 1 second bin; 8-seconds ahead for the 8-second bin), and use interpolation, then I rest my case. But I'm not sure how we can go as far as 1 second... It means that AI commands may be delayed.
Well, on the other hand, we'd probably have AI commands cause temporary accelerations of physics simulation, I guess.
btw, you can use quadratic polynomials with little computational overhead: newV = V + dt * A, newP = P + V * dt + 0.5 * A * A - easily vectorizable and minimal data dependencies)
Cool!
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

chuck_starchaser wrote:
klauss wrote:
chuck_starchaser wrote:...But Errol thinks it is too incorrect.
I'm aware of its incorrectness; but as I mentioned in my last post, I think interpolation, even Hermite, would be just as incorrect, since what we really do is extrapolate, in most cases, anyways.
If we do this, then we won't have two quaternions to interpolate between; so, problem solved :mrgreen:
Ehm... I thought the idea was to simulate one frame ahead in physics, and show stuff with a small lag (one physics frame lag) so we're always interpolating rather than extrapolating.
But you're only thinking about the fastest simulation bin; or are you saying that a unit simulated every 8 seconds would compute frames 8 seconds ahead?
(This might be okay, as long as collisions are computed as often and as interpolated as necessary... just asking.)
I've been thinking about the matter, and I'm thinking exactly what you're thinking ;)
chuck_starchaser wrote:But if what we end up doing is always compute a frame ahead (1 second ahead for the 1 second bin; 8-seconds ahead for the 8-second bin), and use interpolation, then I rest my case. But I'm not sure how we can go as far as 1 second... It means that AI commands may be delayed.
Game-changing AI commands would have to re-schedule the unit. Store the interpolated state as current physical state, and then re-schedule. I've been down that road before.
chuck_starchaser wrote:I think the discrepancies will be negligible. What we all seem to be forgetting is that the reason we want the extra precision rk4 affords us is NOT (immediately) perceptual, but rather for cumulative error reasons. The step to step differences between RK4 and a much simpler simulation should be far beyond our ability to perceive. And if it's not; our stepping is wrong, in the first place.
IOW, I agree about the incorrectnes; but I doubt its import.
Let us not forget that presently the way simulation steps are assigned is horrendously buggy. Sometimes ships in plain sight are obviously using very long steps, even when there's hardly any ships around. I'm assuming we will have a functional way of prioritizing simulation frequencies. And if we do so, I doubt even the inclusion of acceleration will have a very noticeable effect.
You are right, yet you are not.
You are right, I doubt it's that important.
But you are wrong, in that solving this particular problem will be of immense importance for multiplayer - when an update fails to reach the client, such an error resolution and hiding framework will pay off big time (I know how "high" multiplayer is on your list ;) - it is pretty much equally high on mine - but it won't hurt thinking about the subject when we're in such a related matter as physics and state inter/extrapolation).
You are also wrong in that very far units may still be visible to the player - as is the case with the asteroids in an asteroid field (granted, those would probably be close to static, but... well...) or the particles in an accretion disc around a singularity. With long time steps (necessary for asteroids since they're too many to simulate at a higher priority, for instance) come bigger errors, and having error concealment may pay off.

So it's harmless brainstorming.
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

klauss wrote:I've been thinking about the matter, and I'm thinking exactly what you're thinking ;)
Alright; so at each frequency bin we are computing a full delta_t ahead of time, then; and both collision and AI events will cause rescheduling.
So, we're back to interpolating two data points, with the benefit of knowing position, speed and acceleration for each; not sure what's a good algorithm to use.
So, I'm going to modify the physics to,
  • rename "current" to "next"
  • bring back the availability of "previous"
  • and, for clarity, I will NOT have a state named "current", since there is no such thing, so it would be just confusing us; I will have a third state called "scratch", as in scratch-pad, instead, which has obsolete data (prior to "previous") until it is used to compute the 'next-next' (which won't become next until the pointers to the states are rotated, at the end of the physics run on the current frequency bin). The rotation will be done by rotation of three pointers, to avoid copying.
You are right, yet you are not.
You are right, I doubt it's that important.
But you are wrong, in that solving this particular problem will be of immense importance for multiplayer - when an update fails to reach the client, such an error resolution and hiding framework will pay off big time (I know how "high" multiplayer is on your list ;) - it is pretty much equally high on mine - but it won't hurt thinking about the subject when we're in such a related matter as physics and state inter/extrapolation).
Touche. I hadn't thought of that. Crazy latencies, lost packets...
Okay, so we absolutely need interpolation.
You are also wrong in that very far units may still be visible to the player - as is the case with the asteroids in an asteroid field (granted, those would probably be close to static,
Exactly. I don't thin cemetary furniture would look as static. Unless some mod wants asteroids acting crazy; but then they'd have their own controller, rather than expect the physics engine to support craziness.
but... well...) or the particles in an accretion disc around a singularity. With long time steps (necessary for asteroids since they're too many to simulate at a higher priority, for instance) come bigger errors, and having error concealment may pay off.
Makes sense.
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

But you're only thinking about the fastest simulation bin; or are you saying that a unit simulated every 8 seconds would compute frames 8 seconds ahead?
(This might be okay, as long as collisions are computed as often and as interpolated as necessary... just asking.)
That is the idea. When a collision occurs, you would need to stop both objects along their interpolated path and compute the effects of the collision.
So, we're back to interpolating two data points, with the benefit of knowing position, speed and acceleration for each; not sure what's a good algorithm to use.
You have 3 options. I tried to explain them in this thread, but I get the feeling everybody skimmed over the math :oops: , so here is the summary:
x0, v0, and a0 are position, velocity, and acceleration at the "previous" state(I heartily approve of the renaming. I was thinking exactly the same thing), and x1, v1, and a1 are position, velocity, and acceleration in the previous state. dt=t1-t0 where t1 is time at the "next" state and t0 is the time associated with the "previous" state. The time elapsed since t0 is given by t

Option 1: Linear interpolation
x=x0+t*(-x0 + x1)/dt
v will appear to be constant over the interval
a will appear to be zero over the interval
Note that when t=dt, the object should be at the "next" position, and if you plug it in in the above equation, you will see that this is the case.

Option 2: Cubic interpolation
C1=(-2*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
C2=(6*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
x=x0+v0*t+C1*t^2+C2*t^3
v=v0+2*C1*t+3*C2*t^2
a=2*C1+6*C2*t
Note that when t=dt, x=x1 and v=v1, but a!=a1

Option 3: 5th order interpolation
C1=(-3*(3*a0*Power(dt,2) - a1*Power(dt,2) + 12*dt*v0 + 8*dt*v1 + 20*x0 -
20*x1))/Power(dt,3)
C2=(12*(3*a0*Power(dt,2) - 2*a1*Power(dt,2) + 16*dt*v0 + 14*dt*v1 + 30*x0 -
30*x1))/Power(dt,4)
C3=(-60*(a0*Power(dt,2) - a1*Power(dt,2) + 6*dt*v0 + 6*dt*v1 + 12*x0 -
12*x1))/Power(dt,5)
x=x0+v0*t+.5*a0*t^2+C1*t^3+C2*t^4+C3*t^5
v=v0+a0*t+3*C1*t^2+4*C2*t^3+5*C3*t^4
a=a0+6*C1*t+12*C2*t^2+20*C2*t^3

Note that when t=dt, x=x1,v=v1, and a=a1

You CAN get quadratic and quartic solutions, but you have to give up one of your endpoint matching requirements from the cubic and 5th order interpolations respectively.

EDIT: For what its worth, I strongly recommend having all three interpolation routines available, but I think that cubic interpolation gives you a good mix of smoothness with small computational expense. Remember the constants C1, C2, and C3 need only be calculated once per physics step. Things with small physics step will probably be better off with linear interpolation since computing C1 and C2 for the cubic could get computationally expensive with a small physics step. However, things with big physics steps will probably want a better interpolation to keep their trajectory looking smooth.

not sure how perceptible it will be, since it depends how rk4 error compares to said second-order integration (integration that is imprecise, btw, you can use quadratic polynomials with little computational overhead: newV = V + dt * A, newP = P + V * dt + 0.5 * A * A - easily vectorizable and minimal data dependencies)
If you do a direct run-time comparison, rk4 will, in most cases out perform the quadratic polynomial extrapolation (what I call midpoint method).

In other words, even though rk4 is about 5 times as computationally expensive, it gets much better precision meaning that you can take much longer time steps and get the same accuracy (or better). This is ESPECIALLY noticeable for things with large time steps where the higher order terms that midpoint method neglects become more important. If you don't want planets to eat up computation, but you want them to still remain in nice stable orbits, you can have rk4 do 1 minute time steps, or you can have midpoint method do 6 second time steps. Their accuracy will be comparable and the rk4 will take less cpu cycles per second.

I could write a program to prove this by writing two adaptive time step routines that require error<10^-6 and see which finishes first. Rk4 will take bigger steps and still stay within the error range. Midpoint method will take smaller step to stay within the error range. I would test them on planetary orbits first. However, I don't think that is necessary, because it is well documented that, given equivalent error requirements, rk4 will outperform midpoint method so long as the acceleration vector varies non-linearly over the time interval.

Anyway, both algorithms are in EJphysics.h. I would use rk4 for most objects and just make their physics time step bigger if you don't need that much accuracy or their acceleration vectors are slowly varying. I have not tested it to be sure, but if my ballpark calculations are correct, rk4 should be able to do track the trajectory of a planet like earth with a 1 day time step without any noticeable decay over at least 100 orbits (years).

Hmmm, I just thought of something. It may be a good idea to have the interpolated position of an object handy in addition to its "previous" and "next" positions. The reason is as follows. Imagine a planet going on a 1-day physics time step and using cubic interpolation. It is near the end of the planet's physics step when a player flies nearby where the planet used to be nearly 1 day ago. Earth's orbital speed is 30 km/s. Even 1 hour ago, that is 108000 km away from where the planet currently sits. However, when the player's ship calculates the force on it due to nearby gravitational objects, it needs to make sure it grabs the most recently interpolated position of the planet rather than its previous or next state. Otherwise, the player's ship will experience gravity due to an object that is not there anymore.
@errol: oh... now I get it. I indeed was thinking you were going to use <axis,angle> for state as well, hence the misunderstanding.
I hadn't even thought of the problem until chuck mentioned it to me, so I was just keeping the quaternions around. But that is how I "think" chuck intended to handle it as well. We don't exactly speak the same language all the time. I speak physicsease and he speaks computerease
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

Errol_Summerlin wrote:
But you're only thinking about the fastest simulation bin; or are you saying that a unit simulated every 8 seconds would compute frames 8 seconds ahead?
(This might be okay, as long as collisions are computed as often and as interpolated as necessary... just asking.)
That is the idea. When a collision occurs, you would need to stop both objects along their interpolated path and compute the effects of the collision.
Gottcha.
Same goes for any changes by the AI: If the AI was turning to the right, and it now decides to nose-up instead, it will have to notify the physics to please reschedule the unit as to simulate the next frame ASAP.
And I haven't delved into the multiplayer code; but probably the server is trying to predict players' positions, to compensate for internet lag, and to reduce bandwidth; and probably the same applies when user input changes (--i.e.: probably clients only send packets to the server about user input, and only when it changes; and the server would simulate their trajectories, and should switch to fast stepping after it receives a packet for a unit, then slow down the simulation steps; but I'm just guessing).
So, we're back to interpolating two data points, with the benefit of knowing position, speed and acceleration for each; not sure what's a good algorithm to use.
You have 3 options. I tried to explain them in this thread, but I get the feeling everybody skimmed over the math :oops: , so here is the summary:
x0, v0, and a0 are position, velocity, and acceleration at the "previous" state(I heartily approve of the renaming. I was thinking exactly the same thing), and x1, v1, and a1 are position, velocity, and acceleration in the previous state. dt=t1-t0 where t1 is time at the "next" state and t0 is the time associated with the "previous" state. The time elapsed since t0 is given by t

Option 1: Linear interpolation
x=x0+t*(-x0 + x1)/dt
v will appear to be constant over the interval
a will appear to be zero over the interval
Note that when t=dt, the object should be at the "next" position, and if you plug it in in the above equation, you will see that this is the case.

Option 2: Cubic interpolation
C1=(-2*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
C2=(6*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
x=x0+v0*t+C1*t^2+C2*t^3
v=v0+2*C1*t+3*C2*t^2
a=2*C1+6*C2*t
Note that when t=dt, x=x1 and v=v1, but a!=a1

Option 3: 5th order interpolation
C1=(-3*(3*a0*Power(dt,2) - a1*Power(dt,2) + 12*dt*v0 + 8*dt*v1 + 20*x0 -
20*x1))/Power(dt,3)
C2=(12*(3*a0*Power(dt,2) - 2*a1*Power(dt,2) + 16*dt*v0 + 14*dt*v1 + 30*x0 -
30*x1))/Power(dt,4)
C3=(-60*(a0*Power(dt,2) - a1*Power(dt,2) + 6*dt*v0 + 6*dt*v1 + 12*x0 -
12*x1))/Power(dt,5)
x=x0+v0*t+.5*a0*t^2+C1*t^3+C2*t^4+C3*t^5
v=v0+a0*t+3*C1*t^2+4*C2*t^3+5*C3*t^4
a=a0+6*C1*t+12*C2*t^2+20*C2*t^3

Note that when t=dt, x=x1,v=v1, and a=a1

You CAN get quadratic and quartic solutions, but you have to give up one of your endpoint matching requirements from the cubic and 5th order interpolations respectively.

EDIT: For what its worth, I strongly recommend having all three interpolation routines available, but I think that cubic interpolation gives you a good mix of smoothness with small computational expense. Remember the constants C1, C2, and C3 need only be calculated once per physics step. Things with small physics step will probably be better off with linear interpolation since computing C1 and C2 for the cubic could get computationally expensive with a small physics step. However, things with big physics steps will probably want a better interpolation to keep their trajectory looking smooth.
Great stuff!

Yeah, at first glance, the cubic stepping looks the most tempting.
By the way, for interpolation purposes, we only need to interpolate position, so it boils down to,

C1=(-2*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
C2=(6*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
x=x0+v0*t+C1*t^2+C2*t^3

<finito>
Well, we'll need to interpolate velocity when a collision occurs.

And as you say, C1 and C2 are per-physics-step, so the only interpolation code is,

x=x0+v0*t+C1*t^2+C2*t^3



EUREKA!

We don't need to keep around the previous physics frame, then; what we need is a new struct to which we only copy X0 and V0, and compute C1 and C2.
In fact, we may not even need the two frames we got now, if we use a different structure for output... (I need to check the rk4 code to see if there would be any aliasing problems reading from, and writing to, the same state) ...though we will need two ouput structures and a flipper flag, for concurrent update/use from separate threads (namely the graphics thread, as I think collisions will be in the same thread as physics). ((On the other hand, if we unify interpolations for collisions AND graphics, then we have no concurrency requirements for this structure. The way this would work is, we get a display list from the graphics thread, periodically updated. Then, whenever we compute interpolations for collisions, we add to our list of units to interpolate the list from the graphics thread, IF the time to the next graphics frame is near. If we do this, then it is the list of final, interpolated positions that would have concurrency requirements.))

..........
Anyway, both algorithms are in EJphysics.h. I would use rk4 for most objects and just make their physics time step bigger if you don't need that much accuracy or their acceleration vectors are slowly varying. I have not tested it to be sure, but if my ballpark calculations are correct, rk4 should be able to do track the trajectory of a planet like earth with a 1 day time step without any noticeable decay over at least 100 orbits (years).
I'm beginning to think I'd better change my scheduler code to have 16 simulation frequency bins, then, instead of 8. With 8 bins I only have steps from 1/16th of a second to 8 seconds. With another 8 we could have up to a little over one hour per step.
I don't think having steps longer than an hour would make any difference in terms of computational load, but I could be wrong... If we do want to have younger stars with proto-planetary disks, with all planetessimals simulated, as I think Klauss hinted at, then for sure; but then I'm not sure where the initial conditions for such a vast and chaotic set of objects would come from.
Hmmm, I just thought of something. It may be a good idea to have the interpolated position of an object handy in addition to its "previous" and "next" positions. The reason is as follows. Imagine a planet going on a 1-day physics time step and using cubic interpolation. It is near the end of the planet's physics step when a player flies nearby where the planet used to be nearly 1 day ago. Earth's orbital speed is 30 km/s. Even 1 hour ago, that is 108000 km away from where the planet currently sits. However, when the player's ship calculates the force on it due to nearby gravitational objects, it needs to make sure it grabs the most recently interpolated position of the planet rather than its previous or next state. Otherwise, the player's ship will experience gravity due to an object that is not there anymore.
WOW. power( Indeed!, 77.7 )...

(Another factor probably we all tend to forget is that in vegastrike, all speed vectors are "absolute" (relative to the local star), so even for two units of low speeds relative to each other, simulation steps may involve very long travel distances for both. We'll eventually going to have to bite the bullet and implement local, moving frames of reference; but that's going to be a real code challenge.)

Ouch...
You know what this implies?
I'm thinking it implies interpolations for *everything* all the time.
You see, in my mind, we were only going to interpolate units we needed to interpolate, at any given time.
That would have included,
  • a) Units that need interpolation for collisions, which is a small subset, at each basic step, because the first layer or two of the collisions filter are predictive eliminators: Namely, if something is too far from any other units, it says "I won't be colliding with anything for at least 6 seconds, so check with me again then."
  • b) Units that need interpolation for display, which is also a small subset because it ignores firstly all units outside the visual frustrum, and secondly it elliminates possibly all units whose | velocity - player_velocity | * physics_dt / distance < ~ 0.5 * sin( 1 pixel width view angle ).
Well, no big deal; what it implies is adding all gravitational attractors (units of mass > some_mass_TBD), which will probably add maybe a a couple of dozen units to the interpolation list, unconditionally; but that's nothing to write home about. What IS worth writing home about is the fact that physics will have to check this "current" interpolation data, which I'm not sure if it's synchronous or asynchronous...
Well, I think this settles the question, and interpolations will HAVE TO be in the same thread as physics and collisions.

I'm a bit worried that, in time, we're going to keep discovering other cases needing interpolations...
Once we implement more sophisticated sensors, for instance, a small ship 100,000 km away would, according to physics and graphics requirements, only need one physics step every 8-seconds, and no interpolations at all, ever; BUT, if we have a very sensitive heat sensor, that same ship, when it turns such that its engines are pointed towards us, should cause a momentary blip on our screens. But we wouldn't know its rotational position without interpolations.
Another example: Ships behind us may not need graphical interpolations as far as drawing them on the screen, but may need them as far as their blips on our sensors are concerned.

@errol: oh... now I get it. I indeed was thinking you were going to use <axis,angle> for state as well, hence the misunderstanding.
I hadn't even thought of the problem until chuck mentioned it to me, so I was just keeping the quaternions around. But that is how I "think" chuck intended to handle it as well. We don't exactly speak the same language all the time. I speak physicsease and he speaks computerease
I'm not sure what I was thinking, either; or IF I was thinking; can't remember.
I think I'm thinking now, though: Angle and axis is what we need for rotational position, IF we want to interpolate. But for
rotational speed and acceleration we could continue to use quaternions, I suppose; though I'm not sure I understand how the math for expressing rotational speed with quaternions works... Do we have the same problem with having a rotational speed limit beyond which it would have to become a non-unit quaternion or something?
Furthermore, if we're going to do cubic interpolation, and for that we only need rotational X0, V0, C1 and C2, then might as well add these to the output struct; so,

Code: Select all

struct sim_output
{
    vec3<double> X0, V0, C1, C2;
    quat<double> rX0, rV0, rC1, rC2;
};
So, now the problem boils down to making sure that the formula,
rx=rx0+rv0*t+rC1*t^2+rC2*t^3
is implemented in such a way that it works when interpolated angles >> 360 deg.

EDIT:
Well, I guess with quaternions it would go something like,
t/=dt;
rx=rx0+rv0^t+rC1^(t^2)+rC2^(t^3)



If that's not possible, then maybe we'd have to,

Code: Select all

struct sim_output
{
    vec3<double> X0, V0, C1, C2;
    angle_n_axis rX0, rV0, rC1, rC2;
};
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

chuck_starchaser wrote:EDIT:
Well, I guess with quaternions it would go something like,
t/=dt;
rx=rx0+rv0^t+rC1^(t^2)+rC2^(t^3)



If that's not possible, then maybe we'd have to,

Code: Select all

struct sim_output
{
    vec3<double> X0, V0, C1, C2;
    angle_n_axis rX0, rV0, rC1, rC2;
};
Uhhh... Never mind; I forgot that angular speed and acceleration were vectors.

Code: Select all

    vec3< double > acceleration;
    vec3< double > angular_acceleration;
    vec3< double > velocity;
    vec3< double > angular_velocity;
    vec3< double > position;
    quat< double > orientation; //orientation in 3D space represented as a quat<double>
    vec3< double > tensor_of_inertia;
    double mass;
    double inverse_mass;
    double dm_dt; //rate of mass gain (<0 for loss)
    double time; //current time
I'm almost done with the changes.
No problems reading from and writing to the same, single state; so what I have now is two output states, instead.
Currently it is declared as,

Code: Select all

    struct OutputState
    {
        vec3<double> X0, V0, C1, C2;
        quat<double> rX0, rV0, rC1, rC2;
    };
But that's wrong, because rX0 is a quaternion, while rV0 is a vector of rotation.
Then again, I'm not sure at all what rC1 and rC2 are...
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

Errol, there's problems with interpolation. Check the results first: I kept the test program like you had it, generally, with a printout every 1000 iterations, but fleshed it out to include two consecutive physical iterations, and three interpolations between them at t=0, t=dt/2 and t=dt. The last interpolation should equal the physics simulation that follows it. But they don't: At t=0 the interpolation is just right; at t=dt/2 it also seems quite right to me; but at t=dt it always goes way off (dt = 0.1, btw):

Code: Select all

iter=999:
x=0.8633, y=-0.504635, z=0, energy = -0.500014
interpolations:
x=0.8633, y=-0.504635, z=0, @t=0
x=0.886315, y=-0.460317, z=0, @t=0.05
x=11.8171, y=-5.40523, z=0, @t=0.1
iter=1000:
x=0.909369, y=-0.415924, z=0, energy = -0.500014

iter=1999:
x=0.494254, y=-0.869254, z=0, energy = -0.500028
interpolations:
x=0.494254, y=-0.869254, z=0, @t=0
x=0.536378, y=-0.842431, z=0, @t=0.05
x=7.51826, y=-10.5984, z=0, @t=0.1
iter=2000:
x=0.578572, y=-0.815563, z=0, energy = -0.500028

iter=2999:
x=-0.00364066, y=-0.99991, z=0, energy = -0.500042
interpolations:
x=-0.00364066, y=-0.99991, z=0, @t=0
x=0.0462452, y=-0.997596, z=0, @t=0.05
x=1.24997, y=-12.9336, z=0, @t=0.1
iter=3000:
x=0.0962143, y=-0.995277, z=0, energy = -0.500042

iter=3999:
x=-0.496933, y=-0.867661, z=0, energy = -0.500056
interpolations:
x=-0.496933, y=-0.867661, z=0, @t=0
x=-0.452411, y=-0.890283, z=0, @t=0.05
x=-5.29984, y=-11.8635, z=0, @t=0.1
iter=4000:
x=-0.407814, y=-0.912943, z=0, energy = -0.500056

iter=4999:
x=-0.858435, y=-0.512652, z=0, energy = -0.500069
interpolations:
x=-0.858435, y=-0.512652, z=0, @t=0
x=-0.830717, y=-0.554196, z=0, @t=0.05
x=-10.4346, y=-7.74225, z=0, @t=0.1
iter=5000:
x=-0.802954, y=-0.595808, z=0, energy = -0.500069

iter=5999:
x=-0.999377, y=-0.0301945, z=0, energy = -0.500083
interpolations:
x=-0.999377, y=-0.0301945, z=0, @t=0
x=-0.995376, y=-0.0799756, z=0, @t=0.05
x=-12.8828, y=-1.68693, z=0, @t=0.1
iter=6000:
x=-0.991367, y=-0.12984, z=0, energy = -0.500083

iter=6999:
x=-0.889678, y=0.456162, z=0, energy = -0.500097
interpolations:
x=-0.889678, y=0.456162, z=0, @t=0
x=-0.910214, y=0.410637, z=0, @t=0.05
x=-12.0954, y=4.74395, z=0, @t=0.1
iter=7000:
x=-0.930785, y=0.365036, z=0, energy = -0.500097

iter=7999:
x=-0.562211, y=0.826725, z=0, energy = -0.500111
interpolations:
x=-0.562211, y=0.826725, z=0, @t=0
x=-0.602054, y=0.796611, z=0, @t=0.05
x=-8.34203, y=9.96016, z=0, @t=0.1
iter=8000:
x=-0.641963, y=0.766446, z=0, energy = -0.500111

iter=8999:
x=-0.101368, y=0.994598, z=0, energy = -0.500125
interpolations:
x=-0.101368, y=0.994598, z=0, @t=0
x=-0.150739, y=0.987056, z=0, @t=0.05
x=-2.60117, y=12.7287, z=0, @t=0.1
iter=9000:
x=-0.200192, y=0.979502, z=0, energy = -0.500125
Here's the relevant code in rk4Step():

Code: Select all

    ......................
    helper( *this, 0.5*dt, d2 );
    helper( d2, 0.5*dt, d3 );
    helper( d3, dt, d4 );
    //C1=(-2*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
    //C2=(6*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
    inv_dt = 1.0/dt;
    temp = inv_dt*inv_dt;
    ostate.X0 = position;
    ostate.V0 = velocity;
    ostate.C1 = position*3.0 + velocity*(dt*2.0);
    ostate.C2 = position*2.0 + velocity*dt;
    position  += ( ( 1.0/6.0 )*dt*( velocity+2.0*( d2.velocity+d3.velocity )+d4.velocity ) );
    velocity  += ( (1.0/6.0)*dt*( acceleration+2.0*( d2.acceleration+d3.acceleration )+d4.acceleration ) );
    ostate.C1 += position*3.0 + velocity*dt;
    ostate.C2 += position*2.0 + velocity*dt;
    ostate.C1 *= (-2.0*temp);
    ostate.C2 *= (6.0*temp*inv_dt);
    .............................................
Nameley, after the calls of helper(), the two comment lines are the formulas for C1 and C2, copied from your post. Then,
Next two lines: Inverse powers of dt.
Next two lines: Copying of position and velocity.
Next line: x0 and v0 terms for C1
Next line: x0 and v0 terms for C2
Next two lines: position and velocity are updated, becoming x1 and v1
Next line: x1 and v1 terms for C1
Next line: x1 and v1 terms for C2
Next two lines: finalizing C1 and C2

I can't see any errors in the code; but maybe I'm blind.
And here's the interpolation function:

Code: Select all

    vec3< double > get_interpolated_position( double t )
    {
        OutputState state = get_output_state();
        return vec3< double >( state.X0 + state.V0*t + state.C1*(t*t) + state.C2*(t*t*t) );
    }
Which comes from
you wrote:x=x0+v0*t+C1*t^2+C2*t^3
The relevant code in the testing routine:

Code: Select all

    for( size_t i = 0; i < 10000; ++i )
    {
        vec3< double > ipos; //interpolated position
        object.rk4Step( dt );
        //object.MidpointStep(object.get_state(),dt);
        object.AdvanceState();
        if( (i+1)%1000 == 0 )
        {
            std::cout
                << "\niter=" << i << ":\n"
                << "x=" << object.position.get_x()
                << ", y=" << object.position.get_y()
                << ", z=" << object.position.get_z()
                << ", ";
            ::fflush( NULL );
            energy = object.velocity.lengthSquared()/2-1
                     /object.position.length();
            std::cout
                << "energy = " << energy<< "\n";
        }
        if( i%1000 == 0 )
        {
            std::cout << "interpolations:\n";
            ipos = object.get_interpolated_position( 0.0 );
            std::cout
                << "x=" << ipos.get_x()
                << ", y=" << ipos.get_y()
                << ", z=" << ipos.get_z()
                << ", @t=" << 0.0 << "\n";
            ::fflush( NULL );
            ipos = object.get_interpolated_position( 0.5*dt );
            std::cout
                << "x=" << ipos.get_x()
                << ", y=" << ipos.get_y()
                << ", z=" << ipos.get_z()
                << ", @t=" << 0.5*dt << "\n";
            ::fflush( NULL );
            ipos = object.get_interpolated_position( dt );
            std::cout
                << "x=" << ipos.get_x()
                << ", y=" << ipos.get_y()
                << ", z=" << ipos.get_z()
                << ", @t=" << dt << "\n";
            ::fflush( NULL );
            std::cout
                << "iter=" << i << ":\n"
                << "x=" << object.position.get_x()
                << ", y=" << object.position.get_y()
                << ", z=" << object.position.get_z()
                << ", ";
            ::fflush( NULL );
            energy = object.velocity.lengthSquared()/2-1
                     /object.position.length();
            std::cout
                << "energy = " << energy<< "\n";
        }
    }
I suspect some typo, either in your math or in my code.

By the way, the problem remains with rotational interpolations; I need your help: Can't compute rC1 and rC2, the way things stand, because rotational speed and acceleration are vectors, whereas orientation is a quaternion. Furthermore, this quaternion is multiplied by a temporary quaternion at each step, that is obtained by conversion from rotational speed. I think I need to convert everything to vectors of rotation and forget quaternions, but I don't know how exactly.

On another topic: We were discussing gyroscopes by PM, a while ago, and you seemed to say that gyroscopic effects should result from the physics modeling as it stands. I have a nagging suspicion that you're wrong: Right now, RigidBody has no way to inject torque. It is a passive recepient of torque. When you try to turn the axis of a gyroscope, it resists with like a torque of its own making. Of course this is not so, I realize. But in any case, I don't know how the appearance of such a torque could emerge from RigidBody.
Then again, I don't know how an application of torque to an axis of rotation would be expressed.
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

LOL. It was my mistake; I forgot the negative signs in the x1 terms.
Fixed code in rk4Step():

Code: Select all

    helper( *this, 0.5*dt, d2 );
    helper( d2, 0.5*dt, d3 );
    helper( d3, dt, d4 );
    //C1=(-2*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
    //C2=(6*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
    inv_dt = 1.0/dt;
    temp = inv_dt*inv_dt;
    ostate.X0 = position;
    ostate.V0 = velocity;
    ostate.C1 = position*3.0 + velocity*(dt*2.0);
    ostate.C2 = position*2.0 + velocity*dt;
    position  += ( ( 1.0/6.0 )*dt*( velocity+2.0*( d2.velocity+d3.velocity )+d4.velocity ) );
    velocity  += ( (1.0/6.0)*dt*( acceleration+2.0*( d2.acceleration+d3.acceleration )+d4.acceleration ) );
    ostate.C1 -= position*3.0; //<--fixed
    ostate.C2 -= position*2.0; //<--fixed
    ostate.C1 += velocity*dt; //<--fixed
    ostate.C2 += velocity*dt; //<--fixed
    ostate.C1 *= (-2.0*temp);
    ostate.C2 *= (6.0*temp*inv_dt);
Now there are still discrepancies between interpoation at dt and the next physics frame, but rather minor:

Code: Select all

iter=999:
x=0.8633, y=-0.504635, z=0, energy = -0.500014
interpolations:
x=0.8633, y=-0.504635, z=0, @t=0
x=0.886315, y=-0.460317, z=0, @t=0.05
x=0.904665, y=-0.414138, z=0, @t=0.1
iter=1000:
x=0.909369, y=-0.415924, z=0, energy = -0.500014

iter=1999:
x=0.494254, y=-0.869254, z=0, energy = -0.500028
interpolations:
x=0.494254, y=-0.869254, z=0, @t=0
x=0.536378, y=-0.842431, z=0, @t=0.05
x=0.575396, y=-0.81166, z=0, @t=0.1
iter=2000:
x=0.578572, y=-0.815563, z=0, energy = -0.500028

iter=2999:
x=-0.00364066, y=-0.99991, z=0, energy = -0.500042
interpolations:
x=-0.00364066, y=-0.99991, z=0, @t=0
x=0.0462452, y=-0.997596, z=0, @t=0.05
x=0.0954005, y=-0.990311, z=0, @t=0.1
iter=3000:
x=0.0962143, y=-0.995277, z=0, energy = -0.500042

iter=3999:
x=-0.496933, y=-0.867661, z=0, energy = -0.500056
interpolations:
x=-0.496933, y=-0.867661, z=0, @t=0
x=-0.452411, y=-0.890283, z=0, @t=0.05
x=-0.406069, y=-0.908223, z=0, @t=0.1
iter=4000:
x=-0.407814, y=-0.912943, z=0, energy = -0.500056

iter=4999:
x=-0.858435, y=-0.512652, z=0, energy = -0.500069
interpolations:
x=-0.858435, y=-0.512652, z=0, @t=0
x=-0.830717, y=-0.554196, z=0, @t=0.05
x=-0.799119, y=-0.592549, z=0, @t=0.1
iter=5000:
x=-0.802954, y=-0.595808, z=0, energy = -0.500069

iter=5999:
x=-0.999377, y=-0.0301945, z=0, energy = -0.500083
interpolations:
x=-0.999377, y=-0.0301945, z=0, @t=0
x=-0.995376, y=-0.0799756, z=0, @t=0.05
x=-0.98643, y=-0.128858, z=0, @t=0.1
iter=6000:
x=-0.991367, y=-0.12984, z=0, energy = -0.500083
Though minor, the discrepancies (0.5%) seem larger than I would attribute to numerical precision issues, considering all our computations are done in double precision.
But anyhow, check my previous post for remaining problems --r.e.: quaternions vs vectors of rotation, etc.

EDIT:
Committing for now: r25.
EDIT2:
And here's the updated files compressed:
http://wcjunction.com/temp_images/RK4.tar.gz
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

Hey... chuck... if you have the test program, could you turn it into a benchmark?

And measure the frame rate with N simulated items with various time steps, interpolating them ALL at top speed?

That could give us valuable performance information, so we can decide the best architecture. Perhaps interpolation isn't too heavy and we can just interpolate all units in every collision/physics frame, and have the graphics thread do mere linear interpolation in between.
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

klauss wrote:Hey... chuck... if you have the test program, could you turn it into a benchmark?

And measure the frame rate with N simulated items with various time steps, interpolating them ALL at top speed?
This test program that Errol wrote simulates a particle going around an attractor; that's all. It simulates 10,000 steps, and it's ridiculously instantaneous. You press Enter and the results are there. Zero time. :D
Try it. Just uncompress the file, and type ./EJphysics
But I wouldn't know how to integrate it into VS rapidly; I need to do the Unit refactoring thing, --to separate all physics and collision concerns.
That could give us valuable performance information, so we can decide the best architecture. Perhaps interpolation isn't too heavy and we can just interpolate all units in every collision/physics frame, and have the graphics thread do mere linear interpolation in between.
Sounds like a cool idea. It would simplify the code a lot.
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

chuck_starchaser wrote:
klauss wrote:Hey... chuck... if you have the test program, could you turn it into a benchmark?

And measure the frame rate with N simulated items with various time steps, interpolating them ALL at top speed?
This test program that Errol wrote simulates a particle going around an attractor; that's all. It simulates 10,000 steps, and it's ridiculously instantaneous. You press Enter and the results are there. Zero time. :D
Try it. Just uncompress the file, and type ./EJphysics
But I wouldn't know how to integrate it into VS rapidly; I need to do the Unit refactoring thing, --to separate all physics and collision concerns.
No need to integrate it.
Have it simulate N particles around M attractors (make it configurable - a #define somewhere perhaps), and time it.
Count interpolation steps per physics step (that's the achievable FPS).
You'll probably have to port the timing functions from VS, getting portable timing functions that work consistently across platforms is a chore (one that has been done in lin_time.h/cpp).
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

We don't need to keep around the previous physics frame, then; what we need is a new struct to which we only copy X0 and V0, and compute C1 and C2.
In fact, we may not even need the two frames we got now, if we use a different structure for output... (I need to check the rk4 code to see if there would be any aliasing problems reading from, and writing to, the same state) ...though we will need two ouput structures and a flipper flag, for concurrent update/use from separate threads (namely the graphics thread, as I think collisions will be in the same thread as physics). ((On the other hand, if we unify interpolations for collisions AND graphics, then we have no concurrency requirements for this structure. The way this would work is, we get a display list from the graphics thread, periodically updated. Then, whenever we compute interpolations for collisions, we add to our list of units to interpolate the list from the graphics thread, IF the time to the next graphics frame is near. If we do this, then it is the list of final, interpolated positions that would have concurrency requirements.))
Once you have computed C1 and C2, then yes, you no longer need the previous state. You just need the interpolation constants and a range of times when those constants are valid. However, if there is something that interrupts the interpolated trajectory before it reaches the next state, you will need to query for the forces again before you calculate the next physics step.
I don't think having steps longer than an hour would make any difference in terms of computational load, but I could be wrong...
Yeah, it will really depend on how many objects exerting gravitational influence there are in the system. However, if you start doing proto-planetary disks or ENTIRE asteroid fields and what not, you will need to use some tricks to avoid an N-body simulation. One that I am sure you have already thought of is staggering the physics steps of the objects such that only a few objects get a new physics step during the 1/16th second time step. However, even with perfect staggering, getting the force from ALL the other objects in the system is run-time O(N^2) and there needs to be a way reduce this so that we don't turn the simulation into a gravitational simulation, but still get values that look right.

There are two potential ways to do this that I can think of.

One is to use a sphere of influence type calculation. Every gravitational object has a sphere of influence with radius r such that r=someconstant*sqrt(GM). The value of someconstant is somewhat a matter of choice, but it should clearly be large enough so that the sun would have influence on all the objects being simulated in the system. For Sol, assuming you aren't going to simulate the Ort cloud, the farthest object is pluto. So, someconstant should be at least 1300 or so. You determine this value of r for all objects in the system. If the computed value of r is < the spatial extent of the object, then the object exerts no gravitational force on other objects in the system, but it can still be influenced by gravitational effects of others. For other objects, you have determined a distance around them in which objects are affected by their gravity and you ignore the object's gravitational influence outside of this distance. In this example, earth would have a sphere of influence of of around .17 AU. Objects outside of this distance could safely ignore earth's gravity and players would not notice. Jupiter would have a sphere of influence of around 3 AU with the constant I chose above, but the constant can be varied to taste. Also as a reference point, the largest objects in the asteroid belt would have around a .002 AU sphere of influence.

This has several benefits and a couple drawbacks. For benefits, it saves computation time for most objects. It also simplifies orbit stability issues with things like the asteroid belt which consist of many many small objects. The objects in the asteroid belt will only rarely experience gravitation from other asteroids and will, for the most part, travel in a simple low eccentricity orbit around the sun described by their initial conditions. However, occasionally, asteroids will find themselves nearby enough to exert gravitational force.

The downsides are a lack of L3,L4, and L5 lagrange points. The L4 and L5 can be represented accurately using a value of someconstant about 10 times bigger than the nominal value I gave above. Another factor of 2 above that should model the L3 point as well, but it will need to be determined if those values are preferred. The other downside is that it will need to be implemented cleverly to avoid a variety of problems. The first is how to figure out which sphere's of influence and object is in. Naively, to do this, one would need to go through all the existing gravitational bodies and check the distance between them. However, that isn't really any faster than going through and computing the force from each object individually and adding it to a running total as before. I don't have a definitive answer to this yet. I am thinking each object would keep a list of which objects are influencing it gravitationally that would be updated less frequently than a physics time step. So, each physics time step, the object would only need to go through the objects on its list and not through every object in the system. However this can have problems if the object travels through a sphere of influence faster than the frequency with which sphere's of influence are checked.

Anyway, the other possibility is to do something like a PIC simulation. What you do, is every once in a while, you calculate the gravitational potential at a fixed number of points in space. You remember at least a few of these past potential grids. This allows you to interpolate the gravitational potential spatially and extrapolate it temporaly to compute the gravitational potential at any time and place. Taking the gradient of the potential is also easy since you have the interpolation constants. So you can get the force anywhere you want. The downsides of this are that you need a large number of grid points. So this only works if you have a WHOLE LOT of gravitational objects....more than the number of grid points you will need. Dynamic grids around large objects(where you need the grid spacing more fine) are a possibility, but they are very complex.

Bottom line: I don't know how to model the gravitational forces on an object such that it "looks" right, but isn't as computationally expensive as what is actually right. I am curious as to what you guys have in mind.
but then I'm not sure where the initial conditions for such a vast and chaotic set of objects would come from.
I can do this. I just need to know how we are going to handle gravity if there are large numbers of objects in the system.

On interpolations:

The nice thing about an object on an interpolated path is that you can call for its position at any time with minimal computational effort. So, objects like planets that may use very long time steps and may not need to constantly interpolated (ship isn't looking at them and it is out of radar range) but still need a current position for gravitational calculations can essentially treat their getCurrentPosition(time) as a position variable that is always up-to-date and accurate. So, whenever you need it, and for whatever purpose, it is ready to go with a 6 flop calculation. There is no need to store the "current" position. You just have the previous position(including x0,v0, and t0), the next position, and the interpolating constants C1 and C2. Then when you need the position, you can compute it real quick given a time whether it is for display purposes, radar blips, or accurate gravitation, or anything we might think of later down the line. This also allows for the potential of light transit delay for objects that are far away by displaying position from (time-d/c) instead.
I think I'm thinking now, though: Angle and axis is what we need for rotational position, IF we want to interpolate. But for
rotational speed and acceleration we could continue to use quaternions, I suppose; though I'm not sure I understand how the math for expressing rotational speed with quaternions works... Do we have the same problem with having a rotational speed limit beyond which it would have to become a non-unit quaternion or something?
I will need some time to work it out. First, rV0 is not a quaternion. It is a vector. The direction of the vector is the rotational axis and the magnitude is the amount of rotation about that axis. If you write out the equations like a physicist would (ignoring quaternions), you can use exactly the same formula as before except you replace previous quantities with their rotational equivalent. However, then the result you get for the rotation is a vector...not a quaternion. See technically, a rotation can be represented as a single non-unit vector whose direction is the axis of rotation and magnitude is the amount of rotation in radians. The trick then is translating the resulting vector into a rotation quaternion.

delta_theta_vector=rv0*t+rC1*t^2+rC2*t^3


Then, I think delta_theta_vector becomes a quaternion by just constructing the quaternion (cos(mag(delta_theta_vector)/2),unit(delta_theta_vector)*sin(mag(delta_theta_vector)/2))

Then you apply that rotation quaternion to the previous state's quaternion to get the orientation at time(t).

The problem is that C1 and C2 depend on x1-x0 which needs to be in vector form. Anyway, I will write up a rotational interpolation test program eventually to figure that out. Until then, we could just use linear interpolation for rotations. That is straightforward. (I explained it in response to Klauss's question).
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

@Klauss: Hmmm... I'm not sure what exactly you want to evaluate the performance of. Cubic interpolations take only 3 multiplications and 3 additions. Well, this is vector operations, so that would be 9 multiplications and 9 additions; but say 20-20 once we add rotational interpolation. Still, probably not much more at all than interpolations take right now; and infinitely less than an rk4 step; --ridiculously less. Haven't counted rk4 ops, let me see...
helper():
muls: 3 3 3 6 3 4 16 1 4 3 = 46
adds: 3 3 3 2 2 16 1 1 = 31
divs: 1 1 3 = 5
sqrt: 1 1 = 2
trig: 2 = 2
And helper() is called 3 times from rk4Step() :shock:
And rk4Step is pretty heavy, too.
But it runs 10,000 times in zero time, like I said.

Okay, hold on, let me make it run a million times and time it....

3 seconds for a million rk4 steps in my crappy machine.

Hold on, let me try 10,000,000...

27 seconds for 10 mega rk4 steps; so, 2.7 microseconds.
Now, let me add 10 interpolations per step:

32 seconds.
Okay, so the interpolation time is not completely trivial; and right now it's only interpolating position; no rotation.

@Errol: Missed your post; reading now.
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

On another topic: We were discussing gyroscopes by PM, a while ago, and you seemed to say that gyroscopic effects should result from the physics modeling as it stands. I have a nagging suspicion that you're wrong: Right now, RigidBody has no way to inject torque. It is a passive recepient of torque. When you try to turn the axis of a gyroscope, it resists with like a torque of its own making. Of course this is not so, I realize. But in any case, I don't know how the appearance of such a torque could emerge from RigidBody.
Then again, I don't know how an application of torque to an axis of rotation would be expressed.
It doesn't need to apply any torque. If the moment of inertia is a tensor, then torque about one axis can produce angular acceleration about an entirely different one.

consider an INVERSE moment of inertia {{0,0,0},{0,0,0},{1,0,0}}.
That is:

Code: Select all

0 0 0
0 0 0
1 0 0
Apply a torque of {0,0,1}. When you do the matrix multiplication(torque)*(inverse_moment_of_inertia), you get the vector {1,0,0} for the angular acceleration. You applied a torque about the z-axis and got an angular acceleration about the x-axis. The off-diagonal terms can do this.


On the minor discrepancy. The error is mine I believe. Try this instead for the interpolation function with the constants just as you calculated them.

x=x0+v0*t+C1*(t^2)/2+C2*(t^3)/6

Or just, divide the formula for C1 by 2 and the formula for C2 by 3 and leave the interpolating function as is. I apologize. I confused myself when I switched from using "a" and "j" to using C1 and C2

Just to make it absolutely clear

Option 2: Cubic interpolation
C1=(-1*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
C2=(2*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
x=x0+v0*t+C1*t^2+C2*t^3
v=v0+2*C1*t+3*C2*t^2
a=2*C1+6*C2*t
Note that when t=dt, x=x1 and v=v1, but a!=a1

The above should work.
Errol_Summerlin
Merchant
Merchant
Posts: 46
Joined: Thu Mar 11, 2010 4:10 am

Re: rk4 algorithm sent to you chuck

Post by Errol_Summerlin »

Haven't counted rk4 ops, let me see...
You are forgetting the big one for rk4. getForceTorqueandMassLoss will need to cycle through all the gravitational bodies in the actual game. It will need to query for player input such as thrust and possibly compute drag (if using atmospheric braking). And rk4 calls it 5 times. This is where interpolation will really save time. That function is going to be more expensive than all the math in an rk4 step.
32 seconds.
Okay, so the interpolation time is not completely trivial; and right now it's only interpolating position; no rotation.

@Errol: Missed your post; reading now.
I missed yours too. The curse of the long-winded poster. Interpolations are not trivial for the test case because the test case is a very simple getForceTorqueandMassLoss function. If that becomes more computationally expensive, the run time for rk4 goes up a lot, but the interpolation is still the same speed.

Try making an N-Body gravitational simulation for a test run and watch your rk4 run-time go as N^2. With 10 objects, rk4 is now eating up 500 times more time than the interpolation is.
chuck_starchaser
Elite
Elite
Posts: 8014
Joined: Fri Sep 05, 2003 4:03 am
Location: Montreal
Contact:

Re: rk4 algorithm sent to you chuck

Post by chuck_starchaser »

Errol_Summerlin wrote:The downsides are a lack of L3,L4, and L5 lagrange points. The L4 and L5 can be represented accurately using a value of someconstant about 10 times bigger than the nominal value I gave above. Another factor of 2 above that should model the L3 point as well, but it will need to be determined if those values are preferred.
Golden post! L4 and L5 would be nice to have. To hell with L3, I'd say.
However, that isn't really any faster than going through and computing the force from each object individually and adding it to a running total as before.
Exactly.
I don't have a definitive answer to this yet. I am thinking each object would keep a list of which objects are influencing it gravitationally that would be updated less frequently than a physics time step. So, each physics time step, the object would only need to go through the objects on its list and not through every object in the system.
Exactly^2.
However this can have problems if the object travels through a sphere of influence faster than the frequency with which sphere's of influence are checked.
I think collision detection feedback would cause the simulation to accelerate as the objects approach each other, in many cases; and if not, probably the effect doesn't warrant consideration.
That is, if two small asteroids pass each other at high speed, the gravitational attraction is too short-lived to integrate up to too much of a difference anyways; unless they pass each other really close, in which case collision detection would have triggered a temporary acceleration of the physics. Problem solved.
but then I'm not sure where the initial conditions for such a vast and chaotic set of objects would come from.
I can do this. I just need to know how we are going to handle gravity if there are large numbers of objects in the system.
First we need to answer how many gravitational bodies we intend to simulate, I guess; and only Klauss knows that, because my answer would have been 42, whereas he speaks in the thousands. Me, I don't understand. Suppose we wanted to model Saturn's rings as all the individual chunks of ice: What for? If you were there, you'd spend months sitting on one ice chunk, looking around taking photos, comparing to pictures to those of a month earlier, to try to tell if any of the chunks nearby moved. It probably takes years for any two chunks that orbit each other to complete a revolution. We should simulate things that make a difference to the player's experience.
In a protoplanetary disk, same principle applies: Most stuff moves slowly relative to neigboring stuff; or else you get something zooming by you so fast it doesn't even impress your retina. Collisions probably outnumber graviton exchanges :D, and such collisions would best be faked than simulated (just have light flashes happen randomly in the distance). So, I don't understand the motivation for having thousands of units, much less for having more than 42 gravitational attractors.
As for producing data for thousands of chunks of ice, that would make heavy use of rand(), I suppose.
On interpolations:

The nice thing about an object on an interpolated path is that you can call for its position at any time with minimal computational effort. So, objects like planets that may use very long time steps and may not need to constantly interpolated (ship isn't looking at them and it is out of radar range) but still need a current position for gravitational calculations can essentially treat their getCurrentPosition(time) as a position variable that is always up-to-date and accurate. So, whenever you need it, and for whatever purpose, it is ready to go with a 6 flop calculation. There is no need to store the "current" position. You just have the previous position(including x0,v0, and t0), the next position, and the interpolating constants C1 and C2. Then when you need the position, you can compute it real quick given a time whether it is for display purposes, radar blips, or accurate gravitation, or anything we might think of later down the line. This also allows for the potential of light transit delay for objects that are far away by displaying position from (time-d/c) instead.
Exactly. Well, interpolations are not 6 flops, exactly. Remember that X0, V0, C1 and C2 are vectors, not scalars. Secondly, there's also the rotational interpolation, and I'm not sure how the computation goes for that; but I think it's roughly 20 multiplications and 20 additions; that's a whopping 40 double precision flops per interpolation. At 9 and 9, (position interpolation only), they added 5 seconds to a test where there were 100 million interpolations, so that's about 50 nano-seconds in my machine. At 20-20 it would be roughly 100 nanoseconds. So, 10 interpolations per microsecond. If we have 1000 units to interpolate 60 times per second, that's 60,000, so 6 milliseconds per second, or 0.6% of CPU time.
I think I'm thinking now, though: Angle and axis is what we need for rotational position, IF we want to interpolate. But for
rotational speed and acceleration we could continue to use quaternions, I suppose; though I'm not sure I understand how the math for expressing rotational speed with quaternions works... Do we have the same problem with having a rotational speed limit beyond which it would have to become a non-unit quaternion or something?
I will need some time to work it out. First, rV0 is not a quaternion. It is a vector. The direction of the vector is the rotational axis and the magnitude is the amount of rotation about that axis. If you write out the equations like a physicist would (ignoring quaternions), you can use exactly the same formula as before except you replace previous quantities with their rotational equivalent. However, then the result you get for the rotation is a vector...not a quaternion. See technically, a rotation can be represented as a single non-unit vector whose direction is the axis of rotation and magnitude is the amount of rotation in radians. The trick then is translating the resulting vector into a rotation quaternion.

delta_theta_vector=rv0*t+rC1*t^2+rC2*t^3


Then, I think delta_theta_vector becomes a quaternion by just constructing the quaternion (cos(mag(delta_theta_vector)/2),unit(delta_theta_vector)*sin(mag(delta_theta_vector)/2))

Then you apply that rotation quaternion to the previous state's quaternion to get the orientation at time(t).

The problem is that C1 and C2 depend on x1-x0 which needs to be in vector form. Anyway, I will write up a rotational interpolation test program eventually to figure that out. Until then, we could just use linear interpolation for rotations. That is straightforward. (I explained it in response to Klauss's question).
Can't find it.
Anyways, this stuff of converting between vectors of rotation, angle and axis, and quaternions, I think has got out of hand. Maybe we should just standardize on vectors of rotation for everything, and kiss quaternions good-bye. The point of using quaternions was that they were supposedly easy to interpolate; but by now we know better.
You are forgetting the big one for rk4. getForceTorqueandMassLoss will need to cycle through all the gravitational bodies in the actual game. It will need to query for player input such as thrust and possibly compute drag (if using atmospheric braking). And rk4 calls it 5 times. This is where interpolation will really save time. That function is going to be more expensive than all the math in an rk4 step.
Touche.
I hadn't forgotten, exactly; I just didn't consider it part of rk4; --dismissed it as "external code"; but I'd forgotten the fact that it gets called multiple times.
I guess the question of why we need so many attractors comes up again. I don't think we do. I think the only attractors we need are the star, the planets, and the moons. Period. The gravitation of asteroids is too feeble to have any effect perceivable to the player.
Try making an N-Body gravitational simulation for a test run and watch your rk4 run-time go as N^2. With 10 objects, rk4 is now eating up 500 times more time than the interpolation is.
It had better not... I think rk4 is heavy enough without gravity :D
Just kidding.
I think if we have 42 attractors max, and then knock it down to 5 using sphere of influence, we should be alright. Perhaps only one of the 5 attractors (the strongest) should be fully computed for rk4 sub-steps.
Remember that we're only going to simulate a system while the player remains in it; so it won't be years, --except in the case of masochists and maniacs.

EDIT:
Missed the second-last post. Reading now.

EDIT2:
C1=(-1*(2*dt*v0 + dt*v1 + 3*x0 - 3*x1))/Power(dt,2)
C2=(2*(dt*v0 + dt*v1 + 2*x0 - 2*x1))/Power(dt,3)
x=x0+v0*t+C1*t^2+C2*t^3
Works PERFECTLY now :D

EDIT3:
R.E. gyroscopes: Gottcha.

EDIT4:
Damn! Tensor of inertia is a matrix, not a vector. Couldn't it be expressed more succinctly as three quantities?, --namely, the moment of inertia around each axis? I would envision us reading just 3 quantities from units.csv, then expanding that to a matrix, if necessary (in RigidBody's constructor); and maybe the inverse tensor of inertia as well, if it speeds up calculations.
klauss
Elite
Elite
Posts: 7243
Joined: Mon Apr 18, 2005 2:40 pm
Location: LS87, Buenos Aires, República Argentina

Re: rk4 algorithm sent to you chuck

Post by klauss »

chuck_starchaser wrote:@Klauss: Hmmm... I'm not sure what exactly you want to evaluate the performance of. Cubic interpolations take only 3 multiplications and 3 additions. Well, this is vector operations, so that would be 9 multiplications and 9 additions; but say 20-20 once we add rotational interpolation. Still, probably not much more at all than interpolations take right now; and infinitely less than an rk4 step; --ridiculously less. Haven't counted rk4 ops, let me see...
helper():
muls: 3 3 3 6 3 4 16 1 4 3 = 46
adds: 3 3 3 2 2 16 1 1 = 31
divs: 1 1 3 = 5
sqrt: 1 1 = 2
trig: 2 = 2
And helper() is called 3 times from rk4Step() :shock:
And rk4Step is pretty heavy, too.
But it runs 10,000 times in zero time, like I said.
You're ignoring the impact of memory management and bandwidth... if you run it a million times for only one particle, the working set is probably small enough to fit on registers.

A benchmark closer to reality would be simulating multiple particles and multiple attractors, to tax branching and memory just as well as arithmetic.

But I concur, interpolation is very simple when done in vectorized fashion - I just wanted to confirm that given that you already have a test app running the rk4 and interpolation algorithm.
Oíd mortales, el grito sagrado...
Call me "Menes, lord of Cats"
Wing Commander Universe
Post Reply