★ ★ ★ ★ ★ 小木虫: 金币+0.5, 给个红包,谢谢回帖 fegg7502: 金币+4, 多谢交流 2012-07-09 08:07:15
这个其实论坛里有很多例子,参考就能写出来。
给你写了个CODE: function parafit
clear all;
t=[0 2 4 8 24 48];
y=[0.69 0.645 0.635 0.62 0.61 0.61];
y0=0.69;
% Nonlinear least square estimate using lsqnonlin()
beta0=0.5;
lb=[0];ub=[inf];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@Func,beta0,lb,ub,[],t,y,y0);
ci = nlparci(beta,residual,jacobian);
beta;
% result
fprintf('\n Estimated Parameters by Lsqnonlin():\n')
fprintf('\t k = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf(' The sum of the residual squares is: %.1e\n\n',sum(residual.^2))
% plot of fit results
tspan = [0 max(t)];
[tt yc] = ode45(@ModelEqs,tspan,y0,[],beta);
tc=linspace(0,max(t),200);
yca = spline(tt,yc,tc);
plot(t,y,'ro',tc,yca,'r-');
hold on
xlabel('Time');
ylabel('Concentration');
hold off
% =======================================
function f1 = Func(beta,t,y,y0) % Define objective function
tspan =t;
[tt yy] = ode45(@ModelEqs,tspan,y0,[],beta);
yc= spline(tt,yy,t);
f1=y-yc;
% ==================================
function dydt = ModelEqs(t,y,beta) % Model equations
dydt = -beta*(0.7-y).^(1/3)*(y-0.61);