实验数据的数值积分问题

实验数据的数值积分问题

用不同温度下固体铅的定压热容数据，求298K 时的绝对熵。数据见下表：

由热力学方法可知，绝对熵与定压热容有如下关系：

dT T

C S T

P T ⎰

=

采用MA TLAB 来计算，这里将最小二乘B 样条拟合法和三次样条插值函数计算作比较，命令行程序如下：

function Entropy_Int clear all clc

% 读入数据

T = [5. 10. 15. 20. 25. 30. 50. 70. 100. 150. 200. 250. 298.];

Cp = [0.305 2.8 7.0 10.8 14.1 16.5 21.4 23.3 24.5 25.4 25.8 26.2 26.5]; CpT = Cp./T;

T = [0 T]; CpT = [0 CpT];

% 用最小二乘B 样条拟合法和三次样条插值函数计算 knots = 3; K = 3; % 三次B 样条 sp = spap2(knots,K,T,CpT);

cs = csapi(T,CpT); % 生成三次样条插值函数cs % 绘制浓度拟合曲线 Ti = linspace(T(1),T(end),299); CpT_B = fnval(sp,Ti);

CpT_c = fnval(cs,Ti);

plot(T,CpT,'ro',Ti,CpT_B,'b-',Ti,CpT_c,'Bl-.') xlabel('T')

ylabel('Cp/T')

legend('实验值','B 样条拟合','立方样条拟合') % 进行数值积分 pp = fnint(cs); s = fnval(pp,T)

T

C p /T

s =

0 0.0433 0.8733 2.7820 5.3352 8.1104 10.9071 20.7236 28.2599 36.8085 46.9207 54.2919 60.0880 64.7233

在解决本问题时，需首先要考虑积分下限，即0=T 时的情形，因此时的T C P 为0/0。为此，可考虑先求出接近0K 时的P C 数值。可有多种方法此数值，如通过P C T ~图的样条插值曲线读出，或采用回归拟合以及插值计算的方法。这里采用上述方法得到的数值是

05.0=T ，003.0=P C 。由于本题的温度数据间隔不等，直接采用样条插值积分是一较好

的方法。运用IMSL 数学库中的样条插值CSINT 和依据样条插值系数的积分CSITG ，可以求得绝对熵，下面给出调用的主程序和计算结果。

C To obtain the absolute entropy of solid lead at 298 K from

C the experimental data of the heat capacity at constant pressure CP of C solid lead at various temperatures up to and including 298 K

C by the cubic spline interpolant and integration from IMSL library.

C the absolute entropy in CRC Handbook of Chemistry and Physics is 64.8 J K-1 mol-1. INTEGER NDATA, I, NINTV, NOUT PARAMETER (NDATA=14) REAL A, B, BREAK(NDATA), CSCOEF(4,NDATA), CSITG, ERROR,

& EXACT, F, FDATA(NDATA), FI, FLOAT, VALUE, X, XDATA(NDATA) INTRINSIC FLOAT EXTERNAL CSINT, CSITG, UMACH DATA (XDATA(I),I=1,NDATA)/0.05,5.,10.,15.,20.,25.,30.,50., & 70.,100.,150.,200.,250.,298./ DATA (FDATA(I),I=1,NDATA)/0.003,0.305,2.8,7.0,10.8,14.1,16.5, & 21.4,23.3,24.5,25.4,25.8,26.2,26.5/ C Compute cubic spline interpolant DO 10 I=1,NDATA FDATA(I)=FDATA(I)/XDATA(I) 10 CONTINUE

C cubic spline interpolant

CALL CSINT (NDATA, XDATA, FDATA, BREAK, CSCOEF) C cubic spline integration A = 0.0 B = 298. NINTV = NDATA - 1 VALUE = CSITG(A,B,NINTV,BREAK,CSCOEF) EXACT = 64.8 ERROR = EXACT - VALUE C Get output unit number CALL UMACH (2, NOUT) C Print the result WRITE (NOUT,99999) A, B, VALUE, EXACT, ERROR 99999 FORMAT (' On the closed interval (', F4.1, ',', F5.1, & ') we have :', /, 1X, 'Computed Integral = ', F10.5, /, & 1X, 'Exact Integral = ', F7.2, /, 1X, 'Error ' & , ' = ', F10.5, /, /) END

Output

On the closed interval (0.0,298.0) we have : Computed Integral = 64.82892 Exact Integral = 64.8 Error = -0.02892

查阅CRC 化学和物理手册，可得到固体铅298K 时的绝对熵为64.81

1--⋅⋅mol K J ，与计算

结果之间的差别很小。

本题可用多种方法来解决，如Simpson 数值积分法，或采用⎰

∞

-=

T

P T T d C S ln )ln (的形式，或先

用数据点拟合方程得回归参数，然后对回归方程进行积分（回归方程可用y C P =，

，也可在此基础上探索更好的回归方程形式）。

拟合方程 >> % 读入数据

>> T = [5. 10. 15. 20. 25. 30. 50. 70. 100. 150. 200. 250. 298.];

>> Cp = [0.305 2.8 7.0 10.8 14.1 16.5 21.4 23.3 24.5 25.4 25.8 26.2 26.5]; >> figure(1); >> subplot(2,1,1); >> plot(T,Cp,'o');

>> aa=polyfit(T,Cp,3); >> f_x=vpa(poly2sym(aa,'x'),2) f_x =

.56e-5*x^3-.31e-2*x^2+.52*x+.35