北航数值分析大作业二

1、北京航空航天大学2014年研究生数值分析作业第二题院系：XXXXXXXXX学院学号：XXXXXXXXX姓名：XXXXEmail:XXXXXXXXXXXXXXX2014年11月算法设计方案1、矩阵A的输入与存储根据题目所给的条件输入需要求解的10阶矩阵A，即为源程序中的void assignment()子程序。2、将矩阵A拟上三角化利用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵，由子程序void Hessenberg()完成。矩阵A的拟上三角化算法：记A，并记。对r=1，2，n-2执行（1）若全为零，则令，转（5）；否则转（2）。（2）计算，（若，则取），（3）令（4）计算，（5）继续3、对

2、拟上三角化后的矩阵进行QR分解为了了解普通的QR分解过程及结果，下述程序中用void QRfenjie()子程序来对拟上三角化过后的A阵进行QR分解，并输出Q阵R阵RQ阵。其中QRfenjie()既用于基本的QR分解和输出Q阵、R阵、RQ阵，又用于带双步位移的QR方法。m为维数，m<=10。当n=0时，子程序进行基本QR分解；n0时，子程序进行带双步位移的QR分解。4、对拟上三角化后的矩阵进行带双步位移的QR分解。子程序void QRfa()实现对拟上三角化后的A阵进行带双步位移的QR分解，得出特征值并输出，并用子程序void vector()对其中的实数特征值进行求解，得出对应的特征向

3、量。带双步位移的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

4、）记，计算：（为m阶单位矩阵）（对作QR分解）（10）置k=k+1，转（3）。（11）A的全部特征值以计算完毕，停止计算。注：对作QR分解与的计算算法为：记，。对于r=1，2，m-1执行（1）若全为零，则令，转（5）；否则转（2）。（2）计算，（若，则取），（3）令（4）计算，（5）继续源程序#include<stdio.h>#include<iostream>#include<stdlib.h> #include<math.h> #include<float.h>#include<iomanip>#include<

5、time.h>#define N 10 /*定义矩阵的阶数*/#define E 1.0e-12 /*定义全局变量相对误差限*/#define L 20 /*迭代次数*/FILE *fp;struct C /*定义结构体*/double R;double I;void shuchu(double aNN) /*输出一个10*10的矩阵*/int i,j;double b=0;for(i=0;i<N;i+)for(j=0;j<N;j+)if(fabs(aij)<E)fprintf(fp,"%+.12e",b);elsefprintf(fp,"

6、%+.12e",aij);if(j+1)%3=0)fprintf(fp,"n");fprintf(fp,"n");int sign(double x) /*符号函数*/if(x<0) return -1;if(x = 0) return 0;if(x>0) return 1;void assignment(double aNN) /*输入矩阵A*/int i,j;for(i=0;i<N;i+)for(j=0;j<N;j+)if(i=j)aij=1.52*cos(i+1)+1.2*(j+1);elseaij=sin(0.5

7、*(i+1)+0.2*(j+1);void Hessenberg(double aNN) /*将矩阵A拟上三角化*/int r,i,j,b;double d,c,h,t;double pN,qN,uN,wN;for(r=0;r<N-2;r+)b=0;for(i=r+2;i<N;i+)if(fabs(air)>E)b+;if(b=0) continue;elsed=0;for(i=r+1;i<N;i+)d=d+air*air;d=sqrt(d);if(ar+1r=0) c=d;else c=(-1)*sign(ar+1r)*d;h=c*c-c*ar+1r;for(i=0;

8、i<=r;i+)ui=0;ur+1=ar+1r-c;for(i=r+2;i<N;i+)ui=air;for(i=0;i<N;i+)pi=0;for(j=r+1;j<N;j+)pi=pi+aji*uj;pi=pi/h;for(i=0;i<N;i+)qi=0;for(j=r+1;j<N;j+)qi+=aij*uj;qi=qi/h;t=0;for (i=r+1;i<N;i+)t+=pi*ui;t=t/h;for(i=0;i<N;i+)wi=qi-t*ui;for(i=0;i<N;i+)for(j=0;j<N;j+)aij=aij-wi*uj

9、-ui*pj;void QRfenjie(double aNN,int m,int n) /*将m维(m10)矩阵QR分解*/ /*当n=0时用于基本QR分解；n0时用于带双步位移的QR方法*/int r,i,j,k;double d,c,h,t,s,s1;double pN,qN,uN,wN,w1N,vN,bNN=0,gNN=0;if(n=0) goto loop1;s=am-2m-2+am-1m-1;s1=am-2m-2*am-1m-1-am-1m-2*am-2m-1;for(i=0;i<m;i+)for(j=0;j<m;j+)for(r=0;r<m;r+)bij+=ai

10、r*arj;bij=bij-s*aij+s1*(i=j);goto loop2;loop1: for(i=0;i<m;i+)for(j=0;j<m;j+)bij=aij;gij=(i=j);loop2:for(r=0;r<m-1;r+)k=0;for(i=r+1;i<m;i+)if (fabs(bir)>E)k+;if(k=0) continue;d=0;for(i=r;i<m;i+)d +=bir*bir;d=sqrt(d);if(brr=0) c=d;else c=(-1)*sign(brr)*d;h=c*c-c*brr;for(i=0;i<r;i

11、+)ui=0;ur=brr-c;for(i=r+1;i<m;i+)ui=bir;for(i=0;i<m;i+)vi=0;for(j=r;j<m;j+)vi +=bji*uj;vi=vi/h;for(i=0;i<m;i+)for(j=0;j<m;j+)bij=bij-ui*vj;for(i=0;i<m;i+)pi=0;for(j=r;j<m;j+)pi +=aji*uj;pi=pi/h;for(i=0;i<m;i+)qi=0;for(j=r;j<m;j+)qi +=aij*uj;qi=qi/h;t=0;for(i=r;i<m;i+)t+

12、=pi*ui;t=t/h;for(i=0;i<m;i+)wi=qi-t*ui;for(i=0;i<m;i+)for(j=0;j<m;j+)aij=aij-wi*uj-ui*pj;if(n!=0) goto loop3;for(i=0;i<m;i+)w1i=0;for(j=r;j<m;j+)w1i+=gij*uj;for(i=0;i<m;i+)for(j=0;j<m;j+)gij=gij-w1i*uj/h;loop3:;if (n!=0) goto loop4;fprintf(fp, "Q阵：n");shuchu(g);fprintf

13、(fp, "R阵：n");shuchu(b);fprintf(fp, "RQ阵：n");shuchu(a);loop4:;void vector(double d) /*列主元高斯法求属于特征值d的特征向量，并输出*/ int i,j,k,t;double aNN,xN=0,b,m;assignment(a);for(i=0;i<N;i+)aii-=d;for(k=0;k<N-1;k+)t=k;for(i=k+1;i<N;i+)if(fabs(aik)>fabs(atk)t=i;for(j=k;j<N;j+)b=akj;ak

14、j=atj;atj=b;for(i=k+1;i<N;i+)m=aik/akk;for(j=k+1;j<N;j+)aij-=m*akj;xN-1=1;for(k=N-2;k>=0;k-)for(j=N-1;j>k;j-)xk-=akj*xj;xk=xk/akk;b=0;for(i=0;i<N;i+) /*将特征向量归一化*/b+=xi*xi;b=sqrt(b);for(i=0;i<N;i+)xi=xi/b;fprintf(fp,"属于实特征值%+.12e的特征向量为:n",d);for(i=0;i<N;i+)fprintf(fp,&q

15、uot;%+.12en",xi);void QRfa(double aNN) /*带双步位移的QR方法求矩阵的所有特征值*/ /*并输出所有实特征值的特征向量*/int i,k=1,m=N;double b,c,d;struct C yN;for (i=0;i<N;i+)yi.R=0;yi.I=0;loop2:if(fabs(am-1m-2)<=E)ym-1.R=am-1m-1;m-=1;goto loop4;else goto loop1;loop4:if(m=1)y0.R=a00;goto loop3;elseif(m>1) goto loop2;else go

16、to loop3;loop1: b=am-2m-2+am-1m-1;c=am-2m-2*am-1m-1-am-1m-2*am-2m-1;d=b*b-4*c;if (d>=0)d=sqrt(d);ym-1.R=(b+d)/2;ym-2.R=(b-d)/2;elsed=sqrt(fabs(d);ym-1.R=b/2;ym-2.R=b/2;ym-1.I=d/2;ym-2.I=(-1)*d/2;if (m=2) goto loop3;if (fabs(am-2m-3)<=E)m=m-2;goto loop4;if (k=L)fprintf(fp, "未得到A的全部特征值n&quo

17、t;);goto loop5;QRfenjie(a,m,1);k+;goto loop2;loop3:fprintf(fp,"A的全部特征值已计算完毕n");for (i=0;i<N;i+)fprintf(fp,"%+.12e+i*(%+.12e)n",yi.R,yi.I);for (i=0;i<10;i+)if (yi.I=0)vector(yi.R);loop5:;int main()double ANN;if(fp=fopen("SY1405402_谭勇洋","w")=NULL)printf(&q

18、uot;cannot open filen");exit(0);assignment(A);Hessenberg(A);fprintf(fp, " 数值分析第二次大作业n");fprintf(fp, " SY1405402 谭勇洋n");fprintf(fp, "拟上三角化后所得的矩阵A(n-1)：n");shuchu(A);QRfenjie(A, N, 0);assignment(A);Hessenberg(A);QRfa(A);fclose(fp);计算结果数值分析第二次大作业拟上三角化后所得的矩阵A(n-1)：-8.9

19、45216982281e-001 -9.933136491826e-002 -1.099831758877e+000 -7.665038709077e-001 +1.707601141456e-001 -1.934882558889e+000 -8.390208705246e-002 +9.132565113143e-001 -6.407977009188e-001 +1.946733678685e-001 -2.347878362416e+000 +2.372057921598e+000 +1.827998552316e+000 +3.266556884714e-001 +2.0823605

20、83635e-001 +2.088987009941e+000 +1.847861910289e-001 -1.263015266080e+000 +6.790694668499e-001 -4.672150886500e-001 +0.000000000000e+000 +1.735954469946e+000 -1.165023367477e+000 -1.246744443518e+000 -6.298225489084e-001 -1.984820180992e+000 +2.975750060800e-001 +6.339300596595e-001 -1.308518928772e

21、-001 +3.040301036095e-001 +0.000000000000e+000 +0.000000000000e+000 -1.292937563924e+000 -1.126239225902e+000 +1.190782911924e+000 -1.308772983895e+000 +1.860151662666e-001 +4.236733936881e-001 -1.019600826545e-001 +1.943660914505e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +

22、1.577711153032e+000 +8.169358328160e-001 +4.461531723828e-001 -4.365092541609e-002 -4.665979167188e-001 +2.941231566184e-001 -1.034421113665e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -7.728975134989e-001 -1.601028244046e+000 -2.912685474827e-001 -2.4343

23、37858321e-001 +6.736286084510e-001 +2.624772904937e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -7.296773946362e-001 -7.965456279817e-003 +9.710739102007e-001 -1.298967368574e-001 +2.780242081241e-002 +0.000000000000e+000 +0.0000000000

24、00e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +7.945539612976e-001 -4.525143454606e-001 +5.048901527575e-001 -1.211210193512e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+00

25、0 +0.000000000000e+000 +7.039911373514e-001 +1.267535523498e-001 -3.714696735513e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -4.919586872214e-001 +4.081509766399e-001 Q阵：-

26、3.560272500571e-001 +4.439873952785e-001 -6.935939248834e-001 +6.597513287486e-002 +3.701042887335e-001 +1.873680253022e-001 -1.616846253861e-002 +1.142209065421e-001 +4.846147534278e-002 -5.435281546003e-002 -9.344755733655e-001 -1.691554235406e-001 +2.642533895704e-001 -2.513596481179e-002 -1.4100

27、65879813e-001 -7.138562494120e-002 +6.160046789174e-003 -4.351719447171e-002 -1.846341016476e-002 +2.070796067082e-002 +0.000000000000e+000 -8.799213803066e-001 -4.007708665994e-001 +3.812160145536e-002 +2.138528196492e-001 +1.082645668876e-001 -9.342424307230e-003 +6.599886483483e-002 +2.8001899631

28、80e-002 -3.140602039975e-002 +0.000000000000e+000 +0.000000000000e+000 -5.371036454452e-001 -1.260096502461e-001 -7.068831837949e-001 -3.578648243181e-001 +3.088106413322e-002 -2.181569912327e-001 -9.255932185744e-002 +1.038115266702e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+00

29、0 +9.887789321494e-001 -1.266092216426e-001 -6.409685206678e-002 +5.531080075231e-003 -3.907390568777e-002 -1.657821824708e-002 +1.859359069583e-002 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +5.307477730505e-001 -6.851618290904e-001 +5.912435352117e-002 -4.1

30、76796180699e-001 -1.772124834679e-001 +1.987557610042e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -5.886032010032e-001 -9.581355780502e-002 +6.768677853804e-001 +2.871804513252e-001 -3.220922591440e-001 +0.000000000000e+000 +0.0000000

31、00000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -9.929517580142e-001 -9.993700299902e-002 -4.240112211755e-002 +4.755572027992e-002 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e

32、+000 +0.000000000000e+000 +5.375789043480e-001 -5.611563234389e-001 +6.293746914713e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +7.464002047925e-001 +6.654973585866e-001 R

33、阵：+2.512509079248e+000 -2.181265513645e+000 -1.316649918648e+000 -3.235549645996e-002 -2.553867638933e-001 -1.263236417243e+000 -1.428067524844e-001 +8.551127106196e-001 -4.064323860885e-001 +3.672920640297e-001 +0.000000000000e+000 -1.972851789714e+000 +2.276016623827e-001 +7.014634531795e-001 +5.9

34、47854062317e-001 +5.340587633903e-001 -3.303516655914e-001 +6.131164883827e-002 -2.842350072136e-001 -1.020581006877e-001 +0.000000000000e+000 +0.000000000000e+000 +2.407240343439e+000 +1.722528346094e+000 -4.505704077291e-001 +3.392449531044e+000 -1.121444617851e-001 -1.448802456850e+000 +7.3110455

35、93936e-001 -4.847285804893e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +1.595615664669e+000 +6.397406649436e-001 +3.502375215214e-001 -6.543701621428e-002 -3.985833699409e-001 +2.393366716784e-001 -9.059576583919e-002 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e

36、+000 +0.000000000000e+000 -1.456242593458e+000 -1.416208323337e+000 -2.740258546573e-001 +2.820474174287e-001 +3.146375562716e-002 +2.179593898341e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +1.239676225669e+000 +1.437893689789e-001 -

37、1.965876544047e-001 -5.501588373900e-001 -1.563867947060e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -8.001939216933e-001 +3.239240677887e-001 -4.348133701956e-001 +1.296863513999e-001 +0.000000000000e+000 +0.0000

38、00000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +1.309558711582e+000 -4.522300337931e-001 -2.541297741346e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.0000000000

39、00e+000 +0.000000000000e+000 +0.000000000000e+000 -6.591084569144e-001 +4.900007103780e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +6.373330490059e-00

40、2 RQ阵：+1.143817643298e+000 +2.643043667321e+000 -1.774014655111e+000 -7.804541680122e-002 +3.406398561353e-001 +1.461454600932e+000 -9.542620239375e-001 +4.390636920106e-001 +7.812046500198e-001 -3.242676232229e-001 +1.843581807358e+000 +1.334470111481e-001 -9.893074658739e-001 +5.579861874654e-001

41、+3.915081963850e-002 -2.951491456489e-001 +1.302107099262e-002 -6.809912257186e-001 -1.407766316366e-001 +4.534362093238e-003 +0.000000000000e+000 -2.118182245729e+000 -1.889928052623e+000 -5.708018640629e-001 +1.154750215959e+000 -2.585301662617e+000 +1.678124198417e+000 -1.154349560089e+000 -1.428

42、583810585e+000 +8.738827597003e-001 +0.000000000000e+000 +0.000000000000e+000 -8.570109902230e-001 +4.314991197035e-001 -1.023023164209e+000 -8.134730259318e-001 +4.755641447902e-001 -3.951755882759e-001 -4.241791598500e-001 +3.396131374910e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000

43、000e+000 -1.439901996509e+000 -5.672756725063e-001 +1.224964946486e+000 -3.455910828759e-001 +4.516704416505e-001 +3.294865537961e-001 -4.202787739855e-002 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +6.579553960773e-001 -9.340137131103e-001 +2.547201414468e-0

44、01 -6.965685047393e-001 +2.194090527527e-002 -2.596005659042e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +4.709967037320e-001 -2.449715460025e-001 -8.077439833310e-001 +9.726139591002e-002 +8.578610393813e-002 +0.000000000000e+000 +0.

45、000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 -1.300328624888e+000 -3.739826989665e-001 +8.562468804236e-003 -3.914678236392e-001 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000

46、000000e+000 +0.000000000000e+000 -3.543228021145e-001 +7.355995090042e-001 -8.873200325455e-002 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +0.000000000000e+000 +4.757055182990e-002 +4.241434606534

47、e-002 A的全部特征值已计算完毕+3.389613438816e+000+i*（+0.000000000000e+000）-2.336865932239e+000+i*（-8.934379210213e-001）-2.336865932239e+000+i*（+8.934379210213e-001）-1.493147080915e+000+i*（+0.000000000000e+000）+1.590313458807e+000+i*（+0.000000000000e+000）-9.891143464723e-001+i*（-1.084758631502e-001）-9.891143464

48、723e-001+i*（+1.084758631502e-001）+9.432879572769e-001+i*（+0.000000000000e+000）+4.954990923633e-002+i*（+0.000000000000e+000）+6.489488202111e-001+i*（+0.000000000000e+000）属于实特征值+3.389613438816e+000的特征向量为:-1.048719993204e-001-2.176769763196e-001-4.746940122415e-001-2.593836246507e-001-3.046652485206e-001-2.594517466617e-001+8.686641827337e-002+4.052581266927e-001+5.096282896431e-001+2.395146921660e-001属于实特征值-1.493147080915e+000的特征向量为:-5.613409816979e-001+7.781923574579e-001+1.436371665877e-002-2.776019037479e-001+3