反幂法求矩阵特征值.docx

问题描述
用幂法与反幂法求解矩阵特征值
求n阶方阵A的特征值和特征向量，是实际计算中常常碰到的问题，如：机械、结构或电磁振动中的固有值问题等。对于n阶矩阵A，若存在数和n维向量x满足
Ax=x（1）
则称为矩阵A的特征值，x为相应的特征向量。
由线性代数知识可知，特征值是代数方程
|I-A|=+a+…+a+a=0（2）
的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需求解方程（2）的根，就能得到特征值，再解齐次方程组
（I-A）x=0（3）
的解，就可得到相应的特征向量。
上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程，自然就不必说求解方程（2）与（3）的困难了。幂法与反幂法是一种计算矩阵主特征值及对应特征向量的迭代方法,特别是用于大型稀疏矩阵。
这里用幂法与反幂法求解带状稀疏矩阵A[501][501]的特征值。

算法设计
1．幂法
（1）取初始向量u（例如取u=(1,1,…1)）,置精度要求，置k=1.
（2）计算
v=Au,m=max(v),u=v/m
（3）若|m-m|<，则停止计算（m作为绝对值最大特征值，u作为相应的特征向量）否则置k=k+1，转（2）
2.反幂法
（1）取初始向量u（例如取u=(1,1,…1)）,置精度要求，置k=1.
（2）对A作LU分解，即A=LU
（3）解线性方程组Ly=u,Uv=y
（4）计算
m=max(v),u=v/m
（5）若|m-m|<，则停止计算（1/m作为绝对值最小特征值，u作为相应的特征向量）;否则置k=k+1，转（3）.

程序框图
主程序

开始
设置二维数组Q[5][501]存放矩阵A的带内元素
平移求最大特征值和最小特征值
结束
fanmifa(A)反幂法求按模最小特征值
mifa(A)幂法求按模最大特征值
计算u[1]到u[39]

子程序
(1).幂法迭代程序框图

开始
计算幂法初始向量
fabs(b-c)/fabs(b)<=1e-12？

Return(b)并
printf("幂法成功！")
返回
>
b=sgn(h)*u[l];
找到Max(fabs(u[i]))
迭代求u[i]=u[i]+A[i-j+2][j]*y[j];
y[i]=u[i]/fabs(h);
c=b;

(2).反幂法迭代程序框图
开始
计算反幂法初始向量

Return(1/b)并
printf("反幂法成功")
返回
>
LU分解前保存A[i][j]和y[i]的值
c=b;
LU分解
fabs(b-c)/fabs(b)<=1e-12？
解线性方程组求出u[i]
迭代求b=b+y[i]*u[i];

结果显示

计算结果如下：
矩阵的按模最大特征值为:-1.070011361487e+001
矩阵的按模最小特征值为:-5.557910794230e-003
矩阵最大的特征值为:9.724634101479e+000
矩阵最小的特征值为:-1.070011361487e+001

与各最接近的（用表示）的值如下：
v[1]=-1.018293403315e+001u[1]=-1.018949492196e+001
v[2]=-9.585707425068e+000u[2]=-9.678876229054e+000
v[3]=-9.172672423928e+000u[3]=-9.168257536145e+000
v[4]=-8.652284007898e+000u[4]=-8.657638843237e+000
v[5]=-8.093483808675e+000u[5]=-8.147020150328e+000
v[6]=-7.659405407692e+000u[6]=-7.636401457419e+000
v[7]=-7.119684648691e+000u[7]=-7.125782764510e+000
v[8]=-6.611764339397e+000u[8]=-6.615164071601e+000
v[9]=-6.066103226595e+000u[9]=-6.104545378693e+000
v[10]=-5.585101052628e+000u[10]=-5.593926685784e+000
v[11]=-5.114083529812e+000u[11]=-5.083307992875e+000
v[12]=-4.578872176865e+000u[12]=-4.572689299966e+000
v[13]=-4.096470926260e+000u[13]=-4.062070607058e+000
v[14]=-3.5542