数值分析报告第二大题

《数值分析》计算实习报告第二大题学号：DY1305姓名：指导老师：一、题目要求已知10*10阶的矩阵A，使用带双步位移的QR分解法求其全部特征值，及每个实特征值对应的特征向量。此外还有拟上三角化产生的矩阵(1)nA，及其进行QR分解结束后所得的矩阵Q、R以及乘积矩阵RQ。二、算法设计方案先用课本第61到62页的算法将矩阵进行拟上三角化，得到拟上三角化矩阵(1)nA。再用课本第63到64页的带双步位移的QR算法求矩阵(1)nA的全部特征值，同时可以得到矩阵A进行QR分解后的矩阵Q、R以及乘积矩阵RQ。由于矩阵A与矩阵(1)nA相似，则他们的特征值相同，即可得矩阵A的全部特征值。对于各实特征值对应的特征向量，根据()0iAIx，利用实特征值选主元Gauss消元法求解i对应的特征向量。三、运算结果及分析以下计算结果截图。图1（a）对A进行拟上三角化得到的矩阵A(n-1)图1（b）对A进行拟上三角化得到的矩阵A(n-1)图2对A进行QR得到的矩阵Q图3对A进行QR得到的矩阵R图4对A进行QR得到的乘积矩阵RQ图5矩阵A的全部特征值图6各实特征值对应的特征向量四、源程序#includeiostream.h#includestdlib.h#includemath.h#includefloat.h#includeiomanip.h#defineEpsilon1e-12//精度#defineN11#definen10#defineL1000//迭代次数doublea[N][N];doublelambda[N][2];doublefuzhi();//赋值函数voidNSSJH();//拟上三角化函数voidQR();//QR分解函数voidDoubleQR();//带双步位移的QR分解函数voidTZXL();//求取实特征根对应特征向量的函数voidmain(){voidNSSJH();voidQR();DoubleQR();TZXL();}doublefuzhi()//赋值函数，存储行列均为1-10{inti,j;for(i=1;iN;i++){for(j=1;jN;j++){if(i==j){a[i][j]=1.5*cos(i+1.2*j);}elsea[i][j]=sin(0.5*i+0.2*j);}}returna[i][j];}voidNSSJH()//拟上三角化函数{inti,j,r;doublecr,dr,hr,tr,temp;doublepr[N]={0},qr[N]={0},ur[N],wr[N]={0};fuzhi();for(r=1;r=n-2;r++){dr=0;cr=0;hr=0;tr=0;temp=0;for(i=r+1;i=n;i++){temp=temp+a[i][r]*a[i][r];}dr=sqrt(temp);//求解drif(a[r+1][r]0)cr=-dr;elsecr=dr;//求解crhr=cr*cr-cr*a[r+1][r];//求解hrfor(i=1;i=r;i++){ur[i]=0;}ur[r+1]=a[r+1][r]-cr;for(i=r+2;i=n;i++)//生成ur[N]{ur[i]=a[i][r];}for(i=1;i=n;i++){temp=0;for(pr[i]=0.0,j=1;j=n;j++)//求解pr[N]{pr[i]=pr[i]+a[j][i]*ur[j]/hr;}}for(i=1;i=n;i++)//求解qr[N]{temp=0;for(qr[i]=0.0,j=1;j=n;j++){qr[i]=qr[i]+a[i][j]*ur[j]/hr;}}for(i=1;i=n;i++)//求解tr{tr=tr+ur[i]*pr[i]/hr;}for(i=1;i=n;i++)//求解wr[N]{wr[i]=qr[i]-ur[i]*tr;}for(i=1;i=n;i++)//求解A(r+1){for(j=1;j=n;j++){a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];}}}cout******对A进行拟上三角化得到的矩阵A(n-1)******endl;for(i=1;i=n;i++){for(j=1;j=n;j++){coutA[i][j]=setprecision(12)setiosflags(ios::scientific)a[i][j]'\t';if(j%2==0)coutendl;}}}voidQR()//QR分解函数{inti,j,k,r;doubledr,cr,hr,temp,sum;doubleQ[N][N]={0.0},wr[N]={0.0},ur[N]={0.0},pr[N]={0.0},A[N][N]={0.0},b[N]={0.0};for(i=1;i=n;i++){for(j=1;j=n;j++)//生成Q[N]{if(i==j)Q[i][j]=1;elseQ[i][j]=0;}}NSSJH();for(r=1;rn;r++){sum=0;for(i=r+1;i=n;i++){sum=sum+a[i][r]*a[i][r];}if(sum==0)continue;//如果a[i][r]（i=r+1,……n）全为零，r++，继续循环else{sum=sum+a[r][r]*a[r][r];dr=sqrt(sum);//求解drif(a[r][r]0)cr=-dr;elsecr=dr;//求解crhr=cr*cr-cr*a[r][r];//求解hrfor(i=1;ir;i++){ur[i]=0;}ur[r]=a[r][r]-cr;for(i=r+1;i=n;i++)//生成ur[N]{ur[i]=a[i][r];}for(i=1;i=n;i++){temp=0;for(j=1;j=n;j++){temp=temp+Q[i][j]*ur[j];}wr[i]=temp;//求解wr[N]}for(i=1;i=n;i++){for(j=1;j=n;j++){Q[i][j]=Q[i][j]-wr[i]*ur[j]/hr;}}//求解Q[N]，即为所求矩阵Qfor(i=1;i=n;i++){temp=0;for(j=1;j=n;j++){temp=temp+a[j][i]*ur[j]/hr;}pr[i]=temp;//求解pr[N]}for(i=1;i=n;i++){for(j=1;j=n;j++){a[i][j]=a[i][j]-ur[i]*pr[j];}}//求解a[N][N]，即为所求R矩阵for(i=1;i=n;i++){for(j=1;j=n;j++)//求解A[N]，即为所求乘积矩阵RQ{temp=0;for(k=1;k=n;k++){temp+=a[i][k]*Q[k][j];}A[i][j]=temp;}}}}cout******对A进行QR得到的矩阵Q:******endl;//输出Q、R、RQfor(i=1;i=n;i++){for(j=1;j=n;j++){coutQ[i][j]=setprecision(12)setiosflags(ios::scientific)Q[i][j]’\t’;if(j%3==0)coutendl;}}cout******对A进行QR得到的矩阵R:******endl;for(i=1;i=n;i++){for(j=1;j=n;j++){coutR[i][j]=setprecision(12)setiosflags(ios::scientific)a[i][j]’\t’;if(j%3==0)coutendl;}}cout******对A进行QR得到的乘积矩阵RQ:******endl;for(i=1;i=n;i++){for(j=1;j=n;j++){coutQR[i][j]=setprecision(12)setiosflags(ios::scientific)A[i][j]’\t’;if(j%3==0)coutendl;}}}voidDoubleQR()//带双步位移的QR分解函数{inti,j,k,m,r,l;doubleb,c,d,s,t,M[N][N]={0.0},I[N][N]={0.0};doublecr,dr,hr,tr,temp,sum;doublepr[N]={0.0},qr[N]={0.0},ur[N]={0.0},wr[N]={0.0},vr[N]={0.0};k=1;m=10;NSSJH();//先进行拟上三角化for(i=1;i=m;i++){for(j=1;j=m;j++){if(i==j)I[i][j]=1;elseI[i][j]=0;}}//定义单位矩阵loop1:{//降一阶的判断if(fabs(a[m][m-1])=Epsilon){lambda[m][0]=a[m][m];lambda[m][1]=0;m=m-1;gotoloop2;}elsegotoloop3;}loop2:{//判断m=2与否if(m==2){c=a[m-1][m-1]+a[m][m];b=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];d=c*c-4*b;if(d=0){lambda[m][0]=(c+sqrt(d))/2;lambda[m][1]=0;lambda[m-1][0]=(c-sqrt(d))/2;lambda[m-1][1]=0;}else{lambda[m][0]=c/2;lambda[m][1]=sqrt(-d)/2;lambda[m-1][0]=c/2;lambda[m-1][1]=-sqrt(-d)/2;}gotoloopjs;}elseif(m==1){lambda[1][0]=a[1][1];lambda[1][1]=0;gotoloopjs;}elseif(m==0)gotoloopjs;elsegotoloop1;}loop3:{//降两阶的判断if(fabs(a[m-1][m-2])=Epsilon){c=a[m-1][m-1]+a[m][m];b=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];d=c*c-4*b;if(d=0){lambda[m][0]=(c+sqrt(d))/2;lambda[m][1]=0;lambda[m-1][0]=(c-sqrt(d))/2;lambda[m-1][1]=0;}else{lambda[m][0]=c/2;lambda[m][1]=sqrt(-d)/2;lambda[m-1][0]=c/2;lambda[m-1][1]=-sqrt(-d)/2;}m=m-2;gotoloop2;}elsegotoloopPD;}loopPD:{//判断迭代次数是否超出规定if(k==L){cout******计算终止，未得到A的全部特征值******\nendl;gotoloopjs;}elsegotoloopqr;}loopqr:{s=a[m-1][m-1]+a[m][m];//带双步位移的QR分解t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];for(i=1;i=m;i++){for(j=1;j=m;j++){temp=0;for(l=1;l=m;l++){temp=temp+a[i][l]*a[l][j];}M[i][j]=temp-s*a[i][j]+t*I[i][j];}//生成M[N][N]}for(r=1;rm;r++){sum=0;for(i=r+