值分析——帶雙步位移的QR分解求特征值算法

1、數(shù) 值 分 析（B） 大 作 業(yè)（二）1、算法設(shè)計：矩陣的擬上三角化：對實矩陣A進行相似變換化為擬上三角矩陣，其變換矩陣采用Householder矩陣，變換過程如下：若，則；否則，。當時，得，令又是對稱正交矩陣，于是成立，因而與 相似。矩陣的QR分解：矩陣的QR分解過程與擬上三角化過程相似，在這里不再重復(fù)其原理。 求全部特征值矩陣擬上三角化后利用帶雙步位移的QR方法，采用書本Page 63頁具體算法實現(xiàn)。為了使編程方便，并體會goto語句使用的靈活性，程序中的跳轉(zhuǎn)均使用goto Loop*實現(xiàn)。求A的相應(yīng)于實特征值的特征向量求實特征值對應(yīng)的特征向量，即是求解線性方程組，。因此，為得到全部實特征

2、值對應(yīng)的特征向量，解線性方程組的過程要循環(huán)n次（n為矩陣階數(shù)）。線性方程組的求解采用列主元素Gauss消去法。#include <stdio.h>#include <math.h>#define ERR1.0e-12/誤差限#define N10/矩陣行列數(shù) #define L1.0e5/最大迭代次數(shù)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)/顯示特征值對應(yīng)特征向量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);/顯示特征值對應(yīng)特征向量main()double eigenvalueN2;Init_A(