0 个投票

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

0 个投票

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 个评论
显示 2更早的评论 隐藏 2更早的评论

Carey n'eville
Carey n'eville 2020-11-23
编辑：Carey n'eville 2020-11-23
%I couldn't explain myself clearly sorry for this.
% I actually search for labeling Z=0.3 on graph (in this problem Z),
% I need to put a dot to graph
Stephan
Stephan 2020-11-23
编辑：Stephan 2020-11-23
I think there is a mark at the position of Z=0.3 in the graph - maybe i misunderstood. If you only want the second point where Z=0.3 use
scatter(tZe(2),Ze(2),'or')
The t-value is:
tZe(2)
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 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by