北航数值分析大作业第二题(fortran).doc

《北航数值分析大作业第二题(fortran).doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵A(n-1)，给定精度水平ε=10-12和最大迭代次数L；2、记A1=A(n-1)=aij(1)n×n，令k=1，m=n；3、如果m≤2，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果am,m-1(k)≤ε，则

2、可得一个特征值，置m=m-1；转3，否则转5；5、如果am-1,m-2(k)≤ε，则可得两个特征值，置m=m-2；转3，否则转6；6、记Ak=aij1m×m(1≤i,j≤m)，计算s=am-1,m-1(k)+am,m(k)t=am-1,m-1(k)am,m(k)-am,m-1(k)am-1,m(k)Mk=Ak2-sAk+tI(I是m阶单位矩阵)Mk=QkRk(对Mk作QR分解)Ak+1=QkTAkQk7、k=k+1,转38、A的全部特征值已经求出，停止计算。二、计算源程序（FORTRAN）PROGRAMSY1415215_2PAR

3、AMETER(N=10)DIMENSIONA(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部REAL(8)A,A1,A2,C,Q,R,CME=1E-12!精度水平L=1000!迭代最大次数OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')DOI=1,NDOJ=1,NIF(I==J)THENA(I,J)=1.52*COS(I+1.2*J)ELSEA(I,J)=SIN(0.5*I+0.2*J)ENDIFENDDOENDD

4、OA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DOI=1,NDOJ=1,NWRITE(*,"(2X,E20.13,2X,)")A(I,J)WRITE(1,"(2X,E20.13,2X,)")A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALLHESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（

5、n-1）为：')")DOI=1,NDOJ=1,NWRITE(*,"(2X,E20.13,2X,)")A(I,J)WRITE(1,"(2X,E20.13,2X,)")A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵A2=ACALLQRD(A2,N,Q,R)WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')")WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')")DOI=1,NDOJ

6、=1,NWRITE(*,"(2X,E20.13,2X,)")Q(I,J)WRITE(1,"(2X,E20.13,2X,)")Q(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOWRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')")WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')")DOI=1,NDOJ=1,NWRITE(*,"(2X,E20.13,2X,)")R(I,J)WRITE(1,"(2X,E20.13,2X,)")R(I,

7、J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用带双步位移的QR方法求解矩阵A（n-1）的特征值K=1M=NDOWHILE(K<=L)IF(M<=2)THENIF(M==1)THENC(1,M)=A(M,M)ELSEIF(M==2)THENCALLCALCUS(A,N,M,C)ENDIFEXITELSEIF(ABS(A(M,M-1))

8、LSECALLCALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DOJ=1,NWRITE(*,"(E20.13,'+',E20.13,'i')