0
$\begingroup$

I have a system of coupled ODEs which I am solving using Mathematica. I am solving these using NDSolve:

solved=NDSolve[{Cs'[t] == DL*Stuff, 
  Cx1'[t] == DL*Stuff, 
  Cx2'[t] == DL*Stuff, Cs[0] == 10, Cx1[0] == 25, 
  Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Where DL is a constant. Generally following this I will plot them as follows:

Plot[Evaluate[{Cs[t], Cx1[t], Cx2[t]} /. solved], {t, 0, 1000}, 
 PlotRange -> All]

However, what I would like to do is make a 3D plot which shows how the solution to the system changes for different values of the constant DL. How can I do this?

Constraint: The equations mentioned above cannot be solved analytically, or decoupled or whatever.

Full code after attempting changes:

(First, I will Define Constants)

Subscript[Y, sx1] = 1/.14;
U1max = .5;
Subscript[KM, 1] = Subscript[KM, 2] = 10; 
Subscript[Y, x1x2] = 2 ;
U2max = .11;

(Next, I will define the prior algebraic expressions)

Clear[Cs]

U1 = (U1max*Cs[t])/(Subscript[KM, 1] + Cs[t]);

U2 = (U2max*Cs[t])/(Subscript[KM, 2] + Cs[t]);

(Now I put in the remaining algebraic expressions)

rgx1 = U1*Cx1[t];

rgx2 = U2*Cx2[t];

(Now we can enter the coupled ODEs)

solved[DL_] = 
 NDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Subscript[Y, sx1]*rgx1, 
   Cx1'[t] == -DL*Cx1[t] + rgx1 - Subscript[Y, x1x2]*rgx2, 
   Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
   Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Plot3D[Evaluate[Cx2[t] /. solved], {t, 0, 1}, {DL, .01, .1}, 
 PlotRange -> All]

I get two errors:

(1) NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

(2) ReplaceAll::reps: {solved} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

$\endgroup$
1
$\begingroup$

Here you go:

solved[DL_] := 
 NDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Subscript[Y, sx1]*rgx1, 
  Cx1'[t] == -DL*Cx1[t] + rgx1 - Subscript[Y, x1x2]*rgx2, 
  Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
  Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Plot3D[Cx2[tt] /. solved[DL]], {tt, 0, 1000}, {DL, 0.01, 0.1}, 
  PlotRange -> All]

enter image description here

Some important points or it won't work:

(1) the value tt for time has to have a different name to the value t for time inside the NDSolve

(2) the SetDelayed in the definition of solved is important

(3) it's horribly inefficient

But it works nicely.

Comment:

Usually the above approach isn't very useful; e.g. you would find it more useful to plot the final state as a function of DL.

Plot[Evaluate[{Cs[1000], Cx1[1000], Cx2[1000]} /. solved[DL]], {DL, .01, .1}]
$\endgroup$
  • $\begingroup$ I have tried to include the full code. My apologies if it is difficult to read. After attempting to take these steps I have received two errors: (1) NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >> (2) ReplaceAll::reps: {solved} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >> $\endgroup$ – Nick L. Mar 15 '15 at 8:55
  • $\begingroup$ I suggest you edit your code so it can easily be copied and pasted into Mathematica. You need to use SetDelayed and solved[DL] as shown in my edits. $\endgroup$ – djp Mar 15 '15 at 10:24
  • $\begingroup$ Thank you very much. The code works very nicely. I also appreciate the tips which you gave me. I would up-vote you but apparently my reputation is too low. $\endgroup$ – Nick L. Mar 15 '15 at 18:52
0
$\begingroup$

Rewriting the code:

Ysx1 = 1/.14;
U1max = .5;
KM1 = KM2 = 10;
Yx1x2 = 2;
U2max = .11;
Clear[Cs]
U1 = (U1max*Cs[t])/(KM1 + Cs[t]);
U2 = (U2max*Cs[t])/(KM2 + Cs[t]);
rgx1 = U1*Cx1[t];
rgx2 = U2*Cx2[t];
sol = ParametricNDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Ysx1*rgx1, 
   Cx1'[t] == -DL*Cx1[t] + rgx1 - Yx1x2*rgx2, 
   Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
   Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}, {DL}]

Assuming the desired 3D plot relates to {DL, t, Cx2}:>

f[DL_, t_] := (Cx2[DL] /. sol)[t]
Plot3D[Evaluate[f[x, y]], {y, 0, 1000}, {x, 0.01, 0.1}, 
PlotRange -> All, AxesLabel -> {"t", "DL", "Cx2"}, Mesh -> False]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for the code. I am very appreciative that you took the time to do this. This will be very helpful to me. $\endgroup$ – Nick L. Mar 15 '15 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.