北航数值分析大作业第二次的

数值分析第二次大作业姓名：路子威学号：S201451401题目：使用带双步位移的QR分解法求矩阵10*10[]ijAa的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)()1.5cos(1.2)(){ijijijijija(i,j=1,2,……,10)要求：1.用幂法、反幂法和QR方法求矩阵的特征值时，要求迭代精度水平为e=10-12。2.打印以下内容：(1)全部源程序；(2)矩阵A经过拟上三角化后所得的矩阵A(n-1)；(3)对矩阵A(n-1)实行QR方法迭代结束后所得的矩阵；(4)矩阵A的全部特征值λi=(RI,Ii);(i=1,2,…,10)(5)A的相应于实特征值的特征向量。3.采用e型输出实型数，并且至少显示12位有效数字。算法：1）输入需要求解的矩阵A2）对上述生成的矩阵进行拟上三角化3）对拟上三角化后的矩阵进行QR分解（带双步位移）。4）用函数characteristic（）求解矩阵A（n-1）即A的所有特征值。5）用函数characteristicvector（）求解矩阵A（n-1）即A的所有特征向量。6）输出A的所有特征值λ、A的所有实特征值对应的特征向量、拟上三角矩阵A（n-1）、及其Q、R和R*Q说明：为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化。拟上三角化之后矩阵中包含有大量的0元素，在进行QR迭代时可以减少大量的0元素乘法运算，节省计算时间。且上三角矩阵进行QR分解之后，通过RQ得到新的矩阵A仍然是上三角矩阵，因此可以在整个迭代过程中使用。而带双步位移的QR分解法可以加快收敛速度，使矩阵在更短的迭代次数下完成收敛。源程序：#includestdio.h#includemath.h#includeconio.h#definee0.000000000001/*设置精度水平*/#defineN10/*设置矩阵阶数*/#defineL2500/*设置迭代次数*//*子程序*//*sgn函数*//*保证ci与aii符号相反*/doublesgn(doublea){doublex;if(ae)x=1;elsex=-1;returnx;}/*求根函数*//*计算双步位移法中的两个位移量，即右下角二阶子式的两个特征值*/doublepig(doubleb,doublec){doublem;m=b*b-4*c;returnm;}/*矩阵的拟上三角化*//*只要对A1进行拟上三角化就能保证每个矩阵Ak都为拟上三角矩阵，大大简化计算量*/voidhessenberg(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];}}}/*求解矩阵的所有特征值*/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;lL;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=pig(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=pig(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];}}}/*带双步位移的QR分解中M的分解*/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];}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