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

数值分析（B）大作业（二）1、算法设计： 矩阵的拟上三角化：对实矩阵A进行相似变换化为拟上三角矩阵 A（ n-l）,其变换矩阵采用 Householder矩阵，变换过程如下：若 a（r）= 0（i = r + 2,n），则 H = I ；irr否则，S =（0，0，a（r） ，a（r）T,rr+1,rn,rc = -sgn（a（r） ）lls II ，rr+1,r r 2u =s -c e = （0,0, a（r） c ,a（r）a（r），r r r r+1r+l,rr r+2, rnrH =I-2uut/|u |2，A（r+1） = H A（r）H。rr当 r 二 n 2 时，得 A（n1）= H A（n2）H= .二 H H AH H ，n2n2n211n2令P=HHH 又H是对称正交矩阵，于是An-1=PTAP成立，因而An-1与 a相12n-2r似。 矩阵的 QR 分解：矩阵的 QR 分解过程与拟上三角化过程相似，在这里不再重复其原理。 求全部特征值矩阵拟上三角化后利用带双步位移的 QR 方法，采用书本 Page 63 页具体算 法实现。为了使编程方便，并体会 goto 语句使用的灵活性，程序中的跳转均使 用 goto Loop*实现。 求 A 的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组（41-人”=，LL（i =n）。因此，为得到全部实特征值对应的特征向量，解线性方程组的过 程要循环n次（n为矩阵阶数）。线性方程组的求解采用列主元素Gauss消去法。#include #include #define ERR 1.0e-12 /误差限#define N10 /矩阵行列数#define L 1.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); elseAij=sin(0.5*(i+1)+0.2*(j+1); /*void Display_A() int i,j;for(i=0;iN;i+)for(j=0;j=0) return 1;elsereturn -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;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; 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);*(*(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;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;/QR 分解void QR_Reslove(double mkN,int m)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+)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;for(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/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);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); /求一元二次方程的根 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 分解while(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(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 Display_Eigenvector(double Eigenvector,double eigenvalue) / 显示特征值对应特征向 量int i;pri ntf( When 2=%8nEfigenvector=n,eigenvalue);for(i=0;iN;i+)printf(%8.4f,*(Eigenvector+i);printf(n);void Get_Eigenvector(double value2)int i,j;double eigenvalue;double EigenvectorN;利用Gauss消元法求解特征向量/初始化矩阵/Gauss 消元法求解方程组之选主元/Gauss 消元法求解方程组之消元/Gauss 消元法求解方程组之回代for(i=0;iN;i+) if(*(*(value+i)+1)=0) Init_A();eigenvalue=*(*(value+i);for(j=0;jN;j+)Ajj-=eigenvalue;for(j=0;jN-1;j+) Select_Principal_Element(j);Eliminant(j);Back_Substitution(Eigenvector);Display_Eigenvector(Eigenvector,eigenvalue); /显示特征值对应特征向量main()double eigenvalueN2;Init_A();/初始化矩阵On_To_The_Triangle();/矩阵上三角化QR_With_Double_Step_Displacement(eigenvalue); /带双步位移 QR 分解求特征值 printf(Contactme:*nn);Display_Eigenvalue(eigenvalue);/显示特征值Get_Eigenvector(eigenvalue);/利用 Gauss 消元法求解特征向量return 0;bdDOCUIENTS1 .00011.酗醐1 .00015I5J5I5Contact me 二 731363227(?qq.conx.1-3.389ESX.2= -2.323E+0.B93OX?= -2.223E； -R.H93B kg 1.577515= -1.4845-9.9805+0.113917= -0.9805 -0.1139 入 fjw 9.9356 kQ= 9.6361When X -3.3830When X =1,F7?SEis(enuector=0.6217 -0.1115-2.483S-1.3R69 -3.81568.1173-1.2392-9.68E02.6919iJhen A = -1.840Eigenuector=1?.3078 24；8190,4134-8.5720B.H929-0,0?83When X =6.53561.0Milen X=- 0.B565Eigenuector= -S.ltlSH -4,88639.5W52-0.b788 -9.6B43-3.U457 15.748? -7.395B -7.1897continue-H.3403 -0.378SEigenvector-0.43S7 -0.9064-0.B374-1.6&79-12.25717.2412-5.398228.4101-12丄652When X=- R,fi361 teigenuector=4.74503.1579-1.9800-31.87527.7940-10.042616.7076 13.1052-1.0829 -1.2fc93一丄.07240.36251.68292丄12?SETTINGSYADII町STRAIO班桌ffihTAJ?邛 OF NTJIERICAL AIAJY. .PS: Forfurtherdetails,pleasecontactme: