Message Boards Message Boards

1
|
18831 Views
|
9 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Solving Multiple ODE's with NDSolve

  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
9 Replies
Sort By:
Cool. Thanks for the help.
POSTED BY: Josh Wofford
Bill Simpson
Bill Simpson
Posted 11 years 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
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
Bill Simpson
Bill Simpson
Posted 11 years 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
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
Bill Simpson
Bill Simpson
Posted 11 years 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
 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
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
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
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