u= 0.8576 0.6934 0.5623 1.0000
index = 0
k= 1001
修改 M0=0 m= 2.6820
u=
0.8577 0.6937 0.5624 1.0000
index = 1
k=
7
总结以上,幂法如下：
U [1 1 1 1] [1 2 3 4] [3 5 6 7]
m0 0.0001 0.001
大型稀疏矩阵。反幂法是计算海森伯格阵或三角阵的对应一个给定近似特
征值的特征向量的有效方ຫໍສະໝຸດ 之一。二．算法设计及流程图
1、幂法算法
（1）取初始向量 u (0) （例如取 u (0) =(1,1,…1) T ）,置精度要求 ，置 k=1.
（2）计算
v (k ) =Au (k 1) ,m =max(v (k ) ), u (k ) = v (k ) / m
index 1 0 1 1 1 1 1 0 1
k 49 1001 10 9 7 9 7 1001 7
反幂法结果显示：在 m0 为 0 时
M0=0.001 U=[1 1 1 1]
M0=0.1 u=[1 1 1 1]
M0=0 u=[1 3 5 7]
M0=0.1 u=[1 3 5 7]
M0=0.5 u=[1 3 5 7]
的第一式改为求解线性方程组
A v (k ) = u (k 1)
（3）
但由于在反幂法中，每一步迭代都需求解线性方程组（3）式，迭代做了大量的
重复计算，为了节省工作量，可事先把矩阵 A 作 LU 分解，即 A=LU
所以线性方程组（3）改为
Ly (k ) =u (k 1) ,Uv (k ) =y (k)
四、算法程序设计代码
求
3.设计程序并进行计算；
4.对结果进行解释说明；
对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算 法设计和流程图，理论依据，程序及结果进行阐述该问题。
一．问题的分析：
求 n 阶方阵 A 的特征值和特征向量，是实际计算中常常碰到的问题，如：
采
机械、结构或电磁振动中的固有值问题等。对于 n 阶矩阵 A，若存在数 和
n 维向量 x 满足 用
Ax= x
方
则称 为矩阵 A 的特征值，x 为相应的特征向量。
法 由高等代数知识可知，特征值是代数方程
（1）
及
|
I-A|=
n
+a 1
n1 +…+a n1
+a n
=0
（2）
结
的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需
果
求解方程（2）的根，就能得到特征值 ，再解齐次方程组
0.2675 0.5776 0.6344 1.3130
0.5776 1.1503 0.7641 0.1367
0.6344 0.7641 0.0257 0.4193
1.3130 0.1367 0.4193 1.2248
>> u=[1 1 1 1]'; >> [m,u,index,k]=pow(A,u) m=
幂法程序，在 matlab 中建立一个 M 文件并保存。
%pow.m function [m,u,index,k]=pow(A,u,ep,it_max) if nargin<4
it_max=1000; end if nargin<3
ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0; I=eye(n); T=A-m0*I; while k<=it_max
则按 A 1 的特征值绝对值的大小排序，有
| 1 |>| 1 |≥…≥| 1 |
n
n 1
1
对 A 1 实行幂法，就可得 A 1 的绝对值最大的特征值 1/ n 和相应的特征向量， 即 A 的绝对值最小的特征值和相应的特征向量。
由于用 A 1 代替 A 作幂法计算，因此该方法称为反幂法，反幂法的迭代格
幂法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特
征值。稍微修改该方法，也可以用来确定其他特征值。幂法的一个很有用的特
性是它不仅可以生成特征值，而且可以生成相应的特征向量。实际上，幂法经
常用来求通过其他方法确定的特征值的特征向量。
1、幂法的迭代格式与收敛性质
设 n 阶矩阵 A 的特征值 1 ， 2 ,…, n 是按绝对值大小编号的， x i (i=1,2,…,n)为对应 i 的特征向量，且 1 为单根，即
0 0.0001 0.001
0 0.0001 0.001
0
m 2.6813 2.6814 2.6815 2.6813 2.6805 2.6814 2.6819 2.6914 2.692
u [0.8576 0.6934 0.5623 1.0000] [0.5876 0.6934 0.5623 1.0000] [0.8576 0.6935 0.5623 1.0000] [0.8576 0.6934 0.5623 1.0000] [0.8576 0.6934 0.5622 1.0000] [0.8576 0.6934 0.5623 1.0000] [0.8577 0.6937 0.5624 1.0000] [0.8576 0.6934 0.5623 1.0000] [0.8577 0.6937 0.5624 1.0000]
v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1)<ep;
index=1; break; end m=m+m0; m1=m; k=k+1; end
在 matlab 输入面板，输入 A=rand（4）；%产生一个 4 维随机矩阵
B=A+A’；
| 1 |>| 2 |≥…≥| n | 则计算最大特征值与特征向量的迭代格式为
v (k ) =Au (k 1) ,m k =max(v (k ) ), u (k ) = v (k ) / m k 其中 max(v (k ) )表示向量 v (k ) 绝对值的最大分量。
（1）
2、对于幂法的定理
按式（1）计算出 m 和 u (k ) 满足 k
B=A+A’；
u=[1 1 1 1]’;%设立初始向量 [m,u,index,k]=pow_inv(B，u,ep,it_max)%最多可省略 2 个参数 程序结束。 在 M 文件中可以通过改变 m0 的值改变原点位移，从而达到原点位移加速。
【结果显示】
%在 M0=1e-4 >>B=rand(4); >>A=B+B’ A=
lim
k
m
k
=
1
,
lim u (k ) = x1
k
max( x1 )
（二）反幂法算法的理论依据及推导
反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对 幂法的修改，可以给出更快的收敛性。 1、反幂法的迭代格式与收敛性质
设 A 是非奇异矩阵，则零不是特征值，并设特征值为 | 1 |≥| 2 |≥…≥| n1|>| n |
k=0;m1=0
v=invA*u [vmax,i]=max(abs(v))
m1=m;k=k+1
m=v(i);u=v/m
abs(m-m1)< 1e-6
index=1;break; 输出：m，u，index
输入 A； [m,u,inde x=]po结w(束A,1e -6)
三、算法的理论依据及其推导
（一）幂法算法的理论依据及推导
题
幂法和反幂法求矩阵特征值
目
具 随机产生一对称矩阵，对不同的原点位移和初值(至少取 3 个)分别使用幂
体 法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征
内 值及特征向量，并比较不同的原点位移和初值说明收敛。
容
1.认真读题，了解问题的数学原形；
要
2.选择合适问题求解的数值计算方法；
7
修改 M0=0 m=
2.6814
u=
0.8576 0.6934 0.5623 1.0000 index =
1 k=
9
修改 U=[3 5 6 7] 修改 M0=1e-4
m= 2.6819
u= 0.8577 0.6937 0.5624 1.0000
index = 1
k= 7
修改 M0=1e-3 m= 2.6814
index = 1
k= 10
修改 U=[1 2 3 4] 修改 M0=1e-4 m=
2.6813
u= 0.8576 0.6934 0.5623 1.0000
index = 1
k= 9
修改 M0=1e-3 m=
2.6805
u=
0.8576 0.6934 0.5622 1.0000
index =
1
k=
（3）解线性方程组
Ly (k ) =u (k 1) ,Uv (k ) =y (k )
（4）计算
m =max(v (k ) ), u (k ) = v (k ) / m
k
k
（5）若|m k =m k1 |< ，则停止计算（1/m k 作为绝对值最小特征值 n ，u (k ) 作为
相应的特征向量）;否则置 k=k+1，转（3）.
说
（ I-A）x=0
（3）
明
的解，就可得到相应的特征向量。
上述方法对于 n 很小时是可以的。但当 n 稍大时，计算工作量将以惊
人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程，