我是matlab初学者,刚刚接触。
有问题如下:
我已经用ode45将微分方程的图画出来了,但是如果我想要具体某个t时间下,各物质的浓度,如何做到,是不是有什么函数,我看了有限的例子,好像都是把图画出了,求极值之类的,没有具体求某个t下y的具体值的。
很菜的问题,希望高手能释疑。谢谢!
附例:参考书上的,比如我想求出t=3050s时,C(1)、C(2)、C(3)、C(4)的具体值,如何实现?
function BatchReactor
% 间歇反应器中的连串-平行复杂反应系统
%
% Author: Huang Huajiang
% Copyright 2003 UNILAB Research Center,
% East China University of Science and Technology, Shanghai, PRC
% $Revision: 1.0 $ $Date: 2002/07/16 $
clear all
clc
T = 224.6 + 273.15; % Reactor temperature, Kelvin
R = 8.31434; % Gas constant, kJ/kmol K
% Arrhenius constant, 1/s
k0 = [5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8];
% Activation energy, kJ/kmol
Ea = [124670 150386 77954 111528];
% 初始浓度C0(i), kmol/m^3
C0 = [1 0 0 0 0];
tspan = [0 1e4];
[t,C] = ode45(@MassEquations, tspan, C0,[],k0,Ea,R,T)
% 绘图
plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--');
xlabel('Time (s)');
ylabel('Concentration (kmol/m^3)');
legend('A','B','C','D')
CBmax = max(C(:,2)); % CBmax: the maximum concentration of B, kmol/m^3
yBmax = CBmax/C0(1) % yBmax: the maximum yield of B
index = find(C(:,2)==CBmax);
t_opt = t(index) % t_opt: the optimum batch time, s
% ------------------------------------------------------------------
function dCdt = MassEquations(t,C,k0,Ea,R,T)
% Reaction rate constants, 1/s
k = k0.*exp(-Ea/(R*T));
k(5) = 2.16667E-04;
% Reaction rates, kmoles/m3 s
rA = -(k(1)+k(2))*C(1);
rB = k(1)*C(1)-k(3)*C(2);
rC = k(2)*C(1)-k(4)*C(3);
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances
dCdt = [rA; rB; rC; rD; rE];