3楼: Originally posted by 月只蓝 at 2014-07-28 20:03:15
需要指出的是，截图中的计算步骤，在输入y1=... 时，”1./1.5003“这一项的点除是多余的。

function feixianxingnihe_3
clear all;clc
format long


xdata=[0 50.65 293.77 790.14 992.74 5065 10130 15195];
ydata=[0.5650 0.4013 0.2502 0.2324 0.2307 0.1926 0.1812 0.1730];
xspan=xdata;    %x的数据，在此输入
Texp=ydata;    %T的数据，在此输入

x0=[0.1 0.1 0.1 1];
k0=x0;
lb=-[1 1 1 1]*1e9;
ub=[1 1 1 1]*1e9;


%-------------------------------------------------------------------------

% 使用函数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参数 c1 = %.16f',k(1))
fprintf('\n\t参数 c2 = %.16f',k(2))
fprintf('\n\t参数 c3 = %.16f',k(3))
fprintf('\n\t参数 c4 = %.16f',k(4))
fprintf('\n\t残差平方和 = %.6e',resnorm)
y=KineticsEqs(xspan,k);
R2=1-sum((Texp-y).^2)./sum((Texp-mean(y)).^2);
%fprintf('\n\t相关系数之平方R^2 = %.6e',R2);

mm=max(xspan)-min(xspan);
xspan1=min(xspan):0.01:max(xspan)+0.05*mm;
figure(3)
plot(xspan1,KineticsEqs(xspan1,k),'b',xspan,Texp,'or'),legend('拟合曲线','实验数据','Location','Best'),...
    axis([-10 16000 0.1 0.6])


%-------------------------------------------------------------------------

function f = ObjFunc(k,xspan,Texp)
f=KineticsEqs(xspan,k)-Texp;

%------------------------------------------------------------------------
function xt = KineticsEqs(x,n)
k1=n(1);
k2=n(2);
k3=n(3);
k4=n(4);
xt=k1+(k2-k1)./(1+(k3*x).^k4).^(1-1/k4);