藍影

關於 eigenvalue 的問題，並沒有找到太多的相關文獻，
我只知道呆呆的列出下面的式子
 

然後再展開來，以此為列，取得  λ  的1元2次方程式後，去求 λ 的解。
在此先補充一些 1元 n 次方程式的觀念，
如果是 n 次方程式，那麼可能存在的實數根，最多為 n 個
所以，如果行列式為 n 階，展開來後即為 1元 n 次方程式，
所存在的實數根最多為 n 個。

在此，我使用的方法是 "暴力破解法"，
將要求輸入 λ 的上界(up)及下界(down)
以及每次 λ 之增量(step)，還有所允許之最小誤差值(error_small)

函數說明：
1.EIGENVALUE_BAD    使用暴力法求 eigenvalue (up.down.step.error_small 由使用者輸入)
2.EIGENVALUE_BAD2    使用暴力法求 eigenvalue (up.down.step.error_small 已先設定好)

至於 eigenvector，希望有人能提供如何解決重根之問題
這裡 eigenvector 不附上 code.

// =======================================
// author  : edison.shih.
// filename: matrix.cpp
// date    : 2010/2/7
// mail    : forget72@hotmail.com
// ** all copyright reverve **
// =======================================
// ==============================
// eigenvalue 暴力法
unsigned EIGENVALUE_BAD(double **a,
                        unsigned dim,
                        double low,
                        double up,
                        double step,
                        double error_small,
                        double *ans)
{
    unsigned x=0;
    unsigned cnt=0;
    double det1 = 0.0;
    double det2 = 0.0;
    double **tmp = (double**)malloc(sizeof(double*)*dim);
    for(x=0; x<dim; x++) tmp[x] = (double*)malloc(sizeof(double)*dim);
    CopyMatrix(tmp, a, dim, dim);
   
    for(x=0; x<dim; x++) ans[x] = 0.0;
   
   
    for(double i=low; i<=up; i=i+step){
        for(x=0; x<dim; x++) tmp[x][x] = a[x][x] - i;
        if(i==low) det1 = MDET(tmp, dim);
        det2 = MDET(tmp, dim);
        if(det2==0.0) {
            ans[cnt++] = i;
            if(cnt==dim) break;
        }
        else if((det2 < error_small || det1 < error_small) && (det1*det2 < 0.0))
        {       
            if(fabs(det1) < fabs(det2)) ans[cnt++] = i-step; // 上一個值比較好
            else ans[cnt++] = i; // 這次值比較好           
            if(cnt==dim) break;
        }
        det1 = det2;
    }
   
    for(x=0; x<dim; x++) free(tmp[x]);
    free(tmp);
    return cnt;
}

// eigenvalue 預設暴力法
unsigned EIGENVALUE_BAD2(double **a,
                         unsigned dim,
                         double *ans)
{
    unsigned x=0;
    unsigned cnt=0;
    double det1 = 0.0;
    double det2 = 0.0;
    const double low = -1000.0;
    const double up = 1000.0;
    const double step = 0.01;
    const double error_small = 0.0001;
   
   
    double **tmp = (double**)malloc(sizeof(double*)*dim);
    for(x=0; x<dim; x++) tmp[x] = (double*)malloc(sizeof(double)*dim);
    CopyMatrix(tmp, a, dim, dim);
   
    for(x=0; x<dim; x++) ans[x] = 0.0;
   
   
    for(double i=low; i<=up; i=i+step){
        for(x=0; x<dim; x++) tmp[x][x] = a[x][x] - i;
        if(i==low) det1 = MDET(tmp, dim);
        det2 = MDET(tmp, dim);
        if(det2==0.0) {
            ans[cnt++] = i;
            if(cnt==dim) break;
        }
        else if((det2 < error_small || det1 < error_small) && (det1*det2 < 0.0))
        {       
            if(fabs(det1) < fabs(det2)) ans[cnt++] = i-step; // 上一個值比較好
            else ans[cnt++] = i; // 這次值比較好           
            if(cnt==dim) break;
        }
        det1 = det2;
    }
   
    for(x=0; x<dim; x++) free(tmp[x]);
    free(tmp);
    return cnt;
}