如何拟合安托因方程

abingchem

安托因方程是化学工程里一个常用方程，形式如下：
lgP=a-b/(T+c)
有组数据如下（从网上摘抄）
T=[ 293.6500  295.6500  297.6500  299.6500  301.6500  303.9000]
P=[ 750        5940        6350        6600        6940        7300]
用cftool拟合，总是不理想，是默认的初值问题，还是其他问题？求高手指教，最好给出过程，多谢！

附上网上拟合好的结果：
General model:
f(x) = a-b./(c+x)
Coefficients (with 95% confidence bounds):
a = 12.95 (-337.1, 363)
b = 8065 (-6.098e+005, 6.259e+005)
c = 583.9 (-3.322e+004, 3.439e+004)

Goodness of fit:
SSE: 3.999e-005
R-square: 0.9949
Adjusted R-square: 0.9915
RMSE: 0.003651
______________

abingchem

补充另外一组数据，供试验：

温度/C        压力/mmHg
-67.4        0.102
-60.4        0.212
-54.5        0.378
-48.1        0.702
-44.4        0.982
-44        1
-40        2
-30        4
-20        8
-10        15.5
0        29.6
10        54.7
20        96
30        160
40        260.5
50        406
60        625
64.7        760
70        927
80        1341
90        1897
100        2621
110        3561
120        4751
130        6242
140        8071
150        10336
160        13027
170        16292
180        20089
190        24615
200        29787
210        35770
220        42573
230        50414
240        59660

shamohu

1楼问题参考结果：

参数        最佳估算
----------        -------------
a        8.91267919808239
b        0.549416729222603
c        -293.409633846151

abingchem

shamohu 发表于 2012-1-18 09:09
1楼问题参考结果：

参数        最佳估算

版主，这是怎么做出来的呢？需要过程啊。。。。多谢

stats01

clear,clc
%question 1
T=[293.65  295.65  297.65  299.65  301.65  303.9]';
P=[750 5940 6350  6600  6940 7300]';
P=log10(P);
fx1=@(b,x)b(1)+b(2)./(x+b(3));
b1=[9  0.5  -293];
for l=1:5
    b1=lsqcurvefit(fx1,b1,T,P);
    b1=nlinfit(T,P,fx1,b1);
end
b1
figure(1),clf
plot(T,P,'o','markerfacecolor','k')
hold on
x11=linspace(min(T),max(T),200);
y11=fx1(b1,x11);
plot(x11,y11,'r-','linewidth',2)
legend('data','fit','location','best')
axis tight

%question 2
x2=[-67.4, -60.4, -54.5, -48.1, -44.4, -44.0, -40.0, -30.0, -20.0, -10.0, 0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 64.7, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 170.0, 180.0, 190.0, 200.0, 210.0, 220.0, 230.0, 240.0]';
y2=[0.102, 0.212, 0.378, 0.702, 0.982, 1.0, 2.0, 4.0, 8.0, 15.5, 29.6, 54.7, 96.0, 160.0, 260.5, 406.0, 625.0, 760.0, 927.0, 1341.0, 1897.0, 2621.0, 3561.0, 4751.0, 6242.0, 8071.0, 1.034*10^4, 1.303*10^4, 1.629*10^4, 2.009*10^4, 2.461*10^4, 2.979*10^4, 3.577*10^4, 4.257*10^4, 5.041*10^4, 5.966*10^4]';
fx2=@(b,x)b(1)+b(2)*exp(b(3)*x+b(4)*x.^2);
b2=[-105  110  0.04  -.0001];
for l=1:5
    b2=lsqcurvefit(fx2,b2,x2,y2);
    b2=nlinfit(x2,y2,fx2,b2);
end
b2
figure(2),clf
plot(x2,y2,'o','markerfacecolor','k')
hold on
x21=linspace(min(x2),max(x2),200);
y21=fx2(b2,x21);
plot(x21,y21,'r-','linewidth',2)
legend('data','fit','location','best')
axis tight