24小时热门版块排行榜    

查看: 1125  |  回复: 11
【悬赏金币】回答本帖问题,作者ykai_1991将赠送您 100 个金币

ykai_1991

金虫 (初入文坛)

[求助] 1stOpt编写常微分方程组计算反应速率问题已有2人参与

还请各位大神半忙查看下,我也是来学习的哈。另外烦请dingd帮忙查看下代码是否合理并计算下结果。

反应过程和微分方程见附图。

在编写中,数字代表不同的链长,物质是随着反应增长的。

比较麻烦的是微分方程中的多组分的加和,并且每个物质都有个方程,都写出来比较费时间,不知道有没有简便点的方法。
CODE:
本人这么写代码:([b]需要求每步反应速率k值(即为ka和kc),参数很多,不知能否计算出来;并且提供三组数据[/b])
Parameters ka(10), kc1(10), kc2(10), kc3(10), kc4(10), kc5(10), kc6(10), kc7(10), kc8(10), kc9(10), kc10(10) ;       [b]这里的kc有方法简化下么?不然就要写10个了[/b]
Variable t, P(10), Q(10), F;   [b]P1~P10也可以这么定义变量么,直接用P(10)?[/b]
Plot t[x], P1, P2, Q1, Q2, F;   [b]不同组分的含量不一样,在一张图中会画出来看不出来?如果不同物质分开作图,怎么写?[/b]
Plot t[x], P3, P4, Q3, Q4;
ODEFunction
           P1'=-ka1*P1*F-P1*Sum(j=1:10)(kc1[j]*Q[j]);
           P2'=-ka2*P2*F-P2*Sum(j=1:10)(kc2[j]*Q[j])+Sum(j=1:(2-1))(kc2[j]*P2*Q[2-j]);   [b]式子中的数字2都和P2的2有关系,如果P2变成P3,那么其他的2也都为3,那么针对P3'到P10',有没有更为方便的方法?下面再写一个P3'的,要一直写到P10'的
           P3'=-ka3*P3*F-P3*Sum(j=1:10)(kc3[j]*Q[j])+Sum(j=1:(3-1))(kc3[j]*P3*Q[3-j]);[/b]
           同样类似递推编写的方法也适用于Q1'到Q10',因为都只是角标不一样,加和的物质产生变化而已,在matlab中就可以定义出来比较方便,针对Q1'如下:
           Q1'=-Q1*Sum(j=1:10)(kc1[j]*P[j])+ka1*P1*F;   [b]同理,Q2'=-Q2*Sum(j=1:10)(kc2[j]*P[j])+ka2*P2*F; 也要一直写到Q10'[/b]
           F'=-F*Sum(j=1:10)(ka[j]*P[j]);

Data t, P1, Q1, P2;
10, 30, 40, 50, 60, 75, 110, 150, 230;
10, 30, 40, 50, 60, 75, 110, 150, 230;
6.288097876, 4.91765691, 3.878352807, 3.298605506, 2.850889397, 2.241683203, 1.781706784, 0.910037534, 0.693987349;
0.257670552, 0.312407281, 0.244993608, 0.189709656, 0.133296065, 0.0803651, 0.033400217, 0.010654615, 0.005606342;
0.236493005, 0.607274759, 0.615323542, 0.623870854, 0.61234489, 0.562467845, 0.538393919, 0.347783178, 0.305460229;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 120, 170, 230, 305;
5.882896195, 4.593250328, 3.701714247, 2.910582239, 2.571710069, 2.416218035, 1.920578469, 1.68785903, 1.416030469, 1.287574203;
0.255794223, 0.297297319, 0.219253251, 0.161355927, 0.121348782, 0.081854014, 0.039094897, 0.016390052, 0.008676848, 0.005710866;
0.203988176, 0.534247028, 0.637609417, 0.634725031, 0.638037303, 0.569498973, 0.479824153, 0.477902906, 0.424218287, 0.411113563;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 130, 190, 237;
5.879605129, 4.718556629, 3.450891407, 3.045516702, 2.681262558, 2.428241655, 1.569209161, 1.234288027, 1.237674719;
0.22471112, 0.268468006, 0.177155403, 0.1332699, 0.094942111, 0.071844393, 0.019781466, 0.008520798, 0.00503865;
0.226510416, 0.569687824, 0.596106833, 0.599362637, 0.580304305, 0.666507432, 0.464066174, 0.389844442, 0.423282218;


最后再附上编写过的微分方程组部分的matlab程序,便于您对照。
function dydt=KineticsEqs(t, y, ka, kc)
m = 10;      
n = 10;
for i=1:m    %-----------------------将P赋值给y-----------------------------
    P(i) = y(i);
end
for j=1:n    %-----------------------将Q赋值给y-----------------------------
    Q(j) = y(m+j);
end
F = y(m+n+1);%-----------------------将F赋值给y-----------------------------
%----------------------------以下为动力学微分方程----------------------------
for j=1:n
    polynomial_1(j)=kc(1,j)*Q(j);
end
r(1) = -(  ka(1)*P(1)*F+P(1)*sum(polynomial_1)  ); %-------微分方程1--------

for i=1:m
    for j=1:n
    polynomial_2(j)=kc(i,j)*Q(j);
    end
    polysum_2(i)=sum(polynomial_2);
end
for i=2:m
    for j=1:i-1
         polynomial_3(j)=kc(i,j)*P(j)*Q(i-j);
    end
    polysum_3(i)=sum(polynomial_3);
end
r(i) = -(  ka(i)*P(i)*F+P(i)*polysum_2(i)-polysum_3(i)  );

for i=1:m
    for j=1:n
        polynomial_4(j)=kc(i,j)*P(j);
    end
    r(i+n) = -(  Q(i)*sum(polynomial_4)-ka(i)*P(i)*F  );%-----微分方程3-----
end

for j=1:n
        polynomial_5(j)=ka(j)*P(j);
end
r(m+n+1) = -(  F*sum(polynomial_5)  );             %-------微分方程4--------

dydt = r';

1stOpt编写常微分方程组计算反应速率问题
反应方程


1stOpt编写常微分方程组计算反应速率问题-1
微分方程组

[ Last edited by dbb627 on 2015-1-15 at 23:45 ]
回复此楼

» 猜你喜欢

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

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

ykai_1991

金虫 (初入文坛)

本人重新写下代码,有些错误的地方哈
2楼2015-01-16 13:27:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

ben_ladeng: 屏蔽内容 2015-01-17 13:10:30
ben_ladeng: 取消屏蔽 2015-01-17 13:11:09
本人重新写下代码,有些错误的地方哈
  对于两种下角标的,这里的kc有方法简化下么?不然就要写10个了
  P1~P10也可以这样直接用P(1:10)定义变量么?

Parameters ka(1:10), kc1(1:10), kc2(1:10), kc3(1:10), kc4(1:10), kc5(1:10), kc6(1:10), kc7(1:10), kc8(1:10), kc9(1:10), kc10(1:10) ;      
Variable t, P(1:10), Q(1:10), F;      
Plot t[x], P1, P2, Q1, Q2, F;   
Plot t[x], P3, P4, Q3, Q4;
Plot t[x], P5, P6, Q5, Q6;
Plot t[x], P7, P8, Q7, Q8;
ODEFunction
           P1'=-ka1*P1*F-P1*Sum(j=1:10)(kc1[j]*Q[j]);
           P2'=-ka2*P2*F-P2*Sum(j=1:10)(kc2[j]*Q[j])+Sum(j=2-1))(kc2[j]*P2*Q[2-j]);
           P3'=-ka3*P3*F-P3*Sum(j=1:10)(kc3[j]*Q[j])+Sum(j=3-1))(kc3[j]*P3*Q[3-j]);
           P4'=-ka4*P4*F-P4*Sum(j=1:10)(kc4[j]*Q[j])+Sum(j=4-1))(kc4[j]*P4*Q[4-j]);
           P5'=-ka5*P5*F-P5*Sum(j=1:10)(kc5[j]*Q[j])+Sum(j=5-1))(kc5[j]*P5*Q[5-j]);
           P6'=-ka6*P6*F-P6*Sum(j=1:10)(kc6[j]*Q[j])+Sum(j=6-1))(kc6[j]*P6*Q[6-j]);
           P7'=-ka7*P7*F-P7*Sum(j=1:10)(kc7[j]*Q[j])+Sum(j=7-1))(kc7[j]*P7*Q[7-j]);
           P8'=-ka8*P8*F-P8*Sum(j=1:10)(kc8[j]*Q[j])+Sum(j=8-1))(kc8[j]*P8*Q[8-j]);
           P9'=-ka9*P9*F-P9*Sum(j=1:10)(kc9[j]*Q[j])+Sum(j=9-1))(kc9[j]*P9*Q[9-j]);
           P10'=-ka10*P10*F-P10*Sum(j=1:10)(kc10[j]*Q[j])+Sum(j=10-1))(kc10[j]*P10*Q[10-j]);
           Q1'=-Q1*Sum(j=1:10)(kc1[j]*P[j])+ka1*P1*F;
           Q2'=-Q2*Sum(j=1:10)(kc2[j]*P[j])+ka2*P2*F;
           Q3'=-Q3*Sum(j=1:10)(kc3[j]*P[j])+ka3*P3*F;
           Q4'=-Q4*Sum(j=1:10)(kc4[j]*P[j])+ka4*P4*F;
           Q5'=-Q5*Sum(j=1:10)(kc5[j]*P[j])+ka5*P5*F;
           Q6'=-Q6*Sum(j=1:10)(kc6[j]*P[j])+ka6*P6*F;
           Q7'=-Q7*Sum(j=1:10)(kc7[j]*P[j])+ka7*P7*F;
           Q8'=-Q8*Sum(j=1:10)(kc8[j]*P[j])+ka8*P8*F;
           Q9'=-Q9*Sum(j=1:10)(kc9[j]*P[j])+ka9*P9*F;
           Q10'=-Q10*Sum(j=1:10)(kc10[j]*P[j])+ka10*P10*F;
           F'=-F*Sum(j=1:10)(ka[j]*P[j]);

Data t, P1, Q1, P2;
3楼2015-01-16 13:32:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

ben_ladeng: 屏蔽内容 2015-01-17 13:11:02
本帖内容被屏蔽

4楼2015-01-16 20:58:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

除了时间t是min的单位,其他的物质都是mol/L的单位。
5楼2015-01-16 20:59:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

数据修改,补充了0时刻的数据

Data t, P1, Q1, P2, F;
0, 10, 30, 40, 50, 60, 75, 110, 150, 230;
6.569095641, 6.288097876, 4.91765691, 3.878352807, 3.298605506, 2.850889397, 2.241683203, 1.781706784, 0.910037534, 0.693987349;
0, 0.257670552, 0.312407281, 0.244993608, 0.189709656, 0.133296065, 0.0803651, 0.033400217, 0.010654615, 0.005606342;
0, 0.236493005, 0.607274759, 0.615323542, 0.623870854, 0.61234489, 0.562467845, 0.538393919, 0.347783178, 0.305460229;
5.255465, 3.817388062, 2.897843059, 2.434791186, 2.049634402, 1.61835643, 1.169440852, 0.623491234, 0.329066778, 0.158982809;

Data t, P1, Q1, P2, F;
0, 10, 30, 50, 60, 70, 85, 120, 170, 230, 305;
6.735750308, 5.882896195, 4.593250328, 3.701714247, 2.910582239, 2.571710069, 2.416218035, 1.920578469, 1.68785903, 1.416030469, 1.287574203;
0, 0.255794223, 0.297297319, 0.219253251, 0.161355927, 0.121348782, 0.081854014, 0.039094897, 0.016390052, 0.008676848, 0.005710866;
0, 0.203988176, 0.534247028, 0.637609417, 0.634725031, 0.638037303, 0.569498973, 0.479824153, 0.477902906, 0.424218287, 0.411113563;
5.051834547, 4.008120065, 2.876234251, 2.127807567, 1.792536461, 1.454135282, 1.123105548, 0.750702561, 0.437818947, 0.251544912, 0.126883404;

Data t, P1, Q1, P2, F;
0, 10, 30, 50, 60, 70, 85, 130, 190, 237;
6.911080972, 5.879605129, 4.718556629, 3.450891407, 3.045516702, 2.681262558, 2.428241655, 1.569209161, 1.234288027, 1.237674719;
0, 0.22471112, 0.268468006, 0.177155403, 0.1332699, 0.094942111, 0.071844393, 0.019781466, 0.008520798, 0.00503865;
0, 0.226510416, 0.569687824, 0.596106833, 0.599362637, 0.580304305, 0.666507432, 0.464066174, 0.389844442, 0.423282218;
4.837603126, 3.249855677, 2.492115606, 1.848010961, 1.524554737, 1.258887802, 1.025540648, 0.527930211, 0.318209995, 0.194771833;
6楼2015-01-17 10:44:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
光定义的参数就有110个,再加上P3至 P10, Q3至Q10都没有实际数据,又要增加16个参数,也即126个参数,还是微分方程问题,计算太费时间了,一般老点的PC很吃力了,没法耗几个甚至十几个小时替你跑。
7楼2015-01-18 20:24:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

CelestialCYJ

木虫 (小有名气)

【答案】应助回帖

你P3~10,Q2~10在0时刻的值是不是都等于0???
8楼2015-01-19 17:13:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

引用回帖:
8楼: Originally posted by CelestialCYJ at 2015-01-19 17:13:07
你P3~10,Q2~10在0时刻的值是不是都等于0???

是的是的。
9楼2015-01-20 13:06:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ykai_1991

金虫 (初入文坛)

引用回帖:
7楼: Originally posted by dingd at 2015-01-18 20:24:40
光定义的参数就有110个,再加上P3至 P10, Q3至Q10都没有实际数据,又要增加16个参数,也即126个参数,还是微分方程问题,计算太费时间了,一般老点的PC很吃力了,没法耗几个甚至十几个小时替你跑。

哦哦,这样哈,明白了。。那我再去修改,就不用那么多参数了。我把定义的参数都用公式表达出来,这样参数就减少了很多很多。那我弄好之后再重新发帖向您求助吧。谢谢啦!~
10楼2015-01-20 13:10:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ykai_1991 的主题更新
不应助 确定回帖应助 (注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 申博雅思求助 +3 刘嘉研究究 2024-06-24 3/150 2024-06-24 02:15 by fuqijun
[考博] 对象没有,还非常想读博,难以抉择 +10 pvrw0224 2024-06-23 30/1500 2024-06-24 00:47 by pvrw0224
[论文投稿] Pattern Recognition期刊,二审审稿邀请发出一个月了,但有一个审稿人一直未接受审稿 5+4 PLVS_VLTRA 2024-06-19 9/450 2024-06-24 00:16 by holypower
[硕博家园] 联培博士文章第一单位署名问题交流~ +4 橙成成c 2024-06-23 8/400 2024-06-23 23:43 by 风云无泪
[基金申请] 江苏省333人才工程 出校后被刷的概率大不大? +8 maxbirdzhang 2024-06-19 12/600 2024-06-23 20:39 by 查拉图斯特拉
[教师之家] 有没有今年的影响因子? +5 jurkat.1640 2024-06-22 7/350 2024-06-23 18:38 by lyfbangong
[论文投稿] OSA期刊审稿逾期 +3 Thomas_Squid 2024-06-22 3/150 2024-06-23 15:20 by wspglt
[基金申请] 工材口青年基金上会可能性 +10 今晚推荐22 2024-06-19 13/650 2024-06-23 15:11 by Pickfoot
[基金申请] 刚刚收到科研之友邮件 +27 olivermiaoer 2024-06-19 40/2000 2024-06-23 11:45 by qq632458
[硕博家园] 夏至,要不要硕博联谊 +3 我是王小帅 2024-06-21 5/250 2024-06-22 22:51 by 我是王小帅
[基金申请] 自然基金变化?? +13 今晚推荐22 2024-06-21 22/1100 2024-06-22 11:47 by Nano-Firefly
[博后之家] 在国内某高校做全职博士后2年,现在找到新的单位,出站或退站对新工作有什么影响? +10 nxplfcc 2024-06-20 10/500 2024-06-22 07:52 by 徐长安
[有机交流] 锌粉还原反应 20+4 明宝栋 2024-06-17 6/300 2024-06-22 02:29 by hwalle
[基金申请] 面上博士后 +11 jsqy 2024-06-19 13/650 2024-06-21 17:12 by sizhouyi
[有机交流] 跑板能跑开,过柱过不纯怎么办 +8 小胡在努力 2024-06-18 10/500 2024-06-21 08:22 by hptianyan
[论文投稿] 第一篇论文投稿进程记录 +4 慎独的小花卷 2024-06-20 9/450 2024-06-20 20:37 by 慎独的小花卷
[有机交流] 想要用氢化钠拔掉吲哚N上的氢取代酰氯 50+3 光敏剂 2024-06-19 4/200 2024-06-20 18:41 by HF111001
[海外博后] 墨尔本大学博后offer要不要接 +3 kyxblmm 2024-06-18 3/150 2024-06-19 22:39 by blake1111
[基金申请] 基金得中 +4 woaini0218 2024-06-18 4/200 2024-06-19 17:27 by mengzl
[基金申请] 博后基金刷到的BUG,图片来的更直观 +15 carolloo 2024-06-17 16/800 2024-06-19 09:42 by msjy
信息提示
请填处理意见