北航数值分析报告大作业二51639.doc

北京航空航天大学 2014年研究生 《数值分析》作业 第二题 院系：XXXXXXXXX学院 学号：XXXXXXXXX 姓名：XXXX Email:XXXXXXXXXXXXXXX 2014年11月 算法设计方案 1、矩阵A的输入与题目void assignment()子程序。 2、将矩阵A拟上三角化 利用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵，由子程序void Hessenberg()完成。 矩阵A的拟上三角化算法： 记A，并记。对r=1，2，，n-2执行 （1）若全为零，则令，转（5）；否则转（2）。 （2）计算，（若，则取）， （3）令 （4）计算，，，， （5）继续 3、对拟上三角化后的矩阵进行QR分解 为了了解普通的QR分解过程及结果，下述程序中用void QRfenjie()子程序来对拟上三角化过后的A阵进行QR分解，并输出Q阵R阵RQ阵。其中QRfenjie()既用于基本的QR分解和输出Q阵、R阵、RQ阵，又用于带双步位移的QR方法。m为维数，m=10。当n=0时，子程序进行基本QR分解；n≠0时，子程序进行带双步位移的QR分解。 4、对拟上三角化后的矩阵进行带双步位移的QR分解。 子程序void QRfa()实现对拟上三角化后的A阵进行带双步位移的QR分解，得出特征值并输出，并用子程序void vector()对其中的实数特征值进行求解，得出对应的特征向量。 带双步位移的QR方法具体算法如下： （1）使用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵；给定精度水平和迭代最大次数L。 （2）记，令k=1，m=n。 （3）如果，则得到A的一个特征值，置m=n-1，转（4）；否则转（5）。 （4）如果m=1，则得到A的一个特征值，转（11）；如果m=0，则直接转（11）；如果m1，则转（3）。 （5）求二阶子阵的两个特征值和，即计算二次方程的两个根和。 （6）如果m=2，则得到A的两个特征值和，转（11）；否则转（7）。 （7）如果，则得到A的两个特征值和，置m=m-2，转（4）；否则转（8）。 （8）如果k=L，则计算终止，未得到A的全部特征值，否则转（9）。 （9）记，计算： （为m阶单位矩阵） （对作QR分解） （10）置k=k+1，转（3）。 （11）A的全部特征值以计算完毕，停止计算。 ：作QR分解与的计算算法为： 记，，。对于r=1，2，，m-1执行 （1）若全为零，则令，转（5）；否则转（2）。 （2）计算，（若，则取）， （3）令 （4）计算，，，，，， （5）继续 源程序 #includestdio.h #includeiostream #includestdlib.h #includemath.h #includefloat.h #includeiomanip #includetime.h #define N 10 /*定义矩阵的阶数*/ #define E 1.0e-12 /*定义全局变量相对误差限*/ #define L 20 /*迭代次数*/ FILE *fp; struct C /*定义结构体*/ { double R; double I; }; void shuchu(double a[N][N]) /*输出一个10*10的矩阵*/ { int i,j; double b=0; for(i=0;iN;i++) { for(j=0;jN;j++) { if(fabs(a[i][j])E) fprintf(fp,%+.12e,b); else fprintf(fp,%+.12e,a[i][j]); if((j+1)%3==0) fprintf(fp,\n); } fprintf(fp,\n); } } int sign(double x) /*符号函数*/ { if(x0) return -1; if(x == 0) return 0; if(x0) return 1; } void assignment(double a[N][N]) /*输入矩阵A*/ { int i,j; for(i=0;iN;i++) { for(j=0;jN;j++) { if(i==j) a[i][j]=1.52*cos((i+1)+1.2*(j+1)); else a[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } } } void Hessenberg(double a[N]