北京化工大学2009一一2010学年第一学期 将表中所列实验数据变形后，代入多元线性回归子程序（其中调用高斯消去法子程序解正 规方程组)，通过计算确定方程中的各个系数、加、值和复相关系数.(5分) 通过计算的b、b、b:值得到安托尼方程中的参数A,B和C,A=b1,C=-b,=brAC. 《计算化学》期末考试试卷标准答案 (5 第一都分笔试都分（共30分）】 2程序框图(20分) 一、判断题，正确的在括号中划“√”，错误的划“x”(每题1.5分，共15分) 1.×2.×3.4.5.￥6.￥7.89.×10.￥ 开始 二、填空题（每空0.5分，共15分） 输入：数据点数N, 1.(1)函数在区间1ab1内连续，(2)函数在a,b两点之值Ra)与b)洞号. 2.利用矩阵的初等变换将系数矩阵A约化为一个上三角形矩阵，然后再进行求解.(1)选 苯酚的饱和蒸汽压P(I)和温度T(I)(4分) 主元，(2)消元，(3)回代。 3.函数关系，相关关系。插值法，回归分析。 ngp=y,T=x1, 1gp=x2(4分) 4.差商：中心差商. 5.x)函数形式已知，但其积分不能表示成初等函数的闭合形式：x)函数形式未知，但其 离散数据表己给出。 调用多元线性回归子程序计算回归方程的系数b、b、b:值 6.(1)用插值程序求任意点的函数值，2.用Simpson求积程序计算[a,b]区间中离散点下 及复相关系数R(其中调用高斯消去法子程序解正规方程组) (4分) 的面积。 7.实对称矩阵的本征值和本征向量：量子化学中的分子轨道相关的性质。 8。初值。边值。扩散、反应、传质、传热。 A=b1,C=b2,B=-bo+A*C (4 9.ChemOffice:Origin:Gaussian:Materials Studio. 10.线性方程组的系数矩阵为主对角线占优势。 输出：参数A,B和C(4分) 11.三角形：四面体：不用求导，甚至在没有目标函数表达式的情况下也可使用：色谱过 程的顺序优化。 结束 第二部分上机题（共70分） A 3.源程序(20分) 1.计算原理（化学原理和计算方法）(15分) program main (10分) 苯酚的饱和蕊汽压(p)和温度(T)关系遵从安托尼方程 DOUBLE PRECISION XV.S,A,X,BO,QT,U,QE,SE,R,F,X1,AA,BB,CC dimension XV(7)S(7),XI(6,7) B COMMON /AX/A(6,7),X(6,7) Ig(p'/kPa)=A- (T/K)+C OPEN(5,FILE-2010TAI.TXT,STATUS-OLD) OPEN(6,FILE-2010TA2.TXT,STATUS-UNKNOWN) READ(5,*)M,N 将两边进行变形T1gp=(AC-B)+AT-C1gp(省略单位) M1=M+1 READ(5,)((X(I,J),J=1,MI),I=1,N) 令ng,T=x,lg即=，b0=4AC-B,b=A,b=-C DO151=1N X1(LMI=X(L,I(LOGIO(X(1,2))) 得 y=bo+bx+b2xz (5分) X1L,1=XL,1) 15 X1(L,2=L0G10X1,2)
1 北京化工大学 2009——2010 学年第一学期 《计算化学》期末考试试卷标准答案 第一部分 笔试部分（共 30 分） 一、判断题，正确的在括号中划“3”，错误的划“2”（每题 1.5 分，共 15 分） 1. 2 2. 2 3. 3 4. 3 5. 2 6. 2 7. 3 8. 3 9. 2 10. 2 二、填空题（每空 0.5 分，共 15 分） 1．（1）函数在区间[a,b]内连续，（2）函数在 a、b 两点之值 f(a)与 f(b)同号。 2．利用矩阵的初等变换将系数矩阵 A 约化为一个上三角形矩阵，然后再进行求解。（1）选 主元，（2）消元，（3）回代。 3．函数关系，相关关系。插值法，回归分析。 4．差商；中心差商。 5．f(x)函数形式已知，但其积分不能表示成初等函数的闭合形式； f(x)函数形式未知，但其 离散数据表已给出。 6．（1）用插值程序求任意点的函数值，2．用 Simpson 求积程序计算[a,b]区间中离散点下 的面积。 7．实对称矩阵的本征值和本征向量；量子化学中的分子轨道相关的性质。 8．初值，边值。扩散、反应、传质、传热。 9．ChemOffice；Origin；Gaussian；Materials Studio。 10．线性方程组的系数矩阵为主对角线占优势。 11．三角形；四面体；不用求导，甚至在没有目标函数表达式的情况下也可使用；色谱过 程的顺序优化。 第二部分 上机题（共 70 分） A 1.计算原理（化学原理和计算方法）（15 分） 苯酚的饱和蒸汽压（p*）和温度（T）关系遵从安托尼方程 T C B p A + = − ( / K) lg( / kPa) * 将两边进行变形 * * T lg p = (AC − B) + AT − C lg p （省略单位） 令 Tlgp=y, T=x1, lgp=x2, b0=AC-B, b1=A, b2= -C 得 0 1 1 2 2 y = b + b x + b x （5 分） 将表中所列实验数据变形后，代入多元线性回归子程序（其中调用高斯消去法子程序解正 规方程组），通过计算确定方程中的各个系数 b0、b1、b2值和复相关系数。（5 分） 通过计算的 b0、b1、b2值得到安托尼方程中的参数 A，B 和 C，A＝b1，C＝-b2，B= b0-A*C。 （5 分） 2.程序框图（20 分） 3.源程序（20 分） program main （10 分） DOUBLE PRECISION XV,S,A,X,B0,QT,U,QE,SE,R,F,X1,AA,BB,CC dimension XV(7),S(7),X1(6,7) COMMON /AX/A(6,7),X(6,7) OPEN(5,FILE='2010TA1.TXT',STATUS='OLD') OPEN(6,FILE='2010TA2.TXT',STATUS='UNKNOWN') READ(5,*) M,N M1=M+1 READ(5,*) ((X(I,J),J=1,M1),I=1,N) DO 15 I=1,N X1(I,M1)=X(I,1)*(LOG10(X(I,2))) X1(I,1)=X(I,1) 15 X1(I,2)=LOG10(X(I,2)) 开始 输入：数据点数 N, 苯酚的饱和蒸汽压 P（I）和温度 T（I）（4 分） 调用多元线性回归子程序计算回归方程的系数 b0、b1、b2 值 及复相关系数 R（其中调用高斯消去法子程序解正规方程组） （4 分） 输出：参数 A，B 和 C （4 分） 结束 Tlgp=y, T=x1, lgp=x2（4 分） A＝b1，C＝b2，B= - b0+A*C（4 分）
D0181=1,N DO 10 K=I,N D018J=1,M1 10 XV(I)=XV(I+X(K,I) X(IJEXI() 20 XVOEXVUVDBLEONY CALLMLR(M,N,XV,S,BO,QT,U,QE,SE,R,F) D040=1,M AA=A1,MI)） S1=0 CC=-1*A2,M1) D030K=1,N BB=AA*CC-B0 30 SI=s+XK,IXV》*2 WRITE(6,20)M,N S(I)=DSQRT(S(1)) FORMAT(/4X.'INPUT DATA:,'Xnl'1X.'Xn2.....1X.'Xnm,Yn(m=1.'& D050=1,M 2m-1:2:T/7X.XU/T.IX.X2/gp.IX.Tgp-(AC+B ALIO(I)=0 +A*T-C1gp/1X,661H-V6X,XI',8X,'X2',8X,& D050K=1,N Y,7X,Ycale',7X,'DIFF/1X,66(IH-)) 50 ALIO(D=ALIO(D)+(X(K,I}XV(D))(X(K,MI)-XV(MI)) D040=1,N D060=1,MI YCALC=BO DO60 K=LN D030J=1,M 60 X(K,I(X(K,I)-XV(I)VS(I) YCALC=YCALC+A(J,MIYX(LJ) D0701=1,M DIFF=YCALC-X(I,M1) A(LI)=1.0 WRITEI6 50(X(DEL MILYCALC DIFF D070Jj=+1,M1 FORMAT(3X,I0G10.4) A(L,J)-0 WRITE(6,60)(XVI),1=1,M1) D070K=1,N 60 FORMAT(/IX,'XV,10G10.4) 70 A(IJ=A(IJ)+X(K,I)*X(K.J) WRTE(6,70)(SD,=1,M1) D080I=L,M-1 之 FORMAT(/IX,'SV,10G10.4) D080J=+1,M WRITE(6.S0) 80 AUIFA(LI 80 FORMAT(/5X,Regression equation is:'/1x,60(1H-)) WRITE(6,11)M WRITE(6,90)BO,(A(I,M1),I=1,M) 11 FORMAT(4X,'Normal cquation A(m,m+1).2x,'m=,i2/1x,60(IH-)) 90 FORMAT(IX,TIgp-,G10.4,+,G10.4,"T,+,& DO90ELM G10.4,*gp) 90 WRITE(6,33)(A(1J)JJ=I,MI) WRITE(6,100)QT,U,QE,SE,U/QT,R,F,AA,BB,CC CALL GS(A,M,MI,1.D-6) 100 FORMAT(5X'Q=,G18.6W5X,U=,G18.6/5X,Qe-,G18.6/5X& WRITEI6 22) Sc=,G18.65X,N=,G18.6,5X,R=,G18.6/5X,T=,& WRITE(6,33)(A(I,MI).I=1,M) G18.6,/5X,A=,G18.6/5X,'B=,G18.6,5X,C=,G18.6) B0=0.D0 END J=0.D0 D01001=1,M SUBROUTINE MLR(M,N,XV,S,BO,QT,U,QE,SE,R,F) (5分) A(LMIA(LMIYS(MIVS(D) DOUBLE PRECISION XV,S,A,X,BO,QT,U,QE,SE,R,F BO=B0+A(LMI)*XV(I) DIMENSION XV(M+1),S(M+1),ALIO(50) 100 【U=U+A(L,MI)ALI0I) COMMON /AX/A(6,7),X(6,7) BO=XV(MI)-BO MI=M+I QT=S(MI)S(MI) D0201=1,MI QE-QT-U XV(I)=0.DO SE-DSORT(QE/(N-M-1))
2 DO 18 I=1,N DO 18 J=1,M1 18 X(I,J)=X1(I,J) CALL MLR(M,N,XV,S,B0,QT,U,QE,SE,R,F) AA=A(1,M1) CC=-1*A(2,M1) BB=AA*CC -B0 WRITE(6,20) M,N 20 FORMAT(/4X,'INPUT DATA:','Xn1,'1X,'Xn2,....',1X,'Xnm,Yn(m=1,' & ,i2,',n=1,',i2,')'/7X,'X1/T',1X,'X2/lgp',1X,'Tlgp=(AC+B)& +A*T-Clgp'/1X,66(1H-)/6X,'X1',8X,'X2',8X, & 'Y',7X,'Ycalc',7X,'DIFF'/1X,66(1H-)) DO 40 I=1,N YCALC=B0 DO 30 J=1,M 30 YCALC=YCALC+A(J,M1)*X(I,J) DIFF=YCALC-X(I,M1) 40 WRITE(6,50) (X(I,J),J=1,M1),YCALC,DIFF 50 FORMAT(3X,10G10.4) WRITE(6,60) (XV(I),I=1,M1) 60 FORMAT(/1X,'Xv',10G10.4) WRITE(6,70) (S(I),I=1,M1) 70 FORMAT(/1X,'Sv',10G10.4) WRITE(6,80) 80 FORMAT(/5X,'Regression equation is:',/1x,60(1H-)) WRITE(6,90) B0,(A(I,M1),I=1,M) 90 FORMAT(1X,'Tlgp=',G10.4,'+',G10.4,'*T','+', & G10.4,'*lgp') WRITE(6,100) QT,U,QE,SE,U/QT,R,F,AA,BB,CC 100 FORMAT(/5X,'Qt=',G18.6,/5X,'U=',G18.6,/5X,'Qe=',G18.6,/5X, & 'Se=',G18.6,/5X,'N=',G18.6,/5X,'R=',G18.6,/5X,'F=', & G18.6,/5X,'A=',G18.6,/5X,'B=',G18.6,/5X,'C=',G18.6) END SUBROUTINE MLR(M,N,XV,S,B0,QT,U,QE,SE,R,F) （5 分） DOUBLE PRECISION XV,S,A,X,B0,QT,U,QE,SE,R,F DIMENSION XV(M+1),S(M+1),ALI0(50) COMMON /AX/A(6,7),X(6,7) M1=M+1 DO 20 I=1,M1 XV(I)=0.D0 DO 10 K=1,N 10 XV(I)=XV(I)+X(K,I) 20 XV(I)=XV(I)/DBLE(N) DO 40 I=1,M1 S(I)=0 DO 30 K=1,N 30 S(I)=S(I)+(X(K,I)-XV(I))**2 40 S(I)=DSQRT(S(I)) DO 50 I=1,M ALI0(I)=0 DO 50 K=1,N 50 ALI0(I)=ALI0(I)+(X(K,I)-XV(I))*(X(K,M1)-XV(M1)) DO 60 I=1,M1 DO 60 K=1,N 60 X(K,I)=(X(K,I)-XV(I))/S(I) DO 70 I=1,M A(I,I)=1.0 DO 70 J=I+1,M1 A(I,J)=0 DO 70 K=1,N 70 A(I,J)=A(I,J)+X(K,I)*X(K,J) DO 80 I=1,M-1 DO 80 J=I+1,M 80 A(J,I)=A(I,J) WRITE(6,11) M 11 FORMAT(4X,'Normal equation A(m,m+1)',2x,'m=',i2/1x,60(1H-)) DO 90 I=1,M 90 WRITE(6,33) (A(I,J),J=1,M1) CALL GS(A,M,M1,1.D-6) WRITE(6,22) WRITE(6,33) (A(I,M1),I=1,M) B0=0.D0 U=0.D0 DO 100 I=1,M A(I,M1)=A(I,M1)*S(M1)/S(I) B0=B0+A(I,M1)*XV(I) 100 U=U+A(I,M1)*ALI0(I) B0=XV(M1)-B0 QT=S(M1)*S(M1) QE=QT-U SE=DSQRT(QE/(N-M-1))
R=DSQRT(1.DO-QE/QT) Regression equation is (5分) F=N-M-】)LU/OEM FORMAT(/4X,'Roots of mormal equation ,5X,'A(i,m+1),i=1,m'& T1gp=2077. +6.081 *T+110.2 /1x,601H- FORMAT(1X,10G12.4) R= .999999 RETURN Fs 686329. END 6.08138 B= 1407.22 SUBROUTINE GS(A.N.M,EPS) (5分) = -110.160 (10分) DOUBLE PRECISION EPS,A DIMENSION A(6,7) 1.计算原理（化学原理和计算方法）(15分) DO 50K=LN 放射性衰变系列 P=0.0 D0101=K,N 2U4→20Th点→2Ra—→… IFOPGE ABS(AOLKINGOTO I0 P=ABS(A(LK)) 衰变方式服从连串反应动力学， L=I 其动力学方程式为一个常微分方程组： CONTINUE IF(P.LT.EPS)STOPNO SOLUTION IF(L.EQ.K)GOTO 30 d凹=-k,U叮 DO20 EK M (10分) TMP=A(L,J) A(LJA(KJ) d=k心]-kT A(K JETMP CONTINUE P=L0A(KK) 利用RngK法解此方程组。可符2过U和20Th这两种核素在10a内每释2a DO40 K+LM A(KJ)=A(K.J)*P 的质量变化 (5分) DO40 BK+IN A(1,JFA(IJ)-A(L,KA(K,J) 2.程序框图(20分》 CONTINUE 50 CONTINUE 开始 DOMK=2N I=M-K DO 60 J=I+1,N 输入：[U,Th初值234,230 A(L,M)=A(L,M)-A(LJ)*A(J,M) 时间间隔H, 计算精度PS (5分) 60 CONTINUE RETURN 计算：积分步长H0(5分) END 4.运行结果。(15分) 调用变步长Runge--Kuta法子程序解动力学方程组计算时刻t的U, [Th值 (其中调用计算动力学方程组的子程序计算DY值)(5分)
3 R=DSQRT(1.D0-QE/QT) F=(N-M-1)*U/QE/M 22 FORMAT(/4X,'Roots of mormal equation',5X,'A(i,m+1),i=1,m' & /1x,60(1H-)) 33 FORMAT(1X,10G12.4) RETURN END SUBROUTINE GS(A,N,M,EPS) （5 分） DOUBLE PRECISION EPS,A DIMENSION A(6,7) DO 50 K=1,N P=0.0 DO 10 I=K,N IF(P.GE.ABS(A(I,K)))GOTO 10 P=ABS(A(I,K)) L=I 10 CONTINUE IF(P.LT.EPS)STOP'NO SOLUTION' IF(L.EQ.K)GOTO 30 DO 20 J=K,M TMP=A(L,J) A(L,J)=A(K,J) A(K,J)=TMP 20 CONTINUE 30 P=1.0/A(K,K) DO 40 J=K+1,M A(K,J)=A(K,J)*P DO 40 I=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J) 40 CONTINUE 50 CONTINUE DO 60 K=2,N I=M-K DO 60 J=I+1,N A(I,M)=A(I,M)-A(I,J)*A(J,M) 60 CONTINUE RETURN END 4.运行结果。（15 分） Regression equation is: （5 分） ------------------------------------------------------------ Tlgp=-2077. + 6.081 *T+ 110.2 *lgp R= .999999 F= 686329. A= 6.08138 B= 1407.22 C= -110.160 （10 分） B 1.计算原理（化学原理和计算方法）（15 分） 放射性衰变系列 234 92 U⎯⎯→k1 23090Th⎯⎯→k2 22684 Ra ⎯⎯→…… 衰变方式服从连串反应动力学。 其动力学方程式为一个常微分方程组： [U] [Th] d d[Th] [U] dd[U] 1 2 1 k k t k t = − = − （10 分） 利用 Runge-Kutta 法解此方程组，可得 U 23492 和 Th 23090 这两种核素在 10 a 内每隔 2 a 的质量变化。 （5 分） 2.程序框图（20 分） 开始 输入：[U]，[Th]初值 234，230 时间间隔 H， 计算精度 EPS （5 分） 计算：积分步长 H0 （5 分） 调用变步长Runge-Kutta法子程序解动力学方程组计算时刻t的[U]， [Th]值 （其中调用计算动力学方程组的子程序计算 DY 值） （5 分）
D060=,N 60 YO(J)Y2(J) WRITE(6,22YT0,H,(Y2(I),I=1,N) 80 CONTINUE 列表输出：两种核素在10a内每隔2a的质量变化U,[Th) WRITE(6,44) (5分) 11 FORMAT(/IX,N=,I2,2X,T0=,F6.1,2X,T1=,F8.1,2X,& H0=,G8.3,1XEPs=,G10.3 结束 1X,661H-/7X,/a,5X,H/a,10X,", 12X.1h/1X.661H-》 3.源程序(20分) 22 F0RMAT1X2F8.2,2X,2G15.6) PROGRAMMAIN (10分) 33 FORMAT(/IX,'FAILED TO FIND STABLE SOLUTION IN MAX IT/) IMP鬥LICIT REAL◆8A-H.O-Z 44 FORMAT(IX66H-1 DIMENSION YO(10),YI(10),Y2(10) END OPEN(6,FILE-2010TESTB.TXT,STATUS-UNKNOWN) DATA Y0/234.0230.0.80.0/ SUBROUTINE FCYDY (5分) N=2 IMPLICIT REAL*8(A-H,O-Z) 10=0.0 REAL◆8KI,K2 T1=10.0 DIMENSION YU0)DY(10) H0=2.0 DATA K12.8E-6/,K2/8.7E-6d EPs=1.0E-6 DY(IE-KIY(1) H=Ho DY2)=k1*Y1yK2*Y(2) WRITE(6,11)N,TO,TI,HO,EPS RETURN WRITE(6,22)TO,H,(YO(I),I=1,N) END M=IDINT(CTI-TOVHO) D08010=1,M SUBROUTINE RK4(XO,YO,Y,H,N,M) (5分 K=l IMPLICIT REAL8(A-HO-Z) CALL RK4(TO,YO,Y1,HO,N,K) DIMENSION YO(N),Y(N) D030KK=1,10 DIMENSION Y2(20)DY(20)RK(20) K=K+K X=X0 H=HO/K DOSELN CALLRK4(TO,Y0,Y2,H,N,K) 5 Y(I=YO(I) ES=0.0 D050L=1,M D010J=1N CALL F(YDY) 10 ES=ES+ABS((Y2(J)-Y1(J)/Y2(J)) D0101=1,N IF(ES.LT.EPS)GOTO 50 RK(IFDY(D) D020J=1,N 10 Y2(I=Y(IH+0.5*H*DY(I) Y10=Y2U) CALL F(Y2,DY) 30 CONTINUE D0201=1,N WRITE(6,33) RK(IFRK(I)+2*DY(I) STOP 20 Y2(I=Y(IH+0.5H+DY(D) T0=T0+H0 CALL F(Y2,DY)
4 3.源程序（20 分） PROGRAM MAIN （10 分） IMPLICIT REAL *8(A-H,O-Z) DIMENSION Y0(10),Y1(10),Y2(10) OPEN(6,FILE='2010TESTB.TXT',STATUS='UNKNOWN') DATA Y0/234.0,230.0,8*0.0/ N=2 T0=0.0 T1=10.0 H0=2.0 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) DO 80 I0=1,M K=1 CALL RK4(T0,Y0,Y1,H0,N,K) DO 30 KK=1,10 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//& 1X,66(1H-)/7X,'t/a',5X,'H/a',10X,'U', & 12X,'Th' /1X,66(1H-)) 22 FORMAT(1X,2F8.2,2X,2G15.6) 33 FORMAT(/1X,'FAILED TO FIND STABLE SOLUTION IN MAX IT'/) 44 FORMAT(1X,66(1H-)) END SUBROUTINE F(Y,DY) （5 分） IMPLICIT REAL *8(A-H,O-Z) REAL *8 K1,K2 DIMENSION Y(10),DY(10) DATA K1/2.8E-6/,K2/8.7E-6/ DY(1)=-K1*Y(1) DY(2)=k1*Y(1)- K2*Y(2) RETURN END SUBROUTINE RK4(X0,Y0,Y,H,N,M) （5 分） 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) 列表输出：两种核素在 10 a 内每隔 2 a 的质量变化[U]，[Th] （5 分） 结束
D030=1,N RK(IERK(IH2.DY(I) 2. 程序框图(20分) 30 Y2(I)=Y(I)+HDY(I) CALL F(Y2.DY) 开始 D0401=1,N RK(I=RK(D+DY(D) Y(I=Y(I)+RK(I*H/6 输入：x1,x11,r1,r2 (4分)N) CONTINUE RETuRN X2=1-x1;x22=1-x11 (4分) END 4运行结果(15分，每列7.5分) 调用GAUSS消去法计算A,B(4分) N=2T0= 0T1= 10.0H0=2.00 EPS=.100E-05 常数A,B及x1=0200,x2=0.800代入原方程组计算出相应的常数 t H/a U r11,22(4分) ,00 2.00 234.000 230.000 2.00 1.00 233999 229.997 4.00 229.995 3. 结束 233.99 源程序(20分 6.00 1.00 233.996 229.992 PROGRAM MAIN (10分) 8.00 100 233.995 229.989 DIMENSION A(2 3 10.00 1.00 233.993 229.987 DOUBLE PRECISION A,EPS,X11,X21,R11,R21,X12,X22,R12,R22 INTEGER I,N,M OPEN(6,FILE-2010TESTC.Txt.STATUS-UNKNOWN) 计算原理（化学原理和计算方法）(15分) N=2 按照玛格物斯方程，二组分液态混合物中两个组分的活度系数，上与浓度，2间存在 M=N+1 如下关系： EPS=1.0D-06 X11=0.486 In=[A+2(B-A)xx (1) X21=1-X11 R11=0.813 R21=0.782 n5=[B+2(A-B)x] X12=0.200 22=1-X12 将上述方程组(1)变形为： A(1,M)=LOG(R11) A(2,M)=LOG(R21) lnr=[x(1-2x】A+2xx2B (2) A41,1=X21*2)(1-2*X11) Inr =2x,XA+[x(1-2x2)]B AM2,1=2*X21*X112 AM1,2=2*X11*X21**2 将=0.486， 1=0.813,=0.782,=1-1 代入(2)。买用高临消去法星方程阁 A(2,2=(X11*2)*(1-2X21) 可得到4，B的值.然后将x1=0.200,=1-1代入(1)，可求得相应的丙酮和氯仿的活度系数 WRITE(6,YORIGINAL LINEAR EQUATIONS: n,乃的值. WRITE(6.11(A(JEI.N+11EI.N)
5 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 4.运行结果（15 分，每列 7.5 分） N= 2 T0= .0 T1= 10.0 H0=2.00 EPS= .100E-05 ------------------------------------------------------------------ t/a H/a U Th ------------------------------------------------------------------ .00 2.00 234.000 230.000 2.00 1.00 233.999 229.997 4.00 1.00 233.997 229.995 6.00 1.00 233.996 229.992 8.00 1.00 233.995 229.989 10.00 1.00 233.993 229.987 ------------------------------------------------------------------ C 1． 计算原理（化学原理和计算方法）（15 分） 按照玛格勒斯方程，二组分液态混合物中两个组分的活度系数 r1，r2与浓度 x1，x2间存在 如下关系： ㏑ 2 1 12 r A B Ax x =+ − [ 2( ) ] （1） ㏑ 2 2 21 r B A Bx x =+ − [ 2( ) ] 将上述方程组（1）变形为： r x x A x x B r x x A x x B ln 2 [ (1 2 )] ln [ (1 2 )] 2 2 2 1 2 2 2 1 2 1 1 2 2 1 2 = + − = − + （2） 将 x1＝0.486， r1＝0.813， r2＝0.782，x2＝1-x1 代入（2），采用高斯消去法解方程组， 可得到 A，B 的值。然后将 x1＝0.200，x2＝1-x1代入（1）,可求得相应的丙酮和氯仿的活度系数 r1，r2的值。 2． 程序框图（20 分） 3． 源程序（20 分） PROGRAM MAIN （10 分） DIMENSION A(2,3) DOUBLE PRECISION A,EPS,X11,X21,R11,R21,X12,X22,R12,R22 INTEGER I,N,M OPEN(6,FILE='2010TESTC.Txt',STATUS='UNKNOWN') N=2 M=N+1 EPS=1.0D-06 X11=0.486 X21=1-X11 R11=0.813 R21=0.782 X12=0.200 X22=1-X12 A(1,M)=LOG(R11) A(2,M)=LOG(R21) A(1,1)=(X21**2)*(1-2*X11) A(2,1)=2*X21*X11**2 A(1,2)=2*X11*X21**2 A(2,2)=(X11**2)*(1-2*X21) WRITE(6,*)'ORIGINAL LINEAR EQUATIONS:' WRITE(6,11)((A(I,J),J=1,N+1),I=1,N) 开始 输入：x1,x11, r1,r2 （4 分）N） X2=1-x1; x22=1-x11 （4 分） 调用 GAUSS 消去法计算 A，B （4 分） ） 输出：r11，r22 （4 分） 结束 常数 A，B 及 x1=0.200， x2=0.800 代入原方程组计算出相应的常数 r11，r22 （4 分）
11 FORMAT(1X,3F10.4) 4 运行结果(15分) CALL GS(A,N,M,EPS) ORIGINAL LINEAR EOUATIONS (5分) R12=EXPA1.3+24AM23-AL3)}◆X12122◆*2) n74 2568 -2070 R22=EXP(A(2,3+2(AM1,3A2,3)X22)*(X12+2》 2428 -.0066 -2459 WRITE(6.33A(LMLEI.N) A=-1034D+01 (5分) 33F0RMAT5X,A=,D12.4,/5X,B=,D12.40 B=-.7764D+00 WRITE(6,44)R12,R22 WHEN X1=0.200 (5分) 44 FORMAT/5X,WHEN X1=0.200,/5X,'r1=,D12.4/5X,2=,D12.40 rl= 5511D+00 END r2=9536D+00 SUBROUTINE GS(A.N,M,EPS) (10分) D DOUBLE PRECISION A,EPS 【,计算原理（化学原理和计算方法）(15分) DIMENSION A(N,M) 某物质恒压下无相变时，其摩尔恒压热容才Cm与温度T之间的关系可表示为： D050K=1,N P=00 D010=K,N Cn小mol.K=b。+bT+b,T++h,lhT IF(P.GE.ABS(A(I,K)))GOTO 10 P=ABS(A(LK)) bo=12.4,6=2.25,b2=0.34,b3=80.6,b=24.3(5分) 10 CONTINUE 系统的拾变AH-CdT (1)(5分) IF(P.LT.EPS)STOPNO SOLUTION IF(L.EQ.K)GOTO 30 DO20 J=K,M TMP=A(LJ) 由于被积函数Ca小molK=+bT+h,T'+ T +b,lnT有确定的形式，所 A(LJ=A(KJ) A(K,J)=TMP 以可直接利用Sim即son求积分的方法计算可得到(1)式右边的积分值，其中调用自编子程序 20 CONTINUE 30 P=1.0/A(K.K) 计算函数Cp,m的位.然后代入(1)式，可得从1OK升温到50K系统的格变△H.(5分) D040j=K+1,M AKJ=AKJ◆P D0401=K+1,N 2.程序框图(20分)》 A(I,JFA(IJ)-A(LKYA(KJ) 开始 40 CONTINUE 50 CONTINUE 输入：bm=12.4,b1=2.25,b2=0.34, b=80.6,b4=24.3(5分)N D060K=2,N =MK D060J-+1,N A=10K,B=50K,(I=1,N) (5分) A(LM=A(LM)A(LJ)A(J.M) 60 CONTINUE 调用Simpson法子程序计算积分S RETURN END (用自编子程序计算函数C,m值)(5分)
6 11 FORMAT(1X,3F10.4) CALL GS(A,N,M,EPS) R12=EXP((A(1,3)+2*(A(2,3)-A(1,3))*X12)*(X22**2)) R22=EXP((A(2,3)+2*(A(1,3)-A(2,3))*X22)*(X12**2)) WRITE(6,33)(A(I,M),I=1,N) 33 FORMAT(/5X,'A=',D12.4,/5X,'B=',D12.4/) WRITE(6,44) R12,R22 44 FORMAT(/5X,'WHEN X1=0.200', /5X,'r1=',D12.4,/5X,'r2=',D12.4/) END SUBROUTINE GS(A,N,M,EPS) （10 分） DOUBLE PRECISION A,EPS DIMENSION A(N,M) DO 50 K=1,N P=0.0 DO 10 I=K,N IF(P.GE.ABS(A(I,K)))GOTO 10 P=ABS(A(I,K)) L=I 10 CONTINUE IF(P.LT.EPS)STOP'NO SOLUTION' IF(L.EQ.K)GOTO 30 DO 20 J=K,M TMP=A(L,J) A(L,J)=A(K,J) A(K,J)=TMP 20 CONTINUE 30 P=1.0/A(K,K) DO 40 J=K+1,M A(K,J)=A(K,J)*P DO 40 I=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J) 40 CONTINUE 50 CONTINUE DO 60 K=2,N I=M-K DO 60 J=I+1,N A(I,M)=A(I,M)-A(I,J)*A(J,M) 60 CONTINUE RETURN END 4． 运行结果（15 分） ORIGINAL LINEAR EQUATIONS: （5 分） .0074 .2568 -.2070 .2428 -.0066 -.2459 A= -.1034D+01 （5 分） B= -.7764D+00 WHEN X1=0.200 （5 分） r1= .5511D+00 r2= .9536D+00 D 1. 计算原理（化学原理和计算方法）（15 分） 某物质恒压下无相变时，其摩尔恒压热容才 Cp,m与温度 T 之间的关系可表示为： b T T b Cp /J mol K b b T b T ln 4 2 3 0 1 2 -1 -1 ,m ⋅ ⋅ = + + + + b0＝12.4， b1＝2.25， b2＝0.34， b3＝80.6， b4＝24.3 （5 分） 系统的焓变 H C T T T p d 21 Δ = ∫ ,m （1）（5 分） 由于被积函数 b T T b Cp /J mol K b b T b T ln 4 2 3 0 1 2 -1 -1 ,m ⋅ ⋅ = + + + + 有确定的形式，所 以可直接利用 Simpson 求积分的方法计算可得到（1）式右边的积分值，其中调用自编子程序 计算函数Cp,m 的值，然后代入（1）式，可得从 10K 升温到 50K 系统的焓变 ΔH 。（5 分） 2. 程序框图（20 分） 开始 输入：b0＝12.4， b1＝2.25， b2＝0.34， b3＝80.6， b4＝24.3（5 分）N A=10K, B=50K,（I=1，N） （5 分） 调用 Simpson 法子程序计算积分 S （用自编子程序计算函数Cp,m 值） （5 分）
IF(D.LT.EPS)GOTO 50 N-N+N 输出：从1OK升温到50K系统的焓变△H (5分) H=0.5*H T1=12 结束 S1=S2 3.源程序(20分) G0T010 PROGRAM MAIN 50 RETURN DOUBLE PRECISION EPS,TA,TB,S2 (10分) END EXTERNALF 4.运行结果(15分) OPEN(6,FILE-2010TESTD.Txr',STATUS-'UNKNOWN) DH10K->50K=20603.032227Jmol-1 EPS=1.0E-6 TA=10.0 1.计算原理（化学原理和计算方法）(15分) TB=50.0 CALL SIMP(TA,TB,EPS,S2,F) 反应： 2A—→B 10 WRITE(6.1)52 FORMAT('DH(10K->50K-,F12.6,3X,Jmol-1) 已知反应物A和B的初始浓度分别为C。,cm,A的转化率随反应时间()的变化数 END 据。 组分A的速率方程为 FUNCTION F(T) (5分 DOUBLE PRECISION T,F,BO,B1,B2,B3,B4 DATA B0W12.46/,B12.25/,B2/0.34/B3/80.6/.B4/24.3/ v,=de=k,c以 (1)(4分) F=B0+BI*T+B2*T*T+B3/T+B4*LOG(T) dt RETURN END (1)式中，VA为组分A的化学反应速率，ca为反应物A浓度，1为反应时间，n为A的 SUBROUTINE SIMP(A,B,EPS,S2,F) (5分) 反应级数，ka为A的速率常数，dC仙为反应物A浓度随时间的变化率。 DOUBLE PRECISION ABHEPS S2AI H=B-A (1)式两边取对数： T1=O.5H(F(A+F(B》 Inv =Ink +nlnc N=1 s1=0 10 S=0 y=Inv,x Inc,a=Ink,b=n (4分) D020K=1,N A1=A+K-0.5H 20 S=S+F(AI) 得y=a+br T2=0.5*(T1+H*s) S2=T2+(T2-T1)3.0 对应的反应物的浓度©A(,)可由反应关系式求出 cA()= CAd1-0A) (3分 IF(ABS(S2).LE.1.0)D=ABS(S2-S1) 将计算所得同时间1的反应物浓度ca数据通过插值和中心差商法求出反应速率V。·计算 IF(ABS(S2).GT.1.0)THEN 1nV,和1ca后，用一元线性回归子程序计算该反应的级数m,反应的表观速率常数kA·(4分) D=ABS((S2-S1)/S2) END IF
7 3. 源程序（20 分） PROGRAM MAIN DOUBLE PRECISION EPS,TA,TB,S2 （10 分） EXTERNAL F OPEN(6,FILE='2010TESTD.Txt',STATUS='UNKNOWN') EPS=1.0E-6 TA=10.0 TB=50.0 CALL SIMP(TA,TB,EPS,S2,F) 10 WRITE(6,11) S2 11 FORMAT('DH(10K-->50K)=',F12.6,3X,'J mol-1') END FUNCTION F(T) （5 分） DOUBLE PRECISION T,F,B0,B1,B2,B3,B4 DATA B0/12.46/,B1/2.25/,B2/0.34/,B3/80.6/,B4/24.3/ F=B0+B1*T+B2*T*T+B3/T+B4*LOG(T) RETURN END SUBROUTINE SIMP(A,B,EPS,S2,F) （5 分） DOUBLE PRECISION A,B,H,EPS,S2,A1 H=B-A T1=0.5*H*(F(A)+F(B)) N=1 S1=0 10 S=0 DO 20 K=1,N A1=A+(K-0.5)*H 20 S=S+F(A1) T2=0.5*(T1+H*S) S2=T2+(T2-T1)/3.0 IF(ABS(S2).LE.1.0)D=ABS(S2-S1) IF(ABS(S2).GT.1.0)THEN D=ABS((S2-S1)/S2) END IF IF(D.LT.EPS)GOTO 50 N=N+N H=0.5*H T1=T2 S1=S2 GOTO 10 50 RETURN END 4. 运行结果（15 分） DH(10K-->50K)=20603.032227 J mol-1 E 1．计算原理（化学原理和计算方法）（15 分） 反应： 2A ⎯⎯→B 已知反应物 A 和 B 的初始浓度分别为 cA0， cB0， A 的转化率随反应时间（t）的变化数 据。 组分 A 的速率方程为 n k c t c v A A A A dd = − = （1）（4 分） （1）式中， Av 为组分 A 的化学反应速率，cA为反应物 A 浓度，t 为反应时间，n 为 A 的 反应级数，kA为 A 的速率常数，dcA/dt 为反应物 A 浓度随时间的变化率。 （1）式两边取对数： A A A ln v = ln k + nln c 令 y v x c a k b n = ln , 1 = ln A , = ln , = （4 分） 得 y = a + bx 对应的反应物的浓度 cA（ti）可由反应关系式求出 cA（ti）＝ cA0(1-αA) （3 分） 将计算所得同时间 t 的反应物浓度 cA数据通过插值和中心差商法求出反应速率 Av 。计算 ln Av 和 lncA后，用一元线性回归子程序计算该反应的级数 n，反应的表观速率常数 kA。（4 分） 输出：从 10K 升温到 50K 系统的焓变ΔH （5 分） 结束
2.程序框图(20分) D0100I=1,N X(ILOG(Y(D)) 开始 YIDLOG(ABS(V(DD) 100 CONTINUE 输入：数据点数N, 反应温度TO, 时间间隔H CALL PK(N,X.YA,B,R) =B 时间t与XA的实验数据T),XA(① (2分) SS=EXP(A) WRITE(6,11)A,B,R,S,SS 1 FORMAT/5XA=,D18.10M5X'B=,D18.10,5X,R=D18.10M& CA(D)=CAO (1-xA(D))(I=1.N)( 5X,=,D18.10,5Xk=,D18.10W END X(I=T(1).Y(1)=CA(I)(I=1.N)( SUBROUTINE CF(N,H,X,Y,Z) (4分) 调用中心差商子程序计算反应速率R(wW)(dcd山r)(4分) DIMENSION X(N),Y(N),Z(N) DOUBLE PRECISION X,Y.Z,H,T,YA,YB T=X(1)+H X(I)=In[CA(I)].Y(1)=In[-R(I)](I=I,N)( CALL LGRGZ(X,YN,T,YA) T=X(1H7H 调用一元线性回归子程序计算A,B CALLLGRG2(X,YN,T,YB) Z1=4*YA-3Y(1-YB2*H (反应级数N=B,速率常数KS=EXPA) (4分)N) D0101=2N-1 T-X(D)-H 输出：A,B,RN,KS (2分) CALL LGRG2(X,Y,N,T,YA) TX(IH+H 结束 CALLLGRG2(X,YN,T,YB) Z(I(YB-YAM2+H) 10 CONTINUE 3.源程序(20分) T=X(N)-H PROGRAM MAIN (10分) CALL LGRG2(X,Y,N,T,YA) DIMENSION X(5).Y(5)V(5) T=XN-2◆H DOUBLE PRECISION X,Y,V,A,B,R,S,SS,CA,H CALLLGRG2(X,Y.N,T,YB) OPEN(6,FILE-2010TESTE.TXT,STATUS-UNKNOWN) Z(N3Y(N)-4*YA+YBM(2*H) N=5 REIΠRN 1.0E-6 END CA=2.00 DATA X00306090120/ SUBROUTINE PK(N,X,Y,A,B,R) (3分) DATA Y0,0.46,0.63,0.72,0.77/ DOUBLE PRECISION X,Y,A,B,R DO 50 I=1,N DIMENSION X(N),Y(N) Y(IFCA*(1-Y(D)) sX=0.0 0 CONTINUE SY三0.0 CALL CF(N,H,X,Y,V) SXX=0.0
8 2.程序框图（20 分） 3.源程序（20 分） PROGRAM MAIN （10 分） DIMENSION X(5),Y(5),V(5) DOUBLE PRECISION X,Y,V,A,B,R,S,SS,CA,H OPEN(6,FILE='2010TESTE.TXT',STATUS='UNKNOWN') N=5 H=1.0E-6 CA=2.00 DATA X/0.0,3.0,6.0,9.0,12.0/ DATA Y/0.,0.46,0.63,0.72,0.77/ DO 50 I=1,N Y(I)=CA*(1-Y(I)) 50 CONTINUE CALL CF(N,H,X,Y,V) DO 100 I=1,N X(I)=LOG(Y(I)) Y(I)=LOG(ABS(V(I))) 100 CONTINUE CALL PK(N,X,Y,A,B,R) S=B SS=EXP(A) WRITE(6,11)A,B,R,S,SS 11 FORMAT(/5X,'A=',D18.10//5X,'B=',D18.10, 5X,'R=',D18.10// & 5X,'n=',D18.10,5X,'k=',D18.10/) END SUBROUTINE CF(N,H,X,Y,Z) （4 分） DIMENSION X(N),Y(N),Z(N) DOUBLE PRECISION X,Y,Z,H,T,YA,YB T=X(1)+H CALL LGRG2(X,Y,N,T,YA) T=X(1)+2*H CALL LGRG2(X,Y,N,T,YB) Z(1)=(4*YA-3*Y(1)-YB)/(2*H) DO 10 I=2,N-1 T=X(I)-H CALL LGRG2(X,Y,N,T,YA) T=X(I)+H CALL LGRG2(X,Y,N,T,YB) Z(I)=(YB-YA)/(2*H) 10 CONTINUE T=X(N)-H CALL LGRG2(X,Y,N,T,YA) T=X(N)-2*H CALL LGRG2(X,Y,N,T,YB) Z(N)=(3*Y(N)-4*YA+YB)/(2*H) RETURN END SUBROUTINE PK(N,X,Y,A,B,R) （3 分） DOUBLE PRECISION X,Y,A,B,R DIMENSION X(N),Y(N) SX=0.0 SY=0.0 SXX=0.0 开始 输入：数据点数 N, 反应温度 TO, 时间间隔 H 时间 t 与 xA 的实验数据 T(I), XA(I) （I（2 分） CA(I)=CA0（1－xA(I)）（I=1,N） （2 分） 调用中心差商子程序计算反应速率 R(W)（dcA/dt） （4 分） X(I)=ln[CA(I)], Y(I)=ln[-R(I)] (I=1,N) （4 分） 调用一元线性回归子程序计算 A，B （反应级数 N＝B, 速率常数 KS=EXP(A)） （4 分）N） 输出： A，B，R, N, KS （2 分） 结束 X(I)=T(I),Y(I)=CA(I) (I=1,N) （2 分）
SYY-0.0 5Y0.0 10 IF(LABS(K-M).NE.I)THEN IF(ELTX(LTHEN YEY() SY-SY.YI K司 G0T010 SXY-SxYX1-Y1 ENDIF D-5XX-5XSX/N B-DXY/DXX R=DXY/SQRT(DXX*DYY) 00 D10-KM DO 20 J=K,M (3分) 20 END ONTINUE IF(N.LEO)RETURN Z-YGUITHEN 30 END URN IF(N.EQ2)THEN B=.1980802416D+01 R=.9706538910D+00 ,1980802416D+01 k=.1317716592D+00 且SE IF(T.GE.XN-DYTHEN
9 SYY=0.0 SXY=0.0 DO 10 I=1,N X1=X(I) Y1=Y(I) SX=SX+X1 SY=SY+Y1 SXX=SXX+X1*X1 SYY=SYY+Y1*Y1 SXY=SXY+X1*Y1 10 CONTINUE DXX=SXX-SX*SX/N DYY=SYY-SY*SY/N DXY=SXY-SX*SY/N B=DXY/DXX A=(SY-B*SX)/N R=DXY/SQRT(DXX*DYY) RETURN END SUBROUTINE LGRG2(X,Y,N,T,Z) （3 分） 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 ELSE K=1 M=N 10 IF(IABS(K-M).NE.1)THEN L=(K+M)/2 IF(T.LT.X(L))THEN M=L ELSE K=L ENDIF GOTO 10 END IF IF(ABS(T-X(K)).LT.ABS(T-X(M))) THEN K=K-1 ELSE M=M+1 END IF END IF Z=0.0 DO 30 I=K,M S=1.0 DO 20 J=K,M IF(J.NE.I)THEN S=S*(T-X(J))/(X(I)-X(J)) END IF 20 CONTINUE Z=Z+S*Y(I) 30 CONTINUE RETURN END 4. 运行结果（15 分）(每个 3 分) A= -.2026684709D+01 B= .1980802416D+01 R= .9706538910D+00 n= .1980802416D+01 k= .1317716592D+00