北航数值分析大作业第二题答案.docx

目标：使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： (i,j=1,2,10)开始输入矩阵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的实特征值对应的一个特征向量的求法：由于，因此在经过选主元的高斯消元以后，即的最后一行必然为零，左上方变为n-1阶单位矩阵，右上方变为n-1阶向量，然后令，则。这样即求出所有A所有实特征值对应的一个特征向量。#include#include#include#define N 10#define E 1.0e-12#define MAX 10000/以下是符号函数double sgn(double a)double z;if(aE) z=1;else z=-1;return z; /以下是矩阵的拟三角分解 void nishangsanjiaodiv(double ANN) int i,j,k;int m=0;double d,c,h,t;double uN,pN,qN,wN;for(i=0;iN-2;i+) for(j=i+2;jN;j+) if(Aji=E) m=m+1;if(m=(N-2-i) continue; for(j=i+1,d=0;jN;j+) d=d+Aji*Aji;d=sqrt(d);c=-1*sgn(Ai+1i)*d;h=c*c-c*Ai+1i;for(j=i+2;jN;j+) uj=Aji;for(j=0;ji+2;j+) uj=0;ui+1=Ai+1i-c;for(j=0;jN;j+) for(k=i+1,pj=0;kN;k+) pj=Akj*uk+pj;pj=pj/h;for(j=0;jN;j+) for(k=i+1,qj=0;kN;k+) qj=Ajk*uk+qj;qj=qj/h;for(j=0,t=0;jN;j+) t=t+pj*uj;t=t/h;for(j=0;jN;j+) wj=qj-t*uj;for(j=0;jN;j+)for(k=0;kN;k+) Ajk=Ajk-wj*uk-uj*pk; /以下是矩阵的QR分解 void qrdiv(double ANN,double QNN,double RNN) int i,j,k;/int m=0;double d,c,h;double uN,wN,pN;for(i=0;iN;i+)for(j=0;jN;j+) if (i=j) Qij=1; else Qij=0;for(i=0;iN;i+)for(j=0;jN;j+) Rij=Aij;for(i=0;iN-1;i+)/for(j=i+1;jN;j+) if(Rji=E) m=m+1;/if(m=(N-1-i) continue; for(j=i,d=0;jN;j+) d=d+Rji*Rji;d=sqrt(d);c=-1*sgn(Rii)*d;h=c*c-c*Rii;for(j=i+1;jN;j+) uj=Rji;for(j=0;ji;j+) uj=0;ui=Rii-c;for(j=0;jN;j+) for(k=0,wj=0;kN;k+) wj=Qjk*uk+wj;for(j=0;jN;j+) for(k=0;kN;k+) Qjk=Qjk-wj*uk/h;for(j=0;jN;j+) for(k=i,pj=0;kN;k+) pj=Rkj*uk+pj;pj=pj/h;for(j=0;jN;j+)for(k=0;kN;k+) Rjk=Rjk-uj*pk;/矩阵的QR分解 /以下是二次多项式求根double root(double b,double c) double m;m=b*b-4*c;return m; /二次多项式求根 /以下是求解矩阵的所有特征值void characteristic(double ANN,double chaRN,double chaIN) int k=0,m=N-1; int i,j; int L; double s,t,x; double MNN,BNN; int f=0; double d,c,h; double uN,wN,pN; double QNN,RNN; for(L=0;LMAX;L+) next: if (m=0) chaR0=A00;chaI0=0;break; if(fabs(Amm-1)=E) x=sqrt(x);chaRm=s/2+x/2;chaRm-1=s/2-x/2;chaIm=0;chaIm-1=0; else x=sqrt(fabs(x);chaRm=s/2;chaRm-1=s/2;chaIm=x/2;chaIm-1=-x/2; break; if(fabs(Am-1m-2)=E) x=sqrt(x);chaRm=s/2+x/2;chaRm-1=s/2-x/2;chaIm=0;chaIm-1=0; else x=sqrt(fabs(x);chaRm=s/2;chaRm-1=s/2;chaIm=x/2;chaIm-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,Mij=0;k=m;k+) Mij=Aik*Akj+Mij; Mij=Mij-s*Aij+t; else for(k=0,Mij=0;k=m;k+) Mij=Aik*Akj+Mij; Mij=Mij-s*Aij; /以下是M的QR分解 for(i=0;i=m;i+) for(j=0;j=m;j+) if (i=j) Qij=1; else Qij=0; for(i=0;i=m;i+) for(j=0;j=m;j+) Rij=Mij; for(i=0;im;i+) for(j=i+1;j=m;j+) if(Rji=E) f=f+1; if(f=(m-i) continue; for(j=i,d=0;j=m;j+) d=d+Rji*Rji; d=sqrt(d); c=-1*sgn(Rii)*d; h=c*c-c*Rii; for(j=i+1;j=m;j+) uj=Rji; for(j=0;ji;j+) uj=0; ui=Rii-c; for(j=0;j=m;j+) for(k=0,wj=0;k=m;k+) wj=Qjk*uk+wj; for(j=0;j=m;j+) for(k=0;k=m;k+) Qjk=Qjk-wj*uk/h; for(j=0;j=m;j+) for(k=i,pj=0;k=m;k+) pj=Rkj*uk+pj;pj=pj/h; for(j=0;j=m;j+) for(k=0;k=m;k+) Rjk=Rjk-uj*pk; for(j=0;j=m;j+) for(k=0;k=m;k+) Mjk=Qjk; /以上是M的QR分解 for(i=0;i=m;i+) for(j=0;j=m;j+) for(k=0,Bij=0;k=m;k+) Bij=Mki*Akj+Bij; for(i=0;i=m;i+) for(j=0;j=m;j+) for(k=0,Aij=0;k=m;k+) Aij=Bik*Mkj+Aij; /以下是求矩阵的所有特征值的特征向量void eigenvector(double VNN,double TN) double ANN,baozNN,guodN;double c;int i,j,k,m,t;int W=0;for(i=0;iN;i+) for(j=0;jN;j+) baozij=Vij; for(t=0;t6;t+) for(i=0;iN;i+) for(j=0;jN;j+) Aij=baozij;for(i=0;iN;i+) Aii=Aii-Tt;for(i=0;iN-1;i+) for(j=i;jE) k=j;break; for(j=i;jN;j+) guodj=Aij;Aij=Akj;Akj=guodj;for(j=i;jE) for(m=i;mN;m+) Ajm=Ajm/c;for(j=0;jN;j+) c=Aji;if(j!=i)for(m=i;m=0;i-) Vti=AiN-1; /以下是主函数void main()double aNN,bNN,chaRN,chaIN; double qNN,rNN,qrNN;double shiyanN;double f,g;int i,j,k;for(i=0;iN;i+)for(j=0;jN;j+) if(i!=j) aij=sin(0.5*(i+1)+0.2*(j+1);else aij=1.5*cos(i+1)+1.2*(j+1);nishangsanjiaodiv(a);printf(矩阵A的拟上三角分解：n);for(i=0;iN;i+)for(j=0;jN-5;j+) if (fabs(aij)E) aij=0; printf(%22.11e,aij); printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+) if (fabs(aij)E) aij=0; printf(%22.11e,aij); printf(n);printf(n); qrdiv(a,q,r);printf(n);printf(n);printf(n);printf(拟上三角矩阵A的QR分解：n);printf(上三角矩阵R：n); for(i=0;iN;i+)for(j=0;jN-5;j+) if (fabs(rij)E) rij=0;printf(%22.12e,rij);printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+) if (fabs(rij)E) rij=0;printf(%22.12e,rij);printf(n);printf(n);printf(正交矩阵Q：n); for(i=0;iN;i+)for(j=0;jN-5;j+) if (fabs(qij)E) qij=0;printf(%22.12e,qij); printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+) if (fabs(qij)E) qij=0;printf(%22.12e,qij); printf(n);printf(n); for(i=0;iN;i+)for(j=0;jN;j+) for(k=0,qrij=0;kN;k+) qrij=qrij+rik*qkj;printf(n);printf(n);printf(n);printf(R*Q：n); for(i=0;iN;i+)for(j=0;jN-5;j+) if (fabs(qrij)E) qrij=0;printf(%22.12e,qrij); printf(n);printf(n);printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+) if (fabs(qrij)E) qrij=0;printf(%22.12e,qrij); printf(n);printf(n);printf(n);printf(n);characteristic(a,chaR,chaI);for(i=1;iN;i+) if (i3) f=chaRi;g=chaIi;chaRi=chaR7+i;chaIi=chaI7+i;chaR7+i=f;chaI7+i=g;if (i=5) f=chaRi;g=chaIi;chaRi=chaR7;chaIi=chaI7;chaR7=f;chaI7=g; printf(矩阵A所有特征值：n);for(j=0;jN;j+) if(fabs(chaIj)E) printf(%2d =%18.11e +%18.11ein,j+1,chaRj,chaIj);else printf(%2d =%19.11e %19.11ein,j+1,chaRj,chaIj);printf(n);printf(n);printf(n);for(i=0;iN;i+)for(j