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);