《数值分析第二次大作业(12页).doc》由会员分享，可在线阅读，更多相关《数值分析第二次大作业(12页).doc（12页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-一、二、三、 数值分析第二次大作业-第 12 页四、 题目分析使用带双步位移的QR分解法求矩阵的全部特使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： (i,j=1,2,10)五、 算法设计1. 按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为initA()；2. 为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化：对应的函数为nssj()； 3. 对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为tezhengzhi()；这个过程中包含着两个子程序：
2、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;
3、/定义矩阵的特征值数组,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
4、;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
5、+) 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;
6、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)
7、 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
8、); 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
9、_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+) fo
10、r(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
11、(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
12、; 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 =
13、 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
14、=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;
15、 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+=
16、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各个特征值的实部
17、及虚部为：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() ;
18、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.1033492859
19、94e+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.278839
20、089990e+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.
21、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
22、.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+0
23、00 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.701801374364
24、e-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
25、+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-
26、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.22715
27、2128950e-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.396
28、212144889e-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 -
29、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+00
30、0 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-0
31、01 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+000
32、0.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+00
33、0 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）特征值
34、-1.484039822259e+000对应的特征向量为：（-2.484038225520e+000 -5.099127549892e-001 -1.482753404281e+000 -1.485855328801e+000 -1.482839987797e+000 -1.484361479621e+000 -1.484058559054e+000 -1.484024423030e+000 -1.483903316970e+000 5.643548274400e-002）特征值 9.355889078188e-001对应的特征向量为：（1.935588258075e+000 1.901412
35、405301e-003 9.266810627717e-001 9.227337801051e-001 1.117265920074e+000 1.152411293463e+000 1.036047182172e+000 8.504471461628e-001 1.918379516879e+000 4.231871259198e-001）特征值 6.360627875746e-001对应的特征向量为：（1.636059572970e+000 -2.923128459354e-001 6.320831305311e-001 5.930763385190e-001 6.697581334639
36、e-001 8.041619000235e-001 7.582805783539e-001 4.869294230054e-001 1.670372957726e+000 3.319746006653e-001）特征值 5.650488993501e-002对应的特征向量为：（1.056412644529e+000 8.987792003496e-001 -2.984952851574e-001 1.924677455004e-002 1.079212926993e+000 3.452840316522e-001 9.781700222144e-001 3.387720159241e-001 -5.750939079567e-001 7.823951199462e-001）八、 学习体会前期的课程学习中对于这部分用QR分解法求取矩阵特征值的算法还有比较多的疑问和怀疑，经过这次作业自己完成了这一算法的程序并且成功的计算出了题目中给定矩阵的特征值及其相应特征向量等等相关量，使我对QR分解、高斯消元等算法有了更加清晰的了解，使我了解到在工程实践中如何求取相关量，同时也更加增添了我学习数值分析这门课的兴趣和信心，在今后的学习中我定会继续努力，争取能够早日有更大的进步。