matlab拟合Sanchez-lacombe方程
使用Sanchez-lacombe方程研究CO2-polymer体系,要先计算纯物质的特征参数,写了个方程用matlab进行拟合。
function y=myfun(beta,x) %输入温度压力密度矩阵,输入同长度0向量,回归特征参数
y = (x(:,3)./beta(3)).^2+x(:,2)./beta(2)+x(:,1)./beta(1).*(log(1-x(:,3)./beta(3))+(1-8.314*beta(1)*beta(3)./44./beta(2)).*x(:,3)./beta(3));
x=[313 8 0.279493558
313 10 0.631751848
313 12 0.719372707
313 14 0.764467548
313 16 0.795924865
313 18 0.820411847
313 20 0.840618695
313 22 0.857927248
313 24 0.873133677
313 26 0.890709896
313 28 0.899118864
333 8 0.191945948
333 10 0.290807573
333 12 0.436262106
333 14 0.563094769
333 16 0.638855172
333 18 0.688325991
333 20 0.724637681
333 22 0.753238927
333 24 0.7768197
333 26 0.797003268
333 28 0.814663951
353 8 0.160513644
353 10 0.221931246
353 12 0.29731819
353 14 0.384260682
353 16 0.469483568
353 18 0.540102619
353 20 0.594848611
353 22 0.637592451
353 24 0.672043011
353 26 0.700672646
353 28 0.725058005
373 8 0.141392718
373 10 0.188757597
373 12 0.242371361
373 14 0.301731941
373 16 0.364258915
373 18 0.425477599
373 20 0.481255113
373 22 0.529661017
373 24 0.5709066
373 26 0.606023877
373 28 0.636213259
393 8 0.127831467
393 10 0.167442484
393 12 0.210499726
393 14 0.256673511
393 16 0.305026842
393 18 0.353969771
393 20 0.401638686
393 22 0.446388715
393 24 0.487258198
393 26 0.523944252
393 28 0.556637907
413 8 0.117430158
413 10 0.151987233
413 12 0.188665006
413 14 0.227205598
413 16 0.267115421
413 18 0.307654443
413 20 0.347910796
413 22 0.386981928
413 24 0.424124184
413 26 0.458862938
413 28 0.490990327
433 8 0.109060769
433 10 0.14000112
433 12 0.172327629
433 14 0.20584179
433 16 0.240228698
433 18 0.275065328
433 20 0.309837335
433 22 0.344044588
433 24 0.377187689
433 26 0.408914332
433 28 0.438962293
453 8 0.102104371
453 10 0.13029316
453 12 0.159420983
453 14 0.18933298
453 16 0.219809206
453 18 0.250576326
453 20 0.281325606
453 22 0.311720698
453 24 0.341460083
453 26 0.370274373
453 28 0.397962432
473 8 0.096181591
473 10 0.122186652
473 12 0.148840532
473 14 0.176022249
473 16 0.203566485
473 18 0.231288741
473 20 0.258966723
473 22 0.286393447
473 24 0.313361745
473 26 0.339662376
473 28 0.365176746
493 8 0.091049804
493 10 0.115267132
493 12 0.139934511
493 14 0.164953895
493 16 0.19020447
493 18 0.21554976
493 20 0.240830383
493 22 0.265900872
493 24 0.290604748
493 26 0.3148119
493 28 0.338409475];
y=[0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0];
beta0=[600,400,1.5];
[beta]=lsqcurvefit(@myfun,beta0,x,y]
拟合结果和使用1stopt差距很大。不知道是我方法不对还是什么原因 返回小木虫查看更多
Matlab里的“log()”函数在1stOpt里是“ln()”,不知楼主在用1stOpt的时候代码对否?
这一点我是注意到了,在1stopt中也改成ln了,拟合结果相差几个数量级
,
个人觉得你的公式有问题,左边y值全等于0,相当于求解超越方程:
结果1:
beta3: 1.93947598320175E82
beta2: 1.55545140991017E236
beta1: 3.05530394179693E82
结果2:
beta3: -2.90986946492629E161
beta2: 2.60007550366803E163
beta1: 2.87013631764659E139
结果非常大,且不稳定。