24小时热门版块排行榜    

查看: 315  |  回复: 2

咔咔小男人

金虫 (小有名气)

[求助] matlab求模型参数有错误,求指导 已有1人参与

以下的程序是我照葫芦画瓢改的,不知道哪里错了。。。。。。本人刚开始接触matlab,求大神指导啊!
想用方程dy/dx=c*(A-y)^a*(B-y)进行拟合,求出c和a,其中A=0.9,B=8.5。
CODE:
function KineticsEst1_Int_modified_by_Yuezhilan
clear all;clc
format long
tspan=[0,1,2,3,4,5,7,10,15,20,30,45,90,120];
yexp=[0,0.3425,0.61125,0.6125,0.67,0.8025,0.85,0.8375,0.93,0.955,0.95125,0.95375,0.95625,0.96125]';

k0=[1 1];   %%%请注意这里,初值的选取
y0=0.2;
lb=[0 0];
ub=[1000 10];   
yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t待拟合参数 k1 = %.6f\n',k(1))
fprintf('\t待拟合参数 k2 = %.6f\n',k(2))
fprintf(' \t残差平方和= %.6f\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\t决定系数R-Square = %.6f',R2);
figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');
yr=y-yexp;
figure(2)
plot(tspan(2:end),yr,'r*',[-1 15],[0 0]),axis([-1 15 -0.5 0.5]);
figure(3)
plot(yexp,y,'ro',[21 29],[21 29],'b-');

%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k)
ysim = Xsim(2:end);
size(ysim);
size(yexp);
f=ysim-yexp;
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
beta(1)=k(1);
beta(2)=k(2);
dydt =beta(1)*(0.9-y)^beta(2)*(7.5-y);

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

一个人知识的界限就是一个人所能达到的最大范围的界限。对一种事物的理解、认识深度决定了其利用这种事物上所能实现的价值的有限性。
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
咔咔小男人: 金币+20, ★★★★★最佳答案 2015-06-02 18:30:50
CODE:
function KineticsEst1_Int_modified_by_Yuezhilan
clear all;clc
format long
tspan=[0,1,2,3,4,5,7,10,15,20,30,45,90,120];
yexp=[0.3425,0.61125,0.6125,0.67,0.8025,0.85,0.8375,0.93,0.955,0.95125,0.95375,0.95625,0.96125]';

k0=[1 1];   %%%请注意这里,初值的选取
y0=0.2;
lb=[-inf -inf];
ub=[+inf +inf];   
yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t待拟合参数 k1 = %.6f\n',k(1))
fprintf('\t待拟合参数 k2 = %.6f\n',k(2))
fprintf(' \t残差平方和= %.6f\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
[ttt(2:end) y xexp]
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\t决定系数R-Square = %.6f',R2);
figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');



%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k);
ysim = Xsim(2:end);
size(ysim);
size(yexp);
f=ysim-yexp;
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
beta(1)=k(1);
beta(2)=k(2);
dydt =beta(1)*(0.9-y).^beta(2)*(7.5-y);

初值k0的取值很重要,自己再改改吧。
matlab求模型参数有错误,求指导
附图1.png

MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
2楼2015-06-01 18:52:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

咔咔小男人

金虫 (小有名气)

引用回帖:
2楼: Originally posted by 月只蓝 at 2015-06-01 18:52:09
function KineticsEst1_Int_modified_by_Yuezhilan
clear all;clc
format long
tspan=;
yexp=';

k0=;   %%%请注意这里,初值的选取
y0=0.2;
lb=;
ub=;   
yy=;

= ...
    lsqnonlin(@ObjFunc,k0,l ...

没想到作者来了,感谢您
一个人知识的界限就是一个人所能达到的最大范围的界限。对一种事物的理解、认识深度决定了其利用这种事物上所能实现的价值的有限性。
3楼2015-06-02 18:31:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 咔咔小男人 的主题更新
信息提示
请填处理意见