# Solving a system of kinetic equations in the steady state

Posted 9 years ago
10950 Views
|
11 Replies
|
7 Total Likes
|
 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}];varsvars /. solution /. t -> 31I 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.
11 Replies
Sort By:
Posted 9 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 plotIn[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= 10In[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 -> 31Out[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 9 years ago
 Thanks Bill, another few tricks to add to my growing notes.
Posted 9 years ago
 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 9 years ago
 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 9 years ago
 Suggestion:Go back to yourNDSolve[{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.
Anonymous User
Anonymous User
Posted 9 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 9 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? (repost because my name did not appear)
Posted 9 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 thisand 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 9 years ago
 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 -> 31Tried quitting the kernel, still getting the same result with the protected value. By using the == assignment instead of =, I get the overdetermination error.
Posted 9 years ago
 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 9 years ago
 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.