function ode_fit33333333333
clear all;clc

format long

tspan=[0,2,5,7,10,15,20,30,40];
t1=length(tspan)
yexp=[1.530,1.428,1.424,1.326,1.203,1.017,0.910,0.775]';
t2=length(yexp)



y0=1.686;

k0=[5 -10 100];
lb=[-1 -1 -1]*1e9;
ub=[1 1 1]*1e9;        

yy=[y0 yexp'];
options = optimset('MaxFunEvals',1000);
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t反应速率常数 k1 = %.4f\n',k(1))
fprintf('\t反应速率常数 KA = %.4f\n',k(2))
fprintf('\t反应速率常数 KB = %.4f\n',k(3))
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) ;
Xsim;
ysim = Xsim(2:end);
f=ysim-yexp;
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
y0=1.686;
CB=y;
CB0=0.4;
CA0=y0;
k1=k(1);
KA=k(2);
KB=k(3);
dydt=(k1*KA*KB*CA0*CB0^3)./(1+KA*CA0+KB.^(1/3)*CB).^4;