Simulation of differential equations with multiple variables. What is wrong with my code? Monod-like kinetic, Biotechnology, Differential equations, ODE15s

My problem goes like this:

I have to simulate the following Equations

I would do this with a function summing up all the equatiobs and ODE15s.

The code of my function looks like this:

function [ f ] = dXdt( ~, x )
%% Variables%%
Cx = x(1);
cp = x(2);
mu = x(3);
c_s = x(4);
nu = x(5);
%% Constants%%
R=8.3144621; %kJ/(K*kmol)
Kcm=0.03;   % 1/h
Kd0=1.836*10^(-3);% 1/h
Ks=60; % g/L
Ksp=50; % g/L
Kss= 500; % g/L
Kssp=620; % g/L
Pm= 70; % g/L
PIm= 70; % g/L
Cp0= 0; % g/L
Cs0= 220; % g/L
Cx0=0.2; % g/L
E=3.461*10^(4);% kJ/kmol
Ed=1.777*10^(5); % kJ/kmol
Ep=3.496*10^(4);
Y_x=0.41;
Y_p=0.42;
mu_m=0.414; %1/h
nu_m=3.974; %1/h
nu_0=2.511833179;
mu_0=0.241719745;
T=308.15; % K

%% Equations%%

dCxdt=mu*Cx-Kd0*Cx;
dCpdt=Cx*nu;
mu=mu_m*(c_s/(Ks+c_s+(c_s^2)/Kss))*(1-cp/Pm);
dc_sdt=(-1)*((1/Y_x)*dCxdt+(1/Y_p)*dCpdt+Kcm*Cx);
nu=nu_m*(c_s/(Ksp+c_s+(c_s^2)/Kssp))*(1-cp/PIm);

f=[dCxdt; dCpdt; mu; dc_sdt; nu];
end

When executing this function, using

   [t,x]=ode15s(@dXdt,[0,72],[Cx0 Cp0 mu_0 Cs0 nu_0]);
I get a time error:
Warning: Failure at t=4.971802e+00.  Unable to meet integration tolerances without reducing
the step size below the smallest value allowed (1.421085e-14) at time t. 
> In ode15s at 731

But I can't spot the bug in my code, or the transcription/translation of the equations.

Or is this even the proper way to solve this?

If it works, the plots of t,cs ... should look like those on the right.

Thanks in advance for your help !

Answer 1

[Torsten]

The equations for mu and nu are algebraic equations, but at the moment you solve

dmu/dt=mu_m*(c_s/(Ks+c_s+(c_s^2)/Kss))*(1-cp/Pm);
dnu/dt=nu_m*(c_s/(Ksp+c_s+(c_s^2)/Kssp))*(1-cp/PIm);

Best wishes

Torsten.

[Biotechnewbie]

Ah! Thanks a million!

Finally.

[Anna Gardarsdottir]

Hi

Can I ask, 6 years later..

What did you actually do to solve this problem? As I understand it, you were attempting to solve mu and nu as differential equations. But I can't see what needs to be changed in the code to fix this.

[darova]

HEre is what you should change

% dCxdt=mu*Cx-Kd0*Cx;
% dCpdt=Cx*nu;
% mu=mu_m*(c_s/(Ks+c_s+(c_s^2)/Kss))*(1-cp/Pm);
% dc_sdt=(-1)*((1/Y_x)*dCxdt+(1/Y_p)*dCpdt+Kcm*Cx);
% nu=nu_m*(c_s/(Ksp+c_s+(c_s^2)/Kssp))*(1-cp/PIm);
% f=[dCxdt; dCpdt; mu; dc_sdt; nu];
mu=mu_m*(c_s/(Ks+c_s+(c_s^2)/Kss))*(1-cp/Pm);
nu=nu_m*(c_s/(Ksp+c_s+(c_s^2)/Kssp))*(1-cp/PIm);
dCxdt=mu*Cx-Kd0*Cx;
dCpdt=Cx*nu;
dc_sdt=(-1)*((1/Y_x)*dCxdt+(1/Y_p)*dCpdt+Kcm*Cx);
f=[dCxdt; dCpdt; dc_sdt];

Remember about the ode45 call int main code

[Anna Gardarsdottir]

Ahhh right.. thank you!!