ODE Function time output

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

 Risposta accettata

Stephan
Stephan il 23 Nov 2020
Modificato: Stephan il 24 Nov 2020
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 Commenti

Carey n'eville
Carey n'eville il 23 Nov 2020
Modificato: Carey n'eville il 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 il 23 Nov 2020
Modificato: Stephan il 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 il 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 il 24 Nov 2020
The last one is working properly thank you so much!!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Programming in Centro assistenza e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by