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

《数值分析B》大作业二 SY1103120 朱舜杰 算法设计方案： 1.矩阵A的存储 首先输入需要求解的矩阵A。即为下述程序中的void assignment()子程序。 2.将矩阵A拟上三角化 使用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵，由子程序void Hessen()完成。 矩阵A的拟上三角化的算法如下： 记A，并记。对于r=1，2，，n-2执行 （a）若全为零，则令，转（e）；否则转（b）。 （b）计算，（若，则取）， （c）令 （d）计算，，，， （e）继续 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 xiangliang()对其中的实数特征值进行求解，得出对应的特征向量。 带双步位移的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分解） 【对作QR分解与的计算算法如下： 记，，。对于r=1，2，，m-1执行 （a）若全为零，则令，转（e）；否则转（b）。 （b）计算，（若，则取）， （c）令 （d）计算，，，，，， （e）继续】 （10）置k=k+1，转（3）。 （11）A的全部特征值以计算完毕，停止计算。  源程序 #includestdio.h #includeiostream.h #includestdlib.h #includemath.h #includefloat.h #includeiomanip.h #includetime.h #define E 1.0e-12 /*定义全局变量相对误差限*/ #define L 20 /*迭代次数*/ FILE *fp; struct C /*定义结构体*/ { double R; double I; }; void shuchu(double a[][10]) /*输出一个10*10的矩阵*/ { int i,j; double b=0; for(i=0;i10;i++) { for(j=0;j10;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 sgn(double x) /*符号函数*/ { if(x0) return -1; if(x==0) return 0; if(x0) return 1; } void assignment(double a[10][10]) /*输入矩阵A*/ { int i,j; for(i=0;i=9;i++) { for(j=0;j=9;j++) { if(i==j) a[i][j]=1.5*cos((i+1)+1.2*(j+1)); else a[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } } } void Hessen(double a[10][10]) /*将一个10*10的矩阵拟上三角化*/ {