Euler showing exponential growth instead of steady state: Monod growth equation
2 次查看(过去 30 天)
显示 更早的评论
I have a euler model that should be able to reach a steady state, however I am having some difficulty with one of the functions r_x, continually growing with time, instea of getting to the point of a steady state that it should be, and I'm pretty sure it due to how I'm referrring to other functions within the line.
Full Euler loop is:
for i=1:nt
time(i+1) = time(i) + dt;
% Functions to calculate parameter values for each timestep
r_x(i) = Mu_max * f_I * (P(i)/(K_s + P(i))) * X(i);
r_p(i) = r_x(i) * 1/Y_xp;
P_out(i) = P(i) * Overflow / V; %units = mg
X_out(i) = X(i) * Overflow / V; %units = mg
% Rates of change for P, Algae, and Volume
dP(i) = P_in - P_out(i) - r_p(i);
dX(i) = r_x(i) - X_out(i);
% P, Algae and Volume equations
P(i+1) = P(i) + dt * dP(i); % mass P in mg
X(i+1) = X(i) + dt * dX(i); % mass X in mg
% Calculation for Algae concentration at each timestep
X_con(i)= X(i)*1000/V; % Algae concentration vector in mg/L
end
I've narrowed the issue down to the r_x(i) line, in particular how it requires X(i). I have a feeling matlab doesnt like how I have loops of variables that feed off each other with each timetep, since the value for X(i) is dependent on dX(i), which uses r_x(i).... so you can see my issue.
This is my first time attempting this kind of project, not sure what the best way to address this is?
回答(1 个)
Bjorn Gustavsson
2022-10-13
Instead of using this explicit integratin of your coupled (?) odes with Euler's method, you should learn to implement the ODEs in a function and then use any of the built-in tools for ode-integration, for example ode45.
function dPdtdXdt = your_PX_ode(t,PX,other_parameters)
P = PX(1);
X = PX(2);
Mu_max = other_parameters(1);
f_I = other_parameters(2);
K_s = other_parameters(3);
Y_xp = other_parameters(4);
Overflow = other_parameters(5);
r_x = Mu_max * f_I * (P/(K_s + P)) * X;
r_p(i) = r_x * 1/Y_xp;
P_out = P * Overflow / V; %units = mg
X_out = X * Overflow / V; %units = mg
% Rates of change for P, Algae, and Volume
dP = P_in - P_out - r_p;
dX = r_x - X_out;
dPdtdXdt = [dP;
dX];
end
That ode-function you can use with the odeNN-suite of ode-integrating functions. Perhaps something like this:
Mu_m = 12; % You know what values these
f_I = 2^.5*pi;% parameters should have
K_s = exp(3); % the rest of us dont.
Y_xp = exp(-2);
O_f = 12;
other_p = [Mu_m,f_I,K_s,Y_xp,O_f];
PX0 = [0.12,1.32];
t_out = linspace(0,100,201);
[t_o,PX] = ode45(@(t,xp) your_PX_ode(t,xp,other_p),t_out,px0);
That way you separate the problems of the ODE-integration into two separate parts, the definition of the ODE-equations, and then the integration. That is a way preferable structure, since it becomes far easier to isolate where you have your problems.
HTH
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Dialog Boxes 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!