1
|
18305 Views
|
9 Replies
|
2 Total Likes
View groups...
Share
GROUPS:

# Solving Multiple ODE's with NDSolve

Posted 11 years ago
 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?
9 Replies
Sort By:
Posted 11 years ago
 Cool. Thanks for the help.
Posted 11 years ago
 In[15]:= vars /. solution /. t -> 31Out[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 11 years 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 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 11 years 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 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 11 years ago
 Clear["Global'*"]  nut = 46 khi = 100 ro2 = 50 kmrs = 180 kisu = .05 kmp = 1 kres = 1kccc1 = 30k23 = 16kvp = .04kc = 21nc = 5km1 = 24nm1 = 3km2 = 50nm2 = 4kv1 = 29nv1 = 9kv2 = 200nv2 = 4k32 = 18n32 = 9a = .5odec = 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 11 years 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 11 years 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,