c++ - 特征多项式的 Souriau 方法

Question

有谁知道求任何 n × n 矩阵的特征多项式的 Souriau 方法？我找到了第一个系数，很明显，但是我怎样才能找到其他系数呢？在我需要反转矩阵但我知道怎么做之后。

#include <iostream> 
#include <fstream>
using namespace std;

double trace(double a[5][5],int n){ 
        int i;
        double trace=0;
        for(i=0;i<n;i++) 
            trace+=a[i][i]; 
        return trace;
}
double prod(double a[5][5],double b[5][5],int n) {
    double c[5][5];
    int i,j,k;
    cout << "\nProd:\n"; 
    for(i=0;i<n;++i){ 
        for(j=0;j<n;++j){ 
            c[i][j]=0; 
            for(k=0;k<n;++k) 
                c[i][j]=c[i][j]+(a[i][k]*b[k][j]);
                cout << c[i][j] << " "; 
        } 
        cout << "\n"; 
    } 
    return c[i][j]; 
}
double theta(double a[5][5], int n){
    int i;
    double theta[5];
    theta[1]=-trace(a,n);
    for(i=0;i<n;i++)
        cout << "Theta[" << i+1 << "]=" << theta[i+1] << "\n";
    return theta[i+1];
}

int main(){ 
    ifstream f("a.txt"); 
    ifstream g("b.txt");
    double a[5][5],b[5][5];
    int i,j,n; 
    f >> n; 
    g >> n; 
    for(i=0;i<n;++i) 
        for(j=0;j<n;++j) 
            f >> a[i][j]; 
    cout << "Matrix A:"<<endl;
    for(i=0;i<n;++i){
        for(j=0;j<n;++j) 
            cout << a[i][j] << " ";
            cout << endl;
    }
    cout << endl;
    for(i=0;i<n;++i) 
        for(j=0;j<n;++j) 
            g >> b[i][j]; 
    cout << "Matrix B:" << endl;
    for(i=0;i<n;++i){ 
        for(j=0;j<n;++j) 
            cout << b[i][j] << " ";
            cout << endl;
    }

        cout << endl;
        cout << "Trace = ";
        cout << trace(a,n);
        cout << endl;
        prod(a,b,n);
        cout << endl;
        theta(a,n);
    }

Answer 1

取自https://math.stackexchange.com/a/405975/115115经过 Ｊ． Ｍ.

C=A;
for k=1,…,n

    if k>1
        C=A*(C+c[n−k+1]*I);

    c[n−k]=−tr(C)/k;

end for

如果你能读懂德语，在 https://de.wikipedia.org/wiki/Algorithmus_von_Faddejew-Leverrier 有一个带有扩展伪代码算法的维基页面。 (添加 2017: 或同样好的英文版 https://en.wikipedia.org/wiki/Faddeev%E2%80%93LeVerrier_algorithm )

如果您想直接计算逆矩阵，那么您必须像在 wiki 页面中那样，使用通过 C=AB 与上述矩阵 C 相关的矩阵 B。正如可以在 wiki 页面中看到的那样，这给出了一个稍微复杂的算法。但是，那么最后一个矩阵B满足AB=-c[0]*I，所以如果有逆矩阵，可以直接计算。