[5661] | 1 | |
---|
| 2 | #include <stdio.h> |
---|
| 3 | #include <math.h> |
---|
[5662] | 4 | #include "matrix.h" |
---|
[5661] | 5 | |
---|
[5662] | 6 | |
---|
| 7 | class Vector |
---|
[5661] | 8 | { |
---|
[5662] | 9 | public: |
---|
| 10 | Vector (float x, float y, float z) { this->x=x; this->y = y; this->z = z; }; |
---|
| 11 | float x, y, z; |
---|
| 12 | inline Vector cross (const Vector& v) const { return Vector(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x ); } |
---|
| 13 | }; |
---|
[5661] | 14 | |
---|
[5662] | 15 | void Matrix::eigVl(const Matrix& mat) |
---|
| 16 | { |
---|
| 17 | |
---|
[5661] | 18 | float eigValue[3]; |
---|
| 19 | float eigVc[9]; |
---|
| 20 | |
---|
| 21 | float a = 0; |
---|
| 22 | float b = 0; |
---|
| 23 | |
---|
| 24 | float c[3]; |
---|
| 25 | |
---|
| 26 | // c[0] is the determinante of mat |
---|
[5662] | 27 | c[0] = this->m11 * this->m22 * this->m33 + |
---|
| 28 | 2* this->m12 * this->m13 * this->m23 - |
---|
| 29 | this->m11 * this->m23 * this->m23 - |
---|
| 30 | this->m22 * this->m13 * this->m13 - |
---|
| 31 | this->m33 * this->m12 * this->m12; |
---|
[5661] | 32 | |
---|
| 33 | // c[1] is the trace of a |
---|
[5662] | 34 | c[1] = this->m11 * this->m22 - |
---|
| 35 | this->m12 * this->m12 + |
---|
| 36 | this->m11 * this->m33 - |
---|
| 37 | this->m13 * this->m13 + |
---|
| 38 | this->m22 * this->m33 - |
---|
| 39 | this->m23 * this->m23; |
---|
[5661] | 40 | |
---|
| 41 | // c[2] is the sum of the diagonal elements |
---|
[5662] | 42 | c[2] = this->m11 + |
---|
| 43 | this->m22 + |
---|
| 44 | this->m33; |
---|
[5661] | 45 | |
---|
| 46 | |
---|
| 47 | // Computing the roots: |
---|
| 48 | a = (3.0*c[1] - c[2]*c[2]) / 3.0; |
---|
| 49 | b = (-2.0*c[2]*c[2]*c[2] + 9.0*c[1]*c[2] - 27.0*c[0]) / 27.0; |
---|
| 50 | |
---|
| 51 | float Q = b*b/4.0 + a*a*a/27.0; |
---|
| 52 | |
---|
[5662] | 53 | // 3 distinct Roots |
---|
[5661] | 54 | if (Q < 0) |
---|
| 55 | { |
---|
| 56 | printf("good\n"); |
---|
| 57 | float psi = atan2(sqrt(-Q), -b/2.0); |
---|
| 58 | float p = sqrt((b/2.0)*(b/2.0) - Q); |
---|
| 59 | |
---|
| 60 | eigValue[0] = c[2]/3.0 + 2 * pow(p, 1/3.0) * cos(psi/3.0); |
---|
| 61 | eigValue[1] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
| 62 | + sqrt(3.0) * sin(psi/3.0)); |
---|
| 63 | eigValue[2] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
| 64 | - sqrt(3.0) * sin(psi/3.0)); |
---|
| 65 | |
---|
| 66 | } |
---|
[5662] | 67 | // 2 Distinct Roots |
---|
[5661] | 68 | else if (Q == 0) |
---|
| 69 | { |
---|
| 70 | eigValue[0] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
| 71 | eigValue[1] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
| 72 | eigValue[2] = c[2]/3.0 + 2* pow(b/2.0, 1.0/3.0); |
---|
| 73 | } |
---|
[5662] | 74 | // 1 Root (not calculating anything.) |
---|
[5661] | 75 | else if (Q > 0) |
---|
| 76 | { |
---|
| 77 | printf("A is multiple of Identity matrix (lambda * I3))\n"); |
---|
| 78 | eigValue[0] = eigValue[1] = eigValue[2] = 1; |
---|
| 79 | } |
---|
| 80 | |
---|
[5662] | 81 | Matrix M; |
---|
[5661] | 82 | |
---|
[5662] | 83 | float u11, u12, u13, u22, u23, u33; |
---|
| 84 | |
---|
| 85 | this->debug(); |
---|
| 86 | |
---|
[5661] | 87 | printf("%f %f %f\n", eigValue[0], eigValue[1], eigValue[2]); |
---|
| 88 | |
---|
[5662] | 89 | } |
---|
[5661] | 90 | |
---|
| 91 | |
---|
[5662] | 92 | void Matrix::debug() const |
---|
| 93 | { |
---|
| 94 | printf("input: | %f | %f | %f |\n", this->m11, this->m12, this->m13 ); |
---|
| 95 | printf(" | %f | %f | %f |\n", this->m21, this->m22, this->m23 ); |
---|
| 96 | printf(" | %f | %f | %f |\n", this->m31, this->m32, this->m33 ); |
---|
[5661] | 97 | |
---|
| 98 | } |
---|
[5662] | 99 | |
---|