24小时热门版块排行榜    

查看: 341  |  回复: 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的回帖
相关版块跳转 我要订阅楼主 咔咔小男人 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 303求调剂 +3 蓝山月 2026-03-25 3/150 2026-03-25 10:04 by 雾散后相遇lc
[考研] 287求调剂 +9 晨昏线与星海 2026-03-19 10/500 2026-03-25 09:51 by 雾散后相遇lc
[考研] 0805 316求调剂 +6 大雪深藏 2026-03-18 6/300 2026-03-24 22:13 by peike
[考研] 08工学调剂 +12 用户573181 2026-03-20 17/850 2026-03-24 21:45 by 番茄炖了西红柿
[考研] 材料调剂 +6 匹克i 2026-03-23 6/300 2026-03-24 21:09 by greychen00
[有机交流] 有机合成求助 20+3 FENGSHUJEI 2026-03-23 5/250 2026-03-24 19:31 by 88817753
[考研] 0854 考研调剂 招生了!AI 方向 +5 pk3725069 2026-03-19 17/850 2026-03-24 17:30 by zhouxuan..
[考研] 求调剂 +5 林之夕 2026-03-24 5/250 2026-03-24 17:16 by dick_runner
[考博] 申博26年 +4 八6八68 2026-03-19 4/200 2026-03-24 15:49 by 小Ben呵呵
[考研] 341求调剂(一志愿湖南大学070300) +5 番茄头--- 2026-03-22 6/300 2026-03-23 23:45 by Txy@872106
[考研] 0854电子信息求调剂 324 +3 Promise-jyl 2026-03-23 3/150 2026-03-23 13:43 by wangkm
[考研] 276求调剂 +3 YNRYG 2026-03-21 4/200 2026-03-23 08:31 by 醉在风里
[考研] 293求调剂 +3 涛涛Wjt 2026-03-22 5/250 2026-03-22 22:21 by jiangpengfei
[考研] 269专硕求调剂 +6 金恩贝 2026-03-21 6/300 2026-03-22 14:31 by ColorlessPI
[考研] 材料与化工(0856)304求 B区 调剂 +3 邱gl 2026-03-21 3/150 2026-03-21 13:47 by lature00
[考研] 南昌大学材料专硕311分求调剂 +6 77chaselx 2026-03-20 6/300 2026-03-21 07:24 by JourneyLucky
[考研] 317求调剂 +5 申子申申 2026-03-19 9/450 2026-03-20 22:26 by JourneyLucky
[考研] 中南大学化学学硕337求调剂 +3 niko- 2026-03-19 6/300 2026-03-20 21:58 by luoyongfeng
[考研] 290求调剂 +7 ^O^乜 2026-03-19 7/350 2026-03-20 21:43 by JourneyLucky
[考研] 一志愿南理工085701环境302求调剂院校 +3 葵梓卫队 2026-03-20 3/150 2026-03-20 19:28 by zhukairuo
信息提示
请填处理意见