run-time error M6021:MATH - log:DOMAIN error

clchien · 发表于 2014-5-4 00:59:29

以下的程式碼執行後會出現以下錯誤訊息
取對數時的定義域有問題, 煩請各位先進指點迷津, 小弟感恩不盡 !!!
run-time error M6021:MATH
- log:DOMAIN error

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

001     PROGRAM NGDE
002      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
003 
004      PARAMETER (MAX=42)
005      PARAMETER (MAX1=42)
006 
007        DIMENSION V(MAX1+1),XN(MAX1),X(MAX1,MAX1,MAX1),COLF(MAX1,MAX1), &
008       XN1S(MAX1),ZETA(MAX1),DP(MAX1),XM(MAX1)
009 
010      PI=3.1415926
011      XKB=1.38D-23
012      R=8.314
013      XNA=6.02225D23
014 
015      T=1773.D0
016      COOLRATE=1000.D0
017      P=1.
018      XMW=26.95D-3
019      RHO=2700.D0
020      A=948.D0
021      B=0.202
022      C=13.07
023      D=36373.D0
024      AMU=1.966D-6
025      BMU=147.47
026 
027      Q=10.**(12./(MAX-2))
028      V1=XMW/(RHO*XNA)
029 
030      V(1)=V1
031 
032      DO 100 I=1,MAX-1
033          XN(I)=0.
034          XN1S(I)=0.
035100     CONTINUE
036 
037      DO 101 I=2,MAX
038          V(I)=V(1)*(Q**(I-1))
039 
040 
041          DP(I)=(6.*V(I)/PI)**0.333333
042          XM(I)=RHO*V(I)
043101     CONTINUE
044 
045      DP(1)=(6.*V(1)/PI)**0.333333
046        XM(1)=RHO*V(1)
047 
048        DO 102 K=2,MAX-1
049        DO 102 I=2,MAX-1
050        DO 102 J=2,MAX-1
051 
052          IF (V(K) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K+1)) THEN
053                   X(I,J,K)=(V(K+1)-V(I)-V(J))/(V(K+1)-V(K))
054 
055        ELSE IF (V(K-1) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K)) THEN
056                   X(I,J,K)=(V(I)+V(J)-V(K-1))/(V(K)-V(K-1))
057 
058                ELSE
059                   X(I,J,K)=0.D0
060          END IF
061102     CONTINUE
062 
063      TIME=0.D0
064      STEP=1.D-4
065      S=1.001
066      TEMP1=6.*V1/PI
067      S1=PI*(TEMP1**0.66666666667)
068      XMASS=RHO*V1
069      PS=EXP(C-D/T)*101325.*P
070      XNS=PS/(XKB*T)
071      XN(1)=S*XNS
072 
073 
074200   FORMAT(  'Time',5X,'Temperature',5X,'XN(1)',5X,'XJK',5X,'S',5X,'Diameter',5X,'Particle volume',5X,'Total number')
075      WRITE(*,200)
076 
077      ICOUNTER=1
078      ICOUNTER2=1
079 
080        DO WHILE (T .GT. 1600)
081        SIGMA=(A-B*T)*1.D-3
082        PS=EXP(C-D/T)*101325.*P
083        XNS=PS/(XKB*T)
084        S=XN(1)/XNS
085        THETA=(S1*SIGMA)/(XKB*T)
086        AA=(2.*SIGMA)/(PI*XMASS)
087        BB=THETA-(4.*(THETA**3.))/(27.*(DLOG(DABS(S))**2.))
088        XJK=(XNS**2.)*S*V1*(AA**0.5)*EXP(BB)
089        CC=0.6666667*THETA/DLOG(DABS(S))
090        XKSTAR=CC**3
091        DPSTAR=4.*SIGMA*V1/(XKB*T*DLOG(DABS(S)))
092        VSTAR=PI*(DPSTAR**3.)/6.
093 
094        XNTOT=0.D0
095          DO I=2,MAX-1
096            XNTOT=XNTOT+XN(I)
097          END DO
098 
099        VTOT=0.D0
100        VTOTAL=0.D0
101          DO I=2,MAX-1
102            VTOT=VTOT+XN(I)*V(I)
103 
104          END DO
105            VTOTAL=VTOT+XN(1)*V(1)
106 
107        VAV=VTOT/XNTOT
108        DPAV=(6.*VAV/PI)**0.3333333
109 
110          IF (XNTOT .GT. 100.) STEP=1.D-5
111        IF (TIME .GT. 0.14) STEP=5.D-8
112 
113201     FORMAT(8(1X,E14.8))
114202     FORMAT(8(1X,E14.8))
115203     FORMAT(/'time=',E14.8)
116204     FORMAT('N',I2,2(2X,E14.8))
117        OPEN(11,FILE='output1.DAT',STATUS='UNKNOWN',ACCESS='APPEND')
118        OPEN(12,FILE='output2.DAT',STATUS='UNKNOWN',ACCESS='APPEND')
119 
120 
121        ICOUNTER=ICOUNTER+1
122        ICOUNTER2=ICOUNTER2+1
123 
124          IF (ICOUNTER .EQ. 400) THEN
125                WRITE(*,201)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
126          WRITE(11,202)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
127          ICOUNTER=1
128        END IF
129 
130          IF (ICOUNTER2 .EQ. 40000) THEN
131          WRITE(*,203)TIME
132          WRITE(12,203)TIME
133 
134           DO I=1,MAX-1
135             WRITE(*,204)I,DP(I),XN(I)
136             WRITE(12,204)I,DP(I),XN(I)
137           END DO
138          ICOUNTER2=1
139        END IF
140 
141        XKMIN=1.D-9
142        XMU=AMU*(T**1.5)/(BMU+T)
143        XLAMBDA=(XMU/(P*101325.))*SQRT(PI*R*T/(2.*0.04))
144        DO 103 I=1,MAX-1
145          DO 103 J=1,MAX-1
146          XKN1=(2.*XLAMBDA)/DP(I)
147          XKN2=(2.*XLAMBDA)/DP(J)
148          D1=(XKB*T)/(3.*PI*XMU*DP(I))*((5.+4.*XKN1+6.*(XKN1**2)+18.*(XKN1**3))/(5.-XKN1+(8.+PI)*(XKN1**2)))
149          D2=(XKB*T)/(3.*PI*XMU*DP(J))*((5.+4.*XKN2+6.*(XKN2**2)+18.*(XKN2**3))/(5.-XKN2+(8.+PI)*(XKN2**2)))
150          C1=SQRT((8.*XKB*T)/(PI*XM(I)))
151          C2=SQRT((8.*XKB*T)/(PI*XM(J)))
152          XL1=(8.*D1)/(PI*C1)
153          XL2=(8.*D2)/(PI*C2)
154          G1=((DP(I)+XL1)**3)-((DP(I)**2+XL1**2)**1.5)/(3.*DP(I)*XL1)-DP(I)
155          G2=((DP(J)+XL2)**3)-((DP(J)**2+XL2**2)**1.5)/(3.*DP(J)*XL2)-DP(J)
156          COLF(I,J)=2.*PI*(D1+D2)*(DP(I)+DP(J))/((DP(I)+DP(J))/(DP(I)+DP(J)+2.*SQRT(G1**2+G2**2))+(8.*(D1+D2))/(SQRT(C1**2+C2**2)*(DP(I)+DP(J))))
157 
158 
159          IF (COLF(I,J) .LT. XKMIN) THEN
160            XKMIN=COLF(I,J)
161          END IF
162 
163103       CONTINUE
164 
165 
166        DO 104 K=2,MAX-1
167          IF (VSTAR .LT. V1) THEN
168               ZETA(2)=VSTAR/V(2)
169              ELSE IF (V(K-1) .LE. VSTAR .AND. VSTAR .LT. V(K)) THEN
170               ZETA(K)=VSTAR/V(K)
171            ELSE
172               ZETA(K)=0.D0
173             END IF
174104       CONTINUE
175 
176 
177        DO 105 I=1,MAX-1
178          DP(I)=(6.*V(I)/PI)**0.333333
179          XN1S(I)=XNS*EXP(4.*SIGMA*XMW/(R*T*RHO*DP(I)))
180105       CONTINUE
181 
182        T1=0.D0
183        T2=0.D0
184        T3=0.D0
185        T4=0.D0
186 
187        DO 106 K=2,MAX-1
188          SUM1=0.D0
189          SUM2=0.D0
190 
191            IF (STEP .NE. 1.D-4) THEN
192            IF (XN(1) .GT. XN1S(K-1)) THEN
193 
194              IF (K .EQ. 2) THEN
195                ADDTERM=0.D0
196                ELSE
197                ADDTERM=(V(1)/(V(K)-V(K-1)))*COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
198                T1=T1+COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
199              END IF
200              END IF
201 
202            IF (XN(1) .LT. XN1S(K+1)) THEN
203                ADDTERM=-(V(1)/(V(K+1)-V(K)))*COLF(1,K+1)*(XN(1)-XN1S(K+1))*XN(K+1)
204                T2=T2+COLF(1,K+1)*(-XN(1)+XN1S(K+1))*XN(K+1)
205              END IF
206 
207            IF (XN(1) .LT. XN1S(K)) THEN
208                SUBTERM=-(V(1)/(V(K)-V(K-1)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
209                T3=T3+COLF(1,K)*(-XN(1)+XN1S(K))*XN(K)
210              END IF
211 
212            IF (XN(1) .GT. XN1S(K)) THEN
213                SUBTERM=(V(1)/(V(K+1)-V(K)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
214                T4=T4+COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
215              END IF
216          END IF
217          DO I=2,MAX-1
218              SUM2=SUM2+COLF(K,I)*XN(I)
219             DO J=2,K
220               SUM1=SUM1+X(I,J,K)*COLF(I,J)*XN(I)*XN(J)
221             END DO
222          END DO
223          XN(K)=XN(K)+STEP*(0.5*SUM1-XN(K)*SUM2+XJK*ZETA(K)+ADDTERM-SUBTERM)
224 
225106       CONTINUE
226        XN(1)=XN(1)-XJK*XKSTAR*STEP-STEP*(T1+T4-T3-T2)
227        T=T-STEP*COOLRATE
228        TIME=TIME+STEP
229 
230 
231        CLOSE(11)
232          CLOSE(12)
233 
234      END DO
235 
236        STOP
237        END

楚香饭 · 发表于 2014-5-4 05:54:38

楼主，你好。

建议长代码直接上传 *.f90 或 *.for 或 *.f 源代码文件。这样大家好告诉你是第几行的问题。（复制粘贴后，有时候行数不好对应）

根据我的调试，你的问题是，主程序中，Addterm 没有初值。

在第191行，IF (STEP .NE. 1.D-4) THEN 中，STEP ==1.D-4，所以 NE 条件不满足，直接跳到 ENDIF 执行，即217行。虽然在195行对Addterm有给初值，但没有执行。

217行的循环执行完毕后，执行223行的
XN(K)=XN(K)+STEP*(0.5*SUM1-XN(K)*SUM2+XJK*ZETA(K)+ADDTERM-SUBTERM)
时，Addterm 没有初值，于是出现了极其不正常的数据。然后到后面就错了。

clchien · 发表于 2014-5-12 18:06:22

非常感恩chuxf大人的協助,

目前程式已經可以跑了, 我還在繼續測試結果是否合理

[数值问题] run-time error M6021:MATH - log:DOMAIN error

新人勋章

爱心勋章

水王勋章

元老勋章

热心勋章