数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序Word格式文档下载.docx
《数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序Word格式文档下载.docx(22页珍藏版)》请在冰豆网上搜索。
在本论文中,我们主要讨论矩阵的特征值和特征向量的计算,我们知道,有很多现实中的问题都可以用到矩阵特征值与特征向量计算的知识,比如,在物理、力学和工程技术方面有很多的应用,并且发挥着极其重要的作用.因为这些问题都可归结为求矩阵特征值的问题,具体到一些具体问题,如振动问题,物理中某些临界值的确定问题以及一些理论物理中的问题.
在本论文中,我们主要介绍求矩阵的特征值与特征向量的一些原理和方法,原理涉
及高得代数中矩阵的相关定理,方法主要介绍冥法及反冥法并利用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:
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越
小收敛越快,但r越小收敛越快,但r|2|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)
对于任意初始向量
v0ixi(且i,,r不全为零)
则由幂法有
rn
vkAKv01kixii(i)kxi
i1ir1i
r
1k(ixik)
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的关系.
由v0ixi
(2)考查迭代向量序列:
vAKv0k1(1x1k)
kmax(Ak1v0)max(1k1(1x1k1))
1x1k
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|确定.
个子程序求矩阵按模最
实现幂法,每迭代一次主要是计算一次矩阵乘向量(Au),可编大特征值如下:
%这个函数用于使用幂法求矩阵特征向量和特征值%A--矩阵,v--初始向量,e--精度function[t,p]=pm(A,v,e)
u=v./max(abs(v));
%
old=0;
%记录上一次迭代得到的特征值
while
(1)v=A*u;
if(abs(max(v)-old)<
e)
break;
end
old=max(v);
p=u;
t=max(v);
例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.000.000.00
0.000.000.000.00
0.000.000.00
0.000.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,,
vk(ApI)uk1
kmax(vk)(4.2.2)
与定理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)
(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)<
eps%判断第一次迭代后是否满足终止条件
return
%终止条件
%这两步保证取出来的按模最大特征值
%是原值,而非其绝对值。
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);
同样,取一个矩阵进行测试:
计算结果为:
A
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
1.00
0.58
0.54
0.51
0.47
0.43
0.39
0.35
0.30
0.26
0.22
0.18
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:
forj=1:
if(i==j)
A(i,j)=2;
%设置主对角线上的值为2
elseif(i==j-1||i==j+1)%设置主对角线傍边的两条斜线上的的值为-1
A(i,j)=-1;
%这个函数用于使用幂法求矩阵特征向量和特征值
%A--矩阵,v--初始向量,e--精度function[t,p]=pm(A,v,e)u=v./max(abs(v));
%old=0;
u=v./max(abs(v));
%这个程序用于求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);
%.2f'
p);
]\n'
);
end%使用反幂法求矩阵按模最小特征值
function[s,y]=fpm(A,x0,eps)%s为按模最小特征值,y是对应特征向量k=1;
eps
%终止条件.
k=k+1;
y=x./max(abs(x));
60-60阶矩阵的特征值和特征向量
[m,index]=max(abs(x));
s=1/x(index);
%这个程序用于使用反幂法求
[t,p]=fpm(A,v,e);
%d阶:
fprintf(
'
p);
);
使用幂法求矩阵最大特征值和特征向量结果:
60阶:
3.754011[0.58-0.941.00-0.810.54-0.290.13-0.050.01-0.00