数值分析大作业QR分解.doc

题目:设计求取实数矩阵A的所有特征值及其特征向量的数值算法，并以矩阵为例进行具体的求解计算。一、 算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。二、 算法设计：1、 原始实矩阵拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A进行拟上三角化，得到拟上三角化矩阵A记A(1)=A，并记A(r)的第r列到第n列的元素。对于r=1,2,n-2执行(1) 若全为零，则令A(r+1)= A(r)，转(5)；否则转(2)。(2) 计算令(3) 令，(4) 计算(5) 继续算法执行完后，就得到与原矩阵A相似的拟上三角矩阵A(n-1)。2、 拟上三角矩阵QR分解的求原矩阵的全部特征值记Ak是对拟上三角矩阵进行QR算法，产生的矩阵序列，A0是起始拟上三角矩阵，对于k=1,2,n-1执行：(1) 分解选取旋转矩阵，使，从而使第一列次对角元；再选取旋转矩阵，使，从而使第一列次对角元如此进行下去，最多经过n-1步，必然变为上三角矩阵，即必为正交矩阵，且满足(2) 计算(3) 上述两过程反复进行直到变为近似舒尔矩阵，对角线元素即为A的近似特征值。3、 带位移的QR方法求实矩阵全部特征值一般情况下，QR分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR算法所得矩阵序列的收敛。（A是原始矩阵的拟上三角矩阵）(1) 分解，其中即位移量，一般取A的某一特征值的近似值；实际计算通常取为的右下角元素，或取为右下角矩阵特征之中接近者。(2) 计算。(3) 重复过程（1）（2）直到变为近似舒尔矩阵，对角线元素即为A的近似特征值。4、 反幂法求实矩阵的特征向量记通过QR分解得到A的某一特征值的近似值为p，反幂法步骤如下：(1) 三角分解。(2) 对k=1，2，做当时，令 当时，解解求绝对值最大的元素令当或小于规定的误差界时，停止计算，即为所求的特征向量，即为对应特征值的更为精确的取值。三、 程序设计：1、 对生成的矩阵A进行拟上三角化利用hohd函数对A进行householder变换，得到A的拟上三角矩阵。2、 对拟上三角化后的矩阵进行QR分解和带位移的QR分解，求矩阵的全部特征值在绝对误差界为的条件下，利用qrfenjie函数对拟上三角矩阵进行QR分解，利用dp_qrfenjie函数进行带位移的QR分解，比较两者的收敛速度。3、 反幂法求更精确的特征值和特征向量利用vect函数和（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为。4、 输出相关结果。四、 源程序：1、 hohd函数functionA=hohd(a)n=length(a);for i=1:n-2 b=a(i+1:n,i); if b(1)=0 c=-sqrt(sum(b.2); else c=sqrt(sum(b.2); end rho=sqrt(2*c*(c-b(1); u1=(b-c*eye(n-i,1)/rho; H1=eye(n-i)-2*u1*u1; H=eye(n-i); H(i+1:n,i+1:n)=H1; a=H*a*(H);endA=a2、 qrfenjie函数function A=qrfenjie(a)n=length(a);A=a;for i=1:100 theta(n-1)=0; c(n-1)=0; s(n-1)=0; Q=1; for j=1:n-1 theta(j)=atan(A(j+1,j)/A(j,j); c(j)=cos(theta(j); s(j)=sin(theta(j); P=eye(n); P(j,j)=c(j); P(j+1,j)=-s(j); P(j,j+1)=s(j); P(j+1,j+1)=c(j); invP=eye(n); invP(j,j)=c(j); invP(j+1,j)=s(j); invP(j,j+1)=-s(j); invP(j+1,j+1)=c(j); Q=Q*invP; R=P*A; A=R; end A=R*Q; tor=max(abs(diag(A)-diag(a); if tor5*10-5 a=A; else break; endendi,Ak=A,lanmda=diag(A)3、 dp_qrfenjie函数function A=dp_qrfenjie(a)n=length(a);A=a;for i=1:100 A=a-a(n,n)*eye(n); theta(n-1)=0; c(n-1)=0; s(n-1)=0; Q=1; for j=1:n-1 theta(j)=atan(A(j+1,j)/A(j,j); c(j)=cos(theta(j); s(j)=sin(theta(j); P=eye(n); P(j,j)=c(j); P(j+1,j)=-s(j); P(j,j+1)=s(j); P(j+1,j+1)=c(j); invP=eye(n); invP(j,j)=c(j); invP(j+1,j)=s(j); invP(j,j+1)=-s(j); invP(j+1,j+1)=c(j); Q=Q*invP; R=P*A; A=R; end A=R*Q+a(n,n)*eye(n); tor=max(abs(diag(A)-diag(a); if tor5*10-5 a=A; else break; endendi,Ak=A,lanmda=diag(A)4、 vect函数function z=vect(a0,p)n=length(a0);Z=zeros(n,n+2);for x=1:na=a0-p(x)*eye(n);r=zeros(n,n);r(1,1:n)=a(1,1:n);l=eye(n);l(2:n,1)=a(2:n,1)/a(1,1);for i=2:n for j=i:n M=0; for k=1:i-1 M=l(i,k)*r(k,j)+M; end r(i,j)=a(i,j)-M; end if ip2 m=p1; else m=-p2; end z=y/m; tor=abs(m-mo); if tor clear,clc,a=2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1;hohd(a);得到原始实数矩阵的拟上三角矩阵A。A = 2.0000 -1.0000 0.0000 0.0000 -1.0000 1.0000 5.0000 -0.0000 0.0000 5.0000 -2.2000 -0.40000.0000 -0.0000 -0.4000 2.20002、 执行命令： clear,clc,A=2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2;qrfenjie(A);得到迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。i = 38Ak = -5.9068 0.0315 -0.0000 0.0000 0.0315 4.8972 -0.0000 -0.0000 0.0000 -0.0000 2.2138 0.00070.0000 0.0000 0.0007 1.7958lanmda = -5.9068 4.8972 2.2138 1.79583、 执行命令： clear,clc,A=2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2;dp_qrfenjie(A);得到用带位移的QR分解法所用的迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。i = 8Ak = -5.9068 -0.0056 -0.0000 -0.0000 -0.0056 4.8973 0.0000 0.0000 -0.0000 0.0000 1.7958 0.0000 0.0000 0.0000 0.0000 2.2138lanmda = -5.9068 4.8973 1.7958 2.21384、 执行命令： clear,clc,a=2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1,p=-5.9068 4.8972 2.2138 1.7958;vect(a,p);得到Z矩阵，Z矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。Z = Columns 1 through 4 5.000000000000000 -5.906847942119164 0.112419521552950 1.000000000000000 5.000000000000000 4.897302326740487 0.345148654584838 0.505132030845888 5.000000000000000 2.213757601733807 -0.282974649968150 -0.697610774659887 5.000000000000000 1.795788013644870 1.000000000000000 -0.324049010823671 Columns 5 through 6 0.675655845769651 -0.888884062644966 0.510541849590699 1.000000000000000 1