我在做一個動力學的方程擬合，用的是最小二乘法，下面是程序，不知道哪里有問題，結(jié)果偏差很大，求各位大神幫助！
原方程 K=k0*exp(-E/T)*PH^a*PS^b*P，知K T PH PS P,求 k0 E a b
function nihe
global T0 PH0 PS0 P0 K0
T0=[5014.5891 5014.5891 5014.5891 5180.8691 5180.8691 5180.8691 5347.1491 5347.1491 5347.1491];
PH0=[788732.39 592592.59 396039.6 790123.46 394366.19 594059.41 792079.21 591549.29 395061.73];
PS0=[20.19 13.28 7.1 17.71 10.1 10.65 14.2 15.15 8.85];
P0=[24220.91 29411.10 20587.77 20587.77 29411.10 24220.91 29411.10 20587.77 24220.91];
K0=[1047 1122 1233 1188 1225 1148 1815 1185 1125];
y=log(K0);
x1=T0
x2=log(PH0);
x3=log(PS0);
x4=log(P0);
y=y';
X=[ones(size(y)) x1' x2' x3' x4'];
[b bint]=regress(y,X,0.1);
k=exp(b(1));
E=b(2)
a=b(3);
c=b(4);
beta0=[k E a c];
lb=[];ub=[];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
     lsqnonlin(@ObjFunc,beta0,lb,ub);
ci=nlparci(beta,resid,jacobian)

% 殘差關于擬合值的殘差圖
K0c = Rate(beta,PH0,PS0,T0,P0);
plot(K0c,resid,'*')
xlabel('反應速率擬合值, torr s^-^1')
ylabel('殘差R, torr s^-^1')
refline(0,0)

% 參數(shù)辨識結(jié)果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tE = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\ta = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tc = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm)
% ----------------------------------------------------------------
function f=ObjFunc(beta)
global T0 PH0 PS0 P0 K0
f=K0-Rate(beta,T0,PH0,PS0,P0);
% ------------------------------------------------------------------
function K=Rate(beta,T,PH,PS,P)
K=beta(1).*exp(-beta(2)./T).*PH.^beta(3).*PS.^beta(4).*P;



運行結(jié)果如下
ci =

  1.0e+004 *

   -0.0000    0.0000
   -2.4270    3.5058
   -0.0000    0.0002
   -0.0002    0.0001

Estimated Parameters:
        k = 0.0001 ± 0.0015
        E = 5393.89 ± 29664.08
        a = 0.64 ± 1.05
        c = 0.64 ± 1.05
        The sum of the squares is: 2.2e+005 返回小木蟲查看更多