数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序.docx
《数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序.docx》由会员分享，可在线阅读，更多相关《数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序.docx（22页珍藏版）》请在冰豆网上搜索。
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序
矩阵的特征值与特征向量的计算
摘要
物理，力学，工程技术中的很多问题在数学上都归结于求矩阵特征值的问题，例如振动问题（桥梁的振动，机械的振动，电磁振动等）、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分，本文使用幂法和反幂法分别求矩阵的按模最大，按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法，它最大的优点是方法简单，对于稀疏矩阵比较合适，但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值，及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词：
矩阵；特征值；特征向量；冥法；反冥法
THECALCULATIONSOFEIGENVALUEANDEIGENVECTOROFMATRIX
ABSTRACT
Physics,mechanics,engineeringtechnologyinalotofproblemsinmathematicsareattributedtomatrixeigenvalueproblem,suchasvibration（vibrationofthebridge,mechanicalvibration,electromagneticvibration,etc.）inphysics,somecriticalvaluesdetermineproblemsandtheoreticalphysicsinsomeoftheproblems.Matrixeigenvaluecalculationisaveryimportantpartinmatrixcomputation.Inthispaper,weusethepowermethodandinversepowermethodtocalculatethemaximumofthematrix,accordingtotheminimumcharacteristicvectorandthecorrespondingcharacteristicvalue.
Powermethodisaniterativemethodtocalculatetheeigenvaluesofamatrix.Ithastheadvantagethatthemethodissimpleandsuitableforsparsematrices,butsometimestheconvergencerateisveryslow.Thebasicideaistotakeanon-zeroinitialvector.Constructavectorsequencefromthematrixofthematrix.Thentheeigenvaluesandeigenvectorsareobtainedbyusingtheconstructedvectorsequence.
Theinversepowermethodisusedtocalculatetheminimumfeaturevectorsandtheireigenvaluesofthematrix,andtocalculatetheeigenvaluesofthematrix.Inthispaper,weusetheinversepowermethodtocalculatetheminimumeigenvalueofamatrixanditscorrespondingeigenvalues.Thebasicideaofcalculatingtheminimumcharacteristicvectorofamatrixistotransformittothemaximumcharacteristicvectorofthemodulusoftheinversematrix.Then,accordingtothemodel,theminimumfeaturevectoroftheoriginalmatrixisintroduced.
Keywords:
Matrix；Eigenvalue；Eigenvector；Iterationmethods;
1引言1
2相关定理。
1
3符号说明2
4冥法及反冥法2
4.1冥法3
4.2反冥法8
5QR算法14
参考文献18
附录19
1引言
在本论文中，我们主要讨论矩阵的特征值和特征向量的计算，我们知道，有很多现实中的问题都可以用到矩阵特征值与特征向量计算的知识，比如，在物理、力学和工程技术方面有很多的应用，并且发挥着极其重要的作用.因为这些问题都可归结为求矩阵特征值的问题，具体到一些具体问题，如振动问题，物理中某些临界值的确定问题以及一些理论物理中的问题.
在本论文中，我们主要介绍求矩阵的特征值与特征向量的一些原理和方法，原理涉
及高得代数中矩阵的相关定理，方法主要介绍冥法及反冥法并利用MATLAB算法的程序
来求解相关问题，加以验证.
2相关定理
定理2.1如果i（i1,2,...,n）是矩阵A的特征值，则有
nn
1iaiitrA
i1i1
2odetA12n.
定理2.2设A与B为相似矩阵（BT1AT），则
1oA与B有相同的特征值；
2o若x是B的一个特征向量，则Tx是A的特征向量
定理2.3设A（aij）nn，则A的每一个特征值必属于下述某个圆盘之中：
aijaij（i1,2,...,n）.
j1
ji
定义2.1设A是n阶是对称矩阵，对于任意非零向量x，称R（x）（Ax,x）为对应于向量xx,x
的Rayleigh商.
定理2.4设ARnn为对称矩阵（其特征值次序记作12n，对应的特征向量
x1,x2,xn组成规范化正交组，即（x,xj）ij），则
3符号说明
A:
n阶矩阵
B:
n阶矩阵
I：
n阶单位阵
i（i1,2,...,n）:
矩阵特征值
x:
实数域上的n维向量
vi（i0,1,n）：
实数域上的n维向量
uk（k0,1,,n,）：
实属上的规范化向量
4冥法及反冥法
4.1冥法
幂法是一种计算矩阵ARnn的主特征值的一种迭代法，它最大优点是方法简单，
适合于计算大型稀疏矩阵的主特征值
设A（aij）Rnn，其特征值为i，对应特征向量为xi（i1,,n）,即
Axiixi（i1,,n）
4.1.1）
且｛xi,,xn｝线性无关.设A特征值满足：
（即1为强占优）
|1||2||n|
称｛vk｝为迭代向量.
面来分折1及x1与｛vk｝关系.
且设10）
且有
（4.1.3）
n
vkAvk1Av0iixi
i1
vk1k（1x12
（2）kx2n（n）kxn）
1（a1x1k）
由假设（4.1.1）式,则
kim（i）0（i2,,n）
lkimk0
且收敛速度由比值r||确定.且有
1
这说明，当k充分大时，有vk/1k1x1，或vk/1k越来越接近特征向量1x1.下面考虑主特征值1的计算.
用（vk）i表示vk的第i个分量，考虑相邻迭代向量的分量的比值.
从而是
说明相邻迭代向量分量的比值收敛到主特征1,且收敛速度由比值r|2|来度量,r越
1
小收敛越快,但r越小收敛越快,但r|2|1,而接近于1时,收敛可能很慢.
1
定理4.1
（1）设ARnnn个线性无关的特征向量：
（2）设A特征值满足
|1||2||n|
（3）幂法：
v00（且10）
vkAvk1（k1,2,）则
（1）lim（vk1）i1x1；
k（vk）i11
又设A有n个线性无关的特征向量，x1,x2,,xn,其中
Axi1xi（i1,,r）,Axiixi（ir1,,n）
对于任意初始向量
n
v0ixi（且i,,r不全为零）
i1
则由幂法有
rn
vkAKv01kixii（i）kxi
i1ir1i
r
1k（ixik）
i1
且有
r
limvkkixi,（设1,r不全为零）
k1ki1ii1r
（k0当k）
由此，当k充分大时，vk/1k接近于与1对应的特征向量的某个线性组合
应用幂法计算A的主特征值1及对应的特征向量时，如果11（或11），迭代向
量的各个不等于零的分量将随k而趋于无究（或趋于零），这样电算时就可能溢出
为此，就南非要将迭代向量加以规范化.
设有非零向量
其中maxv（）表示向量v绝对值最大的元素，即如果有草药（v）i0max（v）i,则
max（v）（v）i0
1in
其中i0为所有绝对值最大的分量中最小指标
显然有下面性性质：
设t为实数,vrn，则
max（tv）tmax（v）
在定理4.1条件下幂法可改进为：
任取初始向量u0v00（且10）.
迭代：
规范化：
v1Av0v1Au0，u1
max（v1）max（Av0）
A2v0v2A2v0
v2Au1,u22k
max（Av0）max（v2）max（Akv0）（4.1．6）
面考查{uk},{vk}与计算1及x1的关系.
n
由v0ixi
i1
（2）考查迭代向量序列:
vAKv0k1（1x1k）
kmax（Ak1v0）max（1k1（1x1k1））
1x1k
1
max（1x1k1）
max（1x1k）
于是，kmax（vk）11,（当k）
max（1x1k1
定理（改进幂法）
（1）设ARnn有n个线性无关特征向量;
（2）设A特征值满足
且Axiixi（i1,2,,n）
（3）{uk},{vk}由改进幂法得到（（4.1.7）式）,则有
且收敛速度由比值r|2|确定.
个子程序求矩阵按模最
1
实现幂法，每迭代一次主要是计算一次矩阵乘向量（Au），可编大特征值如下：
%这个函数用于使用幂法求矩阵特征向量和特征值%A--矩阵，v--初始向量，e--精度function[t,p]=pm（A,v,e）
u=v./max（abs（v））;%
old=0;%记录上一次迭代得到的特征值
while
（1）v=A*u;
u=v./max（abs（v））;
if（abs（max（v）-old）break;
end
old=max（v）;
end
p=u;
t=max（v）;
end
例1.为检验以上代码的正确性，我们使用以上代码计算以下矩阵的最大特征值和特征向量
1,3,3
A2,1,5
4,3,1
结果为：
例2.利用你所编制的子程序求如下矩阵（从60到70阶）
21
121A
121
12
按模最大、按模最小的特征值及对应特征向量。
解：
代码见附录，运行得到的结果如下：
以上仅给出特征值的计算结果。
特征向量见附录，这里给出70阶的特征向量：
[0.58
-0.94
1.00-0.810.54-0.29
0.13
-0.05
0.01
-0.00
0.00
-0.00
0.00-0.000.000.00
0.000.000.000.00
0.000.000.00
0.00
0.00
0.000.00
0.000.000.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.000.00
0.000.000.00
0.00
0.00
0.00
0.00
0.00
0.000.000.000.000.000.000.000.000.00-0.000.00-0.00
0.00-0.000.01-0.050.13-0.290.54-0.811.00-0.940.58]
4.2反冥法
（1）反幂法可用来计算矩阵按模最小的特征值及对应的特征向量
设ARnn为非厅异矩阵，A特征值满足
2
对应特征向量x1,x1,,xn为线性无关，则A1特征求值为
特征向量为x1,x2,xn.
因此计算A的按模最小的特征值n的部题就是计算A1按模最大的特征值部题对于A1应用幂法迭代（称为反幂法），可求矩阵A1的主特征值1/n.反幂法迭代公式：
文案大全
任取初始向量v0u00（且n0），
k1，2，
4.2.1）
vkAuk1kmax（vk）
ukvk/uk
其中迭代向量vk可通过解方程组求得
Avkuk1
如果ARnn有n个线性无关特征向量且A特征值满足:
1n1n0
则由反幂法（2.11）构造的向量序列{uk},{vk}满足
，现要求对应
|jp|《|ip|,（ji）
|jp|
|ip,（ij）
说明是（ApI）1的主特征值.jp
现对（ApI）1应用幂法得到反幂法计算公式：
取初始向量u0v00（且0）,
j
k1,2,,
1
vk（ApI）uk1
kmax（vk）（4.2.2）
ukvk/uk
与定理8证明类似,可得下述结果.
定理10
（1）设ARnn有n个线性无关特征向量即Axiixi（i1,,n）.
（2）取pj（为A特征值j一个近似值），设（ApI）1存在且
|jp|《～|ip|（ji）
则由反幂法迭代公式（2，12）构造向量序列{uk},{vk}满足：
xj
（a）uk（当k）
max（xj）
1
（b）max（v）（当k）
kkp
反幂法迭代公式中vk是以通过解方程组
（ApI）vkuk1
求得.为了节省计算量，可先将（ApI）进行三角分解.
P（ApI）LU
其中p为置换阵，于是每次迭代求vk相当于求解两个三角形方程组
Lykpuk1
Uvkyk
可按下述方法取v0u0，即选u0使
Uv1L1Pu0（1,,1）T
回代求解即求得v1.
反幂法计算公式：
1．分解计算
P（ApI）LU，且保存L,U及P信息
2．反幂法迭代
（1）Uv1（1,,1）T求v1
u1v1/1,1max（v1）
（2）k2,
1）LykPuk1求yk
Uvkyk求vk
2）kmax（vk）
3）ukvk/k
对于计算对称三对角阵，或计算Hessenberg阵对应于一个给定的近似特征值的特征向量，反幂法是一个有效方法.
使用Matlab编写一个使用反幂法求矩阵最小特征值和特征向量的程序如下：
function[s,y]=fpm（A,x0,eps）%s为按模最小特征值，y是对应特征向量文案大全
k=1;
r=0;%r相当于0?
y=x0./max（abs（x0））;%规范化初始向量
[L,U]=lu（A）;
z=L\y;x=U\z;u=max（x）;
s=1/u;%按模最小为A-1按模最大的倒数.
ifabs（u-r）return
end
%终止条件
%这两步保证取出来的按模最大特征值
%是原值，而非其绝对值。
whileabs（u-r）>epsk=k+1;
r=u;y=x./max（abs（x））;z=L\y;
x=U\z;
u=max（x）;
end[m,index]=max（abs（x））;s=1/x（index）;
end
同样，取一个矩阵进行测试：
1,3,3
A2,1,5
4,3,1
计算结果为：
例2.利用你所编制的子程序求如下矩阵（从60到70阶）
21
121
A
121
12按模最小的特征值及对应特征向量。
代码见附录，程序结果如下图：
同样只给出70阶时的特征值，具体结果见附录
[0.040.090.130.180.220.260.300.350.390.430.470.510.540.58
0.62
0.65
0.68
0.72
0.75
0.77
0.80
0.83
0.85
0.87
0.89
0.91
0.93
0.95
0.96
0.97
0.98
0.99
0.99
1.00
1.00
1.00
1.00
0.99
0.99
0.98
0.97
0.96
0.95
0.93
0.91
0.89
0.87
0.85
0.83
0.80
0.77
0.75
0.72
0.68
0.65
0.62
0.58
0.54
0.51
0.47
0.43
0.39
0.35
0.30
0.26
0.22
0.18
0.13
0.09
0.04]
参考文献
[1]姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：
高等教育出版社，2005：
1-202.
[2]王建卫，曲中水凌滨编著.MATLAB7.X程序设计[M].北京：
中国水利水电出版社，2007：
55-80.
[3]李庆扬，王能超，易大义编著.数值分析（第四版）[M].武汉：
华中科技大学出版
社，2006：
219-245.
附录
%这个函数用来生成老师要求记算的那个矩阵，n是指定阶数
functionA=createMatrix（n）
A=zeros（n）;%先全部初始化为0
fori=1:
n
forj=1:
n
if（i==j）
A（i,j）=2;%设置主对角线上的值为2
elseif（i==j-1||i==j+1）%设置主对角线傍边的两条斜线上的的值为-1
A（i,j）=-1;
end
end
end
end
%这个函数用于使用幂法求矩阵特征向量和特征值
%A--矩阵，v--初始向量，e--精度function[t,p]=pm（A,v,e）u=v./max（abs（v））;%old=0;%记录上一次迭代得到的特征值
while
（1）v=A*u;u=v./max（abs（v））;
if（abs（max（v）-old）break;
end
old=max（v）;
end
p=u;
t=max（v）;
end
%这个程序用于求60-60阶矩阵的特征值和特征向量
clc
clear
e=0.01;
fori=60:
70
A=createMatrix（i）;%生成要计算的矩阵
v=ones（i,1）;%生成初始微量
v
（1）=1;
[t,p]=pm（A,v,e）;%计算
fprintf（'%d阶特征值：
%f\n',i,t）;%输出特征值%以下三句代码为输出特征值和特征微量
%fprintf（'%d阶:
%f[',i,t）;
%fprintf（'%.2f',p）;
%fprintf（']\n'）;
end%使用反幂法求矩阵按模最小特征值
function[s,y]=fpm（A,x0,eps）%s为按模最小特征值，y是对应特征向量k=1;
r=0;%r相当于0?
y=x0./max（abs（x0））;%规范化初始向量
[L,U]=lu（A）;
z=L\y;x=U\z;u=max（x）;
s=1/u;%按模最小为A-1按模最大的倒数.
ifabs（u-r）return
end
whileabs（u-r）>eps
%终止条件.
k=k+1;
r=u;
y=x./max（abs（x））;
z=L\y;
x=U\z;u=max（x）;
60-60阶矩阵的特征值和特征向量
end
[m,index]=max（abs（x））;
s=1/x（index）;
end
%这两步保证取出来的按模最大特征值
%是原值，而非其绝对值。
%这个程序用于使用反幂法求
clc
clear
e=0.01;
fori=60:
70
A=createMatrix（i）;
v=ones（i,1）;
v
（1）=1;
[t,p]=fpm（A,v,e）;
%fprintf（'%d阶特征值：
%f\n',i,t）;
fprintf（'%d阶:
%f[',i,t）;
fprintf（
fprintf（
'%.2f',p）;
']\n'）;
end
使用幂法求矩阵最大特征值和特征向量结果：
60阶:
3.754011[0.58-0.941.00-0.810.54-0.290.13-0.050.01-0.00
升级会员 