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

数值分析（B）大作业（二）1、算法设计：①矩阵的拟上三角化：对实矩阵A进行相似变换化为拟上三角矩阵(1)An，其变换矩阵采用Householder矩阵，变换过程如下：若()a0(2,,)ririrn，则rHI；否则，(r)(r)Trr+1,rn,rs=(0,,0,a,a)，(r)rr+1,rr2c=-sgn(a)||s||，()()()rrrr+11,2,u=s-ce(0,,0,,,,)rrrrrrrrnracaa，2Trrr2H=I-2uu/ru，(1)()rrrrAHAH。当2rn时，得(1)(2)(1)222112nnnnnnAHAHHHAHH，令12n-2P=HHH又rH是对称正交矩阵，于是n-1TA=PAP成立，因而n-1A与A相似。②矩阵的QR分解：矩阵的QR分解过程与拟上三角化过程相似，在这里不再重复其原理。③求全部特征值矩阵拟上三角化后利用带双步位移的QR方法，采用书本Page63页具体算法实现。为了使编程方便，并体会goto语句使用的灵活性，程序中的跳转均使用gotoLoop*实现。④求A的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组(λI-A)=0iix，(1,,)in。因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n次（n为矩阵阶数）。线性方程组的求解采用列主元素Gauss消去法。#includestdio.h#includemath.h#defineERR1.0e-12//误差限#defineN10//矩阵行列数#defineL1.0e5//最大迭代次数doubleA[N][N]={0};voidInit_A()//初始化矩阵{inti,j;for(i=0;iN;i++)for(j=0;jN;j++){if(i==j)A[i][j]=1.5*cos((i+1)+1.2*(j+1));elseA[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}/*voidDisplay_A(){inti,j;for(i=0;iN;i++){for(j=0;jN;j++)printf(%8.4f,A[i][j]);printf(\n);}}*/intSgn(doublea){if(a=0)return1;elsereturn-1;}voidOn_To_The_Triangle()//矩阵拟上三角化{inti,j,r,flag=0;doublecr,dr,hr,tr,temp;doubleur[N],pr[N],qr[N],wr[N];for(r=1;r=N-2;r++){flag=0;for(i=r+2;i=N;i++)if(A[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r+1;i=N;i++)dr+=A[i-1][r-1]*A[i-1][r-1];dr=sqrt(dr);if(0==A[r][r-1])cr=dr;elsecr=-Sgn(A[r][r-1])*dr;hr=cr*cr-cr*A[r][r-1];for(i=1;i=r;i++)ur[i-1]=0;ur[r]=A[r][r-1]-cr;for(i=r+2;i=N;i++)ur[i-1]=A[i-1][r-1];for(i=1;i=N;i++){temp=0;for(j=1;j=N;j++)temp+=A[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i=N;i++){temp=0;for(j=1;j=N;j++)temp+=A[i-1][j-1]*ur[j-1];qr[i-1]=temp/hr;}temp=0;for(i=1;i=N;i++){temp+=pr[i-1]*ur[i-1];tr=temp/hr;}for(i=1;i=N;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i=N;i++)for(j=1;j=N;j++)A[i-1][j-1]=A[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}}voidGet_Roots(doubleeigenvalue[][2],intm,doubless,doublett)//求一元二次方程的根{doublediscriminant=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;}}voidGet_Mk(doublemk[][N],intm,doubless,doublett)//获取Mk,用于带双步位移的QR分解{inti,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)+=A[i][k]*A[k][j];*(*(mk+i)+j)-=ss*A[i][j];if(j==i)*(*(mk+i)+j)+=tt;}}voidQR_Reslove(doublemk[][N],intm)//QR分解{inti,j,r,flag=0;doublecr,dr,hr,tr,temp;doubleur[N],vr[N],pr[N],qr[N],wr[N];doubleB[N][N],C[N][N];for(i=0;im;i++)for(j=0;jm;j++){B[i][j]=*(*(mk+i)+j);C[i][j]=A[i][j];}for(r=1;r=m-1;r++){flag=0;for(i=r+1;i=m;i++)if(B[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r;i=m;i++)dr+=B[i-1][r-1]*B[i-1][r-1];dr=sqrt(dr);if(0==B[r-1][r-1])cr=dr;elsecr=-Sgn(B[r-1][r-1])*dr;hr=cr*cr-cr*B[r-1][r-1];for(i=1;ir;i++)ur[i-1]=0;ur[r-1]=B[r-1][r-1]-cr;for(i=r+1;i=m;i++)ur[i-1]=B[i-1][r-1];for(i=1;i=m;i++){temp=0;for(j=1;j=m;j++)temp+=B[j-1][i-1]*ur[j-1];vr[i-1]=temp/hr;}for(i=0;im;i++)for(j=0;jm;j++)B[i][j]-=ur[i]*vr[j];for(i=1;i=m;i++){temp=0;for(j=1;j=m;j++)temp+=C[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i=m;i++){temp=0;for(j=1;j=m;j++)temp+=C[i-1][j-1]*ur[j-1];qr[i-1]=temp/hr;}temp=0;for(i=1;i=m;i++){temp+=pr[i-1]*ur[i-1];tr=temp/hr;}for(i=1;i=m;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i=m;i++)for(j=1;j=m;j++)C[i-1][j-1]=C[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}for(i=0;im;i++)for(j=0;jm;j++)A[i][j]=C[i][j];}voidDisplay_Eigenvalue(doublevalue[][2])//显示特征值{inti;for(i=0;iN;i++){printf(λ%d=%8.4f,i+1,*(*(value+i)));if(*(*(value+i)+1)0)printf(+%8.4f,*(*(value+i)+1));elseif(*(*(value+i)+1)0)printf(%8.4f,*(*(value+i)+1));printf(\n);}printf(\n);}intQR_With_Double_Step_Displacement(doubleeigenvalue[][2])//带双步位移QR分解求特征值{inti,j,k=1,m=N;doubles,t;doubleMk[N][N];for(i=0;iN;i++)for(j=0;j2;j++)eigenvalue[i][j]=0;do{k++;if(m==1){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}elseif(m==2){s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-1][m-2];Get_Roots(eigenvalue,m,s,t);//求一元二次方程的根m=0;continue;}elseif(m==0)return0;elseif(fabs(A[m-1][m-2])=ERR){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}else{s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];if(fabs(A[m-2][m-3])=ERR){Get_Roots(eigenvalue,m,s,t);//求一元二次方程的根m-=2;continue;}else{Get_Mk(Mk,m,s,t);//获取Mk,用于带双步位移的QR分解QR_Reslove(Mk,m);//QR分解}}}while(kL);printf(Errer\n);//迭代过程不成功，报错}voidSelect_Principal_Element(intk)//Gauss消元法求解方程组之选主元{inti,max=k;doubletrans[N];for(i=k+1;iN;i++)if(fabs(A[i][k])fabs(A[max][k]))max=i;if(k==max)return;else{for(i=k;iN;i++){trans[i]=A[k][i];A[k][i]=A[max][i];A[max][i]=trans[i];}}}voidEliminant(intk)//Gauss消元法求解方程组之消元{inti,j;doublerate;for(i=k+1;iN;i++){rate=A[i][k]/A[k][k];for(j=k;jN;j++)A[i][j]=A[i][j]-rate*A[k][j];}}voidBack_Substitution(doubleEigenvector[])//Gauss消元法求解方程组之回代{inti,j,k;doubleTemp_M[N][N+1];for(i=0;