2
$\begingroup$

I would like some help in a model I'm trying to run. The physical problem is a one that is fed-batch fermentation with a simple chemical reaction: S-->X+P. My independent variable is time(t), and my dependent variables are concentrations(Cs,Cx,Cp), Volume, and flow rate(flow)(μ value is going to be constant for this problem). You can see more details from the picture uploaded.

What I would like to get after running the code is the t vs. Cs,Cx,Cp graphs. The code I input into mathematica is as follows:

rx[cx]:=0.035*cx[t]
volume[t]:=2+flow[t]*t
flow[t]:=(0.035*30*Exp[0.035*t])/(0.67*0.79)
rp[qp,cx]:=qp[ypx]*cx[t]
qp[ypx]:=ypx*0.035+0.05
ypx[cp,cx]:=cp[t]/(cx[t]-30)
rs[rx,rp,ypx,cx]:=-((rx[cx]/0.67)+rp[qp,cx]*ypx[cp,cx]+0.05*cx[t])

sol1=NDSolve[{D[cx[t]*volume[t],t]==rx[cx]*volume[t],
D[cp[t]*volume[t],t]==rp[qp,cx]*volume[t],
D[cs[t]*volume[t],t]==rs[rx,rp,ypx,cx]*volume[t]+flow[t]*0.79,
cx[0]==30.1,cs[0]==79/100,cp[0]==0},{cx,cs,cp},{t,0,30}]

During evaluation of In[192]:= NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

Out[192]= NDSolve[{(1.98375 E^(0.035 t)+0.0694313 E^(0.035 t) t) cx[t]+(2+1.98375 E^(0.035 t) t) (cx^\[Prime])[t]==0.035 (2+1.98375 E^(0.035 t) t) cx[t],(1.98375 E^(0.035 t)+0.0694313 E^(0.035 t) t) cp[t]+(2+1.98375 E^(0.035 t) t) (cp^\[Prime])[t]==(2+1.98375 E^(0.035 t) t) (0.05 +0.035 ypx) cx[t],(1.98375 E^(0.035 t)+0.0694313 E^(0.035 t) t) cs[t]+(2+1.98375 E^(0.035 t) t) (cs^\[Prime])[t]==1.56716 E^(0.035 t)+(2+1.98375 E^(0.035 t) t) (-0.102239 cx[t]-((0.05 +0.035 ypx) cp[t] cx[t])/(-30+cx[t])),cx[0]==30.1,cs[0]==79/100,cp[0]==0},{cx,cs,cp},{t,0,30}]

However, as you can see, I get a non-numerical value error. I tried to follow/simulate similar demonstration projects online, but I can't seem to figure out how to fix this issue. Any help would be appreciated.

Problem Statement

$\endgroup$
7
  • $\begingroup$ Your equation have the symbol ypx via rs[rx, rp, ypx, cx]. Perhap you want rs[rx, rp, ypx[cp,cx], cx]? $\endgroup$ Jun 26, 2021 at 17:08
  • $\begingroup$ Also, your definition of the functions in the first part is not how one generally does it. i.e., instead of volume[t] := 2 + flow[t]*t, you probably want to use t as a pattern: volume[t_] := 2 + flow[t]*t $\endgroup$ Jun 26, 2021 at 17:11
  • $\begingroup$ It is also a good idea to view the equations that you are passing to NDSolve one-by-one to see if they have a form that depends only on the unknown functions and their independent variables. $\endgroup$ Jun 26, 2021 at 17:13
  • $\begingroup$ Hey Craig, thanks for the help. I tried your first two suggestions. After turning the variable t into t_ in the function inputs, the error turned into "Function Cx has no arguments". After doing the ypx etc, now the whole system gives me "There are fewer dependent variables, {cp[t],cx[t]}, than equations, so the system is overdetermined". I have no idea how to approach to this currently. $\endgroup$ Jun 26, 2021 at 19:11
  • $\begingroup$ ypx is undefined in qp[ypx] := ypx*0.035 + 0.05, as can be seen by examine the returned expression for NDSolve`.. $\endgroup$
    – bbgodfrey
    Jun 26, 2021 at 19:25

2 Answers 2

4
$\begingroup$

First, use Clear["Global*"]` to clean up all your definitions. Second, define functions without arguments step by step

rx = 0.035*cx[t];
flow = (0.035*30*Exp[0.035*t])/(0.67*0.79); volume = 2 + flow*t;
ypx = cp[t]/(cx[t] - 30); qp = ypx*0.035 + 0.05;
rp = qp*cx[t];
rs = -((rx/0.67) + rp*ypx + 0.05*cx[t]);

Finally, use explicit Euler method to pass through point At t == 0.0034835532562812544 where step size is effectively zero; singularity or stiff system suspected. So, we use

sol1 = NDSolve[{D[cx[t]*volume, t] == rx*volume, 
   D[cp[t]*volume, t] == rp*volume, 
   D[cs[t]*volume, t] == rs*volume + flow*0.79, cx[0] == 30.1, 
   cs[0] == 79/100, cp[0] == 0}, {cx, cs, cp}, {t, 0, 30}, 
  Method -> "ExplicitEuler", StartingStepSize -> 0.001]

Visualization

{Plot[cx[t] /. sol1[[1]], {t, 0, 30}, FrameLabel -> {"t", "cx"}, 
  Frame -> True, Axes -> False, PlotRange -> All], 
 Plot[cs[t] /. sol1[[1]], {t, 0, 30}, FrameLabel -> {"t", "cs"}, 
  Frame -> True, Axes -> False, PlotRange -> All], 
 Plot[cp[t] /. sol1[[1]], {t, 0, 30}, FrameLabel -> {"t", "cp"}, 
  Frame -> True, Axes -> False, PlotRange -> All]}

Figure 1

$\endgroup$
4
  • $\begingroup$ Thanks a lot for the help! $\endgroup$ Jun 27, 2021 at 13:02
  • 2
    $\begingroup$ @IlkerKaygusuz, It is considered good manners to accept Alex's answer and upvote it. $\endgroup$ Jun 27, 2021 at 13:25
  • $\begingroup$ Since ${C_s}$ is a species concentration, it should always be positive. I suspect the balance is off in the OP and not your approach. $\endgroup$
    – Tim Laska
    Jun 28, 2021 at 14:58
  • $\begingroup$ @TimLaska Thank you very much for your remarks. $\endgroup$ Jun 28, 2021 at 16:27
2
$\begingroup$

I could not follow the derivation in the OP, but I suspect the species balance is in error. This may explain why Alex Trounev's answer has a negative species concentration.

My understanding of the system is that it looks something like this (Note my interpretation of equation 2 in the OP):

Variable batch reactor

This should in turn lead to a system of equations as shown below:

$$\begin{array}{l} \frac{{d{C_S}(t)V(t)}}{{dt}} = F{C_{S,in}} - r(t)V(t)\\ \frac{{d{C_P}(t)V(t)}}{{dt}} = r(t)V(t)\\ \frac{{d{C_X}(t)V(t)}}{{dt}} = r(t)V(t)\\ \frac{{dV(t)}}{{dt}} = F \end{array}$$

The above system of equations can easily be solved with NDSolve as shown below:

tmax = 30;
parms = {μ -> 0.035, csin -> 0.79, cs0 -> 79/100, cp0 -> 0, 
   cx0 -> 30.1, f -> (0.035*30*Exp[0.035*t])/(0.67*0.79), v0 -> 2};
With[{r = (μ cs[t] cx[t] v[t])/csin}, 
  peqns = {D[v[t] cs[t], t] == f csin - r,
    D[v[t] cp[t], t] == r,
    D[v[t] cx[t], t] ==  r,
    D[v[t], t] == f}];
pics = {cs[0] == cs0, cp[0] == cp0, cx[0] == cx0, v[0] == v0};
eqns = Join[peqns, pics] /. parms;
{csfun, cpfun, cxfun, vfun} = 
  NDSolveValue[eqns, {cs, cp, cx, v}, {t, 0, tmax}];
Plot[{csfun[t], cpfun[t], cxfun[t]}, {t, 0, tmax}, 
 PlotLegends -> Automatic]
Plot[{csfun[t], cpfun[t]}, {t, 0, tmax}, PlotRange -> All, 
 PlotLegends -> Automatic]

NDSolve solutions

Now, all the species concentrations are well behaved and consistent with the design goals of Exponential-fed-batch culture

"If the feed rate of the growth-limiting substrate is increased in proportion to the exponential growth rate of the cells, it is possible to maintain the cells' specific growth rate for a long time while keeping the substrate concentration in the culture liquid at a constant level."

$\endgroup$
1
  • $\begingroup$ It is nice picture and code (+1). I understand, that original post has some issues :) $\endgroup$ Jun 28, 2021 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.