北京化工大学2003-2004学年第一学期 《计算化学》期末试卷答案 9 1.计算原理（化学原理和计算方法）(20分) 反应 A+B C+D B+C->D+E 其动力学方程式为一个常微分方程组： dIA]=-[A][B]+(CID] At d[B]-[A][B]+k.(CID]-&,[BJC] dIC]=IA]IB]-&(CID]-&,[B]C] d d[D]-k[A]B]-&[CID]+[BC] dt d回=kBC] 利用Run e-Kutta法解此方程组，可得不同时刻[A],B,[C],D,E值 2.程序框图 30分) 开始 输入：[A,B,[C,D,E初值 时间间隔H,计算精度EPS 计算：积分步长H0 调用变步长Runge--Kuta积分法子程序解动力学方程组计算时刻t 的[A,B,[C],D,E值 (其中调用计算动力学方程组的子程序计算DY值) 列表输出：0-200min内每隔10min所对应的[A,B,[C],D],E]值
北京化工大学 2003-2004 学年第一学期 《计算化学》期末试卷答案 A 1．计算原理（化学原理和计算方法）（20 分） 反应 k1 k2 A+B C+D B C D+E 3k + ⎯⎯→ 其动力学方程式为一个常微分方程组： 1 2 1 23 1 23 1 23 3 d[A] [A][B] [C][D] d d[B] [A][B] [C][D] [B][C] d d[C] [A][B] [C][D] [B][C] d d[D] [A][B] [C][D] [B][C] d d[E] [B][C] d k k t kkk t kkk t kkk t k t =− + =− + − = −− = −+ = 利用 Runge-Kutta 法解此方程组，可得不同时刻[A]，[B]，[C]，[D]，[E]值。 2．程序框图（30 分） 开始 输入：[A]，[B]，[C]，[D]，[E]初值 时间间隔 H， 计算精度 EPS 计算：积分步长 H0 调用变步长 Runge-Kutta 积分法子程序解动力学方程组计算时刻 t 的[A]，[B]，[C]，[D] ，[E]值 （其中调用计算动力学方程组的子程序计算 DY 值） 列表输出：0-200min 内每隔 10min 所对应的[A]，[B]，[C]，[D] ，[E]值
结柬 3.源程序(30分) PROGRAM MAIN IMPLICIT REAL*8(A-H.O-Z) DIMENSION YO(10).Y1(10),Y2(10) OPEN(6 FILE-TESTA TXT STATUS-UNKNOWNY DATA Y0/1.0.1.0.0.0.0.0.0.0,5*0.0/ T0=0.0 T1=100.0 H0=10 EPS-1.0E-6 H-HO WRITE(6.11)N.TO.TI.HO.EPS WRITE(6,22)T0,H,(Y00,I-1,N) M=IDINT((T1-TOVHO+10) D08010=1M CALLRK4(T0.Y0.Y1.H0.N.K) D030KK=1,100 K-K+K H=HO/K CALLRK4(T0.Y0.Y2.H.N.K) Es-00 D010J=1,N 10 ES=ES+ABS((Y2(-Y10VY20) IF(ES.LT.EPS)GOTO 50 D020J=1N YI()-Y2() CONTINUE WRITE(6,33) STOP 50 T0=T0+H0 D060J=1,N 60 Y0J)=Y2) WRITE(6.22)TO.H.(Y2(D)I=1.N) 80 CONTINUE WRITE(6.44) FORMATUIX N=122X TO-F6 12X TI=F8 L2X H0=G8.3,1X'EPS=,G10.320XY,e12,3,4N& 1X96(1H-)/5X,t/min,3X,'H/min,6X,'[A/molL-1',& 6X.[Bl/molL-I'.6X.'[C]/molL-I'6X.'[DVmolL-I'.6X.'[E]/molL-1'& 1X,961H-》
3．源程序（30 分） PROGRAM MAIN IMPLICIT REAL *8(A-H,O-Z) DIMENSION Y0(10),Y1(10),Y2(10) OPEN(6,FILE='TESTA.TXT',STATUS='UNKNOWN') DATA Y0/1.0,1.0,0.0,0.0,0.0,5*0.0/ N=5 T0=0.0 T1=100.0 H0=10 EPS=1.0E-6 H=H0 WRITE(6,11)N,T0,T1,H0,EPS WRITE(6,22)T0,H,(Y0(I),I=1,N) M=IDINT((T1-T0)/H0+10) DO 80 I0=1,M K=1 CALL RK4(T0,Y0,Y1,H0,N,K) DO 30 KK=1,100 K=K+K H=H0/K CALL RK4(T0,Y0,Y2,H,N,K) ES=0.0 DO 10 J=1,N 10 ES=ES+ABS((Y2(J)-Y1(J))/Y2(J)) IF(ES.LT.EPS)GOTO 50 DO 20 J=1,N 20 Y1(J)=Y2(J) 30 CONTINUE WRITE(6,33) STOP 50 T0=T0+H0 DO 60 J=1,N 60 Y0(J)=Y2(J) WRITE(6,22)T0,H,(Y2(I),I=1,N) 80 CONTINUE WRITE(6,44) 11 FORMAT(/1X,'N=',I2,2X,'T0=',F6.1,2X,'T1=',F8.1,2X, & 'H0=',G8.3,1X,'EPS=',G10.3//20X,'Y(I),I=1,2,3,4,...,N'//& 1X,96(1H-)/5X,'t/min',3X,'H/min',6X,'[A]/molL-1', & 6X,'[B]/molL-1',6X,'[C]/molL-1',6X,'[D]/molL-1',6X,'[E]/molL-1'& /1X,96(1H-)) 结束
2 F0RMAT1X2F8.2,2X5G15.6) 33 FORMAT(/IX.'FAILED TO FIND STABLE SOLUTION IN MAX IT/) ) SUBROUTINE F(Y,DY) IMPLICIT REAL8(A-H.O-ZY REAL*8K1K2K3 DIMENSION Y(10).DY(10) DATA K1/1.25E-3,K22.6E-4/,K3/8.4E-3 DY(1=-K1*Y(1)*Y(2)+K2*Y(3)Y(4) DY(2=K1*Y(1)*Y(2+K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(3)=K1*Y(1)*Y(2K2*Y(3)*Y(4)K3*Y(2)*Y(3) DY)-K1Y(Y2)-K2-Y(3)Y(4)+K3Y()Y() DY(5)=K3*Y(2)*Y(3) RETURN END DIMENSION YO(N).Y(N) DIMENSION Y2(20).DY(20),RK(20) X=XO D05=1,N 5 Y()-YOND) D050L=1,M CALL F(YDY) DO I0ILN RK(D=DY(D 10 Y2=Y+0.5*H*DY CALLF(Y2.DY DO 20 I=1.N RKT=RKI①+2◆DYD 20 Y20Y(D+0.5*HDY(I) CALLF(Y2,DY) DO30I=1.N RK=RKI+2*DY④ 30 Y2I=Y(D+H◆DYID CALL F(Y2.DY) D0401=1N RK(I)-RK(I)+DY() 40 Y(I)-Y(I)+RKHW 50 CONTINUE RETURN
22 FORMAT(1X,2F8.2,2X,5G15.6) 33 FORMAT(/1X,'FAILED TO FIND STABLE SOLUTION IN MAX IT'/) 44 FORMAT(1X,96(1H-)) END SUBROUTINE F(Y,DY) IMPLICIT REAL *8(A-H,O-Z) REAL *8 K1,K2,K3 DIMENSION Y(10),DY(10) DATA K1/1.25E-3/,K2/2.6E-4/,K3/8.4E-3/ DY(1)=-K1*Y(1)*Y(2)+K2*Y(3)*Y(4) DY(2)=-K1*Y(1)*Y(2)+K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(3)=K1*Y(1)*Y(2)-K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(4)=K1*Y(1)*Y(2) -K2*Y(3)*Y(4)+K3*Y(2)*Y(3) DY(5)=K3*Y(2)*Y(3) RETURN END SUBROUTINE RK4(X0,Y0,Y,H,N,M) IMPLICIT REAL *8(A-H,O-Z) DIMENSION Y0(N),Y(N) DIMENSION Y2(20),DY(20),RK(20) X=X0 DO 5 I=1,N 5 Y(I)=Y0(I) DO 50 L=1,M CALL F(Y,DY) DO 10 I=1,N RK(I)=DY(I) 10 Y2(I)=Y(I)+0.5*H*DY(I) CALL F(Y2,DY) DO 20 I=1,N RK(I)=RK(I)+2*DY(I) 20 Y2(I)=Y(I)+0.5*H*DY(I) CALL F(Y2,DY) DO 30 I=1,N RK(I)=RK(I)+2*DY(I) 30 Y2(I)=Y(I)+H*DY(I) CALL F(Y2,DY) DO 40 I=1,N RK(I)=RK(I)+DY(I) 40 Y(I)=Y(I)+RK(I)*H/6 50 CONTINUE RETURN
END 0T1= 100.0H0=10.0 EPS=.100E-05 Y0.=1,2,3,4N Vmin H/min ICymolL- [DVmolL-I 00 10.00 100000 100000 000000 000000 .000000 1000 987657 987154 118414E-01 12456E01 502113E03 975627 2324520E-0 262947E.01 192133E0 250 3195 022 41366 500 95250 94546 4045660 545352F-0 7039300 50.0 50 94140 .93087 .480618E-01 691260E-01 105321E-01 500 93061 91608 548602E-01 839168E-0 145283E-01 5.00 92011 13821 87118 7 288113E 10000 5.0 890355 85622 755091E0 143780 341355E-0 10.0 5.00 09 84133 793524E-01 15866 96583E-0 5.0 g634 110 电 858139E.0 57034IE4 831 2168 2990E-0 16000 500 8379s7 6900 930259E-01 231000 689868E-0 170.00 5.00 830101 5509 945967E.01 244901 750024E.0 5.00 52436 9%5449E-0 25858 810196E-D 1900 72796 79927E0 370231E.0 20.0 922E 28520 1.计算原理（化学原理和计算方法）(20分) 合成甲醇的反应是可逆 放热和物质的量减少的反应。由于反应的可逆性 通常产品产率较低。设反应过平衡时生成甲醇的物质的量为xmol,有 C0+2H2 ±CH3OH 起始物质的量mol 0
END 4．运行结果。（20 分） N= 5 T0= .0 T1= 100.0 H0=10.0 EPS= .100E-05 Y(I),I=1,2,3,4,...,N --------------------------------------------------------------------------------------------------------------------------------- t/min H/min [A]/molL-1 [B]/molL-1 [C]/molL-1 [D]/molL-1 [E]/molL-1 --------------------------------------------------------------------------------------------------------------------------------- .00 10.00 1.00000 1.00000 .000000 .000000 .000000 10.00 2.50 .987657 .987154 .118414E-01 .128456E-01 .502113E-03 20.00 2.50 .975627 .973705 .224520E-01 .262947E-01 .192133E-02 30.00 2.50 .963910 .959773 .319535E-01 .402268E-01 .413663E-02 40.00 5.00 .952504 .945465 .404566E-01 .545352E-01 .703930E-02 50.00 5.00 .941406 .930874 .480618E-01 .691260E-01 .105321E-01 60.00 5.00 .930611 .916083 .548602E-01 .839168E-01 .145283E-01 70.00 5.00 .920115 .901164 .609337E-01 .988358E-01 .189511E-01 80.00 5.00 .909911 .886179 .663565E-01 .113821 .237321E-01 90.00 5.00 .899994 .871182 .711950E-01 .128818 .288113E-01 100.00 5.00 .890355 .856220 .755091E-01 .143780 .341355E-01 110.00 5.00 .880989 .841331 .793524E-01 .158669 .396583E-01 120.00 5.00 .871888 .826549 .827730E-01 .173451 .453388E-01 130.00 5.00 .863045 .811904 .858139E-01 .188096 .511411E-01 140.00 5.00 .854452 .797418 .885139E-01 .202582 .570341E-01 150.00 5.00 .846102 .783111 .909076E-01 .216889 .629905E-01 160.00 5.00 .837987 .769000 .930259E-01 .231000 .689868E-01 170.00 5.00 .830101 .755099 .948967E-01 .244901 .750024E-01 180.00 5.00 .822436 .741416 .965449E-01 .258584 .810196E-01 190.00 5.00 .814984 .727961 .979927E-01 .272039 .870231E-01 200.00 5.00 .807740 .714740 .992602E-01 .285260 .929998E-01 --------------------------------------------------------------------------------------------------------------------------------- B 1.计算原理（化学原理和计算方法）（20 分） 合成甲醇的反应是可逆、放热和物质的量减少的反应。由于反应的可逆性， 通常产品产率较低。设反应过平衡时生成甲醇的物质的量为 x mol，有 CO+2H2 CH3OH 起始物质的量/mol CO,0 n H2 ,0 n 0
平衡物质的量/mol nco0x”,0-2x 令起始反应物的总物质的量为1mol,即：no0+m,0=1 反应达平衡时系统中总物质的量为：no0-x+nH,0-2x+x=1-2x 故平衡时系统中各物质的分压为： Pco=pYeo =p-("co 1-2x P4=m4=p-(0-3 1-2x Pcuon PYanon p-(1-2x 式中：y为物质的量分数，p为系统总压。 代入平衡关系式得： K= -22p (1) (c 1-2x 1-2x (1)式经整理后得 1+K pimo(incog K p'nco(r0 41+Kp2) 41+Kp） (2) 令S-1+K,Pa4ne+n） Q=KPrcoa( 41+Kp) 41+K,p) (2)式可写为x2-x2+5x-Q=0 用D表示原料气体中C0与H2的体积比，有 nco0+n40=1 (3) 解(3)得：aoo=D1+D) n4.0=1/1+D) 令：fx)=x3-x2+5Sx-Q 调用Newton--Raphson法求解方程f(x)=0 取初值x。=0.00
平衡物质的量/mol CO,0 n -x H2 ,0 n -2x x 令起始反应物的总物质的量为 1mol，即： 1 nCO,0 + nH2 ,0 = 反应达平衡时系统中总物质的量为： n x n 2x x 1 2x CO,0 − + H2 ,0 − + = − 故平衡时系统中各物质的分压为： ) 1 2 ( ) 1 2 ( ) 1 2 ( CH OH CH OH H ,0 H H CO,0 CO CO 3 3 2 2 2 x x p py p x n x p py p x n x p py p − = = ⋅ − − = = ⋅ − − = = ⋅ 式中：y 为物质的量分数，p 为系统总压。 代入平衡关系式得： CO,0 H ,0 2 2 2 ( ) 1 2 ( )( ) 12 12 p x p x K n x n x p p x x − = − − − − （1） （1） 式经整理后得 22 2 2 22 3 2 H ,0 CO,0 H ,0 CO,0 H ,0 2 2 1 (4 ) ( ) 0 4(1 ) 4(1 ) p p p p K pn n n K pn n xx x Kp Kp ⎡ ⎤ + + −+ − = ⎢ ⎥ + + ⎣ ⎦ （2） 令 22 2 2 2 H ,0 CO,0 H ,0 CO,0 H ,0 2 2 1 (4 ) ( ) 4(1 ) 4(1 ) p p p p K p n n n K pn n S Q Kp Kp + + = = + + （2） 式可写为 0 3 2 x − x + Sx − Q = 用 D 表示原料气体中 CO 与 H2 的体积比，有 n n D n n = + = CO,0 H ,0 CO,0 H ,0 2 2 / 1 （3） 解（3）得： 1/(1 ) /(1 ) H ,0 CO,0 2 n D n D D = + = + 令： f x = x − x + Sx − Q 3 2 ( ) 调用 Newton-Raphson 法求解方程 f (x) = 0 取初值 x0 = 0.00
fx)的导数f"(x)=3x2-2x+s 迭代格式为：=x-化山 f(x) 2.程序框图(30分) 开始☐ 输入：温度T,压力P,平衡常数KP,起始物质的量NCO,NH 体积比D,EPS,PO 计算S,Q 调用Newton-Raphson法子程序求解方程f(x)=0 (其中调用子程序计算f(x)和f'(x)值)得x 计算平衡时各组分的摩尔分数YC,YH,YP 输出：T,P,D,x,YC,YH,YP值 结束 3.源程序(30分) PROGRAM MAIN REALTP.KPNC.NH.D.EPS.S.Q.X0.YC.YH.YPX EXTERNAL F.G COMMON S.O OPEN(6,FILE-TESTB.TXT,STATUS-UNKNOWN) DATA TP/365.3E07/ DATA EPS/IE-6/ WRITEyPlease input volum rate VCO:VH2 READ(,)D T-273.15+T KP=3925/T-9.84*L0G10T-0.00347*T+9.8
f (x) 的导数 f ′(x) = 3x − 2x + S 2 迭代格式为： ( ) ( ) 1 1 1 − − − ′ = − n n n n f x f x x x 2.程序框图（30 分） 3.源程序（30 分） PROGRAM MAIN REAL T,P,KP,NC,NH,D,EPS,S,Q,X0,YC,YH,YP,X EXTERNAL F,G COMMON S,Q OPEN(6,FILE='TESTB.TXT',STATUS='UNKNOWN') DATA T,P/365,3E07/ DATA EPS/1E-6/ WRITE(*,*)'Please input volum rate VCO:VH2:' READ(*,*)D T=273.15+T KP=3925/T-9.84*LOG10(T)-0.00347*T+9.8 开始 输入：温度 T，压力 P，平衡常数 KP，起始物质的量 NCO，NH， 体积比 D，EPS，P0 调用 Newton-Raphson 法子程序求解方程 f (x) = 0 （其中调用子程序计算 f (x) 和 f ′(x) 值）得 x 输出：T，P，D，x，YC，YH，YP 值 结束 计算 S，Q 计算平衡时各组分的摩尔分数 YC，YH，YP
KP=10**KP NC=D/(+D) **2*NH*(4*NC+NH)(4*(1+KP*P**2)) Q=KP*P**2*NC*NH*NH/(4*(1+KP*P**2) X0=0.0 CALL NEWN(XO EPS EGX) YC=ONC.XVL-2*X) YH-(NH-2*X)1-2*X YP=X1-2◆X) WRITE(6.11)T.P.D.x.YC,YH,YP FORMAT(/5X.T-.F10.4.K/5X.P-.G12.4.Pa/5X.VCO:VH2-.F12.4& /5X,X=,F6.4 RESULT:/5X,yCO)尸，F6.4& /5X.y(H2)-.F6.4/5X.y(CH3OH)-.F6.4) END FUNCTION F(X) REALXSO COMMON S O F-X3-X2*X-Q RETURN END FUNCTION G(X) REALX.S.Q COMMON S,Q G=3*X**2-2*X+S RETURN ENd SUBROUTINE NEWN(X0.EPS.F.GX) REAL XO,EPS K=1 10 FX=F(XO) GX-G(X0 X=X0-FX/GX IF((ABS(X-X0)/X).GT.EPS)THEN XO=X K=K+1 IF(K.GT.100)THEN WRITE(,YNO RESOLUTION G0T030 END IF
KP=10**KP NC=D/(1+D) NH=1/(1+D) S=(1+KP*P**2*NH*(4*NC+NH))/(4*(1+KP*P**2)) Q=(KP*P**2*NC*NH*NH)/(4*(1+KP*P**2)) X0=0.0 CALL NEWN(X0,EPS,F,G,X) YC=(NC-X)/(1-2*X) YH=(NH-2*X)/(1-2*X) YP=X/(1-2*X) WRITE(6,11)T,P,D,x,YC,YH,YP 11 FORMAT(/5X,'T=',F10.4,'K'/5X,'P=',G12.4,'Pa'/5X,'VCO:VH2=',F12.4 & /5X,'X=',F6.4//'RESULT:'/5X,'y(CO)=',F6.4 & /5X,'y(H2)=',F6.4 /5X,'y(CH3OH)=',F6.4 ) END FUNCTION F(X) REAL X,S,Q COMMON S,Q F=X**3-X**2+S*X-Q RETURN END FUNCTION G(X) REAL X,S,Q COMMON S,Q G=3*X**2-2*X+S RETURN END SUBROUTINE NEWN(X0,EPS,F,G,X) REAL X0,EPS K=1 10 FX=F(X0) GX=G(X0) X=X0-FX/GX IF((ABS(X-X0)/X).GT.EPS)THEN X0=X K=K+1 IF(K.GT.100) THEN WRITE(*,*)'NO RESOLUTION' GOTO 30 END IF
G0T010 ENDIE 30 RETURN END 4.运行结果。(20分) T=638.1500K 3000E+08 VCO:VH2- 1.5000 X=.1431 RESULT: yC0F.6401 y(H2=.1595 y(CH30H=.2004 1.计算原理（化学原理和计算方法）(20分) 1mol物质在标准压力下(p)从温度T变化至2的熵变化为： △S=S喷-号，若在该温度范围内无相变化，则 等喷空n (1) 本题中，已知周态铅在15K时的标准熵15K)=24JK1ml1,由(1) 式有： (2) 此题中，已知恒压摩尔热容C,a与温度的关系：C~T 可推导出二、T关系，可根据离散点求积的方法即：用S1m即sOm求积方法 T (其中调用Lagrange插值法子程序计算任意点的函数值)计算(2)式右端积分 项，然后通过(2)式求出500K时固态铅的标准痛(500K)
GOTO 10 ENDIF 30 RETURN END 4.运行结果。（20 分） T= 638.1500K P= .3000E+08Pa VCO:VH2= 1.5000 X= .1431 RESULT: y(CO)= .6401 y(H2)= .1595 y(CH3OH)= .2004 C 1.计算原理（化学原理和计算方法）（20 分） 1mol 物质在标准压力下（ ）从温度 T1 变化至 T2 的熵变化为： ，若在该温度范围内无相变化，则 （1） 本题中，已知固态铅在 15K 时的标准熵 ，由（1） 式有： （2） 此题中，已知恒压摩尔热容Cp,m 与温度的关系：Cp,m ~ T 可推导出 T T Cp m ~ , 关系，可根据离散点求积的方法即：用 Simpson 求积方法 （其中调用 Lagrange 插值法子程序计算任意点的函数值）计算（2）式右端积分 项，然后通过（2）式求出 500K 时固态铅的标准熵
2.程序框图(30分) 开始 输入：数据点数N,EPS,S1, T(I),Cp(D)(I=1,N) N) X(I=T(D):Y()-Cp(IT()(I=IN):A=TIB=T2 调用Simpson法子程序计算AB面积S (用插值法计算F),S2-S1+S 输出：S2 结束☐ 3.源程序(30分) PROGRAM MAIN DOUBLE PRECISION EPS.A.B.X.Y.S.KP1.KP DIMENSION X(9).Y(9) N=9 EPS=1 0E.6 OPEN(6.FILE-TESTC.Txl'.STATUS-UNKNOWN) DATA X/15,30,50,100,150,200,30,400,550/ DATA Y/7.2,16.4,21.1,24.3,252,25.8,26.5,27.4,28.9/ A=15 B=500 61=24 D010=1,N X=x 10 Y(I)-Y(I/XI CALLSIMP(N.A.B,X,Y,EPS.S) 92=S+91 WRITE(6 20)A BS WRITE(6.0)B.S2 20 FORMAT(IX.'S(.F6.2.K->.F6.2.K)-.G12.6) 30 FORMAT(1X,'S,F7.2,'K)尸，G12.6,'JK-1mol-1D END SUBROUTINESIMP(M.A.B.X.Y.EPS.S2) DOUBLEPRECISIONA.B.X.Y.EPS.S2.FA.FB.X0.Y0 DIMENSION X(M),Y(M) H=B-A
2.程序框图（30 分） 3.源程序（30 分） PROGRAM MAIN DOUBLE PRECISION EPS,A,B,X,Y,S,KP1,KP2 DIMENSION X(9),Y(9) N=9 EPS=1.0E-6 OPEN(6,FILE='TESTC.Txt',STATUS='UNKNOWN') DATA X/15,30,50,100,150,200,300,400,550/ DATA Y/7.2,16.4,21.1,24.3,25.2,25.8,26.5,27.4,28.9/ A=15 B=500 S1=2.4 DO 10 I=1,N XI=X(I) 10 Y(I)=Y(I)/XI CALL SIMP(N,A,B,X,Y,EPS,S) S2=S+S1 WRITE(6,20) A,B,S WRITE(6,30) B,S2 20 FORMAT(1X,'S(',F6.2,'K-->',F6.2,'K)=',G12.6) 30 FORMAT(1X,'S(',F7.2,'K)=',G12.6,'JK-1mol-1') END SUBROUTINE SIMP(M,A,B,X,Y,EPS,S2) DOUBLE PRECISION A,B,X,Y,EPS,S2,FA,FB,X0,Y0 DIMENSION X(M),Y(M) H=B-A 开始 输入：数据点数 N,EPS,S1 ， T(I),Cp(I) （I=1，N） N） X(I)=T(I)；Y(I)=Cp(I)/T(I) (I=1,N)；A=T1,B=T2 调用 Simpson 法子程序计算 A-B 面积 S （用插值法计算 F），S2=S1+S 输出：S2 结束
N=1 CALLLGRG2(X.Y.M.A.FA) CALLLGRG2(X.Y.M.B.FB) T1=0.5*H(FA+FB) 10 S=0 DO20K=IN X0=A+K-0.5)*H CALLLGRG2(X.Y.M.X0.Y0) 20 S-S+Y0 T2-0.5*(T1+H*S S2=T2+T2.T1/3.0 IFON EO UGOTO 40 30 D=ABS(S2-S1) IF(ABS(S2).GT.1.0)THEN D=ABS((S2-S1VS2) EnD IF IF(D.LT.EPS)GOTO 50 40 N=N+N H=0.5*H T1=2 S1=S2 G0T010 50 Return END SUBROUTINE LGRG2(X,Y.N,T.Z) DIMENSION XON)YON) DOUBLE PRECISION X.Y.T.Z.S Z=00 IF(N.LE.0)RETURN IF(N.EQ.I)THEN Z=Y(1) RETURN END IF IF(N.EQ.2)THEN Z-(Y1)*(T-X2Y2)*(T-X1)MX1-X2) RETURN END IF IFT.LE X(2)THEN K=1 M-3 ELSE IF(T.GE.X(N-1))THEN K=N-2 M-N
N=1 CALL LGRG2(X,Y,M,A,FA) CALL LGRG2(X,Y,M,B,FB) T1=0.5*H*(FA+FB) 10 S=0 DO 20 K=1,N X0=A+(K-0.5)*H CALL LGRG2(X,Y,M,X0,Y0) 20 S=S+Y0 T2=0.5*(T1+H*S) S2=T2+(T2-T1)/3.0 IF(N.EQ.1)GOTO 40 30 D=ABS(S2-S1) IF(ABS(S2).GT.1.0)THEN D=ABS((S2-S1)/S2) END IF IF(D.LT.EPS)GOTO 50 40 N=N+N H=0.5*H T1=T2 S1=S2 GOTO 10 50 RETURN END SUBROUTINE LGRG2(X,Y,N,T,Z) DIMENSION X(N),Y(N) DOUBLE PRECISION X,Y,T,Z,S Z=0.0 IF(N.LE.0)RETURN IF(N.EQ.1)THEN Z=Y(1) RETURN END IF IF(N.EQ.2)THEN Z=(Y(1)*(T-X(2))-Y(2)*(T-X(1)))/(X(1)-X(2)) RETURN END IF IF(T.LE.X(2))THEN K=1 M=3 ELSE IF(T.GE.X(N-1))THEN K=N-2 M=N