24小时热门版块排行榜

南方科技大学公共卫生及应急管理学院2025级博士研究生招生报考通知

返回列表

当前主题已经存档。

jesschen

铜虫 (小有名气)

应助: 0 (幼儿园)
金币: 136.5
帖子: 56
在线: 29.9小时
虫号: 580497
注册: 2008-07-11
专业: 化学工程

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

由于时间紧迫，来不及去从头学习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 ]

回复此楼