值分析——带双步位移的QR分解求特征值算法

1、数 值 分 析（B） 大 作 业（二）1、算法设计：矩阵的拟上三角化：对实矩阵A进行相似变换化为拟上三角矩阵，其变换矩阵采用Householder矩阵，变换过程如下：若，则；否则，。当时，得，令又是对称正交矩阵，于是成立，因而与 相似。矩阵的QR分解：矩阵的QR分解过程与拟上三角化过程相似，在这里不再重复其原理。 求全部特征值矩阵拟上三角化后利用带双步位移的QR方法，采用书本Page 63页具体算法实现。为了使编程方便，并体会goto语句使用的灵活性，程序中的跳转均使用goto Loop*实现。求A的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组，。因此，为得到全部实特征

2、值对应的特征向量，解线性方程组的过程要循环n次（n为矩阵阶数）。线性方程组的求解采用列主元素Gauss消去法。#include <stdio.h>#include <math.h>#define ERR1.0e-12/误差限#define N10/矩阵行列数 #define L1.0e5/最大迭代次数double ANN=0;void Init_A()/初始化矩阵int i,j;for(i=0;i<N;i+)for(j=0;j<N;j+)if(i=j)Aij=1.5*cos(i+1)+1.2*(j+1);else Aij=sin(0.5*(i+1)+0.2*

3、(j+1);/*void Display_A()int i,j;for(i=0;i<N;i+)for(j=0;j<N;j+)printf("%8.4f",Aij);printf("n");*/int Sgn(double a)if(a>=0)return 1;else return -1;void On_To_The_Triangle()/矩阵拟上三角化int i,j,r,flag=0;double cr,dr,hr,tr,temp;double urN,prN,qrN,wrN;for(r=1;r<=N-2;r+)flag=0;f

4、or(i=r+2;i<=N;i+)if(Ai-1r-1!=0)flag=1;break;if(0=flag)continue;dr=0;for(i=r+1;i<=N;i+)dr+=Ai-1r-1*Ai-1r-1;dr=sqrt(dr);if(0=Arr-1)cr=dr;else cr=-Sgn(Arr-1)*dr;hr=cr*cr-cr*Arr-1;for(i=1;i<=r;i+)uri-1=0;urr=Arr-1-cr;for(i=r+2;i<=N;i+)uri-1=Ai-1r-1;for(i=1;i<=N;i+)temp=0;for(j=1;j<=N;j

5、+)temp+=Aj-1i-1*urj-1;pri-1=temp/hr;for(i=1;i<=N;i+)temp=0;for(j=1;j<=N;j+)temp+=Ai-1j-1*urj-1;qri-1=temp/hr;temp=0;for(i=1;i<=N;i+)temp+=pri-1*uri-1;tr=temp/hr;for(i=1;i<=N;i+)wri-1=qri-1-tr*uri-1;for(i=1;i<=N;i+)for(j=1;j<=N;j+)Ai-1j-1=Ai-1j-1-wri-1*urj-1-uri-1*prj-1;void Get_Roo

6、ts(double eigenvalue2,int m,double ss,double tt)/求一元二次方程的根double discriminant=ss*ss-4*tt;/if(discriminant<0)*(*(eigenvalue+m-2)=0.5*ss;*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);*(*(eigenvalue+m-1)=0.5*ss;*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);else*(*(eigenvalue+m-2)=0.5*(ss+sqrt(dis

7、criminant);*(*(eigenvalue+m-2)+1)=0;*(*(eigenvalue+m-1)=0.5*(ss-sqrt(discriminant);*(*(eigenvalue+m-1)+1)=0;void Get_Mk(double mkN,int m,double ss,double tt)/获取Mk,用于带双步位移的QR分解int i,j,k;for(i=0;i<m;i+)for(j=0;j<m;j+)*(*(mk+i)+j)=0;for(i=0;i<m;i+)for(j=0;j<m;j+)for(k=0;k<m;k+)*(*(mk+i)+

8、j)+=Aik*Akj;*(*(mk+i)+j)-=ss*Aij;if(j=i)*(*(mk+i)+j)+=tt;void QR_Reslove(double mkN,int m)/QR分解int i,j,r,flag=0;double cr,dr,hr,tr,temp;double urN,vrN,prN,qrN,wrN;double BNN,CNN;for(i=0;i<m;i+)for(j=0;j<m;j+)Bij=*(*(mk+i)+j);Cij=Aij;for(r=1;r<=m-1;r+)flag=0;for(i=r+1;i<=m;i+)if(Bi-1r-1!=

9、0)flag=1;break;if(0=flag)continue;dr=0;for(i=r;i<=m;i+)dr+=Bi-1r-1*Bi-1r-1;dr=sqrt(dr);if(0=Br-1r-1)cr=dr;else cr=-Sgn(Br-1r-1)*dr;hr=cr*cr-cr*Br-1r-1;for(i=1;i<r;i+)uri-1=0;urr-1=Br-1r-1-cr;for(i=r+1;i<=m;i+)uri-1=Bi-1r-1;for(i=1;i<=m;i+)temp=0;for(j=1;j<=m;j+)temp+=Bj-1i-1*urj-1;vri

10、-1=temp/hr;for(i=0;i<m;i+)for(j=0;j<m;j+)Bij-=uri*vrj;for(i=1;i<=m;i+)temp=0;for(j=1;j<=m;j+)temp+=Cj-1i-1*urj-1;pri-1=temp/hr;for(i=1;i<=m;i+)temp=0;for(j=1;j<=m;j+)temp+=Ci-1j-1*urj-1;qri-1=temp/hr;temp=0;for(i=1;i<=m;i+)temp+=pri-1*uri-1;tr=temp/hr;for(i=1;i<=m;i+)wri-1=qr

11、i-1-tr*uri-1;for(i=1;i<=m;i+)for(j=1;j<=m;j+)Ci-1j-1=Ci-1j-1-wri-1*urj-1-uri-1*prj-1;for(i=0;i<m;i+)for(j=0;j<m;j+)Aij=Cij;void Display_Eigenvalue(double value2)/显示特征值int i;for(i=0;i<N;i+)printf("%d=%8.4f",i+1,*(*(value+i);if(*(*(value+i)+1)>0)printf("+%8.4f",*(

12、*(value+i)+1);else if(*(*(value+i)+1)<0)printf("%8.4f",*(*(value+i)+1);printf("n");printf("n");int QR_With_Double_Step_Displacement(double eigenvalue2)/带双步位移QR分解求特征值int i,j,k=1,m=N;double s,t;double MkNN;for(i=0;i<N;i+)for(j=0;j<2;j+)eigenvalueij=0;dok+;if(m=1)

13、eigenvaluem-10=Am-1m-1;m-;continue;else if(m=2)s=Am-2m-2+Am-1m-1;t=Am-2m-2*Am-1m-1-Am-1m-2*Am-1m-2;Get_Roots(eigenvalue,m,s,t);/求一元二次方程的根m=0;continue;else if(m=0)return 0;else if(fabs(Am-1m-2)<=ERR)eigenvaluem-10=Am-1m-1;m-;continue;elses=Am-2m-2+Am-1m-1;t=Am-2m-2*Am-1m-1-Am-1m-2*Am-2m-1;if(fabs(

14、Am-2m-3)<=ERR)Get_Roots(eigenvalue,m,s,t);/求一元二次方程的根m-=2;continue;elseGet_Mk(Mk,m,s,t);/获取Mk,用于带双步位移的QR分解QR_Reslove(Mk,m);/QR分解while(k<L);printf("Errern");/迭代过程不成功，报错void Select_Principal_Element(int k)/Gauss消元法求解方程组之选主元int i,max=k;double transN;for(i=k+1;i<N;i+)if(fabs(Aik)>fa

15、bs(Amaxk)max=i;if(k=max)return;elsefor(i=k;i<N;i+)transi=Aki;Aki=Amaxi;Amaxi=transi;void Eliminant(int k)/Gauss消元法求解方程组之消元int i,j;double rate;for(i=k+1;i<N;i+)rate=Aik/Akk;for(j=k;j<N;j+)Aij=Aij-rate*Akj;void Back_Substitution(double Eigenvector)/Gauss消元法求解方程组之回代int i,j,k;double Temp_MNN+1;

16、for(i=0;i<N-1;i+)for(j=i;j<N;j+)Temp_Mij=Aij;Temp_MiN=0;Temp_MN-1N-1=1;Temp_MN-1N=1;for(k=N-1;k>=0;k-)*(Eigenvector+k)=Temp_MkN;for(i=k;i<N-1;i+)*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_Mki+1;*(Eigenvector+k)=*(Eigenvector+k)/Temp_Mkk;void Display_Eigenvector(double Eige

17、nvector,double eigenvalue)/显示特征值对应特征向量int i;printf("When =%8.4fnEigenvector=n",eigenvalue);for(i=0;i<N;i+)printf("%8.4f",*(Eigenvector+i);printf("n");void Get_Eigenvector(double value2)/利用Gauss消元法求解特征向量int i,j;double eigenvalue;double EigenvectorN;for(i=0;i<N;i+)i

18、f(*(*(value+i)+1)=0)Init_A();/初始化矩阵eigenvalue=*(*(value+i);for(j=0;j<N;j+)Ajj-=eigenvalue;for(j=0;j<N-1;j+)Select_Principal_Element(j);/Gauss消元法求解方程组之选主元Eliminant(j);/Gauss消元法求解方程组之消元Back_Substitution(Eigenvector);/Gauss消元法求解方程组之回代Display_Eigenvector(Eigenvector,eigenvalue);/显示特征值对应特征向量main()double eigenvalueN2;Init_A(