wu :: forums
« wu :: forums - Moment of Inertia by Random Iterations »

Welcome, Guest. Please Login or Register.
Mar 29th, 2024, 5:31am

RIDDLES SITE WRITE MATH! Home Home Help Help Search Search Members Members Login Login Register Register
   wu :: forums
   riddles
   general problem-solving / chatting / whatever
(Moderators: towr, ThudnBlunder, Icarus, SMQ, william wu, Eigenray, Grimbal)
   Moment of Inertia by Random Iterations
« Previous topic | Next topic »
Pages: 1  Reply Reply Notify of replies Notify of replies Send Topic Send Topic Print Print
   Author  Topic: Moment of Inertia by Random Iterations  (Read 4874 times)
Christoph
Guest

Email

Moment of Inertia by Random Iterations  
« on: Jul 5th, 2005, 6:26pm »
Quote Quote Modify Modify Remove Remove

I am in the middle of writing a program (in C) that will calculate the moment of inertia for a 3D rigid-body ("object"). This object is non-symmetrical and not of any standard shape (i.e. sphere, cube, etc.). I have the exact x,y,z-coordinates of every point in or on this object (actually, it is a protein structure in PDB format. So, I have the x,y,z-coordinates of every atom in the structure). The program (or algorithm) assumes that the object's center of mass (CM) is already at the origin.
 
There are various approaches and algorithms out there. However, it is important, for me, the be able to code my own version. I have to calculate the moment of inertia for several thousand structures.
 
I have settled on finding the eigenvalues of a 3x3 matrix (the inertia tensor matrix). It is quite simple for me to sit down with pencil-and-paper and calculate these values (or use my TI-86 or Mathematica, etc.), however, I need to teach the computer to do this for me. It appears that this is not as simple as I had originally thought.
 
My approach is use random iterations that will converge upon the approximate eigenvalues for a given 3x3 matrix.
 
Below are the steps I have taken so far:
 
First, assume the following 3x3 matrix:
| I11 I12 I13 |
| I21 I22 I23 |
| I31 I32 I33 |
 
1. For each atom/point in the structure/object, calculate and sum the products, etc. of the matrix (i.e. I11 += y^2 + z^2, etc.).
2. Calculate the determinate, cofactor, and trace for the resultant inertia tensor matrix (using the sum of the values from step #1).
3. Iterate randomly through, say one million calculations (any large value should work here), converging upon an approximate eigenvalue whose difference is less than some pre-defined value. Use the following equation during the iterations:
 
x = [ ( (2*r^3) - (t*r^2) + d ) / ( (3*r^2) - (2*t*r) + c ) ],
 
where r = random value at given iteration,
t = trace, d = determinate, and c = cofactor
 
The value of x will be calculated again and again until it converges upon an eigenvalue (i.e. its difference is less than a given value).
 
Okay. So, I have coded the process above and it works. It is extremely fast and its approximate eigenvalues (the three of them) are extremely close to the actual values (as calculated by Mathematica).
 
I know why each step is done and I understand the overall process (I should. I coded it.). Here is what I don't understand, and why I am posting here: Can anyone tell me what that equation above is for or does? I found that equation on a website (can't remember which). The website presented it as "the one to use" but gave no reason why. A simple plot of the iterative values generated by this equation (during a sample run of the algorithm) produces a standard convergence x,y-plot. The equation appears to be an integral of some sort (maybe a harmonic series?). I just can't seem to figure out how it was integrated or why.
 
Any help with the above will be much appreciated. Also, I apologise for the mess that my post is. I hope the above makes sense.
 
Thank you!
IP Logged
Icarus
wu::riddles Moderator
Uberpuzzler
*****



Boldly going where even angels fear to tread.

   


Gender: male
Posts: 4863
Re: Moment of Inertia by Random Iterations  
« Reply #1 on: Jul 5th, 2005, 7:46pm »
Quote Quote Modify Modify

I'm not really understanding what it is you are trying to do. You say that you are trying to calculate the "moments of inertia". This is just summing up the inertias as you have described to get the inertia tensor. The diagonal elements of this tensor are the moments of inertias about your three coordinate axes.  
 
Since you are talking about the eigenvalues, I assume your problem is really to identify the principle axes of your molecules, and moments about them, since the principle axes are the eigenvectors, and the principle moments are the associated eigenvalues.
 
But this is not a difficult task requiring an iterative technique! For a 3x3 matrix, determining the eigenvectors is straightforward, if not exactly nice to do by hand.  But your program should not find it too hard. Let A be your inertia tensor. Find the eigenvalues as roots of the 3rd-degree polynomial det(A - tI) = 0. For each eigenvalue t, find vector v = (v1, v2, v3) for which (A-tI)v = 0. (Assume first that v1 = 0, and solve for v2, v3. If this requires dividing by 0 at any time, or once you've solved for v2, you find two different values for v3, assume v1 = 1 instead. If t is a root of multiplicity > 1, assume in turn that (v1, v2) = (0,0), then (0,1), (1, 0), and finally (1,1) until for one of the assumptions, division by zero is not required, and all three equations give the same result for v3.
 
(Test for multiplicity of a root: Let p(t) = det(A-tI) = at^3 + bt^2 + ct + d, and suppose that t is a root. The multiplicity of t is 1 if p'(t) = 3at^2 + 2bt + c <>0.
Otherwise, the multiplicity is 2 if p''(t) = 6at + 2b <> 0. If both are zero, the multiplicity is 3.)
 
The properties of the inertia matrix insure that real eigenvalues and eigenvectors always exist.
 
I hope this helps. I really don't see any need to pull in approximation techniques for values that are fairly straightforward to calculate directly (even if they are a bit messy).
IP Logged

"Pi goes on and on and on ...
And e is just as cursed.
I wonder: Which is larger
When their digits are reversed? " - Anonymous
Icarus
wu::riddles Moderator
Uberpuzzler
*****



Boldly going where even angels fear to tread.

   


Gender: male
Posts: 4863
Re: Moment of Inertia by Random Iterations  
« Reply #2 on: Jul 5th, 2005, 8:14pm »
Quote Quote Modify Modify

Looking back over your post again - I'll try to more directly answer your question. One problem I have is what you mean by "the cofactor". An arbitrary 3x3 matrix has 9 cofactors, not 1. They are the determinants of the 2x2 matrices that result from removing a column and row from the original matrix. So when you say "c is the cofactor", I have no idea which cofactor you are refering to.
 
Despite this, it looks to me like you are attempting to find a root of the characteristic polynomial p(t)=det(A - tI) by the Newton-Raphson method. The formula you give is the iteration formula tn+1 = tn - p(tn)/p'(tn).
 
Once you have one root t0, you can divide p(t) by (t - t0), then find the remaining 2 roots by using the quadratic formula.
 
------------------------------------------------
 
One other thing: you say that you assume that the c.g. is at zero. Unless you have specifically taken steps to position your molecule so that its c.g. is zero, you cannot assume this!
IP Logged

"Pi goes on and on and on ...
And e is just as cursed.
I wonder: Which is larger
When their digits are reversed? " - Anonymous
Chrsitoph
Guest

Email

Re: Moment of Inertia by Random Iterations  
« Reply #3 on: Jul 6th, 2005, 10:43am »
Quote Quote Modify Modify Remove Remove

Icarus,
 
Thank you, very much indeed, for your help with this. You have certainly pointed me in the right direction. I am actually a little embarrassed that I missed some of the obvious you pointed out.
 
I re-read my post and it is obvious I wasn't clear enough about what I am trying to do here. I should point out that my little program already works and produces the results I need. I am just trying to figure out, more completely, why it works.
 
The part I am especially interested in figuring out was the iterative process of finding the three eigenvalues of the 3x3 matrix. I am completely able to calculate these by hand (or with my TI-86 or Mathematica, etc.); I was just trying to "teach" the computer to do the same thing.
 
From what I have read, on various websites, there is no exact method for teaching a computer to find these eigenvalues (see http://en.wikipedia.org/wiki/Eigenvalue_algorithm).
 
Anyway, I am happy to know that my little program works, is very fast, and produces nearly exact eigenvalues.
 
When I wrote that I am "assuming" my rigid-body (or structure) has its center of mass (CM) at the origin, I meant that I actually _know_ that it is. I have written another little program that first shifts my x,y,z-coordinates relative to the CM before I run my "inertia" program.
 
The next step (after the CM program/function) is the following:
FOREACH ATOM (or "point" in the structure):
// Moments of inertia (sum them):
I11 += (y^2+z^2);
I22 += (x^2+z^2);
I33 += (x^2+y^2);
// Products of inertia (sum them):
I12 = I21 += x*y;
I13 = I31 += x*z;
I23 = I32 += y*z;
END
 
t = trace = I11 + I22 + I33 // sum of main diagonal
c = cofactor = (I11*I22) - (I12*I21) + (I11*I33) - (I13*I31) + (I22*I33) - (I23*I32);
d = determinate = (I11*I22*I33) - (I11*I32*I23) + (I12*I31*I23) - (I12*I21*I33) + (I13*I21*I32) - (I13*I31*I22);
 
// ITERATIVE LOOP
while( (size < 3) && (count < 250000) ) {
   x = rand(1000000); // arbitrarly chosen value
   oldx = x - 500000; // for +/- values
   dx = 100; // set initial difference high
 
   while( dx > 0.01 ) {
      newx = ( (2*oldx^3) - (t*oldx^2) + d ) / ( (3*oldx^2) - (2*t*oldx) + c );
      dx = oldx - newx;
 
      if( dx < 0 ) { dx = -1*dx; } // make positive again
    }
 // Setup a hashing array to capture the unique values (there should only be three)
 
 oldx = newx;
}
 
The above is _only_ pseudocode and missing some minor stuff (in order to preserve clairty and space).
 
As you can see, I only count six cofactors (and I only mentioned one in my previous post as an example of one).
 
Now, thanks to your help, I can see the Newton Method here. My function, f(x) = x^3 - t*x^2, has as it derivative, f'(x) = 3*x^2 - 2*x*t. But where does the "extra" value of 2 come in (multiplied by f(x))? Also, why am I adding the deteriminate in the numerator and adding the cofactor in the denominator? And, more generally, why is my function: f(x) = x^3 - t*x^2? What is the significance of this function? Since this is a general algorithm (it works for _any_ 3D rigid-body), wouldn't that make this function a general one and thus have some general significance?
 
I would like to point out again that I fully understand all the math here (I took linear algebra as an undergrad) and have no problem with the prgramming aspect. I also understand the theory behind Newton's Method.
 
Finally, the entire reason why I am coding this program/algorithm is to answer a more general question about protein structures (my "rigid bodies"). The goal is to have the program:
1. Take two structures;
2. Calculate their individual moments of inertia;
3. Take the three resultant eigenvalues;
4. Take the ratio of the largest eigenvalue with the smallest one (max/min); and then
5. Take the ratio of these values for both structures (structureA_ratio / structureB_ratio).
This final, and single, value should represent the overall shapes of the two structures. That is, is one bigger and/or more spherical than the other, or is one smaller and/or more cylindrical than other, etc.?
IP Logged
Icarus
wu::riddles Moderator
Uberpuzzler
*****



Boldly going where even angels fear to tread.

   


Gender: male
Posts: 4863
Re: Moment of Inertia by Random Iterations  
« Reply #4 on: Jul 6th, 2005, 5:15pm »
Quote Quote Modify Modify

on Jul 6th, 2005, 10:43am, Christoph wrote:
I am actually a little embarrassed that I missed some of the obvious you pointed out.

 
Happens all the time. I was embarrassed after stating in my first post that the process didn't require iteration to realize you are using iteration in the Newton-Raphson method to locate a root of the characteristic polynomial. For a general cubic polynomial, this is not a bad way to find your first root.
 
Quote:
From what I have read, on various websites, there is no exact method for teaching a computer to find these eigenvalues

 
To find eigenvalues of an NxN matrix A, you just find the roots of the characteristic polynomial det(xI - A). Note that the characteristic polynomial is of degree N. What that site is refering to is the fact that for N>4, there is no formula like the quadratic formula that will tell you the roots of ANY polymial (not just characteristic polynomials). So it talks about other techniques people have used to try and find roots.
 
However, there are two things that should be noted:  
(1) the problem indicated is for N>=5. Your matrix has N = 3, so that does not apply.
(2) the site is mistaken: There are algorithms that will find the roots of any polynomial. It is a well-known result that there are no algorithms that involve only standard operations and fractional powers (which is what the author of that entry was thinking of when he said no formula exists). However, at least two methods have been found that use certain transcendental functions to find the roots (similar to the method I give below for cubics, but using more esoteric functions than trigonometric).
 
Quote:
When I wrote that I am "assuming" my rigid-body (or structure) has its center of mass (CM) at the origin, I meant that I actually _know_ that it is.

 
Good. I hoped that this was the case, but the way you described it, I wasn't sure. Since I'm not sure what you know, I thought I had better raise a flag just in case.
 
Quote:
t = trace = I11 + I22 + I33 // sum of main diagonal
c = cofactor = (I11*I22) - (I12*I21) + (I11*I33) - (I13*I31) + (I22*I33) - (I23*I32);
d = determinate = (I11*I22*I33) - (I11*I32*I23) + (I12*I31*I23) - (I12*I21*I33) + (I13*I21*I32) - (I13*I31*I22);
 
...
 
As you can see, I only count six cofactors (and I only mentioned one in my previous post as an example of one).

 
I have never heard the value c called the "cofactor" before. That term is normally defined to mean something different (and not the individual terms in that sum either - cofactors are the determinants of submatrices of your original matrix). Your value c is simply the coefficient of the x term in the characteristic polynomial.
 
Quote:
My function, f(x) = x^3 - t*x^2

 
No - your function, whose roots you are trying to find, is the characteristic polynomial f(x) = det(xI - A) (A is your inertia matrix, I is the identity matrix). When expanded,
f(x) = x3 - tx2 + cx - d.
f'(x) = 3x2 - 2tx + c
 
For the Newton-Raphson method, the formula for the next value of x is x - f(x)/f'(x) = (f'(x)x - f(x))/f'(x). The numerator simplifies to 2x3 - tx2 + d, so the whole iteration becomes:
 
newx = (2*oldx^3 - t*oldx^2 + d)/(3*oldx^2 - 2*t*oldx + c).
 
Which is exactly what you have.
 
(If you are not already doing so, once you have found the first root x1, I would strongly recommend dividing f(x) by (x - x1) and using the quadratic formula on the result to find the other two. This is much faster and more stable than attempting to find them by additional iterations.)
 
As for the significance of the characteristic polynomial (assuming your uncertainty was not just caused by your having the wrong function), recall that a matrix A is invertible if and only if det(A) is non-zero. Also recall that eigenvalues x and the corresponding eigenvectors v of A are the numbers and non-zero vectors which satisfy the equation Av = xv.
 
We can rewrite the last equation as xv - Av = 0, or (xI - A)v = 0. Such a non-zero vector v exists if and only if (xI - A) is not invertible (if it had an inverse B, then multiplying both sides by B would give v = B0 = 0). So an eigenvector exists for x if and only if det(xI - A) = 0. Therefore the eigenvalues of A are the roots of the polynomial det(xI - A).
 
And yes, the ratio of the largest to smallest eigenvalue of the inertia tensor about the CG is indeed a good measure, though not perfect, of how far from spherical an object is. (To see that it isn't perfect, note that a cube will have the same value as a sphere: 1.)
 
--------------------------------------------------
An alternative method to finding the roots of f(x) = x3 - tx2 + cx - d.
 
1) Make the substitution x = y + t/3. The y2 term will cancel out when simplified, leaving a polynomial: y3 - Ay + B.
 
2) Make the substitution y = 2*sqrt(A/3)*z, and divide out the resultant polynomial by 2(A/3)3/2. You will be left with a polynomial of the form:
4z3 - 3z - C.
 
You want to solve the equation 4z3 - 3z = C.
 
3) A formula from trigonometry: 4 cos3 t - 3 cos t = cos 3t.
Set t = (1/3)cos-1 C.
The 3 roots of the z polynomial are cos t, cos(t + 2[pi]/3), cos(t+4[pi]/3).
 
4) Back out your substitutions to find the corresponding values of x.
 
Note that if |C| > 1 or C is complex, then you will need complex values for t, and will end up with complex roots. But for the characteristic polynomial of an inertia tensor, you are guaranteed 3 real eigenvalues. So |C| <= 1 for what you are doing.
« Last Edit: Jul 6th, 2005, 5:20pm by Icarus » IP Logged

"Pi goes on and on and on ...
And e is just as cursed.
I wonder: Which is larger
When their digits are reversed? " - Anonymous
SWF
Uberpuzzler
*****





   


Posts: 879
Re: Moment of Inertia by Random Iterations  
« Reply #5 on: Jul 6th, 2005, 5:38pm »
Quote Quote Modify Modify

While I was typing, I see Icarus made much of what I wrote below obsolete, but I'll leave it posted anyway...
 
I agree with Icarus.  That is Newton's Method for solving the
characteristic equation for p(x)=0, which is in terms of t, c, and d:  
 
  p(x) = x3 -t*x2 + c*x +d
 
In the equation Icarus gave, x(n+1) = x(n) - p(x(n))/p'(x(n)), just get a common denominator (p') to get the form you are using.
 
That program looks inefficient, and using random starting points does not seem like a good idea.  What if the three principal values are similar but different? You might have to get pretty lucky to pick a random number that converges to the middle value. If you have some identical values you have to run through the whole thing. The factoring method works well after you get one solution. Once you know two eigenvalues, you can get the third easily, since the sum of the three equals t, and their product equals d.
 
I like to use Jacobi transformations for finding principal axes, but I also care about the eigenvectors. Basically, you just keep doing coordinate transformations until you diagonalize the matrix.
 
I admit that using a numerical method to solve the for the first root of a cubic is usually my preferred method because it can be programmed very concisely.
 
If you want an exact solution to p(x)=0, try this (sometimes it leads to numerical inaccuracy):
 
Let q = (t*t-3*c)/9    p = (9*t*c - 27d - 2*t3)/27
    w= sqrt(q3-(p/2)2)
    y= atan2( w, -p/2)
    xj = t/3 + 2*sqrt(q)*cos( (y+2*[pi]*j)/3 )
 
The three eigenvalues are obtained by using j= 0, 1, 2.
« Last Edit: Jul 6th, 2005, 5:40pm by SWF » IP Logged
Christoph
Guest

Email

Re: Moment of Inertia by Random Iterations  
« Reply #6 on: Jul 6th, 2005, 6:08pm »
Quote Quote Modify Modify Remove Remove

Perfect! Thanks to both of you, Icarus and SWF. That is exactly what I needed to know.
 
I agree that the program or algorithm I have written (stealing a bunch of information from everywhere) is nowhere near as efficient as it could be. Unfortunately, I am not really a programmer; I am a 1st-year graduate student in Biophysics. I have only taken one computer programming class in my life (Intro to C). Also, although I have had plenty of math, I am not a mathematician.
 
Now I am going to see if I can implement some of your suggestions and see if the program can be made faster. This would be nice, as I have to run this program on several thousand structures. If I get stuck, again, I hope you don't mind if I post here again?
 
All the best!
IP Logged
Icarus
wu::riddles Moderator
Uberpuzzler
*****



Boldly going where even angels fear to tread.

   


Gender: male
Posts: 4863
Re: Moment of Inertia by Random Iterations  
« Reply #7 on: Jul 7th, 2005, 5:19pm »
Quote Quote Modify Modify

on Jul 6th, 2005, 5:38pm, SWF wrote:
I admit that using a numerical method to solve the for the first root of a cubic is usually my preferred method because it can be programmed very concisely.

 
And as you indicated in stating your own alternative method, Newton's method tends to be more stable for finding a root. My own alternative was given not to suggest a better approach, but to demonstrate that you can program a computer to solve these directly (unlike what Christoph was thinking because of the Wikipedia entry), at least for degree 3.
 
But even the non-numeric methods still give only approximations. Those approximations are just hidden inside pre-programmed functions. And the computational load of finding the arctan and cosine may be greater than with Newton's method.
 
So, really, using Newton's method is probably the best way to find the first root. Particularly if you can make a good guess for your starting value. I suggest starting off with oldx=t/3 (this is the average of the three roots).
« Last Edit: Jul 7th, 2005, 5:20pm by Icarus » IP Logged

"Pi goes on and on and on ...
And e is just as cursed.
I wonder: Which is larger
When their digits are reversed? " - Anonymous
Pages: 1  Reply Reply Notify of replies Notify of replies Send Topic Send Topic Print Print

« Previous topic | Next topic »

Powered by YaBB 1 Gold - SP 1.4!
Forum software copyright © 2000-2004 Yet another Bulletin Board