matlab求解方程组,谢谢了!

&& 查看话题
fortran编程过程出错了 麻烦给解决一下 谢谢!!
程序运行过程出错了 不知道怎么办&&麻烦高手给看看。。我用的是FTN95编译器打开.f90的文件的~~& & (复件也有程序)
& & & & PROGRAM POINTEHL
& & & & DIMENSION THETA(15),EALFA(15),EBETA(15)
& & & & COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
& & & & COMMON /COMEK/EK,EAL,EBE
& & & & DATA N,PAI,Z,EAL,EBE,E1,EDA0,RX,RY,X0,XE,W0,US/65,3..68,1.0,1.0,2.21E11,0.,0.05,-2.5,1.5,39.24,1.5/
& & & & DATA THETA/10.,20.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90./
& & & & DATA EALFA/6.612,3.778,2.731,2.397,2.136,1.926,1.754,1.611,1.486,1.378,1.284,1.202,1.128,1.061,1.0/
& & & & DATA EBETA/0.319,0.408,0.493,0.53,0.567,0.604,0.641,0.678,0.717,0.759,0.802,0.846,0.893,0.944,1.0/
& & & & EK=RX/RY
& & & & AA=0.5*(1./RX+1./RY)
& & & & BB=0.5*ABS(1./RX-1./RY)
& & & & CC=ACOS(BB/AA)*180.0/PAI
& & & & DO I=1,15
& & & & IF(CC.LT.THETA(I))THEN
& & & & WRITE(*,*)I
& & & & EAL=EALFA(I-1)+(CC-THETA(I))*(EALFA(I)-EALFA(I-1))/(THETA(I)-THETA(I-1))
& & & & EBE=EBETA(I-1)+(CC-THETA(I))*(EBETA(I)-EBETA(I-1))/(THETA(I)-THETA(I-1))& & & &
& & & & GOTO 1
& & & & ENDIF
& & & & ENDDO
1& & & & EA=EAL*(1.5*W0/AA/E1)**(1./3.0)
& & & & EB=EBE*(1.5*W0/AA/E1)**(1./3.0)
& & & & PH=1.5*W0/(EA*EB*PAI)
& & & & OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
& & & & OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
& & & & OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
& & & & WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
& & & & WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
& & & & H00=0.0
& & & & MM=N-1
& & & & U=EDA0*US/(2.*E1*RX)
& & & & A1=ALOG(EDA0)+9.67
& & & & A2=5.1E-9*PH
& & & & A3=0.59/(PH*1.E-9)
& & & & B=PAI*PH*RX/E1
& & & & W=2.*PAI*PH/(3.*E1)*(B/RX)**2
& & & & ALFA=Z*5.1E-9*A1
& & & & G=ALFA*E1
& & & & AHM=1.0-EXP(-0.68*1.03)
& & & & HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)*AHM
& & & & ENDA=12.*U*(E1/PH)*(RX/B)**3
& & & & UTL=EDA0*US*RX/(B*B*2.E7)
& & & & W0=2.0*PAI*EA*EB*PH/3.0
& & & & WRITE(*,*)'& && && && && &Wait please'
& & & & CALL SUBAK(MM)
& & & & CALL MULTI(N,X0,XE,H00)
& & & & STOP
& & & & END
& & & & SUBROUTINE MULTI(N,X0,XE,H00)
& & & & DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500)
& & & & COMMON /COMEK/EK,EAL,EBE
& & & & DATA MK,G00/200,2.0943951/
& & & & G0=G00*EAL*EBE
& & & & NX=N
& & & & NY=N
& & & & NN=(N+1)/2
& & & & CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
& & & & CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
& & & & M=0
14& & & & KK=15
& & & & CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
& & & & M=M+1
& & & & CALL ERP(N,ER,P,POLD)
& & & & WRITE(*,*)'ER=',ER
& & & & IF(M.LT.MK.AND.ER.GT.1.E-5)GOTO 14
& & & & CALL OUTPUT(N,DX,X,Y,H,P)
& & & & RETURN
& & & & END
& & & & SUBROUTINE ERP(N,ER,P,POLD)
& & & & DIMENSION P(N,N),POLD(N,N)
& & & & ER=0.0
& & & & SUM=0.0
& & & & NN=(N+1)/2
& & & & DO 10 I=1,N
& & & & DO 10 J=1,NN
& & & & ER=ER+ABS(P(I,J)-POLD(I,J))
& & & & SUM=SUM+P(I,J)
10& & & & CONTINUE
& & & & ER=ER/SUM
& & & & DO I=1,N
& & & & DO J=1,N
& & & & POLD(I,J)=P(I,J)
& & & & ENDDO
& & & & ENDDO
& & & & RETURN
& & & & END
& & & & SUBROUTINE INITI(N,DX,X0,XE,X,Y,P,POLD)
& & & & DIMENSION X(N),Y(N),P(N,N),POLD(N,N)
& & & & NN=(N+1)/2
& & & & DX=(XE-X0)/(N-1.)
& & & & Y0=-0.5*(XE-X0)
& & & & DO 5 I=1,N
& & & & X(I)=X0+(I-1)*DX
& & & & Y(I)=Y0+(I-1)*DX
5& & & & CONTINUE
& & & & DO 10 I=1,N
& & & & D=1.-X(I)*X(I)
& & & & DO 10 J=1,NN
& & & & C=D-Y(J)*Y(J)
& & & & IF(C.LE.0.0)P(I,J)=0.0
10& & & & IF(C.GT.0.0)P(I,J)=SQRT(C)
& & & & DO 20 I=1,N
& & & & DO 20 J=NN+1,N
& & & & JJ=N-J+1
20& & & & P(I,J)=P(I,JJ)
& & & & DO I=1,N
& & & & DO J=1,N
& & & & POLD(I,J)=P(I,J)
& & & & ENDDO
& & & & ENDDO
& & & & RETURN
& & & & END
& & & & SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
& & & & DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
& & & & DIMENSION W(150,150),P0(150,150)
& & & & COMMON /COM1/ENDA,A1,A2,A3,Z,HM0/COMAK/AK(0:65,0:65)
& & & & COMMON /COMEK/EK,EAL,EBE
& & & & DATA NW,PAI,PAI1/150,3..2026423/
& & & & NN=(N+1)/2
& & & & CALL VI(NW,N,DX,P,W)
& & & & HMIN=1.E3
& & & & DO 30 I=1,N
& & & & DO 30 J=1,NN
& & & & RAD=X(I)*X(I)+EK*Y(J)*Y(J)
& & & & W1=0.5*RAD
& & & & H0=W1+W(I,J)
& & & & IF(H0.LT.HMIN)HMIN=H0
30& & & & H(I,J)=H0
& & & & IF(KK.EQ.0)THEN
& & & & KG1=0
& & & & H01=-HMIN+HM0
& & & & DH=0.005*HM0
& & & & H02=-HMIN
& & & & H00=0.5*(H01+H02)
& & & & ENDIF
& & & & W1=0.0
& & & & DO 32 I=1,N
& & & & DO 32 J=1,N
32& & & & W1=W1+P(I,J)
& & & & W1=DX*DX*W1/G0
& & & & DW=1.-W1
& & & & IF(KK.EQ.0)THEN
& & & & KK=1
& & & & GOTO 50
& & & & ENDIF
& & & & IF(DW.LT.0.0)THEN
& & & & KG1=1
& & & & H00=AMIN1(H01,H00+DH)
& & & & ENDIF
& & & & IF(DW.GT.0.0)THEN
& & & & KG2=2
& & & & H00=AMAX1(H02,H00-DH)
& & & & ENDIF
50& & & & DO 60 I=1,N
& & & & DO 60 J=1,NN
& & & & H(I,J)=H00+H(I,J)
& & & & IF(P(I,J).LT.0.0)P(I,J)=0.0
& & & & EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
& & & & EDA(I,J)=EDA1
55& & & & RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60& & & & EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA1)
& & & & DO 70 J=NN+1,N
& & & & JJ=N-J+1
& & & & DO 70 I=1,N
& & & & H(I,J)=H(I,JJ)
& & & & RO(I,J)=RO(I,JJ)
& & & & EDA(I,J)=EDA(I,JJ)
70& & & & EPS(I,J)=EPS(I,JJ)
& & & & RETURN
& & & & END
& & & & SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
& & & & DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
& & & & DIMENSION D(70),A(350),B(210),ID(70)
& & & & COMMON /COM1/ENDA,A1,A2,A3,Z,C3/COMAK/AK(0:65,0:65)
& & & & DATA KG1,PAI1,C1,C2/0,0..31,0.31/
& & & & IF(KG1.NE.0)GOTO 2
& & & & KG1=1
& & & & AK00=AK(0,0)
& & & & AK10=AK(1,0)
& & & & AK20=AK(2,0)
& & & & BK00=AK00-AK10
& & & & BK10=AK10-0.25*(AK00+2.*AK(1,1)+AK(2,0))
& & & & BK20=AK20-0.25*(AK10+2.*AK(2,1)+AK(3,0))
2& & & & NN=(N+1)/2
& & & & MM=N-1
& & & & DX1=1./DX
& & & & DX2=DX*DX
& & & & DX3=1./DX2
& & & & DX4=0.3*DX2
& & & & DO 100 K=1,KK
& & & & PMAX=0.0
& & & & DO 70 J=2,NN
& & & & J0=J-1
& & & & J1=J+1
& & & & JJ=N-J+1
& & & & IA=1
8& & & & MM=N-IA
& & & & IF(P(MM,J0).GT.1.E-6)GOTO 20
& & & & IF(P(MM,J).GT.1.E-6)GOTO 20
& & & & IF(P(MM,J1).GT.1.E-6)GOTO 20
& & & & IA=IA+1
& & & & IF(IA.LT.N)GOTO 8
& & & & GOTO 70
20& & & & IF(MM.LT.N-1)MM=MM+1
& & & & D2=0.5*(EPS(1,J)+EPS(2,J))
& & & & DO 50 I=2,MM
& & & & I0=I-1
& & & & I1=I+1
& & & & II=5*I0
& & & & D1=D2
& & & & D2=0.5*(EPS(I1,J)+EPS(I,J))
& & & & D4=0.5*(EPS(I,J0)+EPS(I,J))
& & & & D5=0.5*(EPS(I,J1)+EPS(I,J))
& & & & P1=P(I0,JJ)
& & & & P2=P(I1,JJ)
& & & & P3=P(I,JJ)
& & & & P4=P(I,JJ+1)
& & & & P5=P(I,JJ-1)
& & & & D3=D1+D2+D4+D5
& & & & IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
& & & & IF(H(I,J).LE.0.0)THEN
& & & & ID(I)=2
& & & & A(II+1)=0.0
& & & & A(II+2)=0.0
& & & & A(II+3)=1.0
& & & & A(II+4)=0.0
& & & & A(II+5)=1.0
& & & & A(II-4)=0.0
& & & & GOTO 50
& & & & ENDIF
& & & & IF(D1.GE.DX4)GOTO 30
& & & & IF(D2.GE.DX4)GOTO 30
& & & & IF(D4.GE.DX4)GOTO 30
& & & & IF(D5.GE.DX4)GOTO 30
& & & & ID(I)=1
& & & & IF(J.EQ.NN)P5=P4
& & & & A(II+1)=PAI1*(RO(I0,J)*BK10-RO(I,J)*BK20)
& & & & A(II+2)=DX3*(D1+0.25*D3)+PAI1*(RO(I0,J)*BK00-RO(I,J)*BK10)
& & & & A(II+3)=-1.25*DX3*D3+PAI1*(RO(I0,J)*BK10-RO(I,J)*BK00)
& & & & A(II+4)=DX3*(D2+0.25*D3)+PAI1*(RO(I0,J)*BK20-RO(I,J)*BK10)
& & & & A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
& & & & GOTO 50
30& & & & ID(I)=0
& & & & P4=P(I,J0)
& & & & IF(J.EQ.NN)P5=P4
& & & & A(II+1)=PAI1*(RO(I0,J)*AK10-RO(I,J)*AK20)
& & & & A(II+2)=DX3*D1+PAI1*(RO(I0,J)*AK00-RO(I,J)*AK10)
& & & & A(II+3)=-DX3*D3+PAI1*(RO(I0,J)*AK10-RO(I,J)*AK00)
& & & & A(II+4)=DX3*D2+PAI1*(RO(I0,J)*AK20-RO(I,J)*AK10)
& & & & A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
50& & & & CONTINUE
& & & & CALL TRA4(MM,D,A,B)
& & & & DO 60 I=2,MM
& & & & IF(ID(I).EQ.2)GOTO 60
& & & & IF(ID(I).EQ.0)GOTO 52
& & & & DD=D(I+1)
& & & & IF(I.EQ.MM)DD=0
& & & & P(I,J)=P(I,J)+C2*(D(I)-0.25*(D(I-1)+DD))
& & & & IF(J0.NE.1)P(I,J0)=P(I,J0)-0.25*C2*D(I)
& & & & IF(P(I,J0).LT.0.)P(I,J0)=0.0
& & & & IF(J1.GE.NN)GOTO 54
& & & & P(I,J1)=P(I,J1)-0.25*C2*D(I)
& & & & GOTO 54
52& & & & P(I,J)=P(I,J)+C1*D(I)
54& & & & IF(P(I,J).LT.0.0)P(I,J)=0.0
& & & & IF(PMAX.LT.P(I,J))PMAX=P(I,J)
60& & & & CONTINUE
70& & & & CONTINUE
& & & & DO 80 J=1,NN
& & & & JJ=N+1-J
& & & & DO 80 I=1,N
80& & & & P(I,JJ)=P(I,J)
& & & & CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
100& & & & CONTINUE
& & & & RETURN
& & & & END
& & & & SUBROUTINE TRA4(N,D,A,B)
& & & & DIMENSION D(N),A(5,N),B(3,N)
& & & & C=1./A(3,N)
& & & & B(1,N)=-A(1,N)*C
& & & & B(2,N)=-A(2,N)*C
& & & & B(3,N)=A(5,N)*C
& & & & DO 10 I=1,N-2
& & & & IN=N-I
& & & & IN1=IN+1
& & & & C=1./(A(3,IN)+A(4,IN)*B(2,IN1))
& & & & B(1,IN)=-A(1,IN)*C
& & & & B(2,IN)=-(A(2,IN)+A(4,IN)*B(1,IN1))*C
10& & & & B(3,IN)=(A(5,IN)-A(4,IN)*B(3,IN1))*C
& & & & D(1)=0.0
& & & & D(2)=B(3,2)
& & & & DO 20 I=3,N
20& & & & D(I)=B(1,I)*D(I-2)+B(2,I)*D(I-1)+B(3,I)
& & & & RETURN
& & & & END
& & & & SUBROUTINE VI(NW,N,DX,P,V)
& & & & DIMENSION P(N,N),V(NW,NW)
& & & & COMMON /COMAK/AK(0:65,0:65)
& & & & PAI1=0.2026423
& & & & DO 40 I=1,N
& & & & DO 40 J=1,N
& & & & H0=0.0
& & & & DO 30 K=1,N
& & & & IK=IABS(I-K)
& & & & DO 30 L=1,N
& & & & JL=IABS(J-L)
30& & & & H0=H0+AK(IK,JL)*P(K,L)
40& & & & V(I,J)=H0*DX*PAI1
& & & & RETURN
& & & & END
& & & & SUBROUTINE SUBAK(MM)
& & & & COMMON /COMAK/AK(0:65,0:65)
& & & & S(X,Y)=X+SQRT(X**2+Y**2)
& & & & DO 10 I=0,MM
& & & & XP=I+0.5
& & & & XM=I-0.5
& & & & DO 10 J=0,I
& & & & YP=J+0.5
& & & & YM=J-0.5
& & & & A1=S(YP,XP)/S(YM,XP)
& & & & A2=S(XM,YM)/S(XP,YM)
& & & & A3=S(YM,XM)/S(YP,XM)
& & & & A4=S(XP,YP)/S(XM,YP)
& & & & AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
10& & & & AK(J,I)=AK(I,J)
& & & & RETURN
& & & & END
& & & & SUBROUTINE OUTPUT(N,DX,X,Y,H,P)
& & & & DIMENSION X(N),Y(N),H(N,N),P(N,N)
& & & & NN=(N+1)/2
& & & & A=0.0
& & & & WRITE(8,110)A,(Y(I),I=1,N)
& & & & DO I=1,N
& & & & WRITE(8,110)X(I),(H(I,J),J=1,N)
& & & & ENDDO
& & & & WRITE(10,110)A,(Y(I),I=1,N)
& & & & DO I=1,N
& & & & WRITE(10,110)X(I),(P(I,J),J=1,N)
& & & & ENDDO
110 FORMAT(66(E12.6,1X))
& & RETURN
怎么错了?错误信息是什么?你期待什么? : Originally posted by pippi6 at
怎么错了?错误信息是什么?你期待什么? 直接不能运行 没有输出结果。。几个DAT文件都是空的。错误信息发到附件里了 麻烦看看,期待有结果,谢谢!
D0)FLF7U2802PMM}AX2{6WE.jpg
O05V}WLK6YAUI~~AW6BQU_U.jpg
Y4CBXDGY`CI6513HB{85[WL.jpg : Originally posted by 黑曼巴1991 at
直接不能运行 没有输出结果。。几个DAT文件都是空的。错误信息发到附件里了 麻烦看看,期待有结果,谢谢!
D0)FLF7U2802PMM}AX2{6WE.jpg
O05V}WLK6YAUI~~AW6BQU_U.jpg
Y4CBXDGY`CI6513HB{85[WL.jpg
... 把几个do loop修改了一下,把标号去掉了。可以运算了。没出现你说的问题
PROGRAM POINTEHL
&&DIMENSION THETA(15),EALFA(15),EBETA(15)
&&COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
&&COMMON /COMEK/EK,EAL,EBE
&&DATA N,PAI,Z,EAL,EBE,E1,EDA0,RX,RY,X0,XE,W0,US/65,3..68,1.0,1.0,2.21E11,0.,0.05,-2.5,1.5,39.24,1.5/
&&DATA THETA/10.,20.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90./
&&DATA EALFA/6.612,3.778,2.731,2.397,2.136,1.926,1.754,1.611,1.486,1.378,1.284,1.202,1.128,1.061,1.0/
&&DATA EBETA/0.319,0.408,0.493,0.53,0.567,0.604,0.641,0.678,0.717,0.759,0.802,0.846,0.893,0.944,1.0/
&&EK=RX/RY
&&AA=0.5*(1./RX+1./RY)
&&BB=0.5*ABS(1./RX-1./RY)
&&CC=ACOS(BB/AA)*180.0/PAI
&&DO I=1,15
& &&&IF(CC.LT.THETA(I))THEN
& && &&&WRITE(*,*)I
& && &&&EAL=EALFA(I-1)+(CC-THETA(I))*(EALFA(I)-EALFA(I-1))/(THETA(I)-THETA(I-1))
& && &&&EBE=EBETA(I-1)+(CC-THETA(I))*(EBETA(I)-EBETA(I-1))/(THETA(I)-THETA(I-1))& && &
& && &&&GOTO 1
& &&&ENDIF
1 EA=EAL*(1.5*W0/AA/E1)**(1./3.0)
&&EB=EBE*(1.5*W0/AA/E1)**(1./3.0)
&&PH=1.5*W0/(EA*EB*PAI)
&&OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
&&OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
&&OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
&&WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
&&WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
&&U=EDA0*US/(2.*E1*RX)
&&A1=ALOG(EDA0)+9.67
&&A2=5.1E-9*PH
&&A3=0.59/(PH*1.E-9)
&&B=PAI*PH*RX/E1
&&W=2.*PAI*PH/(3.*E1)*(B/RX)**2
&&ALFA=Z*5.1E-9*A1
&&G=ALFA*E1
&&AHM=1.0-EXP(-0.68*1.03)
&&HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)*AHM
&&ENDA=12.*U*(E1/PH)*(RX/B)**3
&&UTL=EDA0*US*RX/(B*B*2.E7)
&&W0=2.0*PAI*EA*EB*PH/3.0
&&WRITE(*,*)'& && && && && &Wait please'
&&CALL SUBAK(MM)
&&CALL MULTI(N,X0,XE,H00)
END PROGRAM POINTEHL
SUBROUTINE MULTI(N,X0,XE,H00)
&&DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500)
&&COMMON /COMEK/EK,EAL,EBE
&&DATA MK,G00/200,2.0943951/
&&G0=G00*EAL*EBE
&&NN=(N+1)/2
&&CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
&&CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
&&CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
&&CALL ERP(N,ER,P,POLD)
&&WRITE(*,*)'ER=',ER
&&IF(M.LT.MK.AND.ER.GT.1.E-5)GOTO 14
&&CALL OUTPUT(N,DX,X,Y,H,P)
END SUBROUTINE MULTI
SUBROUTINE ERP(N,ER,P,POLD)
&&DIMENSION P(N,N),POLD(N,N)
&&NN=(N+1)/2
&&DO&&I=1,N
& &&&DO&&J=1,NN
& && &&&ER=ER+ABS(P(I,J)-POLD(I,J))
& && &&&SUM=SUM+P(I,J)
& &&&end DO
&&ER=ER/SUM
&&DO I=1,N
& &&&DO J=1,N
& && &&&POLD(I,J)=P(I,J)
& &&&ENDDO
END SUBROUTINE ERP
SUBROUTINE INITI(N,DX,X0,XE,X,Y,P,POLD)
&&DIMENSION X(N),Y(N),P(N,N),POLD(N,N)
&&NN=(N+1)/2
&&DX=(XE-X0)/(N-1.)
&&Y0=-0.5*(XE-X0)
&&DO&&I=1,N
& &&&X(I)=X0+(I-1)*DX
& &&&Y(I)=Y0+(I-1)*DX
&&DO I=1,N
& &&&D=1.-X(I)*X(I)
& &&&DO&&J=1,NN
& && &&&C=D-Y(J)*Y(J)
& && &&&IF(C.LE.0.0)P(I,J)=0.0
& && &&&IF(C.GT.0.0)P(I,J)=SQRT(C)
& &&&end DO
&&DO&&I=1,N
& &&&DO&&J=NN+1,N
& && &&&JJ=N-J+1
& && &&&P(I,J)=P(I,JJ)
& &&&ENDDO
&&DO I=1,N
& &&&DO J=1,N
& && &&&POLD(I,J)=P(I,J)
& &&&ENDDO
END SUBROUTINE INITI
SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
&&DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
&&DIMENSION W(150,150),P0(150,150)
&&COMMON /COM1/ENDA,A1,A2,A3,Z,HM0/COMAK/AK(0:65,0:65)
&&COMMON /COMEK/EK,EAL,EBE
&&DATA NW,PAI,PAI1/150,3..2026423/
&&NN=(N+1)/2
&&CALL VI(NW,N,DX,P,W)
&&HMIN=1.E3
&&DO&&I=1,N
& &&&DO&&J=1,NN
& && &&&RAD=X(I)*X(I)+EK*Y(J)*Y(J)
& && &&&W1=0.5*RAD
& && &&&H0=W1+W(I,J)
& && &&&IF(H0.LT.HMIN)HMIN=H0
30& && &H(I,J)=H0
& &&&ENDDO
&&IF(KK.EQ.0)THEN
& &&&KG1=0
& &&&H01=-HMIN+HM0
& &&&DH=0.005*HM0
& &&&H02=-HMIN
& &&&H00=0.5*(H01+H02)
&&DO&&I=1,N
& &&&DO&&J=1,N
32& && &W1=W1+P(I,J)
& &&&ENDDO
&&W1=DX*DX*W1/G0
&&DW=1.-W1
&&IF(KK.EQ.0)THEN
& &&&GOTO 50
&&IF(DW.LT.0.0)THEN
& &&&KG1=1
& &&&H00=AMIN1(H01,H00+DH)
&&IF(DW.GT.0.0)THEN
& &&&KG2=2
& &&&H00=AMAX1(H02,H00-DH)
50 continue
&&DO&&I=1,N
& &&&DO&&J=1,NN
& && &&&H(I,J)=H00+H(I,J)
& && &&&IF(P(I,J).LT.0.0)P(I,J)=0.0
& && &&&EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
& && &&&EDA(I,J)=EDA1
55& && &RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60& && &EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA1)
& &&&end DO
&&DO&&J=NN+1,N
& &&&JJ=N-J+1
& &&&DO&&I=1,N
& && &&&H(I,J)=H(I,JJ)
& && &&&RO(I,J)=RO(I,JJ)
& && &&&EDA(I,J)=EDA(I,JJ)
70& && &EPS(I,J)=EPS(I,JJ)
& &&&end DO
END SUBROUTINE HREE
SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
&&DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
&&DIMENSION D(70),A(350),B(210),ID(70)
&&COMMON /COM1/ENDA,A1,A2,A3,Z,C3/COMAK/AK(0:65,0:65)
&&DATA KG1,PAI1,C1,C2/0,0..31,0.31/
&&IF(KG1.NE.0)GOTO 2
&&AK00=AK(0,0)
&&AK10=AK(1,0)
&&AK20=AK(2,0)
&&BK00=AK00-AK10
&&BK10=AK10-0.25*(AK00+2.*AK(1,1)+AK(2,0))
&&BK20=AK20-0.25*(AK10+2.*AK(2,1)+AK(3,0))
2 NN=(N+1)/2
&&DX1=1./DX
&&DX2=DX*DX
&&DX3=1./DX2
&&DX4=0.3*DX2
&&DO&&K=1,KK
& &&&PMAX=0.0
& &&&DO&&J=2,NN
& && &&&J0=J-1
& && &&&J1=J+1
& && &&&JJ=N-J+1
& && &&&IA=1
8& && & MM=N-IA
& && &&&IF(P(MM,J0).GT.1.E-6)GOTO 20
& && &&&IF(P(MM,J).GT.1.E-6)GOTO 20
& && &&&IF(P(MM,J1).GT.1.E-6)GOTO 20
& && &&&IA=IA+1
& && &&&IF(IA.LT.N)GOTO 8
& && &&&cycle
20& && &IF(MM.LT.N-1)MM=MM+1
& && &&&D2=0.5*(EPS(1,J)+EPS(2,J))
& && &&&DO&&I=2,MM
& && && &&&I0=I-1
& && && &&&I1=I+1
& && && &&&II=5*I0
& && && &&&D1=D2
& && && &&&D2=0.5*(EPS(I1,J)+EPS(I,J))
& && && &&&D4=0.5*(EPS(I,J0)+EPS(I,J))
& && && &&&D5=0.5*(EPS(I,J1)+EPS(I,J))
& && && &&&P1=P(I0,JJ)
& && && &&&P2=P(I1,JJ)
& && && &&&P3=P(I,JJ)
& && && &&&P4=P(I,JJ+1)
& && && &&&P5=P(I,JJ-1)
& && && &&&D3=D1+D2+D4+D5
& && && &&&IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
& && && &&&IF(H(I,J).LE.0.0)THEN
& && && && &&&ID(I)=2
& && && && &&&A(II+1)=0.0
& && && && &&&A(II+2)=0.0
& && && && &&&A(II+3)=1.0
& && && && &&&A(II+4)=0.0
& && && && &&&A(II+5)=1.0
& && && && &&&A(II-4)=0.0
& && && && &&&GOTO 50
& && && &&&ENDIF
& && && &&&IF(D1.GE.DX4)GOTO 30
& && && &&&IF(D2.GE.DX4)GOTO 30
& && && &&&IF(D4.GE.DX4)GOTO 30
& && && &&&IF(D5.GE.DX4)GOTO 30
& && && &&&ID(I)=1
& && && &&&IF(J.EQ.NN)P5=P4
& && && &&&A(II+1)=PAI1*(RO(I0,J)*BK10-RO(I,J)*BK20)
& && && &&&A(II+2)=DX3*(D1+0.25*D3)+PAI1*(RO(I0,J)*BK00-RO(I,J)*BK10)
& && && &&&A(II+3)=-1.25*DX3*D3+PAI1*(RO(I0,J)*BK10-RO(I,J)*BK00)
& && && &&&A(II+4)=DX3*(D2+0.25*D3)+PAI1*(RO(I0,J)*BK20-RO(I,J)*BK10)
& && && &&&A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
& && && &&&GOTO 50
30& && && &ID(I)=0
& && && &&&P4=P(I,J0)
& && && &&&IF(J.EQ.NN)P5=P4
& && && &&&A(II+1)=PAI1*(RO(I0,J)*AK10-RO(I,J)*AK20)
& && && &&&A(II+2)=DX3*D1+PAI1*(RO(I0,J)*AK00-RO(I,J)*AK10)
& && && &&&A(II+3)=-DX3*D3+PAI1*(RO(I0,J)*AK10-RO(I,J)*AK00)
& && && &&&A(II+4)=DX3*D2+PAI1*(RO(I0,J)*AK20-RO(I,J)*AK10)
& && && &&&A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
50& && && &CONTINUE
& && &&&end DO
& && &&&CALL TRA4(MM,D,A,B)
& && &&&DO&&I=2,MM
& && && &&&IF(ID(I).EQ.2)GOTO 60
& && && &&&IF(ID(I).EQ.0)GOTO 52
& && && &&&DD=D(I+1)
& && && &&&IF(I.EQ.MM)DD=0
& && && &&&P(I,J)=P(I,J)+C2*(D(I)-0.25*(D(I-1)+DD))
& && && &&&IF(J0.NE.1)P(I,J0)=P(I,J0)-0.25*C2*D(I)
& && && &&&IF(P(I,J0).LT.0.)P(I,J0)=0.0
& && && &&&IF(J1.GE.NN)GOTO 54
& && && &&&P(I,J1)=P(I,J1)-0.25*C2*D(I)
& && && &&&GOTO 54
52& && && &P(I,J)=P(I,J)+C1*D(I)
54& && && &IF(P(I,J).LT.0.0)P(I,J)=0.0
& && && &&&IF(PMAX.LT.P(I,J))PMAX=P(I,J)
60& && && &CONTINUE
& && &&&end DO
70& && &CONTINUE
& &&&end do
& &&&DO J=1,NN
& && &&&JJ=N+1-J
& && &&&DO I=1,N
80& && && &P(I,JJ)=P(I,J)
& && &&&end do
& &&&end do
& &&&CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
END SUBROUTINE ITER
SUBROUTINE TRA4(N,D,A,B)
&&DIMENSION D(N),A(5,N),B(3,N)
&&C=1./A(3,N)
&&B(1,N)=-A(1,N)*C
&&B(2,N)=-A(2,N)*C
&&B(3,N)=A(5,N)*C
&&DO&&I=1,N-2
& &&&IN=N-I
& &&&IN1=IN+1
& &&&C=1./(A(3,IN)+A(4,IN)*B(2,IN1))
& &&&B(1,IN)=-A(1,IN)*C
& &&&B(2,IN)=-(A(2,IN)+A(4,IN)*B(1,IN1))*C
10& &B(3,IN)=(A(5,IN)-A(4,IN)*B(3,IN1))*C
&&D(1)=0.0
&&D(2)=B(3,2)
&&DO&&I=3,N
20& &D(I)=B(1,I)*D(I-2)+B(2,I)*D(I-1)+B(3,I)
END SUBROUTINE TRA4
SUBROUTINE VI(NW,N,DX,P,V)
&&DIMENSION P(N,N),V(NW,NW)
&&COMMON /COMAK/AK(0:65,0:65)
&&PAI1=0.2026423
&&DO&&I=1,N
& &&&DO&&J=1,N
& && &&&H0=0.0
& && &&&DO&&K=1,N
& && && &&&IK=IABS(I-K)
& && && &&&DO&&L=1,N
& && && && &&&JL=IABS(J-L)
30& && && && &H0=H0+AK(IK,JL)*P(K,L)
& && && &&&end DO
& && &&&end DO
40& && &V(I,J)=H0*DX*PAI1
& &&&end DO
END SUBROUTINE VI
SUBROUTINE SUBAK(MM)
&&COMMON /COMAK/AK(0:65,0:65)
&&S(X,Y)=X+SQRT(X**2+Y**2)
&&DO&&I=0,MM
& &&&XP=I+0.5
& &&&XM=I-0.5
& &&&DO J=0,I
& && &&&YP=J+0.5
& && &&&YM=J-0.5
& && &&&A1=S(YP,XP)/S(YM,XP)
& && &&&A2=S(XM,YM)/S(XP,YM)
& && &&&A3=S(YM,XM)/S(YP,XM)
& && &&&A4=S(XP,YP)/S(XM,YP)
& && &&&AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
10& && &AK(J,I)=AK(I,J)
& &&&end DO
END SUBROUTINE SUBAK
SUBROUTINE OUTPUT(N,DX,X,Y,H,P)
&&DIMENSION X(N),Y(N),H(N,N),P(N,N)
&&NN=(N+1)/2
&&WRITE(8,110)A,(Y(I),I=1,N)
&&DO I=1,N
& &&&WRITE(8,110)X(I),(H(I,J),J=1,N)
&&WRITE(10,110)A,(Y(I),I=1,N)
&&DO I=1,N
& &&&WRITE(10,110)X(I),(P(I,J),J=1,N)
110 FORMAT(66(E12.6,1X))
END SUBROUTINE OUTPUT : Originally posted by pippi6 at
把几个do loop修改了一下,把标号去掉了。可以运算了。没出现你说的问题
PROGRAM POINTEHL
&&DIMENSION THETA(15),EALFA(15),EBETA(15)
&&COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
&&COMMON /COMEK/EK,EAL,EBE
... 我用你的试了一下 还是一样的问题啊。。:hand:&&奇怪了&&是不是我用的FTN95的原因?麻烦你给我发一下修改后的.90文件?、万分感谢!!!您还未登陆,请登录后操作!
悬赏20爱心点
分享到微博
请选择登录方式
梦见枣树上长苹果了,怎么解,谢谢!~~
了,怎么解,谢谢!~~
也许是由于晚上我和朋友一起吃苹果时说到绿苹果和红苹果了,结果晚上我就梦见,在我老家的院子里的枣树上长了苹果了,一棵结的是绿苹果,一棵结的是红苹果,但树身子是枣树,我还说是我爸嫁接的呢,请问是什么意思呢?请专业人士指点,谢谢!!~~~~~
,卖多一点钱,并通过这一个梦赞扬你父亲的聪明与勤劳的本性。
您的举报已经提交成功,我们将尽快处理,谢谢!
大家还关注

我要回帖

更多关于 结构力学求解器 的文章

 

随机推荐