24小时热门版块排行榜    

查看: 2913  |  回复: 5

wtb1987

新虫 (初入文坛)

[交流] 【求助】催化反应动力学 已有5人参与

大家好,有研究催化反应动力学的吗?大家可以交流一下,从实验到数值拟合个都可以,我做这方面也很长时间了,但是有很多东西搞不明白。如,下面是对某个反应进行的数值模拟,使用的软件是matlab,但是求出来的动力学数值有很大的置信区间,想问大家一下,这是我做的实验数据的问题?还是程序的问题?希望做动力学的高手给指点一下。万分感谢!!!!
function HZSM_8
clear all;clc
expdata=[4.28 0.0724 0.0614 0.0075 0.0013 0.0021;
         8.56 0.1126 0.0859 0.0164 0.0039 0.0064;
         16.69 0.1529 0.1034 0.0278 0.0079 0.0138;
         33.76 0.1928 0.1134 0.0412 0.0132 0.0251;
         69.72 0.2305 0.1158 0.0555 0.0192 0.04;];
yexp=expdata(:,2:6);
x0=[0 0 0 0 0];
k0=[1.862 0.578 1.321 0.8 9.813 1.018 0.2409 0.626 4.5301 18.174 40.311 314.1];
options=optimset('Algorithm','Levenberg-Marquardt','LargeScale','off');
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@objfunc,k0,[],[],options,x0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('\t k1=%.4f±%.4f\n',k(1),ci(1,2)-k(1))
fprintf('\t k2=%.4f±%.4f\n',k(2),ci(2,2)-k(2))
fprintf('\t k3=%.4f±%.4f\n',k(3),ci(3,2)-k(3))
fprintf('\t k4=%.4f±%.4f\n',k(4),ci(4,2)-k(4))
fprintf('\t k5=%.4f±%.4f\n',k(5),ci(5,2)-k(5))
fprintf('\t k6=%.4f±%.4f\n',k(6),ci(6,2)-k(6))
fprintf('\t k7=%.4f±%.4f\n',k(7),ci(7,2)-k(7))
fprintf('\t k8=%.4f±%.4f\n',k(8),ci(8,2)-k(8))
fprintf('\t k9=%.4f±%.4f\n',k(9),ci(9,2)-k(9))
fprintf('\t k10=%.4f±%.4f\n',k(10),ci(10,2)-k(10))
fprintf('\t k11=%.4f±%.4f\n',k(11),ci(11,2)-k(11))
fprintf('\t k12=%.4f±%.4f\n',k(12),ci(12,2)-k(12))
fprintf('\t the sum of the squares is:%.1e\n\n',resnorm)
function f=objfunc(k,x0,yexp)
%使用四阶龙格库塔法求取常微分方程组的值
tspan=[0 4.28 8.56 16.69 33.76 69.72];
[t,x]=ode45(@equations,tspan,x0,[],k)
y=x(2:6,;
%目标函数的定义
f1=y(:,1)-yexp(:,1);
f2=y(:,2)-yexp(:,2);
f3=y(:,3)-yexp(:,3);
f4=y(:,4)-yexp(:,4);
f5=y(:,5)-yexp(:,5);
f=[f1;f2;f3;f4;f5];
function dxdt=equations(t,x,k)
P=1;
PT=(1-x(1))*P/(6+x(1));
Ppx=x(2)*P/(6+x(2));
PO=x(4)*P/(6+x(1));
PMx=x(3)*P/(6+x(1));
PTMB=x(5)*P/(6+x(5));
PM=[0.5-(x(1)+x(5))]/(6+x(1));
Z=(1+k(7)*PT+k(8)*PM+k(9)*Ppx+k(10)*PMx+k(11)*PO+k(12)*PTMB)^2;
dx1dt=k(1)*PT*PM/Z;
dx2dt=[k(1)*PT*PM-k(2)*(Ppx-PMx/2.2)-k(4)*Ppx*PM]/Z;
dx3dt=[k(2)*(Ppx-PMx/2.2)-k(3)*(PMx-PO/0.4)-k(5)*PMx*PM]/Z;
dx4dt=[k(3)*(PMx-PO/0.4)-k(6)*PO*PM]/Z;
dx5dt=(k(4)*Ppx*PM+k(5)*PMx*PM+k(6)*PO*PM)/Z;
dxdt=[dx1dt;dx2dt;dx3dt;dx4dt;dx5dt];

结果是:
k1=1.9174±58287.8297
         k2=0.5759±17506.6005
         k3=1.4271±43384.6026
         k4=0.8909±27081.9850
         k5=9.4426±287049.8962
         k6=2.5515±77564.2240
         k7=-0.0804±227820.0664
         k8=1.4310±276004.5168
         k9=5.7924±106309.8775
         k10=17.1832±259330.6709
         k11=37.7388±551649.7103
         k12=314.3976±4536110.8346
         the sum of the squares is:6.2e-007
回复此楼

» 猜你喜欢

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

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zt_chem

木虫 (正式写手)

能不能将速率常数k写成k=A*EXP(-Ea/(R*T))代入常微分方程组中来直接估计指前因子和活化能呢?
2楼2011-01-10 22:33:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wtb1987

新虫 (初入文坛)

我认为可以这么做,但是这样的话更增加了方程组的非线性,本来这样已经12个未知参数了,如果直接代入求解将会变得更加复杂。
3楼2011-01-12 09:31:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dongyan-608

银虫 (小有名气)

有道理
相信自己
4楼2011-09-06 09:45:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cyzeng

金虫 (小有名气)


小木虫: 金币+0.5, 给个红包,谢谢回帖
我也是碰到这样的问题,该怎么办啊
5楼2012-09-07 09:57:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

spiderone

木虫 (小有名气)


小木虫: 金币+0.5, 给个红包,谢谢回帖
楼上的诸位 此问题怎么解决的啊 !请指教~~~
6楼2013-12-02 16:35:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wtb1987 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-09-02 00:19 by xiangy672
[基金申请] 学科评审组评审是指会评吗? +4 瞬息宇宙 2026-08-31 4/200 2026-09-01 14:58 by jiaoxg
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 基金不中,共勉 +12 eulota 2026-08-26 12/600 2026-08-31 08:42 by ZJTJZ
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 看板上这么多中的,有点像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 北京莱茵编辑
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
信息提示
请填处理意见