程序运行出现NAN

昵称

自己编的程序代码如下：      

[Fortran] 纯文本查看 复制代码

program ex
                implicit none
        integer E,P,I,J,K,PP,NE,NP,NDT,L,MDT,ICOUNT,IPRINT,IJK(20,3),MP(3)
        REAL MU,DT,DF,TXX,TYY,AREA,MPP,D,TIME,DFMAX,OLDVAL,SUMGST,X(22),Y(22),H0(22),H(22),BB(22),B(3),C(3),G(22,22),ST(22,22)
        NE=20
        NP=18
        NDT=20
        DT=0.5
        TXX=10.0
        TYY=10.0
        MU=0.005
        OPEN(1,FILE="XZB.txt")
        READ(1,*)(X(I),I=1,NP)
        CLOSE(1)
        OPEN(1,FILE="YZB.txt")
        READ(1,*)(Y(I),I=1,NP)
        CLOSE(1)
        OPEN(1,FILE="IJKDY.txt")
        READ(1,*)((IJK(E,P),P=1,3),E=1,NE)
        CLOSE(1)
        WRITE(*,5)(X(I),I=1,NP)
   5    FORMAT(1X,"X(NP)=",/,10(9F5.1,/))
        WRITE(*,6)(Y(I),I=1,NP)
   6    FORMAT(1X,"Y(NP)=",/,10(9F5.1,/))
        WRITE(*,10)((IJK(E,P),P=1,3),E=1,NE)
   10   FORMAT(1X,"IJK(NE,3)=",/,100(4(3I4,2X),/))

                DO 20 P=1,NP
                H0(P)=16.
                20   H(P)=16.
                DO 25 P=16,18
                H0(P)=11.
   25   H(P)=11.
        DO 30 P=1,NP
                DO 30 PP=1,NP
                ST(P,PP)=0.0
   30   G(P,PP)=0.0
        DO 100 E=1,NE
                I=IJK(E,1)
                J=IJK(E,2)
                K=IJK(E,3)
                B(1)=Y(J)-Y(K)
                B(2)=Y(K)-Y(I)
                B(3)=Y(I)-Y(J)
                C(1)=X(K)-X(J)
                C(2)=X(I)-X(J)
                C(3)=X(J)-X(I)
                AREA=ABS(B(2)*C(3)-B(3)*C(2))/2.
                MP(1)=I
                MP(2)=J
                MP(3)=K
                DO 50 P=1,3
                L=MP(P)
                G(L,I)=G(L,I)+(TXX*B(P)*B(1)+TYY*C(P)*C(1))/(4.*AREA)
                G(L,J)=G(L,J)+(TXX*B(P)*B(2)+TYY*C(P)*C(2))/(4.*AREA)
                G(L,K)=G(L,K)+(TXX*B(P)*B(3)+TYY*C(P)*C(3))/(4.*AREA)
                DO 40 PP=1,3
                MPP=MP(PP)
                IF(L.EQ.MPP)GO TO 36
                D=12.0
                GO TO 38
   36   D=6.0
   38   ST(L,MPP)=ST(L,MPP)+MU*AREA/D
   40   CONTINUE
   50   CONTINUE
   100  CONTINUE
        WRITE(*,110)
   110  FORMAT(/,5X,"TIME",31X,"HEAD")
        IPRINT=2
                ICOUNT=1
                TIME=DT
                DO 500 MDT=1,NDT
                DO 150 L=1,NP
                BB(L)=0.0
                DO 150 P=1,NP
   150  BB(L)=BB(L)+ST(L,P)*H0(P)/DT

   200  DFMAX=0.0
        DO 400 L=1,NP
                IF(L.EQ.1.OR.L.EQ.2.OR.L.EQ.3.OR.L.EQ.16.OR.L.EQ.17.OR.L.EQ.18) GOTO 400
                OLDVAL=H(L)
                SUMGST=0.0
                DO 300 P=1,NP
                IF(P.EQ.L) GO TO 300
                SUMGST=SUMGST+(G(L,P)+ST(L,P)/DT)*H(P)
   300  CONTINUE
        H(L)=(BB(L)-SUMGST)/(G(L,L)+ST(L,L)/DT)
                DF=ABS(OLDVAL-H(L))
                IF(DF.GT.DFMAX) DFMAX=DF
   400  CONTINUE
        IF(DFMAX.GT.0.01) GO TO 490
                DO 450 L=1,NP
   450  H0(L)=H(L)
        IF(ICOUNT.NE.IPRINT) GO TO 490
                WRITE(*,460)TIME,(H(I),I=1,NP)
   460  FORMAT(1X,F8.2,1X,9F7.2,/,10X,9F7.2)
        ICOUNT=0
   490  TIME=TIME+DT
        ICOUNT=ICOUNT+1
   500  CONTINUE
        STOP
                END 

运行结果全部为NAN  求大神指导怎么回事？

楚香饭

请上传输入文件。

昵称

输入文件已上传。。。多谢大神指导

楚香饭

X数组前3个数为0（数据文件）
所以
C(1)=X(K)-X(J)
C(2)=X(I)-X(J)
C(3)=X(J)-X(I)
这里的 C2 和 C3 为 0

AREA=ABS(B(2)*C(3)-B(3)*C(2))/2.
所以 AREA 为 0

这里0为分母，开始出现NaN，以后就很多很多 NaN 了
G(L,I)=G(L,I)+(TXX*B(P)*B(1)+TYY*C(P)*C(1))/(4.*AREA)

昵称

chuxf 发表于 2014-6-30 14:59
X数组前3个数为0（数据文件）
所以
C(1)=X(K)-X(J)

谢谢您，我马上修改一下。真是非常感谢

昵称

chuxf 发表于 2014-6-30 14:59
X数组前3个数为0（数据文件）
所以
C(1)=X(K)-X(J)

您好。好像不是这里的原因。I,J,K分别为二位数组IJK的第E行的三个数，这三个数在X数组里的顺序不是排在一起的，例如：第一行I,J,K分别为1,2,5.能否麻烦一下再找找其他的问题？多谢

楚香饭

i j 分别为 1和2
所以x1 x2都为0
c2  c3 为零

昵称

chuxf 发表于 2014-6-30 15:26
i j 分别为 1和2
所以x1 x2都为0
c2  c3 为零

经重新核查公示，发现应为C（2）=X(I)-X(K)。这次倒是没有NAN的问题了，但是却没有输出任何数值了 不知道这个原因能否帮忙解决一下？

昵称

chuxf 发表于 2014-6-30 15:26
i j 分别为 1和2
所以x1 x2都为0
c2  c3 为零

              WRITE(*,460)TIME,(H(I),I=1,NP)
096
   460  FORMAT(1X,F8.2,1X,9F7.2,/,10X,9F7.2)
097
        ICOUNT=0
098
   490  TIME=TIME+DT
099
        ICOUNT=ICOUNT+1
就是这一部分输出不来了，请问是怎么回事呢？

楚香饭

因为 IF(ICOUNT.NE.IPRINT) GO TO 490 总是满足，于是每次都会跳到 490，就不输出了