回 帖 发 新 帖 刷新版面

主题:求救 NRTL 三元回归参数的问题....

希望高手能帮下忙, 我做了很久了 现在很头疼啊...
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 ;
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];
x=nlinfit(x1,r1exp,@fun552,beta0)%最小二乘法拟合


function r1cal = fun552(x,x1,x2,x3)
 R=8.315;
 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;
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;
%x1=[x1 x1];
%x2=[x2 x2];
%x3=[x3 x3];
T=[364 360.2 358.5 356.3 355.3 354.5 354.3 354.1 353.9 353.7 353.6];

%T=[T T];


for n=1:11%11个数据点
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(n)=(tao(3)*G(3)*x2(n)+tao(5)*G(5)*x3(n))/(x1(n)+G(3)*x2(n)+G(5)*x3(n));
B(n)=-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(n)=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(n)=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(n)+B(n)+C(n)+D(n));%组分1的活度系数计算值
end
现在 我只将r1进行与r1实验值拟合 结果 可行 但是 我想把r1和r2 的实验值一块儿 与计算值 都进行 拟合回归得到参数 老是报错。
程序如下

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 ;
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];
beta=nlinfit(xx,rexp,@fun552,beta0);%最小二乘法拟合


function rcal = fun552(x,x1,x2,x3)
 R=8.315;
 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;

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];

for n=1:11
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(n)=(tao(3)*G(3)*x2(n)+tao(5)*G(5)*x3(n))/(x1(n)+G(3)*x2(n)+G(5)*x3(n));
B(n)=-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(n)=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(n)=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(n)+B(n)+C(n)+D(n));%组分1的活度系数计算值

A1(n)=(tao(1)*G(1)*x1(n)+tao(6)*G(6)*x3(n))/(G(1)*x1(n)+G(6)*x3(n)+x2(n));
B1(n)=-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(n)=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(n)=x3(n)*G(4)/(G(2)*x1(n)+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(n)+B1(n)+C1(n)+D1(n)); %组分2的活度系数计算值

rcal=[r1cal r2cal];
end

回复列表 (共1个回复)

沙发


或者谁有现成的程序 发给我 参考下 谢谢
chaiqing123@126.com
邮箱 感激不尽!

我来回复

您尚未登录,请登录后再回复。点此登录或注册