北航数值分析第二次大作业--QR分解

0《数值分析A》计算实习题目二姓名学号联系方式班级指导教师2012年10月1一、算法设计方案整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。2.对经过拟上三角化处理的矩阵进行QR分解。3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。2二、源程序代码#includestdio.h#includemath.h#includestring.hinti,j,k,l,m;//定义外部变量doubled,h,b,c,t,s;doubleA[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10];doubleX[10][10],Y[10][10],Qt[10][10],M[10][10];doubleU[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0};doubleepsilon=1e-12;voidmain(){voidQuasiuppertriangular(doubleA[][10]);voidQRdecomposition(doubleA[][10]);voidDoublestepsQR(doubleA[][10]);inti,j;for(i=0;i10;i++){for(j=0;j10;j++){A[i][j]=sin(0.5*(i+1)+0.2*(j+1));Q[i][j]=0;AA[i][j]=A[i][j];}A[i][i]=1.5*cos(2.2*(i+1));AA[i][i]=A[i][i];3}Quasiuppertriangular(A);//调用拟上三角化函数printf(\nA经过拟上三角化矩阵为：\n\n);for(i=0;i10;i++)//输出拟上三角化矩阵{for(j=0;j10;j++){printf(%.12e,A[i][j]);//输出拟上三角化矩阵}printf(\n\n);}QRdecomposition(A);//调用QR分解函数printf(进行QR分解后，R矩阵为：\n\n);//输出R矩阵for(i=0;i10;i++){for(j=0;j10;j++){printf(%.12e,R[i][j]);}printf(\n\n);}printf(Q矩阵为:\n\n);//输出Q矩阵for(i=0;i10;i++){for(j=0;j10;j++){printf(%.12e,Q[i][j]);4}printf(\n\n);}printf(RQ矩阵为：\n\n);//输出RQ矩阵for(i=0;i10;i++){for(j=0;j10;j++){printf(%.12e,RQ[i][j]);}printf(\n\n);}DoublestepsQR(A);//调用双步位移函数printf(\n\n特征值实部依次为：\n\n);//输出特征值实部for(j=0;j10;j++){printf(%.12e,Re[j]);}printf(\n\n特征值虚部依次为：\n\n);//输出特征值虚部for(j=0;j10;j++){printf(%.12e,Im[j]);}//按行输出特征向量printf(\n\n按行输出实特征根相应特征向量为：\n\n);for(i=0;i10;i++){5if(i==1||i==2||i==5||i==6){continue;}for(j=0;j10;j++){printf(%.12e,X[i][j]);}printf(\n\n);}getchar();}//拟上三角化函数voidQuasiuppertriangular(doubleA[][10]){for(j=0;j8;j++){for(i=0;i10;i++){U[i]=0;P[i]=0;T[i]=0;W[i]=0;}m=0;for(i=j+2;i10;i++){if(A[i][j]!=0)6{m=m+1;}}if(m==0){continue;}d=0;for(i=j+1;i10;i++){d=d+pow(A[i][j],2);}d=sqrt(d);c=-d;if(A[j+1][j]=0){c=d;}h=c*(c-A[j+1][j]);U[j+1]=A[j+1][j]-c;for(i=j+2;i10;i++){U[i]=A[i][j];}for(i=0;i10;i++){7for(k=0;k10;k++){P[i]=P[i]+U[k]*A[k][i];}P[i]=P[i]/h;}t=0;for(i=0;i10;i++){for(k=0;k10;k++){T[i]=T[i]+U[k]*A[i][k];}T[i]=T[i]/h;t=t+P[i]*U[i];}t=t/h;for(i=0;i10;i++){W[i]=T[i]-t*U[i];for(k=0;k10;k++){A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];if(abs(A[i][k])1e-12){A[i][k]=0;}8}}}}//QR分解函数voidQRdecomposition(doubleA[][10]){for(i=0;i10;i++){for(j=0;j10;j++){RQ[i][j]=0;Q[i][j]=0;R[i][j]=A[i][j];}Q[i][i]=1;}for(j=0;j9;j++){for(i=0;i10;i++){U[i]=0;P[i]=0;W[i]=0;}m=0;for(i=j+1;i10;i++){if(R[i][j]!=0)9{m=m+1;}}if(m==0){continue;}d=0;for(i=j;i10;i++){d=d+pow(R[i][j],2);}d=sqrt(d);c=-d;if(R[j][j]=0){c=d;}h=c*(c-R[j][j]);U[j]=R[j][j]-c;for(i=j+1;i10;i++){U[i]=R[i][j];}10for(i=0;i10;i++){for(k=0;k10;k++){W[i]=W[i]+U[k]*Q[i][k];}}for(i=0;i10;i++){for(k=0;k10;k++){Q[i][k]=Q[i][k]-((W[i]*U[k])/h);}}for(i=0;i10;i++){for(k=0;k10;k++){P[i]=P[i]+U[k]*R[k][i];}P[i]=P[i]/h;}for(i=0;i10;i++){for(k=0;k10;k++){R[i][k]=R[i][k]-U[i]*P[k];if(abs(R[i][k])epsilon){11R[i][k]=0;}}}}for(i=0;i10;i++)//计算A(n+1)=RQ{for(j=0;j10;j++){for(k=0;k10;k++){RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];}}}}//双步位移法计算特征值特征向量函数voidDoublestepsQR(doubleA[][10]){intL=1000,m=9;//定义最大循环次数for(i=0;iL;i++){for(;m-1;){if(abs(A[m][m-1])=epsilon){Re[m]=A[m][m];m=m-1;//降阶if(m==0)//4{12Re[0]=A[0][0];break;}if(m==-1){break;}if(m1){continue;}}b=-A[m][m]-A[m-1][m-1];//5c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];if(m==1)//6{if((b*b-4*c)=0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)0){Re[m]=-b/2;Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2;Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-1;//循环出口条件break;}if((m1)&&(abs(A[m-1][m-2])epsilon))//813{if(i==L-1){printf(Noresults!\n);m=0;//循环出口条件break;}break;}if((m1)&&(abs(A[m-1][m-2])=epsilon))//7{if((b*b-4*c)0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)0){Re[m]=-b/2;Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2;Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-2;//降阶if(m0){continue;}if(m==0){Re[0]=A[0][0];break;14}}}if(m=0){break;}s=A[m-1][m-1]+A[m][m];//9t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];for(j=0;j10;j++){for(k=0;k10;k++){Qt[j][k]=0;Q[j][k]=0;M[j][k]=0;X[j][k]=0;Y[j][k]=0;}}for(j=0;jm+1;j++){for(k=0;km+1;k++){for(l=0;lm+1;l++){M[j][k]=M[j][k]+A[j][l]*A[l][k];}15}}for(j=0;jm+1;j++){for(k=0;km+1;k++){M[j][k]=M[j][k]-s*A[j][k];}M[j][j]=M[j][j]+t;}//调用QR分解函数对M矩阵进行分解并传递参数矩阵QQRdecomposition(M);for(j=0;j10;j++){for(k=0;k10;k++){Qt[j][k]=Q[k][j];}}for(j=0;jm+1;j++){for(k=0;km+1;k++){for(l=0;lm+1;l++){X[j][k]=X[j][k]+Qt[j][l]*A[l][k];}}}16for(j=0;jm+1;j++){for(k=0;km+1;k++){for(l=0;lm+1;l++){Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];}}}for(j=0;j10;j++){for(k=0;k10;k++){A[j][k]=Y[j][k];}}}//应用列主元高斯消元法计算实部特征向量for(l=0;l10;l++){if(l==1||l==2||l==5||l==6){continue;}for(k=0;k10;k++){for(m=0;m10;m++){A[k][m]=AA[k][m];17}A[k][k]=A[k][k]-Re[l];}for(j=0;j9;j++){m=j;for(i=j+1;i10;i++){if(abs(A[i][j])abs(A[m