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

 Réponse acceptée

Stephan
Stephan le 23 Nov 2020
Modifié(e) : Stephan le 24 Nov 2020

0 votes

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 commentaires

Carey n'eville
Carey n'eville le 23 Nov 2020
Modifié(e) : Carey n'eville le 23 Nov 2020
%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 le 23 Nov 2020
Modifié(e) : Stephan le 23 Nov 2020
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 le 24 Nov 2020
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 le 24 Nov 2020
The last one is working properly thank you so much!!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Centre d'aide et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by