ODE Function time output

1 次查看(过去 30 天)
Dear all,
I wrote the code below now, I have initial concentration 6mg/L but i need to calculate when it reaches (0.05*6) mg/L. Could you help me please? I need the t value when concentration is equal to (0.05*6) mg/L
Z0=6;
tspan = [0 24];
[tZ,Z] = ode45(@ConcDCE,tspan,Z0);
function dZdt=ConcDCE(t,Z)
k1=1.26;
k2=0.74;
k3=0.22;
X0=1;
Y0=4;
Z0=6;
dZdt=k2*(Y0*exp(-k2*t)+((X0*k1)/(k2-k1)*(exp(-k1*t)-exp(-k2*t))))-k3*(Z0*exp(-k3*t)+((Y0*k2)/(k3-k2)*(exp(-k2*t)-exp(-k3*t)))+X0*k1*k2*((exp(-k1*t))/((k2-k1)*(k3-k1))-(exp(-k2*t))/((k2-k1)*(k3-k2))+(exp(-k3*t))/((k3-k1)*(k3-k2))));
end

采纳的回答

Stephan
Stephan 2020-11-23
编辑:Stephan 2020-11-24
Use events:
Z0 = 6;
Zt=0.05*6;
tspan = [0 24];
opts = odeset('Events',@(tZ,Z)EventsFcn(tZ,Z,Zt));
[tZ,Z,tZe,Ze,iZe] = ode45(@ConcDCE,tspan,Z0,opts);
plot(tZ,Z)
hold on
scatter(tZe,Ze,'or')
function dZdt=ConcDCE(t,~)
k1=1.26;
k2=0.74;
k3=0.22;
X0=1;
Y0=4;
Z0=6;
dZdt=k2*(Y0*exp(-k2*t)+((X0*k1)/(k2-k1)*(exp(-k1*t)-exp(-k2*t))))-k3*(Z0*exp(-k3*t)+((Y0*k2)/(k3-k2)*(exp(-k2*t)-exp(-k3*t)))+X0*k1*k2*((exp(-k1*t))/((k2-k1)*(k3-k1))-(exp(-k2*t))/((k2-k1)*(k3-k2))+(exp(-k3*t))/((k3-k1)*(k3-k2))));
end
function [Conc,isterminal,direction] = EventsFcn(~,Z,Zt)
Conc = Z - Zt; % The value that we want to be zero
isterminal = 0; % Halt integration
direction = 0; % The zero can be approached from either direction
end
  4 个评论
Stephan
Stephan 2020-11-24
I edited my code in my answer. For some reason I missed Z0 as initial value - accidentally I used Zt as initial value. This is corrected now. Sorry for this.
Carey n'eville
Carey n'eville 2020-11-24
The last one is working properly thank you so much!!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by