Refactoring joints

[Julio Jerez]

Did you check out? Is actually running with a simple Euler, I am writing a mid point euler just for reference, because I know the mid point is some time actually worse, but it can be used as the base for a second crack at the predictor corrector integrator.
Predictor correction is good and bad. Good because is adapt to any degree of acuracy, but that flexibility is its downfall. This is because prediction correction is based on averaging the derivative a one point with the derivative at the end point to calculate a new point. This repeat until the average does not changes.

The problem is that since it never takes derivative at intermediate points, is possible as it usually is with non linear functions, that the derivative are smooth at the end but not at the middle, example of this are collision or constraint forces.
But anyway for a free body, it should be ok.

What pencil are you talking about?

[JoeJ]

Did you check out?

Yes, but i see it gains energy. Don't you use RK4 by default?
I would have no complaints against an integrator that looses energy for such rare 'phenomena'

What pencil are you talking about?

You remember, the example of the pencil without gravity.
I would assume if you apply an impulse to the very end, it would only spin but not move.
And if you apply it at the com, it would only move but not spin.
But laws of physics claim linear velocity is the same in both cases.
It makes sense angular motion is faster on one end and slower on the other relative to com velocity, so it cancels each other out, but still this won't get into my brain without a doupt.
I'd really like to see an experiment proofing this.

[Julio Jerez]

Oh The engine use RK 4 only for islands that has joints. free bodies has a special predictor corrector that is faster, and according to teh book more stable.

The test, yes it gain energy, like I said I commented out the predictor corrector integrator, because is smooth the derivative too much.

to do s quick test of how no linera teh derivatoev are, I simpel placxe in a loop and devide teh dt by 8 and still gain energy, of couse tshi does nothong becaus eis just a small Euler step that gain less enery per step, but take more steps.
I am now writing the predictor corrector again.

bu tshi time I am doing the inegration wit the step, tah was a different tha before where I try to reuse the same velocity integration and calcuale teh average velocity.
bu thsi is too slow and does not seem to work.
anyway this is a simple task.

an I knwo thsi si kind of liek a pedante trhong to do, but I liek to have the thong doing eh right thing and applity option to remove thenm rather that not have then and have to wriet special code to get the funtionality.

I do see that this is a tall call, because integrating bodies with high speed with large difference in the inerta meatris is a extreme non linear.
you can see how that demo is quite trivial, the difference in inertial is very small and the angular velocity is not that high, and it takes a few seconds to blow up.

We may need to have this as an option. also the use can controlly by setting the inertia to be uniform. or very close to uniform by having at least two axis of symmetry.

as for the pencil thing.
if you apply an impulse of a force to a body, the body will get a change of momentum.
it like using a level to move an object.

it easy to demonstrate that and yes in fact it will move the center the same regardless where the apply point is. but that only numerically, in practice the center move different base of where the point of action is because most objects aren't ideal rigid bodies, so the impulse will flex long thin bodies.
I know is sound estrange but it is a fact. I made an image to demonstrate it.

Impulse.png (5.3 KiB) Viewed 3823 times

you can see that in both cases, the total impulse or force applied to the rigid body
is: P = [p, r x p ]

I assume you agree that a rotational impulse applied to the body the right will not move the center of the body. only the linear part can do that.

[JoeJ]

I assume you agree that a rotational impulse will not move the center of the body. onley teh linear part can do that.

yes yes, i agree. Also i remember i did this correctly in my own phyiscs sim ages ago, and the assumption hitting exactly at the end would cause no com movement at all is just stupid (major progress ). But still i wonder the rotational motion happens for free and takes no extra energy to get started.
Hmmm... i think i have it: It takes no energy because there is no resistance assuming vacuum. That's it! I did not thought about air resistance but still assumed it to be there. Got it now

[Julio Jerez]

Ok I do not have any bug in the predictor corrector integration, the problem is that the equation is quite nasty, basically it comes down to integrate this expression.

I * a + w × (I * w) = T

This is valid in all frame of reference but calculation of derivative of inertia is a complex problem on any arbitrary frame that is not aligned to the principal axis.
Choosing the principal axis as the frame of reference, the inertia is a constant, and the second expression become a skew cross product, that in effect change the inertia matrix.
This makes easier to extrapolate one step to get the value of the derivative and use the average coeficintes.

Here is the huge problem, the second expression has three carasteristic form. Two of them with a positive coefficient and one with a negative coefficient.
The two with positive coeffiecint has a solution of the form of and exponential decay, and these are the rotation in which the maximum magnitude of w match either the larger or the smaller principal axis.
These integration is very stable, since the errors in angular velocity on the other two axis decay rapidly.

The third characteristic equation has a negative coefficient and the solution to this is a sin wave of a relatively short period. And that's the problem with predictor corrector, it is incapable to average such fast changes, this is why it almost work for some cases, and for others it simple smooth out the error and align the body axis to the axis with max spin.
It is in fact quite stable, and does not violate the conservation of angular momentum or energy, but is not write.

This is a case for using implicit derivative, basically find the effect of a small change of the angular velocity on the other two axis, and use that the expand the above expression.

This is very similar to what the engine does in the joint jacobian already, which is why the system is stable for linked bodies. Is not a big problem, just need to write in paper.

The one problem with implicit derivatives is that it act the opposite way than explicit method.
With explicit methods in general the integration is above the ideal plot, this is because it takes derivatives and advance the as straight lines, taken more only average them at different sub steps.
The result is that the area under the curve will always be larger that the ideal plot area.

Implicit integration, at its core uses series expansion to predict the derivative at some time in the future, and uses that derivative as if it was the derivation at time t.
What this means is that it the derivative intersect the ideal plot in the future.therefore the integration will always be a smaller area that the ideal curve.

What this means is that it will rotate for a while and it will gradually stops as if there was some damping.
This is inevitable, but given the option between gaining and losing energy, I take losing, but this does not means it can not be annoying in some cases

[JoeJ]

With explicit methods in general the integration is above the ideal plot

Sounds like if you could use a second but different approach for approximization which tends to result below the plot and alternate both methods feeding each other to get a very accurate result after a small number of iterations. (math god has spoken )

[Julio Jerez]

ye the explicit method get better as the step is smaller, but its never correct, you can see that is newton every where a spring damper is used, you never get a perpetual sin oscillation even when teh damper is zero, and teh time step is very small, but this is not really a problem since we can make an options.
I put the expression in Mathematica to get the partial derivatives and this is what I get
Mathematica script

Code: Select all

Wxx = wx * Ix + (wy * wz * Iz - wz * wy * Iy - Tx) * dt
Wyy = wy * Iy + (wz * wx * Ix - wx * wz * Iz - Ty) * dt
Wzz = wz * Iz + (wx * wy * Iy - wy * wx * Ix - Tz) * dt

CreateDocument [{TextCell["Wx ="], Wxx,
               TextCell["dwx/dwx ="], D[Wxx, {wx, 1}],
               TextCell["dwx/dwy ="], D[Wxx, {wy, 1}],
               TextCell["dwx/dwz ="], D[Wxx, {wz, 1}]}]
            
CreateDocument [{TextCell["Wy ="], Wyy,
               TextCell["dwy/dwx ="], D[Wyy, {wx, 1}],
               TextCell["dwy/dwy ="], D[Wyy, {wy, 1}],
               TextCell["dwy/dwz ="], D[Wyy, {wz, 1}]}]

CreateDocument [{TextCell["Wz ="], Wzz,
               TextCell["dwz/dwx ="], D[Wzz, {wx, 1}],
               TextCell["dwz/dwy ="], D[Wzz, {wy, 1}],
               TextCell["dwz/dwz ="], D[Wzz, {wz, 1}]}]

the derivative is:

Code: Select all

dWx/dwx = Ix
dWx/dwy = dt * (-Iy * wz + Iz * wz)
dWx/dwz = dt * (-Iy * wy + Iz * wy)

dWy/dwx = dt * (Ix * wz - Iz * wz)
dWy/dwy = Iy
dWy/dwz = dt * (Ix * wx - Iz * wx)

dWz/dwx = dt * (-Ix * wy + Iy * wy)
dWz/dwy = dt * (-Ix * wx + Iy * wx)
dWz/dwz = Iz


ass you can see the derivative for a mass is in fact a correction of the diagonal matrix, now I only need to pass this and find the common expressions.
The correction is no symmetric, so it nee fo complete matrix solver. to calculate the inverse

take for example the first row reduces to

Code: Select all

dWx/dwx = Ix
dWx/dwy = (Iz  - Iy) * wz * dt
dWx/dwz = (Iz  - Iy) * wy * dt

you can see what happens, when a body has symmetry along y and z axis the Iz == Iyy
and the of diagonal values are zero.
when the Iz is diffrent that Iy you can see that there is a off diagoila coupling effect that is a funtion o fteh angula velocity of the on the oteh two axis.

also we can see that we can step that and any fration of the time step. and instead fo calculation omega in one step wie can do it is say tow or three base of hwo big is wy * dt or wz * dt
this seems promising

[Julio Jerez]

Ok Joe I commited the code now, and as predicted the work is now reversed.
The the total energy decay frame by frame.
also as predicted the larger the sub step the fastest the decay so I am doing two sub steps.

doing too many substeps is a mistake because is result into lots of noise added by small floats added to large floats, and also this requires a general Gaussian elimination matrix solver per step.

would you please check it pout and tell me what you think?

Notice that the lost of momentum only happen if the Inertia is asymmetric, for symmetric inertia the matrix coupling are all zeros, and only teh main diagonal is left, so it become a expensive simple Euler

we can add a check that only do it for Asymmetric inertia, but that just an optimization.

[Julio Jerez]

I made an image so that I can explain teh difference between explicit and implicit integration

integration.png (11.75 KiB) Viewed 3813 times

you can see hwo the secund method lose energy by design. of teh tow method the explicit integration is far more accurate, because is use teh actual derivatives.
teh implicit method has to use and approximation of the derivative by expanding teh funtion using Taylor or use forward numerical differences. This on itself is already and approximation, since the series has to be truncated to a first order.

there is a method is tha very accurate which teh newton method which truncate the series at the derivation and the secund derivation is call the Hessian and teh firs derivation teh gradient, bu this method is too hard to use with vectors, is better for optimizing scalar energy functions.
I do not think I will go there.

also this video is a very good and clear explanation of the effect.
https://www.youtube.com/watch?v=cnF4G-InNqU

[Julio Jerez]

Hey Joe now that we have the expression, it is actually very intitive
after all the math hankypanky is all come down to this.
The angular momentum at time t must be equal to the angular momentum at t + dt
which is basically this equation

Code: Select all

// dW/dt * I + w x (w * I) = T
// dW * I + w x (w * I) * dt = T * dt
// since c * (a x b) =  (c * a) x (c * b)
//  we get
// dw * I + (w * dt) x ((w * dt) * I) = T * dt
// dw * I + dw x (dw * I) = dT

//dw * I + dw x (dw * I) - dT = 0


in the last equation solving for dw that make the expression zero at any dt, is the solution.
notice that if dT is zero, it has a trivial solution which is dw equal zero. which of course is the case for a non spinning body. for spinning bodies is has a non trivial solution which requires the inertia matrix changing with time.
what happens is that the Taylor expansion of the expression makes an artificial inertial matrix that is not diagonal. and that's the explanation of why it works. after that it becomes very intuitive.

I applied the final adjustment an now is near perfect. it runs for a long time and is does not gain or losses angular momentum or energy in the numerical sense, if you let it go for very long time it degenerates but it take hours and there is nothing we can be done about that because that's the nature of numerical integration, maybe running at smaller time step.

On this joe

I know this video, or a similar one in a space station.
Interestingly all this stuff - also the flywheels - looks totally logical to me and i expect this behavior

you will be suprrice hwo many of the so called adavace engines that invent brand new Laws of physics aren't even aware of that problem.
These engine talk in the language of:
some engine do 1 impulse pass, and 4 position correction,
multiply velocity by 0.7,
long distance constraint, and stuff like that
scale the Inertia matrix by 10.
clipping impulse, adaptive forces
add joint n time in some specific order
and stuff like that.

Ok now to back to the two dof joints.

[JoeJ]

So conversation of angular momentum did the trick?
That's a really useful law of physics. Recently i used it for automatic UV mapping, which works like a soft body simulation of triangles in 2D. Results are really good, better than Microsofts UV Atlas, and no need for tweaking parameters.

Gyroscopic looks great now

[Julio Jerez]

I now completed the cork screw joint, and it is derived from the slider.

I could be a sub classed from a hinge but slider is betters because there is another joint, sliding contact, which is used for wheels with suspensions. the cockscreew does not have practical usage other that is is very simple an can ge used as the model to make other joints.

Now I will do the Sliding Contact next, which should be a copy, paste and edit of the cork screw.
hopefully by the end of this week we can cover all the basics joint and get back to the active Character

[Julio Jerez]

the slider contact was piece for cake. Later I will make a spring so that people can see that this is a wheel.
I added a funtion to trace Energy ans angal momentum
and the are the results

on the beginning

E = 550.845581 L (0.007644 0.032097 55.084538)
E = 550.845642 L (0.007911 0.031870 55.084541)
E = 550.845581 L (0.008244 0.031671 55.084541)
E = 550.845520 L (0.008644 0.031515 55.084530)

after few minutes
E = 552.705872 L (-0.048777 0.198817 55.177235)
E = 552.738647 L (-0.160179 0.045880 55.178997)
E = 552.777527 L (-0.251570 -0.140124 55.180443)
E = 552.823669 L (-0.314742 -0.357644 55.181438)
E = 552.877869 L (-0.341150 -0.603136 55.181854)

as you can see it does gain some energy but this is after few minutes of simulations, there is not way that I know this can be resolved, if there was, I will be inventing a way to solve non linear equations by using linear approximation.
I am actually very please with those results.

[Julio Jerez]

Ok I think I now completed all of the 1 and 2 standard joints.
they now inherited from some joint, for example the corkscrew and the sliding contact inherit from the slider so they has all the properties of the slider.
the universal inherit from the hinge, so it has all the features of a the hinge, plus its own.

My mistake was that I was too greedy, I wanted to inhered from a common joint and that makes it too difficult, instead is like evolution, joints has common heritage and is a leaf joint needs a functionality from some other joint that is on a separate branch, they copy the functionality (converge evolution).
this is the case of the sliding contact, is a slider but is also a hinge, but since it has more in common with a slider, it is derived from slider and copy, paste and edited the hinge functionality.

This seem to be a lot more consistent and simpler. because now similar Joints have the same reference,
for example before if I was using a hinge and I wanted to change to a slider contact to add suspension, they were using different frames, so lot high level setup had to changed, now the joint can be a slider contact and making a is just and option.
And example of that is a tank vehicle, some wheels has suspension and some do not, now they all can be slider contacts.
This does not solve everything but I like better.

now to the 3dof joints, and these are there one lineage.

Ah Joe, if you sync, I committed two Universal joints spinning on the two axis.
you will see that one bends and the other does not.

this is where the physics get strange, the joint on the right, is the one that is correct.
while the one of the left does not bends but is wrong.

here is the reasoning and why the change I made to move the derivatives to the joint payoff.

Remember the Gyro, if the object has one angular velocity and the gravity apply a torque to try to toppled the disk, the joint generates a side torque that gives a processing angular velocity.

here we have the pivot at the center of mass but the wheel has a precessing velocity, therefore there move be a conservative torque acting on it to bend to the side.
you can see that effect in this video, at 2:30
https://www.youtube.com/watch?v=NeXIV-wMVUk
see that when the professor hangs a weight, the wheel changes its velocity and bends to match the conservation of angular momentum.
then he removes the weights but the wheel keeps the bend angle regardless.
what happen is that the is the orientation that matches the main angular velocity and the processing velocity.
The moral is that each processing velocity has a unique side tilt angle. and that's inevitable.

This is what was different before. I was using this angle as an contraint error violation and apply a torque to fix it, but that generates a perpetual non conservative torque in the system that try to rotate the entire system.

now the joint does what the derivative say they has to do, and the error violation is added by the joint implementation.

in the case of the universal joint now we have the option
joint1->SetHardMiddleAxis(0);

which is applied to the left joint and that one use the angular error to generate a torque that keeps the wheel up.
But remember that joint is actually applying a torque to the two bodies, if the parent body was has a non zero mass, that system will be spinning around in way that would not make sense, because it getting a torque from the thin air. the one on the right will be acting as a body that will be moving in place, because after the momentum align all torque will be conservatives and produce zero work.

This is what explain all the problems we were having when using joint that rotate in three axis to control the balancing gait, is was clear that the joint was adding these weird torques that were coming form the conversion of error violations to joint accelerations.
This is new information, we did no have before.

anyway I will continue with the 3dof joints now.

[JoeJ]

Did some side by side comparision test with old July version again to ensure your improvements have no bad side-effects.

All fine

(I updated my joint test on github because your ball socket does some angular stuff now, so i made it independent of this)