冪法求解矩陣特徵值及特徵向量_演算法

  冪法求解矩陣特徵值及特徵向量



【演算法原理】
冪法是通過求矩陣特徵向量來求出特徵值的一種迭代法.其基本思想是:若我們求某個n階方陣A的特徵值和特徵向量,先任取一個初始向量X(0),構造如下序列:
       X(0)  ,X(1)  =AX(0)  ,X(2)  =AX(1) ,…, X(K)  =AX(K+1)  ,…  ⑴
       當k增大時,序列的收斂情況與絕對值最大的特徵值有密切關係,分析這一序列的極限,即可求出按模最大的特徵值和特徵向量.
       假定矩陣A有n個線性無關的特徵向量.n個特徵值按模由大到小排列:
                         │λ1│>=│λ2│>=…>=│λn│              ⑵
    其相應的特徵向量為:
                         V1 ,V2 , …,Vn                                       ⑶
     它們構成n維空間的一組基.任取的初始向量X(0)由它們的線性組合給出
                         X(0)=a1V1+a2V2+…+anVn                         ⑷
     由此知,構造的向量序列有
          X(k)  =AX(k-1) = A2X(k-2) =…=AkX(0)  = a1λ1kV1+a2 λ2kV2+…+anλnkVn                    ⑸
   下面按模最大特徵值λ1是單根的情況討論:
      
 由此公式(5)可寫成
                   X(k) = λ1k (a1V1+a2 (λ2/λ1)kV2+…+an(λn/λ1)kVn  )       ⑹ 
   若a1≠0，由於|λi/λ1 |<1 (i≥2),故k充分大時，
                  X(k) = λ1k (a1V1+εk)
    其中εk為一可以忽略的小量，這說明X(k)與特徵向量V1相差一個常數因子，即使a1=0，由於計算過程的舍入誤差，必將引入在方向上的微小分量，這一分量隨著迭代過程的進展而逐漸成為主導，其收斂情況最終也將與相同。
特徵值按下屬方法求得：
            λ1 ≈Xj(k+1)/ Xj(k)                                        ⑺
其中Xj(k+1)， Xj(k)分別為X(k+1)，X(k)的第j各分量。
    實際計算時，為了避免計算過程中出現絕對值過大或過小的數參加運算，通常在每步迭代時，將向量“歸一化”即用的按模最大的分量         max  |Xj(k)| 1≤j≤n
去除X(k)的各個分量,得到歸一化的向量Y(k),並令X(k+1) = AY(k)


由此得到下列選代公式 :
          Y(k) = X(k)/║ X(k)║∞
           X(k+1) = AY(k)           k=0,1,2,…             ⑻
當k充分大時,或當║ X(k)- X(k+1)║<ε時,
          Y(k)≈V1
          max  |Xj(k)| ≈ λ1                                ⑼
          1≤j≤n
【演算法描述】

設矩陣 有 個線性無關的特徵向量，主特徵值 滿足

                   ，則 ，下式構造的向量序列 
                       
           
     
              有     ，  

【原始碼】

//////////////////////////////////////////////////////
//冪法求矩陣特徵值      // 
//author:zhiyong fang           //
//date:15/11/2004              //
////////////////////////////////////////////////////
#include <iostream.h>
#include <math.h>
#define N 3
void matrixx(double A[N][N],double x[N],double v[N])
{
    for(int i=0;i<N;i++)
          {
            v[i]=0;
            for(int j=0;j<N;j++)
                v[i]+=A[i][j]*x[j];
          }
    
}
double slove(double v[N])
{
    double max;
    for(int i=0;i<N-1;i++) max=v[i]>v[i+1]?v[i]:v[i+1];
    return max;
}
void main()    
{
    //data input
double A[N][N]={1.0,1.0,0.5,1.0,1.0,0.25,0.5,0.25,2.0};
    double x[N]={1,1,1};
    double v[N]={0,0,0};
    double u[N]={0,0,0};
    double p[N]={0,0,0};
    double e=1e-10,delta=1;
    int k=0;
    while(delta>=e)
    {
                for(int q=0;q<N;q++) p[q]=v[q];
        matrixx(A,x,v);
        for(int i=0;i<N;i++) u[i]=v[i]/(slove(v));
        delta=fabs(slove(v)-slove(p));
        k++;
        for(int l=0;l<N;l++) x[l]=u[l];
    }
    cout << "迭代次數" << k << endl;
    cout << "矩陣的特徵值" << slove(v) << endl;
    cout << "(" ;
    for(int i=0;i<N;i++) cout << u[i] << " " ;
    cout  << ")" << endl;
} 
[checkdata]
未經處理資料
double A[N][N]={1.0,1.0,0.5,1.0,1.0,0.25,0.5,0.25,2.0};
    double x[N]={1,1,1};
    double v[N]={0,0,0};
    double u[N]={0,0,0};
    double p[N]={0,0,0};
    double e=1e-10,delta=1;
輸出結果
迭代次數41
矩陣的特徵值2.53653
(0.748221 0.649661 1 )

精度滿足要求，程式設計及演算法合理