function ode_fit
format long
clear all
clc

x0 = [6.7 4.86 3.17];
k0 = [1.5 1.5 1.5 6.8 2/3 0.1 0.1 0.01];   
lb = -[1 1 1 1 1 1 1 1 ]*1e9;
ub = [1 1 1 1 1 1 1 1 ]*1e9;

data=...
    [

88        6.7        3.17        4.86
92        7.5        2.87        5.17                                   
96        8.4        2.33        5.72                 
100        9.1        1.87        6.18
104        9.3        1.36        6.73               
108        9.1        0.95        7.13
112        8.6        0.52        7.56
];
tspan = [data(:,1)'-88];
yexp = [data(2:end,2) data(2:end,4) data(2:end,3)];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk6 = %.9f \n',k(6))
fprintf('\tk7 = %.9f \n',k(7))
fprintf('\tk8 = %.9f \n',k(8))

figure(1)

ts=0:((max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,4) data(:,3)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro',ts,ys(:,3),'k',tspan,yy(:,3),'ko'),
legend('X的计算值','X的实验值','P的计算值','P的实验值','S的计算值','S的实验值','Location','best');



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

size(ysim(:,1));
size(ysim(:,2));
size(yexp(:,1));
size(yexp(:,2));

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

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
umax=k(1);Ks=k(2);Ki=k(3);Pmax=k(4);a=k(5);Yxs=k(6);Yps=k(7);m=k(8);
X=C(1);P=C(2);S=C(3);

dXdt=umax*(S./(S+Ks)).*(1./(1+S/Ki)).*(1-P/Pmax).*X;
dPdt=a*dXdt;
dSdt=-Yxs*dXdt-Yps*dPdt-m*X;

dCdt = [dXdt; dPdt;dSdt];

ConstStr f1=umax*(S/(S+Ks))*(1/(1+S/Ki))*(1-P/Pmax)*X,
         f2=a*f1+b*X;
Variable t,X,S,P;
ODEFunction X'=f1;
            P'=f2;
            S'=-Yxs*f1-Yps*f2-m*X;
Data;
//t         X           S           P
88        6.7        3.17        4.86
92        7.5        2.87        5.17
96        8.4        2.33        5.72
100        9.1        1.87        6.18
104        9.3        1.36        6.73
108        9.1        0.95        7.13
112        8.6        0.52        7.56

6楼: Originally posted by 月只蓝 at 2014-03-22 15:39:25
k1~k8分别对应于1楼中除去β以外的8个参数。

5楼: Originally posted by 月只蓝 at 2014-03-22 15:36:01
1.  8个参数的初值，你自己再调一下；
2.  原问题没有给出t=0时的初值，代码编写时以t=88为零时刻。

代码和计算结果如下：

function ode_fit
format long
clear all
clc

x0 = ;
k0 = ;   
lb = -*1e ...

6楼: Originally posted by 月只蓝 at 2014-03-22 15:39:25
k1~k8分别对应于1楼中除去β以外的8个参数。