非线性动力学matlab代码大全,matlab lsqnonlin函数进行非线性动力学参数拟合，求高手指点....

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