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

实用文档数值分析第二题目录数值分析第二题...........................................11引言：.............................................11.1矩阵的拟上三角化.............................11.2矩阵的特征值求解.............................21.3矩阵的特征向量求解...........................42算法的程序实现.....................................52.1主程序.......................................52.2子程序的实现.................................73计算结果...........................................83.2矩阵Q、R以及乘积RQ..........................93.3各实特征值及其相对应的特征向量..............104实验结论..........................................12附录源代码.........................................12《数值分析》计算实习题目二11引言1.1矩阵的拟上三角化为了减少计算量，对矩阵A利用Householder矩阵进行相似变换，把A化为上三角矩阵𝐴(n−1)。对A拟上三角化，得到拟上三角矩阵)1(nA，具体算法如下：记AA)1(，并记)(rA的第r列至第n列的元素为nrrjniarij,,1,;,,2,1)(。对于2,,2,1nr执行1.若nrriarir,,3,2)(全为零，则令)()1(rrAA,转5；否则转2。2.计算nririrrad12)(rrrrrrrrrrdcadac则取,0sgn)(,1)(,1若)(,12rrrrrracch3.令nTrnrrrrrrrrrRaacau)()(,2)(,1,,,,0,,0。4.计算rrTrrhuAp/)(rrrrhuAq/)(rrTrrhupt/rrrrutqTrrTrrrrpuuAA)()1(《数值分析》计算实习题目二25.继续。1.2矩阵的特征值求解使用带双步位移的QR方法计算矩阵)1(nA的全部特征值，也是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，计算《数值分析》计算实习题目二3kkTkkkkkkkkkkmmkmmkmmkmmkmmkmmQAQAQRMRQMItIsAAMaaaataas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/《数值分析》计算实习题目二4rrrrutqTrrTrrrrpuuCC15.继续。此算法执行完后，就得到mkCA1。1.3矩阵的特征向量求解用列主元素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对应于实特征值的特征向量。《数值分析》计算实习题目二52算法的程序实现2.1主程序#includeiostream#includematrix-tool.h#includesimi-trian.h#includeQR_decompose.h#includeGauss_by_max.h#definethrehold(1e-12)usingnamespacestd;intn=10;externdouble*a;externdouble*Mk;externdouble*a_lambda;externdouble*b_lambda;double*spec_u;//高斯求解方程组的右端常数列向量double*G_b;intmain(){simi_Initial();cout处理之前的A矩阵为：endl;Print_Matrix(a,n);cout拟上三角化后的A矩阵为：endl;pseudo_triangle(a,n);《数值分析》计算实习题目二6BI_STEP_QR(a,n);/******************特征向量求取代码段*****************///初始化spec_u[10]spec_u=newdouble[n];G_b=newdouble[n];chari;simi_Initial();//为G_b赋值for(i=0;in;i++){(*(G_b+i))=0;}for(i=0;in;i++){if((*(b_lambda+i))==0){Gauss_root(a,G_b,n,(*(a_lambda+i)),spec_u);cout属于(*(a_lambda+i))i(*(b_lambda+i))的特征向量为：endl;Print_Vector(spec_u,n);}}system(pause);return0;}《数值分析》计算实习题目二72.2子程序的实现2.2.1“matrix-tool.h”voidMatrix_M_Vector(double*A,double*b,double*y,intn);函数实现N*N的矩阵乘以N*1的向量，得到N*1向量YvoidMatrix_Converted_M_Vector(double*A,double*b,double*y,intn);函数实现转置后N*N矩阵乘以N*1向量,得到N*1向量YvoidColumn_M_Row(double*a,double*b,double*A,intn);函数实现列乘以行得到一个矩阵AvoidRow_M_Column(double*a,double*b,double*y,intn);函数实现行乘以列得到一个数值yvoidPrint_Matrix(double*a,intn);函数实现将矩阵a显示出来voidPrint_Vector(double*b,intn);函数实现将向量b显示出来2.2.2“Gauss_by_max.h”voidGauss_root(double*a,double*b,intn,doublelambda,double*u)这个函数用于使用高斯列主元消元法实现方程组的求解。M_switch(double*a,double*b,intn,intk,intflag)这个函数用于实现按照列主元的绝对值大小重新排列方程组2.2.3“QR_decompose.h”voidBI_STEP_QR(double*a,intn)这个函数用于实现矩阵的双步位移QR分解。《数值分析》计算实习题目二82.2.4“simi-trian.h”voidpseudo_triangle(double*A,intn);这个函数主要实现矩阵的拟上三角化voidsimi_Initial();这个函数主要用于按照题目条件实现矩阵的初始化3计算结果3.1拟三角化后的矩阵An)1(如图1所示：图1拟三角化后的矩阵A《数值分析》计算实习题目二93.2矩阵Q、R以及乘积RQ图2QR分解后得到的RQ矩阵图3QR分解后得到的Q矩阵《数值分析》计算实习题目二10图4QR分解后得到的R矩阵3.3各实特征值及其相对应的特征向量得到的特征向量如下图5所示：图5计算得到A的特征值《数值分析》计算实习题目二11得到的A的各特征值对应的特征向量为：图6A的各特征值对应的特征向量图7A的各特征值对应的特征向量《数值分析》计算实习题目二124实验结论本次作业运用带双步位移的QR分解法求出了给定矩阵的全部特征值和对应于实特征值的特征向量，通过本人调试后，达到了题目的设计要求。本次调试过程中，还是发现了比较多的问题。1.当一个数据被多个程序段进行操作的时候，需要在具体程序中，对被操作的数据另行创建副本，以避免这种操作改变了原始数据而影响了其他的程序。我在本次程序设计过程中，就出现这种问题：在对矩阵A进行QR分解求取特征值的程序和后来求取特征向量的程序中，都没有创建原始变量的副本，导致程序之间互相干扰。2.在调试过程中，由于我是使用指针来操作数组的，所以整个程序对二维数组的维数非常敏感。例如我在调试过程中，在完成矩阵A的QR分解过程中，原来是通过语句来(∗(MkC+i∗m+j))=(∗(a+i∗m+j))实现的，最后发现程序在m=10时，程序运行正确，但是当m=8后，程序就出现问题了。最后才知道这句赋值语句出现了问题，应该是(∗(MkC+i∗m+j))=(∗(a+i∗n+j))。3.在调试过程中，我发现我写的函数间的耦合性非常大，所使用的全局变量也特别多。导致整个程序调试起来非常不方便。以后可以通过减少全局变量的个数等其他方式来降低程序的耦合性。附录源代码1.“matrix-tool.c”#includematrix-tool.h#includemath.h#includeiostream#includeiomanip《数值分析》计算实习题目二13usingnamespacestd;voidMatrix_M_Vector(double*A,double*b,double*y,intn){inti,j;doublesum=0;for(i=0;in;i++){for(j=0;jn;j++){sum+=(*(A+i*n+j))*(*(b+j));}(*