隐式QR法求实矩阵的全部特征值matlab实现.docx
《隐式QR法求实矩阵的全部特征值matlab实现.docx》由会员分享，可在线阅读，更多相关《隐式QR法求实矩阵的全部特征值matlab实现.docx（12页珍藏版）》请在冰豆网上搜索。
隐式QR法求实矩阵的全部特征值matlab实现
隐式QR法求实矩阵的全部特征值matlab实现
要求：
用matlab编写通用子程序，利用隐式QR法求实矩阵的全部特征值和特征向量。
思想：
隐式QR法实质上就是将一个矩阵Schur化，之后求解特征值就比较方便。
而隐式QR法还需要用到household变换，以及上hessenberg变换。
最后使用QR迭代，达到Schur化的结果。
步骤：
1.将矩阵A上hessenberg化（算法6.4.1），送而得到一个上hessenberg形矩阵H；
2.可约性判定，也就是判断次对角线元素是否非零，如果次对角线元素非零，则不可约。
3.Schur化，也就是通过QR迭代，将矩阵H变化成为某些次对角线元素变成0，同时还要满足，这些元素之间间隔最大为1，那么，所得到的最重的矩阵H就是一个Schur形矩阵。
4.假如两个等于0的次对角线元素间隔为0，那么该元素的上面一个元素，也就是H的对角线上的元素，即为其中一个特征值；假如两个等于0的次对角线元素间隔为1，那么在这两个元素之间就形成了一个2*2的矩阵，可以求解一个一元二次方程来得到两个共轭的特征值。
实验代码：
详见附录2
实验结果：
（代码相见附录2）
（i）设矩阵A如下：
求x=0.9,1.0,1.1时的特征值和特征向量。
X=0.9：
r是特征值，V是特征向量矩阵。
X=1：
r是特征值，V是特征向量矩阵。
X=1.1：
r是特征值，V是特征向量矩阵。
（ii）
求
的所有根。
附录2
隐式QR迭代：
主程序：
function[r,V]=SchurQR（A）
%向量r用来储存特征值
%Hessenberg分解:
[m,m]=size（A）;
fork=1:
m-2
[v,b]=house（A（k+1:
m,k））;
H1=eye（m-k）-b*v*v';
H2=eye（m）;
fori=k+1:
m
forj=k+1:
m
H2（i,j）=H1（i-k,j-k）;
end
end
ifk==1;
H=H2;
else
H=H*H2;
end
A（k+1:
m,k:
m）=H1*A（k+1:
m,k:
m）;
A（1:
m,k+1:
m）=A（1:
m,k+1:
m）*H1;
end
u=10e-5;
fori=2:
m;
ifabs（A（i,i-1））<=（abs（A（i,i））+abs（A（i-1,i-1）））*u;
A（i,i-1）=0;
end
end
%QR迭代：
H22=A;
x=Ifreducible（H22）;
whilex==1
H22=Francis（H22）;
x=Ifreducible（H22）;
end
[r,V]=EigValue（H22）;
子程序1：
function[r,V]=EigValue（A）%计算A的特征值，特征向量
[n,n]=size（A）;
r=zeros（1,n）;
y=zeros（1,n-1）;%y用来储存次对角线元素
fori=1:
n-1
y（i）=A（i+1,i）;
end
m=0;
fori=1:
n-1
ifabs（y（i）-0）<1e-5
m=m+1;
end
end
ifm==0
x=1;
else
z=zeros（1,m）;%z用来储存值为0的y向量的角标。
j=1;
i=1;
while（iifabs（y（i）-0）<1e-5
z（j）=i;
j=j+1;
end
i=i+1;
end
end
ifz
（1）==2%次对角线第一个等于0的元素的位置不同，需要2分类讨论
p=[1,A（1,1）+A（2,2）,A（1,1）*A（2,2）-A（1,2）*A（2,1）]
r（1:
2）=roots（p）;%求2*2矩阵的特征值
j=1;
whilejifz（j+1）-z（j）==1
r（z（j+1））=A（z（j+1）,z（j+1））;
end
if（z（j+1）-z（j）==2）
p=[1,-（A（z（j+1）-1,z（j+1）-1）+A（z（j+1）,z（j+1）））,A（z（j+1）-1,z（j+1）-1）*A（z（j+1）,z（j+1））-A（z（j+1）-1,z（j+1））*A（z（j+1）,z（j+1）-1）];
r（（z（j+1）-1）:
z（j+1））=roots（p）;
end
j=j+1;
end
ifn-z（m）==1
r（n）=A（n,n）;
else
p=[1,-（A（n-1,n-1）+A（n,n））,A（n-1,n-1）*A（n,n）-A（n-1,n）*A（n,n-1）];
r（n-1:
n）=roots（p）;
end
else
r
（1）=A（1,1）;
j=1;
whilejifz（j+1）-z（j）==1
r（z（j+1））=A（z（j+1）,z（j+1））;
end
if（z（j+1）-z（j）==2）
p=[1,-（A（z（j+1）-1,z（j+1）-1）+A（z（j+1）,z（j+1）））,A（z（j+1）-1,z（j+1）-1）*A（z（j+1）,z（j+1））-A（z（j+1）-1,z（j+1））*A（z（j+1）,z（j+1）-1）];
r（（z（j+1）-1）:
z（j+1））=roots（p）;
end
j=j+1;
end
ifn-z（m）==1
r（n）=A（n,n）;
else
p=[1,-（A（n-1,n-1）+A（n,n））,A（n-1,n-1）*A（n,n）-A（n-1,n）*A（n,n-1）];
r（n-1:
n）=roots（p）;
end
end
子程序2：
function[x]=Ifreducible（A）%判断H22是否已经schur形x=1表示还不是schur形，x=0表示已经是schur形
%y用来储存A次对角线元素
%z用来储存y元素为零的角标
n=size（A）;
y=zeros（1,n-1）;
x=1;
fori=1:
n-1
y（i）=A（i+1,i）;
end
m=0;
fori=1:
n-1
ifabs（y（i）-0）<1e-5
m=m+1;
end
end
ifm==0
x=1;
else
z=zeros（1,m）;
j=1;
i=1;
while（iifabs（y（i）-0）<1e-5
z（j）=i;
j=j+1;
end
i=i+1;
end
i=1;
while（iifz（i+1）-z（i）>2
x=1;
break;
end
i=i+1;
end
ifi>=m
x=0;
end
end
子程序3：
function[H22]=Francis（A）%QR迭代
H22=A;
[q,q]=size（H22）;
p=q-1;
s=H22（p,p）+H22（q,q）;
t=H22（p,p）*H22（q,q）-H22（p,q）*H22（q,p）;
x=H22（1,1）*H22（1,1）+H22（1,2）*H22（2,1）-s*H22（1,1）+t;
y=H22（2,1）*（H22（1,1）+H22（2,2）-s）;
z=H22（2,1）*H22（3,2）;
fork=0:
q-3;
[v,b]=house（[x,y,z]'）;
w=max（1,k）;
H22（k+1:
k+3,w:
q）=（eye（3）-b*v*v'）*H22（k+1:
k+3,w:
q）;
r=min（k+4,q）;
H22（1:
r,k+1:
k+3）=H22（1:
r,k+1:
k+3）*（eye（3）-b*v*v'）;
x=H22（k+2,k+1）;
y=H22（k+3,k+1）;
ifkz=H22（k+4,k+1）;
end
end
[v,b]=house（[x,y]'）;
H22（q-1:
q,q-2:
q）=（eye
（2）-b*v*v'）*H22（q-1:
q,q-2:
q）;
H22（1:
q,q-1:
q）=H22（1:
q,q-1:
q）*（eye
（2）-b*v*v'）;
子程序4：
function[v,b]=house（x）%house变换
n=length（x）;
m=max（abs（x））;
x=x/m;
q=x（2:
n）'*x（2:
n）;
v
（1）=1;
v（2:
n）=x（2:
n）;
ifq==0
b=0;
else
a=（x
（1）^2+q）^（1/2）;
ifx
（1）<=0
v
（1）=x
（1）-a;
else
v
（1）=-q/（x
（1）+a）;
end
b=2*v
（1）^2/（q+v
（1）^2）;
v=v/v
（1）;
end
v=v';
升级会员 