function k1k2k3
format long
clear all
clc
tspan = 0:12:144;
x0 = [0.35  0      99.8];

k0 = [0.5 8 1.5 2 -0.5 0.1 0.5 5 -0.01];   %k的初值，最需要调节

lb = [1 1 1 1 1 1 1 1 1 ]*-inf;
ub = [1 1 1 1 1 1 1 1 1]*inf;

data=...
    [


12  2.75  0.55   88.96
24  4.6   3.72   73.49
36  6.27  8.19   61.57
48  8.07  14.12  47.38
60  9.73 19.37  33.63
72  10.41 22.54  24.65
84  10.7  27.61  11.19
96  10.53 32.49  3.2
108 10.56 35.62  0
120 10.59 37.58  0
132 10.4  38.44  0
144 10.77 39.88  0
];
yexp = data(:,2:4);
options=optimset('MaxFunEvals',1500);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);

[t uu]=ode45(@KineticsEqs,tspan,x0,[],k);
uexp=[0   0.35  0      99.8
12  2.75  0.55   88.96
24  4.6   3.72   73.49
36  6.27  8.19   61.57
48  8.07  14.12  47.38
60  9.73 19.37  33.63
72  10.41 22.54  24.65
84  10.7  27.61  11.19
96  10.53 32.49  3.2
108 10.56 35.62  0
120 10.59 37.58  0
132 10.4  38.44  0
144 10.77 39.88  0
];
figure(1)

size(uexp(:,2));
size(uexp(:,3));
size(uexp(:,4));
size(uu(:,1));
size(uu(:,2));
size(uu(:,3));

R1_Square=1-sum((uexp(2:end,2)-uu(2:end,1)).^2)./sum((uexp(2:end,2)-mean( uu(2:end,1)  )).^2);
R2_Square=1-sum((uexp(2:end,3)-uu(2:end,2)).^2)./sum((uexp(2:end,3)-mean( uu(2:end,2)  )).^2);
R3_Square=1-sum((uexp(2:end,4)-uu(2:end,3)).^2)./sum((uexp(2:end,4)-mean( uu(2:end,3)  )).^2);



plot(tspan,uexp(:,2),'or',tspan,uu(:,1),'-r',tspan,uexp(:,3),'>b',tspan,uu(:,2),'b-',tspan,uexp(:,4),'*k',tspan,uu(:,3),'k-')


fprintf('\n\tk1=%.9f',k(1))
fprintf('\n\tk2=%.9f',k(2))
fprintf('\n\tk3=%.9f',k(3))
fprintf('\n\tk4=%.9f',k(4))
fprintf('\n\tk5=%.9f',k(5))
fprintf('\n\tk6=%.9f',k(6))
fprintf('\n\tk7=%.9f',k(7))
fprintf('\n\tk8=%.9f',k(8))
fprintf('\n\tk9=%.9f',k(9))

fprintf('\n\t方程一的决定系数R1_Square=%.6f',R1_Square);
fprintf('\n\t方程二的决定系数R2_Square=%.6f',R2_Square);
fprintf('\n\t方程三的决定系数R3_Square=%.6f',R3_Square);




function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,2);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);


f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3))];

function dCdt = KineticsEqs(t,u,k)  
x=u(1);
p=u(2);
s=u(3);

a=k(1);
b=k(2);
c=k(3);
d=k(4);
e=k(5);
f=k(6);
g=k(7);
h=k(8);
i=k(9);



dxdt = a*x*(1-x/b).*(1+s/c).^(-1);
dpdt=d*dxdt+e*(s./(s+f)).*x;
dsdt=g*dxdt+h*dpdt+i*x;
dCdt = [dxdt;dpdt;dsdt];

2楼: Originally posted by 月只蓝 at 2015-08-31 16:48:55
要拟合九个参数，这种问题还是推荐用1stOpt，参数个数太多，用MATLAB比较困难。
下面给出MATLAB参考代码，代码中初值k0可以试着调一下，我初步调了一下初值，拟合效果不好。

function k1k2k3
format long
cle ...

2楼: Originally posted by 月只蓝 at 2015-08-31 16:48:55
要拟合九个参数，这种问题还是推荐用1stOpt，参数个数太多，用MATLAB比较困难。
下面给出MATLAB参考代码，代码中初值k0可以试着调一下，我初步调了一下初值，拟合效果不好。

function k1k2k3
format long
cle ...

4楼: Originally posted by 皓小天 at 2015-08-31 22:03:03
我记得已经回复了你这个问题，http://muchong.com/bbs/viewthread.php?tid=9255377&authorid=1864745