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

北航数值分析大作业第二题(fortran)北航数值分析大作业第二题(fortran)北航数值分析大作业第二题(fortran)北航数值分析大作业第二题(fortran)编制仅供参考审核批准生效日期地址：电话：传真:邮编:“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵A(n-1)2、记A13、如果m≤2，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果am,m-1(k)5、如果am-1,m-26、记Aks=t=MMA7、k=k+1,转38、A的全部特征值已经求出，停止计算。二、计算源程序（FORTRAN）PROGRAMSY1415215_2PARAMETER(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)=*COS(I+*J)ELSEA(I,J)=SIN*I+*J)ENDIFENDDOENDDOA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DOI=1,NDOJ=1,NWRITE(*,"(2X,,2X,\)")A(I,J)WRITE(1,"(2X,,2X,\)")A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALLHESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")DOI=1,NDOJ=1,NWRITE(*,"(2X,,2X,\)")A(I,J)WRITE(1,"(2X,,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=1,NWRITE(*,"(2X,,2X,\)")Q(I,J)WRITE(1,"(2X,,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,,2X,\)")R(I,J)WRITE(1,"(2X,,2X,\)")R(I,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))<E)THENC(1,M)=A(M,M)M=M-1ELSEIF(ABS(A(M-1,M-2))<E)THENCALLCALCUS(A,N,M,C)M=M-2ELSECALLCALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DOJ=1,NWRITE(*,",'+',,'i')")C(1,J),C(2,J)WRITE(1,",'+',,'i')")C(1,J),C(2,J)ENDDO!使用反幂法求解A的相应于实特征值的特征向量J=1DOI=1,NIF(C(2,I)==0)THENCR(J)=C(1,I)J=J+1ENDIFENDDOJC=J-1WRITE(*,"('矩阵A的实特征值为：')")WRITE(1,"('矩阵A的实特征值为：')")DOI=1,JCWRITE(*,"")CR(I)WRITE(1,"")CR(I)ENDDODOII=1,JCDOI=1,NDOJ=1,NIF(I==J)THENA(I,J)=A1(I,J)-CR(II)ELSEA(I,J)=A1(I,J)ENDIFENDDOENDDOCALLINPOVERMETHOD(A,N,CM)WRITE(*,"('与实特征值',,'对应的特征向量为：')")CR(II)WRITE(1,"('与实特征值',,'对应的特征向量为：')")CR(II)DOI=1,NWRITE(*,"(2X,,2X,\)")CM(I)WRITE(1,"(2X,,2X,\)")CM(I)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOCLOSE(1)END!***************拟上三角化子函数*************************！SUBROUTINEHESSENBERG(A,N)DIMENSIONA(N,N),P(N),Q(N),W(N),U(N),AT(N,N)REAL(8)A,P,Q,W,U,ATREAL(8)S0,S1,S2,S3,S4,TDOL=1,N-2JUDGE=0DOI=L+2,NIF(A(I,L)/=0)THENJUDGE=1EXITENDIFENDDOIF(JUDGE==0)THENA=ACYCLEELSEIF(JUDGE/=0)THEN!计算DRS0=0DOI=L+1,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L+1,L)==0)THENCR=DRELSECR=-SGN(A(L+1,L))*DRENDIF!计算HRHR=CR**2-CR*A(L+1,L)!给u赋值DOI=1,NIF(I<L+1)THENU(I)=0ELSEIF(I==L+1)THENU(I)=A(I,L)-CRELSEIF(I>L+1)THENU(I)=A(I,L)ENDIFENDDO!计算PDOI=1,NDOJ=1,NAT(I,J)=A(J,I)ENDDOENDDODOI=1,NS1=0DOJ=1,NS1=S1+AT(I,J)*U(J)ENDDOP(I)=S1/HRENDDO!计算QDOI=1,NS2=0DOJ=1,NS2=S2+A(I,J)*U(J)ENDDOQ(I)=S2/HRENDDO!计算TS3=0DOI=1,NS3=S3+P(I)*U(I)ENDDOT=S3/HR!计算WDOI=1,NW(I)=Q(I)-T*U(I)ENDDO!计算A（r+1）DOI=1,NDOJ=1,NA(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J)ENDDOENDDOENDIFENDDORETURNEND!***************符号函数子程序*****************！FUNCTIONSGN(X)REAL(8)XIF(X>0)THENSGN=1ELSEIF(X<0)THENSGN=-1ELSEIF(X==0)THENSGN=0ENDIFEND!*********计算二阶子阵特征值s1，s2子函数*******！SUBROUTINECALCUS(X,N,M,Y)DIMENSIONX(N,N),Y(2,N)REAL(8)A,B,C,D,X,YA=1C=X(M-1,M-1)*X(M,M)-X(M-1,M)*X(M,M-1)B=-(X(M-1,M-1)+X(M,M))D=B**2-4*CIF(D>=0)THENY(1,M)=(-B-SQRT(D))/2Y(1,M-1)=(-B+SQRT(D))/2ELSEIF(D<0)THENY(1,M)=-B/2Y(1,M-1)=-B/2Y(2,M)=-SQRT(-D)/2Y(2,M-1)=-SQRT(-D)/2ENDIFRETURNEND!*********计算Mk，Ak+1子函数************！SUBROUTINECALM(A,M,N)DIMENSIONA(N,N),MK(M,M),X(M,M),QK(M,M),RK(M,M),S1(M,M),S2(M,M),QKT(M,M)REAL(8)A,MK,X,QK,RK,QKTREAL(8)S0,S1,S2DOI=1,MDOJ=1,MIF(I==J)THENX(I,J)=1ELSEX(I,J)=0ENDIFENDDOENDDOS=A(M-1,M-1)+A(M,M)T=A(M-1,M-1)*A(M,M)-A(M,M-1)*A(M-1,M)DOI=1,MDOJ=1,MS0=0DOK=1,MS0=S0+A(I,K)*A(K,J)ENDDOMK(I,J)=S0-S*A(I,J)+T*X(I,J)ENDDOENDDO!对Mk做QR分解CALLQRD(MK,M,QK,RK)DOI=1,MDOJ=1,MQKT(I,J)=QK(J,I)ENDDOENDDODOI=1,MDOJ=1,MS1(I,J)=0DOK=1,MS1(I,J)=S1(I,J)+QKT(I,K)*A(K,J)ENDDOENDDOENDDOA=S1DOI=1,MDOJ=1,MS2(I,J)=0DOK=1,MS2(I,J)=S2(I,J)+A(I,K)*QK(K,J)ENDDOENDDOENDDOA=S2RETURNEND!************QR分解子程序***************！SUBROUTINEQRD(A,N,Q,R)DIMENSIONA(N,N),AT(N,N),Q(N,N),U(N),W(N),P(N),R(N,N)REAL(8)A,AT,Q,U,W,P,RREAL(8)DR,S0,CR,HR,S1,S2DOI=1,NDOJ=1,NIF(I==J)THENQ(I,J)=1ELSEQ(I,J)=0ENDIFENDDOENDDODOL=1,N-1JUDGE=0DOI=L+1,NIF(A(I,L)/=0)THENJUDGE=1EXITENDIF!A(I,L)中有一个不为零,判断条件为真，跳出循环转ENDDOIF(JUDGE==0)THENQ=QA=ACYCLE!A(I,L)全为零，结束本循环，进入下一个ELSEIF(JUDGE/=0)THEN!计算DRS0=0DOI=L,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L,L)==0)THENCR=DRELSECR=-SGN(A(L,L))*DRENDIF!计算HRHR=CR**2-CR*A(L,L)!给u赋值DOI=1,NIF(I<L)THENU(I)=0ELSEIF(I==L)THENU(I)=A(I,L)-CRELSEIF(I>L)THENU(I)=A(I,L)ENDIFENDDO!计算WDOI=1,NS1=0DOJ=1,NS1=S1+Q(I,J)*U(J)ENDDOW(I)=S1ENDDO!计算Q（r+1）DOI=1,NDOJ=1,NQ(I,J)=Q(I,J)-W(I)*U(J)/HRENDDOENDDO!计算PDOI=1,NDOJ=1,NAT(I,J)=A(J,I)ENDDOENDDODOI=1,NS2=0DOJ=1,NS2=S2+AT(I,J)*U(J)ENDDOP(I)=S2/HRENDDO!计算A（r+1）DOI=1,NDOJ=1,NA(I,J)=A(I,J)-U(I)*P(J)ENDDOENDDOENDIFENDDOQ=QR=ARETURNEND!*************运用反幂法求解矩阵A实特征值的特征向量***********！SUBROUTINEINPOVERMETHOD(A,N,Y)DIMENSIONA(N,N),U(N),Y(N),U1(N,N),L1(N,N)REAL(8)E,Z,Z1,Z2,S1,S2,BREAL(8)A,U,Y,U1,L1DOI=1,NU(I)=1ENDDO!任取非零向量UCALLDETA(A,N,U1,L1)Z2=EIZ1=E=1K=1DOWHILE(E>1E-12)S1=0DOI=1,NS1=S1+U(I)**2ENDDOB=SQRT(S1)!1DOI=1,NY(I)=U(I)/BENDDO!2CALLDOOLITTLE(U1,L1,Y,N,U)!3利用DOOLITTLE分解法法求解Au=yS2=0DOI=1,NS2=S2+Y(I)*U(I)ENDDOZ1=Z2Z2=S2!4E=ABS(1/Z2-1/Z1)/ABS(1/Z2)!判断是否满足精度K=K+1ENDDORETURNENDSUBROUTINEDOOLITTLE(U,L,B1,N,X)DIMENSIONB(N),U(N,N),X(N),Y(N),B1(N)REAL(8)L(N,N)REAL(8)B,U,X,Y,B1REAL(8)S1,S2,S3,S4B=B1Y(1)=B(1)DOI=2,NS3=0 DOM=1,I-1 S3=S3+L(I,M)*Y(M) ENDDOY(I)=B(I)-S3ENDDOX(N)=Y(N)/U(N,N)DOI=N-1,1,-1S4=0 DOM=I+1,N S4=S4+U(I,M)*X(M) ENDDOX(I)=(Y(I)-S4)/U(I,I)ENDDORETURNENDSUBROUTINEDETA(A1,N,U,L)DIMENSIONA(N,N),U(N,N),A1(N,N)REAL(8)L(N,N)REAL(8)X,S1,S2REAL(8)A,U,A1X=1A=A1!对矩阵A进行Doolittle分解DOK=1,N DOJ=K,N S1=0 DOM=1,K-1 S1=S1+L(K,M)*U(M,J) ENDDO U(K,J)=A(K,J)-S1 A(K,J)=U(K,J) ENDDO IF(K==N)THEN EXIT ELSE DOI=K+1,N S2=0 DOM=1,K-1 S2=S2+L(I,M)*U(M,K) ENDDO L(I,K)=(A(I,K)-S2)/U(K,K) A(I,K)=L(I,K) ENDDO ENDIFENDDORETURNEND三、计算结果矩阵A为：+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+01+00+00+00+00+00+00+00+00+00+01+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+01+00+00+00+00+00+00+00+00+00+00+01+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+01拟上三角化后矩阵A（n-1）为：+00+01+00+00+01+00+00+00+01+01+01+00+00+01+00+01+00+00+01+01+01+00+01+00+00+00+00+01+01+01+01+00+00+00+00+01+00+00+00+00+00+00+01+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00对矩阵A（n-1）实行QR方法迭代结束后所得Q为：+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00+00