zoukankan      html  css  js  c++  java
  • 幂法求解矩阵特征值及特征向量

    幂法求解矩阵特征值及特征向量





    【算法原理】
    幂法是通过求矩阵特征向量来求出特征值的一种迭代法.其基本思想是:若我们求某个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
    【算法描述】

    设矩阵 有 个线性无关的特征向量,主特征值 满足

                       ,则 ,下式构造的向量序列 
                           
               
         
                  有     ,  

    【源代码】
    #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;
    }

    [TEST]
    原始数据
    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 )

    精度满足要求,程序设计及算法合理
  • 相关阅读:
    在TreeView控件节点中显示图片
    PAT 甲级 1146 Topological Order (25 分)
    PAT 甲级 1146 Topological Order (25 分)
    PAT 甲级 1145 Hashing
    PAT 甲级 1145 Hashing
    PAT 甲级 1144 The Missing Number (20 分)
    PAT 甲级 1144 The Missing Number (20 分)
    PAT 甲级 1151 LCA in a Binary Tree (30 分)
    PAT 甲级 1151 LCA in a Binary Tree (30 分)
    PAT 甲级 1149 Dangerous Goods Packaging
  • 原文地址:https://www.cnblogs.com/erwin/p/701184.html
Copyright © 2011-2022 走看看