24小时热门版块排行榜    

查看: 417  |  回复: 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的回帖
相关版块跳转 我要订阅楼主 咔咔小男人 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 能否申诉? +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[文学芳草园] 梦想 +5 myrtle 2026-08-26 5/250 2026-08-31 14:40 by hahaboy
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 11:51 by xiangy672
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 国社科又开始会评了,不知道这次命运如何 +6 雨打竹帘 2026-08-30 10/500 2026-08-31 09:29 by huixian257
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +12 Kittylucky 2026-08-27 13/650 2026-08-29 00:04 by superceng
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
信息提示
请填处理意见