方法求矩阵全部特征值

1、数值分析课程设计Q昉法求矩阵全部特征值问题复述用QR算法求矩阵特征值：2 3 4 5 6一6 2 114 4 5 6 7(i ) (A) = 2 3 1(ii ) H = 0 3 6 7 81111002890 0 0 1 0_ 要求：(1)根据QR算法原理编制求(i)与(ii )中矩阵全部特征值的程序并输出计算结果(要求误差10,)(2)直接用现有的数学软件求(i), (ii )的全部特征值，并与(1)的结果比 较。问题分析QR方法是目前公认为计算矩阵全部特征值的最有效的方法。它适用于任一 种实矩。QR方法的原理是利用矩阵的正交分解产生一个与矩阵A相似的矩阵迭代序列，这个序列将收敛于一个上三

2、角阵或拟上三角阵，从而求得原矩阵A的全部特征值。QFg一个迭代算法，每一步迭代都要进行 QR分解，再作逆序的矩阵乘法。要是对矩阵A直接用QR方法，计算量太大，效率不高。因此，一个实际的 QR方法计算过程包括两个步骤，首先要对矩阵A用初等相似变换约化为上Hessenberg矩阵，再用 QR7J法求上Hessenberg矩阵的全部特征值。 示意如下：用Householder阵作正交相似变换第一步 A寸）吧2T 上HessenbergB =X1格一止 c用QR变换产生迭代序列，迭 代计算Bk =QkRk弟少B -.Bk = RkQk因此常常用平面旋对B矩阵的约化只需将每列次对角线上的元素约化为 转阵

3、（Givens变换阵）来进行约化。一、Q叫法原理及收敛性 设AWRn也已实现了 QR分解，记A = A = Q1R1其中是正交阵，是上三角阵。因为 QiT=Q/，用对作正交相似变换有Q1 A1Q1 - A2可改写为A2=QAQi=QQiRQ =RQi显然只是的QR分解因子阵的逆序相乘，而且与原矩阵有相同的特征值。对矩阵再重复以上过程并继续下去，可以得到一个与原矩阵A有相同特征值的矩阵序列。其过程如下：对作正交分解A=Q1R作矩阵A2=RQi,A=A对作正交分解A2 -Q2R2作矩阵 A=R2Q2,重复以上过程可得一般的形式为对作正交分解Ak =QkRk (A = A)构成矩阵序列A = RQk

4、(k=1,2。A从矩阵A开始得到一个矩阵序列 一 AM , 2, 3, k ,这个矩阵序列中每一个矩阵都与原矩阵 A( = A)相似，即都有与A相同的特征值。这个矩阵序列AJ在实质上收敛于依次以 储,”,%为对角元的上三角阵。具体可以表示为% x x +3A,2:lkmxI，印其中aii(k)(k.) i a j时aj(k) t 0 (kT 吟i>j时ajk)的极限不一定存在二、用正交相似变换约化矩阵为上 Hessenberg阵用Householder变换可以将一个向量指定的某个分量以下的各分量变为 0。我 们只要求消掉A的次对角线以下的元素，即将 A约化为上Hessenberg阵。为了

5、 使变换前后矩阵的特征值不变，需要用Householder矩阵对A作相似变换，即用 正交阵同时左乘和右乘A时，原来已变为0的元素不再改变。若设是Householder 矩阵，用它对A的第一列元素的变换示意如下:xH 1 AHxxx依次对A的各列进行类似的变换，一共要进行 n-2次变换，最终可以得到一 个与原矩阵A有相同特征值的上Hessenberg阵。x x+mx J-:H n/H nH 2H lAH iH 2 H nH n/ = 工 三 ：*二 xIx以上约化A为上Hessenberg阵的过程可以用一系列Householder矩阵来实现。H kCk = -a ke其中，R/十算公式为其中，对

6、于每一个有Ck = (ak 也Jk)二 k = sgn(a(k)、T.nA,ank ) Rn(k)(k)k 1,k )( -(aiki#1Uk =Ck . 3,，Ck = Ck/|Ck|L/Ck-)2)1/2'k ,j(akiJk),Rk = I一1UkUkT.经过n -2步约化就可得到一个上 Hessenberg阵（A的第n-1列不需要约化）An-1= Hn_/“ H2H1AH1H2Hn_2(记为)B -Va11 xx一C (2)“1a22(3).一灯"22 a33n -2(n -1)(n-1)(n-1)(n-1) n(n-1)an;)、Hessenberg 阵的 QRB法

7、设矩阵AWRn殉，其特征值都是实数。若已用 Householder变换约化为上Hessenberg 阵PAPT对已得到的上Hessenberg阵可用QR变换，经过迭代过程约化为上三角形矩阵以 求出A的特征值。只要A的特征值是实数，将B约化为上三角形矩阵总是可以做到的。由于B矩阵结构上的特点，对B矩阵的约化只需将每列次对角线上的元素约 化为0.可用平面旋转阵（Givens变换阵）来进行约化。n阶方阵R（i, j）为平面旋转阵cos 二(j >i)1R(i, j)=cos?Rj非奇异，显然有是一个平面旋转阵。 有以下几种作用：det(Rj) -1还是一个非对称的正交阵，有RjRij=I ,R

8、j）也1、左乘向量x =（x1, x2，xn）T只改变x的第i个和第j个分量。现构造对 x作变换使Rjx的结果将x的第j个分量约化为0yi = xi cos。 xj sin令 y = Rjx,有 j yj = -xi sin 9 + xj cos 9yk = xk(k =1,2, ,n;k ; i, j)调整角可使yj =0。若记C =cos"S =sina ,按下式选取C = xi /、x i2 x j2 l x i / r , S = x j / v xi2 x j2 l x j / r于是 J 丫 = Cxi + Sxj = r =4x： ”2、 y j = - Sk + Cx

9、 j =0有 Rjx = ("，x-, r, xi由，xj.0, xj由，xn)T。2、对非零的n维向量x连续左乘Ri,i + ,Ri,(i42)， ,可将x的第i+1到第n个分量都约化为零；即RinRi(nf Ri(i 2)Ri(i ” = (Xi，x-,ri,0,，0)T其中ri222XiXi 1Xn3 、用对矩阵A作变换得到的结论是RjT左乘A只改变A的第i,j行；RjT右乘A只改变A的第i,j歹1;用对A作正交相似变换RjTARj改变了 A的i行和j行以及i列和j歹I。用一系列R（i，j）连续左乘矩阵A,可以将矩阵A化为上三角阵。数据结构描述主要数据成员说明double An

10、n存放矩阵Adouble Qnn,存放Q的解式的正交矩阵Qdouble Rnn存放Q的解式的上三角阵Rdouble pnnGivens 矩阵 pdouble InnN阶单位阵double Vnn存放Q矩阵的转置double Tnn初等反射阵Tdouble eps精度double max最大值double det存放行列式的值int count存放迭代次数主要函数成员说明double Det(double Lnn)用高斯列主元方法求行列式int Non_singularMatrix(double Lnn)判断是否是非奇异矩阵void Disp(double Hnn)输出矩阵int IsZero(

11、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)判断是否满足QRS法条件，满足 则进行QR7T法求特征值算法的描述（流程图）主函数流程图:源程序C程序为：#include<stdio.h>#include<math.h&

12、gt;#define n 3#define eps 1E-5double Det(double Lnn)/用高斯列主元方法求行列式double det=1,t;for(int k=0;k<n-1;k+)double max=Lkk;int Ik=k;for(int j=k+1;j<n;j+)if(max<Ljk) max=Ljk; Ik=j;if(max=0) return 0;if(Ik!=k)for(int j=k;j<n;j+)t=LIkj;LIkj=Lkj;Lkj=t;det*=-1;for(int i=k+1;i<n;i+) t=Lik/Lkk;Lik寸

13、for(int j=k+1;j<n;j+) Lij=Lij-t*Lkj;det*=Lkk;if(Ln-1n-1=0) return 0;elseprintf(" 矩阵 A叩 的行列式为：.5fn",det*Ln-1n-1);return det*Ln-1n-1;int Non_singularMatrix(double Lnn)/判断是否是非奇异矩阵printf(" 判断矩阵A叩的行列式是否为0?n");if(Det(L)!=0) return 1;return 0;void Disp(double Hnn)/ 输出矩阵for(int i=0;i&

14、lt;n;i+)for(int j=0;j<n;j+)printf("%.6f ",Hij);printf("n");int IsZero(double a口,int j)判断数组是否全为 0for(int i=0;i<j;i+)if(ai!=0)return 0;return 1;int sgn(double y)/符号函数if(y>0) return 1;else if(y=0) return 0;else return -1;void Hessenberg(double Ann)/将矩阵化为上 Hessenberg 矩阵print

15、f("原矩阵为:n");Disp(A);/输出原矩阵doubleTnn,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;k<n-2;k+)for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)if(i=j) Iij=1;else Iij=0;/ 定义单位阵I叩 double max=fabs(Ak+1k);/求最大值for(i=0;i<n-k-1;i+)if(max<fabs(Ai+k+1k) max=fabs(Ai+k+1k);for(i=

16、0;i<n-k-1;i+)ci=Ai+k+1k/max;/标准化数组if(IsZero(c,n-k-1)判断数组是否全为0continue;/数组为0,则这一步不需要约化 for(i=0,t=0.0;i<n-k-1;i+)t+=ci*ci;vk=sgn(Ak+1k)*sqrt(t);/u0=c0+vk;for(j=1;j<n-k-1;j+)uj=cj;/ 求矩阵 u叩 w=vk*(c0+vk);for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)Rij=Iij-ui*uj/w;for(i=0;i<n;i+)for(j=0;j<n

17、;j+)if(i=j) Tij=1;else Tij=0;/ 定义矩阵T为单位阵)for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)Ti+k+1j+k+1=Rij;/初等反射阵 Tfor(i=0;i<n;i+) /矩P$ T左乘矩阵 Afor(j=0;j<n;j+)for(m=0,s=0.0;m<n;m+) s+=Tim*Amj;Bij=s;)for(i=0;i<n;i+) /矩P$ T右乘矩阵 Afor(j=0;j<n;j+)for(m=0,s=0.0;m<n;m+) s+=Bim*Tmj;Cij=s;)for(i=0

18、;i<n;i+)for(j=0;j<n;j+) Aij=Cij;)printf("原矩阵化成上 Hessenberg 矩阵后为:n");Disp(A);)int IsHessenberg(double Enn)/判断是否是上 Hessenberg 矩阵for(int i=2;i<n;i+)for(int j=0;j+1<i;j+)if(Eij!=0)return 0;return 1;算法求特征值int count=1;void QRAlgorithm(double Ann)/QRint i,j,k,m;double pnn,Qnn,Rnn,Fnn,

19、Vnn,c,s,t,y,max;for(i=0;i<n;i+)/ 赋 Q为单位阵for(int j=0;j<n;j+)if(i!=j) Qij=0；else Qij=1;for(m=0;m<n-1;m+) 对矩阵 A 进行 QR分解for(i=0;i<n;i+)/ 赋 p 为单位阵for(j=0;j<n;j+)if(i!=j) pij=0;else pij=1;c=Amm/sqrt(Amm*Amm+Am+1m*Am+1m);s=Am+1m/sqrt(Amm*Amm+Am+1m*Am+1m);为 Givenspmm=c;pmm+1=s;pm+1m=-s;pm+1m+

20、1=c;/p矩阵for(i=0;i<n;i+)/p 矩阵左乘 W巨阵,p矩阵左乘Q矩阵 for(j=0;j<n;j+)t=0;y=0;for(k=0;k<n;k+) t+=pik*Akj; y+=pik*Qkj;Rij=t;Fij=y;for(i=0;i<n;i+)/A,R口赋新值for(j=0;j<n;j+)Aij=Ri皿 Q皿产Fi皿 for(i=0;i<n;i+)/Q 转置for(j=0;j<n;j+)Vij=Qji；for(i=0;i<n;i+)for(j=0;j<n;j+)Q皿产V皿for(i=0;i<n;i+)/矩P$ A

21、为矩阵R和矩阵Q的乘积for(j=0卜n;j+)t=0;for(k=0;k<n;k+) t+=Rik*Qkj;Aij=t；count+;max=fabs(A10);/求A矩阵下三角的绝对值最大值for(i=1;i<n;i+)for(j=0;j<i;j+)if(fabs(Aij)>max)max=fabs(Aij);if(max<eps)/满足精度条件，输出Q矩阵，R矩阵以及特征值printf(" 在精度.1e下,QR算法迭代次数为：次坨输出矩阵A%d:n",eps,count-1,count)；Disp(A)；/输出 A矩阵printf(&qu

22、ot; 输出矩阵 Q:n");Disp(Q);printf("输出矩阵 R:n");Disp(R);printf(" 输出A矩阵的全部特征值:");for(i=0;i<n;i+)if(i%3=0) printf("n");printf("a%d=%.6忖”,i+1,Aii);printf("n");return;QRAlgorithm(A);判断是否满足QR算法条件，满足则void SeekEigenvalue(double Ann)/进行QR方法求特征值double Znn;for(in

23、t i=0;i<n;i+)for(int j=0;j<n;j+)Zij=A皿;printf(" 判断矩阵A叩 是否是非奇异矩阵？n");if(Non_singularMatrix(Z)=0)printf("矩阵A叩不是非奇异矩阵！不符合分解条件！n");return;elseprintf("矩阵A叩 是非奇异矩阵！符合分解条件!n");printf(" 判断矩阵A 是否是上Hessenberg矩阵？n");if(IsHessenberg(A)=0)printf(" 不是上Hessenberg矩

24、阵！将其化为上Hessenberg矩阵.n");Hessenberg(A);/将矩阵A转化为上Hessenberg矩阵elseprintf(" 是上Hessenberg矩阵！不需要将其再化为上 Hessenberg矩 阵!n");printf("QR 方法求全部特征值:n");QRAlgorithm(A);/ 用 QRTJ法求特征值void main()/*doubleA155=2,3,4,5,6,4,4,5,6,7,03678,0Q2,8,9,0Q0,1,0;用此组数据时#define n 5*/double A133=6,2,1,2,3,1

25、,1,1,1;SeekEigenvalue(A1);MATLAB?:»坛陶41岛3.1乩L 11 ；» 8=忆3#%的6&%5, 6"山施6门萨；心也2/歌9山曼为1,0;» ei£(A)» eig (H)测试数据及运行结果测试数据为：6 2 1A= 2 3 1 H1 1 1程序结果为： 对矩阵A2345644567036780028900010_阵白12口:输出口门口矩阵的全部特征值:al=7.287992a2=2.133074a3=0.578933Press any key to continue3.&QQQ&am

26、p;Q1.0000000.0000021.QQQ0QQ Q.0000001.0000000.0000000.2000000.6000002.QQQQQQ1.0000006.000000 -2.2360680.0000001.000000 -0.000002 0.000000-2.2360683.400000 0.200000-0.000000-0.QQQ&QQ 1.000000-0.000000 0.0QQQ&Q 0.578933-0.000000 0.000000 0.5789337.287992 -0.000004 -0.000004 2.133074 -0.000000

27、0.000000 输出矩阵Q门口;输出矩阵R门口; 7.287992-0.000018&.&&&&&&2.133074-0.0000000.000000,|n|?-LLL 可富列异丕SS其IIF:9?rHe将?0r 牛r上 0 AM e JFUKruz 0 . .T b 4/ 矩否腼衿en化 SB-日AE口行育AEer00 LUJV - - - b . n 2 口郸 口日TH口seA T.J fi s -阵徐为00 矩断阵A矩上阵00阵曹TH矩00 判 矩判不原6.原矩阵化成上Hew wen berg矩阵后为：?,阵lW *C ： Vs

28、er s * £三h'D e skt op I> 色b口云上口由法求特 fiF 值-exe*QR无法求全邰特征值;苣精度10-眄5下,QR算法迭代次数为：11次 输出矩阵口12口:对矩阵HIW：VEers£shDesktopD>el>iig：tQ®法求特征值.eze |n| x判断库审是否是非向昌矩阵？前断矩阵口门 口前强底票否为0?恒制江是非奇异售冷符合分解,像M断矩阵口 口口是否是上H5cnberg矩阵？是上H5EnbcFg矩除不需要将其再化为上H-5cnbcrg矩阵，帆无法求全部特征看年精度10-眄5下,QR算法迭代次数为：22次输出矩阵由23口:13.17234711.2224361.383301-12.287651-2.1843