24小时热门版块排行榜    

查看: 1588  |  回复: 3
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

chaiqing123

金虫 (小有名气)

[求助] MATLAB用于三元回归参数

function文件为下
function rcal = fun552(x,x1,x2)
for n=1:11%m为所取的数据点个数
    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(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));

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+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));
rcal=[r1cal r2cal];
end
script文件如下
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%
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)));

%A=7.23255; B=1750.286;C=235;
%A1=7.16879;B1=1552.601;C1=222.419;
%p1s=10.^(A-B./(t+C))
%p2s=10.^(A1-B1./(t+C1))
p3s=0;
%a为aotonie方程的三个系数
P=101.3%总压
r1exp=y1.*P./(x1.*p1s);
r2exp=y2.*P./(x2.*p2s);
rexp=[r1exp r2exp]
x0=[1,1,1,1,1,1];
xx=[xx1 xx1];
x=nlinfit(xx,rexp,'fun552',x0)

最小二乘法 这个有问题 但不知道怎么改....
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chaiqing123

金虫 (小有名气)

引用回帖:
2楼: Originally posted by 695710829 at 2012-04-01 10:52:29:
x=nlinfit(xx,rexp,@fun552,x0)

chaiqing123@126.com  我的邮箱 如果你有源程序 麻烦 发一份 我参考下
4楼2012-04-01 14:48:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 4 个回答

695710829

新虫 (初入文坛)

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
chaiqing123: 金币+1 2012-04-06 19:46:41
chaiqing123: 金币+1, 有帮助 2012-04-12 09:50:54
x=nlinfit(xx,rexp,@fun552,x0)
2楼2012-04-01 10:52:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chaiqing123

金虫 (小有名气)

引用回帖:
2楼: Originally posted by 695710829 at 2012-04-01 10:52:29:
x=nlinfit(xx,rexp,@fun552,x0)

不是@的问题 我两个程序不是放在一个文件.m文件中的,我现在又把程序修改了下 放在同一文件中 不过 只能成功运行 对r1的拟合,如何r1与r2 放在一块 拟合的话 老是 说矩阵 长度不一致
3楼2012-04-01 14:46:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
信息提示
请填处理意见