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

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

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

3、=0;iN;i+)for(j=0;j=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;for(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;

4、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+)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

5、;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_Roots(double eigenvalue2,int m,double ss,double tt)/求一元二次方程的根double discriminant=ss*ss-4*tt;/if(discriminant0)*(*(eigenvalue+m-2)=0.5*ss;*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);

6、*(*(eigenvalue+m-1)=0.5*ss;*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);else*(*(eigenvalue+m-2)=0.5*(ss+sqrt(discriminant);*(*(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;

7、for(i=0;im;i+)for(j=0;jm;j+)*(*(mk+i)+j)=0;for(i=0;im;i+)for(j=0;jm;j+)for(k=0;km;k+)*(*(mk+i)+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;im;i+)for(j=0;jm;j+

8、)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!=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;ir;i+)uri-1=0;urr-1=Br-1r-1-cr;for(i=r+1;i=m;i+)uri-1=Bi-1r-1;f

9、or(i=1;i=m;i+)temp=0;for(j=1;j=m;j+)temp+=Bj-1i-1*urj-1;vri-1=temp/hr;for(i=0;im;i+)for(j=0;jm;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

10、/hr;for(i=1;i=m;i+)wri-1=qri-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;im;i+)for(j=0;jm;j+)Aij=Cij;void Display_Eigenvalue(double value2)/显示特征值int i;for(i=0;i0)printf(+%8.4f,*(*(value+i)+1);else if(*(*(value+i)+1)0)printf(%8.4f,*(*(value+i)+1);printf(n)

11、;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;iN;i+)for(j=0;j2;j+)eigenvalueij=0;dok+;if(m=1)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

12、);/求一元二次方程的根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(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分解whi

13、le(kL);printf(Errern);/迭代过程不成功，报错void Select_Principal_Element(int k)/Gauss消元法求解方程组之选主元int i,max=k;double transN;for(i=k+1;ifabs(Amaxk)max=i;if(k=max)return;elsefor(i=k;iN;i+)transi=Aki;Aki=Amaxi;Amaxi=transi;void Eliminant(int k)/Gauss消元法求解方程组之消元int i,j;double rate;for(i=k+1;iN;i+)rate=Aik/Akk;for(

14、j=k;jN;j+)Aij=Aij-rate*Akj;void Back_Substitution(double Eigenvector)/Gauss消元法求解方程组之回代int i,j,k;double Temp_MNN+1;for(i=0;iN-1;i+)for(j=i;j=0;k-)*(Eigenvector+k)=Temp_MkN;for(i=k;iN-1;i+)*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_Mki+1;*(Eigenvector+k)=*(Eigenvector+k)/Temp_Mkk;void D

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

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