目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。

已知：sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10)算法：以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法：()[]()()[]()[]()111111I 00000i n n n B A I gause i n Q A I u Bu u λλ-⨯-⨯-=-⨯-⎛⎫⎪-=−−−−→=−−−−−−→= ⎪⎝⎭选主元的消元检查知无重特征值由于=0i A I λ-，因此在经过选主元的高斯消元以后，i A I λ-即B 的最后一行必然为零，左上方变为n-1阶单位矩阵[]()()11I n n -⨯-，右上方变为n-1阶向量[]()11n Q ⨯-，然后令n u 1=-，则()1,2,,1j j u Q j n ==⋅⋅⋅-。

这样即求出所有A所有实特征值对应的一个特征向量。

#include<stdio.h>#include<math.h>#include<conio.h>#define N 10#define E 1.0e-12#define MAX 10000//以下是符号函数double sgn(double a){double z;if(a>E) z=1;else z=-1;return z;}//以下是矩阵的拟三角分解void nishangsanjiaodiv(double A[N][N]){int i,j,k;int m=0;double d,c,h,t;double u[N],p[N],q[N],w[N];for(i=0;i<N-2;i++){for(j=i+2;j<N;j++) if(A[j][i]<=E) m=m+1;if(m==(N-2-i)) continue;for(j=i+1,d=0;j<N;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;j<N;j++) u[j]=A[j][i];for(j=0;j<i+2;j++) u[j]=0;u[i+1]=A[i+1][i]-c;for(j=0;j<N;j++){for(k=i+1,p[j]=0;k<N;k++) p[j]=A[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;j<N;j++){for(k=i+1,q[j]=0;k<N;k++) q[j]=A[j][k]*u[k]+q[j];q[j]=q[j]/h;}for(j=0,t=0;j<N;j++) t=t+p[j]*u[j];t=t/h;for(j=0;j<N;j++) w[j]=q[j]-t*u[j];for(j=0;j<N;j++){for(k=0;k<N;k++) A[j][k]=A[j][k]-w[j]*u[k]-u[j]*p[k];}}}//以下是矩阵的QR分解void qrdiv(double A[N][N],double Q[N][N],double R[N][N]){int i,j,k;//int m=0;double d,c,h;double u[N],w[N],p[N];for(i=0;i<N;i++){for(j=0;j<N;j++) {if (i==j) Q[i][j]=1; else Q[i][j]=0;}}for(i=0;i<N;i++){for(j=0;j<N;j++) R[i][j]=A[i][j];}for(i=0;i<N-1;i++){//for(j=i+1;j<N;j++) if(R[j][i]<=E) m=m+1;//if(m==(N-1-i)) continue;for(j=i,d=0;j<N;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<N;j++) u[j]=R[j][i];for(j=0;j<i;j++) u[j]=0;u[i]=R[i][i]-c;for(j=0;j<N;j++) {for(k=0,w[j]=0;k<N;k++) w[j]=Q[j][k]*u[k]+w[j];}for(j=0;j<N;j++) {for(k=0;k<N;k++) Q[j][k]=Q[j][k]-w[j]*u[k]/h;}for(j=0;j<N;j++){for(k=i,p[j]=0;k<N;k++) p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;j<N;j++){for(k=0;k<N;k++) R[j][k]=R[j][k]-u[j]*p[k];}}}//矩阵的QR分解//以下是二次多项式求根double root(double b,double c){double m;m=b*b-4*c;return m;} //二次多项式求根//以下是求解矩阵的所有特征值void characteristic(double A[N][N],double chaR[N],double chaI[N]){int k=0,m=N-1;int i,j;int L;double s,t,x;double M[N][N],B[N][N];int f=0;double d,c,h;double u[N],w[N],p[N];double Q[N][N],R[N][N];for(L=0;L<MAX;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--;goto next;}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;goto next;}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; else Q[i][j]=0;}}for(i=0;i<=m;i++){for(j=0;j<=m;j++) R[i][j]=M[i][j];}for(i=0;i<m;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;j<i;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];} }}}//以下是求矩阵的所有特征值的特征向量void eigenvector(double V[N][N],double T[N]){double A[N][N],baoz[N][N],guod[N];double c;int i,j,k,m,t;int W=0;for(i=0;i<N;i++) for(j=0;j<N;j++) baoz[i][j]=V[i][j];for(t=0;t<6;t++){for(i=0;i<N;i++) for(j=0;j<N;j++) A[i][j]=baoz[i][j];for(i=0;i<N;i++) A[i][i]=A[i][i]-T[t];for(i=0;i<N-1;i++){for(j=i;j<N;j++) if(fabs(A[j][i])>E) {k=j;break; }for(j=i;j<N;j++) {guod[j]=A[i][j];A[i][j]=A[k][j];A[k][j]=guod[j];}for(j=i;j<N;j++){c=A[j][i];if(fabs(c)>E) for(m=i;m<N;m++) A[j][m]=A[j][m]/c;}for(j=0;j<N;j++){c=A[j][i];if(j!=i) {for(m=i;m<N;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];}}}//以下是主函数void main(){double a[N][N],b[N][N],chaR[N],chaI[N];double q[N][N],r[N][N],qr[N][N];double shiyan[N];double f,g;int i,j,k;for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j) a[i][j]=sin(0.5*(i+1)+0.2*(j+1));else a[i][j]=1.5*cos((i+1)+1.2*(j+1));}}nishangsanjiaodiv(a);printf("矩阵A的拟上三角分解：\n");for(i=0;i<N;i++){for(j=0;j<N-5;j++) {if (fabs(a[i][j])<E) a[i][j]=0; printf("%22.11e",a[i][j]);}printf("\n");}printf("\n");for(i=0;i<N;i++){for(j=N-5;j<N;j++) {if (fabs(a[i][j])<E) a[i][j]=0; printf("%22.11e",a[i][j]);} printf("\n");}printf("\n");qrdiv(a,q,r); printf("\n");printf("\n");printf("\n");printf("拟上三角矩阵A的QR分解：\n");printf("上三角矩阵R：\n");for(i=0;i<N;i++){for(j=0;j<N-5;j++) {if (fabs(r[i][j])<E) r[i][j]=0;printf("%22.12e",r[i][j]);}printf("\n");} printf("\n");for(i=0;i<N;i++){for(j=N-5;j<N;j++) {if (fabs(r[i][j])<E) r[i][j]=0;printf("%22.12e",r[i][j]);}printf("\n");} printf("\n");printf("正交矩阵Q：\n");for(i=0;i<N;i++){for(j=0;j<N-5;j++) {if (fabs(q[i][j])<E) q[i][j]=0;printf("%22.12e",q[i][j]);}printf("\n");} printf("\n");for(i=0;i<N;i++){for(j=N-5;j<N;j++) {if (fabs(q[i][j])<E) q[i][j]=0;printf("%22.12e",q[i][j]);}printf("\n");} printf("\n");for(i=0;i<N;i++){for(j=0;j<N;j++) for(k=0,qr[i][j]=0;k<N;k++) qr[i][j]=qr[i][j]+r[i][k]*q[k][j];} printf("\n");printf("\n");printf("\n");printf("R*Q：\n");for(i=0;i<N;i++){for(j=0;j<N-5;j++) {if (fabs(qr[i][j])<E) qr[i][j]=0;printf("%22.12e",qr[i][j]);}printf("\n");} printf("\n");printf("\n");printf("\n");for(i=0;i<N;i++){for(j=N-5;j<N;j++) {if (fabs(qr[i][j])<E) qr[i][j]=0;printf("%22.12e",qr[i][j]);}printf("\n");} printf("\n");printf("\n");printf("\n");characteristic(a,chaR,chaI);for(i=1;i<N;i++){if (i<3) {f=chaR[i];g=chaI[i];chaR[i]=chaR[7+i];chaI[i]=chaI[7+i];chaR[7+i]=f;chaI[7+i]=g;} if (i==5) {f=chaR[i];g=chaI[i];chaR[i]=chaR[7];chaI[i]=chaI[7];chaR[7]=f;chaI[7]=g;}}printf("矩阵A所有特征值：\n");for(j=0;j<N;j++){if(fabs(chaI[j])<=E) printf("λ%2d =%19.11e\n",j+1,chaR[j]);else if(chaI[j]>E) printf("λ%2d =%18.11e +%18.11ei\n",j+1,chaR[j],chaI[j]);else printf("λ%2d =%19.11e %19.11ei\n",j+1,chaR[j],chaI[j]);}printf("\n");printf("\n");printf("\n");for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j) a[i][j]=sin(0.5*(i+1)+0.2*(j+1));else a[i][j]=1.5*cos((i+1)+1.2*(j+1));}}//重新输入矩阵Aeigenvector(a,chaR);printf("相应实特征值对应的特征向量：\n");for(i=0;i<6;i++){printf("λ%d的一个特征向量为：\n ",(i+1));for(j=0;j<N-5;j++) printf("%22.11e",a[i][j]);printf("\n ");for(j=N-5;j<N;j++) printf("%22.11e",a[i][j]);printf("\n");}getch();}。