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.
ParametricPlot[{kb[t], cb[t]}, {t, 0, 60}, <options>]
. $\endgroup$ – march Jul 27 '15 at 3:41