对含温度系数的化学反应速率方程进行拟合。测试发现用cftool的拟合函数进行拟合比调用lsqnonlin()要好用得多。

调试情况：

 a=8.188704E+013  b= -4799.760247
拟合残差 0.996347

输入校正系数
0.001
 A= 1.000208E+011  B= -842.079002  C= 2.499924
拟合残差 1.000000
输入校正系数，或按0键接受当前值，进行下一组拟合
0


拟合需要经验，只有估测值好的时候，可以使相关度趋近于1.

源代码：

function kneticMultiFit
%  功能：计算含温度系数的反应速率常数  
%       
%  说明：
%      1.拟合方程形式  K=A*exp(-E/RT)*(T/298)^n
%      2.令 x=1/T, a=A, b=-E/R, c=n
%      3.首先调用expl1Fit进行k=a*exp(b*x)拟合获得初值系数 a
%      4.调用expl2Fit进行k=a*exp(b*x)/(x*298)^c拟合得到系数 a b c
%      5.需要调用expl1Fit和expl2Fit文件
%
%  补充说明：根据经验，n>1时温度系数对高温段影响大，也就是1/T坐标轴靠近0的区域。
%  而通过调整系数A使它的初值接近真实值时，而其他参数B和C只需要置零，
%  也可使拟合效果比较理想。通过调整factor因子，使曲线在1/T坐标轴靠近0的区域拟合
%  度增大，便可以使最终残差趋近于1.
%      为了有好的收敛性，将expl2Fit中（第57行左右）ABC的拟合范围做了规定，当拟
%  合情况不理想，出现拟合值在边界上时，需要调整ABC的取值范围。默认为：
%           'Lower',[0 -Inf    0.3],'Upper',[Inf   0   3]
%  非线性拟合收敛性定义：
%      'MaxFunEvals',800,'MaxIter',600,'TolFun',1e-7,'TolX',1e-7
clear all,clc
% % -----测试计算系数---------------------
 T=[1000:50:2000];
 A=1e11;
 En=-7E3;
 R=8.314;
 KN=[A*exp(En/R./T).*(T/298).^2.5;A*exp(En/R./T).*(T/298).^1.5];
 
% %--------------------------------------

% % -----键盘输入T和k
% T=input('input Temperature T in [x1 x2 ... xn] form');
% k=input('input knetic k in [x1 x2 ... xn] form');
% %-----------------------------------------

% %----直接输入 T和k
% T=[];
% k=[];
% %-----------------------------------------

% %----读入文件中的数据
% data_fit
% %-------------------------------------------

fit_results=[];

for i=1:size(KN)

    k=KN(i,:);

% % 取T的倒数用于拟合计算
Tr=1./T;

% 用k=a*exp(b*x)拟合，得到org,g1
[org,g1]=expl1Fit(Tr,k);
fprintf('\n a=%E\t b= %f\n',org.a, org.b);
fprintf('拟合残差 %f \n',g1.rsquare);

% 输入校正系数 factor=1 0.1 10 etc.
factor=input('输入校正系数\n');
% 选用不同的 a系数
af1=org.a*factor;

% 定义其他初值
bf=0;
nf=0;

% % 用k=a*exp(b*x)/(x*298)^c拟合，得到f,f_g
[f1,f_g1]=expl2Fit(Tr,k,af1,bf,nf);
fprintf('\n A= %E\t B= %f\t C= %f\n',f1.a, f1.b, f1.c);
fprintf('拟合残差 %f \n',f_g1.rsquare);

% 提示修改校正系数
factor=input('输入校正系数，或按0键接受当前值，进行下一组拟合\n');

while factor~=0
    close all
    af1=org.a*factor;
    [f1,f_g1]=expl2Fit(Tr,k,af1,bf,nf);
    fprintf('\n A= %E\t B= %f\t C= %f\n',f1.a, f1.b, f1.c);
    fprintf('拟合残差 %f \n',f_g1.rsquare);
    factor=input('输入校正系数，或按0键接受当前值，进行下一组拟合\n');
end
   
    fit_results=[fit_results;f1.a f1.b f1.c f_g1.rsquare];

i=i+1;
close all
end

final_data=[KN fit_results];
save fit final_data

 function [cf_,g] =expl1Fit(Tr,k)

% Set up figure to receive data sets and fits
figure
f_ = clf;
figure(f_);
set(f_,'Units','Pixels','Position',[658 269 688 485]);
% Line handles and text for the legend.
legh_ = [];
legt_ = {};
% Limits of the x-axis.
xlim_ = [Inf -Inf];
% Axes for the plot.
ax_ = axes;
set(ax_,'Units','normalized','OuterPosition',[0 0 1 1]);
set(ax_,'Box','on');
axes(ax_);
hold on;

% --- Plot data that was originally in data set "k vs. Tr"
Tr = Tr(:);
k = k(:);
h_ = line(Tr,k,'Parent',ax_,'Color',[0.333333 0 0.666667],...
    'LineStyle','none', 'LineWidth',1,...
    'Marker','.', 'MarkerSize',12);
xlim_(1) = min(xlim_(1),min(Tr));
xlim_(2) = max(xlim_(2),max(Tr));
legh_(end+1) = h_;
legt_{end+1} = 'k vs. Tr';

% Nudge axis limits beyond data limits
if all(isfinite(xlim_))
    xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
    set(ax_,'XLim',xlim_)
else
    set(ax_, 'XLim',[0.000495, 0.001005]);
end

% --- Create fit "fit 1"
ok_ = isfinite(Tr) & isfinite(k);
if ~all( ok_ )
    warning( 'GenerateMFile:IgnoringNansAndInfs',...
        'Ignoring NaNs and Infs in data.' );
end
st_ = [101312471405283.05 -3686.6509320548034 ];
ft_ = fittype('exp1');

% Fit this model using new data
[cf_,g] = fit(Tr(ok_),k(ok_),ft_,'Startpoint',st_);
% Alternatively uncomment the following lines to use coefficients from the
% original fit. You can use this choice to plot the original fit against new
% data.
%    cv_ = { 159605243069299.56, -4424.3830760570581};
%    cf_ = cfit(ft_,cv_{:});

% Plot this fit
h_ = plot(cf_,'fit',0.95);
set(h_(1),'Color',[1 0 0],...
    'LineStyle','-', 'LineWidth',2,...
    'Marker','none', 'MarkerSize',6);
% Turn off legend created by plot method.
legend off;
% Store line handle and fit name for legend.
legh_(end+1) = h_(1);
legt_{end+1} = 'fit 1';

% --- Finished fitting and plotting data. Clean up.
hold off;
% Display legend
leginfo_ = {'Orientation', 'vertical', 'Location', 'NorthEast'};
h_ = legend(ax_,legh_,legt_,leginfo_{:});
set(h_,'Interpreter','none');
% Remove labels from x- and y-axes.
xlabel(ax_,'');
ylabel(ax_,'');

function [cf_,g]= expl2Fit(Tr,k,af,bf,nf)
% Set up figure to receive data sets and fits
figure
f_ = clf;
figure(f_);
set(f_,'Units','Pixels','Position',[1021 249 688 485]);
% Line handles and text for the legend.
legh_ = [];
legt_ = {};
% Limits of the x-axis.
xlim_ = [Inf -Inf];
% Axes for the plot.
ax_ = axes;
set(ax_,'Units','normalized','OuterPosition',[0 0 1 1]);
set(ax_,'Box','on');
axes(ax_);
hold on;

% --- Plot data that was originally in data set "k vs. Tr"
Tr = Tr(:);
k = k(:);
h_ = line(Tr,k,'Parent',ax_,'Color',[0.333333 0 0.666667],...
    'LineStyle','none', 'LineWidth',1,...
    'Marker','.', 'MarkerSize',12);
xlim_(1) = min(xlim_(1),min(Tr));
xlim_(2) = max(xlim_(2),max(Tr));
legh_(end+1) = h_;
legt_{end+1} = 'k vs. Tr';

% Nudge axis limits beyond data limits
if all(isfinite(xlim_))
    xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
    set(ax_,'XLim',xlim_)
else
    set(ax_, 'XLim',[0.000495, 0.001005]);
end

% --- Create fit "fit 1"
fo_ = fitoptions('method','NonlinearLeastSquares','Lower',[0 -Inf    0.3],'Upper',[Inf   0   3],'MaxFunEvals',800,'MaxIter',600,'TolFun',1e-7,'TolX',1e-7);
ok_ = isfinite(Tr) & isfinite(k);
if ~all( ok_ )
    warning( 'GenerateMFile:IgnoringNansAndInfs',...
        'Ignoring NaNs and Infs in data.' );
end
st_ = [af bf nf ];
set(fo_,'Startpoint',st_);
ft_ = fittype('a*exp(b*x)/(298*x)^c',...
    'dependent',{'y'},'independent',{'x'},...
    'coefficients',{'a', 'b', 'c'});

% Fit this model using new data
[cf_,g] = fit(Tr(ok_),k(ok_),ft_,fo_);
% Alternatively uncomment the following lines to use coefficients from the
% original fit. You can use this choice to plot the original fit against new
% data.
%    cv_ = { 2524104546234.3926, -1993.4420796170032, 1.5508056550249052};
%    cf_ = cfit(ft_,cv_{:});

% Plot this fit
h_ = plot(cf_,'fit',0.95);
set(h_(1),'Color',[1 0 0],...
    'LineStyle','-', 'LineWidth',2,...
    'Marker','none', 'MarkerSize',6);
% Turn off legend created by plot method.
legend off;
% Store line handle and fit name for legend.
legh_(end+1) = h_(1);
legt_{end+1} = 'fit 1';

% --- Finished fitting and plotting data. Clean up.
hold off;
% Display legend
leginfo_ = {'Orientation', 'vertical', 'Location', 'NorthEast'};
h_ = legend(ax_,legh_,legt_,leginfo_{:});
set(h_,'Interpreter','none');
% Remove labels from x- and y-axes.
xlabel(ax_,'');
ylabel(ax_,'');

天乐树

2011.10.09