幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx
《幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx》由会员分享，可在线阅读，更多相关《幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx（21页珍藏版）》请在冰豆网上搜索。
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量
一.幂法
1.幂法简介：
当矩阵A满足一定条件时，在工程中可用幂法计算其主特征值（按模最大）及其特征向量。
矩阵A需要满足的条件为：
（1）
（2）存在n个线性无关的特征向量，设为
1.1计算过程：
不全为0，那么有
可见，当
越小时，收敛越快；且当k充分大时，有
，对应的特征向量即是
。
2算法实现
3matlab程序代码
function[t,y]=lpowerA,x0,eps,N）%t为所求特征值，y是对应特征向量
k=1;
z=0;%z相当于
y=x0./max（abs（x0））;%标准化初始向量
x=A*y;%迭代格式
b=max（x）;%b相当于
ifabs（z-b）t=max（x）;
return;
end
whileabs（z-b）>eps&&kk=k+1;
z=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是否是对应的特征向量。
结果如下：
结果正确，说明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与eig（A）的得到结果比拟，再计算A*y-t*y，验证y是否是对应的特征向量。
设置初始向量为x0=ones（15,1），结果显示如下
可见，结果正确。
得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。
二.反幂法
1.反幂法简介及其理论
在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。
其根本理论如下，与幂法根本一样：
，可知，A和A-1的特征值互为倒数，求A按模最小特征值即求A-1的按模最大特征值,取倒数即为A的按模最小特征值所以算法根本一样，区别就是在计算
2.算法实现
3matlab程序代码
function[s,y]=invpower（A,x0,eps,n）%s为按模最小特征值，y是对应特征向量
k=1;r=0;%r相当于
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）;%是原值，而非其绝对值。
end
4举例验证
同幂法一样，选取一个矩阵A，代入程序，得到结果，并与eig（A）的得到结果比拟，再计算A*y-t*y，验证y是否是对应的特征向量。
可见结果正确，然后利用此程序计算15阶Hilb矩阵，eig（A）的得到结果比拟，再计算A*y-s*y，验证y是否是对应的特征向量。
设置初始向量为x0=ones（15,1），结果显示如下
可见，结果真确。
得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。
三.计算条件数
矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond（A）=‖A‖·‖A^（-1）‖，对应矩阵的3种范数，可以定义3种条件数。
函数cond（A,1）、cond（A）或cond（Ainf）是判断矩阵病态与否的一种度量，条件数越大说明矩阵的病态程度越大.
而如果A为对称矩阵，如Hilb矩阵，
的最大最小特征值，分别为A的最大最小特征值的平方。
所以cond（A）为A的最大最小特征值得比值。
对于本例中的15阶Hilb矩阵来说，利用上面计算结果得其条件数（选择第二种条件数〕为：
3.0934e+017；这与直接利用cond（A）得到的结果：
2.5083e+017在同一数量级，再次说明了上述算得得最大最小特征值的正确性，同时又说明Hilb矩阵是病态矩阵。
四.Aitken商加速法
1.简介与原理
同幂法和反幂法计算最大和最小特征值类似，如果计算最大特征值，那么迭代格式为
；计算最小特征值时，迭代格式为
。
2.算法实现
计算按模最大特征值算法如下：
类似幂法和反幂法可以写出按模最小特征值算法，此处不再赘述。
3.matlab程序代码
function[r,y]=aitken（A,x0,eps,n）%r按模最大特征值,y为对应特征向量
k=1;
a0=0;%a相当于
a1=1;%a1相当于
r0=1;%相当于2中的
y=x0./max（abs（x0））;%标准化初始向量
x=A*y;
a2=max（abs（x））;%a2相当于
r=a0-（a1-a0）^2/（a2-2*a1+a0）;%相当于
if（a2-2*a1+a0）==0%假设上式中分母为0，那么迭代失败，返回
disp"初始向量迭代失败"
return;
end
ifabs（r-r0）return
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）^2/（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,x0,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）^2/（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）^2/（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
end
4.计算Hilb矩阵特征值
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones（15,1），结果如下，并将结果与幂法和反幂法得到结果比拟
这与幂法得到的特征值和特征向量一致，说明算法和代码正确；同理，最小特征值结果如下：
这与反幂法得到的结果一致，说明结果正确。
五，对称矩阵的Rayleigh商加速法
1.简介与原理
原理如下：
2.算法实现
3.Matlab程序代码
function[r,y]=rayleigh（A,x0,eps,n）%r是特征值，y是特征向量
k=1;r0=0;
y=x0./max（abs（x0））;
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
4.计算Hilb矩阵特征值
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones（15,1），结果如下，并将结果与幂法和反幂法得到结果比拟
这与幂法得到结果一致，说明算法和代码正确。
同理，最小特征值如下：
与反幂法得到结果一致，说明代码和算法正确。
升级会员 