Eigenvalues and eigenvectors of a 3x3 matrix (principal stresses)[programming]

[xiss burg]

I need to compute the 3 eigenvalues and 3 eigenvectors of a symmetric 3x3 matrix, namely a stress tensor, computationaly (in C++). More specific details http://en.wikipedia.org/wiki/Principal_stress#Principal_stresses_and_stress_invariants". Basically 2 questions:

1. I am running into trouble in finding the 3 real roots. http://en.wikipedia.org/wiki/Cubic_equation#General formula of roots" we have the general formula for the roots. According to Wikipedia, again, "The characteristic equation has three real roots λ, i.e. not imaginary due to the symmetry of the stress tensor.", then I'm sure these huge formulas can be simplified, and it is indeed a hard work to try simplifying it. Have you ever seen a simplified formula for cubic equations where the discriminant is always greater than zero (then, three real roots).

2. After finding the eigenvalues how to find the eigenvectors? I read about general purpose methods like the Power Method but these are for nxn matrices. For a 3x3 matrix there must be a much simpler technique.


Thanks in advance.

 

[Cyosis]

You can find the eigenvector belonging to eigenvalue \lambda by solving the system of equations Tv=\lambda v. You will get three equations which you can solve for v_1,v_2 and v_3.

 

[xiss burg]

You can find the eigenvector belonging to eigenvalue \lambda by solving the system of equations Tv=\lambda v. You will get three equations which you can solve for v_1,v_2 and v_3.

I know how to do this with pencil and paper, or in my mind perhaps ^^. But how to do it in programming? I don't know methods to compute a solution space (not only a single solution as most usual methods do like Gauss Seidel, Conugate Gradients, etc)..

 

[jrosen13]

Mathematica would just give you the answer, Eigenvalues[matrix], done. To calculate yourself, sure get your cubic polynomial and figure out a way to solve it, numerically perhaps?

 

[Lord Crc]

There's probably an even easier way for 3x3 matrices, but Jacobi (Givens) rotation is fairly easy to code and should be quick. See http://scholar.lib.vt.edu/theses/available/etd-62597-173629/unrestricted/chapter5a.PDF and http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm.

edit: Actually this page is more describing than the other Wikipedia link imho: http://en.wikipedia.org/wiki/Jacobi_rotation

 

[xiss burg]

Mathematica would just give you the answer, Eigenvalues[matrix], done. To calculate yourself, sure get your cubic polynomial and figure out a way to solve it, numerically perhaps?

Actually, I have to implement it myself in C++ code. Since there's a formula for the exact solution of a cubic equation I think its better to use it then. Numerical methods have problems finding all roots, don't they?

There's probably an even easier way for 3x3 matrices, but Jacobi (Givens) rotation is fairly easy to code and should be quick. See http://scholar.lib.vt.edu/theses/available/etd-62597-173629/unrestricted/chapter5a.PDF and http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm.

edit: Actually this page is more describing than the other Wikipedia link imho: http://en.wikipedia.org/wiki/Jacobi_rotation

Uh I think using these methods is a waste for 3x3 matrices. After looking more thoroughly into this problem I got a few conclusions which I am still unsure about. The eigen space of λ (vector space generated by all eigenvectors of λ) is the null space of the matrix (A-λI), or equivalently the solution space of the system (A-λI)x=0. We know that the row space of a matrix is orthogonal to its null space, then we can compute the eigenvector(s) of an eigenvalue by verifying the linear independence of the rows of (A-λI). I could compute the cross product of two pairs of row vectors, like r0Xr1 and r1Xr2 , then
-if these two cross products are zero, the rows of (A-λI) are linearly dependent, then the eigenspace is a plane orthogonal to any of the rows of (A-λI)
-if one of the cross products is different zero, then the cross product is the eigenvector itself
-it should not be possible for both to be different zero because then the eigenspace would be zero dimensional

What you think eh?

Thanks

 

[Lord Crc]

Uh I think using these methods is a waste for 3x3 matrices.

Perhaps. All I know is that it would get you the eigenvalues and the eigenvectors in one go, and the code I got is shorter than the stable cubic equation solver I use (formatted C).

 

[xiss burg]

Now I see that its better to use a numerical method, specially because I don't need a very accurate approximation for my specific application. Then I will give a try to that Jacobi method. Do you guys think this is one of the best choices for symmetrical 3x3 matrices?

Thanks.

 

[jrosen13]

That was my feeling before, that numerically you can solve a cubic polynomial to whatever accuracy you desire, up to some fundamental limit of time cost perhaps, but that you should be able to get really close without too much trouble. However, it is possible to analytically solve it! e.g. http://en.wikipedia.org/wiki/Cubic_function

That should be really fast! And exact too!