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

| 1-λ 2    |
|  3   4-λ | 

然後再展開來,以此為列,取得  λ  的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;
}

arrow
arrow
    全站熱搜

    Edison 發表在 痞客邦 留言(1) 人氣()