qr方法计算中小型矩阵的全部特征值

计算方法课程设计报告学生姓名:学 号：学 院:班 级:题 目:QR方法计算中小型矩阵的全部特征值指导教师： 职称： 教 授 讲 师 实验师 2015年12月31日目 录目录I一、选题背景11.1 QR方法11.2 矩阵的特征值 1二、算法设计12.1 QR方法的理论12.2 基本QR方法 22.3 Householder QR分解 22.4 带原点位移的QR方法 3三、程序设计及功能说明43.1主要程序以及主要功能43.1.1矩阵约化为上海森伯格矩阵43.1.1 QR方法求实矩阵全部特征值4四、结果分析. 6五、总结及心得体会.19参考文献.20源程序.21一、选题背景1.1 QR方法矩阵是高等数学中的常用工具，在很多方面都有重要运用，而矩阵特征值问题在许多领域的研究中有重要的地位，矩阵特征值的一些基本计算方法，研究不同种类矩阵的计算方法和最优计算方法.其中求解矩阵的普通方法包括传统的求法以及初等变换求矩阵的特征值方法；其他的一些优化方法包括幂法、反幂法、Jacobi方法、QR方法.在实际的求解矩阵特征值的问题，根据矩阵的不同特点，选择最快速的方法求解，从而达到最优化解决实际问题。QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为Hessenberg矩阵，然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法，与此正规正交矩阵的通用符号Q有关。目前QR方法主要用来计算：(1) Hessenberg矩阵的全部特征值问题；(2)计算对称三对角矩阵的全部特征值问题。1.2 矩阵的特征值关于计算矩阵A的特征值问题，当n=2,3时，我们还可按行列式展开的办法求，我们还可按行列式展开的办法求det()=0的根，但当n较大时，如果按展开行列式的办法，首先求出det()的系数，再求det()的根，工作量就非常大，用这种办法求矩阵用这种办法求矩阵的特征值是不切实际的的特征值是不切实际的，由此需要研究求A的特征值 及特征向量的数值解法。二、算法设计2.1 QR方法的理论QR方法的理论：对任意一个非奇异矩阵（可逆矩阵）A，可以把它分解成一个正交阵Q和一个上三角阵R的乘积，称为对矩阵A的QR分解，即A=QR。如果规定R的对角元取正实数，这种分解是唯一的。若A是奇异的，则A有零特征值。任取一个不等于A的特征值的实数，则A-I是非奇异的。只要求出A-I的特征值和特征向量就容易求出矩阵A的特征值和特征向量，所以假设A是非奇异的，不失一般性。2.2 基本QR方法 设A=A1，对A1作QR分解，得A1= Q1R1，交换该乘积的次序，得，由于Q1正交矩阵，A1到A2的变换为正交相似变换，于是A1和A2就有相同的特征值。一般的令A1=A，对k=1，2，3，. A(k)=QkRk （QR分解） A(k-1)=RkQk (迭代定义) 这样，可得到一个迭代序列Ak，这就是QR方法的基本过程。2.3 Householder QR分解从QR 算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值，那麽涉及的计算量是很大的，因此应该先对A作预处理。应用中常先对 做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素比三角矩阵多了一条次对角线，其形式为： 上Hessenberg矩阵 下Hessenberg矩阵一般矩阵相似约化到Hessenberg矩阵的方法。定理: 任取非零向量X=(X1,X2,.,Xn)TRn，可以选择一个Householder矩阵P，使Px=1式中e1=(1,0,.,0)T是Rn的单位向量， 证： 作u=x+1，用u做一个Householder矩阵P=I-1uuT, 因为 而证毕定理中=(+x1)，为避免+x1出现两个相似数相减引起有效数字损失，实用中常选 这样可得Householder矩阵的计算公式为u1=x1+=u1 P=I-1uuT2.4 带原点位移的QR方法设A=A1，对A1-s1I进行QR分解A1-s1I=Q1R1；形成矩阵 A2=R1Q1+s1I=1(A-s1I) Q1+s1I=1A1Q1；求得Ak后，将Ak-skI进行QR分解Ak-skI=QkRk，k=3，4， 形成矩阵A(k+1)=RkQk+skI=kAkQk.如果令k=Q1Q2Qk, k=RkR2R1,则有A(k+1)= kAk，并且矩阵(A-s1I)(A-s2I) (A-snI) (A)有QR分解式(A)= kk.也可以首先用正交变换(左变换)将Ak-skI化为上三角矩阵，即PnP2P1(Ak-skI)=Rk当A为Householder矩阵或对称三对角矩阵，Pi可为平面旋转矩阵，则Ak= P(n-1)P2P1(Ak-skI) 12(n-1)+skI.三、程序及功能说明3.1主要程序以及主要功能3.1.1 矩阵约化为上海森伯格矩阵function k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k)*sign(Ak(k+1,k), uk= Ak(k+1:n,k)+ Sk*eye(n-k,1),ck=(norm(uk,2)2)/2,Pk= eye(n-k,n-k)-uk*uk/ck,Uk=eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk,A1=Uk*Ak;Ak=A1,end 程序功能：通过初等反射矩阵正交相似约化实矩阵为上海森伯格矩阵Ak。来实现矩阵约化为上海森伯格矩阵，以方便使用QR方法求矩阵A的全部特征值。3.1.2 QR方法求实矩阵全部特征值function tzg=qr4(A,t,max1)n,n=size(A); k=0;Ak=A;tzg=zeros(n); state=1;for i=1:n;while(k1)b1=abs(Ak(n,n-1); b2=abs(Ak(n,n); b3=abs(Ak(n-1,n-1);b4=min(b2, b3); jd=10(-t); jd1=jd*b4;if(b1=jd1) sk=Ak(n,n); Bk=Ak-sk*eye(n); Qk,Rk=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp(请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量，)disp( Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：)i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;else disp(请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数，) disp( 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.) i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state=1; i=i+1;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1);tzgk=Akdisp(请注意：n阶实对称矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：)lamoda=sort(eig(A)程序功能：QR方法来计算Hessenberg矩阵或对称的矩阵的全部特征值至少含t个有效数字的近似值。四、结果分析实验结果分析：1.用QR方法求下列矩阵的全部近似特征值，其中精度为；解：（1）先将矩阵转化为上海森伯格矩阵，MATLAB程序如下： A=1 2 3;1 4 5;0 2 6; k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)结果如下：k = 1Sk = 1uk = 2 0ck = 2Pk = -1 0 0 1Uk = 1 0 0 0 -1 0 0 0 1Ak = 1 2 3 -1 -4 -5 0 2 6然后用最末元位移QR方法求实矩阵k全部特征值.MATLAB程序如下： A=1 2 3;-1 -4 -5;0 2 6; tzg=qr4(A,5,100)结果如下：请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 1k = 4sk = 4.7505Bk = -7.2733 0.1020 -5.3467 -2.5877 -3.9784 -5.5837 0 0.0000 0Qk = -0.9421 0.3352 0.0000 -0.3352 -0.9421 -0.0000 0 0.0000 -1.0000Rk = 7.7199 1.2374 6.9091 0 3.7824 3.4684 0 0 0.0000At = -2.9375 1.4219 -6.9091 -1.2679 1.1870 -3.4684 0 0.0000 4.7505请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 1tzgk = 4.7505k = 4sk = 4.7505B = -2.9375 1.4219 -1.2679 1.1870请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 2k = 7sk = 0.6930Bk = -3.1366 2.6863 -0.0035 0Qk = -1.0000 0.0011 -0.0011 -1.0000Rk = 3.1366 -2.6863 0 0.0030At = -2.4406 2.6897 -0.0000 0.6900请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 2tzgk = 0.6900k = 7sk = 0.6930B = -2.4406tzgk = -2.4406请注意：n阶实对称矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：lamoda = -2.4406 0.6900 4.7505tzg = -2.4406 0.6900 4.75052.用QR方法求下列矩阵的全部近似特征值，其中精度为；解：MATLAB程序如下： A=1 2 2;2 1 2;2 2 1; tzg=qr4(A,5,100)结果如下：请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 1k = 3sk = -0.9883Bk = 5.9418 0.4553 0.2639 0.4553 0.0231 0.0202 0.2639 0.0202 0Qk = 0.9961 0.0764 0.0441 0.0763 -0.9971 0.0034 0.0442 0.0000 -0.9990Rk = 5.9650 0.4562 0.2644 0 0.0117 0.0000 0 0 0.0117At = 5.0000 0.0009 0.0005 0.0009 -1.0000 0.0000 0.0005 0.0000 -1.0000请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 1tzgk = -1.0000k = 3sk = -0.9883B = 5.0000 0.0009 0.0009 -1.0000请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 2k = 4sk = -1.0000Bk = 6.0000 0.0009 0.0009 0Qk = 1.0000 0.0001 0.0001 -1.0000Rk = 6.0000 0.0009 0 0.0000At = 5.0000 0.0000 0.0000 -1.0000请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 2tzgk = -1.0000k = 4sk = -1.0000B = 5.0000tzgk = 5.0000请注意：n阶实对称矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：lamoda = -1.0000 -1.0000 5.0000tzg = -1.0000 -1.00005.00003.用QR方法求下列矩阵的全部近似特征值，其中精度为；解：（1）先将矩阵转化为上海森伯格矩阵，MATLAB程序如下： A=1 2 3 4;2 3 4 5;3 4 5 6;4 5 6 7; k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)结果如下：k = 2Sk = -0.4549uk = -0.5400 -0.4468ck = 0.2456Pk = -0.1871 -0.9823 -0.9823 0.1871Uk = 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.1871 -0.9823 0 0 -0.9823 0.1871Ak = 1.0000 2.0000 3.0000 4.0000 -5.3852 -7.0564 -8.7277 -10.3989 0.0000 0.4549 0.9097 1.36460.0000 -0.0000 -0.0000 -0.0000（2）然后用最末元位移QR方法求实对称矩阵全部特征值，因为矩阵非奇异，故MATLAB程序如下： A=1 2 3;-5.3852 -7.0564 -8.7277;0.0000 0.4549 0.9097; tzg=qr4(A,5,100)结果如下：请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 1k = 6sk = -1.4987e-004Bk = -4.6826 11.4681 -3.9762 -0.0002 -0.4637 -0.7271 0 0.0000 0Qk = -1.0000 0.0000 0.0000 -0.0000 -1.0000 -0.0000 0 0.0000 -1.0000Rk = 4.6826 -11.4681 3.9763 0 0.4641 0.7269 0 0 0.0000At = -4.6823 11.4683 -3.9763 -0.0000 -0.4642 -0.7269 0 0.0000 -0.0001请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 1tzgk = -1.4992e-004k = 6sk = -1.4987e-004B = -4.6823 11.4683 -0.0000 -0.4642请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i = 2k = 7sk = -0.4642Bk = -4.2181 11.4683 -0.0000 0Qk = -1.0000 0.0000 -0.0000 -1.0000Rk = 4.2181 -11.4683 0 0.0000At = -4.6823 11.4683 -0.0000 -0.4643请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i = 2tzgk = -0.4643k = 7sk = -0.4642B = -4.6823tzgk = -4.6823请注意：n阶矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：lamoda = -4.6823 -0.4643 -0.0001tzg = -4.6823 -0.4643 -0.0001五、总结及心得体会在本次课程设计中，虽然我们在完成过程中，遇到这样那样问题，通过老师和同学帮助，我们算是基本熟悉QR方法的基本原理；得到QR方法求解矩阵特征值的通用程序完成本次课程设计，使我们了解各种求解矩阵特征值计算方法的之间优缺点。参考文献1 李庆扬，王能超，易大义数值分析（第5版）北京：清华大学出版社，2008.2 李庆扬，关治，白峰杉 数值计算原理北京：清华大学出版社，20003 李海，邓樱MATLAB程序设计教程北京：高等教育出版社，20024 王萼芳,石生明高等代数M北京:高等教育出版社2003源程序1 矩阵约化为上海森伯格矩阵function k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k)*sign(Ak(k+1,k), uk= Ak(k+1:n,k)+ Sk*eye(n-k,1),ck=(norm(uk,2)2)/2,Pk= eye(n-k,n-k)-uk*uk/ck,Uk=eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk,A1=Uk*Ak;Ak=A1,end 2 QR方法求实矩阵全部特征值的function tzg=qr4(A,t,max1)n,n=size(A); k=0;Ak=A;tzg=zeros(n); state=1;for i=1:n;while(k1)b1=abs(Ak(n,n-1); b2=abs(Ak(n,n); b3=abs(Ak(n-1,n-1);b4=min(b2, b3); jd=10(-t); jd1=jd*b4;if(b1=jd1) sk=Ak(n,n); Bk=Ak-sk*eye(n); Qk,Rk=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp(请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量，)disp( Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：)i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;else disp(请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数，) disp( 下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.) i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state=1; i=i+1;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1);tzgk=Akdisp(请注意：n阶实对称矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：)lamoda=sort(eig(A)