Message Boards Message Boards

0
|
12091 Views
|
11 Replies
|
7 Total Likes
View groups...
Share
Share this post:

Solving a system of kinetic equations in the steady state

 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.
POSTED BY: Josh Wofford
11 Replies
Sort By:
Thanks Bill, another few tricks to add to my growing notes.
POSTED BY: Josh Wofford
Bill Simpson
Bill Simpson
Posted 11 years ago
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}
POSTED BY: Bill Simpson
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.
POSTED BY: Josh Wofford
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)
POSTED BY: Josh Wofford
Updating Name
Updating Name
Posted 11 years ago
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?
POSTED BY: Updating Name
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.
POSTED BY: David Reiss
Bill Simpson
Bill Simpson
Posted 11 years ago
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.
POSTED BY: Bill Simpson
 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.
POSTED BY: Josh Wofford
 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.
POSTED BY: Josh Wofford
Bill Simpson
Bill Simpson
Posted 11 years ago
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.
Attachments:
POSTED BY: Bill Simpson
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?
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract