A - xE = 0
を満たす x (要するに固有値)を求める方程式を得ることです。

このように一つ小さい行列式にバラしてやることによって再帰的に求めます。
typedef vector<vector<Polynomial> > Matrix; //行列
Polynomial determinant(const Matrix& m) {
int n;
n = m.size() - 1;
if(n == 0)
return m[0][0];
else {
Polynomial det = 0;
Matrix m2(n);
for(int i = 0; i < n; i++)
m2[i].resize(n);
for(int i = 0; i <= n; i++) {
//一番上の行とi列を除いた行列を作る
for(int j = 0; j < i; j++)
for(int k = 1; k <= n; k++)
m2[j][k-1] = m[j][k];
for(int j = i + 1; j <= n; j++)
for(int k = 1; k <= n; k++)
m2[j-1][k-1] = m[j][k];
//一番上の行の項と行列式をかけて足す
if(i % 2 == 0)
det += m[i][0] * determinant(m2);
else
det -= m[i][0] * determinant(m2);
}
return det;
}
}
行列はvectorを要素とするvectorで実現しています。
いかにも効率が悪そうですが、いまさら気にしません。
Matrix m(5); //matrix 5x5
for(int i = 0; i < 5; i++)
m[i].resize(5);
//m = (aij)
for(int i = 0; i < 5; i++) {
for(int j = 0; j < 5; j++) {
//aijという文字列を作って多項式にして各項の要素とする
ostringstream ost;
ost << "a" << i + 1 << j + 1;
m[i][j] = Polynomial(ost.str());
}
}
ostringstream とか使い方がよく分からないのですが、これでいいのかな。
//m = (aij - xδij)
for(int i = 0; i < 5; i++)
m[i][i] -= Polynomial("x");
Polynomial det = determinant(m);
map<int,Polynomial> c = det.getTerms("x");
map<int,Polynomial>::const_iterator p;
for(p = c.begin(); p != c.end(); ++p)
cout << "x^" << p->first << " : " << p->second << endl;
完成したファイルは、
polynomial.cpp
これと前回のヘッダファイルでコンパイルして、実行した結果が次のようです。
(以上略)
x^2 : a11 * a22 * a33 + a11 * a22 * a44 + a11 * a22 * a55 - a11 * a23 * a32
- a11 * a24 * a42 - a11 * a25 * a52 + a11 * a33 * a44 + a11 * a33 * a55 -
a11 * a34 * a43 - a11 * a35 * a53 + a11 * a44 * a55 - a11 * a45 * a54 - a1
2 * a21 * a33 - a12 * a21 * a44 - a12 * a21 * a55 + a12 * a23 * a31 + a12
* a24 * a41 + a12 * a25 * a51 + a13 * a21 * a32 - a13 * a22 * a31 - a13 *
a31 * a44 - a13 * a31 * a55 + a13 * a34 * a41 + a13 * a35 * a51 + a14 * a2
1 * a42 - a14 * a22 * a41 + a14 * a31 * a43 - a14 * a33 * a41 - a14 * a41
* a55 + a14 * a45 * a51 + a15 * a21 * a52 - a15 * a22 * a51 + a15 * a31 *
a53 - a15 * a33 * a51 + a15 * a41 * a54 - a15 * a44 * a51 + a22 * a33 * a4
4 + a22 * a33 * a55 - a22 * a34 * a43 - a22 * a35 * a53 + a22 * a44 * a55
- a22 * a45 * a54 - a23 * a32 * a44 - a23 * a32 * a55 + a23 * a34 * a42 +
a23 * a35 * a52 + a24 * a32 * a43 - a24 * a33 * a42 - a24 * a42 * a55 + a2
4 * a45 * a52 + a25 * a32 * a53 - a25 * a33 * a52 + a25 * a42 * a54 - a25
* a44 * a52 + a33 * a44 * a55 - a33 * a45 * a54 - a34 * a43 * a55 + a34 *
a45 * a53 + a35 * a43 * a54 - a35 * a44 * a53
x^3 : -a11 * a22 - a11 * a33 - a11 * a44 - a11 * a55 + a12 * a21 + a13 * a3
1 + a14 * a41 + a15 * a51 - a22 * a33 - a22 * a44 - a22 * a55 + a23 * a32
+ a24 * a42 + a25 * a52 - a33 * a44 - a33 * a55 + a34 * a43 + a35 * a53 -
a44 * a55 + a45 * a54
x^4 : a11 + a22 + a33 + a44 + a55
x^5 : -1
えっと、しゃれにならないですね。
よく考えたら、xの0乗の項は5!=120項あるんですね。
これでいちおう固有値を求める方程式ができました。