1、 作矩阵 ，重复以上过程可得一般的形式为对作正交分解 构成矩阵序列 （k=1,2） A 从矩阵A开始得到一个矩阵序列 这个矩阵序列中每一个矩阵都与原矩阵相似，即都有与A相同的特征值。这个矩阵序列在实质上收敛于依次以为对角元的上三角阵。具体可以表示为其中 的极限不一定存在。二、 用正交相似变换约化矩阵为上Hessenberg阵 用Householder变换可以将一个向量指定的某个分量以下的各分量变为0。我们只要求消掉A的次对角线以下的元素，即将A约化为上Hessenberg阵。为了使变换前后矩阵的特征值不变，需要用Householder矩阵对A作相似变换，即用正交阵同时左乘和右乘A时，原来已变为
2、0的元素不再改变。若设是Householder矩阵，用它对A的第一列元素的变换示意如下： 依次对A的各列进行类似的变换，一共要进行次变换，最终可以得到一个与原矩阵A有相同特征值的上Hessenberg阵。 以上约化A为上Hessenberg阵的过程可以用一系列Householder矩阵来实现。其中，对于每一个有经过步约化就可得到一个上Hessenberg阵（A的第列不需要约化） 三、 Hessenberg阵的QR算法 设矩阵，其特征值都是实数。若已用Householder变换约化为上Hessenberg阵对已得到的上Hessenberg阵可用QR变换，经过迭代过程约化为上三角形矩阵以求出A的特
3、征值。只要A的特征值是实数，将B约化为上三角形矩阵总是可以做到的。 由于B矩阵结构上的特点，对B矩阵的约化只需将每列次对角线上的元素约化为0.可用平面旋转阵（Givens变换阵）来进行约化。n阶方阵为平面旋转阵 。还是一个非对称的正交阵，有，也是一个平面旋转阵。 有以下几种作用： 1、左乘向量只改变x 的第i个和第j个分量。现构造对x作变换使的结果将x的第j个分量约化为0。令，有调整角可使。若记，按下式选取于是有。 2、对非零的n维向量x连续左乘，可将x的第i+1到第n个分量都约化为零；即 其中 3、用对矩阵A作变换得到的结论是 左乘A只改变A的第i,j 行； 右乘A只改变A的第i,j列； 用
4、对A作正交相似变换改变了A的i行和j行以及i列和j列。用一系列连续左乘矩阵A，可以将矩阵A化为上三角阵。 数据结构描述主要数据成员说明double Ann 存放矩阵Adouble Qnn,存放QR分解式的正交矩阵Qdouble Rnn存放QR分解式的上三角阵Rdouble pnnGivens矩阵pdouble InnN阶单位阵double Vnn存放Q矩阵的转置double Tnn初等反射阵Tdouble eps精度double max最大值double det存放行列式的值int count存放迭代次数主要函数成员double Det（double Lnn）用高斯列主元方法求行列式int N
5、on_singularMatrix（double Lnn）判断是否是非奇异矩阵void Disp（double Hnn）输出矩阵int IsZero（double a,int j）判断数组是否全为0int sgn（double y）符号函数void Hessenberg（double Ann）将矩阵化为上Hessenberg矩阵int IsHessenberg（double Enn）判断是否是上Hessenberg矩阵void QRAlgorithm（double Ann）QR算法求特征值void SeekEigenvalue（double Ann）判断是否满足QR算法条件，满足则进行QR方法
6、求特征值算法的描述（流程图） 源程序C程序为：#includemath.h#define n 3#define eps 1E-5double Det（double Lnn）/用高斯列主元方法求行列式double det=1,t;for（int k=0;kn-1;k+）double max=Lkk;int Ik=k;for（int j=k+1;jn;j+）if（maxLjk）max=Ljk;Ik=j;if（max=0） return 0;if（Ik!=k）for（int j=k;j+）t=LIkj;LIkj=Lkj;Lkj=t;det*=-1;for（int i=k+1;i0） return 1
7、;else if（y=0） return 0;else return -1;void Hessenberg（double Ann）/将矩阵化为上Hessenberg矩阵原矩阵为:Disp（A）;/输出原矩阵double Tnn,Bnn,Cnn,cn-1,vn-1,un-1,Rn-1n-1,In-1n-1,t,w,s;int i,j,k,m;for（k=0;n-2;for（i=0;n-k-1;for（j=0;if（i=j） Iij=1;else Iij=0;/定义单位阵Idouble max=fabs（Ak+1k）;/求最大值fabs（Ai+k+1k）max=fabs（Ai+k+1k）;ci=A
8、i+k+1k/max;/标准化数组if（IsZero（c,n-k-1）/判断数组是否全为0continue;/数组为0，则这一步不需要约化for（i=0,t=0.0;t+=ci*ci;vk=sgn（Ak+1k）*sqrt（t）;/u0=c0+vk;for（j=1;uj=cj;/求矩阵uw=vk*（c0+vk）;Rij=Iij-ui*uj/w;if（i=j） Tij=1;else Tij=0;/定义矩阵T为单位阵Ti+k+1j+k+1=Rij;/初等反射阵Ti+） /矩阵T左乘矩阵Afor（m=0,s=0.0;mm+）s+=Tim*Amj;Bij=s;i+） /矩阵T右乘矩阵As+=Bim*Tm
9、j;Cij=s;Aij=Cij;原矩阵化成上Hessenberg矩阵后为:int IsHessenberg（double Enn）/判断是否是上Hessenberg矩阵for（int i=2;j+1max）max=fabs（Aij）;if（maxeps）/满足精度条件，输出Q矩阵，R矩阵以及特征值在精度%.1e下,QR算法迭代次数为:%d次n输出矩阵A%d:,eps,count-1,count）;/输出A矩阵输出矩阵Q:Disp（Q）;输出矩阵R:Disp（R）;输出A矩阵的全部特征值:if（i%3=0） printf（a%d=%.6ft,i+1,Aii）;return;QRAlgorithm
10、（A）;void SeekEigenvalue（double Ann）/判断是否满足QR算法条件，满足则进行QR方法求特征值double Znn;Zij=Aij;判断矩阵A是否是非奇异矩阵?if（Non_singularMatrix（Z）=0）矩阵A不是非奇异矩阵!不符合分解条件!return;else矩阵A是非奇异矩阵!符合分解条件!判断矩阵A是否是上Hessenberg矩阵?if（IsHessenberg（A）=0）不是上Hessenberg矩阵!将其化为上Hessenberg矩阵.nHessenberg（A）;/将矩阵A转化为上Hessenberg矩阵else 是上Hessenberg矩
11、阵!不需要将其再化为上Hessenberg矩阵!QR方法求全部特征值:/用QR方法求特征值void main（）/*double A155=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;用此组数据时#define n 5*/double A133=6,2,1,2,3,1,1,1,1;SeekEigenvalue（A1）;MATLAB程序为： 测试数据及运行结果测试数据为:程序结果为：对矩阵A对矩阵HMATLAB结果为： 程序中的标识符的说明（类型、含义等）见程序。完成情况说明： 通过编程，完成了对矩阵是否符合QR分解条件的判断，对于不是上Hessenberg阵的矩阵，将其约化为上Hessenberg阵，然后通过QR算法进行递归，最后求得满足精度条件的矩阵A的全部特征值。并且，与MATLAB软件求解的结果进行比较，两者的结果是一致的，这也证明了程序的正确性。
