Group Abstract Group Abstract

Message Boards Message Boards

Solving Multiple ODE's with NDSolve

GROUPS:
  In[1]:= Clear["Global'*"]
  
  In[2]:= nut = 46
  khi = 100
  ro2 = 50
  kmrs = 180
  kisu = .05
  kmp = 1
  kres = 1
 kccc1 = 30
 k23 = 16
 kvp = .04
 kc = 21
 nc = 5
 km1 = 24
 nm1 = 3
 km2 = 50
 nm2 = 4
 kv1 = 29
 nv1 = 9
 kv2 = 200
 nv2 = 4
 k32 = 18
 n32 = 9
 a = .5
 
 Out[2]= 46
 
 Out[3]= 100
 
 Out[4]= 50
 
 Out[5]= 180
 
 Out[6]= 0.05
 
 Out[7]= 1
 
 Out[8]= 1
 
 Out[9]= 30
 
 Out[10]= 16
 
 Out[11]= 0.04
 
 Out[12]= 21
 
 Out[13]= 5
 
 Out[14]= 24
 
 Out[15]= 3
 
 Out[16]= 50
 
 Out[17]= 4
 
 Out[18]= 29
 
 Out[19]= 9
 
 Out[20]= 200
 
 Out[21]= 4
 
 Out[22]= 18
 
 Out[23]= 9
 
 Out[24]= 0.5
 
 In[25]:= odec =
  khi*(1/(1 + (c[t]/kc)^(nc)))*nut -
   kmrs*(1/(1 + (c[t]/km1)^(nm1)))*(1/(1 + (fs[t]/km2)^(nm2)))*c[t] -
   kccc1 (1 - (1/(1 + (c[t]/kv1)^(nv1))))*(1/(1 + (f3[t]/kv2)^(nv2)))*
    c[t] - a*c[t]
 odefm = kmrs*(1/(1 + (c[t]/km1)^(nm1)))*(1/(1 + (fs[t]/km2)^(nm2)))*
    c[t] - kisu*fm[t] - a*fm[t]
 odefs = kisu*fm[t] - kmp*fs[t]*o2[t] - a*fs[t]
 odemp = kmp*fs[t]*o2[t] - a*mp[t]
 odeo2 = ro2 - kres*fs[t]*o2[t] - a*o2[t]
 odef2 = kccc1*(1 - (1/(1 + (c[t]/kv1)^(nv1))))*(1/(1 + (f3[t]/
           kv2)^(nv2)))*c[t] -
   k23*(1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - a*f2[t]
 odef3 = k23*(1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - kvp*f3[t] -
   a*f3[t]
 odevp = kvp*f3[t] - a*vp[t]
 
 Out[25]= -0.5 c[t] + 4600/(1 + c[t]^5/4084101) - (
  30 c[t] (1 - 1/(1 + c[t]^9/14507145975869)))/(
  1 + f3[t]^4/1600000000) - (
  180 c[t])/((1 + c[t]^3/13824) (1 + fs[t]^4/6250000))
 
 Out[26]= -0.55 fm[t] + (
  180 c[t])/((1 + c[t]^3/13824) (1 + fs[t]^4/6250000))
 
 Out[27]= 0.05 fm[t] - 0.5 fs[t] - fs[t] o2[t]
 
Out[28]= -0.5 mp[t] + fs[t] o2[t]

Out[29]= 50 - 0.5 o2[t] - fs[t] o2[t]

Out[30]= -0.5 f2[t] - 16 (1 - 1/(1 + c[t]^9/198359290368)) f2[t] + (
30 c[t] (1 - 1/(1 + c[t]^9/14507145975869)))/(1 + f3[t]^4/1600000000)

Out[31]= 16 (1 - 1/(1 + c[t]^9/198359290368)) f2[t] - 0.54 f3[t]

Out[32]= 0.04 f3[t] - 0.5 vp[t]

In[36]:= vars =
List[c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]]

Out[36]= {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]}

In[38]:= 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}]

During evaluation of In[38]:= NDSolve::deqn: Equation or list of equations expected instead of -0.5 c[t]+4600/(1+c[t]^5/4084101)-(30 c[t] (1-1/(1+Power[<<2>>]/14507145975869)))/(1+f3[t]^4/1600000000)-(180 c[t])/((1+c[t]^3/13824) (1+fs[t]^4/6250000)) in the first argument {-0.5 c[t]+4600/(1+c[<<1>>]^5/4084101)-(30 c[t] (1-1/(1+Times[<<2>>])))/(1+f3[<<1>>]^4/1600000000)-(180 c[t])/((1+c[<<1>>]^3/13824) (1+fs[<<1>>]^4/6250000)),-0.5 f2[t]-16 (1-1/(1+Times[<<2>>])) f2[t]+(30 c[t] (1-1/(1+Times[<<2>>])))/(1+f3[<<1>>]^4/1600000000),16 (1-1/(1+Times[<<2>>])) f2[t]-0.54 f3[t],-0.55 fm[t]+(180 c[t])/((1+c[<<1>>]^3/13824) (1+fs[<<1>>]^4/6250000)),0.05 fm[t]-0.5 fs[t]-fs[t] o2[t],-0.5 mp[t]+fs[t] o2[t],50-0.5 o2[t]-fs[t] o2[t],0.04 f3[t]-0.5 vp[t],c[0]==30,f2[0]==100,f3[0]==100,fm[0]==50,fs[0]==50,mp[0]==0,100,vp[0]==100}. >>

Out[38]= NDSolve[{-0.5 c[t] + 4600/(1 + c[t]^5/4084101) - (
   30 c[t] (1 - 1/(1 + c[t]^9/14507145975869)))/(
   1 + f3[t]^4/1600000000) - (
   180 c[t])/((1 + c[t]^3/13824) (1 + fs[t]^4/6250000)), -0.5 f2[t] -
   16 (1 - 1/(1 + c[t]^9/198359290368)) f2[t] + (
   30 c[t] (1 - 1/(1 + c[t]^9/14507145975869)))/(
   1 + f3[t]^4/1600000000),
  16 (1 - 1/(1 + c[t]^9/198359290368)) f2[t] -
   0.54 f3[t], -0.55 fm[t] + (
   180 c[t])/((1 + c[t]^3/13824) (1 + fs[t]^4/6250000)),
  0.05 fm[t] - 0.5 fs[t] - fs[t] o2[t], -0.5 mp[t] + fs[t] o2[t],
  50 - 0.5 o2[t] - fs[t] o2[t], 0.04 f3[t] - 0.5 vp[t], c[0] == 30,
  f2[0] == 100, f3[0] == 100, fm[0] == 50, fs[0] == 50, mp[0] == 0,
  100, vp[0] == 100}, {c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t],
   vp[t]}, {t, 0, 31}]
Having a little bit of trouble working on a system of ODE's. I think I have the syntax right but it keeps giving me an error so obviously not. Any pointers?
POSTED BY: Josh Wofford
Answer
7 months ago
I am only a student, very new to mathematica, and no expert in anyways. I expect that someone else can definitely can give a better insight than I can.

As far as I understand, NDSolve do not take such complicated equations as eqns.

I hope this helps to give a little clue to the real issue.

Thank you,
POSTED BY: Seokin Yeh
Answer
7 months ago
Josh, could you brush up your code a little bit - why are you giving the outputs? People can run the code and see for themselves. Now it is just cluttering and makes it difficult to copy-paste to notebook. You should post the code that other people can just copy to notebook to check.
POSTED BY: Sam Carrettie
Answer
7 months ago
 Clear["Global'*"]
 
 nut = 46
 khi = 100
 ro2 = 50
 kmrs = 180
 kisu = .05
 kmp = 1
 kres = 1
kccc1 = 30
k23 = 16
kvp = .04
kc = 21
nc = 5
km1 = 24
nm1 = 3
km2 = 50
nm2 = 4
kv1 = 29
nv1 = 9
kv2 = 200
nv2 = 4
k32 = 18
n32 = 9
a = .5

odec = khi*(1/(1 + (c[t]/kc)^(nc)))*nut -
  kmrs*(1/(1 + (c[t]/km1)^(nm1)))*(1/(1 + (fs[t]/km2)^(nm2)))*c[t] -
  kccc1 (1 - (1/(1 + (c[t]/kv1)^(nv1))))*(1/(1 + (f3[t]/kv2)^(nv2)))*
   c[t] - a*c[t]
odefm = kmrs*(1/(1 + (c[t]/km1)^(nm1)))*(1/(1 + (fs[t]/km2)^(nm2)))*
   c[t] - kisu*fm[t] - a*fm[t]
odefs = kisu*fm[t] - kmp*fs[t]*o2[t] - a*fs[t]
odemp = kmp*fs[t]*o2[t] - a*mp[t]
odeo2 = ro2 - kres*fs[t]*o2[t] - a*o2[t]
odef2 = kccc1*(1 - (1/(1 + (c[t]/kv1)^(nv1))))*(1/(1 + (f3[t]/
          kv2)^(nv2)))*c[t] -
  k23*(1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - a*f2[t]
odef3 = k23*(1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - kvp*f3[t] -
  a*f3[t]
odevp = kvp*f3[t] - a*vp[t]

vars = List[c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], vp[t]]

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}]
Sorry about that, rather new to using the forums. Here is the code in input mode only.
POSTED BY: Josh Wofford
Answer
7 months ago
Looks like some = versus == problems.
 In[1]:= Clear["Global'*"];
 nut = 46; khi = 100; ro2 = 50; kmrs = 180; kisu = .05; kmp = 1; kres = 1;
 kccc1 = 30; k23 = 16; kvp = .04; kc = 21; nc = 5; km1 = 24; nm1 = 3;
 km2 = 50; nm2 = 4; kv1 = 29; nv1 = 9; kv2 = 200; nv2 = 4; k32 = 18; n32 = 9; a = .5;
 odec = khi*(1/(1 + (c[t]/kc)^(nc)))*nut - kmrs*(1/(1 + (c[t]/km1)^(nm1)))*
    (1/(1 + (fs[t]/km2)^(nm2)))*c[t] - kccc1 (1 - (1/(1 + (c[t]/kv1)^(nv1))))*
    (1/(1 + (f3[t]/kv2)^(nv2)))*c[t] - a*c[t];
 odefm = kmrs*(1/(1 + (c[t]/km1)^(nm1)))*(1/(1 + (fs[t]/km2)^(nm2)))*c[t] - kisu*fm[t] - a*fm[t];
 odefs = kisu*fm[t] - kmp*fs[t]*o2[t] - a*fs[t];
odemp = kmp*fs[t]*o2[t] - a*mp[t];
odeo2 = ro2 - kres*fs[t]*o2[t] - a*o2[t];
odef2 = kccc1*(1 - (1/(1 + (c[t]/kv1)^(nv1))))*(1/(1 + (f3[t]/kv2)^(nv2)))*c[t] - k23*
   (1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - a*f2[t];
odef3 = k23*(1 - (1/(1 + (c[t]/k32)^(n32))))*f2[t] - kvp*f3[t] - a*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]};
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}]

Out[12]= {{c[t] -> InterpolatingFunction[{{0., 31.}},<>][t], etc, etc, etc}
POSTED BY: Bill Simpson
Answer
7 months ago
Thanks a lot Bill. If you could, would you show me how to take those interpolating functions to plots? When I did this for a single ODE, Mathematica was nice enough to give me a button where it automatically does it for me, but not this time.
POSTED BY: Josh Wofford
Answer
7 months ago
In[1]:= Clear["Global'*"];
   (*blah blah blah*)
   solution = NDSolve[{c'[t] == odec, (*blah blah blah*)}, vars, {t, 0, 31}];
   Plot[c[t] /. solution, {t, 0, 31}]

Out[15]= ...PlotSnipped... And use that last line as an example for all your other plots too.
POSTED BY: Bill Simpson
Answer
7 months ago
Thanks. One (hopefully last) question for now. If I wanted to make an array that output the values of the eight ODE's at t=31, how would I code that?
POSTED BY: Josh Wofford
Answer
7 months ago
In[15]:= vars /. solution /. t -> 31

Out[15]= {{27.2861, 7.58359, 219.506, 1530.97, 54.0141, 99.0828, 0.917194, 17.5605}}

I suggest creating a small separate Mathematica notebook that you leave on your desktop. Each time you discover a trick, like how to plot an InterpolatingFunction, you take a minute or two to create a very small well written simple example and paste it into that desktop notebook with a clear descriptive title in a Section cell above the trick. Then you will be able to refer back to the trick the next time you need it. Gradually that notebook will grow with the tools and techniques that you will use over and over again. (Hint: That was how I was able to remember the answer to your question)
POSTED BY: Bill Simpson
Answer
7 months ago
Cool. Thanks for the help.
POSTED BY: Josh Wofford
Answer
7 months ago