Euler showing exponential growth instead of steady state: Monod growth equation

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?

[Torsten]

We need the full code to test it.

Answer 1

[Bjorn Gustavsson]

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