北航数值分析大作业二.pdf

《北航数值分析大作业二.pdf》由会员分享，可在线阅读，更多相关《北航数值分析大作业二.pdf（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数值分析 B大作业二SY1103120朱舜杰一 算法设计方案：1.矩阵 A 的存储首先输入需要求解的矩阵 A。即为下述程序中的 void assignment()子程序。2.将矩阵 A 拟上三角化使用矩阵的拟上三角化的算法把矩阵 ARnn化为拟上三角矩阵A(n1)，由子程序 void Hessen()完成。矩阵 A 的拟上三角化的算法如下：(r)记A(1)A，并记A(r)aijnn。对于 r=1，2，n-2 执行（a）若ai(,rr)(i r 2,r 3,（b）计算dr)hr cr2crar(r1,r（e）；否则转（b）。,n)全为零，则令A(r1)A(r)，转ir1(an(r)2ir)(r)

2、，cr sign(ar(r1,r)dr（若ar1,r 0，则取cr dr），（c）令ur(0,)(r),0,ar(r1,rcr,ar2,r,(r)Tn,an,r)R（d）计算pr A(r)Turhr，qr A(r)urhr，tr prTurhr，r qrtrur，TTA(r1)A(r)rururpr（e）继续3.对拟上三角化后的矩阵进行 QR 分解为了直观的了解普通的 QR 分解过程及结果，下述程序中用 void QRfenjie()子程序来对拟上三角化过后的 A 阵进行 QR 分解，并输出 Q 阵 R 阵 RQ 阵。其中 QRfenjie()既用于基本的 QR 分解和输出 Q 阵、R 阵、R

3、Q 阵又用于带双步位移的 QR 方法。m 为维数，m10。n=0 时，子程序进行基本 QR 分解；n0 时，子程序进行带双步位移的 QR 分解。4.对拟上三角化后的矩阵进行带双步位移的QR 分解。子程序 void QRfa()实现对拟上三角化后的 A 阵进行带双步位移的 QR 分解，得出特征值并输出，并用子程序 void xiangliang()对其中的实数特征值进行求解，得出对应的特征向量。带双步位移的 QR 方法具体算法如下：（1）使用矩阵的拟上三角化的算法把矩阵 ARnn化为拟上三角矩阵A(n1)；给定精度水平0和迭代最大次数 L。(1)（2）记A1 A(n1)aijnn，令 k=1，m

4、=n。(k)(k)（3）如果am,m1，则得到 A 的一个特征值am,m，置 m=n-1，转（4）；否则转（5）。(k)（4）如果 m=1，则得到 A 的一个特征值a1,1，转（11）；如果 m=0，则直接转（11）；如果 m1，则转（3）。(k)(k)am1,m1am1,m（5）求二阶子阵Dk(k)的两个特征值s1和s2，即计算二次方程(k)am,m1am,m(k)(k)2(ama1,m1m,m)detDk0的两个根s1和s2。（6）如果 m=2，则得到 A 的两个特征值s1和s2，转（11）；否则转（7）。(k)（7）如果am1,m2，则得到 A 的两个特征值s1和s2，置 m=m-2，转

5、（4）；否则转（8）。（8）如果 k=L，则计算终止，未得到 A 的全部特征值，否则转（9）。(k)（9）记Ak aijmm(1i,j m)，计算(k)(k)s am1,m1am,m(k)(k)(k)(k)t am1,m1am,mam,m1am1,mMk Ak2 sAktI（I为 m 阶单位矩阵）Mk QkRk（对Mk作 QR 分解）TAk1 QkAkQk【对Mk作 QR 分解与Ak1的计算算法如下：(1)(r)bB b记B1 Mk，ijrijmmmm，C1 Ak。对于 r=1，2，m-1 执行)（a）若bi(,r,r 2,r(i r 1,m)全为零，则令Br1 Br,Cr1 Cr，转（e）；

6、否则转（b）。（b）计算dr(birm(r)2ir)，cr sign(br(,rr)dr（若br(,rr)0，则取cr dr），hr cr2crbr(,rr)（c）令ur(0,r),0,br(,rr)cr,br(1,r,(r)Tm,bm,r)RTTBr1 Brurvrtr prurhr，qr Crurhr，（d）计算vr BrTurhr，pr CrTurhr，r qrtrur，Cr1 CrrurTurprT（e）继续】（10）置 k=k+1，转（3）。（11）A 的全部特征值以计算完毕，停止计算。二 源程序#include#include#include#include#include#inc

7、lude#include#define E 1.0e-12/*定义全局变量相对误差限*/#define L 20/*迭代次数*/FILE*fp;struct C/*定义结构体*/double R;double I;void shuchu(double a10)/*输出一个 10*10 的矩阵*/int i,j;double b=0;for(i=0;i10;i+)for(j=0;j10;j+)if(fabs(aij)E)fprintf(fp,%+.12e,b);elsefprintf(fp,%+.12e,aij);if(j+1)%3=0)fprintf(fp,n);fprintf(fp,n);i

8、nt sgn(double x)/*符号函数*/if(x0)return 1;void assignment(double a1010)/*输入矩阵 A*/int i,j;for(i=0;i=9;i+)for(j=0;j=9;j+)if(i=j)aij=1.5*cos(i+1)+1.2*(j+1);elseaij=sin(0.5*(i+1)+0.2*(j+1);void Hessen(double a1010)/*将一个 10*10 的矩阵拟上三角化*/int r,i,j,b;double d,c,h,t;double p10,q10,u10,w10;for(r=0;r=7;r+)b=0;fo

9、r(i=r+2;iE)b+;if(b=0)continue;elsed=0;for(i=r+1;i=9;i+)d+=air*air;d=sqrt(d);if(ar+1r=0)c=d;else c=(-1)*sgn(ar+1r)*d;h=c*c-c*ar+1r;for(i=0;i=r;i+)ui=0;ur+1=ar+1r-c;for(i=r+2;i=9;i+)ui=air;for(i=0;i10;i+)pi=0;for(j=r+1;j10;j+)pi+=aji*uj;pi=pi/h;for(i=0;i10;i+)qi=0;for(j=r+1;j10;j+)qi+=aij*uj;qi=qi/h;t

10、=0;for(i=r+1;i10;i+)t+=pi*ui;t=t/h;for(i=0;i10;i+)wi=qi-t*ui;for(i=0;i10;i+)for(j=0;j10;j+)aij=aij-wi*uj-ui*pj;void QRfenjie(double a1010,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 p10,q10,u10,w10,w110,v10,b1010=0,g1010=0;if(n=0)goto

11、 loop1;s=am-2m-2+am-1m-1;s1=am-2m-2*am-1m-1-am-1m-2*am-2m-1;for(i=0;im;i+)for(j=0;jm;j+)for(r=0;rm;r+)bij+=air*arj;bij=bij-s*aij+s1*(i=j);goto loop2;loop1:for(i=0;im;i+)for(j=0;jm;j+)bij=aij;gij=(i=j);loop2:for(r=0;rm-1;r+)k=0;for(i=r+1;iE)k+;if(k=0)continue;d=0;for(i=r;im;i+)d+=bir*bir;d=sqrt(d);if

12、(brr=0)c=d;else c=(-1)*sgn(brr)*d;h=c*c-c*brr;for(i=0;ir;i+)ui=0;ur=brr-c;for(i=r+1;im;i+)ui=bir;for(i=0;im;i+)vi=0;for(j=r;jm;j+)vi+=bji*uj;vi=vi/h;for(i=0;im;i+)for(j=0;jm;j+)bij=bij-ui*vj;for(i=0;im;i+)pi=0;for(j=r;jm;j+)pi+=aji*uj;pi=pi/h;for(i=0;im;i+)qi=0;for(j=r;jm;j+)qi+=aij*uj;qi=qi/h;t=0;f

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

14、(fp,RQ 阵：n);shuchu(a);loop4:;void xiangliang(double d)/*用列主元高斯法求属于特征值 d 的特征向量，并输出*/int i,j,k,t;double a1010,x10=0,b,m;assignment(a);for(i=0;i10;i+)aii-=d;for(k=0;k9;k+)t=k;for(i=k+1;ifabs(atk)t=i;for(j=k;j10;j+)b=akj;akj=atj;atj=b;for(i=k+1;i10;i+)m=aik/akk;for(j=k+1;j=0;k-)for(j=9;jk;j-)xk-=akj*xj;

15、xk/=akk;b=0;for(i=0;i10;i+)/*将特征向量归一化*/b+=xi*xi;b=sqrt(b);for(i=0;i10;i+)xi/=b;fprintf(fp,属于实特征值%+.12e 的特征向量如下:n,d);for(i=0;i10;i+)fprintf(fp,%+.12en,xi);void QRfa(double a1010)/*带双步位移的 QR 方法求矩阵的所有特征值*/*并输出所有实特征值的特征向量*/int i,k=1,m=10;double b,c,d;struct C y10;for(i=0;i10;i+)yi.R=0;yi.I=0;loop2:if(fa

16、bs(am-1m-2)1)goto loop2;else goto 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-=2;goto loop4;if(k=L)fprintf(fp,未得到 A

17、的全部特征值n);goto loop5;QRfenjie(a,m,1);k+;goto loop2;loop3:fprintf(fp,A 的全部特征值已计算完毕n);for(i=0;i10;i+)fprintf(fp,%+.12e+i*%+.12en,yi.R,yi.I);for(i=0;i10;i+)if(yi.I=0)xiangliang(yi.R);loop5:;main()double A1010;if(fp=fopen(SY1103120 朱舜杰,w)=NULL)printf(cannot open filen);exit(0);assignment(A);Hessen(A);fpr

18、intf(fp,数值分析 B 第二题n);fprintf(fp,SY1103120 朱舜杰n);fprintf(fp,A(n-1)n);shuchu(A);QRfenjie(A,10,0);assignment(A);Hessen(A);QRfa(A);fclose(fp);三 程序结果A(n-1):-8.827516758830e-001-9.933136491826e-002-1.103349285994e+000-7.600443585637e-001+1.549101079914e-001-1.946591862872e+000-8.782436382927e-002-9.255889

19、387184e-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-001+0.000000000000e+000+1.728274599968e+000-1.171467

20、642785e+000-1.243839262699e+000-6.399758341743e-001-2.002833079037e+000+2.924947206124e-001-6.412830068395e-001+9.783997621285e-002+2.557763574160e-001+0.000000000000e+000+0.000000000000e+000-1.291669534130e+000-1.111603513396e+000+1.171346824096e+000-1.307356030021e+000+1.803699177750e-001-4.246385

21、358369e-001+7.988955239304e-002+1.608819928069e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+1.560126298527e+000+8.125049397515e-001+4.421756832923e-001-3.588616128137e-002+4.691742313671e-001-2.736595050092e-001-7.359334657750e-002+0.000000000000e+000+0.000000000000e+000+0.000000

22、000000e+000+0.000000000000e+000-7.707773755194e-001-1.583051425742e+000-3.042843176799e-001+2.528712446035e-001-6.709925401449e-001+2.544619929082e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-7.463453456938e-001-2.708365157019e-002-9.486521

23、893682e-001+1.195871081495e-001+1.929265617952e-002+0.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-001+0.000000000000e+000+0.000000000000e+000+0.000000

24、000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+7.013167092107e-001+1.582180688475e-001+3.862594614233e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000

25、000000e+000+4.843807602783e-001+3.992777995177e-001Q 阵：-3.519262579534e-001+4.427591982245e-001-6.955982513607e-001+6.486200753651e-002+3.709718861896e-001+1.855847143605e-001-1.628942319628e-002-1.181053169648e-001-5.255375383724e-002-5.486582943568e-002-9.360277287361e-001-1.664679186545e-001+2.61

26、5299548560e-001-2.438671728934e-002-1.394774360893e-001-6.977585391241e-002+6.124472142963e-003+4.440505443139e-002+1.975907909728e-002+2.062836970533e-002+0.000000000000e+000-8.810520554692e-001-3.989762796959e-001+3.720308728479e-002+2.127794064090e-001+1.064463557221e-001-9.343171079758e-003-6.77

27、4200464527e-002-3.014340698675e-002-3.146955080444e-002+0.000000000000e+000+0.000000000000e+000-5.371806806439e-001-1.234945854205e-001-7.063151608719e-001-3.533456368496e-001+3.101438948264e-002+2.248676491598e-001+1.000601783527e-001+1.044622748702e-001+0.000000000000e+000+0.000000000000e+000+0.00

28、0000000000e+000+9.892235468621e-001-1.239411731211e-001-6.200358589825e-002+5.442272839461e-003+3.945881637235e-002+1.755813350011e-002+1.833059462907e-002+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+5.323610690264e-001-6.733900344896e-001+5.910581205868e-002+4.28

29、5425323867e-001+1.906901343193e-001+1.990794495295e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-6.059761505747e-001-9.165783032819e-002-6.645586508974e-001-2.957110877580e-001-3.087207462557e-001+0.000000000000e+000+0.000000000000e+000+0.00

30、0000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+9.933396625117e-001-9.690440311939e-002-4.311990584470e-002-4.501694411183e-002+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+5.41

31、0088006061e-001-5.817540838226e-001-6.073480580545e-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.221591336735e-001+6.917269588876e-001R 阵：+2.508342744917e+000-2.185646885493e+000-

32、1.314609070786e+000-3.558787493835e-002-2.609857850388e-001-1.283121847090e+000-1.390878610606e-001-8.712897972161e-001+3.849367902971e-001+3.353802899665e-001+0.000000000000e+000-1.961603277854e+000+2.407523727633e-001+7.054714572823e-001+5.957204318279e-001+5.526978774676e-001-3.268209924413e-001-

33、5.769498668364e-002+2.871129330189e-001-8.895128754189e-002+0.000000000000e+000+0.000000000000e+000+2.404534601993e+000+1.706758096328e+000-4.239566704091e-001+3.405332305815e+000-1.050017655852e-001+1.462257102734e+000-6.684487469283e-001-4.027646209664e-001+0.000000000000e+000+0.000000000000e+000+

34、0.000000000000e+000+1.577122080722e+000+6.399535133956e-001+3.468127872427e-001-5.701786649768e-002+4.014788054433e-001-2.222476176311e-001-6.317059236442e-002+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-1.447846997770e+000-1.415724007744e+000-2.806139044665e-001-

35、2.817910521892e-001-4.611434881851e-002+1.996629079956e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+1.231641451542e+000+1.619701003419e-001+1.962638275504e-001+5.350035621760e-001-1.509273424767e-001+0.000000000000e+000+0.000000000000e+000+

36、0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-7.753441914209e-001-3.464514508821e-001+4.312226803504e-001+1.234643696237e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+

37、1.296312940612e+000-4.288053318338e-001+2.737334158165e-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-6.707396440648e-001-4.842320121884e-001+0.000000000000e+000+0.000000000000e+000+

38、0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+7.168323926323e-002RQ 阵：+1.163074414164e+000+2.632670934508e+000-1.772796003272e+000-8.668899138521e-002+3.300503471047e-001+1.455162371214e+000-9.730650448593e

39、-001-4.873031174655e-001-7.756411630489e-001-3.249201979113e-001+1.836115060851e+000+1.144286420080e-001-9.880381403133e-001+5.589725694767e-001+4.694190067101e-002-2.978478237007e-001+1.617130577649e-002+6.936977702522e-001+1.367670571405e-001+1.419099231519e-002+0.000000000000e+000-2.118520153533e

40、+000-1.876189745783e+000-5.407071940597e-001+1.171538359721e+000-2.550323020223e+000+1.691577936540e+000+1.229951613262e+000+1.387947777212e+000+8.667502917242e-001+0.000000000000e+000+0.000000000000e+000-8.471995127808e-001+4.382910468318e-001-1.008632199185e+000-7.959374261495e-001+4.769258865577e

41、-001+4.072683083890e-001+4.096390493527e-001+3.363378940862e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-1.432244342447e+000-5.742284908055e-001+1.213151477723e+000-3.457508625575e-001-4.749853573124e-001-3.176158274191e-001-4.294507015032e-002+0.000000000000e+000+0.000000000000e

42、+000+0.000000000000e+000+0.000000000000e+000+6.556779598004e-001-9.275250974463e-001+2.529079844053e-001+6.905949216976e-001-2.374430675823e-002-2.429781119781e-001+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+4.698400884876e-001-2.730776009527e

43、-001+7.821296259798e-001-9.580964936399e-002+7.846239841323e-002+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+1.287679058937e+000-3.576058900348e-001-4.116725408807e-003+3.914268216423e-001+0.000000000000e+000+0.000000000000e

44、+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000-3.628760503545e-001+7.398980975354e-001+7.241608309576e-002+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e+000+0.000000000000e

45、+000+0.000000000000e+000-5.176670596524e-002+4.958522909877e-002A 的全部特征值已计算完毕+3.383039617436e+000+i*+0.000000000000e+000-2.323496210212e+000+i*-8.930405177200e-001-2.323496210212e+000+i*+8.930405177200e-001-1.484039822259e+000+i*+0.000000000000e+000+1.577548557113e+000+i*+0.000000000000e+000-9.80530

46、9562902e-001+i*-1.139489127430e-001-9.805309562902e-001+i*+1.139489127430e-001+9.355889078188e-001+i*+0.000000000000e+000+5.650488993501e-002+i*+0.000000000000e+000+6.360627875745e-001+i*+0.000000000000e+000属于实特征值+3.383039617436e+000 的特征向量如下:-1.052215734361e-001-2.183641973272e-001-4.730178777758e-0

47、01-2.608874788703e-001-3.058056251403e-001-2.583797832214e-001+8.733793639505e-002+4.054540338526e-001+5.090131137201e-001+2.409250445621e-001属于实特征值-1.484039822259e+000 的特征向量如下:-5.601181168003e-001+7.793415087695e-001+1.337801148572e-002-2.774092878907e-001+3.005575883015e-003-2.534834089601e-003-2.

48、062848490883e-002-1.101348291072e-002-1.224861665252e-002+3.236209304091e-002属于实特征值+1.577548557113e+000 的特征向量如下:+6.249625599173e-002-1.120900451496e-002-2.496667051788e-001-1.313644861499e-001-3.835412561202e-001+8.159443103526e-001-1.245603520675e-001-6.835610547208e-002+2.705874017772e-001+1.00519

49、1081460e-001属于实特征值+9.355889078188e-001 的特征向量如下:+8.056515482686e-002+4.611134688033e-002-1.502437919765e-002-4.812083461841e-002-3.536338920690e-001+2.089190778112e-001-1.557459663954e-001+8.196704272714e-001-3.509825128186e-001+2.885138527827e-002属于实特征值+5.650488993501e-002 的特征向量如下:-2.090175185815e-0

50、01-2.000637966849e-001+3.891760618170e-001-2.779391187071e-002-3.932365479065e-001-1.247038109520e-001+6.448109395581e-001-3.027798531081e-001-2.910952633140e-001+4.094365558118e-002属于实特征值+6.360627875745e-001 的特征向量如下:+1.070374683972e-001+7.123455832557e-002+3.902379171385e-001-4.466555137564e-002-7.