2022年隐式QR法求实矩阵的全部特征值matlab实现 .pdf

《2022年隐式QR法求实矩阵的全部特征值matlab实现 .pdf》由会员分享，可在线阅读，更多相关《2022年隐式QR法求实矩阵的全部特征值matlab实现 .pdf（6页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、隐式Q R法求实矩阵的全部特征值m a tla b实现要求：用 matlab 编写通用子程序，利用隐式QR法求实矩阵的全部特征值和特征向量。思想：隐式 QR法实质上就是将一个矩阵 Schur 化，之后求解特征值就比较方便。而隐式 QR法还需要用到 household 变换，以及上 hessenberg 变换。最后使用 QR迭代，达到Schur 化的结果。步骤：1. 可约性判定， 也就是判断次对角线元素是否非零，如果次对角线元素非零， 则不可约。2.Schur 化，也就是通过QR迭代，将矩阵 H变化成为某些次对角线元素变成0，同时还要满足，这些元素之间间隔最大为1，那么，所得到的最重的矩阵H就是

2、一个 Schur 形矩阵。3. 假如两个等于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 ）求01341xx的所有根。附录 2 隐式 QR

3、迭代：主程序：function r,V=SchurQR(A)% 向量r 用来储存特征值%Hessenberg分解:m,m=size(A);for k=1:m-2 v,b=house(A(k+1:m,k); H1= eye(m-k)-b*v*v; H2=eye(m);for i=k+1:m名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 6 页 - - - - - - - - - for j=k+1:m H2(i,j)=H1(i-k,j-k);endendif k=1; H=

4、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;endu=10e-5;for i=2:m;if abs( A(i,i-1)=(abs(A(i,i)+ abs(A(i-1,i-1)*u; A(i,i-1)=0;endend%QR 迭代：H22=A;x=Ifreducible(H22);while x=1 H22=Francis(H22); x=Ifreducible(H22);endr,V=EigValue(H22);子程序 1：function r,V=EigValue(A)% 计算A

5、的特征值，特征向量n,n=size(A);r=zeros(1,n);y=zeros(1,n-1);%y用来储存次对角线元素for i=1:n-1 y(i)=A(i+1,i);endm=0;for i=1:n-1if abs(y(i)-0)1e-5 m=m+1;endendif m=0 x=1;名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 6 页 - - - - - - - - - elsez=zeros(1,m);%z用来储存值为 0的y向量的角标。j=1;i=1;wh

6、ile(in)if abs(y(i)-0)1e-5 z(j)=i; j=j+1;end i=i+1;endendif z(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;while jmif z(j+1)-z(j)=1 r(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

7、,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;endif n-z(m)=1 r(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);endelse r(1)=A(1,1); j=1;while jmif z(j+1)-z(j)=1 r(z(j+1)=A(z(j+1),z(j+1);end名师资料总结 -

8、- -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 6 页 - - - - - - - - - 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;endif n-z(m)=1 r(n)=A(n,n);

9、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：functionx=Ifreducible(A)%判 断 H22 是 否 已 经 schur形x=1 表 示还 不是schur 形，x=0 表示已经是 schur 形%y用来储存 A次对角线元素%z用来储存 y元素为零的角标n=size(A);y=zeros(1,n-1);x=1;for i=1:n-1 y(i)=A(i+1,i);endm=0;for i=1:n-1if abs(y(i)-0)1e-5 m

10、=m+1;endendif m=0 x=1;elsez=zeros(1,m);j=1;i=1;while(in)if abs(y(i)-0)1e-5名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 6 页 - - - - - - - - - z(j)=i; j=j+1;end i=i+1;endi=1;while(i2 x=1;break ;end i=i+1; endif i=m x=0;endend子程序 3：function H22=Francis(A)%QR 迭代H

11、22=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);for k=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); H

12、22(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);if kq-3 z=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：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - -

13、- - 第 5 页，共 6 页 - - - - - - - - - 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);if q=0 b=0;else a=(x(1)2+q)(1/2);if x(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);endv=v;名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 6 页 - - - - - - - - -