方法求矩阵全部特征值0001

1、数值分析计QF方法求矩阵全部特征值问题复述用QR算法求矩阵特征值:23456162们44567231(ii ) H =0367811 一0028900010 一(i)(A)二要求:（1）根据QR算法原理编制求(i )与）中矩阵全部特征值的程序并输出计算结果（要求误差门0冷（2）直接用现有的数学软件求（i），（ii ）的全部特征值，并与（1）的结果比 较。问题分析QR方法是目前公认为计算矩阵全部特征值的最有效的方法。它适用于任一 种实矩。QR方法的原理是利用矩阵的正交分解产生一个与矩阵 A相似的矩阵迭 代序列，这个序列将收敛于一个上三角阵或拟上三角阵， 从而求得原矩阵A的全 部特征值。QF是一个

2、迭代算法，每一步迭代都要进行 QR分解，再作逆序的矩阵乘法。要是对矩阵A直接用QR方法，计算量太大，效率不高。因此，一个实际的QR方法计算过程包括两个步骤，首先要对矩阵A用初等相似变换约化为上Hessenberg矩阵，再用 QR方法求上Hessenberg矩阵的全部特征值。示意如下：第一步A用Householder阵作正交相似变换不需迭代，对A的每一列计算一次 、上 Hesse nberg车 B =,L x:Iax x第二步B用QR变换产生迭代序列，迭 代计算Bk 二 Qk Rk、YtBk = RkQk0。因此常常用平面旋对B矩阵的约化只需将每列次对角线上的元素约化为 转阵（Give ns变换

3、阵）来进行约化。一、QR方法原理及收敛性 设A Rn n已实现了 QR分解，记其中是正交阵，是上三角阵。因为 q1t，用对作正交相似变换有可改写为 A? =QAQ1 =QQRQ =RQ显然只是的QR分解因子阵的逆序相乘，而且与原矩阵有相同的特征值。对矩阵 再重复以上过程并继续下去，可以得到一个与原矩阵A有相同特征值的矩阵序列。其过程如下：记A =A对作正交分解A hQjR作矩阵A? =RQi，A=A对作正交分解AQ?R?作矩阵A=R2Q2，AA重复以上过程可得一般的形式为对作正交分解Ak二QkRk (A二A)构成矩阵序列A. RQk (k=1,2)Ak 1 A从矩阵A开始得到一个矩阵序列A A

4、 123k ,这个矩阵序列中每一个矩阵都与原矩阵 A( = A)相似，即都有与A相同的特征值。这个矩阵序列AJ在实质上收敛于依次以 a匕,n为对角元的上三角阵。 具体可以表示为人Xx|时=X其中 a'1")(k )i j时a/" ； 0 (kr )i j时aj的极限不一定存在二、用正交相似变换约化矩阵为上Hessenberg阵用Householder变换可以将一个向量指定的某个分量以下的各分量变为 0。我 们只要求消掉A的次对角线以下的元素，即将 A约化为上Hessenberg阵。为了 使变换前后矩阵的特征值不变，需要用Householder矩阵对A作相似变换，即用

5、 正交阵同时左乘和右乘A时，原来已变为0的元素不再改变。若设是Householder 矩阵，用它对A的第一列元素的变换示意如下:xlH 1 AH依次对A的各列进行类似的变换，一共要进行 n-2次变换，最终可以得到一 个与原矩阵A有相同特征值的上Hessenberg阵。nz（k）/ P z （k）S 二 sgn（ak1,k ）（佝丘i 土岀Uk =Ck 0,x *k（ak1,k（k）,Rk = I - 1 _1UkUkT.其中，对于每一个有H nH 2H 2H 1AH 出 2 H n_3H n _2以上约化A为上Hessenberg阵的过程可以用一系列Householder矩阵来实现。H kCk

6、 = -二2其中，Rk计算公式为 CkaJ,rZ Ck = Ck/|Ck|L 仪也）2）1/2,经过n -2步约化就可得到一个上 Hessenberg阵（A的第n -1列不需要约化）代一1 二 Hn_2 H2H1AH1H2 Hn_2（记为）BRin Ri (n 书Ri (i 2) Ri(i 1)X =(X1 , Xj 二，仃，0 ,0)知-C1Xa22a33n -2a（n -1）（n）（n_1） a（n-1） n（n1）X（n -1）ann三、Hessenberg阵的QR算法设矩阵A Rn n，其特征值都是实数。若已用Householder变换约化为上Hessenberg 阵XpaptBxHe

7、ssenberg阵可用QR变换，经过迭代过程约化为上三角形矩阵以对已得到的上求出A的特征值。只要A的特征值是实数，将B约化为上三角形矩阵总是可以做到的。由于B矩阵结构上的特点，对B矩阵的约化只需将每列次对角线上的元素约 化为0.可用平面旋转阵(Give ns变换阵)来进行约化。n阶方阵R(i, j)为平面旋转阵R(i, j)二(j >i)-sin vcostdet(Rj) = 1还是一个非对称的正交阵，有RjRj(T)=l，尺厂也Rj非奇异，显然有是一个平面旋转阵。有以下几种作用：1 、左乘向量(x1, x2/ ,xn)T只改变x的第i个和第j个分量。现构造对x作变换使RjX的结果将x的

8、第j个分量约化为0。令 y =RjX,yi = Xj cos 日 + xj sin 日有 yj = -Xj sin 寸xj cos 寸yk = Xk(k =1,2,n;k = i, j)调整角可使yj = 0。若记C=cos3S=s in二，按下式选取C 二 Xi / i x i$ x j 匚 x i / r , S = x j / v Xi $ x j x j / rV、二 CXi SXj = r 二 Xi2 Xj2yj = _SXi 十 Cx j = 0有 RjX =(X1，x.,r,Xi 1,，Xjj,O,Xj .1,，Xn)T。2、对非零的n维向量x连续左乘R宀，R(七)，,Rin，可

9、将x的第i+1到第n个分量都约化为零；即其中3、用对矩阵A作变换得到的结论是叮左乘A只改变A的第i,j行;叮右乘A只改变A的第i,j列；用对A作正交相似变换R/ARj改变了 A的i行和j行以及i列和j列用一系列R（i，j）连续左乘矩阵A,可以将矩阵A化为上三角阵。数据结构描述主要数据成员说明double An n存放矩阵Adouble Qnn,存放QR分解式的正交矩阵Qdouble Rn n存放QR分解式的上三角阵Rdouble pn nGive ns 矩阵 pdouble lnnN阶单位阵double Vn n存放Q矩阵的转置double Tnn初等反射阵Tdouble eps精度doubl

10、e max最大值double det存放行列式的值int count存放迭代次数主要函数成员说明double Det(double Ln n)用咼斯列主元方法求行列式int Non si ngularMatrix(double Lnn)判断是否是非奇异矩阵void Disp(double H n n)输出矩阵int IsZero(double a,i nt j)判断数组是否全为0int sgn( double y)符号函数void Hesse nberg(double An n)将矩阵化为上Hessenberg矩阵int IsHesse nberg(double En n)判断是否是上Hess

11、enberg矩阵void QRAIgorithm(double Ann)QR 算法求特征值void SeekEige nv alue(double Ann)判断是否满足QR算法条件，满足 则进行QR方法求特征值算法的描述（流程图）主的数溢程團:源程序C程序为：#i nclude<stdio.h>#in clude<math.h>#defi ne n 3#defi ne eps 1E-5double Det(double Lnn)用高斯列主元方法求行列式double det=1,t;for(int k=O;k<n-1;k+)double max=Lkk;int Ik

12、=k;for(i nt j=k+1;j <n ;j+)if(max<Ljk)max=Ljk;Ik=j;if(max=0) retur n 0;if(Ik!=k)for(i nt j=k;j< n;j+)t=LIkj;LIkj=Lkj;Lkj=t;det*=-1;for(i nt i=k+1;i <n ;i+)t=Lik/Lkk;Lik=t;for(i nt j=k+1;j <n ;j+)Lij=Lij-t*Lkj;det*=Lkk;if(L n-1 n-1=0) return 0;elseprintf(" 矩阵 A的行列式为：.5fn",det

13、*Ln-1n-1); return det*L n-1 n-1;int Non_si ngularMatrix(double Ln n)/判断是否是非奇异矩阵printf("判断矩阵A的行列式是否为0?n");if(Det(L)!=0) return 1;return 0;void Disp(double Hnn)/输出矩阵for(i nt i=0;i< n;i+)for(i nt j=0;j< n;j+)prin tf("%.6f ",Hij);prin tf("n");int lsZero(double a,int j

14、)/判断数组是否全为 0for(i nt i=0;i<j;i+)if(ai!=0)return 0;int sgn( double y)符号函数if(y>0) return 1;else if(y=0) retur n 0;else return -1;void Hessenberg(double Ann)/将矩阵化为上 Hessenberg 矩阵printf("原矩阵为:n");Disp(A);输出原矩阵doubleT nn ,B nn ,C nn ,c n-1,v n-1,u n-1,R n-1 n-1,l n-1 n-1,t,w,s;int i,j,k,m

15、;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;/ 定义单位阵 Idouble 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=0;i< n-k-1;i+)ci=Ai+k+1k/max;/标准化数组if(lsZero(c,n-k-1)判断数组是否全为0continue;/数组为0,则这一步不需要约化for(i=0,t=0.

16、0;i< n-k-1;i+)t+=ci*ci;vk=sg n(Ak+1k)*sqrt(t);u0=c0+vk;for(j=1;j< n-k-1;j+)uj=cj;/ 求矩阵 叩w=vk*(c0+vk);for(i=0;i< n-k-1;i+)for(j=0;j< n-k-1;j+)Rij=lij-ui*uj/w;for(i=0;i <n ;i+)for(j=0;jv n;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=Ri

17、j;/初等反射阵 Tfor(i=0;i<n;i+) / 矩阵 T左乘矩阵 A for(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+) / 矩阵 T右乘矩阵 Afor(j=0;j< n;j+)for(m=0,s=0.0;m <n ;m+)s+=Bim*Tmj;Cij=s;for(i=0;i <n ;i+)for(j=0;j <n ;j+)Aij=Cij;printf(”原矩阵化成上 Hessenberg 矩阵后为:n");Disp(A);int IsH

18、essenberg(double Enn)判断是否是上 Hessenberg 矩阵for(i nt i=2;i <n ;i+)for(i nt j=0;j+1<i;j+)if(Eij!=0)return 0;int coun t=1;算法求特征值void QRAIgorithm(double An n )/QRint i,j,k,m;double pn n,Q n n,R nn ,F nn ,V nn ,c,s,t,y,max; for(i=0;i<n;i+)赋 Q为单位阵for(i nt j=0;j< n;j+)if(i!=j) Qij=0;else Qij=1;fo

19、r(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);为 Give nsPmm=c;pmm+1=s;pm+1m=-s;pm+1m+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

20、+=pik*Akj;y+=pik*Qkj;Rij=t;Fij=y;for(i=0;i< n;i+)/A,R赋新值for(j=0;j< n;j+)Aij=Rij;Qij=Fij;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+)Qij=Vij;for(i=0;i<n;i+) 矩阵A为矩阵R和矩阵Q的乘积for(j=0;j< n;j+)t=0;for(k=0;k< n;k+)t+=Rik*Qkj;Aij=t;coun t+;max=fab

21、s(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算法迭代次数为:%d次n输出矩阵A%d:n",eps,count-1,count);Disp(A);输出 A矩阵printf("输出矩阵 Q:n");Disp(Q);printf("输出矩阵 R:n");Disp(R);printf("

22、输出A矩阵的全部特征值:");for(i=0;i< n;i+)if(i%3=0) prin tf("n");prin tf("a%d=%.6ft",i+1,Aii);prin tf("n");return;QRAlgorithm(A);判断是否满足QR算法条件，满足则void SeekEige nvalue(double Ann)/ 进行QR方法求特征值double Zn n;for(int i=0;i<n;i+)for(i nt j=0;j< n;j+)Zij=Aij;printf(" 判断矩阵

23、A是否是非奇异矩阵?n");if(Non_si ngularMatrix(Z)=0)printf("矩阵A不是非奇异矩阵！不符合分解条件!n");return;elseprintf("矩阵A是非奇异矩阵！符合分解条件!n");printf(" 判断矩阵A是否是上Hessenberg矩阵?n");if(lsHesse nberg(A)=0)printf(" 不是上Hessenberg矩阵！将其化为上Hessenberg矩阵.n");Hessenberg(A);/将矩阵A转化为上Hessenberg矩阵els

24、eprintf("是上Hessenberg矩阵！不需要将其再化为上 Hessenberg矩阵!n");printf("QR 方法求全部特征值:n");QRAlgorithm(A);/ 用QR方法求特征值void mai n()/*doubleA155=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;SeekEige nvalue(A1);MATLAB?序为：» A=6J2J1;2Jhl；»

25、; B= 2, 3,4,0 5a 3, 7; 0 36,0, 0, 2 3 ； 00, 1,0;» ei£(A)» eig (H)IM1Vsers.fshDesktopDebug.QBj法求特征值.eze测试数据及运行结果测试数据为：6 2们A= 23 1 H1 1 1.0 05620程序结果为:对矩阵Aa3=0.5789333.0000001.000000-2.2360683.4000000.2000000.0000021.0000000.0000001.0000000.0000000.2000000.600000-0.0000000.0000000.57893

26、3-0.0000000.0000000.578933-0.000000 -0.QQQ&QQ1.000000EB *C： WsersVfsh.DesktopDebug.QBj§求特征值.exe*|n|?e +?阵ss 舸矩He ?.J-牙 g - ? 0 ±1 FIJ 0 阵为腼过be为00 矩否腼场en化00 曰苜疋.0铲SS其.0 ll:9?rHe将 1 是至矩是餌00 不昴列异不義00 曰疋门行番疋er00 LUJV - - - b ，&InT n 2 门制门曰疋门se AI7J fi s 阵壮为00 矩断阵A矩上阵00阵曹疋矩00 判 矩判不原6.2.

27、QQQQQQ1.000000原矩阵化成上He00enhe"矩阵后为:6.000000-2.2360680.000000QR忑法求全部特征值： 在精度 1 ra-rarac_f7 e 騎出疽阵A12EJE1：-0.0000042.1330740.000000度恥-0测下,QH算法迭代次数为：11.刖7.287992-0.000004-0.000000输岀矩阵QH1：1.000000 -0.0000020.000000输岀矩阵7.287992-0.0000180-0000002.133074-0.000000 0.000000输岀A 1 矩阵的全部特征值：al=7.287992a2=2.

28、133074Press any key to continue_对矩阵HMWIftfi n 是否是非苛昌i矩阵？囲断楚阵AEHlffiff列我龛否为即矩阵缸的行列式为：50.&&&&&舗库严0是非首异誓险符合分解卷僧M断矩阵A 1 1是否是上Heswenhei'g矩阵？ 是上HewwEnljei'g矩牢不需要将其再化为上HewgenIjeFy矩阵? Q病法求全部特征每査精贋1忌-卿下®算法迭代次数为：22输出矩阵A23：11J 1 . 1113.17234711.2224361.383301-12.287651-2.1843910.0000036.5518