How to determine steady state of X given code
2 views (last 30 days)
Show older comments
I have the code below, but I need to determine the Steady state of total SA reached, and need a graph
Question:
- The total SA concentration is approximated by adding SA1 and SA2, and division by 140 dL (for a normal person). When is a steady state of total SA reached? Support by a graph (extend the time if needed).
STARTTIME = 0;
STOPTIME = 24;
% initial values
SA1 = 0;
SA2 = 0;
ASA1 = 0;
ASA2 = 0;
y0 = [SA1; SA2; ASA1; ASA2];
% solve ode
[t,y] = ode45(@odefun,[STARTTIME STOPTIME],y0);
% Plot results
yyaxis left
plot(t,y(:,[1,3]))
ylabel("SA1,ASA1")
yyaxis right
plot(t,y(:,[2,4]))
ylabel("SA2,ASA2")
legend(["SA1","ASA1","SA2","ASA2"])
function ddt = odefun(t,y)
SA1 = y(1);
SA2 = y(2);
ASA1 = y(3);
ASA2 = y(4);
% Constants
Vnf = 1;
Cnf = SA1/Vnf;
Qnf = 1;
Vf = 1;
Cf = SA2/Vf;
Qf = 1;
Vb = 1;
Cb = ASA1/Vb;
Pf = 20;
Cfv = Cf/Pf;
Vl = 1;
Cl = ASA2/Vl;
Ql = 1;
vmax = 1;
Km = 1;
Pl = 2;
Clv = Cl/Pl;
Pb = 18;
% Benzene conc in inhaled air
Ci = 0.32;
% alveolar ventilation rate
Qp = 5.74;
% Flows
Jnf = Qnf*(Cb-Cnf);
Jf = Qf*(Cb-Cfv);
Jl = Ql*(Cb - Clv);
Jmetab = vmax*Cl/(Km + Cl);
% Replace squarepulse with " IF t>6h Jresp=0 "
if t<=6
Jresp = Qp*(Ci - (Cb/Pb));
else
Jresp = 0;
end
% Reservoirs Blood, Fat, NonFat and Liver
ddt(1,1) = Jnf;
ddt(2,1) = Jf;
ddt(3,1) = -Jnf - Jf - Jl + Jresp;
ddt(4,1) = Jl - Jmetab;
end
2 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!