Solving a system of kinetic equations in the steady state

Questions
, , , ,
1 
Clear["Global'*"];
nut = 46; khi = 16; ro2 = 3400; kmrs = 4; kisu = 1.5; kmp = .5; kres \
= 25;
kccc1 = 38; k23 = 16; kvp = .04; kc = 21; nc = 5;
km2 = 15; nm2 = 5; kv1 = 17; nv1 = 6; kv2 = 200; nv2 = 4; k32 = 18; \
n32 = 9; a = .5;
odec = -(a c[t]) - (
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1) - (
   kmrs c[t])/((fm[t]/km2)^nm2 + 1) + (khi nut)/((c[t]/kc)^nc + 1);
odefm = -(a fm[t]) + (kmrs c[t])/((fm[t]/km2)^nm2 + 1) - kisu fm[t] - 
   kmp fm[t] o2[t];
odefs = kisu fm[t] - a fs[t];
odemp = kmp fm[t] o2[t] - a mp[t];
odeo2 = -(kmp fm[t] o2[t]) - kres fs[t] o2[t] + 
   ro2/((o2[t]/245)^9 + 1);
odef2 = -(a f2[t]) - k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) + (
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1);
odef3 = -(a f3[t]) + k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) - 
   kvp f3[t];
odevp = kvp f3[t] - a vp[t];
vars = {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]};
solution = 
  NDSolve[{Derivative[1][c][t] == odec, Derivative[1][f2][t] == odef2,
     Derivative[1][f3][t] == odef3, Derivative[1][fm][t] == odefm, 
    Derivative[1][fs][t] == odefs, Derivative[1][mp][t] == odemp, 
    Derivative[1][o2][t] == odeo2, Derivative[1][vp][t] == odevp, 
    c[0] == 30, f2[0] == 100, f3[0] == 100, fm[0] == 50, fs[0] == 50, 
    mp[0] == 0, o2[0] == 100, vp[0] == 100}, vars, {t, 0, 31}];
vars
vars /. solution /. t -> 31

I am trying to solve these equations for the steady state conditions (i.e. c’=odec=0 and so on). When I try to add in these conditions, be it together or separately, I get that 0 is protected or that the system has more equations than dependent variables. I might also be using NDSolve wrong in this case and I might be better off using DSolve or something. Any hints or suggestions would be much appreciated.

2

If you’re trying to find the steady state solution, why not set all the time derivatives equal to zero and use NSolve to solve the resulting equations?

3

Often “0 is protected” is because you are trying to assign a value to something that has already been made into zero.

When that happens I click Evaluation->Quit the Kernal->Local->Yes I really mean to do this
and then try evaluating the notebook fresh.

Attached is a notebook with a few changes that should not make any difference for your error message.
It evaluates just fine. This again makes me think this was a previously assigned value.

4 
Clear["Global'*"];
nut = 46; khi = 16; ro2 = 3400; kmrs = 4; kisu = 1.5; kmp = .5; kres \
= 25;
kccc1 = 38; k23 = 16; kvp = .04; kc = 21; nc = 5;
km2 = 15; nm2 = 5; kv1 = 17; nv1 = 6; kv2 = 200; nv2 = 4; k32 = 18; \
n32 = 9; a = .5;
odec = -(a c[t]) - (
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1) - (
   kmrs c[t])/((fm[t]/km2)^nm2 + 1) + (khi nut)/((c[t]/kc)^nc + 1);
odefm = -(a fm[t]) + (kmrs c[t])/((fm[t]/km2)^nm2 + 1) - kisu fm[t] - 
   kmp fm[t] o2[t];
odefs = kisu fm[t] - a fs[t];
odemp = kmp fm[t] o2[t] - a mp[t];
odeo2 = -(kmp fm[t] o2[t]) - kres fs[t] o2[t] + 
   ro2/((o2[t]/245)^9 + 1);
odef2 = -(a f2[t]) - k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) + (
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1);
odef3 = -(a f3[t]) + k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) - 
   kvp f3[t];
odevp = kvp f3[t] - a vp[t];
vars = {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]};
solution = 
 NSolve[{Derivative[1][c][t] == odec == 0, 
   Derivative[1][f2][t] == odef2 == 0, 
   Derivative[1][f3][t] == odef3 == 0, 
   Derivative[1][fm][t] == odefm == 0, 
   Derivative[1][fs][t] == odefs == 0, 
   Derivative[1][mp][t] == odemp == 0, 
   Derivative[1][o2][t] == odeo2 == 0, 
   Derivative[1][vp][t] == odevp == 0}, vars]

My attempt at NSolve, I tried to run it and quit the kernel after 5 minutes.

5 
Clear["Global'*"];
nut = 46; khi = 16; ro2 = 3400; kmrs = 4; kisu = 1.5; kmp = .5; kres \
= 25;
kccc1 = 38; k23 = 16; kvp = .04; kc = 21; nc = 5;
km2 = 15; nm2 = 5; kv1 = 17; nv1 = 6; kv2 = 200; nv2 = 4; k32 = 18; \
n32 = 9; a = .5;
odec = -a c[t] - 
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1))/((f3[t]/kv2)^nv2 + 1) - 
   kmrs c[t]/((fm[t]/km2)^nm2 + 1) + khi nut/((c[t]/kc)^nc + 1);
odefm = -a fm[t] + kmrs c[t]/((fm[t]/km2)^nm2 + 1) - kisu fm[t] - 
   kmp fm[t] o2[t];
odefs = kisu fm[t] - a fs[t];
odemp = kmp fm[t] o2[t] - a mp[t];
odeo2 = -kmp fm[t] o2[t] - kres fs[t] o2[t] + ro2/((o2[t]/245)^9 + 1);
odef2 = -a f2[t] - k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) + 
   kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1))/((f3[t]/kv2)^nv2 + 1);
odef3 = -a f3[t] + k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) - kvp f3[t];
odevp = kvp f3[t] - a vp[t];
vars = {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]};
solution = 
 NDSolve[{c'[t] == odec = 0, f2'[t] == odef2 = 0, f3'[t] == odef3 = 0,
    fm'[t] == odefm = 0, fs'[t] == odefs = 0, mp'[t] == odemp = 0, 
   o2'[t] == odeo2 = 0, vp'[t] == odevp = 0, c[0] == 30, f2[0] == 100,
    f3[0] == 100, fm[0] == 50, fs[0] == 50, mp[0] == 0, o2[0] == 100, 
   vp[0] == 100}, vars, {t, 0, 31}]
vars /. solution /. t -> 31

Tried quitting the kernel, still getting the same result with the protected value. By using the == assignment instead of =, I get the overdetermination error.

6

Suggestion:

Go back to your

NDSolve[{c'[t] == odec, f2'[t] == odef2,...

form.

Take the InterpolatingFunction solution for each of your c[ t ] and f2[ t ] and … functions from that.

Look at http://reference.wolfram.com/mathematica/tutorial/ApproximateFunctionsAndInterpolation.html down where it starts “If you differentiate an approximate function…” and study that information to see if you can determine how to then solve for where the derivative of your InterpolatingFunction will equal zero.

I think that is what you are really trying to accomplish.

7

The equations and conditions in NDSolve are always expressed using == rather than =.

In your first argument of NDSolve you have a bunch of things that say xxxxx==yyyyy=0. I am not sure that those =0’s are there for. If you remove them as in

solution = 
 NDSolve[{c'[t] == odec, f2'[t] == odef2, f3'[t] == odef3, 
   fm'[t] == odefm, fs'[t] == odefs, mp'[t] == odemp, o2'[t] == odeo2,
    vp'[t] == odevp, c[0] == 30, f2[0] == 100, f3[0] == 100, 
   fm[0] == 50, fs[0] == 50, mp[0] == 0, o2[0] == 100, vp[0] == 100}, 
  vars, {t, 0, 31}]

then your NDSolve evaluates promplty.

8

I went back to my previous format, and am trying to follow the documentation you sent me, but I guess I need to assign the interpolating functions from the list that NDSolve provides. I tried flattening the list and assigning like that but it doesn’t appear to work. Any suggestions?

9

I went back to my previous format, and am trying to follow the documentation you sent me, but I guess I need to assign the interpolating functions from the list that NDSolve provides. I tried flattening the list and assigning like that but it doesn’t appear to work. Any suggestions? (repost because my name did not appear)

10

Yes I can do that. What I am trying to see is where the derivatives of the equations are zero so I can show the steady state conditions. By plotting out the functions (Plot[c/.solution,{t,0,31}]) and so on, I can see that the functions are relatively flat past around t=24. What I wanted to verify is that the derivatives past this point are indeed zero.

11

Math folks, and Mathematica make a big distinction between really close to zero and exactly zero. Even for something exponentially decreasing towards zero it is not going to be exactly zero.

In[1]:= ...CodeFromYourPreviousPostSnipped...
vars = {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]};
solution = NDSolve[{c'[t] == odec, f2'[t] == odef2, f3'[t] == odef3, fm'[t] == odefm, 
    fs'[t] == odefs, mp'[t] == odemp, o2'[t] == odeo2, vp'[t] == odevp, 
    c[0] == 30, f2[0] == 100, f3[0] == 100, fm[0] == 50, fs[0] == 50, 
    mp[0] == 0, o2[0] == 100, vp[0] == 100}, vars, {t, 0, 31}];
allf = vars /. solution[[1]]; (* functions *)
allfprime = Map[D[#, t] &, allf]; (* derivatives *)

c[t] is first in the list, look at a plot

In[17]:= Plot[allf[[1]], {t, 0, 31}]

Out[17]= ...PlotSnipped...

Look at a plot of c'[t]

In[18]:= Plot[allfprime[[1]], {t, 0, 31}]

Out[18]= ...PlotSnipped...

Beyond t == 20 it looks close to zero, but how close?

In[19]:= Table[allfprime[[1]], {t, 0, 31}]

Out[19]= {-947.944, 0.726472, 0.0672147, -0.00443547, -0.00923631, -0.00695448, 
-0.00470561, -0.00311015, -0.00204019, -0.00133365, -0.00087007, 
-0.000566914, -0.000369095, -0.000240173, -0.000156233, -0.000101607, 
-0.0000660714, -0.0000429592, -0.000027931, -0.0000181579, -0.0000118052, 
-7.67559*10^-6, -4.9897*10^-6, -3.24209*10^-6, -2.10926*10^-6, 
-1.37057*10^-6, -9.01694*10^-7, -8.27389*10^-7, -3.74807*10^-7, 
-2.41653*10^-7, -1.55533*10^-7, -1.03173*10^-7}

So if you are going to solve that for exactly zero you have a problem, unless you get the solution near t==3 that you don' t want.

Let's compromise and, after looking at that plot of c'[t], find where c'[t]== -.001, which looks like it should be around t= 10

In[20]:= FindRoot[allfprime[[1]] == -.001, {t, 10}]

Out[20]= {t -> 9.67451}

So you might be able to see how to use that to "solve" for t for other values of your derivatives.

None of your derivatives are exactly zero, even at t==31, but most are darned small and getting smaller.

In[21]:= allfprime /. t -> 31

Out[21]= {-1.03173*10^-7, 1.68236*10^-8, -5.48649*10^-7, -4.21631*10^-7,
 -9.11946*10^-6, 0.0000527061, 6.29529*10^-7, -7.2189*10^-6}
12

Thanks Bill, another few tricks to add to my growing notes.