0
$\begingroup$

I try to plot a phase diagram with ParametricPlot but everytime Mathematica does nothing.

I start by ;

sol = Solve[{α A k^(α - 1) - ρ - δ == 0,  A k^α - δ k - c == 0}, {c, k}] // FullSimplify


{csol, ksol} = {sol[[1, 1, 2]], sol[[1, 2, 2]]}


Clear[dec, dek, α, A, σ, ρ, δ];

dec = c'[t] == -(c[t]/σ) (A α k[t] ^(α - 1) - δ - ρ);

dek = k'[t] == -(A k[t]^α - c[t] - δ k[t]);

Clear[α, σ, ρ, δ];
α = 0.3; σ = 1; ρ = 0.04; δ = 0.02; A = 1;

Clear[nsba, cb, kb, cf, kf];
nsba = NDSolve[  {dec, dek, c[0] == csol (1 - 10^-5), k[0] == 
                                    ksol (1 - 10^-5)}, {c[t], k[t]}, {t, 0, 255}]

cb[t_] := Evaluate[nsba[[1, 1, 2]]]
kb[t_] := Evaluate[nsba[[1, 2, 2]]]

Clear[tende];
tende = 164;
cf[t_] := cb[tende - t]
kf[t_] := kb[tende - t]

After all, I try this last command ;

tra = ParametricPlot[{kf[t], cf[t]}, {t, 0, tende},
                     AxesLabel -> {k, c}, PlotStyle -> Thickness[0.01]];

Mathematica stays actionless. What could be the eventual problem ?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ is the first Solve[ ] working for you? What is returning as solutions? $\endgroup$ – Dr. belisarius Jul 21 '15 at 22:32
  • 2
    $\begingroup$ The semicolon after your plot command suppresses output. $\endgroup$ – Simon Rochester Jul 21 '15 at 22:33
  • $\begingroup$ @belisarius yes the command Solve[ ] works fine. It gives me the right solution for k and c that I have also solved by hand. $\endgroup$ – optimal control Jul 21 '15 at 22:38
  • 1
    $\begingroup$ Please post the solutions here. Thanks! $\endgroup$ – Dr. belisarius Jul 21 '15 at 22:41
  • $\begingroup$ Optimal, unfortunately very little works in your code. Even trying to understand what you want to do by reading the code is pretty hard. Maybe you could describe the problem you are trying to solve in your question as well. As @belisarius hinted in his comment, the first Solve doesn't seem to work for me as written (it throws multiple errors). $\endgroup$ – MarcoB Jul 21 '15 at 23:20
2
$\begingroup$

There seem to be problems with your use of both Solve and NDSolve. However, ignoring the warnings Mathematica issues and plunging ahead, I get

sol = 
  Solve[{α A k^(α - 1) - ρ - δ == 0, A k^α - δ k - c == 0}, {c, k}] // FullSimplify;
{csol, ksol} = {sol[[1, 1, 2]], sol[[1, 2, 2]]};
dec = c'[t] == -(c[t]/σ) (A α k[t]^(α - 1) - δ - ρ);
dek = k'[t] == -(A k[t]^α - c[t] - δ k[t]);
α = 0.3; σ = 1; ρ = 0.04; δ = 0.02; A = 1;
nsba = 
  NDSolve[{dec, dek, c[0] == csol (1 - 10^-5), k[0] == ksol (1 - 10^-5)}, 
    {c[t], k[t]}, {t, 0, 255}];

At this point I arrive at the point in the code that is affecting your plot. I fix that by defining

cb = Head @ nsba[[1, 1, 2]];
kb = Head @ nsba[[1, 2, 2]];

Now

With[{tende = 164.},
  With[{cf = cb[tende - #] &, kf = kb[tende - #] &},
    ParametricPlot[{kf[t], cf[t]}, {t, 0, tende},
      AxesLabel -> {k, c},
      AspectRatio -> 1,
      PlotStyle -> Thickness[0.01]]]]

produces

plot

I have no idea whether or not this is a meaningful plot, but at least it shows you what your code is generating.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks so much, it is very helpful :) I think the problem is possible to be solved analytically. $\endgroup$ – optimal control Jul 22 '15 at 14:32

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.