1、隱式QR法求實(shí)矩陣的全部特征值matlab實(shí)現(xiàn)要求：用matlab編寫(xiě)通用子程序，利用隱式QR法求實(shí)矩陣的全部特征值和特征向量。思想：隱式QR法實(shí)質(zhì)上就是將一個(gè)矩陣Schur化，之后求解特征值就比較方便。而隱式QR法還需要用到household變換，以及上hessenberg變換。最后使用QR迭代，達(dá)到Schur化的結(jié)果。步驟：1.可約性判定，也就是判斷次對(duì)角線元素是否非零，如果次對(duì)角線元素非零，則不可約。2.Schur化，也就是通過(guò)QR迭代，將矩陣H變化成為某些次對(duì)角線元素變成0,同時(shí)還要滿足，這些元素之間間隔最大為1，那么，所得到的最重的矩陣H就是一個(gè)Schur形矩陣。3.假如兩個(gè)等于0的

2、次對(duì)角線元素間隔為0，那么該元素的上面一個(gè)元素，也就是H的對(duì)角線上的元素，即為其中一個(gè)特征值；假如兩個(gè)等于0的次對(duì)角線元素間隔為1，那么在這兩個(gè)元素之間就形成了一個(gè)2*2的矩陣，可以求解一個(gè)一元二次方程來(lái)得到兩個(gè)共軛的特征值。實(shí)驗(yàn)代碼：詳見(jiàn)附錄2實(shí)驗(yàn)結(jié)果：（代碼相見(jiàn)附錄2）（i）設(shè)矩陣A如下：求x=,時(shí)的特征值和特征向量。X=：r是特征值，V是特征向量矩陣。X=1：r是特征值，V是特征向量矩陣。X=：r是特征值，V是特征向量矩陣。（ii）求x41x310的所有根。附錄2隱式QR迭代：主程序：functionr,V=SchurQR(A)%向量r用來(lái)儲(chǔ)存特征值%Hessenberg分解:m,m=s

3、ize(A);fork=1:m-2v,b=house(A(k+1:m,k);H1=eye(m-k)-b*v*v'H2=eye(m);fori=k+1:mforj=k+1:mH2(i,j)=H1(i-k,j-k);endendifk=1;H=H2;elseH=H*H2;endA(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;endu=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;endend%Q嵌代：H22=A;x=Ifr

4、educible(H22);whilex=1H22=Francis(H22);x=Ifreducible(H22);endr,V=EigValue(H22);子程序1：functionr,V=EigValue(A)%算A的特征值，特征向量n,n=size(A);r=zeros(1,n);y=zeros(1,n-1);%擁來(lái)儲(chǔ)存次對(duì)角線元素fori=1:n-1y(i)=A(i+1,i);endm=0;fori=1:n-1ifabs(y(i)-0)<1e-5m=m+1;endendifm=0x=1;elsez=zeros(1,m);%和來(lái)儲(chǔ)存值為0的y向量的角標(biāo)。j=1;i=1;while(

5、i<n)ifabs(y(i)-0)<1e-5z(j)=i;j=j+1;endi=i+1;endendifz(1)=2%次對(duì)角線第一個(gè)等于0的元素的位置不同，需要2分類(lèi)討論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<mifz(j+1)-z(j)=1r(z(j+1)=A(z(j+1),z(j+1);endif(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

6、)*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);endj=j+1;endifn-z(m)=1r(n)=A(n,n);elsep=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);endelser(1)=A(1,1);j=1;whilej<mifz(j+1)-z(j)=1r(z(j+1)=A(z(j+1),z(j+1);endif(z(j+1)-z(j)=2)p=1,-(A

7、(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);endj=j+1;endifn-z(m)=1r(n)=A(n,n);elsep=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);endend子程序2：%|J斷H22M否已經(jīng)schur形x=1表示還不是schur形，x=0

8、表示已經(jīng)是functionx=Ifreducible(A)schur形%擁來(lái)儲(chǔ)存A次對(duì)角線元素%和來(lái)儲(chǔ)存y元素為零的角標(biāo)n=size(A);y=zeros(1,n-1);x=1;fori=1:n-1y(i)=A(i+1,i);endm=0;fori=1:n-1ifabs(y(i)-0)<1e-5m=m+1;endendifm=0x=1;elsez=zeros(1,m);j=1;i=1;while(i<n)ifabs(y(i)-0)<1e-5z(j)=i;j=j+1;endi=i+1;endi=1;while(i<m)ifz(i+1)-z(i)>2x=1;break

9、;endi=i+1;endifi>=mx=0;endend子程序3：functionH22=Francis(A)%Q嵌代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

10、+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);ifk<q-3z=H22(k+4,k+1);endendv,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：functionv,b=house(x)%hous吸換n=length(x);m=max(abs(x);x=