Matlab实现隐式QR方法求解矩阵特征值特征向量.docx
《Matlab实现隐式QR方法求解矩阵特征值特征向量.docx》由会员分享，可在线阅读，更多相关《Matlab实现隐式QR方法求解矩阵特征值特征向量.docx（4页珍藏版）》请在冰点文库上搜索。
以下为4个matlab子程序和一个主程序：
实现了用隐式QR方法求方阵的特征值及特征向量
实质上是将一个矩阵Hessenberg化，其中要用到Householder变换，然后使用QR迭代得到特征值
最后用反幂法求解特征向量。
子程序1:
function[v,b]=house（x）
%householder变换H：
满足H=I-bvv'且Hx只有第一个分量非零
n=length（x）;
e=norm（x,inf）;
x=x/e;
s=x（2:
n）'*x（2:
n）;
v=zeros（n,1）;
v（2:
n）=x（2:
n）;
ifs==0
b=0;
else
a=sqrt（x
（1）^2+s）;
ifx
（1）<=0
v
（1）=x
（1）-a;
else
v
（1）=-s/（x
（1）+a）;
end
b=2*v
（1）^2/（s+v
（1）^2）;
v=v/v
（1）;
end
子程序2:
function[A,d]=Hessen（A）
%将A上Hessenberg化，次对角线以下用来存储Householder变换
%次对角线及以上用来存A的上Hessenberg化
n=length（A）;
d=zeros（n-1,1）;
fork=1:
n-2
[v,beta]=house（A（k+1:
n,k））;
A（k+1:
n,k:
n）=A（k+1:
n,k:
n）-（beta*v）*（v'*A（k+1:
n,k:
n））;
A（1:
n,k+1:
n）=A（1:
n,k+1:
n）-（A（1:
n,k+1:
n）*v）*（beta*v'）;
A（k+2:
n,k）=v（2:
n-k）;
d（k）=beta;
end
子程序3:
functionv=ipow（A,t）
%反幂法计算对应于A的特征值为t的特征向量v
n=length（A）;
[L,U,P]=lu（A-t*eye（n））;
v=ones（n,1）;
fori=2:
n
v（i）=v（i）-L（i,1:
i-1）*v（1:
i-1）;
end
v（n）=v（n）/U（n,n）;
fori=n-1:
-1:
1
v（i）=v（i）-U（i,i+1:
n）*v（i+1:
n）;
v（i）=v（i）/U（i,i）;
end
v=v/norm（v）;
子程序4:
functionH=QR2shift（H）
%双重步位移的QR方法，输入H并返回其一次双重步位移迭代后的矩阵
%调用格式:
循环H=QR2shift（H）
n=length（H）;
m=n-1;
s=H（m,m）+H（n,n）;
t=H（m,m）*H（n,n）-H（m,n）*H（n,m）;
x=H（1,1）*H（1,1）+H（1,2）*H（2,1）-s*H（1,1）+t;
y=H（2,1）*（H（1,1）+H（2,2）-s）;
z=H（2,1）*H（3,2）;
fork=0:
n-3
[v,b]=house（[x,y,z]'）;
q=max（1,k）;
H（k+1:
k+3,q:
n）=H（k+1:
k+3,q:
n）-（b*v）*（v'*H（k+1:
k+3,q:
n））;
r=min（k+4,n）;
H（1:
r,k+1:
k+3）=H（1:
r,k+1:
k+3）-（H（1:
r,k+1:
k+3）*v）*（b*v'）;
x=H（k+2,k+1）;
y=H（k+3,k+1）;
ifkz=H（k+4,k+1）;
end
end
[v,b]=house（[x,y]'）;
H（n-1:
n,n-2:
n）=H（n-1:
n,n-2:
n）-（b*v）*（v'*H（n-1:
n,n-2:
n））;
H（1:
n,n-1:
n）=H（1:
n,n-1:
n）-（H（1:
n,n-1:
n）*v）*（b*v'）;
主程序：
function[eigval,V]=QReig（A）
%计算A的特征值以及对应的特征向量
%调用格式[eigval,V]=QReig（A）,其中eigval为特征值，V的列为特征向量
u=1e-15;%机器精度
n=length（A）;
H=A;
H=Hessen（H）;%将A上Hessenberg化
fori=1:
n-2
H（i+2:
n,i）=zeros（n-i-1,1）;
end
while
（1）
i1=0;
i2=0;
fori=2:
n
if（abs（H（i,i-1））<=（abs（H（i,i））+abs（H（i-1,i-1）））*u）;
H（i,i-1）=0;
end
end
%将次对角线上足够小的元素置零
fori=n:
-1:
2
if（（i~=2&&H（i,i-1）~=0）&&H（i-1,i-2）~=0）
i2=i;
forj=i2:
-1:
2
if（H（j,j-1）==0）
i1=j;
H（i1:
i2,i1:
i2）=QR2shift（H（i1:
i2,i1:
i2））;
break;
end
end
if（i1==0）
H（1:
i2,1:
i2）=QR2shift（H（1:
i2,1:
i2））;
end
break;
end
end
if（i==2）
break;
end
end
%以下求复特征值
fori=2:
n
if（H（i,i-1）~=0）
p=H（i,i）+H（i-1,i-1）;
q=H（i,i）*H（i-1,i-1）-H（i,i-1）*H（i-1,i）;
delta=p^2-4*q;
H（i-1,i-1）=（p+sqrt（delta））/2;
H（i,i）=（p-sqrt（delta））/2;
end
end
eigval=diag（H）;
eigval=sort（eigval）;
%解得特征值后调用ipow解特征向量
V=zeros（n）;
fori=1:
n
V（1:
n,i）=ipow（A,eigval（i））;
end
升级会员 