Changeset 19030

Timestamp:

01/21/15 10:23:41 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added eigen values and vector calculations for 2x2 symetric matrices

Location:

issm/trunk-jpl/src/c/shared/Matrix

Files:

• MatrixUtils.cpp (modified) (1 diff)
• matrix.h (modified) (1 diff)

issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

@@ -333 +333 @@
 
 }/*}}}*/
+void Matrix2x2Eigen(IssmDouble* plambda1,IssmDouble* plambda2,IssmDouble* pvx, IssmDouble* pvy,IssmDouble a11, IssmDouble a21,IssmDouble a22){/*{{{*/
+        /*From symetric matrix (a11,a21;a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
+
+        /*Output*/
+        IssmDouble lambda1,lambda2;
+        IssmDouble vx,vy;
+
+        /*To get the eigen values, we must solve the following equation:
+         *     | a11 - lambda    a21        |
+         * det |                            | = 0
+         *     | a21             a22-lambda |
+         *
+         * We have to solve the following polynom:
+         *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
+
+        /*Compute polynom determinant*/
+        IssmDouble b=-a11-a22;
+        IssmDouble delta=b*b - 4*(a11*a22-a21*a21);
+
+        /*Compute norm of M to avoid round off errors*/
+        IssmDouble normM=a11*a11 + a22*a22 + a21*a21;
+
+        /*1: normM too small: eigen values = 0*/
+        if(normM<1.e-30){
+                lambda1=0.;
+                lambda2=0.;
+                vx=1.;
+                vy=0.;
+        }
+        /*2: delta is small -> double root*/
+        else if (delta < 1.e-5*normM){
+                lambda1=-b/2.;
+                lambda2=-b/2.;
+                vx=1.;
+                vy=0.;
+        }
+        /*3: general case -> two roots*/
+        else{
+                delta   = sqrt(delta);
+                lambda1 = (-b-delta)/2.;
+                lambda2 = (-b+delta)/2.;
+
+                /*Now, one must find the eigen vectors. For that we use the following property of the inner product
+                 *    <Ax,y> = <x,tAy>
+                 * Here, M'(M-lambda*Id) is symmetrical, which gives:
+                 *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
+                 * And we have the following:
+                 *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
+                 * We have shown that
+                 *    Im(M') ⊥ Ker(M')
+                 *
+                 * To find the eigen vectors of M, we only have to find two vectors
+                 * of the image of M' and take their perpendicular as long as they are
+                 * not 0.
+                 * To do that, we take the images (1,0) and (0,1):
+                 *  x1 = (a11 - lambda)      x2 = a21
+                 *  y1 = a21                 y2 = (a22-lambda)
+                 *
+                 * We take the vector that has the larger norm and take its perpendicular.*/
+
+                IssmDouble norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
+                IssmDouble norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
+
+                if(norm2<norm1){
+                        norm1=sqrt(norm1);
+                        vx = - a21/norm1;
+                        vy = (a11-lambda1)/norm1;
+                }
+                else{
+                        norm2=sqrt(norm2);
+                        vx = - (a22-lambda1)/norm2;
+                        vy = a21/norm2;
+                }
+        }
+
+        /*Assign output*/
+        *plambda1 = lambda1;
+        *plambda2 = lambda2;
+        *pvx      = vx;
+        *pvy      = vy;
+
+
+}/*}}}*/
 
 void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){/*{{{*/

issm/trunk-jpl/src/c/shared/Matrix/matrix.h

@@ -14 +14 @@
 void Matrix2x2Invert(IssmDouble* Ainv, IssmDouble* A);
 void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A);
+void Matrix2x2Eigen(IssmDouble* plambda1,IssmDouble* plambda2,IssmDouble* pvx, IssmDouble* pvy,IssmDouble a11, IssmDouble a21,IssmDouble a22);
 
 void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A);