How to determine steady state of X given code

2 次查看（过去 30 天）

显示更早的评论

Justin Lee 2021-8-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/897332-how-to-determine-steady-state-of-x-given-code

评论： Justin Lee 2021-8-16

在 MATLAB Online 中打开

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