|
|
???????,???!????????}??!??????????????????W?,?F????????3??,??6??k1 k2 k3 k4 k5 k6!
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC
???????????????a!???????RUN,??????????!
function ODE_parafit2
%-------data
yexp=[5.14 12.76 45.41 21.52 3.31;
5.36 13.20 46.94 20.90 3.17;
5.51 13.47 48.68 20.01 3.08;
5.78 14.03 48.16 19.34 3.04;
6.07 14.33 47.80 18.65 2.94]';
x0=[93.7 6.3 0 0 0];
k0 = [1 1 1 1 1 1];
lb = [0 0 0 0 0 0];
ub = [1 1 1 1 1 1]*1e5;
tspan = [0 38.46 50 62.5 100 166.666 250 500];
%-----------------
[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);
fprintf('\n\n??ú???fmincon()?????????????:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',fval)
[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 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',resnorm)
%% plot figure
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i, =x(end,;
end
Ypre=fy(;
Yexp=yexp(;
plot(Ypre,'or',Yexp,'b.')
legend('prediction','experiment')
%----------fmincon fun-------------------
function f = ObjFuncFmincon(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i)=sum((yexp(i,-x(end,).^2);
end
f = sum(fy) ;
% -------------------lsqnonlin fun-------------------
function f = ObjFunc(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);
fy(i,=yexp(i,-x(end,;
end
f =fy(;
%-----------odefun--------------------------------
function dYdt = KineticsEqs(t,y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);
DEC=y(1);PA=y(2);DPC=y(3);EA=y(4);EPC=y(5);
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA;
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA;
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC;
dy1dt=-R1+R3;
dy2dt=-R1-R2;
dy3dt=R2+R3;
dy4dt=R1+R2;
dy5dt=R1-R2-R3;
dYdt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt]; |
|
回复此楼
