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

北航数值分析大作业第二题(fortran) 北航数值分析大作业第二题(fortran) PAGE PAGE PAGE 24 北航数值分析大作业第二题(fortran) PAGE “数值分析“计算实习大作业第二题 ——SY1415215 孔维鹏 一、计算说明 本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。 在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示： 1、计算出A拟上三角化后的矩阵A(n-1) 2、记A1 3、如果m≤2，则可直接计算出最后1或2个特征值，转8，否则转4； 4、如果am,m-1 5、如果am-1,m-2 6、记Ak s= t= M M A 7、k=k+1,转3 8、A的全部特征值已经求出，停止计算。 二、计算源程序（FORTRAN） PROGRAM SY1415215_2 PARAMETER (N=10) DIMENSION A(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,CM E=1E-12 !精度水平 L=1000 !迭代最大次数 OPEN(1,FILE=数值分析大作业第二题计算结果.TXT) DO I=1,N DO J=1,N IF(I==J) THEN A(I,J)=*COS(I+*J) ELSE A(I,J)=SIN*I+*J) ENDIF ENDDO ENDDO A1=A WRITE(*,(矩阵A为：)) WRITE(1,(矩阵A为：)) DO I=1,N DO J=1,N WRITE(*,(2X,,2X,\)) A(I,J) WRITE(1,(2X,,2X,\)) A(I,J) ENDDO WRITE(*,( )) WRITE(1,( )) ENDDO !使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1） CALL HESSENBERG(A,N) WRITE(*,(拟上三角化后矩阵A（n-1）为：)) WRITE(1,(拟上三角化后矩阵A（n-1）为：)) DO I=1,N DO J=1,N WRITE(*,(2X,,2X,\)) A(I,J) WRITE(1,(2X,,2X,\)) A(I,J) ENDDO WRITE(*,()) WRITE(1,()) ENDDO !计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵 A2=A CALL QRD(A2,N,Q,R) WRITE(*,(对矩阵A（n-1）实行QR方法迭代结束后所得Q为：)) WRITE(1,(对矩阵A（n-1）实行QR方法迭代结束后所得Q为：)) DO I=1,N DO J=1,N WRITE(*,(2X,,2X,\)) Q(I,J) WRITE(1,(2X,,2X,\)) Q(I,J) ENDDO WRITE(*,()) WRITE(1,()) ENDDO WRITE(*,(对矩阵A（n-1）实行QR方法迭代结束后所得R为：)) WRITE(1,(对矩阵A（n-1）实行QR方法迭代结束后所得R为：)) DO I=1,N DO J=1,N WRITE(*,(2X,,2X,\)) R(I,J) WRITE(1,(2X,,2X,\)) R(I,J) ENDDO WRITE(*,()) WRITE(1,()) ENDDO !使用带双步位移的QR方法求解矩阵A（n-1）的特征值 K=1 M=N DO WHILE(K=L) IF(M=2) THEN IF(M==1) THEN C(1,M)=A(M,M) ELSE IF(M==2) THEN CALL CALCUS(A,N,M,C) ENDIF EXIT ELSE IF(ABS(A(M,M-1))E) THEN C(1,M)=A(M,M) M=M-1 ELSE IF(ABS(A(M-1,M-2))E) THEN CALL CALCUS(A,N,M,C)