プログラミング雑感  新JavaScript入門  JavaScript,Neo-Generation  掲示板  表紙
2. 多項式クラス(2)  4. 多項式クラス(4) 
プログラミング雑感
Written 7/15/01
3. 多項式クラス(3)
やりたかったことは、ある正方行列Aに対して、
    A - xE = 0
 
を満たす x (要するに固有値)を求める方程式を得ることです。
なんかの本にこれは難しいと書いてあったのですが、 ここまでくれば簡単なのでやってみましょ。
行列式
行列式は割り算を使わない素朴な方法で求めましょう。 というか、私がすぐに思いつく方法は、順列を使う方法と、再帰的に求める方法です。 ここでは比較的簡単と思われる再帰的に求める方法を採りましょう。
ピカピカの1年生が線形代数で習うと思うのですが、これで正しかったかな?

このように一つ小さい行列式にバラしてやることによって再帰的に求めます。

    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で実現しています。 いかにも効率が悪そうですが、いまさら気にしません。
結果
さて、すべて準備できたので計算してみましょう。
行列Aはa11,a12,...などを項にします。
    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 とか使い方がよく分からないのですが、これでいいのかな。
あと、対角成分で引き算して、行列式を求め、xのべきごとに表示します。
    //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項あるんですね。 これでいちおう固有値を求める方程式ができました。
first, prev, next, exit