24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 398  |  回复: 4
当前主题已经存档。
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

jesschen

铜虫 (小有名气)

[交流] 【求助】请高手帮忙修改一个动力学的程序,多谢

由于时间紧迫,来不及去从头学习matlab了。所以希望各位高手帮帮忙编一下程序。需要金币的话可以相送。
我的动力学方程式子是:2A——>B+C,是个可逆反应,正反应速率常数为k1,逆反应为k2。现在知道A,B,C各个时间的浓度,不过C的似乎有些问题,不过不要紧的。目前没有做到平衡,所以不能由平衡时的数值来算常数。
我想用的是采用定步长的龙格库塔方数值积分和Powell发来估算模型的参数k1,和k2。
希望高手帮忙编一下程序,我自己仿照例子编了半天,不是这里有问题就是那里有问题。我自己编的如下:
数据是随便写的
function KineticsEst1_int

% 动力学方程为rA=dCA/dt=-k1*CA^2+k2*CB*CC
clear all;
clc
global CAm CBm CCm
t=[0 10 30 60 90 150 210 270 330];
CAm=[10 9 8 7 6 5 4 3 2];
CBm=[0 0.5 1.5 2 2.5 3 3.5 4 4.5]
CCm = [有问题,不列了];

% 非线性拟合
beta0=[0.0000053 10];
tspan = [0 10 30 60 90 150 210 270 330];
CA0 = 9.672;
CB0 = 0;
CC0 = 0;
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@OptObjFunc,beta0,[0 0],[],[],tspan,[CA0 CB0 CC0])
ci = nlparci(beta,resid,jacobian)

% 拟合效果图(实验与拟合的比较)
[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1)  tspan(end)],[CA0 CB0 CC0],[],beta);
plot(tspan,CAm,'bo',t4plot,CA4plot,'k-')
legend('Exp','Model')
xlabel('时间t, min')
ylabel('浓度C_A, mol/L')

%残差关于拟合值的残差图
[t Y]=ode45(@KineticsEqs,tspan,CA0,[],beta);
Figure, plot(CAc,resid, '*')
xlabel('浓度拟合值(mol/L)'), ylabel('残差R(mol/L)'), refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n'),
fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tk2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2))
% ------------------------------------------------------------------
function f = OptObjFunc(beta,tspan,CA0,CAm)
global CAm CBm CCm
[t Y] = ode45(@KineticsEqs,tspan,CA0,[],beta);
if length(t)==length(tspan)
    f = Y - [CAm CBm CCm];
else
    f=100000;
end
% ------------------------------------------------------------------
function dy = KineticsEqs(t,Y,beta)    这个地方一直看不懂例子为什么这么编
CA=Y(1);
CB=Y(2);
CC=Y(3);
dCAdt = -beta(1)*CA^2 +beta(2)*CB*CC;            % k1= beta(1), k2= beta(2)
dCBdt = beta(1)*CA*CB-beta(2)*CC;
dCCdt = beta(2)*CC;
dy=[dCAdt;dCBdt;dCCdt];
多谢各位。

[ Last edited by jesschen on 2009-6-19 at 18:19 ]
回复此楼

» 猜你喜欢

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

ljling

木虫 (小有名气)


小木虫(金币+0.5):给个红包,谢谢回帖交流
我也想学格,一起想想办法。。。。
3楼2009-06-08 16:57:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 5 个回答

youol

金虫 (小有名气)

sunxiao(金币+0,VIP+0):此话怎讲?请解释,以免误会,呵呵 6-10 04:53
kuhailangyu(金币+0,VIP+0):是楼主有点为难人吧? 6-10 16:23
版主有点为男人呢!
2010,该丰收的一年!
2楼2009-06-01 18:13:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

youol

金虫 (小有名气)


小木虫(金币+0.5):给个红包,谢谢回帖交流
是有点难的意思!打字打错了!sorry!
2010,该丰收的一年!
4楼2009-06-10 10:27:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

youol

金虫 (小有名气)


sunxiao(金币+1,VIP+0):谢谢,原来是这个意思,呵呵,那应该是楼主吧,而不应该是版主,嘿嘿 6-11 04:58
呵呵  就是这个意思!
2010,该丰收的一年!
5楼2009-06-10 22:57:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 283分求调剂 +6 试试看呗 2026-04-04 6/300 2026-04-05 08:42 by barlinike
[考研] 材料调剂 +5 dxy调剂 2026-04-04 5/250 2026-04-05 08:34 by qlm5820
[考研] 材料工程310专硕调剂 +12 捞捞我…. 2026-04-04 13/650 2026-04-05 08:30 by qlm5820
[考研] 复试调剂 +9 呼呼?~+123456 2026-04-05 9/450 2026-04-05 08:26 by 无际的草原
[考研] 324求调剂 +9 想上学求调 2026-04-03 9/450 2026-04-04 23:57 by 果冻大王
[考研] 求调剂 +7 xzghyuj 2026-04-04 7/350 2026-04-04 22:25 by oooqiao
[考研] 349求调剂 +11 zwjjjjjj 2026-03-31 11/550 2026-04-04 19:52 by 蓝云思雨
[考研] 化学308分调剂 +23 你好明天你好 2026-03-30 24/1200 2026-04-04 18:29 by macy2011
[考研] 一志愿0817化学工程与技术,求调剂 +24 我不是只因 2026-04-02 28/1400 2026-04-04 15:15 by dongzh2009
[考研] 材料295 +13 小英11 2026-04-03 14/700 2026-04-04 09:02 by 来看流星雨10
[考研] 0856调剂 +8 曲听筠 2026-03-30 8/400 2026-04-04 08:46 by tianyyysss
[考研] 考研调剂 +3 15615482637 2026-04-03 3/150 2026-04-03 22:50 by ms629
[考研] 求调剂 +4 压力??大 2026-04-03 4/200 2026-04-03 21:36 by 啵啵啵0119
[考研] 343求调剂085601 +6 要努力学习x 2026-03-29 7/350 2026-04-03 19:49 by 百灵童888
[考研] 材料科学与工程339求调剂 +12 hyz0119 2026-03-31 13/650 2026-04-03 18:33 by ls刘帅
[考研] 求调剂不挑专业 +3 xrh030412 2026-04-01 3/150 2026-04-03 14:40 by 氮气气气
[考研] 316求调剂 +14 舟自梗 2026-04-01 18/900 2026-04-03 10:28 by linyelide
[考研] 材料工程322分 +8 哈哈哈吼吼吼哈 2026-04-01 8/400 2026-04-02 11:53 by 3041
[考研] 261求B区调剂 +5 明仔· 2026-04-01 7/350 2026-04-02 11:17 by 邹尉尉
[考研] 环境工程 085701,267求调剂 +15 minht 2026-03-29 16/800 2026-04-01 10:13 by li_sujuan99
信息提示
请填处理意见