幂法反幂法求解矩阵最大最小特征值和对应的特征向量.docx
《幂法反幂法求解矩阵最大最小特征值和对应的特征向量.docx》由会员分享，可在线阅读，更多相关《幂法反幂法求解矩阵最大最小特征值和对应的特征向量.docx（22页珍藏版）》请在冰豆网上搜索。
幂法反幂法求解矩阵最大最小特征值和对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量
—一.幂法
1.幕法简介：
当矩阵A满足一定条件时，在工程中可用幕法计算其主特征值（按模最大）
及其特征向量。
矩阵A需要满足的条件为：
（1）|「|「2卜…门|-0,'i为A的特征值
⑵存在n个线性无关的特征向量，设为X2，…，Xn
1.1计算过程：
n
对任意向量x（0），有x（0）-7rUi,〉i不全为0，则有
i4
x（k—Ax（k）Ak1x（0）
nn
k1k1
八Aam=aA5
「》（八2
/n、k1
（^）3nUn
1
i=1iT
1u1
可见，当丨厘丨越小时，收敛越快；且当k充分大时，有
x（k1）
■k1:
-
x（k⑴
[x（k）
二k：
U
11U1
x（k）
；（k^）
1，对应的特征向量即是x
2算法实现
（1）.输入矩阵A，初始向量X，误差限
最大迭代次数N
（k）
⑵.k「1,」
0;y（k）
max（abs（x（k））
⑶.计算xAy,…max（x）;
⑷若「-」：
；,输出］,y,否则,转（5）
（5）若kN,置k-k•1/：
--,转3,否则输出失败信息，停机.
3matlab程序代码
function[t,y]=lpowerA,xO,eps,N）%t为所求特征值，y是对应特征向量
k=1;
z=0;
y=x0./max（abs（x0））;x=A*y;
%z相当于■
%规范化初始向量
%迭代格式
b=max（x）;
%b相当于：
ifabs（z-b）%判断第一次迭代后是否满足要求
end
whileabs（z-b）>eps&&kz=b;
y=x./max（abs（x））;
x=A*y;
b=max（x）;
end
[m,index]=max（abs（x））;%这两步保证取出来的按模最大特征值
t=x（index）;end
%是原值，而非其绝对值。
4举例验证
选取一个矩阵A，代入程序，得到结果，并与eig（A）的得到结果比较，再计算A*y-t*y，验证y是否是对应的特征向量。
结果如下：
»A=[l10.5;110.25M5,0.25.2]
A=
1.0000
1.00000.5000
i.0000
1.00000.2500
0.5000
0.25002.0000
»
;eps=le-5:
X二20;
>>lx,y]=lpower（A,xO,eps,X）
ans=
2.5365
0.7482
0.6497
L0000
^0.0166
1.4801
2.5365
»A#y-t*y
ans-
1.0e-004+
-0.1603
-0,1684
0
结果正确，表明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与
eig（A）的得到结果比较，再计算A*y-t*y，验证y是否是对应的特征向量。
设
置初始向量为x0=ones（15,1）,结果显示如下»A二hilb（15）;
»x0=ones（15,1）:
】；eps-le^6:
»N=30;
>>Lt,y]-1power（A,xO,eps,X}
t二
L8459
y=
V-
»eig（A）
>>A*y-T*y
1.0000
ans=
ans=
0.6228
0.4706
0.3838
-0.00DO
0.0000
L0e-007*
a
0.3264
0.0000
-0“5516
0.0000
F643G
0.2851
0.0000
-0“6465
0.2537
0.0000
池6245
0.2290
0.0000
-①5953
0.2089
0.0000
-0.5650
0.1922
0.0000
75359
0.1780
0,0000
-0.5Q阴
0.1659
0,0004
-0”4834
0,0056
-0.4603
0.1554
a0572
-0.4390
0.1462
a.4266
-0.4195
0.1381
1.8459
-0.401b
jn.rnrs
可见，结果正确。
得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。
二•反幂法
1■反幕法简介及其理论
在工程计算中，可以利用反幕法计算矩阵按模最小特征值及其对应特征向量。
其基本理论如下，与幕法基本相同：
1
由Ax二’X=x=A‘（’X）,则A’xx，可知，A和A-1的特征值互为倒数，
Z
求A按模最小特征值即求A-1的按模最大特征值，取倒数即为A的按模最小特征值所以算法基本相同，区别就是在计算x（k1｝时，不是令x（k1）=Ay（k）,而是x（kA-1y（k）具体计算时，变换为
人严1）y（k）；对a做LU分解，来计算x（k1）
2■算法实现
（1）.输入矩阵A,初始向量x,误差限；,最大迭代次数N,
（2）•置k-1,，°“0,y,
max（abs（x））
（3）.作三角分解A二LU
（4）.解方程组LUx=y（Lz二y,Ux二z）,
（5）.：
—max（x）,■-:
一■-
（6）若|:
-咒
0卜；,输出〔,y,停机,否则转（7）,
转（4）;
max（abs（x））
（7）.若k：
：
N,置k・k1,'，屮
否则输出失败信息，停机.
3matlab程序代码
function[s,y]=invpower（A,xO,eps,n）%s为按模最小特征值，y是对应特征向量k=1;r=0;%r相当于-o
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&&kk=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）;%是原值，而非其绝对值。
end4举例验证
同幕法一样，选取一个矩阵A，代入程序，得到结果，并与eig（A）的得到结果比较，再计算A*y-t*y，验证y是否是对应的特征向量。
»A二[1J,0.5；1?
1,0.25；0.5,0.25,2]:
xO=[l;l：
j]:
eps=le~6;
n=30;
_s,y]=invpower（A,xOjeps,n）
s=
-0.0166
»eig（A）
ans二
-0.0166
1.4801
2.5365
ans=
y二
L0000
-0.9517
-0.1300
1.0e-016宕
-0.1388
-0.0694
-0.1908
可见结果正确，然后利用此程序计算15阶Hilb矩阵，eig（A）的得到结果比
较，再计算A*y-s*y，验证y是否是对应的特征向量。
设置初始向量为
x0=ones（15,1）,结果显示如下
»A=liilb（15）:
x0=ones（1571）;eps=le~6;n=30;
>>[s,y]=invpower（A?
eps^n）
-5.9673e-018
0.0000
-0.0000
0.0000
-0.0006
0.0050
-0.0245
0.0699
-0.0968
-0.0421
0.4497
-0.9084
L0000
-Q.6527
0.2379
-0.0375
»eig（A）
ans=
-0.0000
0.0000
0.0000
0.0000
0.0000
0+0000
0+0000
0,0000
0”0000
0.0000
0,0004
0.0056
0.0572
0.4266
1・8459
>>A*vs*y
已nw=
L0e-017*
0.3903
-0.5638
0.0000
-0.5208
-0.4741
63540
0.4537
-0.2746
CL9724
-0.6556
-0.6288
66618
0.0659
0.1419
可见，结果真确。
得到了量。
15阶Hilb矩阵的按模最大特征值和对应的特征向
三•计算条件数
矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond（A）=IIAA（-1）II,对应矩阵的3种范数，可以定义3种条件数。
函数
cond（A,1）、cond（A）或cond（Ainf）是判断矩阵病态与否的一种度量，条件数越大表明矩阵的病态程度越大•
这里我们选择矩阵的2范数，即cond（A）「分别为矩阵AtA的最大和最小特征值
而如果A为对称矩阵，如Hilb矩阵，AtA的最大最小特征值，分别为A的最大最小特征值的平方。
所以cond（A）为A的最大最小特征值得比值。
对于本例中的15阶Hilb矩阵来说，利用上面计算结果得其条件数（选择第二种条件数）
为：
3.0934e+017；这与直接利用cond（A）得到的结果：
2.5083e+017在同一
数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明阵是病态矩阵。
Hilb矩
四.Aitken商加速法
1.简介与原理
若：
ak收敛与a,且lim一二c=0;即：
ak：
线性收敛，k护ak-a
当k充分大时，有玉口玉心
ak-aakd1-a
、，、，（Xn雄—xn*）
yn2=Xn「
Xn42一2Xn*+Xn
1rr（ak*—ak）2“
2aak=-=Ok
ak2一2ak-1'ak
用ak逼近a,这种方法称为Aitken加速法.
同幕法和反幕法计算最大和最小特征值类似，如果计算最大特征值,
为x（k1）=Ay（k）；计算最小特征值时，迭代格式为
则迭代格式
x（k1^A-1y（k）,即y（k^Ax（k1）
2.算法实现
计算按模最大特征值算法如下：
（1）.输入A-（aj）,初始向量x,误差限；,最大迭代次数N,
（2）.置—7十0…"0'厂max（x）'
⑶.计算x=Ay,置max（x）二:
2
—（°-Of）
⑷.计算O。
-丄丄，
1*0
:
:
:
;,输出■,y停机，否则转（6）,
（5）.若卩-町
（6）.若k:
:
N,置-：
：
1=-；o,-'2=1，儿=八o,
k*“mx^xr弹⑶'
否则,输出失败信息，停机.
类似幕法和反幕法可以写出按模最小特征值算法，此处不再赘述。
3.matlab程序代码function[r,y]=aitken（A,xO,eps,n）%r按模最大特征值$为对应特征向量
k=1;
a0=0;%a相当于>0
a仁1;%a1相当于:
1
r0=1;%相当于2中的‘°
y=x0./max（abs（x0））;%规范化初始向量
x=A*y;
a2=max（abs（x））;%a2相当于：
2
r=a0-（a1-a0）A2/（a2-2*a1+a0）;%相当于'
if（a2-2*a1+a0）==0%若上式中分母为0,则迭代失败，返回
disp"初始向量迭代失败"
return;
end
ifabs（r-r0）end
whileabs（r-r0）>eps&&kk=k+1;
a0=a1;
a1=a2;
r0=r;
y=x./max（abs（x））;
x=A*y;%迭代格式
a2=max（abs（x））;
if（a2-2*a1+a0）==0%若分母为0,则迭代失败，返回return;
end
r=a0-（a1-a0）A2/（a2-2*a1+a0）;
[m,index]=max（abs（eig（A）））;%以下代码保证取出来的按模最大特征值aa=eig（A）;%是原值，而非其绝对值。
ifaa（index）>0||aa（index）==0
r=r;
else
r=-r;
end
end
end
类似可得按模最小特征值和特征向量的代码如下：
与上面类似，所不同的只是迭代格式不同•
function[r,y]=invaitken（A,xO,eps,n）
k=1;
a0=0;
a1=1;
r0=1;
y=x0./max（abs（x0））;
[L,U]=lu（A）;%迭代格式的不同
z=L\y;
x=U\z;
a2=max（abs（x））;
r=a0-（a1-a0）A2/（a2-2*a1+a0）;
if（a2-2*a1+a0）==0
disp"初始向量迭代失败"
return;
end
ifabs（r-r0）return
end
whileabs（r-r0）>eps&&kk=k+1;
a=b;
b=c;
r0=r;
y=x./max（abs（x））;
z=L\y;
x=U\z;
a2=max（abs（x））;
if（a2-2*a1+a0）==0
return;
end
r=a0-（a1-a0）A2/（a2-2*a1+a0）;
end
[m,index]=min（abs（eig（A）））;%以下代码保证取出来的按模最大特征值
aa=eig（A）;%是原值，而非其绝对值。
ifaa（index）>0||aa（index）==0
r=1/r;
else
r=-1/r;
end
end4.计算Hilb矩阵特征值
ones（15,1）
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为
结果如下，并将结果与幕法和反幕法得到结果比较
»A二hilb（15）;
xO=ones（1&I）;
»eps=le-6;
»n=30;
y_=aitken（AJxOjeps,n）
1.8459
L0000
0.6229
0.4706
0,3838
0.3264
0.2851
0,2537
0.2290
0.2089
0.1922
0.1780
0.1659
0.1554
0.1462
0.1381
同理，最小
这与幕法得到的特征值和特征向量一致，表明算法和代码正确;特征值结果如下：
>>[r,y]=invaitken&x0fepSjn）
5.9673e-018
0+0000-0,0000
0.0000-0,0006
0.0050-0.0245
0.0699
-0.0968
-Q.0421
0.4497
7+9084
L0000
-①6527
Q.2379-0.0375
这与反幕法得到的结果一致，表明结果正确
五，对称矩阵的Rayleigh商加速法
1.简介与原理
xAx
A为对称矩阵，设x=0,则称R（x）为关于A的Rayleigh商
xx
原理如下：
设Xj=0为•的特征向量，即AXj二iXi
用Xi左乘上式有:
（k）
y（）
x（k1）
即k））
_x（k）
max（x（k））
二Ay（k）
_~（k）\T/（k）、
（y（））（y（））
这称为Rayleigh商加速法
（k）
其中R（y（k）^1?
y（k）>
max（x（k））
2.算法实现
（1）.输入矩阵A,初始向量X,误差限；,最大迭代次数N,
（2）.置匕1,r°「0,y「
x
max（abs（x））
（3）xA*y,r「
（4）
x
max（abs（x））
.若|r_r°卜：
；,输出r,y,停机,否则转（5）,
（5）.若k:
:
N,置k「k1,r^r,y「
否则输出失败信息，停机.
3.Matlab程序代码
function[r,y]=rayleigh（A,xO,eps,n）%r是特征值，y是特征向量
k=1;r0=0;
y=xO./max（abs（xO））;
x=A*y;%迭代格式计算新的x
r=dot（y,x）/dot（y,y）;%Reyleigh商
ifabs（r-r0）return
end
whileabs（r-r0）>eps&&kk=k+1;
r0=r;
y=x./max（abs（x））;
x=A*y;
r=dot（y,x）/dot（y,y）;
end
end
类似得计算按模最小特征值的Rayleigh商加速法，如下：
function[r,y]=invrayleigh（A,x0,eps,n）
k=1;r0=0;
y=x0./max（abs（x0））;
[L,U]=lu（A）;%迭代格式不同
z=L\y;
x=U\z;
r=dot（y,x）/dot（y,y）;
ifabs（r-r0）return
end
whileabs（r-r0）>eps&&kk=k+1;
r0=r;
y=x./max（abs（x））;
z=L\y;
x=U\z;
r=dot（y,x）/dot（y,y）;
end
r=1/r;
end
ones（15,1）
4.计算Hilb矩阵特征值
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为
结果如下，并将结果与幕法和反幕法得到结果比较
»A=hilb（15）:
xO=ones（15>I）;eps=le~6;
n—30■
>>_r,y]=rayleigh（A>x0Peps^n）
8459
v二
L0000
0.6229
0.4706
0.3839
0,3265
0.2852
0.2538
0.2290
0.2089
0.1922
0.1781
0.1660
0.1555
0.1463
这与幕法得到结果一致，表明算法和代码正确
0.1381
同理，最小特征值如下:
»[rfy]=invray1eigh（A,xO,eps,n）
-5.9673e-018
0.0000-0.0000
0.0000
-0.0006
0.0050
-0.0245
0.0699
-0.0968
-0.0421
0.4497
-0.9084
1.0000
-0.6527
0.2379
与反幕法得到结果一致，表明代码和算法正确
70375
根据企业发展战略的要求，有计划地对人力、资源进行合理配置，通过对企业中员工的招聘、培训、使用、考核、评价、激励、调整等一系列过程，调动员工地积极性，发挥员工地潜能，为
企业创造价值，确保企业战略目标的实现。
读书是一种感悟人生的艺术读杜甫的诗使人感悟人生的辛酸，读李白的诗使人领悟官场的腐败，读鲁迅的文章使人认清社会的黑暗，读巴金的文章使人感到未来的希望每一本书都是一个朋友，
教会我们如何去看待人生读书是人生的一门最不缺少的功课，阅读书籍，感悟人生，助我们走好人生的每一步
升级会员 