Çó¶¯Á¦Ñ§²ÎÊý
·½³ÌΪ£ºr1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
t c
0 8.769713879
3 5.909461709
6 5.736047766
9 4.756936544
12 4.505799773
15 4.354160734
18 4.124025932
21 3.994725485
31 3.732949397
41 3.680613046
51 3.537445782
61 3.283577643
Çók1,k2,k3,,k4
ÏÂÃæÊÇдµÄ³ÌÐò£¬´ó¼Ò°ïæ¿´¿´£¬ÎªÊ²Ã´½á¹ûÊÇÖ±ÏßÄØ£¿
function parafit
%
% r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
%
% dCAdt = - r1;
clear all
clc
% t/min CA / mol/L
Kinetics=[0 7.029559717
1 3.960706476
2 3.163862085
3 2.963057756
6 2.820434159
9 2.698109429
12 2.556212898
15 2.447100967
18 2.207246045
21 2.13921916
31 2.000403189
41 1.917172016
51 1.851299168
61 1.636287388];
c=Kinetics( :,1)';
k0 = [0.671 0.02253 0.641 0.04028]; % ²ÎÊý³õÖµ
lb = [0 0 0 0]; % ²ÎÊýÏÂÏÞ
ub = [100 100 100 100]; % ²ÎÊýÉÏÏÞ
x0 = [7.029559717];
yexp = Kinetics;
warning off
% ʹÓú¯Êýfmincon()½øÐвÎÊý¹À¼Æ
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\nʹÓú¯Êýfmincon()¹À¼ÆµÃµ½µÄ²ÎÊýֵΪ:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf(' The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
warning off
% ʹÓú¯Êýlsqnonlin()½øÐвÎÊý¹À¼Æ
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\nʹÓú¯Êýlsqnonlin()¹À¼ÆµÃµ½µÄ²ÎÊýֵΪ:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf(' The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% ÒÔº¯Êýfmincon()¹À¼ÆµÃµ½µÄ½á¹ûΪ³õÖµ£¬Ê¹Óú¯Êýlsqnonlin()½øÐвÎÊý¹À¼Æ
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\nÒÔfmincon()µÄ½á¹ûΪ³õÖµ£¬Ê¹Óú¯Êýlsqnonlin()¹À¼ÆµÃµ½µÄ²ÎÊýֵΪ:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf(' The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('ca-pr','ca-real')
x(:,1)
figure;plot(t,x);
hold on
% ------------------------------------------------------------------
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,2) = x(:,1);
f = sum((y(:,2)-yexp(:,2)).^2) ;
% ------------------------------------------------------------------
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,2) = x(:,1);
f1 = y(:,2) - yexp(:,2);
f = [f1];
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
Kinetics=[0 7.029559717
1 3.960706476
2 3.163862085
3 2.963057756
6 2.820434159
9 2.698109429
12 2.556212898
15 2.447100967
18 2.207246045
21 2.13921916
31 2.000403189
41 1.917172016
51 1.851299168
61 1.636287388];
c=Kinetics(:,1)';
for i=1:14
c(i)=c(:,i);
r1 =(k(1)*k(2)*c(i)^2-0.25*k(3)*k(4)*(9.404-c(i))^2)/(1+k(2)*c(i)+0.5*k(4)*(9.404-c(i)));
dCAdt = - r1;
dxdt = [dCAdt];
end
½á¹ûÈçÏ£º
![Çó´ó¼Ò°ïÎÒ¿´¿´ÏÂÃæÖÐmatlabÖеijÌÐòÎÊÌâ³öÔÚÄÄÀ]()
H8_DJFZ~6NZC}07J0EXL1]L.png |