帮忙看看这个程序？ c c c DIMENSION R(10),RS(10),QS(10),XL(10),KG(10),NTEXT(40),NUM(10),IPAR 1(4),PAR(4),DEV(100,2),JENS(10),NU(100,2) DIMENSION X(5,4),F(5),XB(4),XS(4),XM(4),XE(4),XX(4),XR(4),XK(4) COMMON T(100),NM(100,2),XXX(100,2),GME(100,2),GMC(100,2),GMR(100,2 1),NNY(10,10),Q(10),A(10,10) DATA NC,NP/5,6/ READ(NC,16)NCOMP,NG,NOBS,LAAF,NOIT 16 FORMAT(20I3) C C NCOMP=NUMBER OF COMPONENTS C NG=NUMBER OF DIFFERENT GROUPS C NOBS=NUMBEI OF DATA POINTS C LAAF,LAAF=1 THE EXPERIMENTAL ACTIVITY COEFFICIENTS ARE READ C LAAF=2 THE LOGARITHMS TO THE EXPERIMENTAL ACTIVITY COEFFICIENTS C ARE READ C NOIT,NOIT=1 PARAMETER ESTIMATION IS PERFORMED C NOIT=2 NO PARAMETER ESTIMATION.THE ACTIVITY C COEFFICIENTS ARE CALCULATED BASED ON C THE GIVEN R,Q AND A MATRIX IF(NOIT.EQ.2) GO TO 1301 WRITE(NP,15) 15 FORMAT(1H1,' PARAMETER ESTIMATION') GO TO 1302 1301 WRITE(NP,1303) 1303 FORMAT(1H1,' CALCULATION OF THE ACTIVITY COEFFICIENTS BASED ON 1THE GIVEN R,O AND A MATRIX') 1302 CONTINUE DO 14 I=1,16 READ(NC,13)NTEXT 14 WRITE(NP,13)NTEXT 13 FORMAT(40A2) READ(NC,16)(NUM(I),I=1,NCOMP) C C NUM GIVEN THE NUMBERS ATTANHED TO C THE DIFFERENT COMPONENTS C DO 1304 I=1,NCOMP 1304 READ(NC,2)(NNY(I,K),K=1,NG) 2 FORMAT(40I2) C C NNY(I,K) IS THE MATRIX GIVING THE NUMBER OF GROUPS OF C KIND K IN MOLECULE IC C READ(NC,16)(KG(K),K=1,NG) C C KG(K) IS THE NUMBER ATTACHED TO GROUP K C DO 100 N=1,NOBS 100 READ(NC,101)T(N),(NU(N,I),XXX(N,I),GME(N,I),I=1,2) 101 FORMAT (F7.2,3(I2,2F7.4)) C C N=DATA POINT NUMBER C T(N)=TRMPERATURE IN K C NU(N,I)=NUMBER ATTACHED TO COMPONENT I C (NU(N,I)=NUM(I)) C XXX(N,I)=LIQUID MOLE FRACTION C GME(N,I),LAAF=1 XPERIMENTAL ACTIVITY COEFFICIENT C LAAF=2 LOGARITHMS TO THE EXPERIMENTAL ACTIVITY COEFFICIENT C IF(LAAF-1)820,820,821 821 CONTINUE DO 822 N=1,NOBS DO 822 I=1,2 GGG=GME(N,I) GME(N,I)=EXP(GGG) 822 CONTINUE 820 CONTINUE DO 705 I=1,NG READ(NC,704) R(I),Q(I),(A(I,J),J=1,NG) 704 FORMAT(8F10.4) C C R(I) IS THE GROUP VOLUME OF GROUP I C Q(I) IS THE GROUP AREA OF GROUP I C A(I,J) IS THE GROUP INTERACTION PARAMETER C BETWEEN GROUPS I AND J C 705 CONTINUE DO 7 I=1,NCOMP RS(I)=0. QS(I)=0. DO 8 J=1,NG RS(I)=RS(I)+NNY(I,J)*R(J) QS(I)=QS(I)+NNY(I,J)*Q(J) 8 CONTINUE 7 XL(I)=5.*(RS(I)-QS(I))-RS(I)+1. WRITE(NP,5) 5 FORMAT(1HO,'COMPONENT/GROUPS',/) WRITE(NP,12)(KG(I),I=1,NG) 12 FORMAT(17X,10I3) DO 10 I=1,NCOMP 10 WRITE(NP,6)NUM(I),(NNY(I,J),J=1,NG) 6 FORMAT(10X,I3,4X,10I3) WRITE(NP,3) 3 FORMAT(1HO, 'GROUP NO GROUP R GROUP Q',/) DO 11 I=1,NG 11 WRITE(NP,4)KG(I),R(I),Q(I) 4 FORMAT(4X,I3,5X,2E12.4) WRITE(NP,1400) 1400 FORMAT(1HO,' GROUP INTERACTION PARAMETERS',/) DO 1401 I=1,NG 1401 WRITE(NP,1402)(A(I,J),J=1,NG) 1402 FORMAT(10E12.4) DO 102 N=1,NOBS DO 102 I=1,2 DO 104 J=1,NCOMP IF(NU(N,I)-NUM(J))104,103,104 103 NM(N,I)=J GO TO 102 104 CONTINUE 102 CONTINUE CALL PFAC3(RS,QS,XL,NOBS) IF(NOIT-1)830,830,831 831 CONTINUE CALL PFAC4(NG,NOBS) DO 832 NR=1,NOBS DO 832 I=1,2 GMR(NR,I)=GMC(NR,I)+GMR(NR,I) GMR(NR,I)=EXP(GMR(NR,I)) 832 CONTINUE GO TO 833 830 CONTINUE READ(NC,16)NPAR,KRIT,IDEN C C NPAR IS THE NUMBER OF PARAMETERS TO C BE ESTIMATED C C KRIT DETERMINES THE OBJECTIVE FUNCTION C IE KRIT=1,THE SUM OF THE SQUARED DIFFERENCES BETWEEN C THE EXPERIMENTAL AND CALCULATED AVTIVITY COEFFICIENTS C IS MINIMIZED C IE KRIT=2,THE LOFARITHMS OF THE ACTIVITY COEFFICIENTS C ARE USED C C IDEN IS THE NUMBER OF IDENTICAL PAIRS C OE INTERANTION PARAMETERS C READ(NC,16)(IPAR(I),I=1,NPAR) C C IPAQ IS THE VECTOR INDICATING THE C PARAMETERS TO BE ESTIMATED IN THE C A MATRIX C IF(IDEN)950,950,951 951 IDEN=2*IDEN READ(NC,16)(JENS(J),J=1,IDEN) 950 CONTINUE C C JENS IS THE VECTOR INDICATING THE IDENTICAL C PAIR OF PARAMETERS C NN=NPAR+1 N=NPAR READ(NC,400)(X(I,I),I=1,NPAR) C C X(I,I)=ROW OF INITIAL PARAMETERS C NO INITIAL PARAMETER MUST BE ZERO C SA=1.E-6 C C SA IS THE STANDARD ERROR AS DEFINED BY C NELDER-MEAD C DO 201 J=2,NN DO 201 I=1,N IF(J-I-1)202,203,202 203 X(J,I)=1.1*X(1,I) GO TO 201 202 X(J,I)=X(1,I) 201 CONTINUE WRITE(NP,300)(IPAR(I),I=1,N) 300 FORMAT(1HO,' INITIAL PARAMETERS',413/) DO 204 J=1,NN 204 WRITE(NP,400)(X(J,I),I=1,N) DO 1 J=1,NN DO 21 I=1,N 21 XX(I)=X(J,I) CALL FMIN(NPAR,IPAR,PAR,NOBS,NG,XX,FF,KRIT,JENS,IDEN) 1 F(J)=FF NF=NN C C NF IS THE NUMBER OF CALCULATION OF F C ALFA=1. BETA=1.0/2.0 GAMMA=2. ITER=0 JPR=0 400 FORMAT(8F10.3) C C ESTIMATION OF THE LOWEST VALUE OF F=FB C 25 FB=F(1) DO 98 I=1,N 98 XB(I)=X(1,I) JB=1 DO 31 J=2,NN IF(FB-F(J))31,31,108 108 FB=F(J) JB=1 DO 41 I=1,N 41 XB(I)=X(J,I) 31 CONTINUE C C ESTINATION OF THE HEGHEST VALUE OF F=FS C FS=F(1) DO 51 I=1,N 51 XS(I)=X(1,I) JS=1 DO 61 J=2,NN IF (FS-F(J))111,61,61 111 FS=F(J) JS=J DO 71 I=1,N 71 XS(I)=X(J,I) 61 CONTINUE C C CALCULATION OF THE CENTROID XM(I) OF POINTS C EXCLUDING XS(I) C DO 81 I=1,N 81 XM(I)=-XS(I) DO 9 J=1,NN DO 122 I=1,N 122 XM(I)=XM(I)+X(J,I) 9 CONTINUE DO 121 I=1,N 121 XM(I)=XM(I)/FLOAT(N) C C REFLECTION C DO 131 I=1,N 131 XR(I)=XM(I)+ALFA*(XM(I)-XS(I)) CALL FMIN(NPAR,IPAR,PAR,NOBS,NG,XR,FR,KRIT,JENS,IDEN) NF=NF+1 C C EXPANSION C IF(FR-FB)141,151,151 141 DO 161 I=1,N 161 XE(I)=XM(I)+GAMMA*(XR(I)-XM(I)) CALL FMIN(NPAR,IPAR,PAR,NOBS,NG,XE,FE,KRIT,JENS,IDEN) NF=NF+1 IF (FE-FB)17,18,18 17 DO 19 I=1,N X(JS,I)=XE(I) 19 XS(I)=XE(I) F(JS)=FE C C CALCULATION OF THE HALTING CRITERION C 27 FM=0. DO 20 J=1,NN 20 FM=FM+F(J) FM=FM/FLOAT(NN) FRMS=0. DO 22 J=1,NN 22 FRMS=(F(J)-FM)**2+FRMS RMS=SQRT(FRMS/FLOAT(N)) ITER=ITER+1 JPR=JPR+1 IF(ITER-200) 500,500,23 500 CONTINUE IF(JPR-1)902,902,903 903 CONTINUE IF(JPR-6)901,904,904 904 JPR=1 902 CONTINUE WRITE(NP,107)ITER,NF 107 FORMAT(1HO,'ITERATION',I4,' NUMBER OF CALLS FOR THE SUBROUTIN 1E',I5) WRITE(NP,109) 109 FORMAT(' PARAMETERS') WRITE(NP,400)(X(JS,I),I=1,N) WRITE(NP,106)F(JS),RMS 106 FORMAT(1H,' FMIN=',E14.5,' SD=',E14.5) 901 CONTINUE IF(RMS-SA)23,23,25 C C NEW SIMPLEX C FE GREATER THAN FB C 18 DO 26 I=1,N X(JS,I)=XR(I) 26 XS(I)=XR(I) F(JS)=FR FS=FR GO TO 27 C C NEW SIMPLEX C FR GREATER THAN FB C 151 DO 30 J=1,NN IF(J-JS)28,30,28 28 IF(FR-F(J))18,18,30 30 CONTINUE IF(FR-FS)91,91,32 91 DO 33 I=1,N X(JS,I)=XR(I) 33 XS(I)=XR(I) F(JS)=FR FS=FR 32 DO 34 I=1,N 34 XK(I)=XM(I)+BETA*(XS(I)-XM(I)) CALL FMIN(NPAR,IPAR,PAR,NOBS,NG,XK,FK,KRIT,JENS,IDEN) NF=NF+1 C C NEW SIMPLEX C AFTER CONTRACTION C IF(FK-FS)35,35,36 35 DO 37 I=1,N X(JS,I)=XK(I) 37 XS(I)=XK(I) F(JS)=FK FS=FK GO TO 27 36 DO 38 J=1,NN DO 39 I=1,N 39 X(J,I)=(X(J,I)+XB(I))/2. 38 CONTINUE GO TO 27 23 WRITE(NP,905) 905 FORMAT(1HO,' FINAL PARAMETERS') WRITE(NP,400)(X(JS,I),I=1,N) WRITE(NP,106)F(JS),RMS 833 CONTINUE DO 906 N=1,NOBS DO 906 I=1,2 DEV(N,I)=(GMR(N,I)-GME(N,I))*100./GME(N,I) 906 CONTINUE WRITE(NP,207) 207 FORMAT(1HO,' TEMP NUMBER X GAMEXP GAMCAL DEV 1',/) DO 222 N=1,NOBS WRITE(NP,16)N DO 223 I=1,2 223 WRITE(NP,96)T(N),NU(N,I),XXX(N,I),GME(N,I),GMR(N,I),DEV(N,I) 96 FORMAT(F8.2,I5,F8.4,2F14.4,F6.1) 222 CONTINUE IF(NOIT-1)835,835,836 835 CONTINUE WRITE(NP,1400) DO 112 I=1,NG 112 WRITE(NP,1402)(A(I,J),J=1,NG) WRITE(NP,996) 996 FORMAT(1HO,' THE SUM OF THE SQUARED DIFFERENCES BETWEEN THE') IF(KRIT-1)997,997,998 998 WRITE(NP,995) 995 FORMAT(' LOGARITHMS TO HE') 997 WRITE(NP,999) 999 FORMAT(' EXPERIMENTAD AND CALCULATED ACTIVITY COEFFICIENTS IS 1MINIMIZED') 836 CONTINUE STOP END c 计算活度系数的组合项 SUBROUTINE PFAC3(RS,QS,XL,NOBS) c DIMENSION THETA(2),PHI(2),RS(10),QS(10),XL(10) COMMON T(100),NM(100,2),XXX(100,2),GME(100,2),GMC(100,2),GMR(2 1),NNY(10,10),Q(10),A(10,10) DO 3 N=1,NOBS SQ=0. SR=0. SXL=0. DO 2 I=1,2 J=NM(N,I) SXL=SXL+XL(J)*XXX(N,I) SQ=SQ+QS(J)*XXX(N,I) 2 SR=SR+RS(J)*XXX(N,I) DO 3 I=1,2 J=NM(N,I) THETA(I)=QS(J)/SQ PHI(I)=RS(J)/SR GMC(N,I)=ALOG(PHI(I))+5.*QS(J)*ALOG(THETA(I)/PHI(I))+XL(J)-PHI(I)* 1SXL 3 CONTINUE RETURN END c 计算活度系数剩余项 c SUBROUTINE PFAC4(NG,NOBS) C DIMENSION GMOL(10),ATET(10),ANYK(10),BNYK(10),GK(10,3),P(10,10) COMMON T(100),NM(100,2),XXX(100,2),GME(100,2),GMC(100,2),GMR(100, 12),NNY(10,10),Q(10),A(10,10) C C CALCULATION OF THE PSI MATRIX DO 250 NR=1,NOBS DO 7 I=1,NG DO 7 J=1,NG 7 P(I,J)=EXP(-A(I,J)/T(NR)) C C CALCULATION OF GROUP MOLE FRACTIONS C DO 105 II=1,3 IF(II-2)100,100,101 100 SNYK=0. C C PURE COMPONENT C J=NM(NR,II) DO 12 K=1,NG 12 SNYK=SNYK+FLOAT(NNY(J,K)) DO 13 K=1,NG 13 GMOL(K)=FLOAT(NNY(J,K))/SNYK GO TO 102 101 SNYK=0. C C MIXITURE C DO 2 I=1,2 J=NM(NR,I) DO 2 K=1,NG 2 SNYK=SNYK+FLOAT(NNY(J,K))*XXX(NR,I) DO 3 K=1,NG GNYK=0. DO 4 I=1,2 J=NM(NR,I) 4 GNYK=GNYK+FLOAT(NNY(J,K))*XXX(NR,I) 3 GMOL(K)=GNYK/SNYK C C CALCULATION OF GROUP AREA FRACTIONS C 102 SNYK=0. DO 5 K=1,NG 5 SNYK=SNYK+Q(K)*GMOL(K) DO 6 K=1,NG 6 ATET(K)=Q(K)*GMOL(K)/SNYK C C CALCULATION OF GAMMA K C DO 9 K=1,NG ANYK(K)=0. DO 10 M=1,NG SNYK=0. DO 8 N=1,NG 8 SNYK=SNYK+ATET(N)*P(N,M) 10 ANYK(K)=ATET(M)*P(K,M)/SNYK+ANYK(K) BNYK(K)=0. DO 11 M=1,NG 11 BNYK(K)=BNYK(K)+ATET(M)*P(M,K) BNYK(K)=ALOG(BNYK(K)) 9 GK(K,II)=Q(K)*(1.-BNYK(K)-ANYK(K)) 105 CONTINUE DO 201 I=1,2 J=NM(NR,I) SNYK=0. DO 200 K=1,NG 200 SNYK=SNYK+FLOAT(NNY(J,K))*(GK(K,3)-GK(K,I)) 201 GMR(NR,I)=SNYK 250 CONTINUE RETURN END c 一组参数函数的目标函数 SUBROUTINE FMIN(NPAR,IPAR,PAR,NOBS,NG,XX,FF,KRIT,JENS,IDEN) C DIMENSION IPAR(4),PAR(4),XX(4),JENS(10) COMMON T(100),NM(100,2),XXX(100,2),GME(100,2),GMC(100,2),GMR(100,2 1),NNY(10,10),Q(10),A(10,10) DO 2 I=1,NPAR,2 KI=IPAR(I) KJ=IPAR(I+1) A(KI,KJ)=XX(I) 2 A(KJ,KI)=XX(I+1) IF (IDEN)9,9,8 8 KKI=JENS(1) KKJ=JENS(2) DO 7 J=3,IDEN,2 IKI=JENS(J) IKJ=JENS(J+1) A(IKI,IKJ)=A(KKI,KKJ) A(IKJ,IKI)=A(KKJ,KKI) 7 CONTINUE 9 CONTINUE CALL PFAC4(NG,NOBS) DO 200 NR=1,NOBS DO 200 I=1,2 GMR(NR,I)=GMC(NR,I)+GMR(NR,I) GMR(NR,I)=EXP(GMR(NR,I)) 200 CONTINUE FF=0. DO 3 N=1,NOBS DO 3 I=1,2 IF (KRIT-1)10,10,20 10 FF=FF+(GMR(N,I)-GME(N,I))**2 GO TO 3 20 GCAL=GMR(N,I) GEXP=GME(N,I) FF=FF+(ALOG(GCAL)-ALOG(GEXP))**2 3 CONTINUE RETURN END查看更多4个回答 . 3人已关注
第1天
关注盖德视界
添加小助手
