北航数值分析大作业第二题

数值分析第二次大作业史立峰SY1505327一、方案（1）利用循环结构将sin(0.50.2)()1.5cos(1.2)(){ijijijijija(i,j=1,2,……,10)进行赋值，得到需要变换的矩阵A；（2）然后，对矩阵A利用Householder矩阵进行相似变换，把A化为上三角矩阵A(n-1)。对A拟上三角化，得到拟上三角矩阵A(n-1)，具体算法如下：记A(1)=A，并记A(r)的第r列至第n列的元素为nrrjniarij,,1,;,,2,1)(。对于2,,2,1nr执行1.若nrriarir,,3,2)(全为零，则令A(r+1)=A(r),转5；否则转2。2.计算nririrrad12)(rrrrrrrrrrdcadac则取,0sgn)(,1)(,1若)(,12rrrrrracch3.令nTrnrrrrrrrrrRaacau)()(,2)(,1,,,,0,,0。4.计算rrTrrhuAp/)(rrrrhuAq/)(rrTrrhupt/rrrrutqTrrTrrrrpuuAA)()1(5.继续。（3）使用带双步位移的QR方法计算矩阵A(n-1)的全部特征值，也是A的全部特征值，具体算法如下：1.给定精度水平0和迭代最大次数L。2.记nnijnaAA][)1()1()1(，令nmk,1。3.如果)(1,kmma，则得到A的一个特征值)(kmma，置1:mm（降阶），转4；否则转5。4.如果1m，则得到的一个特征值)(11ka，转11；如果1m，则转3。5.求2阶子阵)()(1,(k),1)(1,1kmmkmmmmkmmkaaaaD的两个特征值1s和2s，即计算二次方程0det)()()(1.12kkmmkmmDaa的两个根1s和2s。6.如果2m，则得到A的两个特征值1s和2s，转11；否则转7。7.如果)(2,1kmma，则得到A的两个特征值1s和2s，置2:mm（降阶），转4；否则转88.如果Lk，则计算终止，未得到A的全部特征值；否则转9。9.记),1(][)()(mjiaAmmkijk，计算kkTkkkkkkkkkkmmkmmkmmkmmkmmkmmQAQAQRMRQMItIsAAMaaaataas12)(,1)(1,)()(1,1)()(1,1)分解作对阶单位矩阵是()m(10.置1:kk，转3。11.A的全部特征值已计算完毕，停止计算。其中，kM的QR分解与1kA的计算用下列算法实现：记kmmrijrmmijkACbBbMB1)()1(1,][,][。对于1,,2,1mr执行1.若mrribrir,,2,1)(全为零，则令rrrrCCBB11,,转5；否则转2。2.计算mririrrbd2)(rrrrrrrrrrdcadbc则取,0sgn)(,1)(若)(2rrrrrrbcch3.令mTrmrrrrrrrrrRbbcbu)()(,1)(,,,,0,,0。4.计算rrTrrhuBv/TrrrrvuBB1rrTrrhuCp/rrrrhuCq/rrTrrhupt/rrrrutqTrrTrrrrpuuCC15.继续。此算法执行完后，就得到mkCA1。（4）计算Q,R，一边求R*Q矩阵。记阶单位阵）n（令,][，并记A1)(r1IQaAAnnrij对于1n,,2,1r执行1.若nrriarir,,3,2)(全为零，则令rrrrAAQQ11转5；否则转2.2.计算mririrrbd2)(rrrrrrrrrrdcadbc则取,0sgn)(,1)(若)(2rrrrrrbcch3.令mTrmrrrrrrrrrRbbcbu)()(,1)(,,,,0,,0。4.计算TrrrrrrTrrTrrrrrrpuAAhuAphuQQuQ1r1r5.继续当此算法执行完毕后，就得到正交矩阵nQQ和上三角矩阵R6.然后计算出QR矩阵（5）用列主元素Gauss消去法计算矩阵A对应于实特征值的特征向量，具体算法如下：记)()1(的实特征值为AIAA1.消元过程对于1,,2,1nk执行(1)选行号ki，使)()(maxkiknikkkiaak。(2)交换)(kkja与),,1,()(nkkjakjik所含的数值。(3)对于nkki,,2,1计算),,2,1(/)()()1()()(nkkjamaaaamkkiikkijkijkkkkikik2.回代过程1,,2,11)(1)(nnkaxaxxkkknkjjkkjkn最终得到的向量Tnxxx)1,,,(21的即为A对应于实特征值的特征向量。二、源程序#include#include#include#definen10#defineEvoidHouseholder(doublea[n][n]);//拟上三角化函数doublesgn(doublea);//符号函数voidQRfenjie(doublea[n][n],doubleQ[n][n],doubleR[n][n]);voidQR(doublea[n][n],doubleL[n][2]);//带双步位移的QR分解voidMxM(doubleM[n][n],doubleA[n][n],doubleB[n][n],intm);//矩阵相乘voidsolve(doublea[n][n],doubles1[2],doubles2[2],intm);//解方程函数voidGauss(doublea[n][n],doublex[n]);//定义列主元高斯消去法函数voidmain(){inti,j,k;doublea[n][n],qr[n][n],q[n][n],r[n][n],L[n][2],x[n];for(i=0;in;i++)//录入要进行变换的矩a，并对q初始赋值。{for(j=0;jn;j++){if(i==j)a[i][j]=*cos(i+1+*(j+1));elsea[i][j]=sin*(i+1)+*(j+1));}}//coutendl;Householder(a);//调用拟上三角化函数coutendlendl对矩阵A拟上三角化前七列的结果:endl;for(i=0;in;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout矩阵A的第i+1行元素为：endl;for(j=0;j7;j++){if(fabs(a[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)a[i][j];//k++;//coutendl;//if(k%3==0)//判断每行是否到达三个元素}coutendl;}cout对矩阵A拟上三角化后三列的结果:endl;for(i=0;in;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout矩阵A的第i+1行元素为：endl;for(j=7;jn;j++){if(fabs(a[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)a[i][j];//k++;//if(k%3==0)//判断每行是否到达三个元素//coutendl;}coutendl;}coutendl;QRfenjie(a,q,r);QR(a,L);coutendl对矩阵A进行QR分解后所得矩阵前七列的结果:endl;for(i=0;in;i++){for(j=0;j7;j++){if(fabs(a[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)a[i][j];}coutendl;}cout对矩阵A进行QR分解后所得矩阵三列的结果:endl;for(i=0;in;i++){for(j=7;jn;j++){if(fabs(a[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)a[i][j];}coutendl;}/*cout对矩阵A进行QR分解后所得矩阵:endl;for(i=0;in;i++){k=0;cout矩阵A的第i+1行元素为：endl;for(j=0;jn;j++){if(fabs(a[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)a[i][j];k++;if(k%3==0)coutendl;}coutendl;}*/for(i=0;in;i++){for(j=0;jn;j++)for(k=0,qr[i][j]=0;kn;k++)qr[i][j]=qr[i][j]+r[i][k]*q[k][j];}coutendlR*Q相乘的前七列:endl;for(i=0;in;i++){//k=0;//coutR*Q的第i+1行元素为：endl;for(j=0;j7;j++){if(fabs(qr[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)qr[i][j];//k++;//if(k%3==0)//coutendl;}coutendl;}coutendlR*Q相乘的后三列:endl;for(i=0;in;i++){for(j=7;jn;j++){if(fabs(qr[i][j])E)a[i][j]=0;coutsetprecision(12)setiosflags(ios::scientific)qr[i][j];}coutendl;}coutendl矩阵A的全部特征值为endl;for(i=0;in;i++){if(L[i][1]==0)coutλi+1=L[i][0]endl;elsecoutλi+1=L[i][0]+L[i][1]iendl;}coutendl实特征值对应的特征向量分别是：endl;for(k=0;kn;k++){if(L[k][1]==0)//判断特征值是否为实特征值{for(i=0;in;i++)//重新录入矩阵A{for(j=0;jn;j++){if(i==j)a[i][j]=*cos(i+1+*(j+1))-L[k][0];elsea[i][j]=sin*(i+1)+*(j+1));}}Gauss(a,x);coutλk+1=L[k][0]对应的特征向量是:endl;for(j=0;jn;j++){coutx[j]endl;}}elsecontinue;}}voidHouseholder(doublea[n][n])//拟上三角化函数{inti,j,k;intm=0;doubled,c,h,t;doubleu[n],p[n