24小时热门版块排行榜    

查看: 2886  |  回复: 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的回帖

wtb1987

新虫 (初入文坛)

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

zt_chem

木虫 (正式写手)

能不能将速率常数k写成k=A*EXP(-Ea/(R*T))代入常微分方程组中来直接估计指前因子和活化能呢?
2楼2011-01-10 22:33:17
已阅   回复此楼   关注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的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 建议基金发布提前给出明确的时间点 +11 kulium 2026-08-21 14/700 2026-08-24 02:12 by 流流伤
[基金申请] 范进中举一文的中心思想 +4 炎黄贵胄 2026-08-22 5/250 2026-08-24 02:09 by 流流伤
[基金申请] 2026国自然函评费到账 +16 羊腰板 2026-08-21 17/850 2026-08-23 23:02 by tianxiaochun
[基金申请] 什么时候开奖? +9 CrisMessi 2026-08-18 10/500 2026-08-23 12:03 by 丶昵称占用
[基金申请] 今天放榜吗? +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 93BebMhtakh前后11位开头都是大写 +7 且听虎啸 2026-08-17 8/400 2026-08-22 21:55 by 医学老男孩
[基金申请] 只有每年这种时候来逛逛小木虫 +24 yaoyewhu2008 2026-08-20 26/1300 2026-08-22 17:43 by kammury
[基金申请] filecode,4个jtjc了 +13 ziyangfang 2026-08-19 16/800 2026-08-22 17:08 by WH3796
[基金申请] 人气不行了 +8 fansofjerry 2026-08-21 8/400 2026-08-22 16:30 by zyqchem
[基金申请] 今天基金会出结果吗?20260819 +16 kkkl_v 2026-08-19 17/850 2026-08-22 16:12 by 阿布Abu
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +16 医学老男孩 2026-08-20 20/1000 2026-08-21 21:13 by Ldrop2023
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 感觉是下周放榜了 +7 angus9576 2026-08-17 12/600 2026-08-21 13:38 by weiyin
[基金申请] 我面上完蛋了 +7 且听虎啸 2026-08-20 8/400 2026-08-21 12:31 by 酷酷墨镜
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
信息提示
请填处理意见