function kinetics1
clear all
clc
a0=[10^11 81 0.3];
lb=[0 0 0];ub=[+inf +inf +inf];
c0=[0.0525564 100];
c16exp=...
[0.0525564 100
0.0403533 73.3083
0.0267938 50.3759
0.015037 42.1053
0.0082474 24.812
0.00415801 13.9098
0.00278195 11.6541
0.00185881 10.1504
0.00183435 9.77444
0.0015879 9.77444];
c14exp=...
[0.0525564 100
0.0493758 89.8496
0.0441658 81.203
0.0378303 72.9323
0.0357841 67.6692
0.0337337 61.6541
0.0276196 54.8872
0.0239832 47.7444
0.018766 40.6015
0.0160359 35.7143
0.0153455 34.5865];
c12exp=...
[0.0526 100
0.0516 98.1203
0.0507 94.7368
0.0484 93.2331
0.0464 89.8496
0.0439 80.4511
0.0429 80.8271
0.0407 77.8195
0.0393 76.3158
0.0381 75.5639
0.0365 72.1805];
[a,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@objlm,a0,lb,ub,[],c0,c16exp,c14exp,c12exp);
ci=nlparci(a,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n'),a
function f=objlm(a,c0,c16exp,c14exp,c12exp)
%160度
tspan1=[0 0.1005 0.2143 0.3360 0.4365 0.5661 0.6667 0.7593 0.8677 0.9603];
[t,c16]=ode45(@kinetic16,tspan1,c0,[],a);
f1=c16(:,1)-c16exp(:,1);
f2=c16(:,2)-c16exp(:,2);
%140度
tspan2=[0 0.1005 0.1984 0.2857 0.3571 0.4471 0.5529 0.6746 0.8042 0.9074 0.9683];
[t,c14]=ode45(@kinetic14,tspan2,c0,[],a);
f3=c14(:,1)-c14exp(:,1);
f4=c14(:,2)-c14exp(:,2);
%120度
tspan3=[0 0.1005 0.2011 0.2937 0.4021 0.5000 0.5979 0.7090 0.8069 0.8862 0.9894];
[t,c12]=ode45(@kinetic12,tspan3,c0,[],a);
f5=c12(:,1)-c12exp(:,1);
f6=c12(:,2)-c12exp(:,2);
f=[f1;f2;f3;f4;f5;f6]; %构造目标函数
function dc16dt=kinetic16(t,c16,a) %160度
T1=433.15;R=0.008314;
dc16dt=...
[-(a(1).*exp(-a(2)./(R.*T1)).*((c16(2)./100).^(a(3))).*c16(1).*0.943)
-(7.*a(1).*exp(-a(2)./(R.*T1)).*((c16(2)./100).^(a(3))).*c16(1).*0.943)];
function dc14dt=kinetic14(t,c14,a) %140度
T2=413.15;R=0.008314;
dc14dt=...
[-(a(1).*exp(-a(2)./(R.*T2)).*((c14(2)./100).^(a(3))).*c14(1).*0.926)
-(7.*a(1).*exp(-a(2)./(R.*T2)).*((c14(2)./100).^(a(3))).*c14(1).*0.926)];
function dc12dt=kinetic12(t,c12,a) %120度
T3=393.15;R=0.008314;
dc12dt=...
[-(a(1).*exp(-a(2)./(R.*T3)).*((c12(2)./100).^(a(3))).*c12(1).*0.907)
-(7.*a(1).*exp(-a(2)./(R.*T3)).*((c12(2)./100).^(a(3))).*c12(1).*0.907)];