题目:设计求取n n ⨯实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵20010-1-24A=0-2131431⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭为例进行具体的求解计算。

一、 算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。

其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。

二、 算法设计：1、 原始实矩阵A Rn n⨯∈拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)rij a i n j r r n ==+。

对于r=1,2,…,n -2执行(1) 若()(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。

(2) 计算 令()2()()1,1,1,sgn(0,sgn()=1)r r r r r rr r r r r c aa a ρ+++=-=，（若则取(3) 令-0=r n r u u ⎛⎫ ⎪⎝⎭，()()()-1,2,1(,,...,)r r r Tn r r r r r r nr u a c a a ρ++=- (4) 计算(r)(r)(r)Tn-r r+1,r r r+2,r nr r Tn-r T n-r n-r n-r n-r r+1r 1u =(a -c ,a ,...,a )ρI H =I -2uu =H H =I -2u u A =HA H⎛⎫⎪⎝⎭(5) 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。

2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)A必然变为上三角矩阵k-1R ，即k-1n-121k-1-1=-1-1-1k-1n-12112n-1R =P P P A =P P P P P PQ （）k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q(2) 计算k k-1k-1A =R Q(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。

3、 带位移的QR 方法求实矩阵A Rn n⨯∈全部特征值一般情况下，QR 分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR 算法所得矩阵序列k A 的收敛。

0A =A （A 是原始矩阵的拟上三角矩阵）(1) 分解k-1-1k-1k-1A -u I=R k Q ，其中-1u k 即位移量，一般取A 的某一特征值的近似值；实际计算通常取为k-1A 的右下角元素(k-1)n,n a ，或取为右下角22⨯矩阵特征之中接近(k-1)n,n a 者。

(2) 计算k k-1k-1-1A =+R u I k Q 。

(3) 重复过程（1）（2）直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。

4、 反幂法求实矩阵A Rn n⨯∈的特征向量记通过QR 分解得到A 的某一特征值的近似值为p ，反幂法步骤如下： (1) 三角分解A -pI =LR 。

(2) 对k=1，2，…，做①当k=1时，令Tu =(1 11) 当k 1≠时，解-1Lu=z k②解Ry =u k③求-1y k 绝对值最大的元素k m ④令k k k z =y /m⑤当k k-1m -m 或k k-1z -z 小于规定的误差界时，停止计算，k z 即为所求的特征向量，k p +1/m 即为对应特征值的更为精确的取值。

三、 程序设计：1、 对生成的矩阵A 进行拟上三角化利用hohd 函数对A 进行householder 变换，得到A 的拟上三角矩阵。

2、 对拟上三角化后的矩阵进行QR 分解和带位移的QR 分解，求矩阵的全部特征值在绝对误差界为-41102⨯的条件下，利用qrfenjie 函数对拟上三角矩阵进行QR 分解，利用dp_qrfenjie 函数进行带位移的QR 分解，比较两者的收敛速度。

3、 反幂法求更精确的特征值和特征向量利用vect 函数和（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为-61102⨯。

4、 输出相关结果。

四、 源程序：1、 hohd 函数function[A]=hohd(a) n=length(a); for i=1:n-2b=a(i+1:n,i); if b(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=a2、qrfenjie函数function A=qrfenjie(a)n=length(a);A=a;for i=1:100theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;for j=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(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)));if tor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'3、dp_qrfenjie函数function A=dp_qrfenjie(a)n=length(a);A=a;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-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(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+a(n,n)*eye(n);tor=max(abs(diag(A)-diag(a)));if tor>5*10^-5a=A;elsebreak;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(2:n,1)=a(2:n,1)/a(1,1);for i=2:nfor j=i:nM=0;for k=1:i-1M=l(i,k)*r(k,j)+M;endr(i,j)=a(i,j)-M;endif i<nfor j=i+1:nN=0;for k=1:i-1N=l(j,k)*r(k,i)+N;endl(j,i)=(a(j,i)-N)/r(i,i);endelsebreak;endendR=r;L=l;u=ones(n,1);mo=0;invr=inv(R);invl=inv(L);for i=1:100y=invr*u;p1=max(y);p2=max(-y);if p1>p2m=p1;elsem=-p2;endz=y/m;tor=abs(m-mo);if tor<5*10^-7break;elsemo=m;endu=invl*z;endZ(x,3:n+2)=z';Z(x,1)=i;Z(x,2)=p(x)+1/m;endZ五、运行结果：1、执行命令：>> 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.00000.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.00000.0315 4.8972 -0.0000 -0.00000.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.00000.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矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。