3
$\begingroup$

I have the following differential equation

r = 1.4;
b = 0.36;
k = 0.4;
\[Epsilon] = 0.1;
a = 10;
sol[t_] = 
 NDSolve[{x'[t] == -2 x [t] + r (k + 1) y[t] + 2( k + 1) p[t],
   y'[t] == -(2 k + 1)/(k + 1) x[t] + r k y[t] + (2 k + 1) p[t],
   z'[t] == b p[t],
   \[Epsilon]  p'[t] == 1/(k + 1) x[t] - z[t] - (1 + a) p[t] + a,
   x[0] == 0.1, y[0] == 0.5, z[0] == 0.2, p[0] == 1.3}, {x, y, z, 
   p}, {t, 0, 500}]

ParametricPlot[{x[t], y[t]} /. sol[t_], {t, 200, 400}, 
 PlotRange -> Full]

With these parameter values a limit cycle oscillation should appear. However I am not getting any limit cycle. Please help me to rectify the problem.

$\endgroup$
2
  • $\begingroup$ If I change your code to sol=NDSolve[...{x[t],y[t],z[t],p}...][[1]]; ParametricPlot3D[ {x[t],y[t],z[t]}/.sol, {t,0,500}, PlotRange->Full] then I do get a visible plot. That is as close to your code I can get and not have what I think are errors. If I restrict that to a 2D plot with ParametricPlot[ {x[t],y[t]}/.sol, {t,0,500}, PlotRange->Full] then I get a very different plot. That doesn't address your limit cycle issue, but it should get you to the point where you can actually see a plot. $\endgroup$
    – Bill
    May 9 '18 at 3:04
  • $\begingroup$ This looks like a linear system, but limit cycles are a nonlinear phenomenon (physics.stackexchange.com/questions/235687/…). Why do you think there is a limit cycle to be found? $\endgroup$
    – Chris K
    May 9 '18 at 14:12
4
$\begingroup$

With $r=3.25$, $b=2.36$, $k=0.14$, $ϵ=1000$, and $a=0.01$ your "limit cycle" appears.

r = 3.25; b = 2.36; k = 0.14; ϵ=1000; a = 0.01;

sol[t_] = NDSolve[{x'[t] == -2 x[t] + r (k + 1) y[t] + 2 (k + 1) p[t], 
        y'[t] == -(2 k + 1)/(k + 1) x[t] + r k y[t] + (2 k + 1) p[t], 
        z'[t] == b p[t], ϵ p'[t] == 1/(k + 1) x[t] - z[t] - (1 + a) p[t] + a, 
        x[0] == 0.1, y[0] == 0.5, z[0] == 0.2, p[0] == 1.3}, {x, y, z, p},
        {t, 0, 2000}];

First[ParametricPlot3D[{x[t], y[t], z[t]} /. sol[t_], 
        {t, 0, 2000}, PlotRange -> All, PlotStyle -> {Black, Thickness[0.002]}, 
        LabelStyle -> Directive[Black, Small], PlotPoints -> 1000, 
        BoxRatios -> {1, 1, 1}, AspectRatio -> 1, PlotTheme -> "Detailed"]]

Plot[x[t] /. sol[t_], {t, 0, 500}, PlotStyle -> {Blue, Thickness[0.002]}, 
        AxesStyle -> Directive[Black, Small, Arrowheads[0.03]], 
        LabelStyle -> Directive[Black, Small]]

Phase portrait and time series enter image description here

Calling it the limit cycle is an abuse of language. It is a center. Only in piecewise linear systems is it possible to see limit cycles.

$\endgroup$
2
$\begingroup$

If you rationalize all parameters, you can solve the dgls analytically.

rat = Rationalize[{x'[t] == -2 x[t] + r (k + 1) y[t] + 2 (k + 1) p[t],
        y'[t] == -(2 k + 1)/(k + 1) x[t] + r k y[t] + (2 k + 1) p[t], 
        z'[t] == b p[t], \[Epsilon] p'[t] == 
          1/(k + 1) x[t] - z[t] - (1 + a) p[t] + a, x[0] == 0.1, 
        y[0] == 0.5, z[0] == 0.2, p[0] == 1.3}, 0];

sol = DSolve[rat, {x, y, z, p}, t];

You need very high workingPrecision for higher t.

plxyzp = Plot[
           Evaluate[{x[t], y[t], z[t], p[t]} /. First@sol], {t, 0, 20}, 
           WorkingPrecision -> 1000, PlotRange -> All, 
           PlotStyle -> {Red, Blue, Green, Magenta}]

enter image description here

$\endgroup$
2
  • $\begingroup$ What does dgls mean? $\endgroup$
    – Chris K
    May 9 '18 at 14:10
  • 1
    $\begingroup$ Have a guess! (just: differential equations. My mixed german-english abreviation) $\endgroup$
    – Akku14
    May 10 '18 at 15:03

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.