CODE:
function feixianxingnihe_0508
clear all;clc
format long
data=[2.42e-6 353.15
2.64e-6 363.15
2.71e-6 373.15
];
xspan=data(:,1);
Texp=data(:,2);
k0=[1e-1 -1e-5 1e-5];
lb=-[0 1 0]*1;
ub=[1 0 1]*1;
%-------------------------------------------------------------------------
% 使用函数lsqnonlin()进行参数估计
OPTIONS=optimset('MaxFunEvals',1000);
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,xspan,Texp);
ci = nlparci(k,residual,jacobian);
%residual;
fprintf('\n\n拟合结果:\n')
fprintf('\n\t参数 a = %.16f',k(1))
fprintf('\n\t参数 b = %.16f',k(2))
fprintf('\n\t参数 c = %.16f',k(3))
y=KineticsEqs(xspan,k);
R2=1-sum((Texp-y).^2)./sum((Texp-mean(y)).^2);
fprintf('\n\t相关系数之平方R^2 = %.16f',R2);
figure(1)
mm=max(xspan)-min(xspan);
xspan1=min(xspan)-0.1*mm:1e-8:max(xspan)+0.1*mm;
plot(xspan1,KineticsEqs(xspan1,k),'b-',xspan,Texp,'or'),legend('计算值','实验值','Location','Best')
%-------------------------------------------------------------------------
function f = ObjFunc(k,xspan,Texp)
f=KineticsEqs(xspan,k)-Texp;
%------------------------------------------------------------------------
function T = KineticsEqs(x,k)
R=8.314;p=101325;
V=x;
a=k(1);b=k(2);c=k(3);
T=(p+a./(V.*(V+b)+c*(V-b))).*(V-b)/R;