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

 Akzeptierte Antwort

Stephan
Stephan am 23 Nov. 2020
Bearbeitet: Stephan am 24 Nov. 2020

0 Stimmen

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 Kommentare

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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Hilfe-Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by