反應(yīng)動力學(xué)，ode45求助！

天字：一：號

clear all
clc
DATA1=[
1    0.1578    0.3888    0.6673    0.7643;
4 0.1202    0.2235    1.1420    1.1303;
7 0.1136    0.2287    1.1510    1.1419;
10 0.1027    0.2172    1.2187    1.1922;
30  0.0701    0.2235    1.3194    1.2827;
60   0.0594    0.2472    1.2912    1.2790;
90   0.0538    0.3108    1.1387    1.1946];


k0 = [0.5  0.5];         % 參數(shù)初值                                                                                                
lb = [0  0 ];                   % 參數(shù)下限
ub = [+inf  +inf  ];    % 參數(shù)上限
x0 = [0.0625 0.09375 0 0]';                                                                                                                                          
yexp = DATA1(:,2:5);                  % yexp: 實驗數(shù)據(jù)[x1x2x3x4]
% 使用函數(shù)fmincon()進(jìn)行參數(shù)估計
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函數(shù)fmincon()估計得到的參數(shù)值為:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;


% 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函數(shù)lsqnonlin()估計得到的參數(shù)值為:\n')
Output


% 以函數(shù)fmincon()估計得到的結(jié)果為初值，使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計
k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的結(jié)果為初值，使用函數(shù)lsqnonlin()估計得到的參數(shù)值為:\n')
Output




% ------------------------------------------------------------------
function f = ObjFunc4Fmincon(k,x0,yexp)
tspan = [1 4 7 10 30 60 90];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,1) = x(:,1);
y(:,2:4) = x(:,2:4);
f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)   ...
    + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2);


% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = [1 4 7 10 30 60 90];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,1) = x(:,1);
y(:,2:4) = x(:,2:4);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f = [f1; f2; f3; f4];


% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
dxdt =1.17855*(1.6649E-05^2)*(k(1)*x(1)*x(2)-k(2)*x(3)*x(4));