1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
| double[] eigenvaluesOfSymmetric3x3( double[][] M ) {
double p = Math.pow(M[0][1],2)+Math.pow(M[0][2],2)+Math.pow(M[1][2],2);
if (p==0) { // M is diagonal.
return new double[] {M[0][0],M[1][1],M[2][2]};
}
double q = trace(M)/3;
p = Math.pow(M[0][0]-q,2) + Math.pow(M[1][1]-q,2) + Math.pow(M[2][2]-q,2) + 2*p;
p = Math.sqrt(p / 6);
double[][] B = { // B = (M - q*I);
{ (M[0][0]-q), M[0][1], M[0][2] },
{ M[1][0], (M[1][1]-q), M[1][2] },
{ M[2][0], M[2][1], (M[2][2]-q) }
};
double r = det(B) / (2*p*p*p);
double phi = 0;
if (r<=-1) {phi = Math.PI/3;} else if (r>=1) {phi=0;} else {phi=Math.acos(r)/3;}
double eig1 = q + 2 * p * Math.cos(phi);
double eig2 = q + 2 * p * Math.cos(phi + 2*Math.PI/3);
double eig3 = q + 2 * p * Math.cos(phi + 4*Math.PI/3);
return new double[] {eig1,eig2,eig3};
}
double trace( double[][] M ) {
return M[0][0]+M[1][1]+M[2][2];
}
double det( double[][] M ) { // Sarrus formula
return M[0][0]*M[1][1]*M[2][2] + M[0][1]*M[1][2]*M[2][0] + M[0][2]*M[1][0]*M[2][1]
- M[0][2]*M[1][1]*M[2][0] - M[0][1]*M[1][0]*M[2][2] - M[0][0]*M[1][2]*M[2][1];
} |
Partager