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

数值分析第二次 SY1107214 王宽一、 题目分析使用带双步位移的QR分解法求矩阵的全部特使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： (i,j=1,2,10)二、 算法设计1. 按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为initA()；2. 为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化：对应的函数为nssj()； 3. 对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为tezhengzhi()；这个过程中包含着两个子程序：QR()和qmk()，分别用来对矩阵Mk进行QR分解并得到Ak+1和计算mk的值；4. 使用带原点平移的反幂法求出其对应的特征向量，对应的函数为：fmifa()，这个过程中求解线性方程时用到列主元的高斯消元法，对应的函数为：gauss(double );根据数值分析课本的相关知识，步骤3中带双步位移的QR分解法的流程图如下：本次作业所使用的编译环境为：Visual c+6.0三、 程序源代码#include#includedouble A1010,rr1010,qq1010,rq1010,uk1010;double xy10,x10=0;/反幂法中double tzr10=0,tzi10=0;/定义矩阵的特征值数组,r实部、i虚部/对矩阵A进行初始化，赋值void initA() int i,j; for(i = 1; i = 10; i+) for(j= 1; j= 10;j+) Ai-1j-1=sin(0.5*i+0.2*j); for(i= 1; i= 10;i+) Ai- 1i- 1=cos(i+1.2*i)*1.5;/定义函数对矩阵A进行拟上三角化void nssj() int r,i,j; int sgn(double a); double d,c,h,t,tr,u10,p10,q10,w10; double sum=0,tmp=0; /开始循环 for (r=0;r8;r+) t= 0;for(i= r+2;i10;i+) t= t+(Air=0);if (t=8- r) continue; else sum=0; for (i=r+1;i10;i+) sum=sum+ Air*Air; d= sqrt(sum); c= -1*sgn(Ar+1r)*d; h= c*c- c*Ar+1r; /step (3) for (i=0;i10;i+)ui= 0; for (i=r+2;i10;i+)ui= Air; ur+1= Ar+1r- c; /step (4) for(i= 0;i 10; i+) tmp= 0; sum= 0; for(j= 0;j 10; j+) tmp= tmp+Aji* uj; sum= sum+Aij* uj; pi= tmp/h; qi= sum/h; sum=0; for(i= 0;i 10; i+) sum= sum+ pi*ui; tr= sum/h; for(i= 0;i 10; i+) wi= qi- tr*ui; for(i= 0;i 10;i+) for(j= 0;j 10;j+) Aij= Aij- wi*uj- ui*pj; /以上的A存在一定的误差，消除误差 for(i=0; i10; i+) for(j= 0; j 10; j+) if (-1.0e-12Aij & Aij= 0) return 1; else return -1;/计算A的特征值的函数（课本P63 11步）int tezhengzhi() int k=1,m=9,n=9; int flag=0; double ms ,mt,mk1010;/step5 double detdk,fb,temp,s1_r,s1_i,s2_r,s2_i; void qmk(double mk1010,int m,double ms,double mt);/声明求mk 的函数 void QR (double mk1010,int m);/声明对mk进行qr分解的函数 while (1) if (-1e-12 Amm-1 & Amm-11e-12) tzrn= Amm; n-;m-; flag=1; /step 4 if(flag) flag= 0; if (m=0) tzrn= Amm; return 0; else if(m=-1) return 0; else continue; /step 5 else detdk = Am-1m-1 * Amm - Am-1m * Amm-1; fb= -1.0 * (Am - 1m - 1 + Amm); temp = fb * fb - 4.0 * detdk; if(temp 0) temp = -temp; s1_r = -0.5 *fb ; s1_i = 0.5* sqrt(temp); s2_r = -0.5 *fb; s2_i = -0.5 * sqrt(temp); else s1_r = (-1.0 *fb + sqrt(temp)/ 2; s1_i = 0; s2_r = (-1.0 *fb - sqrt(temp)/ 2; s2_i = 0; /step 6if(m=1) tzr1= s2_r; tzi1= s2_i; tzr0= s1_r; tzi0= s1_i; return 0;else if (-1e-12Am-1m-2 & Am-1m-21e-12) tzrn= s2_r; tzin= s2_i; n-;m-; tzrn= s1_r; tzin= s1_i; n-;m-;if (m=0)tzrn= Amm; return 0;else if(m=-1) return 0;else continue; else if (k= 500) return 1; elsems = Am-1m-1+ Amm; mt = Am-1m-1* Amm-Am-1m*Amm-1; qmk(mk,m,ms,mt);QR(mk,m);k+;/计算mk的子程序； void qmk(double mk1010, int m, double ms, double mt) int i,j,k; double tmp_sum; for(i = 0; i = m; i+) for(j = 0; j = m; j+) tmp_sum = 0; for(k = 0; k = m; k+) tmp_sum+= Aik * Akj; mkij = tmp_sum - ms * Aij; for(i = 0; i = m; i+) mkii+= mt;/对mk做QR分解并得到A（k+1）的函数void QR(double mk1010,int m)inttmp_ir;intr, i, j;double ur10, vr10, pr10, qr10, wr10, tmp_dr, tmp_vr, tmp_pr, tmp_qr, tmp_tr, dr, cr, hr, tr;for(r=0; rm;r+) /sp 1 tmp_ir = 0; for(i=r+1;i=m;i+) tmp_ir=tmp_ir+(mkir=0); if (tmp_ir=m-r) continue; else /sp 2 tmp_dr = 0; for(i = r; i = m; i+) tmp_dr = tmp_dr + mkir * mkir; dr = sqrt(tmp_dr); cr = -1.0 * sgn(mkrr) * dr; hr = cr * cr - cr * mkrr;/sp 3 for(i = 0; i r; i+) uri = 0; for(i = r + 1; i = m; i+) uri = mkir; urr = mkrr-cr;/sp 4 for(i=0; i=m;i+) tmp_vr = 0; tmp_pr = 0; tmp_qr = 0; for(j = 0; j = m; j+) tmp_vr = tmp_vr + mkji * urj; tmp_pr = tmp_pr + Aji * urj; tmp_qr = tmp_qr + Aij * urj; vri = tmp_vr / hr; pri = tmp_pr / hr; qri = tmp_qr / hr; tmp_tr = 0; for(i = 0; i = m; i+) tmp_tr = tmp_tr + pri * uri; tr = tmp_tr / hr; for(i = 0; i = m; i+) wri = qri - tr * uri; for(i = 0; i = m; i+) for(j = 0; j v1)?v1:v2); double max(double v1,double v2,double v3)double temp;temp= (v1v2)?v1:v2); return (tempv3)?temp:v3);/构造函数translation 用于进行原点平移void translation (double x) int i=0,j=0; for(i=0;i10;i+)for(j=0;j10;j+)Aji= Aji-x;/列主元的高斯消元法求解特征向量void gauss(double b10) int i,j,k,mk; double mm,f,mik; for(k = 0;k9;k+) /*选主元并消元*/ mm=Akk;/记录最大组员mk = k;/记录最大组员的行号for(i=k;i10;i+)/*选择第K列主元素 if(fabs(mm)fabs(Aik) mm = Aik; mk = i; if(mk!=k)/* 将第K列主元素换行到对角线上*/ for(j=k;j10;j+) f = Akj; Akj=Amkj; Amkj=f; f = bk; bk=bmk; bmk=f; for(i=k+1;i10;i+) /*将第K列对角线以下消元为零*/ mik = Aik/Akk; for(j=k+1;j=0;k-) f= 0;for (j=k+1;j10;j+)f= f+ Akj*xj;xk= (bk-f)/Akk; /构造反幂法函数void fmifa() intm,k,i; doubleu10,y10; double beita1=1,beita2=0,sum=0, h=0; for ( i=0;i10;i+) ui=1; do if(k!=0) beita1=beita2; for (m=0;m10;m+) sum+=um*um; for (m=0;m10;m+) ym= um/sqrt (sum); sum= 0 ; beita2=0; gauss ( y); for(m=0;m= fabs(beita2)*exp(-12); for (m=0;m10;m+) xym=ym; / 主程序void main()int i,j;initA();nssj() ;printf(经过上三角化后的矩阵为：nn );for (i=0;i10;i+) for (j=0;j10;j+) printf(%.12e ,Aij); printf(nn); /求矩阵特征值tezhengzhi();printf(矩阵A各个特征值的实部及虚部为：nn);for (i=0;i10;i+) printf(第%d个特征值为：（%.12e，%.12ei）nn,i+1,tzri,tzii); printf(nn );printf(经带双步位移的QR分解后的矩阵为：n );for (i=0;i10;i+) for (j=0;j10;j+) if (-1.0e-12Aij & Aij 1.0e-12) Aij=0; printf(%.12e ,Aij); printf(nnn); initA();printf(特征向量如下：nn);for (i=0;i10;i+)if (tzii != 0) continue ;else initA() ;translation (tzri) ;fmifa();printf(特征值 %.12e对应的特征向量为： （,tzri);for (j=0;j9;j+) if (-1.0e-12xyj & xyj 1.0e-12) xyj=0; printf(%.12e,xyj+tzri);if (j= 9) if (-1.0e-12xyj & xyj 1.0e-12) xyj=0; printf(%.12e,xyj+tzr9);printf(）nn );四、 输出结果1. 经过上三角化后的矩阵即:-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000 -7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000 -8.782436382927e-002 -9.255889387184e-001 6.032599440534e-001 1.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+000 3.237804101546e-001 2.205798440320e-001 2.102692662546e+000 1.816138086098e-001 1.278839089990e+000 -6.380578124405e-001 -4.154075603804e-0010.000000000000e+000 1.728274599968e+000 -1.171467642785e+000 -1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+000 2.924947206124e-001 -6.412830068395e-001 9.783997621285e-002 2.557763574160e-0010.000000000000e+000 0.000000000000e+000 -1.291669534130e+000 -1.111603513396e+000 1.171346824096e+000 -1.307356030021e+000 1.803699177750e-001 -4.246385358369e-001 7.988955239304e-002 1.608819928069e-0010.000000000000e+0000.000000000000e+000 0.000000000000e+000 1.560126298527e+000 8.125049397515e-001 4.421756832923e-001 -3.588616128137e-0024.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000 -3.042843176799e-001 2.528712446035e-001 -6.709925401449e-001 2.544619929082e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -7.463453456938e-001 -2.708365157019e-002 -9.486521893682e-001 1.195871081495e-001 1.929265617952e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -7.701801374364e-001 -4.697623990618e-001 4.988259468008e-001 1.137691603776e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 7.013167092107e-001 1.582180688475e-001 3.862594614233e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 4.843807602783e-001 3.992777995177e-0012. 矩阵A各个特征值的实部及虚部为：3. 经带双步位移的QR分解后的矩阵为：3.383039617436e+000 8.948776735202e-001 -8.956760626575e-001 8.465153677934e-002 2.612677876847e-001 1.610398501091e+000 -1.022613073552e+000 9.371886210517e-002 -1.002578700827e+000 -4.086260812135e-0010.000000000000e+000 -2.118477444213e+000 -2.361527904476e+000 -3.455612566063e-002 -4.736641031449e-002 1.816397532103e+000 -2.318977562028e-001 -1.435516012783e-001 -6.537076947730e-001 3.227152128950e-0020.000000000000e+000 3.555130807938e-001 -2.528514976210e+000 -6.375262961974e-001 2.023867361503e-002 1.838634248158e+000 1.868763022455e-001 -2.932577853277e-001 1.987065844256e+000 1.004631189744e+0000.000000000000e+000 0.000000000000e+000 0.000000000000e+000 1.577548557113e+000 -1.396212144889e-002 -6.971943990408e-001 1.555685237633e-001 8.405054265524e-003 -8.154391705514e-002 -1.086119035359e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -1.484039822259e+000 -1.005291633153e-001 4.249889459973e-002 2.623603124719e-002 1.040075636749e-001 -1.180720843108e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -6.948721401107e-001 -5.289206673609e-001 2.679156722087e-001 -5.962368505233e-001 -4.911592278103e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 1.788270336426e-001 -1.266189772470e+000 4.715000043047e-002 2.904780036855e-001 -3.570390801052e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 9.355889078188e-001 1.877411051033e-001 1.360256519262e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 6.360511809728e-001 2.737034413363e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+0000.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 2.457609935855e-005 5.651649653672e-0024. 特征向量如下：特征值 3.383039617436e+000对应的特征向量为：（2.383044269616e+000 4.356301670926e+000 2.966232952005e+000 3.414040144924e+000 3.192820812991e+000 3.460517224709e+000 3.228349157601e+000 3.896462664798e+000 3.720436504998e+000 1.407124666763e-001）特征值 1.577548557113e+000对应的特征向量为：（5.775526303775e-001 2.527748803591e+000 8.788044927043e-001 1.776934466823e+000 1.562614187082e+000 1.952049971107e+000 2.083198669547e+000 1.094012955120e+000 1.161445426852e+000 -1.757901376022e-001）特征值 -1.484039822259e+000对应的特征向量为：（-2.48403