隐式QR法求实矩阵的全部特征值matlab实现.docx
- 文档编号:23306115
- 上传时间:2023-05-16
- 格式:DOCX
- 页数:12
- 大小:59.59KB
隐式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(i ifabs(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; whilej ifz(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; whilej ifz(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(i ifabs(y(i)-0)<1e-5 z(j)=i; j=j+1; end i=i+1; end i=1; while(i ifz(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); ifk z=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';
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 隐式 QR 求实 矩阵 全部 特征值 matlab 实现