Solving Multiple ODE's with NDSolve

Questions
,
1 
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?

2

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,

3

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.

4 
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.

5

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}
6

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.

7 
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.
8

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?

9 
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)

10

Cool. Thanks for the help.