北航数值分析大作业第二题(fortran)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北航数值分析大作业第二题(f o r t r a n)
-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN
“数值分析“计算实习大作业第二题
——SY1415215
孔维鹏一、计算说明
本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。
在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：
1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；
2、记，令k=1，m=n；
3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；
4、如果，则可得一个特征值，置m=m-1；转3，否则转5；
5、如果，则可得两个特征值，置m=m-2；转3，否则转6；
6、记，计算
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)) C(1,M)=A(M,M) M=M-1 ELSE IF(ABS(A(M-1,M-2)) CALL CALCUS(A,N,M,C) M=M-2 ELSE CALL CALM(A,M,N) ENDIF K=K+1 ENDDO WRITE(*,"('矩阵A的全部特征值为：')") WRITE(1,"('矩阵A的全部特征值为：')") DO J=1,N WRITE(*,",'+',,'i')") C(1,J),C(2,J) WRITE(1,",'+',,'i')") C(1,J),C(2,J) ENDDO !使用反幂法求解A的相应于实特征值的特征向量J=1 DO I=1,N IF(C(2,I)==0)THEN CR(J)=C(1,I) J=J+1 ENDIF ENDDO JC=J-1 WRITE(*,"('矩阵A的实特征值为：')")