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

数值分析第二次大作业姓名:袁二凯 学号:SY0917145 联系方式:13488854452 6421011362011年11月8日题目:使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： (i,j=1,2,10)一、 算法的设计方案：（一）、总体方案设计：构造矩阵 A，先利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A(n-1)，然后进行带双步位移的QR分解求解矩阵的全部特征值，然后对实特征值利用列主元高斯消元法求解其对应的特征向量。题目要求输出矩阵A(n-1)QR分解得到的矩阵Q、R以及乘积矩阵RQ，可用单独的QR分解法进行计算输出。（二）具体算法设计：1、对矩阵赋值：即为下述程序中的double fuzhi（）子程序。2、对上述生成的矩阵进行拟上三角化：即为下述程序中的void NSSJH（）子程序，输出矩阵A(n-1)。3、对拟上三角化后的矩阵进行QR分解：程序中用void QR（）子程序来对拟上三角化过后的A阵进行QR分解，并输出Q阵、R阵和乘积矩阵RQ阵。4、对拟上三角化后的矩阵进行带双步位移的QR分解。子程序void doubleQR（）实现对拟上三角化后的A阵进行带双步位移的QR分解，得出特征值并输出。5、求解每一个实特征值相应的特征向量：用子程序void TZXL（）对其中的实数特征值进行求解，得出对应的特征向量并输出。二、源程序：#include#include#include #include #include#include /*头文件*/ /*定义全局变量*/ #define Epsilon 1e-12 /*定义精度*/ #define N 11 /*取N为11，可实现从1到10的存储，方便计算*/ #define n 10 /*为了方便计算取n*/ #define L 1000 /*规定迭代次数*/ double aNN; double lambdaN2;double fuzhi(); /*赋值函数*/void NSSJH(); /*拟上三角化函数*/void QR(); /*QR分解函数*/void DoubleQR();/*带双步位移的QR分解函数*/void TZXL(); /*求取实特征根对应特征向量的函数*/double fuzhi() /*赋值函数，存储行列均为1-10*/int i,j; for(i=1;iN;i+) for(j=1;jN;j+) if(i=j) aij = 1.5*cos(i+1.2*j); else aij = sin(0.5*i+0.2*j); return aij;void NSSJH() /*拟上三角化函数*/ int i,j,r; double cr,dr,hr,tr,temp;double prN=0,qrN=0,urN,wrN=0;fuzhi(); for(r=1;r=n-2;r+) dr=0;cr=0;hr=0;tr=0;temp=0;for(i=r+1;i0)cr=-dr;else cr=dr;/*求解cr*/hr=cr*cr-cr*ar+1r; /*求解hr*/for(i=1;i=r;i+)uri=0;urr+1=ar+1r-cr;for(i=r+2;i=n;i+)uri=air; /*生成urN*/ for(i=1;i=n;i+) temp=0;for(pri=0.0,j=1;j=n;j+)pri=pri+aji*urj/hr; /*求解prN*/ for(i=1;i=n;i+) temp=0;for(qri=0.0,j=1;j=n;j+)qri=qri+aij*urj/hr;/*求解qrN*/ for(i=1;i=n;i+) tr=tr+uri*pri/hr;/*求解tr*/ for(i=1;i=n;i+) wri=qri-uri*tr;/*求解wrN*/ for(i=1;i=n;i+)for(j=1;j=n;j+)aij=aij-wri*urj-uri*prj; /*求解A(r+1)*/printf(%对A进行拟上三角化得到的矩阵A(n-1):n); for(i=1;i=n;i+)for(j=1;j=n;j+)printf(%13.11e ,aij);printf(n,aij);/*输出拟上三角化后的A（n-1）*/ void QR() /*QR分解函数*/int i,j,k,r; double dr,cr,hr,temp,sum; double QNN=0.0,wrN=0.0,urN=0.0,prN=0.0,ANN=0.0,bN=0.0; for(i=1;i=n;i+) for(j=1;j=n;j+) if(i=j) Qij=1; else Qij=0; /*生成QN*/ NSSJH(); for(r=1;rn;r+) sum=0; for(i=r+1;i0) cr=-dr; else cr=dr; /*求解cr*/ hr=cr*cr-cr*arr; /*求解hr*/ for(i=1;ir;i+) uri=0; urr=arr-cr; for(i=r+1;i=n;i+) uri=air; /*生成urN*/ for(i=1;i=n;i+) temp=0; for(j=1;j=n;j+)temp=temp+Qij*urj;wri=temp; /*求解wrN*/ for(i=1;i=n;i+) for(j=1;j=n;j+) Qij=Qij-wri*urj/hr; /*求解QN，即为所求矩阵Q*/ for(i=1;i=n;i+) temp=0; for(j=1;j=n;j+) temp=temp+aji*urj/hr; pri=temp; /*求解prN*/ for(i=1;i=n;i+) for(j=1;j=n;j+) aij=aij-uri*prj; /*求解aNN，即为所求R矩阵*/ for(i=1;i=n;i+) for(j=1;j=n;j+) temp=0; for(k=1;k=n;k+) temp+=aik*Qkj; Aij=temp; /*求解AN，即为所求乘积矩阵RQ*/ printf(%对A进行QR得到的矩阵Q:n); /*输出Q、R、RQ*/ for(i=1;i=n;i+)for(j=1;j=n;j+)printf(%13.11e ,Qij);printf(n); printf(%对A进行QR得到的矩阵R:n); for(i=1;i=n;i+)for(j=1;j=n;j+)printf(%13.11e ,aij);printf(n); printf(%对A进行QR得到的乘积矩阵RQ:n); for(i=1;i=n;i+)for(j=1;j=n;j+)printf(%13.11e ,Aij);printf(n); void DoubleQR() /*带双步位移的QR分解函数*/int i,j,k,m,r,l;double b,c,d,s,t,MNN=0.0,INN=0.0; double cr,dr,hr,tr,temp,sum; double prN=0.0,qrN=0.0,urN=0.0,wrN=0.0,vrN=0.0; k=1;m=10;NSSJH(); /*先进行拟上三角化*/ for(i=1;i=m;i+) for(j=1;j=m;j+) if(i=j) Iij=1; else Iij=0; /*定义单位矩阵*/loop1: /*降一阶的判断*/ if(fabs(amm-1)=Epsilon) lambdam0=amm; lambdam1=0; m=m-1; goto loop2; else goto loop3; loop2: /*判断m=0) lambdam0=(c+sqrt(d)/2; lambdam1=0;lambdam-10=(c-sqrt(d)/2; lambdam-11=0; else lambdam0=c/2; lambdam1=sqrt(-d)/2;lambdam-10=c/2; lambdam-11=-sqrt(-d)/2; goto loopjs; else if(m=1) lambda10=a11; lambda11=0; goto loopjs; else if(m=0) goto loopjs; else goto loop1; loop3: /*降两阶的判断*/ if(fabs(am-1m-2)=0) lambdam0=(c+sqrt(d)/2; lambdam1=0;lambdam-10=(c-sqrt(d)/2; lambdam-11=0; else lambdam0=c/2; lambdam1=sqrt(-d)/2;lambdam-10=c/2; lambdam-11=-sqrt(-d)/2; m=m-2; goto loop2; else goto loopPD; loopPD: /*判断迭代次数是否超出规定*/ if(k=L) printf(计算终止，未得到A的全部特征值n); goto loopjs; else goto loopqr; loopqr:s=am-1m-1+amm; /*带双步位移的QR分解*/ t=am-1m-1*amm-amm-1*am-1m; for(i=1;i=m;i+) for(j=1;j=m;j+) temp=0; for(l=1;l=m;l+) temp=temp+ail*alj; Mij=temp-s*aij+t*Iij; /*生成MNN*/ for(r=1;rm;r+) sum=0; for(i=r+1;i0) cr=-dr; else cr=dr;/*求解cr*/ hr=cr*cr-cr*Mrr;/*求解hr*/ for(i=1;ir;i+) uri=0.0; urr=Mrr-cr; for(i=r+1;i=m;i+) uri=Mir; /*生成urN*/ for(j=1;j=m;j+) temp=0;for(i=1;i=m;i+) temp+=Mij*uri/hr;vrj=temp;/*求解vrN*/ for(i=1;i=m;i+) for(j=1;j=m;j+) Mij=Mij-uri*vrj; /*求解MN（k+1）*/ for(j=1;j=m;j+) temp=0; for(i=1;i=m;i+) temp+=aij*uri/hr; prj=temp;/*求解prN*/ for(i=1;i=m;i+) temp=0;for(j=1;j=m;j+) temp+=aij*urj/hr; qri=temp;/*求解qrN*/ for(tr=0.0,i=1;i=m;i+) tr+=pri*uri/hr;/*求解tr*/ for(i=1;i=m;i+) wri=qri-tr*uri;/*求解wrN*/ for(i=1;i=m;i+) for(j=1;j=m;j+) aij=aij-wri*urj-uri*prj; /*求解A(k+1)*/ goto loopk; loopk:k=k+1; goto loop1;loopjs: ; printf(%矩阵A的全部特征值:n); for(i=1;i=n;i+) printf(%13.11e+%13.11ein,lambdai0,lambdai1); /*输出求得的特征值*/ void TZXL() /*求解实特征根对应的特征向量，列主元素Gauss消去法*/int i,j,k,l,ik;double xN=0.0,INN=0.0,m,temp;for(i=1;i=n;i+) for(j=1;j=n;j+) if(i=j) Iij=1; else Iij=0; /*定义单位矩阵*/DoubleQR();/*求出特征值*/for(l=1;l=n;l+)if(lambdal1!=0)continue;/*如果不是实根，则l+，继续循环*/elsefor(i=1;in;i+)for(j=1;jn;j+)aij=aij-lambdal0*Iij; /*构造矩阵（A-*I）*/for(k=1;kn;k+)ik=k;for(i=k;i=n;i+)if(aikkaik)ik=i; /*选主元*/for(j=k;j=n;j+) temp=0;temp=aikj;aikj=akj;akj=temp; for(i=k+1;i=n;i+) m=aik/akk; for(j=k+1;j0;k-) temp=0; for(j=k+1;j=n;j+) temp=temp+akj*xj/akk; xk=xk-temp; /*回代过程，xN即为特征值对应的特征向量*/ printf(特征值%13.11e对应的特征向量为：n,lambdal0); for(i=1;i=n;i+) printf(%13.11e ,xi); printf(n); /*输出特征向量*/ void main()void NSSJH();void QR();DoubleQR();TZXL();二、运行结果输出对A进行拟上三角化得到的矩阵A(n-1):-8.82751675883e-001 -9.93313649183e-002 -1.10334928599e+000 -7.60044358564e-001 1.54910107991e-001 -1.94659186287e+000 -8.78243638293e-002 -9.25588938718e-001 6.03259944053e-001 1.51886095647e-001-2.34787836242e+000 2.37237010494e+000 1.81929082221e+000 3.23780410155e-0012.20579844032e-001 2.10269266255e+000 1.81613808610e-001 1.27883908999e+000-6.38057812440e-001 -4.15407560380e-0010.00000000000e+000 1.72827459997e+000 -1.17146764279e+000 -1.24383926270e+000-6.39975834174e-001 -2.00283307904e+000 2.92494720612e-001 -6.41283006839e-001 9.78399762128e-002 2.55776357416e-0010.00000000000e+000 -5.40440307281e-018 -1.29166953413e+000 -1.11160351340e+000 1.17134682410e+000 -1.30735603002e+000 1.80369917775e-001 -4.24638535837e-001 7.98895523930e-002 1.60881992807e-0010.00000000000e+000 -7.99016164618e-017 8.29733316381e-017 1.56012629853e+0008.12504939751e-001 4.42175683292e-001 -3.58861612814e-002 4.69174231367e-001-2.73659505009e-001 -7.35933465775e-0020.00000000000e+000 1.28519026519e-018 -9.51733763689e-017 0.00000000000e+000-7.70777375519e-001 -1.58305142574e+000 -3.04284317680e-001 2.52871244603e-001 -6.70992540145e-001 2.54461992908e-0010.00000000000e+000 1.83527558474e-016 1.57658260703e-017 0.00000000000e+0000.00000000000e+000 -7.46345345694e-001 -2.70836515702e-002 -9.48652189368e-001 1.19587108150e-001 1.92926561795e-0020.00000000000e+000 -4.31880044039e-017 1.12848923378e-016 0.00000000000e+0000.00000000000e+000 0.00000000000e+000 -7.70180137436e-001 -4.69762399062e-001 4.98825946801e-001 1.13769160378e-0010.00000000000e+000 7.14764653623e-017 -8.54810208228e-017 0.00000000000e+0000.00000000000e+000 0.00000000000e+000 0.00000000000e+000 7.01316709211e-0011.58218068848e-001 3.86259461423e-0010.00000000000e+000 7.92254903113e-017 2.72921948474e-017 0.00000000000e+0000.00000000000e+000 0.00000000000e+000 0.00000000000e+000 0.00000000000e+0004.84380760278e-001 3.99277799518e-001对A进行QR得到的矩阵Q:-3.51926257953e-001 4.42759198224e-001 -6.95598251361e-001 6.48620075365e-0023.70971886190e-001 1.85584714361e-001 -1.62894231963e-002 -1.18105316965e-001 -5.25537538372e-002 -5.48658294357e-002-9.36027728736e-001 -1.66467918654e-001 2.61529954856e-001 -2.43867172893e-002 -1.39477436089e-001 -6.97758539124e-002 6.12447214296e-003 4.44050544314e-002 1.97590790973e-002 2.06283697053e-0020.00000000000e+000 -8.81052055469e-001 -3.98976279696e-001 3.72030872848e-0022.12779406409e-001 1.06446355722e-001 -9.34317107976e-003 -6.77420046453e-002 -3.01434069867e-002 -3.14695508044e-0020.00000000000e+000 2.75509484197e-018 -5.37180680644e-001 -1.23494585420e-001-7.06315160872e-001 -3.53345636850e-001 3.10143894826e-002 2.24867649160e-001 1.00060178353e-001 1.04462274870e-0010.00000000000e+000 4.07328114527e-017 3.04286786187e-017 9.89223546862e-001-1.23941173121e-001 -6.20035858983e-002 5.44227283946e-003 3.94588163723e-0021.75581335001e-002 1.83305946291e-0020.00000000000e+000 -6.55173387862e-019 -3.95151900673e-017 4.30563223546e-0175.32361069026e-001 -6.73390034490e-001 5.91058120587e-002 4.28542532387e-001 1.90690134319e-001 1.99079449530e-0010.00000000000e+000 -9.35599774666e-017 1.59243342181e-017 2.46175662403e-017-3.22773883960e-017 -6.05976150575e-001 -9.16578303282e-002 -6.64558650897e-001 -2.95711087758e-001 -3.08720746256e-0010.00000000000e+000 2.20166864990e-017 4.47273055557e-017 -5.82522024871e-017-2.97858484095e-017 -1.51379970365e-016 9.93339662512e-001 -9.69044031194e-002 -4.31199058447e-002 -4.50169441118e-0020.00000000000e+000 -3.64377783058e-017 -3.19016158777e-017 5.08230496082e-0171.68129388097e-017 1.09570105276e-016 3.27463232780e-017 5.41008800606e-001-5.81754083823e-001 -6.07348058054e-0010.00000000000e+000 -4.03881310792e-017 1.53941362314e-017 1.40674401021e-018-2.05037071608e-017 -4.84030334626e-017 1.21453791695e-017 -1.34808250804e-017 -7.22159133673e-001 6.91726958888e-001对A进行QR得到的矩阵R:2.50834274492e+000 -2.18564688549e+000 -1.31460907079e+000 -3.55878749384e-002 -2.60985785039e-001 -1.28312184709e+000 -1.39087861061e-001 -8.71289797216e-001 3.84936790297e-001 3.35380289967e-0010.00000000000e+000 -1.96160327785e+000 2.40752372763e-001 7.05471457282e-0015.95720431828e-001 5.52697877468e-001 -3.26820992441e-001 -5.76949866836e-002 2.87112933019e-001 -8.89512875419e-0020.00000000000e+000 0.00000000000e+000 2.40453460199e+000 1.70675809633e+000-4.23956670409e-001 3.40533230581e+000 -1.05001765585e-001 1.46225710273e+000-6.68448746928e-001 -4.02764620966e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+000 1.57712208072e+0006.39953513396e-001 3.46812787243e-001 -5.70178664977e-002 4.01478805443e-001-2.22247617631e-001 -6.31705923644e-0020.00000000000e+000 0.00000000000e+000 0.00000000000e+000 0.00000000000e+000-1.44784699777e+000 -1.41572400774e+000 -2.80613904467e-001 -2.81791052189e-001 -4.61143488185e-002 1.99662907996e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+000 0.00000000000e+0001.76632666819e-016 1.23164145154e+000 1.61970100342e-001 1.96263827550e-0015.35003562176e-001 -1.50927342477e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+000 0.00000000000e+000-1.55036704937e-017 0.00000000000e+000 -7.75344191421e-001 -3.4645145088