北航数值分析大作业二

数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1算法设计方案首先将矩阵A进行拟上三角化，把矩阵A进行QR分解，计算出RQ。要得出矩阵A的全部特征值，首先对A进行QR的双步位移得出特征值。最后，采用列主元的高斯消元法求解特征向量。1.1A的拟上三角化因为对矩阵进行QR分解并不改变矩阵的结构，因此在进行QR分解前对矩阵A进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。具体算法如下所示，记AA)1(，并记)(rA的第r列至第n列的元素为nrrjniarij,,1,;,,2,1)(。对于2,,2,1nr执行若nrriarir,,3,2)(全为零，则令)()1(rrAA,转5；否则转2。计算nririrrad12)(rrrrrrrrrrdcadac则取,0sgn)(,1)(,1若)(,12rrrrrracch令nTrnrrrrrrrrrRaacau)()(,2)(,1,,,,0,,0。计算rrTrrhuAp/)(rrrrhuAq/)(rrTrrhupt/rrrrutqTrrTrrrrpuuAA)()1(继续。1.2A的QR分解具体算法如下所示，记)1(1nAA，并记nnrijraA)(，令IQ1对于1,,2,1nr执行1.若nrriarir,,3,1)(全为零，则令rrQQ1rrAA1,转5；否则转2。2.计算nririrrad2)(rrrrrrrrrrdcadac则取,0sgn)(,1)(,1若)(,2rrrrrracch3.令nTrnrrrrrrrrrRaacau)()(,2)(,,,,,0,,0。4.计算rrTrrhuAp/)(rrrrhuAq/)(rrTrrhupt/rrrrutqTrrTrrrrpuuAA)()1(5.继续。1.3A的特征值为了加快收敛速度，采用带双步移位的QR分解求解A的全部特征值，具体算法如下，使用矩阵的拟上三角化算法把矩阵A化为你上三角阵；给定经度水平和最大迭代次数L。记，令k=1,m=n。如果，则得到A的一个特征值，置m:=m-1(降阶)，转(5)，否则转(4)。如果，则得到矩阵A的两个特征值,为二阶子阵)()(1,(k),1)(1,1kmmkmmmmkmmkaaaaD的两个特征根，置m:=m-1(降阶)，转(5)；否则转(6)。如果m=1，则得到A的一个特征根，转(9)，如果m=0，转(9)。否则转(3)。如果k=L，则计算失败。否则转(7)。记，计算kkTkkkkkkkkkkmmkmmkmmkmmkmmkmmQAQAQRMRQMItIsAAMaaaataas12)(,1)(1,)()(1,1)()(1,1)分解作对阶单位矩阵是()m(，转(3)。A的全部特征值以计算完毕，停止计算。其中，kM的QR分解与1kA的计算用下列算法实现：记kmmrijrmmijkACbBbMB1)()1(1,][,][。对于1,,2,1mr执行若mrribrir,,2,1)(全为零，则令rrrrCCBB11,,转5；否则转2。计算mririrrbd2)(rrrrrrrrrrdcadbc则取若,0sgn)(,1)()(2rrrrrrbcch令mTrmrrrrrrrrrRbbcbu)()(,1)(,,,,0,,0。计算rrTrrhuBv/TrrrrvuBB1rrTrrhuCp/rrrrhuCq/rrTrrhupt/rrrrutqTrrTrrrrpuuCC1继续。此算法执行完后，就得到mkCA1。1.4A的特征向量记λ为矩阵A的实特征值，x为对应的特征向量。则Ax=λx，即(A-λI)x=0的解即为矩阵A的特征向量。因此对于特征向量的求解可采用列主元素的Gauss消元法。其具体算法如下，消元过程对于1,,2,1nk执行选行号ki，使)()(maxkiknikkkiaak。交换)(kkja与),,1,()(nkkjakjik所含的数值。对于nkki,,2,1计算),,2,1(/)()()1()()(nkkjamaaaamkkiikkijkijkkkkikik回代过程1,,2,11)(1)(nnkaxaxxkkknkjjkkjkn最终得到的向量Tnxxx)1,,,(21的即为A对应于实特征值的特征向量。2.程序运行结果2.1拟上三角化后2.2QR分解Q阵R阵RQ阵2.3QR的双步位移求特征值和特征向量特征值特征向量3源程序#includestdio.h#includemath.h#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;i11;i++){for(j=1;j11;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(x0){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;r9;r++){pingfanghe=0;drpingfang=0;for(i=1;i11;i++){ur[i-1]=0;pr[i-1]=0;qr[i-1]=0;wr[i-1]=0;}tr=0;//置零for(i=r+2;i11;i++){pingfanghe=pingfanghe+a[i-1][r-1]*a[i-1][r-1];}if(pingfanghe==0){continue;//(1)}else{for(i=r+1;i11;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;i11;i++){ur[i-1]=a[i-1][r-1];}ur[r]=ur[r]-cr;//(3)urfor(j=1;j11;j++){for(i=1;i11;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;i11;i++){for(j=1;j11;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;i11;i++){tr=tr+pr[i-1]*ur[i-1];}tr=tr/hr;for(i=1;i11;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i11;i++){for(j=1;j11;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;i11;i++){for(j=1;j11;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;doubledrpingfang;doubledr,cr,hr;doubleur[10],wr[10],pr[10];for(i=1;i11;i++){for(j=1;j11;j++){if(i==j){Q[i-1][j-1]=1;}else{Q[i-1][j-1]=0;}}}for(r=1;r10;r++){pingfanghe=0;drpingfang=0;for(i=1;i11;i++){ur[i-1]=0;wr[i-1]=0;pr[i-1]=0;}for(i=r+1;i11;i++){pingfanghe=pingfanghe+a[i-1][r-1]*a[i-1][r-1];}if(pingfanghe==0){continue;}else{for(i=r;i11;i++){drpingfang=drpingfang+a[i-1][r-1]*a[i-1][r-1];}dr=sqrt(drpingfang);cr=-sgn(a[r-1][r-1])*dr;hr=cr*cr-cr*a[r-1][r-1];for(i=r;i11;i++){ur[i-1]=a[i-1][r-1];}ur[r-1]=ur[r-1]-cr;for(i=1;i11;i++){for(j=1;j11;j++){wr[i-1]=wr[i-1]+Q[i-1][j-1]*ur[j-1];}}for(i=1;i11;i++){for(j=1;j11;j++){Q[i-1][j-1]=Q[i-1][j-1]-wr[i-1]*ur[j-1]/hr;}}for(i=1;i11;i++){for(j=1;j11;j++){pr[i-1]=pr[i-1]+a[j-1][i-1]*ur[j-1]/hr;}}for(i=1;i11;i++){for(j=1;j11;j++){a[i-1][j-1]=a[i-1][j-1]-ur[i-1]*pr[j-1];}}}}printf(正交矩阵Q：\n);for(i=1;i11;i++){for(j=1;j11;j++){printf(%.12e\t,Q[i-1][j-1]);}printf(\n);}printf(上三角矩阵R:\n);for(i=1;i11;i++){for(j=1;j11;j++){printf(%.12e\t,a[i-1][j-1]);}printf(\n);}for(i=1;i11;i++){for(j=1;j11;j++){for(t=1;t11;t++){RQ[i-1][j-1]=RQ[i-1][j-1]+a[i-1][t-1]*Q[t-1][j-1];}}}//计算RQprintf(RQ阵：\n);for(i=1;i11;i++){for(j=1;j11;j++){printf(%.12e\t,RQ[i-1