Physics - Inertia - message from Danske

By: Danske - danske
Angular Momentum (probably Yet Again)  
2006-07-08 22:23

I've been working on my own physics simulation (semi-gamish) and this site came up in some search. I have a few quibbles/suggestions regarding some of the angular momentum pages. 
 
I haven't read every page on the site from start to finish, and there could definitely be more links when concepts are referenced (like intertia matrix), so what I'm saying may be covered somewhere I didn't read. My math notation is also bad, so some creative interpretation may be required. Please cut me a little slack, I'm trying to make a positive contribution, really. ;) 
 
So here are my comments on the following pages: 
https://www.euclideanspace.com/physics/dynamics/energy/rigidbody/index.htm (Inertia Matrix) 
https://www.euclideanspace.com/physics/dynamics/rotation/index.htm (spinning pencil and other principle axis examples) 
 
It seems to me that it's better (or at least enlightening) to consider torque, angular acceleration, and angular velocity starting from angular momentum. 
 
So my attempt at equations would be: 
(with T = torque, L = angular momentum, a = angular acceleration, w = angular velocity, I = moment of inertia; please interpret the correct representations and operations) 
 
L = Iw 
dL/dt = d(Iw)/dt 
= I'w + Iw' 
= I'w + Ia 
 
T = dL/dt 
= I'w + Ia 
 
(all analogous to: 
 
p = mv 
dp/dt = d(mv)/dt 
= m'v + mv' 
= m'v + ma 
 
F = dp/dt 
= m'v + ma 
 
if m' = 0 
F = ma ) 
 
Angular momentum is a vector and torque is a vector expressing its change over time. That's all easy to understand and visualize. The problem is how angular momentum relates to the motion of an object. 
 
So let's just start here: 
 
L = Iw 
 
'w' is a vector representing the axis of rotation and angular velocity about that axis. 'I' is the inertia matrix which combined with w gives us the angular momentum of the object. 
 
The first problem is that of initial conditions. Just as with "p = mv" we must start with two values, calculate the third, and then integrate over time to find subsequent values. It's an unusual problem to be solving for the inertia matrix, so we'll take that as one given and also one that we can calculate or know at any time for any object. At the start of the simulation we could have either the angular momentum given and solve for angular velocity (although given some inertia maxtrices there may be no solution) or be given the initial angular velocity and solve for initial angular momentum. In either case we begin knowing all three. 
 
Let's start by examining an object which begins life with a non-zero angular momentum and experiences no net torque (constant angular momentum). 
 
L = Iw = constant (non-zero) 
dL/dt = I'w + Ia = 0 
 
Now if I is constant over time then angular acceleration is zero and angular velocity is constant over time. The question is if I' can be non-zero. 
 
I believe the answer is yes. Take the following example: 
 
Imagine the point-mass dumbell from the example of inertia matrices in the referenced page above. Lay it in the x-y plane, and rotate it about the z-axis. The angular momentum is pointing straight up along the z-axis (<0,0,something>) and w is likewise pointing straight up. As the dumbell rotates about the z-axis you can see that the inertia matrix will change over time. When it is aligned along the y-axis the Ixx value is at its largest, when it is aligned along the x-axis its Ixx value is near-zero. QED, the inertia matrix changes over time. 
 
Back to: 
 
dL/dt = I'w + Ia = 0 
 
So if the inertia maxtrix does vary over time then I' is non-zero and both I and I' must be recomputed for each time step. Furthermore, as demonstrated in the dumbell example, the rate of change over time of I depends on w (I changes as the object rotates, and rotation is w). 
 
However, for the dumbell lying in the x/y plane and rotating about the z-axis this isn't so bad: 
 
Ixx = ?? (variable) 
Iyy = ?? (variable) 
Izz = lzz (constant) 
Ixy = ?? (variable) 
Ixz = 0 (z = 0, so xz = 0) 
Iyz = 0 (z = 0, so yz = 0) 
 
So Ia = 
 
?? ?? 0 ax 
?? ?? 0 * ay 
0 0 Izz az 
 
and I'w =  
 
?? ?? 0 0 
?? ?? 0 * 0 
0 0 0 wz 
 
(okay, with the simple dumbell example we could solve equations for the ?? in the matrices, but it's not important here.) 
 
With so many zeroes the equation isn't bad. We're left with having to solve: 
 
(from dL/dt = I'w + Ia = 0 ) 
 
ax*Ixx - ay*Ixy = 0 
-ax*Ixy + ay*Iyy = 0 
az*Izz = 0 
 
which produces a = <0, 0, 0>. As the dumbell is already rotating about one of its principle axis there is no angular acceleration. 
 
Because angular momentum is a vector in the global frame of reference but an object's inertia matrix is changing as it rotates, the object's angular velocity may change over time (may just be the axis moving) in order to maintain constant angular momentum. So combined the terms (I'w) and (Ia) determine how an object tumbles when it is not rotating about one of its principle axis. 
 
(For non-rigid objects such as the physics class experiment of a student sitting on a spinning stool with weights in his/her hands changing their moment of inertia, knowing (I'w) lets you solve for the angular acceleration in (Ia); the faster the student pulls in their arms the greater the accompanying acceleration and the greater the initial w the greater the needed acceleration, e.g. if I is changed so that w is doubled, then doubling 2/s to 4/s requires greater acceleration than doubling 1/s to 2/s.) 
 
Now to the example of the inclined dumbell. There is a sign error here, I believe, a value error, and I think it needs a note as well. 
 
If the maxtrix equation for the example is: 
 
[Tx] [ 1 -1 0] [ax] 
[Ty] = [-1 1 0] [ay] 
[Tz] [ 0 0 2] [az] 
 
(the masses need to be 0.5kg each for that to be the matrix, that's the value error) 
(is that a better way to represent matrices in text?) 
 
then I think the solutions are actually: 
 
Tx = ax - ay 
Ty = -ax + ay 
Tz = 2az (unsed in the example) 
Tx = -Ty 
 
so that the acceleration is actually on a line of x = -y. The image(GIF) of the equation also has an (-Ixy) in the upper-right corner of the inertia matrix where it should have an (-Ixz). 
 
One problem I have with the example is that in simulations it seems to me you usually know what the angular momentum at any given moment should be (as you calculate the forces you also calculate the torques and therefore know how angular momentum is changing; integrating torque over time gives you L), and have to calculate angular velocity and angular acceleration from that. The example treats T as an unknown instead of a known. 
 
If we instead proceed through the inclined dumbell example with a known dL/dt then another problem is exposed. Let's say we know we have an external torque producing a change in angular momentum of <0, 1, 0> (along the y-axis). Then the equation gives us: 
 
0 = ax - ay 
1 = -ax + ay 
 
ax = ay 
1 = -ay + ay 
1 = 0 (Ooops!) 
 
The example dumbell with its point masses has a moment of inertia of zero about its long axis. This would mean that however fast it rotates about its long axis it would never generate any angular momentum, or that any torque about that axis produces infinite (or undefined if that's how you prefer to describe dividing by zero) angular acceleration; this is the note I feel should be added. Real objects wouldn't have any axis of rotation about which they had a moment of inertia of zero. So having any angular momentum along this axis is impossible ( L * <1, 1, 0> = 0 ), which is what solving for torque in the example has actually proven. 
 
If we amend the example to give the dumbell a non-zero moment of inertia through its long axis by having four masses of 0.25kg at <1.1, 0.9> <0.9, 1.1> <-1.1, -0.9> and <-0.9, -1.1>, and setting dL/dt to <0, 1, 0> we have the following equation: 
 
[Tx] [ Ixx -Ixy -Ixz] [ax] 
[Ty] = [-Ixy Iyy -Iyz] [ay] 
[Tz] [-Ixz -Iyz -Izz] [az] 
 
[0] [ 1.01 -0.99 0] [ax] 
[1] = [-0.99 1.01 0] [ay] 
[0] [ 0 0 2.02] [az] 
 
ax = 0.9802ay 
ax = 24.75 
ay = 25.25 
az = 0 
 
So with a torque producing a dL/dt of <0, 1, 0> the resulting angular acceleration is <24.75, 25.25, 0>, or very nearly but not quite along the long axis of the dumbell. Why is this? If we look at this revised dumbell's moment of inertia along its long axis it's still very small, only 0.04 by my calculation, compared to the ~4 about any axis perpendicular to the long axis. Any rotation about the long axis provides very little angular momentum, so we need the axis of rotation very close to the long axis if the angular momentum has any significant component parallel to that axis. 
 
To better illustrate that and to provide a sanity check, let's rotate the dumbell so that it's long axis is in line with the x-axis. That gives us: 
 
[Tx] [ 0.04, 0, 0] [ax] 
[Ty] = [ 0, 2, 0] [ay] 
[Tz] [ 0, 0, 2.02] [az] 
 
Hey! It's the principle axes! ;) (Sorry, it's getting later at night as I write this now.) 
 
Tx = 0.04*ax 
Ty = 2*ay 
Tz = 2.02*az 
 
So if Tx (dLx/dt) is non zero it produces a very large angular acceleration along the x axis, as (ax = 25*Tx). 
 
If dL/dt = < 1.414, 1.414, 0 >, at 45-degrees from the long axis again, then: 
 
0.707 = 0.02*ax 
0.707 = 2*ay 
ax = 35.355 
ay = 0.354 
az = 0 
 
The magnitude of the angular acceleration is the same as before, and once again very nearly in line with the long axis, which has been rotated to be the global x-axis now. 
 
If we use this initial I then we can also pretend we started in this orientation with a large L at 45-degrees to the long axis, say L = <100, 100, 0>. With [L] = [I][w] then [w] is along the same vector as [a] was before, that is nearly in line with the long axis but "cheated" just a hair toward the [L] vector, with the same solutions as for a: 
 
100 = 0.02*wx 
100 = 2*wy 
wx = 5000 
wy = 50 
 
L = <100, 100, 0> 
w = <5000, 50, 0> 
 
(because we are not rotating about a principle axis the axis of rotation is not in line with the angular momentum) 
 
Notice that the angular momentum, angular velocity, and long axis of the dumbell are all in the same plane, the x-y plane at the start. I believe that this relation must hold as long as the angular momentum remains constant (net torque is zero). 
 
What does it mean physically for the axis of rotation to be nearly but not quite in line with the dumbell's long axis? Well, if the dumbell is rotating about that axis, and its own long-axis is not in line with that axis, then the long axis must be rotating about that axis, at least at that initial moment in time. Since the axis of rotation has a slight y component the dumbbell is rotating slightly around the y axis, which gives the long axis a slight negative-z velocity at the positive x end. 
 
If the angular momentum hasn't changed (still pointed off at <100, 100, 0>), but the axis of rotation is not in line with that then the axis of rotation must also change. When we solved for the angular velocity we found that we needed a very small component perdendicular to the long axis of the dumbbell in the direction of the angular momentum. In order to keep that the case I think the axis of rotation must itself rotate about the angular momentum vector. If we were already rotating about a principle axis then the axis of rotation would be in line with the angular momentum vector and so it wouldn't move. As this is not the case it must. This changing of the axis of rotation is the a in (dL/dt = I'w + Ia). 
 
As some further examples of this I think wobbling semi-gyroscopes are the best. A "lame duck" thrown football is obviously spinning very nearly along its long axis, while that long axis is also rotating about the angular momentum vector. This is not caused by air resistance and you can toss a football just slightly into the air in front of you with this same behavior. A wobbling frisbee does the same. This ties in to the pencil example. You can't make the pencil rotate about an arbitrary axis, but you can give it an angular momentum pointed in an arbitrary direction about which its axis of rotation will rotate. You can demonstrate this also by tossing a flashlight, mug, or other cylindrical object into the air and trying to almost but not quite spin it about its long axis. Perhaps a bicycle wheel gently tossed and spun in the air will also give you proof that the gyroscopic "wobbling" effects are real and don't depend on external forces. ;)  
 
So all this stems from properly differentiating angular momentum over time, giving us: 
 
dL/dt = I'w + Ia 
 
I don't know if there's an easy way to calculate I and I'. There are those three orthogonal principle axes, but I don't know if you can simply rotate that inertia matrix to find the correct I for any given orientation of the object. I'm not sure the total physical properties of the object needed for the computation of the inertia matrix are completely encoded in those nine numbers. If you have a simple, fairly symmetrical object maybe you can rotate the inertia matrix, but I really doubt a highly asymmetric object can be recontructed just from the matrix. You may be need a specific formula for each object. I'll have to play with some numbers. 
 
Anyway, the inertia matrix is what I was looking for, even though I had come up with the principles I've expounded on here; I didn't have the right mathematics to express it and I thank you for publishing them on your pages. I hope I've helped in return, although it's possible all this was already hidden somewhere on the site that I hadn't read. 
 
Oh, I was also working on the constraint equations for the different joints, trying to model them as simply distances between points on objects that must be maintained. If you could point me at a reference describing the basic joints that would also be appreciated; I don't think I found exactly what I was looking for where I looked, and there might be some more information somewhere that I missed. I mentioned that I'm approaching my effort more as a physics simulation, and my most recent efforts have been in working out using Euler's Method (which I think I independently re-invented from a physics-simulation-application-specific perspective because I had forgotten it from my purely theoretically mathematical classwork. ;) ). I could help write a bit more about that if you'd like. 
 
*yawn* Very late now. Sorry if this has all been covered here already.

By: Martin Baker - martinbaker
RE: Angular Momentum (probably Yet Again)  
2006-07-09 11:06

Thanks, this is great, I appreciate you taking the time to write this up and I definitely want to include this on the website. 
 
I don't think your points are already covered, in fact its made me realise the 'dynamics' pages on the site need a lot of improvement. 
 
I have already linked this message from the 'correspondence' section at the bottom of the pages mentioned (hope that's alright). 
 
The next stage is to incorporate the points into the web pages themselves, so I will: 
* Introduce from L = [I]w as you suggest 
* Make the point much more strongly that, in general, [I] changes with orientation (when in inertial frame coordinates). 
* Link this to the examples as you have done (I need to work through that in more detail to make sure I understand) 
* Make the corrections you suggest. 
 
I think I also need to make some changes to the order of the pages so, for example, the Inertia Matrix is introduced properly before its used. So these changes will take me a few days to make. 
 
> I've been working on my own physics simulation (semi-gamish)  
Sounds interesting, is it open source? or anything you can talk about? 
 
> Oh, I was also working on the constraint equations for the different joints, 
> trying to model them as simply distances between points on objects that 
> must be maintained. If you could point me at a reference describing the 
> basic joints that would also be appreciated 
The only stuff I have on joints is here: 
https://www.euclideanspace.com/physics/kinematics/joints/jointstructures/ 
I would like to improve this. 
 
> my most recent efforts have been in working out using Euler's Method 
> ... I could help write a bit more about that if you'd like 
Yes, I have wanted to have something about Eulers Method and the Lagrangian Method for some time. Anything about either of these topics would be most welcome, but I would feel guilty about imposing on your time too much. 
 
Thanks, Martin

By: Danske - danske
RE: Angular Momentum (probably Yet Again)  
2006-07-09 17:11

 
>> I've been working on my own physics simulation (semi-gamish) 
>Sounds interesting, is it open source? or anything you can talk about? 
 
Right now it's near entirely just scribbling in a notebook. There is this old tabletop game called "Car Wars". It basically simulated driving cars around with guns and other weaponry shooting at each other in a post-apacolyptic future. ;) The time step size was 0.1s/turn, but with all the dice rolling, table lookups, measuring and maneuvering of little playing pieces, and bookkeeping it was very, very far from a realtime simulation. Anything approaching realistic physics was also absent. What began as just having some fun trying to make a more realistic game became writing the mechanics textbook for the class I never took, taking the mechanics section from the physics textbook and expanding it to three dimensions. If I get to heat-induced brake-fade I'll even make it to thermodynamics. 
 
I'm inclined to try for a simulation of the actual physical state of the system instead of just using equations which give the correct behavior; that's just the way I like it to really know the hows and whys of the physics. However, in some cases the physics are just too complex to not use equations (tires would really need FEM to simulate physically), or I can convince myself the equations are not approximations and I understand why they work. In other cases it seems it might just be easier to physically model the system even at the expense of processing time rather than spend the time to derive the proper equations because I'm not concerned about processing time. Yeah, that's not really a game-making perspective, but I do intend to play with it. In a realtime game I'd say you should at least know how much accuracy you're throwing away before you throw it, though. 
 
So right now I'm working on physically modelling the car suspension, along with finally getting started on coding the underlying simulation mathematics, hence my inquiry about joints and having found the site in the first place. The various rods could be modelled as very stiff springs, but very stiff springs results in very high forces and accelerations. If these accelerations are too high for the time step then bad things happen. Which leads us to Euler's Method. 
 
Two appropriate pages would seem to be: 
https://www.euclideanspace.com/maths/differential/index.htm  
https://www.euclideanspace.com/threed/games/options/timestep/index.htm 
 
First off we need to expand the explanation of the Taylor series. 
There are more terms to the Taylor series than in: (x(t0 + h) = x(t0) + h x'(t0) + h^2/2 x''(t0) + O(h^3)). It continues off in an infinite series: 
 
x(t0 + h) = x(t0) + x'(t0)h + x''(t0)/2 h^2 + x'''(t0)/6 h^3 + x''''(t0)/24 h^4 + ... + x('n)(t0)/n! h^n + error 
 
Err...I'm going to write this quickly, so it may not be comprehensive or appropriately split-up for those two pages. 
 
If we're integrating over time, as in a simulation, then we need to be able to calculate the derivative over time for each variable we're attempting to integrate for. For position this derivative is velocity, for velocity this derivative is acceleration. If we could directly calculate the derivative for acceleration than we could use this as well to predict future acceleration. In a physical simulation acceleration (from forces) typically depend on the positions and velocities of all the objects in the system, so we don't have an equation for how acceleration (forces) are changing over time. If we did that would be like having an equation which would tell us that two objects will collide in the simulation before we actually get to that time. This is, however, what we will attempt to numerically emulate. 
 
Rewriting the Taylor series for physical quantities (and in the order usually used in physics texts) we have: 
 
x(t0 + h) = 0.5a(t0)h^2 + v(t0)h + x(t0) 
(the equation for position given constant acceleration) 
v(t0 + h) = a(t0)h + v(t0) 
(the equation for velocity given constant acceleration) 
 
if we had the derivative of acceleration then we could instead have: 
 
x(t0 + h) = 1/6 * a'(t0)h^3 + 0.5a(t0)h^2 + v(t0)h + x(t0) 
v(t0 + h) = 0.5a'(t0)h^2 + a(t0)h + v(t0) 
a(t0 + h) = a'(t0)h + a(t0) 
 
We can estimate a' using our simulation, however. We can calculate the forces in the simulation given positions and velocities off the objects. If only we knew the positions and velocites of the objects in the future we could know the forces in the future. Wait a second, computing positions and velocities of the future is what this simulation is doing. We can use it! 
 
How does this work? We can use our current forces and accelerations to provide us with a prediction of the positions and velocites of all the objects (h) in the future, calculate the forces and accelerations in our predictions, and then use those future accelerations to estimate (a') at (t0). We can then use that estimated (a') to add another term to our Taylor series, improving our estimation of future positions and accelerations, repeating the process until the difference between our prediction and corrected prediction is small. 
 
Here's what it looks like mathematically: 
 
a(t0 + h) = a'(t0)h + a(t0) 
a'(t0) = (a(t0 + h) - a(t0))/h 
(you'll recognize that as the definition of the derivative if h -> 0) 
 
Instead of using the "actual" value of (a(t0 + h)) we use our predicted value of the future. 
 
a'(t0) = (ap(t0 + h) - a(t0))/h 
 
If we put that back into our Taylor series (omitting some of the alebra) and calling our new values "corrected", we have, : 
 
xc(t0 + h) = 0.25*(ap(t0 + h) + a(t0))h^2 + v(t0)h + x(t0) 
vc(t0 + h) = 0.5*(ap(t0 + h) + a(t0))h + v(t0) 
 
So instead of using simply (a(t0)) in the Taylor series without (a'(t0)), we replace it with the average of the current acceleration and our predicted accelerations in the future ((a(t0) + ap(t0))/2). 
 
In terms of the simulation what we're doing is generating the world one time step in the future as if all the forces were constant during that time step, evaluating the forces in the future, and then using those future values to compute new average forces to use over the timestep. If the forces aren't constant over time then our new average forces will result in a more accurate simulation than pretending the forces are constant during each timestep. 
 
As a concrete example let us consider a spring being compressed over time. At (t0) it is producing some force, at (t0 + h) it is producing a greater force. If we integrate position over time with the initial force then we've underestimated the force of the spring during that time step. By "peeking" ahead we see that its force is increasing and that we should use a greater value as an average force over the time step. Doing so gives us a more accurate simulation of the forces from the spring. 
 
In another example which may be more intuitive, consider a car driving off a cliff sometime in the middle of the time step. At (t0) it is still on the ground. At (t0 + h) we find it hanging in mid air. Somewhere in the middle of the timestep the force from the ground when from the weight of the car to zero. Obviously the car has been experiencing an average acceleration somewhere between zero and -9.8m/s^2 over the timestep and our use of a zero acceleration is inaccurate. So we recalculate the future position and velocity of the car at (t0 + h) using an average acceleration of -4.9m/s^2 instead of zero. This simulates the car having driven off the cliff exactly half-way through the time step. Unless we want to instead decrease the time step to more accurately determine when the car went off the cliff then this result of simulating it going off half-way between when we last saw it on the ground and when we first found it hanging in mid-air, this is better than simulating it having just driven off the cliff an infintesimal amount of time before we first found it hanging in the air. On average it will have driven off in the middle of the timestep, which is what using the average acceleration of -4.9m/s^2 over the timestep gives us. 
 
In coding this Modifier Euler method we do need to maintain at least two whole world states, the initial state at (t0) and our predicted state at (t0 + h), in order to compute the average forces to use in our corrected computations. The corrected values (using the computed averages) can then be stored in the memory allocated for the initial state (since we're done with it unless we want to re-correct with another estimate of average forces using the forces we find in our corrected world), so memory usage only doubles (or maybe as bad as triples depending on how we're forced to compute averages of the data). Since we're computing a future world at least twice for each time step the simulation rate will drop, but the accuracy improves more than if we halved the time step. If we're having problems with accuracy, or just want more of it, and can afford the memory usage and the added complexity of computing the average forces, then this may be faster than decreasing the timestep to achieve the same accuracy. 
 
In summary, what this method does is guess and check the world state in the future. We take our initial guess at the future assuming all the forces are constant over the time step, use the future forces in our first guess of the future to provide us with a better guess at average forces over the time step, and use those averages to make a better guess at the future. 
 
If code is your thing, here's a snippet computing a one-dimension mass on a spring: 
 
double x = 10.0; 
double v = 0.0; 
double k = 1000.0; 
double t = 0; 
double t_step = 0.001; 
double t_max = 10.0; 
while ( t < t_max ) { 
double a = ((x * -k) / m); 
double xp = (x + v * t_step + 0.5 * a * t_step * t_step); 
double vp = (v + a * t_step); 
double ap = ((xp * -k) / m); 
x = x + v * t_step + 0.25*(a + ap) * t_step * t_step; 
v = v + 0.5*(a + ap) * t_step; 
t += t_step; 
} 
 
 
Oh, and I played at solving some simple constraints using forces basically using this method. It involves "peeking" into the future with the constraining forces slightly altered, seeing how that change affects the constraint condition, and then rescales the constraining force to satisfy the constraint. You may recognize that as Newton's Method. The problem is that in a system this creates an NxN-dimensional matrix where N is the number of constraining forces. Each constraining force (one of the N) must also be "tweaked" independently to see how it affects *all* constraints (the other N). So the number of iterations goes up with the number of constraints, if using this method. If you consider even the "simple" example of a jointed bar swinging from a pivot on the ceiling then it's still a messy solution. With the complexity of physical laws I'm not sure there's an easy mathematical solution if you want to keep those laws unviolated. But I'll look at the Lagrangian Method and see if I have any insight.