数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1算法设计方案首先将矩阵A进行拟上三角化，把矩阵A进行QR分解，计算出RQ。要得出矩阵A的全部特征值，首先对A进行QR的双步位移得出特征值。最后，采用列主元的高斯消元法求解特征向量。1.1A的拟上三角化因为对矩阵进行QR分解并不改变矩阵的结构，因此在进行QR分解前对矩阵A进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。具体算法如下所示，记，并记的第r列至第n列的元素为。对于执行若全为零，则令,转5；否则转2。计算令。计算继续。1.2A的QR分解具体算法如下所示，记，并记，令对于执行1.若全为零，则令,转5；否则转2。2.计算3.令。4.计算5.继续。1.3A的特征值为了加快收敛速度，采用带双步移位的QR分解求解A的全部特征值，具体算法如下，使用矩阵的拟上三角化算法把矩阵A化为你上三角阵；给定经度水平和最大迭代次数L。记，令k=1,m=n。如果，则得到A的一个特征值，置m:=m-1(降阶)，转(5)，否则转(4)。如果，则得到矩阵A的两个特征值,为二阶子阵的两个特征根，置m:=m-1(降阶)，转(5)；否则转(6)。如果m=1，则得到A的一个特征根，转(9)，如果m=0，转(9)。否则转(3)。如果k=L，则计算失败。否则转(7)。记，计算，转(3)。A的全部特征值以计算完毕，停止计算。其中，的分解与的计算用下列算法实现：记。对于执行若全为零，则令,转5；否则转2。计算令。计算继续。此算法执行完后，就得到。1.4A的特征向量记λ为矩阵A的实特征值，x为对应的特征向量。则Ax=λx，即(A-λI)x=0的解即为矩阵A的特征向量。因此对于特征向量的求解可采用列主元素的Gauss消元法。其具体算法如下，消元过程对于执行选行号，使。交换与所含的数值。对于计算回代过程最终得到的向量的即为对应于实特征值的特征向量。2.程序运行结果2.1拟上三角化后2.2QR分解Q阵R阵RQ阵2.3QR的双步位移求特征值和特征向量特征值特征向量3源程序#include#include#definen11#defineE1e-12#defineL2500doublefuzhi(doublea[10][10]);//给矩阵A赋值doublenishangsanjiaohua(doublea[10][10]);doubleQR(doublea[10][10]);doublesgn(doublex);voidgauss();doubleA[10][10];doubletzR[n],I[n][n];intnR;voidQRshuangbu();voidmain(){fuzhi(A);nishangsanjiaohua(A);QR(A);QRshuangbu();gauss();}doublefuzhi(doublea[10][10])//给矩阵A赋值{inti,j;for(i=1;i<11;i++){for(j=1;j<11;j++){if(i==j){a[i-1][j-1]=1.52*cos(i+1.2*j);}else{a[i-1][j-1]=sin(0.5*i+0.2*j);}}}}doublesgn(doublex)//定义sgn函数{if(x>0){return1;}if(x<=0){return-1;}}doublenishangsanjiaohua(doublea[10][10]){intr,i,j;doublepingfanghe;doubledrpingfang;doubledr,cr,hr;doubleur[10],pr[10],qr[10],tr,wr[10];for(r=1;r<9;r++){pingfanghe=0;drpingfang=0;for(i=1;i<11;i++){ur[i-1]=0;pr[i-1]=0;qr[i-1]=0;wr[i-1]=0;}tr=0;//置零for(i=r+2;i<11;i++){pingfanghe=pingfanghe+a[i-1][r-1]*a[i-1][r-1];}if(pingfanghe==0){continue;//(1)}else{for(i=r+1;i<11;i++){drpingfang=drpingfang+a[i-1][r-1]*a[i-1][r-1];}dr=sqrt(drpingfang);//(2.1)drcr=-sgn(a[r][r-1])*dr;//(2.2)crhr=cr*cr-cr*a[r][r-1];//(2.3)hrfor(i=r+1;i<11;i++){ur[i-1]=a[i-1][r-1];}ur[r]=ur[r]-cr;//(3)urfor(j=1;j<11;j++){for(i=1;i<11;i++){pr[j-1]=pr[j-1]+a[i-1][j-1]*ur[i-1];}pr[j-1]=pr[j-1]/hr;}//(4.1)prfor(i=1;i<11;i++){for(j=1;j<11;j++){qr[i-1]=qr[i-1]+a[i-1][j-1]*ur[j-1];}qr[i-1]=qr[i-1]/hr;}//(4.2)qrfor(i=1;i<11;i++){tr=tr+pr[i-1]*ur[i-1];}tr=tr/hr;for(i=1;i<11;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i<11;i++){for(j=1;j<11;j++){a[i-1][j-1]=a[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}}}}for(i=1;i<11;i++){for(j=1;j<11;j++){printf("A%d,%d=%.12e\t",i,j,a[i-1][j-1]);}printf("\n");}}doubleQR(doublea[10][10]){doubleQ[10][10];doubleRQ[10][10];inti,j;intt;intr;doublepingfanghe;do...
VIP
