24小时热门版块排行榜    

查看: 1107  |  回复: 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个字符以上)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[论文投稿] 编辑是什么意思 15+3 s090604054 2024-06-15 3/150 2024-06-16 10:29 by bobvan
[有机交流] 车间生产,真空度很高,温度很高,但减压蒸馏速度很慢。 10+12 召唤鬼泣lL 2024-06-13 36/1800 2024-06-16 09:20 by ddc805
[基金申请] 博后面上今天有bug可以看到是否资助? +20 lyfbangong 2024-06-12 31/1550 2024-06-15 21:18 by since—2010
[基金申请] 关于博后基金的bug问题 +6 lxr1991 2024-06-14 9/450 2024-06-15 21:17 by since—2010
[考博] 上海交大招收材料化学方向科研助理/“申请考核”博士(请勿回复帖子或站内投条) +3 灵梦and紫 2024-06-12 4/200 2024-06-15 20:58 by 1822836277
[论文投稿] 审稿人含糊拒稿,还需要回复吗?如何回复? 20+3 BruceChum 2024-06-15 17/850 2024-06-15 20:19 by arthas_007
[教师之家] 我们学院常年位居 各学院 倒数第二。专业撤销的话,在编者有什么补偿? +13 河西夜郎 2024-06-09 14/700 2024-06-15 19:44 by LittleBush
[基金申请] 2024国社科通讯评审 +9 qsd10086 2024-06-13 14/700 2024-06-15 15:51 by thesuna
[基金申请] 为什么我的博后基金还在流动站审核中?不会是学院给我卡了吧? +14 王凯12 2024-06-13 26/1300 2024-06-15 15:22 by 好人与坏人
[论文投稿] 投稿时忘记修改一作 +7 gll123456 2024-06-13 11/550 2024-06-15 11:49 by gll123456
[论文投稿] 求机械类四区sci推荐 5+3 迷茫小旷 2024-06-14 4/200 2024-06-15 11:25 by bobvan
[基金申请] 有没有机械的前辈分享一下评上海优都是什么成果啊 +7 wulala800 2024-06-10 7/350 2024-06-15 09:33 by 晓目崇
[食品] 食品博士导师 +6 小李醒yy 2024-06-11 9/450 2024-06-14 23:37 by 小李醒yy
[论文投稿] 审稿问题:为什么荧光激发波长和紫外吸收波长差的大? 10+4 sdawege 2024-06-14 8/400 2024-06-14 22:39 by 东北读书人
[基金申请] 工材E口JQ有消息了吗 +4 babyduck 2024-06-11 4/200 2024-06-14 17:23 by firepick
[基金申请] 75批博后基金 +10 kyukitu 2024-06-13 13/650 2024-06-14 10:31 by kyukitu
[考博] 博导选择 +3 bing85977 2024-06-12 3/150 2024-06-13 15:34 by 我是邱尧
[有机交流] 原料反应完了,怎么知道是产物还是中间体 +6 小胡在努力 2024-06-11 8/400 2024-06-13 13:33 by 091602
[硕博家园] 求助 +6 LYWwrz 2024-06-09 9/450 2024-06-11 13:12 by powerhours
[论文投稿] water research状态咨询 5+3 Flyyawa 2024-06-10 6/300 2024-06-11 09:45 by bobvan
信息提示
请填处理意见