隐式QR法求实矩阵的全部特征值matlab实现Word格式文档下载.docx

隐式QR法求实矩阵的全部特征值matlab实现Word格式文档下载.docx

《隐式QR法求实矩阵的全部特征值matlab实现Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《隐式QR法求实矩阵的全部特征值matlab实现Word格式文档下载.docx（12页珍藏版）》请在冰豆网上搜索。

求x=0.9,1.0,1.1时的特征值和特征向量。

X=0.9：

r是特征值，V是特征向量矩阵。

X=1：

X=1.1：

（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:

H2（i,j）=H1（i-k,j-k）;

end

ifk==1;

H=H2;

else

H=H*H2;

A（k+1:

m,k:

m）=H1*A（k+1:

m）;

A（1:

m,k+1:

m）=A（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;

%QR迭代：

H22=A;

x=Ifreducible（H22）;

whilex==1

H22=Francis（H22）;

x=Ifreducible（H22）;

[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）;

m=0;

ifabs（y（i）-0）<

1e-5

m=m+1;

ifm==0

x=1;

else

z=zeros（1,m）;

%z用来储存值为0的y向量的角标。

j=1;

i=1;

while（i<

n）

z（j）=i;

j=j+1;

i=i+1;

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;

whilej<

ifz（j+1）-z（j）==1

r（z（j+1））=A（z（j+1）,z（j+1））;

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）;

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）;

r

（1）=A（1,1）;

子程序2：

function[x]=Ifreducible（A）%判断H22是否已经schur形x=1表示还不是schur形，x=0表示已经是schur形

%y用来储存A次对角线元素

%z用来储存y元素为零的角标

n=size（A）;

x=1;

m）

ifz（i+1）-z（i）>

2

break;

ifi>

=m

x=0;

子程序3：

function[H22]=Francis（A）%QR迭代

[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:

q）;

r=min（k+4,q）;

H22（1:

r,k+1:

k+3）=H22（1:

k+3）*（eye（3）-b*v*v'

x=H22（k+2,k+1）;

y=H22（k+3,k+1）;

ifk<

q-3

z=H22（k+4,k+1）;

[v,b]=house（[x,y]'

H22（q-1:

q,q-2:

q）=（eye

（2）-b*v*v'

）*H22（q-1:

H22（1:

q,q-1:

q）=H22（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:

ifq==0

b=0;

a=（x

（1）^2+q）^（1/2）;

ifx

（1）<

=0

v

（1）=x

（1）-a;

v

（1）=-q/（x

（1）+a）;

b=2*v

（1）^2/（q+v

（1）^2）;

v=v/v

（1）;

v=v'