实验七、QR 算法一、实验目的1、熟悉matlab 编程并学习QR 算法原理及计算机实现；2、学习用matlab 内置函数eig 和QR 算法求矩阵的特征值，并比对二者差异。

二、实验题目1、课本第277页第1题已知矩阵11261112376111671123456110787445677565,,.0367886109002897591000010A B H ⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪ ⎪⎪⎪=== ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎪⎝⎭（1）用MATLAB 函数“eig ”求矩阵全部特征值；（2）用基本QR 算法求全部特征值（可用MA TLAB 函数“qr ”实现矩阵的QR 分解）。

2、用QR 算法求矩阵特征值：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=111132126)(i ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=010*******876307654465432)(ii根据QR 算法原理编制求（i ）及（ii ）中矩阵全部特征值的程序并输出计算结果（要求误差<10 -5）.三、实验原理与理论基础QR 方法是一种变换方法，是计算一般矩阵（中小型矩阵）全部特征值问题的最有效方法之一。

目前QR 方法主要用来计算上海森伯格矩阵和对称三对角矩阵的全部特征值问题，且QR 方法具有收敛快、算法稳定等特点。

对于一般矩阵n nA ⨯∈（或对称矩阵），首先用豪斯霍尔德方法将A 化为上海森伯格矩阵B （或对称三对角矩阵），然后再用QR 方法计算B 的全部特征值。

1、矩阵的QR 分解设n nA ⨯∈非奇异，则存在正交矩阵P ，使PA=R ，其中R 为上三角矩阵。

用Householder变换构造正交矩阵P ，记(0)AA =，它的第一列记为(0)1a ，不妨设(0)10a ≠，可按公式（3.2）（Th14,约化定理 设12(,,,)0,T n x x x x =≠则存在初等反射矩阵H 使1Hx e σ=-，其中）112121122,sgn(),,().T H I uu x x u x e u x βσσβσσ-⎧=-⎪=⎪⎨=+⎪⎪==+⎩ 找到矩阵111111,,n nTH H I u u β⨯-∈=-使(0)11111,(1,0,,0).T nH a e e σ=-=∈于是(1)1(1)(0)(0)(0)(1)111121(1)(,,,),0n b AH A H a H a H a A σ⎛⎫-===⎪ ⎪⎝⎭其中(1)(1)(1)(1)(1)(1)121(,,,).n n n Aa a a -⨯--=∈一般地，设(1)(1)(1)()0j j j D B AA ---⎛⎫= ⎪ ⎪⎝⎭，其中(1)j D -为(j-1)阶方阵，其对角线以下元素均为0，(1)j A-为(n-j+1)阶方阵，设其第一列为(1)1j a -，可选择(n-j+1)的Householder 矩阵变换()()n j n j j H -⨯-∈，使(1)1111,(1,0,,0),j n j j j H a e e σ--+=-=∈根据j H 构造n*n 阶的变换矩阵j H 为10,0j j j I H H -⎛⎫=⎪⎝⎭于是有()()()(1).0j j j j j j D B A H A A -⎛⎫== ⎪ ⎪⎝⎭它和(1)j A -有类似的形式，只是()j D 为j 阶方阵，其对角线以下元素是0，这样经过n-1步运算得到(1)11,n n H H A A R --==其中(1)n R A -=为上三角矩阵，11n P H H -=为正交矩阵，从而有PA=R 。

2、QR 算法设n nA ⨯∈，且对A 进行QR 分解，即A QR =，其中R 为上三角矩阵，Q 为正交矩阵，于是可得到一个新矩阵TB RQ Q AQ ==。

显然，B 是由A 经过正交相似变换得到，因此B 与A 特征值相同，再对B 进行QR 分解，又可得一新的矩阵，重复这一过程可得到矩阵序列：设1A A =将1A 进行QR 分解111A Q R =作矩阵211111TA R Q Q A Q ==求得k A 后将k A 进行QR 分解k k k A Q R =形成矩阵1Tk k k k k k A R Q Q A Q +==QR 算法，就是利用矩阵的QR 分解，按上述递推法则构造矩阵序列{}k A 的过程。

只要A 为非奇异矩阵，则由QR 算法就完全确定{}k A 。

四、实验内容1、用matlab内置函数eig求矩阵的全部特征值；2、编写求特征值的QR算法程序，并用之求矩阵特征值；3、比较两种方法的结果差异。

（1）QR算法的m文件function qrsf(A,r)[Q,R]=qr(A);t=A(1,1) %tempA=R*Q;for k=1:50[Q,R]=qr(A);t=A(1,1);A=R*Q;if( abs( A(1,1)-t )<r )break;endendn=size(A,1);for i=1:nformat long gdisp( ['特征值λ',num2str(i),'=',num2str( A(i,i) )] ); end%disp('');for i=1:ndisp('特征值');format long g,A(i,i)endformat long g,A,Q,R（2）改进后的QR算法的m文件function qrsf(A,r)[Q,R]=qr(A);%t=A(1,1) %tempt(1)=max(abs(diag(R)));A=R*Q;for k=2:50[Q,R]=qr(A);z=diag(A);t(k)=max(abs(diag(R)));A=R*Q;if( abs( t(k) - t(k-1) ) < r )break;endendn=size(A,1);for i=1:nformat long gdisp( ['特征值λ',num2str(i),'=',num2str( A(i,i) )] ); end%disp('');for i=1:ndisp('特征值');format long g,A(i,i)endformat long g,A,Q,R五、实验结果>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];>> B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0];>> H6=hilb(6);1、eig求矩阵特征值>> eig(A),eig(B),eig(H6)ans = 0.01015004839789240.8431071498550323.8580574559449530.2886853458021ans = 13.17235139810326.551878351915661.59565457314994-0.390788045416488-0.929096277752298ans = 1.08279948406811e-0071.25707571226224e-0050.0006157483541826520.01632152131987580.242360870575211.618899858924342、QR算法求矩阵特征值>> eig(A),qrsf(A,10^-8)ans = 0.01015004839789240.8431071498550323.8580574559449530.2886853458021t = 10特征值λ1=30.2887特征值λ2=3.8581特征值λ3=0.84311特征值λ4=0.01015特征值ans = 30.2886853457915特征值ans = 3.85805737835431 特征值ans = 0.843107227456257 特征值ans = 0.0101500483978911 >> eig(B),qrsf(B,10^-8)ans = 13.17235139810326.551878351915661.59565457314994-0.390788045416488-0.929096277752298t = 2特征值λ1=13.1724特征值λ2=6.5519特征值λ3=1.5957特征值λ4=-0.9291特征值λ5=-0.39079特征值ans = 13.1723513891479 特征值ans = 6.55187836087093 特征值ans = 1.59565457937031 特征值ans = -0.929096283974607 特征值ans = -0.390788045414554 >> eig(H6),qrsf(H6,10^-8)ans = 1.08279948406811e-0071.25707571226224e-0050.0006157483541826520.01632152131987580.242360870575211.61889985892434t = 1特征值λ1=1.6189特征值λ2=0.24236特征值λ3=0.016322特征值λ4=0.00061575特征值λ5=1.2571e-005特征值λ6=1.0828e-007特征值ans = 1.6188998588068 特征值ans = 0.24236087069274 特征值ans = 0.01632152131988 特征值ans = 0.000615748354182639 特征值ans = 1.25707571226506e-005 特征值ans = 1.08279948456401e-0073、改进后的QR算法求特征值>> eig(A),qrsf(A,10^-8)ans = 0.01015004839789240.8431071498550323.8580574559449530.2886853458021特征值λ1=30.2887特征值λ2=3.8581特征值λ3=0.84311特征值λ4=0.01015特征值ans = 30.2886853458019 特征值ans = 3.85805745223919 特征值ans = 0.843107153560957 特征值ans = 0.0101500483978911 >> eig(B),qrsf(B,10^-8)ans = 13.17235139810326.551878351915661.59565457314994-0.390788045416488-0.929096277752298特征值λ1=13.1724特征值λ2=6.5519特征值λ3=1.5957特征值λ4=-0.9291特征值λ5=-0.39079特征值ans = 13.1723513936489 特征值ans = 6.55187835636998 特征值ans = 1.59565456952802 特征值ans = -0.929096274131193 特征值ans = -0.390788045415675 >> eig(H6),qrsf(H6,10^-8)ans = 1.08279948406811e-0071.25707571226224e-0050.0006157483541826520.01632152131987580.242360870575211.61889985892434特征值λ1=1.6189特征值λ2=0.24236特征值λ3=0.016322特征值λ4=0.00061575特征值λ5=1.2571e-005特征值λ6=1.0828e-007特征值ans = 1.61889985892171 特征值ans = 0.242360870577844 特征值ans = 0.0163215213198758 特征值ans = 0.000615748354182638 特征值ans = 1.25707571226506e-005特征值ans = 1.08279948456401e-007六、实验结果分析与小结从实验结果可以看出，用MA TLAB内置函数eig求矩阵特征值与用QR算法求矩阵特征值的结果基本一致，数据只有微小差别。