2
$\begingroup$

Whatever I have tried, I could not achieve to plot numerical values coming from my system of differential equation. It is a optimal control model with two state variables. My differential system is as follows ;

   dek = k'[t] == k[t]^α (r[t] s[t])^β - δ k[t] - c[t] - g k[t];
   des = s'[t] == -r[t] s[t] - (-rr) s[t];
   dec = c'[t] == c[t]/σ (α k[t]^(α - 1) (r[t] s[t])^β - δ - ρ);
   der = r'[t] ==  r[t]/(β -  1) (α  c[t]/k[t] - δ (1 - α)) + r[t]^2; 

I use NDSolve to solve the system (I have calculated initial values in previous stages which are not relevant to my question) ;

soldif[t_] = with[{σ = 1.8, α = 0.4, β = 0.1, ρ= 0.02, δ = 0.05, g = -0.0029411764705882357,  rr = 0.017647058823529415}, 
NDSolve[{dec, dek, der, des, c[0] == 0.130567, k[0] == 0.167055, r[0] == 0.0157895, 
s[0] == 1}, {c, k, r, s}, {t, 0, 60}], 
Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}]

I define ;

  cb[t_] := Evaluate[soldif[t_][[2, 1, 1, 2]]]
  kb[t_] := Evaluate[soldif[t_][[2, 1, 2, 2]]]
  rb[t_] := Evaluate[soldif[t_][[2, 1, 3, 2]]]
  sb[t_] := Evaluate[soldif[t_][[2, 1, 4, 2]]]

Finally, I use the following plot command in order to draw a phase diagram on a plane (c,k) but it does not work ;

  tra1 = Plot[{kb[t], cb[t]}, {t, 0, 60}, AxesLabel -> {k, c}, PlotStyle -> Thickness[0.1]]

I have no idea why it does not work, any hints or suggestions ? Thanks so much.

$\endgroup$
  • 1
    $\begingroup$ If you want a phase plot, you should be doing ParametricPlot[{kb[t], cb[t]}, {t, 0, 60}, <options>]. $\endgroup$ – march Jul 27 '15 at 3:41
9
$\begingroup$

First, you used with when you should have used With and the assignments to cb,... won't work this way. Let me give you an easier alternative that uses NDSolveValue

With[{σ=1.8,α=0.4,β=0.1,ρ=0.02,
      δ=0.05,g=-0.0029411764705882357,rr=0.017647058823529415},
  dek=k'[t]==k[t]^α (r[t] s[t])^β-δ k[t]-c[t]-g k[t];
  des=s'[t]==-r[t] s[t]-(-rr) s[t];
  dec=c'[t]==c[t]/σ (α k[t]^(α-1) (r[t] s[t])^β-δ-ρ);
  der=r'[t]==r[t]/(β-1) (α c[t]/k[t]-δ (1-α))+r[t]^2;
]

{cb, kb, rb, sb} = NDSolveValue[{dec, dek, der, des, c[0] == 0.130567, 
    k[0] == 0.167055, r[0] == 0.0157895, s[0] == 1}, 
    {c, k, r, s}, {t, 0, 60}, 
   Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}]


Plot[{cb[t], rb[t]}, {t, 0, 60}]

Mathematica graphics

Some explanation

OK, the issue of your wrong with needs a bit of explanation, because when you simply replace it with With, it won't work after the next start of Mathematica (or Kernel restart).

What happens when you use With is, that it replaces all occurrences of the defined variables with their value. The important part is that it won't work if the variable is not explicitly in the code. Example:

ClearAll[a];
expr = a + a;
With[{a = 2},
 {expr, a + a}
]
(* {2 a, 4} *)

Note that only the explicit a+a is replaced but not expr which evaluates afterwards to 2a.

Now, why did it work when you used with but it doesn't when you use With? Exactly for this reason. Assume you wanted to calculate the above expression and you mistyped With with boing. Let's see what happens and why it happens:

boing[{a = 2},
 expr
 ]
(* boing[{2}, 4] *)

You see your function boing is still there (because doesn't have a definition) but on the way of evaluating the arguments, you accidentally assigned the global variable a=2 which is then used when expr is evaluated. Therefore, it is kind of the same as doing

a = 1;
expr

The conclusion to all this is that you need to use With where you really want to replace the parameters alpha, etc.. and this is in the assignment of dek, des, ... Therefore, I have put the With around this to make the replacement work.

$\endgroup$
  • 1
    $\begingroup$ Despite the fact that there is a Plot command at the end of the original post, I'm pretty sure the OP wants a phase plot, so probably Plot should be ParametricPlot. $\endgroup$ – march Jul 27 '15 at 3:43
  • $\begingroup$ @halirutan thanks for this clear and helpful answer. I think it could be easy to implement a phase diagram with ParametricPlot as march has mentionned above. $\endgroup$ – user1905 Jul 27 '15 at 10:13

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.