您所在位置：主题：[讨论]跪求，求助：NRTL 用matlab拟合回归三元参数 方程
chaiqing123
[专家分：0] 发布于 2012-04-01 10:15:00
具体代码如下，编了两三个星期了，没编出来。
function NRTL
xx1=[0.05 0.105 0.153 0.248 0.347 0.446 0.474 0.535 0.624 0.740 0.870];%不含离子液体的x1
xx2=1-xx1;
x1=xx1/1.05%离子液体浓度为5%，含离子液体的x1
x2=xx2/1.05
x3=1-x1-x2
global x1 x2 x3 T;
y1=[0.397 0.483 0.521 0.601 0.635 0.676 0.685 0.713 0.745 0.825 0.897];
y2=1-y1;
T=[364 360.2 358.5 356.3 355.3 354.5 354.3 354.1 353.9 353.7 353.6];
t=T-273.15;
a=[7.23255 1750.286 235;7.16879 1552.601 222.419];
p1s=10.^(a(1,1)-a(1,2)./(t+a(1,3)))
p2s=10.^(a(2,1)-a(2,2)./(t+a(2,3)))
p3s=0;
P=101.3%总压
r1exp=y1.*P./(x1.*p1s);%组分1的活度系数实验值
r2exp=y2.*P./(x2.*p2s);
rexp=[r1exp r2exp]%组分1,2的活度系数计算值
beta0=[1,1,1,1,1,1];
xx=[x1 x1];
A=nlinfit(xx,rexp,@fun552,beta0)%最小二乘法拟合
function rcal = fun552(x,x1,x2,x3)
for n=1:11%11个数据点
R=8.315;
global x1 x2 x3 T;
tao12=x(1)./R/T(n);%x(1)=g12-g22
tao21=x(2)./R/T(n); %x(2)=g21-g11
tao13=x(3)./R/T(n); %x(3)=g13-g33
tao31=x(4)./R/T(n); %x(4)=g31-g11
tao23=x(5)./R/T(n); %x(5)=g23-g33
tao32=x(6)./R/T(n); %x(6)=g32-g22
tao=[tao12,tao13,tao21,tao23,tao31,tao32];
G=exp(-0.3*tao);
A=(tao(3)*G(3)*x2(n)+tao(5)*G(5)*x3(n))/(x1(n)+G(3)*x2(n)+G(5)*x3(n));
B=-x1(n)*(x2(n)*tao(3)*G(3)+x3(n)*tao(5)*G(5))/((x1(n)+G(3)*x2(n)+G(5)*x3(n))^2);
C=x2(n)*G(1)*(tao(1)-(x1(n)*tao(1)*G(1)+x3(n)*tao(6)*G(6))/(G(1)*x1(n)+x2(n)+G(6)*x3(n)))/(G(1)*x1(n)+x2(n)+G(6)*x3(n));
D=x3(n)*G(2)*(tao(2)-(x1(n)*tao(2)*G(2)+x2(n)*tao(4)*G(4))/(G(2)*x1(n)+G(4)*x2(n)+x3(n)))/(G(2)*x3(n)+G(4)*x2(n)+x3(n));
r1cal(n)=exp(A+B+C+D);%组分1的活度系数计算值
%A1=(tao12*G12*x1+tao32*G32*x3)/(G12*x1+G32*x3+x2);
%B1=-x2*(x1*tao12*G12+x3*tao32*G32)/((G12*x1+x2+G32*x3)^2);
%C1=x1*G21/(x1+G21x2+G31x3)*(tao21-(x2*tao21*G21+x3*tao31*G31)/(x1+G21*x2+G31*x3));
%D1=x3*G23/(G13*x1+G23*x2+x3)*(tao23-(x1*tao13*G13+x2*tao23*G23)/(G13*x1+G23*x2+x3));
A1=(tao(1)*G(1)*x1(n)+tao(6)*G(6)*x3(n))/(G(1)*x1(n)+G(6)*x3(n)+x2(n));
B1=-x2(n)*(x1(n)*tao(1)*G(1)+x3(n)*tao(6)*G(6))/((G(1)*x1(n)+x2(n)+G(6)*x3(n))^2);
C1=x1(n)*G(3)/(x1(n)+G(3)*x2(n)+G(5)*x3(n))*(tao(3)-(x2(n)*tao(3)*G(3)+x3(n)*tao(5)*G(5))/(x1(n)+G(3)*x2(n)+G(5)*x3(n)));
D1=x3(n)*G(4)/(G(2)*x1+G(4)*x2(n)+x3(n))*(tao(4)-(x1(n)*tao(2)*G(2)+x2(n)*tao(4)*G(4))/(G(2)*x1(n)+G(4)*x2(n)+x3(n)));
r2cal(n)=exp(A1+B1+C1+D1); %组分2的活度系数计算值
rcal=[r1cal r2cal];%组分1,2的活度系数计算值
end
老是提示nlinfit报错。[em10][em10]
function NRTL
xx1=[0.05 0.105 0.153 0.248 0.347 0.446 0.474 0.535 0.624 0.740 0.870];%不含离子液体的x1
xx2=1-xx1;
x1=xx1/1.05%离子液体浓度为5%，含离子液体的x1
x2=xx2/1.05
x3=1-x1-x2
global x1 x2 x3 T;
y1=[0.397 0.483 0.521 0.601 0.635 0.676 0.685 0.713 0.745 0.825 0.897];
y2=1-y1;
T=[364 360.2 358.5 356.3 355.3 354.5 354.3 354.1 353.9 353.7 353.6];
t=T-273.15;
a=[7.23255 1750.286 235;7.16879 1552.601 222.419];
p1s=10.^(a(1,1)-a(1,2)./(t+a(1,3)))
p2s=10.^(a(2,1)-a(2,2)./(t+a(2,3)))
p3s=0;
P=101.3%总压
r1exp=y1.*P./(x1.*p1s);%组分1的活度系数实验值
r2exp=y2.*P./(x2.*p2s);
rexp=[r1exp r2exp]%组分1,2的活度系数计算值
beta0=[1,1,1,1,1,1];
xx=[x1 x1];
A=nlinfit(xx,rexp,@fun552,beta0)%最小二乘法拟合
function rcal = fun552(x,x1,x2,x3)
for n=1:11%11个数据点
R=8.315;
global x1 x2 x3 T;
tao12=x(1)./R/T(n);%x(1)=g12-g22
tao21=x(2)./R/T(n); %x(2)=g21-g11
tao13=x(3)./R/T(n); %x(3)=g13-g33
tao31=x(4)./R/T(n); %x(4)=g31-g11
tao23=x(5)./R/T(n); %x(5)=g23-g33
tao32=x(6)./R/T(n); %x(6)=g32-g22
tao=[tao12,tao13,tao21,tao23,tao31,tao32];
G=exp(-0.3*tao);
A=(tao(3)*G(3)*x2(n)+tao(5)*G(5)*x3(n))/(x1(n)+G(3)*x2(n)+G(5)*x3(n));
B=-x1(n)*(x2(n)*tao(3)*G(3)+x3(n)*tao(5)*G(5))/((x1(n)+G(3)*x2(n)+G(5)*x3(n))^2);
C=x2(n)*G(1)*(tao(1)-(x1(n)*tao(1)*G(1)+x3(n)*tao(6)*G(6))/(G(1)*x1(n)+x2(n)+G(6)*x3(n)))/(G(1)*x1(n)+x2(n)+G(6)*x3(n));
D=x3(n)*G(2)*(tao(2)-(x1(n)*tao(2)*G(2)+x2(n)*tao(4)*G(4))/(G(2)*x1(n)+G(4)*x2(n)+x3(n)))/(G(2)*x3(n)+G(4)*x2(n)+x3(n));
r1cal(n)=exp(A+B+C+D);%组分1的活度系数计算值
%A1=(tao12*G12*x1+tao32*G32*x3)/(G12*x1+G32*x3+x2);
%B1=-x2*(x1*tao12*G12+x3*tao32*G32)/((G12*x1+x2+G32*x3)^2);
%C1=x1*G21/(x1+G21x2+G31x3)*(tao21-(x2*tao21*G21+x3*tao31*G31)/(x1+G21*x2+G31*x3));
%D1=x3*G23/(G13*x1+G23*x2+x3)*(tao23-(x1*tao13*G13+x2*tao23*G23)/(G13*x1+G23*x2+x3));
A1=(tao(1)*G(1)*x1(n)+tao(6)*G(6)*x3(n))/(G(1)*x1(n)+G(6)*x3(n)+x2(n));
B1=-x2(n)*(x1(n)*tao(1)*G(1)+x3(n)*tao(6)*G(6))/((G(1)*x1(n)+x2(n)+G(6)*x3(n))^2);
C1=x1(n)*G(3)/(x1(n)+G(3)*x2(n)+G(5)*x3(n))*(tao(3)-(x2(n)*tao(3)*G(3)+x3(n)*tao(5)*G(5))/(x1(n)+G(3)*x2(n)+G(5)*x3(n)));
D1=x3(n)*G(4)/(G(2)*x1+G(4)*x2(n)+x3(n))*(tao(4)-(x1(n)*tao(2)*G(2)+x2(n)*tao(4)*G(4))/(G(2)*x1(n)+G(4)*x2(n)+x3(n)));
r2cal(n)=exp(A1+B1+C1+D1); %组分2的活度系数计算值
rcal=[r1cal r2cal];%组分1,2的活度系数计算值
end
老是提示nlinfit报错。[em10][em10]