Fortran Coder

查看: 19401|回复: 14
打印 上一主题 下一主题

[数值问题] 程序运行出现NAN

[复制链接]

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
跳转到指定楼层
楼主
发表于 2014-6-30 14:36:16 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
自己编的程序代码如下:      
[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  求大神指导怎么回事?
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

736

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
700 元
贡献
359 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

沙发
发表于 2014-6-30 14:42:39 | 只看该作者
请上传输入文件。

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
板凳
 楼主| 发表于 2014-6-30 14:46:34 | 只看该作者
输入文件已上传。。。多谢大神指导

XZB.txt

85 Bytes, 下载次数: 1

YZB.txt

82 Bytes, 下载次数: 1

IJKDY.txt

251 Bytes, 下载次数: 1

736

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
700 元
贡献
359 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

地板
发表于 2014-6-30 14:59:11 | 只看该作者
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)

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
5#
 楼主| 发表于 2014-6-30 15:00:38 | 只看该作者
chuxf 发表于 2014-6-30 14:59
X数组前3个数为0(数据文件)
所以
C(1)=X(K)-X(J)

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

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
6#
 楼主| 发表于 2014-6-30 15:08:28 | 只看该作者
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.能否麻烦一下再找找其他的问题?多谢

736

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
700 元
贡献
359 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

7#
发表于 2014-6-30 15:26:12 | 只看该作者
本帖最后由 chuxf 于 2014-6-30 15:29 编辑

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

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
8#
 楼主| 发表于 2014-6-30 16:05:10 | 只看该作者
chuxf 发表于 2014-6-30 15:26
i j 分别为 1和2
所以x1 x2都为0
c2  c3 为零

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

9

帖子

1

主题

0

精华

入门

F 币
34 元
贡献
20 点
9#
 楼主| 发表于 2014-6-30 16:15:08 | 只看该作者
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
就是这一部分输出不来了,请问是怎么回事呢?

736

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
700 元
贡献
359 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

10#
发表于 2014-6-30 16:39:39 | 只看该作者
因为 IF(ICOUNT.NE.IPRINT) GO TO 490 总是满足,于是每次都会跳到 490,就不输出了
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-12-23 21:52

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表