资源描述：
《数值分析大作业qr分解》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、题目:设计求取实数矩阵A的所有特征值及其特征向量的数值算法，并以矩阵为例进行具体的求解计算。一、算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。二、算法设计：1、原始实矩阵拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A进行拟上三角化，得到拟上三角化矩阵A’记A(1)=A，并记A(r)的第r列到第n列的元素。对于r=1,2,…,n-2执行(1)若全为零，则令A
2、(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、近似特征值。1、带位移的QR方法求实矩阵全部特征值一般情况下，QR分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR算法所得矩阵序列的收敛。（A是原始矩阵的拟上三角矩阵）(1)分解，其中即位移量，一般取A的某一特征值的近似值；实际计算通常取为的右下角元素，或取为右下角矩阵特征之中接近者。(2)计算。(3)重复过程（1）（2）直到变为近似舒尔矩阵，对角线元素即为A的近似特征值。2、反幂法求实矩阵的特征向量记通过QR分解得到A的某一特征值的近似值为p，反幂法步骤如下：(1)三角分解。(2)对k=1，2，…，做①
4、当时，令当时，解②解③求绝对值最大的元素④令⑤当或小于规定的误差界时，停止计算，即为所求的特征向量，即为对应特征值的更为精确的取值。一、程序设计：1、对生成的矩阵A进行拟上三角化利用hohd函数对A进行householder变换，得到A的拟上三角矩阵。2、对拟上三角化后的矩阵进行QR分解和带位移的QR分解，求矩阵的全部特征值在绝对误差界为的条件下，利用qrfenjie函数对拟上三角矩阵进行QR分解，利用dp_qrfenjie函数进行带位移的QR分解，比较两者的收敛速度。3、反幂法求更精确的特征值和特征向量利用vect函数和（2）得到的特征值
5、的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为。4、输出相关结果。二、源程序：1、hohd函数function[A]=hohd(a)n=length(a);fori=1:n-2b=a(i+1:n,i);ifb(1)>=0c=-sqrt(sum(b.^2));elsec=sqrt(sum(b.^2));endrho=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
6、=a1、qrfenjie函数functionA=qrfenjie(a)n=length(a);A=a;fori=1:100theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;forj=1:n-1theta(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(
7、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;endA=R*Q;tor=max(abs(diag(A)-diag(a)));iftor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'2、dp_qrfenjie函数functionA=dp_qrfenjie(a)n=length(a);A=a;fori=1:100A=a-a(n,n)*eye(n);theta(n-1)=0;c(n-1)=0;s(n-1)
8、=0;Q=1;forj=1:n-1theta(j)=atan(A(j+1,j)/A(j,j));c(j)=cos(theta(j));s(j)=sin(theta(j));P=e
