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

《计算方法》课程设计报告学生姓名:学号：学院:班级:题目:QR方法计算中小型矩阵的全部特征值指导教师：职称：教授讲师实验师2015年12月31日-I-目录目录……………………………………………………………………………I一、选题背景…………………………………………………………………11.1QR方法………………………………………………………………11.2矩阵的特征值………………………………………………………1二、算法设计…………………………………………………………………12.1QR方法的理论………………………………………………………12.2基本QR方法………………………………………………………22.3HouseholderQR分解……………………………………………22.4带原点位移的QR方法……………………………………………3三、程序设计及功能说明……………………………………………………43.1主要程序以及主要功能……………………………………………43.1.1矩阵约化为上海森伯格矩阵………………………………43.1.1QR方法求实矩阵全部特征值………………………………4四、结果分析…………..……………………………………………………6五、总结及心得体会…..……………………………………………………19参考文献…………….…..……………………………………………………20源程序………………..…..……………………………………………………21-1-一、选题背景1.1QR方法矩阵是高等数学中的常用工具，在很多方面都有重要运用，而矩阵特征值问题在许多领域的研究中有重要的地位，矩阵特征值的一些基本计算方法，研究不同种类矩阵的计算方法和最优计算方法.其中求解矩阵的普通方法包括传统的求法以及初等变换求矩阵的特征值方法；其他的一些优化方法包括幂法、反幂法、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.1QR方法的理论QR方法的理论：对任意一个非奇异矩阵（可逆矩阵）A，可以把它分解成一个正交阵Q和一个上三角阵R的乘积，称为对矩阵A的QR分解，即A=QR。如果规定R的对角元取正实数，这种分解是唯一的。若A是奇异的，则A有零特征值。任取一个不等于A的特征值的实数μ，则A-μI是非奇异的。只要求出A-μI的特征值和特征向量就容易求出矩阵A的特征值和特征向量，所以假设A是非奇异的，不失一般性。-2-2.2基本QR方法设A=A1，对A1作QR分解，得A1=Q1R1，，交换该乘积的次序，得11AQQR1Q1A2T，由于Q1正交矩阵，A1到A2的变换为正交相似变换，于是A1和A2就有相同的特征值。一般的令A1=A，对k=1，2，3，…..A(k)=QkRk（QR分解）A(k-1)=RkQk(迭代定义)这样，可得到一个迭代序列{Ak}，这就是QR方法的基本过程。2.3HouseholderQR分解从QR算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值，那麽涉及的计算量是很大的，因此应该先对A作预处理。应用中常先对做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素比三角矩阵多了一条次对角线，其形式为：上Hessenberg矩阵下Hessenberg矩阵一般矩阵相似约化到Hessenberg矩阵的方法。定理:任取非零向量X=(X1,X2,.....,Xn)T∈Rn，可以选择一个Householder矩阵P，使Px=－σｅ1式中e1=(1,0,......,0)T是Rn的单位向量，证：作u=x+σｅ1，用u做一个Householder矩阵P=I－β-1uuT,因为而证毕-3-定理中β=σ(σ+x1)，为避免σ+x1出现两个相似数相减引起有效数字损失，实用中常选这样可得Householder矩阵的计算公式为u1=x1+σβ=σu1P=I－β-1uuT2.4带原点位移的QR方法设A=A1，对A1-s1I进行QR分解A1-s1I=Q1R1；形成矩阵A2=R1Q1+s1I=TQ1(A-s1I)Q1+s1I=TQ1A1Q1；求得Ak后，将Ak-skI进行QR分解Ak-skI=QkRk，k=3，4，…形成矩阵A(k+1)=RkQk+skI=TQkAkQk.如果令Qk=Q1Q2…Qk,Rk=Rk…R2R1,则有A(k+1)=TRkAQk，并且矩阵(A-s1I)(A-s2I)…(A-snI)≡ψ(A)有QR分解式ψ(A)=QkRk.也可以首先用正交变换(左变换)将Ak-skI化为上三角矩阵，即Pn…P2P1(Ak-skI)=Rk当A为Householder矩阵或对称三对角矩阵，Pi可为平面旋转矩阵，则Ak=P(n-1)…P2P1(Ak-skI)TP1TP2…TP(n-1)+skI.-4-三、程序及功能说明3.1主要程序以及主要功能3.1.1矩阵约化为上海森伯格矩阵function[k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)n=size(A);Ak=A;fork=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程序功能：通过初等反射矩阵正交相似约化实矩阵A为上海森伯格矩阵Ak。来实现矩阵约化为上海森伯格矩阵，以方便使用QR方法求矩阵A的全部特征值。3.1.2QR方法求实矩阵全部特征值functiontzg=qr4(A,t,max1)[n,n]=size(A);k=0;Ak=A;tzg=zeros(n);state=1;fori=1:n;while((k=max1)&(state==1)&(n1))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-5-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个有效数字的近似值。-6-四、结果分析实验结果分析：1.用QR方法求下列矩阵的全部近似特征值，其中,620541321A精度为510；解：（1）先将矩阵转化为上海森伯格矩阵，MATLAB程序如下：A=[123;145;026];[k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)结果如下：k=1Sk=1uk=20ck=2Pk=-1001Uk=1000-10001-7-Ak=123-1-4-5026然后用最末元位移QR方法求实矩阵Ak全部特征值.MATLAB程序如下：A=[123;-1-4-5;026];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.27330.1020-5.3467-2.5877-3.9784-5.583700.00000Qk=-0.94210.33520.0000-0.3352-0.9421-0.000000.0000-1.0000-8-Rk=7.71991.23746.909103.78243.4684000.0000At=-2.93751.4219-6.9091-1.26791.1870-3.468400.00004.7505请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数，下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i=1tzgk=4.7505k=4sk=4.7505B=-2.93751.4219-1.26791.1870请注意：下面的i表示求第i个特征值，k是迭代次数，sk是原点位移量，Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵：-9-i=2k=7sk=0.6930Bk=-3.13662.6863-0.00350Qk=-1.00000.0011-0.0011-1.0000Rk=3.1366-2.686300.0030At=-2.44062.6897-0.00000.6900请注意：i表示求第i个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数，下面的矩阵B是m阶矩阵At的（m-1）阶主子矩阵，即At降一阶.i=2tzgk=0.6900-10-k=7sk=0.6930B=-2.4406tzgk=-2.4406请注意：n阶实对称矩阵A的全部真特征值lamoda和至少含t个有效数字的近似特征值tzg如下：lamoda=-2.44060.69004.7505tzg=-2.44060.69004.75052.用QR方法求下列矩阵的全部近似特征值，其中,122212221A精度为510；解：MATLAB程序如下：A=[122;212;221];tzg=qr4(A,5,100)结果如下：-11-请注意：下面的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.94180.45530