北航数值分析第二次大作业概述.docx

1、北航数值分析第二次大作业概述数值分析第二次大作业姓名：李潇 学号：SY1303314 题目:使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： (i,j=1,2,10)1、算法设计：矩阵的拟上三角化：对实矩阵A进行相似变换化为拟上三角矩阵，其变换矩阵采用Householder矩阵，变换过程如下：若，则；否则，。当时，得，令又是对称正交矩阵，于是成立，因而与相似。矩阵的QR分解：矩阵的QR分解过程与拟上三角化过程相似，在这里不再重复其原理。 求全部特征值矩阵拟上三角化后利用带双步位移的QR方法，采用书本Page 63页具体算法实现。为了使编程方便，并体

2、会goto语句使用的灵活性，程序中的跳转均使用goto Loop*实现。求A的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组，。因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n次（n为矩阵阶数）。线性方程组的求解采用列主元素Gauss消去法。2、程序源代码：#include Stdio.h#include Conio.h#include math.h/*/ N ：待求矩阵的阶数 / Epsilon：精度水平 / MaxL ：限制求解的最大迭代次数 /*/#define N 10#define Epsilon 0.000000000001#define M

3、axL 200void Creat_MatrixA(double arrayN+1N+1);void Load_MatrixA(double MatrixAN+1N+1,double arrayAN+1N+1);void Clear_Elements(double arrayN+1N+1);void Initial_to_I(double arrayN+1N+1);void Nishang(double arrayAN+1N+1);void QR_Fenjie(double arrayAN+1N+1,double QN+1N+1,int m);void Solution_namda(doubl

4、e arrayAN+1N+1,double namdaN+12);void Solution_A_Next(double arrayAN+1N+1,int m,double s,double t);void Solution_vector(double arrayN+1N+1,double namdaN+12,double xN+1N+1);void GaussElemitation_select(int m,int flag,double arrayN+1N+1,double xN+1N+1);void Print_test(double arrayN+1N+1);void Print_te

5、stAB(double arrayAN+1N+1,double arrayBN+1N+1);void Print_namda(double namdaN+12);void Print_vector(double xN+1N+1,double namdaN+12);int main(void)/*/ 定义主程序中用到的变量 /*/ MatrixA：用来存储源矩阵A / Q ：用来存储对A(n-1)进行QR分解得到的正交阵 / namda ：用来存储矩阵A的全部特征值。第1列存储实部，第2列存储虚部 /*/ double MatrixAN+1N+1,QN+1N+1,namdaN+12,xN+1N+

6、1;/*/ 对矩阵A进行拟上三角化，再进行QR分解见 /*/ Creat_MatrixA(MatrixA); printf(A经过拟上三角化得到的矩阵A(n-1)为：n);/*对矩阵A进行拟上三角化，并输出矩阵A(n-1)*/ Nishang(MatrixA); Print_test(MatrixA); printf(nA(n-1)经过QR分解得到的矩阵Q、R和RQ分别为：n);/*对矩阵A(n-1)进行QR分解*/ QR_Fenjie(MatrixA,Q,N);/*输出分解得到的Q、R和RQ*/ printf(Q:n); Print_test(Q); printf(R:n); Print_t

7、est(MatrixA); printf(RQ:n); Print_testAB(MatrixA,Q); /*/ 用带双步位移的QR方法求矩阵A的全部特征值 /*/*重新生成矩阵A*/ Creat_MatrixA(MatrixA);/*调用Solution_namda求矩阵A的全部特征值，并输出*/ Solution_namda(MatrixA,namda); Print_namda(namda);/*/ 求实特征值对应的特征向量 /*/*重新生成矩阵A*/ Creat_MatrixA(MatrixA); printf(n实特征值对应的特征向量：);/*调用Solution_vector求矩阵

8、A实特征值对应的特征向量，并输出*/ Solution_vector(MatrixA,namda,x); Print_vector(x,namda); getch(); return 0;/*创建源矩阵A*/void Creat_MatrixA(double arrayN+1N+1) int i,j;/*矩阵第0行和第0列赋为0*/ for(i=0;i=N;i+) arrayi0=0; for(j=0;j=N;j+) array0j=0; for(i=1;i=N;i+) for(j=1;j=N;j+) if(i=j) arrayij=1.5*cos(i+1.2*j); else arrayij

9、=sin(0.5*i+0.2*j); /*做矩阵A的拷贝*/* 求某个特征值对应的特征向量，实际上是求解一个n元线性方程组的根， 因此每求一个特征向量都要重新载入矩阵A。*/void Load_MatrixA(double MatrixAN+1N+1,double arrayAN+1N+1) int i,j; for(i=1;i=N;i+) for(j=1;j=N;j+) arrayAij=MatrixAij;/*输出给定矩阵的所有元素，按2个元素一行输出*/void Print_test(double arrayN+1N+1) int i,j; for(i=1;i=N;i+) for(j=1

10、;j=N;j+) printf(a%d%d= %.12e ,i,j,arrayij); /*2个元素一行输出*/ if(j%2=0) printf(n); /*输出给定矩阵A、B的积AB的所有元素，按2个元素一行输出*/void Print_testAB(double arrayAN+1N+1,double arrayBN+1N+1) int i,j,k; double temp; for(i=1;i=N;i+) for(j=1;j=N;j+) temp=0; for(k=1;k=N;k+) temp+=arrayAik*arrayBjk; printf(%.12e ,temp); /*2个元

11、素一行输出*/ if(j%2=0) printf(n); /*输出求得的所有特征值*/void Print_namda(double namdaN+12) int i; for(i=1;i=N;i+) if(namdai1!=0) printf(%d)= %.12e+%.12ein,i,namdai0,namdai1); else printf(%d)= %.12en,i,namdai0); /*输出矩阵A的实特征值对应的特征向量，按3个元素一行输出*/void Print_vector(double xN+1N+1,double namdaN+12) int i,j; for(i=1;i=N

12、;i+) if(namdai1=0) printf(n(%d)=%.12e对应的特征向量为：n,i,namdai0); for(j=1;j=N;j+) printf(%.12e ,xij); /*3个元素一行输出*/ if(j%3=0) printf(n); else ; /*将低于精度水平的矩阵元素清0*/* 拟上三角化过程后矩阵A的下三角元素的精确值应该为0，但由于计算机中存储的数与精确值都有一定误差， 因此，实际上得到的下三角元素不为0，而是一个绝对值很小的数。为了避免这些数对后面的计算造成影响， 故要对拟上三角化后的矩阵A下三角元素进行清0。*/void Clear_Elements(

13、double arrayN+1N+1) int i,j; for(i=1;i=N;i+) for(j=1;j=i-1;j+) if(fabs(arrayij)Epsilon) arrayij=0; /*将矩阵初Q始化为单位矩阵I*/void Initial_to_I(double arrayN+1N+1) int i,j; for(i=0;i=N;i+) for(j=0;j=N;j+) arrayij=0; for(i=1;i=N;i+) arrayii=1;/*矩阵A的拟上三角化*/void Nishang(double arrayAN+1N+1) int r,i,j,flag; doubl

14、e d_r,c_r,h_r,t_r,temp; double w_rN+1,p_rN+1,q_rN+1,uN+1; for(r=1;r=N-2;r+) /*判断(i=r+1,，n)，Air是否全为0*/ flag=0; for(i=r+2;i=N;i+) if(arrayAir!=0) flag=1; break; else ; /*Air不全为0时，由A(r)求A(r+1)*/ if(flag=1) /*计算d_r、c_r和h_r*/ d_r=0; for(i=r+1;i0)?(-d_r):(d_r); h_r=c_r*(c_r-arrayAr+1r); /*给u_r各元素赋值*/ for(

15、i=1;i=N;i+) if(i=r) ui=0; else if(i=r+1) ui=arrayAr+1r-c_r; else ui=arrayAir; /*计算p_r、q_r各元素*/ for(i=1;i=N;i+) temp=0; for(j=1;j=N;j+) temp+=arrayAji*uj; p_ri=temp/h_r; for(i=1;i=N;i+) temp=0; for(j=1;j=N;j+) temp+=arrayAij*uj; q_ri=temp/h_r; /*计算t_r的值和w_r各元素*/ temp=0; for(i=1;i=N;i+) temp+=p_ri*ui;

16、 t_r=temp/h_r; for(i=1;i=N;i+) w_ri=q_ri-t_r*ui; /*由A(r)计算A(r+1)*/ for(i=1;i=N;i+) for(j=1;j=N;j+) arrayAij=arrayAij-w_ri*uj-ui*p_rj; /endif /*Air全为0时，A(r+1)=A(r)*/ else ; /*调用Clear_Elements将下三角清0*/ Clear_Elements(arrayA);/*矩阵A的QR分解*/void QR_Fenjie(double arrayAN+1N+1,double QN+1N+1,int m) int r,i,j

17、,flag; double d_r,c_r,h_r,temp; double w_rN+1,p_rN+1,uN+1; Initial_to_I(Q); for(r=1;r=m-1;r+) flag=0; for(i=r+1;i=m;i+) if(arrayAir!=0) flag=1; break; else ; if(flag=1) d_r=0; for(i=r;i0)?(-d_r):(d_r); h_r=c_r*(c_r-arrayArr); for(i=1;i=m;i+) if(i=r-1) ui=0; else if(i=r) ui=arrayArr-c_r; else ui=arra

18、yAir; for(i=1;i=m;i+) temp=0; for(j=1;j=m;j+) temp+=Qij*uj; w_ri=temp; for(i=1;i=m;i+) for(j=1;j=m;j+) Qij=Qij-w_ri*uj/h_r; for(i=1;i=m;i+) temp=0; for(j=1;j=m;j+) temp+=arrayAji*uj; p_ri=temp/h_r; for(i=1;i=m;i+) for(j=1;j=m;j+) arrayAij=arrayAij-ui*p_rj; /endif else ; /*/ 带双步位移的QR分解求矩阵的全部特征值 /*/vo

19、id Solution_namda(double arrayAN+1N+1,double namdaN+12) int k,m; double det,s,t,namda12,namda22;/*初始化存放特征值的数组元素*/ for(k=0;k=N;k+) for(m=0;m=1;m+) namdakm=0;/*Loop1：对矩阵A进行拟上三角化*/ Nishang(arrayA);/*Loop2：迭代次数初始化为1，矩阵阶数初始化为N*/ k=1; m=N;Loop3: if(fabs(arrayAmm-1)1) goto Loop3; /*转3*/ else if(m=1) namdam

20、0=arrayAmm; goto Loop11; /*转11结束*/ else if(m=0) goto Loop11; /*转11结束*/Loop5: /*求二阶子阵的两个特征值*/ s=arrayAm-1m-1+arrayAmm; t=arrayAm-1m-1*arrayAmm-arrayAm-1m*arrayAmm-1; det=pow(s,2)-4*t; if(det=0) namda10=(s+sqrt(det)/2; namda20=(s-sqrt(det)/2; else if(det=0) namdam0 =namda10; namdam-10=namda20; m=m-2;

21、else namdam0 =namda10; namdam1 =namda11; namdam-10=namda20; namdam-11=namda21; m=m-2; goto Loop11; else goto Loop7;Loop7: if(fabs(arrayAm-1m-2)=0) namdam0 =namda10; namdam-10=namda20; m=m-2; else namdam0 =namda10; namdam1 =namda11; namdam-10=namda20; namdam-11=namda21; m=m-2; goto Loop4; else goto L

22、oop8;Loop8: if(k=MaxL) goto Loop12; /*计算终止，未得到全部特征值*/ else goto Loop9;Loop9: Solution_A_Next(arrayA,m,s,t); goto Loop10;Loop10: k+; goto Loop3;Loop11: printf(n求解成功: 总共迭代k=%d次n,k); goto Loop13;Loop12: printf(迭代次数超出限制，算法失败！); goto Loop13;Loop13: ;/*由A(k)计算A(k+1)*/void Solution_A_Next(double arrayAN+1N+1,int m,double s,double t) int i,j,k; double QN+1N+1,MN+1N+1; double temp;/*计算矩阵M_k：M_k=AA-SA+TI*/ for(i=1;i=m;i+) for(j=1;j=m;j+)