北航数值分析大作业第二题精解

目标：使用带双步位移的QR分解法求矩阵10*10[]ijAa的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)()1.5cos(1.2)(){ijijijijija(i,j=1,2,……,10)算法：以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A的实特征值对应的一个特征向量的求法：111111I00000innnBAIgauseinQAIuBuu选主元的消元检查知无重特征值由于=0iAI，因此在经过选主元的高斯消元以后，iAI即B的最后一行必然为零，左上方变为n-1阶单位矩阵11Inn，右上方变为n-1阶向量11nQ，然后令nu1，则1,2,,1jjuQjn。开始输入矩阵A用函数“nishangsanjiadiv”将矩阵A拟上三角化即为A（n-1）用函数“characteristic”求解矩阵A（n-1）即A的所有特征值用函数“qrdiv”将A（n-1）QR分解用函数“characteristicvector”求解矩阵A（n-1）即A的所有实特征值对应的特征向量输出A的所有特征值λ、A的所有实特征值对应的特征向量、拟上三角矩阵A（n-1）、及其Q、R和R*Q这样即求出所有A所有实特征值对应的一个特征向量。#includestdio.h#includemath.h#includeconio.h#defineN10#defineE1.0e-12#defineMAX10000//以下是符号函数doublesgn(doublea){doublez;if(aE)z=1;elsez=-1;returnz;}//以下是矩阵的拟三角分解voidnishangsanjiaodiv(doubleA[N][N]){inti,j,k;intm=0;doubled,c,h,t;doubleu[N],p[N],q[N],w[N];for(i=0;iN-2;i++){for(j=i+2;jN;j++)if(A[j][i]=E)m=m+1;if(m==(N-2-i))continue;for(j=i+1,d=0;jN;j++)d=d+A[j][i]*A[j][i];d=sqrt(d);c=-1*sgn(A[i+1][i])*d;h=c*c-c*A[i+1][i];for(j=i+2;jN;j++)u[j]=A[j][i];for(j=0;ji+2;j++)u[j]=0;u[i+1]=A[i+1][i]-c;for(j=0;jN;j++){for(k=i+1,p[j]=0;kN;k++)p[j]=A[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;jN;j++){for(k=i+1,q[j]=0;kN;k++)q[j]=A[j][k]*u[k]+q[j];q[j]=q[j]/h;}for(j=0,t=0;jN;j++)t=t+p[j]*u[j];t=t/h;for(j=0;jN;j++)w[j]=q[j]-t*u[j];for(j=0;jN;j++){for(k=0;kN;k++)A[j][k]=A[j][k]-w[j]*u[k]-u[j]*p[k];}}}//以下是矩阵的QR分解voidqrdiv(doubleA[N][N],doubleQ[N][N],doubleR[N][N]){inti,j,k;//intm=0;doubled,c,h;doubleu[N],w[N],p[N];for(i=0;iN;i++){for(j=0;jN;j++){if(i==j)Q[i][j]=1;elseQ[i][j]=0;}}for(i=0;iN;i++){for(j=0;jN;j++)R[i][j]=A[i][j];}for(i=0;iN-1;i++){//for(j=i+1;jN;j++)if(R[j][i]=E)m=m+1;//if(m==(N-1-i))continue;for(j=i,d=0;jN;j++)d=d+R[j][i]*R[j][i];d=sqrt(d);c=-1*sgn(R[i][i])*d;h=c*c-c*R[i][i];for(j=i+1;jN;j++)u[j]=R[j][i];for(j=0;ji;j++)u[j]=0;u[i]=R[i][i]-c;for(j=0;jN;j++){for(k=0,w[j]=0;kN;k++)w[j]=Q[j][k]*u[k]+w[j];}for(j=0;jN;j++){for(k=0;kN;k++)Q[j][k]=Q[j][k]-w[j]*u[k]/h;}for(j=0;jN;j++){for(k=i,p[j]=0;kN;k++)p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;jN;j++){for(k=0;kN;k++)R[j][k]=R[j][k]-u[j]*p[k];}}}//矩阵的QR分解//以下是二次多项式求根doubleroot(doubleb,doublec){doublem;m=b*b-4*c;returnm;}//二次多项式求根//以下是求解矩阵的所有特征值voidcharacteristic(doubleA[N][N],doublechaR[N],doublechaI[N]){intk=0,m=N-1;inti,j;intL;doubles,t,x;doubleM[N][N],B[N][N];intf=0;doubled,c,h;doubleu[N],w[N],p[N];doubleQ[N][N],R[N][N];for(L=0;LMAX;L++){next:if(m==0){chaR[0]=A[0][0];chaI[0]=0;break;}if(fabs(A[m][m-1])=E){chaR[m]=A[m][m];chaI[m]=0;m--;gotonext;}s=A[m-1][m-1]+A[m][m];t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];if(m==1){x=root(s,t);if(x=E){x=sqrt(x);chaR[m]=s/2+x/2;chaR[m-1]=s/2-x/2;chaI[m]=0;chaI[m-1]=0;}else{x=sqrt(fabs(x));chaR[m]=s/2;chaR[m-1]=s/2;chaI[m]=x/2;chaI[m-1]=-x/2;}break;}if(fabs(A[m-1][m-2])=E){x=root(s,t);if(x=E){x=sqrt(x);chaR[m]=s/2+x/2;chaR[m-1]=s/2-x/2;chaI[m]=0;chaI[m-1]=0;}else{x=sqrt(fabs(x));chaR[m]=s/2;chaR[m-1]=s/2;chaI[m]=x/2;chaI[m-1]=-x/2;}m=m-2;gotonext;}for(i=0;i=m;i++){for(j=0;j=m;j++){if(i==j){for(k=0,M[i][j]=0;k=m;k++)M[i][j]=A[i][k]*A[k][j]+M[i][j];M[i][j]=M[i][j]-s*A[i][j]+t;}else{for(k=0,M[i][j]=0;k=m;k++)M[i][j]=A[i][k]*A[k][j]+M[i][j];M[i][j]=M[i][j]-s*A[i][j];}}}//以下是M的QR分解for(i=0;i=m;i++){for(j=0;j=m;j++){if(i==j)Q[i][j]=1;elseQ[i][j]=0;}}for(i=0;i=m;i++){for(j=0;j=m;j++)R[i][j]=M[i][j];}for(i=0;im;i++){for(j=i+1;j=m;j++)if(R[j][i]=E)f=f+1;if(f==(m-i))continue;for(j=i,d=0;j=m;j++)d=d+R[j][i]*R[j][i];d=sqrt(d);c=-1*sgn(R[i][i])*d;h=c*c-c*R[i][i];for(j=i+1;j=m;j++)u[j]=R[j][i];for(j=0;ji;j++)u[j]=0;u[i]=R[i][i]-c;for(j=0;j=m;j++){for(k=0,w[j]=0;k=m;k++)w[j]=Q[j][k]*u[k]+w[j];}for(j=0;j=m;j++){for(k=0;k=m;k++)Q[j][k]=Q[j][k]-w[j]*u[k]/h;}for(j=0;j=m;j++){for(k=i,p[j]=0;k=m;k++)p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;j=m;j++){for(k=0;k=m;k++)R[j][k]=R[j][k]-u[j]*p[k];}}for(j=0;j=m;j++){for(k=0;k=m;k++)M[j][k]=Q[j][k];}//以上是M的QR分解for(i=0;i=m;i++){for(j=0;j=m;j++){for(k=0,B[i][j]=0;k=m;k++)B[i][j]=M[k][i]*A[k][j]+B[i][j];}}for(i=0;i=m;i++){for(j=0;j=m;j++){for(k=0,A[i][j]=0;k=m;k++)A[i][j]=B[i][k]*M[k][j]+A[i][j];}}}}//以下是求矩阵的所有特征值的特征向量voideigenvector(doubleV[N][N],doubleT[N]){doubleA[N][N],baoz[N][N],guod[N];doublec;inti,j,k,m,t;intW=0;for(i=0;iN;i++)for(j=0;jN;j++)baoz[i][j]=V[i][j];for(t=0;t6;t++){for(i=0;iN;i++)for(j=0;jN;j++)A[i][j]=baoz[i][j];for(i=0;iN;i++)A[i][i]=A[i][i]-T[t];for(i=0;iN-1;i++){for(j=i;jN;j++)if(fabs(A[j][i])E){k=j;break;}for(j=i;jN;j++){guod[j]=A[i][j];A[i][j]=A[k][j];A[k][j]=guod[j];}for(j=i;jN;j++){c=A[j][i];if(fabs(c)E)for(m=i;mN;m++)A[j][m]=A[j][m]/c;}for(j=0;jN;j++){c=A[j][i];if(j!=i){for(m=i;mN;m++)A[j][m]=A[j][m]-A[i][m]*c;}}}V[t][N-1]=-1;for(i=N-2;i=0;i--){V[t][i]=A[i][N-1];}}}//以下是主函数voidmain(){doublea[N][N],b[N][N],chaR[N],chaI[N];doubleq[N][N],r[N][N],qr[N][N];doubleshiyan[N];doublef,g;inti,j,k;for(i=0;iN;i++){for(j=0;jN;j++){if(i!=j)a[i][j]=sin(0.5*(i+1)+0.2*(j+1));elsea[i][j]=1.5*