幂法求解矩阵特征值及特征向量
【算法原理】
幂法是通过求矩阵特征向量来求出特征值的一种迭代法.其基本思想是:若我们求某个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 )
精度满足要求，程序设计及算法合理