数值分析报告(共21页).doc

《数值分析报告(共21页).doc》由会员分享，可在线阅读，更多相关《数值分析报告(共21页).doc（21页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上北航2008级研究生数值分析A计算实习题2008年11月15日目 录4计算5671 引言算法背景：幂法，Jacobi方法及QR方法是求矩阵特征值和特征向量的常用数值方法，它们都是造构造迭代产生的矩阵序列来达到目的的。幂法计算简单，特别适用于高阶稀疏矩阵，但其收敛速度不能令人满意，要想加快幂法的收敛速度可采用反幂法及位移技术。Jacobi方法是古典方法，它收敛快，精度高，便于并行计算且算法稳定。用Jacobi方法求出的特征向量在较好的正交性，不过它的计算量较大，当阶数n增大时收敛速度减慢，因此Jacobi方法适用于求低阶的对称矩阵的全部特征值和特征向量。由J.G.Fra

2、ncis提出的QR算法的基本思想源于LR算法，即对矩阵分解获得一个相似于原矩阵的序列，使其收敛到一个易于求得特征值的形式。LR算法比QR算法收敛速度快，但不稳定。QR方法是60年代发展起来的，被人们称为数值数学，最值得注意的算法之一，它是目前求任意矩阵全部特征值和特征向量最有效的方法。利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换，使其变化为一个近似的上三角矩阵来求全部特征值。这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。但是基本QR算法的收敛往往过慢，因此通常采用带原点位移的QR算法。又由于A是实矩阵，它可能有复特征值，如果采用一般带原点位移的QR算法

3、，变换中带有复位移量，则造成复数运算。故而采用带双步位移的QR算法，可以减少迭代次数加速收敛，避免复数运算，在计算上是比较经济的。2 算法设计方案2.1 题目 试用带双步位移的QR分解法求矩阵A=的全部特征值，并对其中的每一个实特征值求相应的特征向量。要求程序输出矩阵A经过拟上三角化后所得到的矩阵；对矩阵进行QR分解方法结束后所得的矩阵；矩阵A的全部特征值及相应于实特征值的特征向量。采用e型输出数据，并且至少显示12位有效数字。要求迭代的精度水平为=。已知 2.2 算法2.2.1用带双步位移的QR方法求实矩阵全部特征值的算法（1）使用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵；给定精度水平

4、和迭代最大次数L。【矩阵A的拟上三角化的算法如下：记A，并记。对于r=1，2，n-2执行（a）若全为零，则令，转（e）；否则转（b）。（b）计算，（若，则取），（c）令（d）计算，（e）继续】（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

5、）如果k=L，则计算终止，未得到A的全部特征值，否则转（9）。（9）记，计算（为m阶单位矩阵）（对作QR分解）【对作QR分解与的计算算法如下：记，。对于r=1，2，m-1执行（a）若全为零，则令，转（e）；否则转（b）。（b）计算，（若，则取），（c）令（d）计算，（e）继续】（10）置k=k+1，转（3）。（11）A的全部特征值以计算完毕，停止计算。2.2.2 列主元素Gauss消去法求特征值对应的特征向量的算法列主元素Gauss消去法具体算法如下：记，1.消元过程对于1，2，n-1执行（1）选行号，使（2）交换与以及与所含的数值（3）对于i=k+1,k+2, ,n计算，2.回代过程3特殊情

6、况处理 (1)由于拟上三角化和QR分解的算法有相同的部分故将两部分的子程序集成在一起（当m=1时为矩阵的拟上三角化程序和当m=0时为QR分解程序）(2)当输出“未求出矩阵A的全部特征值”时，检查L的大小，增大L的数值再计算。4 计算结果4.1 矩阵A经过拟上三角化后得到的矩阵 （由于一行放不下十个元素，故将矩阵每一行的十个元素分成两行列出）矩阵A经过拟上三角化后所得的矩阵A(n-1)：-8.e-001-9.e-002-1.e+000-7.e-0011.e-001-1.e+000-8.e-002-9.e-0016.e-0011.e-001-2.e+0002.e+0001.e+0003.e-001

7、2.e-0012.e+0001.e-0011.e+000-6.e-001-4.e-001-1.e-0161.e+000-1.e+000-1.e+000-6.e-001-2.e+0002.e-001-6.e-0019.e-0022.e-001-5.e-017-3.e-017-1.e+000-1.e+0001.e+000-1.e+0001.e-001-4.e-0017.e-0021.e-0011.e-0175.e-017-3.e-0171.e+0008.e-0014.e-001-3.e-0024.e-001-2.e-001-7.e-0021.e-016-3.e-0178.e-0170.e+000-

8、7.e-001-1.e+000-3.e-0012.e-001-6.e-0012.e-0011.e-016-2.e-016-3.e-0170.e+0000.e+000-7.e-001-2.e-002-9.e-0011.e-0011.e-0021.e-0167.e-017-8.e-0170.e+0000.e+000-1.e-016-7.e-001-4.e-0014.e-0011.e-001-2.e-017-6.e-0178.e-0170.e+0000.e+000-4.e-0171.e-0167.e-0011.e-0013.e-001-2.e-017-1.e-016-2.e-0170.e+0000.

9、e+0004.e-0171.e-0160.e+0004.e-0013.e-0014.2 对矩阵进行QR分解后所得的矩阵（由于一行放不下十个元素，故将矩阵每一行的十个元素分成两行列出）对拟三角化后的矩阵A(n-1)进行QR分解方法结束后的矩阵：1.e+000-2.e+000-5.e-002-1.e-0013.e-001-4.e-001-4.e-0011.e-0014.e-0029.e-002-3.e+000-5.e-001-2.e-001-6.e-0011.e+0002.e+000-7.e-001-8.e-001-7.e-0012.e+0003.e-0161.e+000-1.e+0007.e-0

10、01-6.e-001-1.e+0002.e-0013.e-0015.e-001-9.e-001-2.e-016-2.e-016-4.e-001-5.e-0011.e+0005.e-001-3.e-001-1.e-001-4.e-0025.e-0011.e-016-9.e-0183.e-0161.e+0003.e-001-1.e-0011.e-001-1.e-001-2.e-0011.e-001-1.e-017-7.e-017-8.e-017-4.e-017-6.e-001-9.e-0012.e-001-5.e-002-8.e-002-3.e-0011.e-0161.e-0161.e-016-4

11、.e-017-2.e-016-4.e-001-1.e+000-2.e-0012.e-001-4.e-0014.e-0174.e-0174.e-017-1.e-017-1.e-016-6.e-017-4.e-0018.e-001-1.e-0011.e-0012.e-0174.e-0172.e-017-1.e-017-7.e-017-3.e-017-1.e-016-1.e-0019.e-0023.e-003-2.e-017-5.e-017-1.e-0181.e-017-8.e-017-3.e-0172.e-017-5.e-016-2.e-0017.e-0014.5 矩阵A的全部特征值矩阵A的全部特

12、征值:6.e-001+0.e+000i5.e-002+0.e+000i9.e-001+0.e+000i-9.e-001+1.e-001i-9.e-001-1.e-001i-1.e+000+0.e+000i1.e+000+0.e+000i-2.e+000+8.e-001i-2.e+000-8.e-001i3.e+000+0.e+000i4.6 矩阵A的相应于实特征值的特征向量（由于一行放不下十个元素，故将矩阵每一行的十个元素分成两行列出）特征值对应的特征向量(按行排列)(不是实特征值的,不求对应的特征向量):-1.e-001-7.e-002-3.e-0014.e-0027.e-001-1.e-0

13、012.e-001-3.e-001-2.e-001-2.e-0022.e-0012.e-001-3.e-0012.e-0023.e-0011.e-001-6.e-0013.e-0012.e-001-4.e-002-8.e-002-4.e-0021.e-0024.e-0023.e-001-2.e-0011.e-001-8.e-0013.e-001-2.e-002不是实特征值,不求对应的特征向量不是实特征值,不求对应的特征向量5.e-001-7.e-001-1.e-0022.e-001-3.e-0032.e-0032.e-0021.e-0021.e-002-3.e-002-6.e-0021.e-0

14、022.e-0011.e-0013.e-001-8.e-0011.e-0016.e-002-2.e-001-1.e-001不是实特征值,不求对应的特征向量不是实特征值,不求对应的特征向量1.e-0012.e-0014.e-0012.e-0013.e-0012.e-001-8.e-002-4.e-001-5.e-001-2.e-0015 结论QR方法是目前求任意矩阵全部特征值和特征向量最有效的方法。利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换，使其变化为一个近似的上三角矩阵（拟上三角矩阵）来求全部特征值，算出特征值后对实特征值利用列主元素Gauss消去法求实特征值对应的特征向

15、量。利用QR分解将矩阵化为一个正交矩阵Q和一个上三角矩阵R左乘的形式。从QR 算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值，那麽涉及的计算量是很大的，因此应该先对A作预处理。应用中常先对 做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素比三角矩阵多了一条次对角线。采用此程序计算题目所给矩阵的特征值时，迭代次数k=20，可见采用带双步位移的QR算法，可以减少迭代次数加速收敛，避免复数运算，在计算上是比较经济的。6 参考文献【1】 冯康等编：数

16、值计算方法，国防工业出版社，1978年。【2】 颜庆津编：数值分析（第三版），北京航空航天大学出版社，2006年。【3】 谭浩强编：C程序设计（第二版），清华大学出版社，2004年。7 附录附录A 全部源程序（包括主程序和每个子程序的功能注释）/带双步位移的QR算法求矩阵特征值和特征向量程序/声明：本程序为北航2008级研究生数值分析A计算实习题第一题。# include # include void input(int m,double n10)/数据输入子程序（将aij的数据输入二维数组）int i,j;for(i=1;im+1;i+)for(j=1;j0)b=1;elseb=-1;ret

17、urn b;void mul(int n1,int n2,int n3, double a10,double m110,double m210)/定义矩阵乘法的子程序int i,j,l;for(i=0;in1;i+)for(j=0;jn3;j+)aij=0;for(l=0;ln2;l+)aij=aij+m1il*m2lj;void equ(double A10,double B10)/定义矩阵赋值的子程序int i,j;for(i=0;i10;i+)for(j=0;j10;j+)Aij=Bij;void sub(int n1,int n2,double m,double a10,double

18、m110,double m210)/定义矩阵减法的子程序int i,j;for(i=0;in1;i+)for(j=0;jn2;j+)aij=0;aij=m1ij-m*m2ij;void add(int n1,int n2,double m,double a10,double m110,double m210)/定义矩阵加法的子程序int i,j;for(i=0;in1;i+)for(j=0;j=0)x00=(a+sqrt(a*a-4*b)/2;x01=0;x10=(a-sqrt(a*a-4*b)/2;x11=0;elsex00=a/2;x01=(sqrt(4*b-a*a)/2;x10=a/2;

19、x11=(-sqrt(4*b-a*a)/2;void qr(int n,int m,double B10,double C10)/矩阵的拟上三角化和qr分解的集成子程序(当m=1时为拟上三角化m=0时为qr分解)double Br1010,BrT1010,Br11010;double Cr1010,CrT1010,Cr11010;double vr1010,vrT1010,ur1010,urT1010,uv1010;double pr1010,prT1010,qr1010,tr1010,wr1010;double wu1010,up1010;double dr,cr,hr;int i,j,r,

20、l;equ(Br,B);equ(Cr,C);for(r=1;rn-m;r+)Br,Cr;for(i=r+1+m;in+1;i+)l=0;if(Bri-1r-1=0)l=1;continue;else break;if(l=1)continue;elsedr=0;for(i=r+m;in+1;i+)dr=dr+Bri-1r-1*Bri-1r-1;dr=sqrt(dr);cr=-sign(Brr+m-1r-1)*dr;hr=cr*cr-cr*Brr+m-1r-1;for(i=1;in+1;i+)if(ir+m)uri-10=0;else if(i=r+m)uri-10=Bri-1r-1-cr;el

21、seuri-10=Bri-1r-1;if(m=0) /m=0时为qr分解for(i=1;in+1;i+)for(j=1;jn+1;j+)BrTi-1j-1=Brj-1i-1;mul(n,n,1,vr,BrT,ur);for(i=1;in+1;i+)vri-10=vri-10/hr;for(i=1;in+1;i+)vrT0i-1=vri-10;mul(n,1,n,uv,ur,vrT);sub(n,n,1,Br1,Br,uv);equ(Br,Br1);for(i=1;in+1;i+)for(j=1;jn+1;j+)CrTi-1j-1=Crj-1i-1;mul(n,n,1,pr,CrT,ur);fo

22、r(i=1;in+1;i+)pri-10=pri-10/hr;mul(n,n,1,qr,Cr,ur);for(i=1;in+1;i+)qri-10=qri-10/hr;for(i=1;in+1;i+)prT0i-1=pri-10;mul(1,n,1,tr,prT,ur);tr00=tr00/hr;sub(n,1,tr00,wr,qr,ur);for(i=1;in+1;i+)urT0i-1=uri-10;mul(n,1,n,wu,wr,urT);mul(n,1,n,up,ur,prT);sub(n,n,1,Cr1,Cr,wu);sub(n,n,1,Cr,Cr1,up);if(m=1) /m=1时

23、为拟上三角化equ(Br,Cr);equ(C,Cr);void gauss(int n,double x1010,double A1010,double b1010 )/列主元素gauss消去法子程序int i,j,k;int ik;double aa;double AK1010,bK1010,m1010;equ(bK,b);equ(AK,A);for(k=1;kn;k+)aa=AKk-1k-1;ik=k;for(i=k+1;iaa)aa=AKi-1k-1;ik=i;elsecontinue;for(j=1;jn+1;j+)aa=AKk-1j-1;AKk-1j-1=AKik-1j-1;AKik

24、-1j-1=aa;aa=bKk-10;bKk-10=bKik-10;bKik-10=aa;for(i=k+1;in+1;i+)mi-1k-1=AKi-1k-1/AKk-1k-1;for(j=k;j0;k-)aa=0;for(j=k+1;jn+1;j+)aa=aa+AKk-1j-1*xj-10;xk-10=(bKk-10-aa)/AKk-1k-1;void main()/主函数double ev102,ss22;double x1010,xx1010,b1010;double A1010,A11010,A222;double M1010,M11010,AA1010,I1010;double ep

25、s=1.0e-12,x1;/定义精度水平=1.0e-12double s,t;int i,j,k,r=1;int n=10,m=n,l=1,L=40;/定义迭代最大次数L=40char c1=+,c2=i;FILE *fp;for(i=1;in+1;i+)/定义单位矩阵for(j=1;jn+1;j+)if(i-1=j-1)Ii-1j-1=1;elseIi-1j-1=0;input(m,A);/调用数据输入子程序equ(A1,A);qr(n,1,A1,A1);/进行拟上三角化if(fp=fopen(data.txt,w)=NULL)printf(cannot write this filen);

26、fprintf(fp,矩阵A经过拟上三角化后所得的矩阵A(n-1)(按行排列)：n);/拟上三角化后矩阵的输出for(i=1;in+1;i+)for(j=1;jn+1;j+)fprintf(fp,%12.11et,A1i-1j-1);if(j=5)fprintf(fp,n);fprintf(fp,nn);for(k=1;k=L;k+) /求解特征值if(fabs(A1m-1m-2)1)continue;elseA200=A1m-2m-2;A201=A1m-2m-1;A210=A1m-1m-2;A211=A1m-1m-1;function(ss,A2);/求解特征方程if(m=2)evl-10=

27、ss00;evl-11=ss01;evl0=ss10;evl1=ss11;l=l+2;m=m-2;break;else if(fabs(A1m-2m-3)1)continue;elses=A1m-2m-2+A1m-1m-1;t=A1m-2m-2*A1m-1m-1-A1m-1m-2*A1m-2m-1;mul(m,m,m,AA,A1,A1);sub(m,m,s,M1,AA,A1);add(m,m,t,M,M1,I);qr(m,0,M,A1);/ 调用qr分解子程序 if(r=1)fprintf(fp,nnn对拟三角化后的矩阵A(n-1)进行QR分解方法结束后的矩阵(按行排列)：n); /qr分解后

28、矩阵的输出for(i=1;im+1;i+)for(j=1;jm+1;j+)fprintf(fp,%12.11et,A1i-1j-1);if(j=5)fprintf(fp,n);fprintf(fp,nn);r=r+1;fprintf(fp,nnn);continue;if(m=0)break;for(i=1;in+1;i+)bi-10=0;equ(A1,A);for(i=1;in+1;i+)/求特征值对应的特征向量if(evi-11=0)sub(n,n,evi-10,A,A1,I);gauss(n,xx,A,b);for(j=1;jn+1;j+)xj-1i-1=xxj-10;for(i=1;i

29、n+1;i+)if(evi-11=0)x1=0;for(j=1;jn+1;j+)x1=x1+xj-1i-1*xj-1i-1;x1=sqrt(x1);for(j=1;j=L)fprintf(fp,未求出矩阵A的全部特征值n);fprintf(fp,矩阵A的全部特征值:n);/特征值输出for(i=1;i0)fprintf(fp,%12.11et+%12.11einn,evi-10,evi-11);elsefprintf(fp,%12.11et%12.11einn,evi-10,evi-11);fprintf(fp,nnn特征值对应的特征向量(按行排列)(不是实特征值的,不求对应的特征向量):n);for(i=1;in+1;i+)if(evi-11=0)for(j=1;jn+1;j+)fprintf(fp,%12.11et,xj-1i-1);if(j=5)fprintf(fp,n);fprintf(fp,nn);elsefprintf(fp,不是实特征值,不求对应的特征向量nn);fclose(fp);专心-专注-专业