(参考)数值分析实验报告7.docVIP

实验七、QR算法 一、实验目的 二、实验题目 （1）用MATLAB函数“eig”求矩阵全部特征值； （2）用基本QR算法求全部特征值（可用MATLAB函数“qr”实现矩阵的QR分解）。 2、用QR算法求矩阵特征值： 根据QR算法原理编制求（i）及（ii）中矩阵全部特征值的程序并输出计算结果（要求误差10 -5）. 三、实验原理与理论基础 （或对称矩阵），首先用豪斯霍尔德方法将A化为上海森伯格矩阵B（或对称三对角矩阵），然后再用QR方法计算B的全部特征值。 1、矩阵的QR分解 设非奇异，则存在正交矩阵P，使PA=R，其中R为上三角矩阵。用Householder变换构造正交矩阵P，记，它的第一列记为，不妨设，可按公式（3.2）（Th14,约化定理 设则存在初等反射矩阵H使，其中） 找到矩阵使 于是 其中 一般地，设，其中为(j-1)阶方阵，其对角线以下元素均为0，为(n-j+1)阶方阵，设其第一列为，可选择(n-j+1)的Householder矩阵变换，使根据构造n*n阶的变换矩阵为于是有它和有类似的形式，只是为j阶方阵，其对角线以下元素是0，这样经过n-1步运算得到其中为上三角矩阵，为正交矩阵，从而有PA=R。 2、QR算法 设，且对A进行QR分解，即，其中R为上三角矩阵，Q为正交矩阵，于是可得到一个新矩阵。显然，B是由A经过正交相似变换得到，因此B与A特征值相同，再对B进行QR分解，又可得一新的矩阵，重复这一过程可得到矩阵序列： 设 将进行QR分解 作矩阵 求得后将进行QR分解 形成矩阵 QR算法，就是利用矩阵的QR分解，按上述递推法则构造矩阵序列的过程。只要A为非奇异矩阵，则由QR算法就完全确定。 四、实验内容 function qrsf(A,r) [Q,R]=qr(A); t=A(1,1) %temp A=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; end end n=size(A,1); for i=1:n format long g disp( [特征值λ,num2str(i),=,num2str( A(i,i) )] ); end %disp(); for i=1:n disp(特征值); format long g,A(i,i) end format long g,A,Q,R （2）改进后的QR算法的m文件 function qrsf(A,r) [Q,R]=qr(A); %t=A(1,1) %temp t(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; end end n=size(A,1); for i=1:n format long g disp( [特征值λ,num2str(i),=,num2str( A(i,i) )] ); end %disp(); for i=1:n disp(特征值); format long g,A(i,i) end format 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.0101500483978924 0.843107149855032 3.85805745594495 30.2886853458021 ans = 13.1723513981032 6.55187835191566 1.59565457314994 -0.390788045416488 -0.929096277752298 ans = 1.08279948406811e-007 1.25707571226224e-005 0.000615748354182652 0