南京工业大学研究生工程应用数学B下第二次作业

专业：安全科学与工程姓名：孙皓学号：652083700020

1、（1）将无限稀释活度系数的估算过程用MATLAB实现，要求将 lnγ1∕x22、-lnγ2∕x12和ln（γ1∕γ2）与x1的关系图画出来。

答：首先将已知列入表1-1，数据的第一行和最后一行包含需要求解的未知数，暂不考虑，将其它数据录入为huoduxishu.m文件。

表1-1 X1 0 0.0819 0.2192 0.3584 0.3831 0.5256 0.8478 0.9872 1 Y1 0 0.1869 0.4065 0.5509 0.5748 0.6786 0.8741 0.9863 1 P（mmHg） 178.08 202.74 236.86 266.04 270.73 293.36 324.66 327.39 327.05 γ1 ？ 1.415 1.343 1.25 1.242 1.158 1.023 1 1 γ2 1 1.008 1.011 1.046 1.048 1.116 1.508 1.968 ？ huoduxishu.m文件内容如下：

X1=[0.0819;0.2192;0.3584;0.3831;0.5256;0.8478;0.9872]; X2=1-X1;

Y1=[0.1869;0.4065;0.5509;0.5748;0.6786;0.8741;0.9863]; P=[202.74;236.86;266.04;270.73;293.36;324.66;327.39]; gama1=[1.415;1.343;1.25;1.242;1.158;1.023;1];

gama2=[1.008;1.011;1.046;1.048;1.116;1.508;1.968]; lng1g2=log(gama1./gama2); lng1x22=log(gama1)./(X2.^2);

lng2x12=((-1)*log(gama2))./(X1.^2);

在matlab命令窗中输入如下：

huoduxishu

plot(X1,lng1g2,'k*-');

hold on

plot(X1,lng2x12,'b*-'); hold on

plot(X1,lng1x22,'r*-'); hold off

legend ('lng1g2','lng1x12','lng1x22')

得到图1-1，

图1-1

lnγ1∕x22曲线和-lnγ2∕x12曲线分别与ln（γ1∕γ2）曲线外推在x1=0和1处有交点。

由图可见，lnγ1∕x22曲线的后半段和-lnγ2∕x12曲线的前半段分别对其拟合计算可能有一定影响，考虑分别取前半段与后半段（4或5个点）进行拟合计算。

（1）编辑huoduxishux1.m文件内容如下：

X1=[0.3831;0.5256;0.8478;0.9872]; gama1=[1.242;1.158;1.023;1]; gama2=[1.048;1.116;1.508;1.968]; lng1g2=log(gama1./gama2);

lng2x12=((-1)*log(gama2))./(X1.^2);

在matlab命令窗中输入如下：

clear all

>> huoduxishux1 plot(X1,lng1g2,'k*-'); hold on

plot(X1,lng2x12,'b*-'); hold off

legend ('lng1g2','lng2x12')

得到图1-2，

图1-2

运行：[P,S]=polyfit(X1,lng2x12,1) 得：P = -0.6032 -0.0820 S = R: [2x2 double] df: 2 normr: 0.0286

由此得，-lnγ2∕x12 的拟合方程为Y= -0.6032X-0.0820。X1=1时，Y=-0.6852，

输入：Y=-0.6852;

gama2wuqiongda=exp(-Y)， 得到：gama2wuqiongda = 1.9842

（2）编辑huoduxishux2.m文件内容如下：

X1=[0.0819;0.2192;0.3584;0.3831]; X2=1-X1;

gama1=[1.415;1.343;1.25;1.242]; gama2=[1.008;1.011;1.046;1.048]; lng1g2=log(gama1./gama2); lng1x22=log(gama1)./(X2.^2);

在matlab命令窗中输入如下：

clear all

>> huoduxishux2 plot(X1,lng1g2,'k*-'); hold on

plot(X1,lng1x22,'r*-'); hold off

legend ('lng1g2','lng1x22')

得到图1-3，

图1-3

运行：[P,S]=polyfit(X1,lng1x22,1) 得：P = 0.4991 0.3717

S = R: [2x2 double]

df: 2 normr: 0.0111

由此得，-lnγ1∕x22的拟合方程为Y= 0.4991X+0.3717。X1=0时，Y=0.3717，

输入：Y=0.3717;

gama1wuqiongda=exp(Y)，

得到：gama1wuqiongda = 1.4502。 综上，γ

∞

1

=1.4502，γ

∞

2

=1.9842。

（3）编辑huoduxishuxhuatu.m文件内容如下：

X1=[0;0.0819;0.2192;0.3584;0.3831;0.5256;0.8478;0.9872;1]; X2=1-X1;

gama1=[1.4502;1.415;1.343;1.25;1.242;1.158;1.023;1;1]; gama2=[1;1.008;1.011;1.046;1.048;1.116;1.508;1.968;1.9842]; lng1g2=log(gama1./gama2) lng1x22=log(gama1)./(X2.^2) lng2x12=((-1)*log(gama2))./(X1.^2)

在matlab命令窗中输入如下：

clear all

huoduxishuhuatu

运行得到：

lng1g2 = 0.3717 0.3392 0.2840 0.1782 0.1698 0.0369

-0.3880 -0.6770 -0.6852

lng1x22 = 0.3717 0.4118 0.4837 0.5421 0.5695 0.6518 0.9816 0 NaN

lng2x12 = NaN -1.1879 -0.2277 -0.3501 -0.3194 -0.3973 -0.5715 -0.6947 -0.6852

再运行：

plot(X1,lng1g2,'k*-'); hold on

plot(X1,lng2x12,'b*-'); hold on

plot(X1,lng1x22,'r*-'); hold off

legend ('lng1g2','lng1x12','lng1x22')

得到图1-4，即lnγ1∕x22、-lnγ2∕x12和ln(γ1∕γ2)与X1的关系图。

图1-4lnγ1∕x22、-lnγ2∕x12和ln(γ1∕γ2)与X1的关系 根据以上结果完成表1-1并添加lnγ1∕x22、-lnγ2∕x12和ln(γ∕γ2)的数据得表1-2。

表1-2 X1 0 Y1 0 P（mmHg） 178.08 202.74 236.86 266.04 270.73 293.36 γ1 1.4502 1.415 1.343 1.25 1.242 1.158 γ2 1 1.008 1.011 1.046 1.048 1.116 ln(γ1∕γ2) 0.3717 0.3392 0.2840 0.1782 0.1698 0.0369 lnγ1∕x22 0.3717 0.4118 0.4837 0.5421 0.5695 0.6518 -lnγ2∕x12 -1.1879 -0.2277 -0.3501 -0.3194 -0.3973 1

0.0819 0.1869 0.2192 0.4065 0.3584 0.5509 0.3831 0.5748 0.5256 0.6786 X1 Y1 P（mmHg） 324.66 327.39 327.05 γ1 1.023 1 1 γ2 1.508 1.968 1.9842 ln(γ1∕γ2) -0.3880 -0.6770 -0.6852 lnγ1∕x22 0.9816 0 -lnγ2∕x12 -0.5715 -0.6947 -0.6852 0.8478 0.8741 0.9872 0.9863 1 1

2、分别采用梯形法求积分进行热力学一致性检验，并将结果与辛普森积分法进行比较。 答：

在matlab命令窗中输入如下：

clear all clc

X1=[0;0.0819;0.2192;0.3584;0.3831;0.5256;0.8478;0.9872;1]; gama1=[1.4502;1.415;1.343;1.25;1.242;1.158;1.023;1;1]; gama2=[1;1.008;1.011;1.046;1.048;1.116;1.508;1.968;1.9842]; lng1g2=(-1).*log(gama1./gama2)

得： lng1g2 = -0.3717 -0.3392 -0.2840 -0.1782 -0.1698 -0.0369 0.3880 0.6770

0.6852

继续输入：

p = polyfit(X1,lng1g2, 2)

得到：

p = 0.8868 0.1613 -0.3629

即二次多项式拟合得到的方程为：

Y=0.8868.*X1.^2+0.1613.*X1-0.3629

输入：

X = 0:.001:1;

Y = 0.8868.*X.^2+0.1613.*X-0.3629; T = trapz(X,Y)

得到：T = 0.0134 输入：

X = 0:.001:1;

Y = 0.8868.*X.^2+0.1613.*X-0.3629; Y=abs(Y); T2 = trapz(X,Y)

得到：T2 = 0.2654

一致性检验：0.0134÷0.2654=0.05048，合5.0%。其值比原方法大，但两结果同样是大于2%，误差主要是由于采用多项式拟合方程导致。

3、在积分方法统一采用采用辛普森法下，采用线性插值描述ln(γ

1

∕γ2)与的x1的关系，并将热力学一致性检验结果与样条插值法进行比较。

答：编写jianyan.m程序，内容如下：

clear all clc

X=[0;0.0819];

Y=[0.3717;0.3392]; X1=[0:.01: 0.0819];

Y1=interp1(X,Y,X1, 'linear')

X=[0.0819;0.2192]; Y=[0.3392;0.2840];

X2=[0.0819:.01:0.2192]; Y2=interp1(X,Y,X2, 'linear') X=[0.2192;0.3584]; Y=[0.2840;0.1782];

X3=[0.2192:.01:0.3584]; Y3=interp1(X,Y,X3, 'linear') X=[0.3584;0.3831]; Y=[0.1782;0.1698];

X4=[0.3584:.01:0.3831]; Y4=interp1(X,Y,X4, 'linear') X=[0.3831;0.5256]; Y=[0.1698;0.0369];

X5=[0.3831:.01:0.5256]; Y5=interp1(X,Y,X5, 'linear') X=[0.5256;0.8478]; Y=[0.0369;-0.3880]; X6=[0.5256:.01:0.8478]; Y6=interp1(X,Y,X6, 'linear') X=[0.8478;0.9872]; Y=[-0.3880;-0.6770]; X7=[0.8478:.01:0.9872]; Y7=interp1(X,Y,X7, 'linear') X=[0.9872;1];

Y=[-0.6770;-0.6852]; X8=[0.9872:.01:1];

Y8=interp1(X,Y,X8, 'linear')

运行：jianyan 得到：

Y1 = 0.3717 0.3677 0.3479 0.3439 0.3400

Y2 = 0.3392 0.3352 0.3151 0.3111 0.3070 0.2869

Y3 = 0.2840 0.27 0.2384 0.2308 0.2232 0.1852

0.3638 0.3312 0.3030 0.2688 0.2156 0.3598 0.3271 0.2990 0.2612 0.2080 0.3558 0.3519 0.3231 0.3191 0.2950 0.2910 0.2536 0.2460 0.2004 0.1928 Y4 = 0.1782 0.1748 0.1714

Y5 = Columns 1 through 14 0.1698 0.1605 0.1511 0.1418 0.1325 0.1232 0.1138 0.1045 0.0952 0.0859 0.0765 0.0672 0.0579 0.0486 Column 15 0.0392

Y6 = Columns 1 through 14 0.0369 0.0237 0.0105 -0.0027 -0.0158 -0.0290 -0.0422 -0.0554 -0.0686 -0.0818 -0.0950 -0.1082 -0.1213 -0.1345

Columns 15 through 28 -0.1477 -0.1609 -0.1741 -0.1873 -0.2005 -0.2137 -0.2268 -0.2400 -0.2532 -0.26 -0.2796 -0.2928 -0.3060 -0.3192

Columns 29 through 33 -0.3323 -0.3455 -0.3587 -0.3719 -0.3851 Y7 = -0.3880 -0.4087 -0.4295 -0.4502 -0.4709 -0.4917 -0.5124 -0.5331 -0.5539 -0.5746 -0.5953 -0.6160 -0.6368 -0.6575 Y8 = -0.6770 -0.6834

输入：

X=[0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800];

Y=[0.3717 0.3677 0.3638 0.3598 0.3558 0.3519 0.3479 0.3439 0.3400]; p1 = polyfit(X, Y, 1)

X=[0.0819 0.0919 0.1019 0.1119 0.1219 0.1319 0.1419 0.1519 0.1619 0.1719 0.1819 0.1919 0.2019 0.2119];

Y=[0.3392 0.3352 0.3312 0.3271 0.3231 0.3191 0.3151 0.3111 0.3070 0.3030 0.2990 0.2950 0.2910 0.2869]; p2 = polyfit(X, Y, 1)

X=[0.2192 0.2292 0.2392 0.2492 0.2592 0.2692 0.2792 0.22 0.2992 0.3092 0.3192 0.3292 0.3392 0.3492];

Y=[0.2840 0.27 0.2688 0.2612 0.2536 0.2460 0.2384 0.2308 0.2232 0.2156 0.2080 0.2004 0.1928 0.1852];

p3 = polyfit(X, Y, 1) X=[0.36 0.37 0.38]; Y=[0.17820.17480.1714]; p4 = polyfit(X, Y, 1)

X=[0.3831 0.3931 0.4031 0.4131 0.4231 0.4331 0.4431 0.4531 0.4631 0.4731 0.4831 0.4931 0.5031 0.5131 0.5231];

Y=[0.16980.16050.15110.14180.13250.12320.11380.10450.09520.08590.07650.06720.05790.04860.0392]; p5 = polyfit(X, Y, 1)

X=[0.5256 0.5356 0.5456 0.5556 0.5656 0.5756 0.5856 0.5956 0.6056 0.6156 0.6256 0.6356 0.56 0.6556 0.6656 0.6756 0.6856 0.6956 0.7056 0.7156 0.7256 0.7356 0.7456 0.7556 0.7656 0.7756 0.7856 0.7956 0.8056 0.8156 0.8256 0.8356 0.8456];

Y=[0.03690.02370.0105-0.0027-0.0158-0.0290-0.0422-0.0554-0.0686-0.0818-0.0950-0.1082-0.1213-0.1345-0.1477-0.1609-0.1741-0.1873-0.2005-0.2137-0.2268-0.2400-0.2532-0.26-0.2796-0.2928-0.3060-0.3192-0.3323-0.3455-0.3587-0.3719-0.3851]; p6 = polyfit(X, Y, 1)

X=[0.8478 0.8578 0.8678 0.8778 0.8878 0.78 0.9078 0.9178 0.9278 0.9378 0.9478 0.9578 0.9678 0.9778];

Y=[-0.3880-0.4087-0.4295-0.4502-0.4709-0.4917-0.5124-0.5331-0.5539-0.5746-0.5953-0.6160-0.6368-0.6575]; p7 = polyfit(X, Y, 1) X=[0.99 1]; Y=[-0.6770-0.6834]; p8 = polyfit(X, Y, 1)

得到

p1 = -0.3965 0.3717

p2 = -0.4021 0.3721 p3 = -0.7600 0.4506 p4 = -0.3400 0.3006 p5 = -0.9326 0.5271 p6 = -1.3187 0.7300 p7 = -2.0731 1.3696 p8 = -0.00 -0.0434

输入：

S1=quad('-0.3965*x+0.3717',0,0.0819)+quad('-0.4021*x+0.3721',0.0819,0.2192)+quad('-0.7600*x+0.4506',0.2192,0.3584)+quad('-0.3400*x+0.3006',0.3584,0.3831)+quad('-0.9326*x+0.5271',0.3831,0.5256)+quad('-1.3187*x+0.7300',0.5256,0.8478)+quad('-2.0731*x+1.3696',0.8478,0.9872)+quad('-0.00*x-0.434',0.9872,1)

输出：S1 = -0.0214 输入：

S=quad('abs(-0.3965*x+0.3717)',0,0.0819)+quad('abs(-0.4021*x+0.3721)',0.0819,0.2192)+quad('abs(-0.7600*x+0.4506)',0.2192,0.3584)+quad('abs(-0.3400*x+0.3006)',0.3584,0.3831)+quad('abs(-0.9326*x+0.5271)',0.3831,0.5256)+quad('abs(-1.3187*x+0.7300)',0.5256,0.8478)+quad('abs(-2.0731*x+1.3696)',0.8478,0.9872)+quad('abs(-0.00*x-0.434)',0.9872,1)

输出：S = 0.2686

0.0214÷0.2686=0.08＞0.02，热力学一致性检验结果与原方法相比较大。