应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求助用matlab三次样条法算一个反应的反应速率
142/2返回列表上一页12
查看: 2499  |  回复: 13
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

justin534

新虫 (初入文坛)
???????:
10?: Originally posted by ????? at 2014-03-29 15:11:44
?????????????????????????????????
?????????????dingd ??1stopt???????????????????...

???????,???!????????}??!??????????????????W?,?F????????3??,??6??k1 k2 k3 k4 k5 k6!
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC
???????????????a!???????RUN,??????????!


function ODE_parafit2
%-------data
yexp=[5.14  12.76  45.41  21.52  3.31;
5.36  13.20  46.94  20.90  3.17;
5.51  13.47  48.68  20.01  3.08;
5.78  14.03  48.16  19.34  3.04;
6.07  14.33  47.80  18.65  2.94]';
x0=[93.7 6.3 0 0 0];
k0 = [1 1 1 1 1 1];
lb = [0 0 0 0 0 0];
ub = [1 1 1 1 1 1]*1e5;
tspan = [0 38.46 50 62.5 100 166.666 250 500];
%-----------------
[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);
fprintf('\n\n??ú???fmincon()?????????????:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',fval)

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n??ú???lsqnonlin()?????????????:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',resnorm)

%%  plot figure
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i,=x(end,;
end
Ypre=fy(;
Yexp=yexp(;
plot(Ypre,'or',Yexp,'b.')
legend('prediction','experiment')
%----------fmincon fun-------------------
function f = ObjFuncFmincon(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i)=sum((yexp(i,-x(end,).^2);
end
f =  sum(fy) ;
% -------------------lsqnonlin fun-------------------
function f = ObjFunc(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666];
for i=1:5
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i,=yexp(i,-x(end,;
end
f =fy(;
%-----------odefun--------------------------------
function dYdt = KineticsEqs(t,y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);
DEC=y(1);PA=y(2);DPC=y(3);EA=y(4);EPC=y(5);
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA;
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA;
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC;
dy1dt=-R1+R3;
dy2dt=-R1-R2;
dy3dt=R2+R3;
dy4dt=R1+R2;
dy5dt=R1-R2-R3;
dYdt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt];
赞一下 回复此楼
11楼2014-04-24 14:50:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

justin534

新虫 (初入文坛)
引用回帖:
10楼: Originally posted by 月只蓝 at 2014-03-29 15:11:44
可能是方程的问题，也可能是实验数据的问题。
你可以请实习版主dingd 用1stopt软件做一下，来验证一下。...

sorry,剛剛的編碼有誤,這個是OK的編碼!還有data是隨意打的,所以只是請你幫我看可不可以RUN!

function ODE_parafit2
%-------data
yexp=[5.14  12.76  45.41  21.52  3.31 1 1;
5.36  13.20  46.94  20.90  3.17 1 1;
5.51  13.47  48.68  20.01  3.08 1 1;
5.78  14.03  48.16  19.34  3.04 1 1;
6.07  14.33  47.80  18.65  2.94 1 1]';
x0=[93.7 6.3 0 0 0];
k0 = [1 1 1 1 1 1];
lb = [0 0 0 0 0 0];
ub = [1 1 1 1 1 1]*1e5;
tspan = [0 38.46 50 62.5 100 166.666 250 500];
%-----------------
[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);
fprintf('\n\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',fval)

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk6 = %16f \n',k(6))
fprintf('The sum of squares is: %.16f \n\n',resnorm)

%%  plot figure
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i,=x(end,;
end
Ypre=fy(;
Yexp=yexp(;
plot(Ypre,'or',Yexp,'b.')
legend('prediction','experiment')
%----------fmincon fun-------------------
function f = ObjFuncFmincon(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666 250 500];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i)=sum((yexp(i,-x(end,).^2);
end
f =  sum(fy) ;
% -------------------lsqnonlin fun-------------------
function f = ObjFunc(k,x0,yexp)
tspan = [0 38.46 50 62.5 100 166.666];
for i=1:8
[t x] = ode45(@KineticsEqs,tspan,x0,[],k);   
fy(i,=yexp(i,-x(end,;
end
f =fy(;
%-----------odefun--------------------------------

function dydt = KineticsEqs(t,y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);
DEC=y(1);PA=y(2);DPC=y(3);EA=y(4);EPC=y(5);
R1=k(1)*DEC*PA-(k(1)/k(2))*EPC*EA;
R2=k(3)*EPC*PA-(k(3)/k(4))*DPC*EA;
R3=k(5)*EPC.^2-(k(5)/k(6))*DPC*DEC;
dy1dt=-R1+R3;
dy2dt=-R1-R2;
dy3dt=R2+R3;
dy4dt=R1+R2;
dy5dt=R1-R2-R3;
dydt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt];
赞一下 回复此楼
12楼2014-04-24 15:11:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

WULONGLOVE

铁虫 (初入文坛)
引用回帖:
2楼: Originally posted by 月只蓝 at 2014-03-26 20:00:43
比较粗糙的办法是求数值导数，根据t CA的数据即可求出，如果需要MATLAB代码参见：http://muchong.com/bbs/viewthread.php?tid=7020495；
精确点的方法是先拟合出动力学方程中的各参数，以第一组数据（即T=325.15K这组 ...

大神， lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp); 中的y0代表啥意思啊
赞一下 回复此楼
13楼2014-12-03 17:15:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
引用回帖:
13楼: Originally posted by WULONGLOVE at 2014-12-03 17:15:12
大神， lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp); 中的y0代表啥意思啊...

预估的初值
回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
14楼2014-12-03 18:09:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 liqianmelody 的主题更新
142/2返回列表上一页12
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的？ +15 爱看书的可乐 2026-08-26 17/850 2026-08-27 03:31 by xxniao123
[硕博家园] 售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急 +3 LIbGuocjEEYw 2026-08-26 4/200 2026-08-27 02:01 by Ie9AyIAvGbvs
[基金申请] 看板上这么多中的，有点像50人群里49个人都是骗子的那种感觉…… +3 a089 2026-08-26 3/150 2026-08-26 23:21 by young90
[基金申请] 我不理解！ +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 怎么查啊 +5 huang1991js 2026-08-26 5/250 2026-08-26 16:09 by Equinoxhua
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 科研孤儿太难了 +21 我4大白菜 2026-08-20 22/1100 2026-08-26 11:29 by nlgza
[基金申请] 怎么看青基中了没有啊 +4 叶九微 2026-08-26 4/200 2026-08-26 10:52 by xiacongxu
[职场人生] 学生撞见辅导员送外卖，第二天全班都沉默了 +3 鲸鱼融金_浙江_� 2026-08-22 3/150 2026-08-26 09:01 by zzuzxg
[基金申请] 在坚冰还盖着北海的时候，我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了！？ +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 今日不放榜？网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[基金申请] 如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[教师之家] 跳槽后在研项目怎么办？ +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 时间戳今天，20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧？ +4 哈哈蛤？ 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
信息提示
关闭
请填处理意见
关闭