| 查看: 474 | 回复: 4 | |||
| 当前主题已经存档。 | |||
[交流]
【求助】请高手帮忙修改一个动力学的程序，多谢
|
|||
|
由于时间紧迫，来不及去从头学习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 ] |
回复此楼» 猜你喜欢
FileCode能看出啥？
已经有15人回复
-
filecode
已经有8人回复
-
确定了，国自然21号放榜
已经有7人回复
-
帮忙看看fileCode
已经有5人回复
-
综述论文作为代表作会不会影响评审专家的印象分？
已经有11人回复
-
关于代码变化问题，想知道的进来
已经有24人回复
-
为什么网上很多人说本周 12号出结果
已经有3人回复
-
静等基金结果
已经有16人回复
-
我的国基提前知道中了，可是同事的操作让我实在接受不了，怎么会有这样的人
已经有11人回复
-
关于Filecode分析方法
已经有3人回复
2楼2009-06-01 18:13:40
3楼2009-06-08 16:57:44
4楼2009-06-10 10:27:41
5楼2009-06-10 22:57:49
100