幂法和反幂法求矩阵特征值课程设计-学位论文.doc

题目
幂法和反幂法求矩阵特征值
具
体
内
容

随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量，并比较不同的原点位移和初值说明收敛。

要求
1.认真读题，了解问题的数学原形；
2.选择合适问题求解的数值计算方法；
3.设计程序并进行计算；
4.对结果进行解释说明；
采用方法
及结果
说明

对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算法设计和流程图，理论依据，程序及结果进行阐述该问题。
一．问题的分析：
求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）的困难了。幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别是用于大型稀疏矩阵。反幂法是计算海森伯格阵或三角阵的对应一个给定近似特征值的特征向量的有效方法之一。

二．算法设计及流程图
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）.

幂法流程图：

开始
输入A；[m,u,index]
=pow(A,1e-6)
k=0;m1=0
v=A*u
[vmax,i]=max(abs(v))
m=v(i);u=v/m
abs(m-m1)< 1e-6
index=1;break;
输出：m，u，index
结束
m1=m;k=k+1

反幂法流程图

开始
输入A；[m ,u,index]
=pow_inv(A,1e-6)
k=0;m1=0
v=invA*u
[vmax,i]=max(abs(v))
m=v(i);u=v/m
abs(m-m1)< 1e-6
index=1;break
输出：m，u，index
结束
m1=m;k=k+1
输入A；[m,u,index]
=pow(A,1e-6)

三、算法的理论依据及其推导
（一）幂法算法的理论依据及推导
幂法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特征值。稍微修改该方法，也可以用来确定其他特征值。幂法的一个很有用的特性是它不仅可以生成特征值，而且可以生成相应的特征向量。实际上，幂法经常用来求通过其他方法确定的特征值的特征向量。
1、幂法的迭代格式与收敛性质
设n阶矩阵A的特征值，,…,是按绝对值大小编号的，x(i=1,2,…,n)为对应的特征向量，且为单根，即
||>||≥…≥||
则计算最大特征值与特征向量的迭代格式为
v=Au,m=max(v), u= v/ m （1）
其中max(v)表示向量v绝对值的最大分量。
2、对于幂法的定理
按式（1）计算出m和u满足
m=, u=
（二）反幂法算法的理论依据及推导
反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对幂法的修改，可以给出更快的收敛性。
1、反幂法的迭代格式与收敛性质
设A是非奇异矩阵，则零不是特征值，并设特征值为
||≥||≥…≥||>||
则按A的特征值绝对值的大小排序，有
||>||≥…≥||
对A实行幂法，就可得A的绝对值最大的特征值1/和相应的特征向量，即A的绝对值最小的特征值和相应的特征向量。

由于用A代替A作幂法计算，因此该方法称为反幂法，反幂法的迭代格式为 v= Au,m=max(v), u= v/ m （2）
2、对于反幂法的定理
按式（2）计算出的m和u满足：
m=, u=
在式（2）中，需要用到A，这给计算带来很大的不方便，因此，把（2）式的第一式改为求解线性方程组
A v= u （3）
但由于在反幂法中，每一步迭代都需求解线性方程组（3）式，迭代做了大量的重复计算，为了节省工作量，可事先把矩阵A作LU分解，即 A=LU
所以线性方程组（3）改为
Ly=u,Uv=y

四、算法程序设计代码
幂法程序，在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
v=T*u;
[vmax,i]=max(abs(v));
m=v(i);
u=v/m;
if abs(m-m1)产生一个4维随机矩阵
B=A+A’；
u=[1 1 1 1]’;%设立初始向量
[m,u,index,k]=pow(B，u,ep,it_max)%最多可省略2个参数
程序结束。
在M文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速。

反幂法程序设计代码：
在matlab中建立一个M文件并保存。
%pow_inv.m
function[m,u,index,k]=pow_inv(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;
invT=inv(T);
while k<=it_max
v=invT*u;
[vmax,i]=max(abs(v));
m=v(i);
u=v/m;
if abs(m-m1)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 =
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 =

2.6813

u =

0.8576
0.6934
0.5623
1.0000

index =

1

k =

49

修改M0=1e-3
m =

2.6814

u =

0.8576
0.6934
0.5623
1.0000

index =

0

k =
1001

修改M0=0 %此时为幂法
m =

2.6815

u =

0.8576
0.6935
0.5623
1.0000

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 =

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

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
m0
m
u
index
k
[1 1 1 1]
0.0001
2.6813
[0.8576 0.6934 0.5623 1.0000]
1
49
0.001
2.6814
[0.5876 0.6934 0.5623 1.0000]
0
1001
0
2.6815
[0.8576 0.6935 0.5623 1.0000]
1
10
[1 2 3 4]
0.0001
2.6813
[0.8576 0.6934 0.5623 1.0000]
1
9
0.001
2.6805
[0.8576 0.6934 0.5622 1.0000]
1
7
0
2.6814
[0.8576 0.6934 0.5623 1.0000]
1
9
[3 5 6 7]
0.0001
2.6819
[0.8577 0.6937 0.5624 1.0000]
1
7
0.001
2.6914
[0.8576 0.6934 0.5623 1.0000]
0
1001
0
2.692
[0.8577 06937 0.5624 1.0000]
1
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]

M0=0 u=[2 3 4 5]

M0=0.1 u=[2 3 4 5]

M0=0.7 u=[2 3 4 5]

综上，反幂法结果如下：
u
m0
m
u
index
k
[1 1 1 1]
0.1
0.3847
[-0.8996 1.0000 0.2726 -0.2364]
1
15

0.001
0.3847
[-0.8996 1.0000 0.2726 -0.2364]
1
16

0
0.3847
[-0.8996 1.0000 0.2726 -0.2364]
1
16
[1 3 5 7]
0.5
0.3847
[-0.8995 1.0000 0.2726 -0.2364]
1
27

0.1
0.3847
[-08996 1.0000 0.2726 -0.2364]
1
17

0
0.3847
[-0.8996 1.0000 0.2726 -0.2364]
1
20
[2 3 4 5]
0.7
0.7091
[-0.6962 -0.4497 0.2196 1.0000]
1
5

0.1
0.3847
[-0.8995 1.0000 0.2726 -0.2364]
1
17

0
0.3847
[-0.8996 1.0000 0.2726 -0.2364]
1
19

五、结果分析
采用幂法和反幂法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。
同时，原点位移m的选取也影响收敛的速度。但原点位移m0的适当选取依赖于对矩阵A的大致了解。

成员
1007024104辛志贤
1007024107张 容
1007024108罗言月