CODE:

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