2Â¥: Originally posted by dingd at 2014-03-12 09:56:13
k1, k2, k3ÒªÇ󶼴óÓÚ0µÄ»°Ð§¹û»áºÜ²îµÄ¡£

3Â¥: Originally posted by 805592836_01 at 2014-03-12 10:02:43
àÅ£¬¶¼Òª´óÓÚ0µÄ¡­¡­...

4Â¥: Originally posted by dingd at 2014-03-12 10:35:53
²»¿ÉÄÜÓкýá¹û£¬¼ì²éÏÂÄ£ÐÍ°É¡£...

function test_123
clear all;clc

format long

tspan=[0,180,300,420,600,900,1200,1800,2400];

yexp=[1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124]';




y0=1.845;

k0=[-1 1 1];
lb=-[inf 0 0];
ub=+[inf inf inf];        

yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\nʹÓú¯Êýlsqnonlin()¹À¼ÆµÃµ½µÄ²ÎÊýֵΪ:\n')
fprintf('\tk1=%.4f\n',k(1)),  %±¾Õ÷¶¯Á¦Ñ§²ÎÊý
fprintf('\tk2=%.4f\n',k(2)),  % AµÄÎü¸½Æ½ºâ³£ÊýKA
fprintf('\tk3=%.4f\n',k(3)),  % H2µÄÎü¸½Æ½ºâ³£ÊýKB
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\tÏà¹ØϵÊý֮ƽ·½R^2 = %.6f',R2);
plot(ts,ys,'b',tspan,yy,'or'),legend('¼ÆËãÖµ','ʵÑéÖµ','Location','best'),


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % Ä¿±êº¯Êý
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
ysim = Xsim(2:end);
f=ysim-yexp;
%----------------------------------------------------------

function dxdt = KineticsEqs(t,y,k)
y0=1.845;
CB=y;
CB0=0.4;
CA0=y0;
k1=k(1);
KA=k(2);
KB=k(3);
dxdt=k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA0+k(3)^(1/6)*CB.^(1/2)).^7;