How to get Limit cycle oscillation from the following differential equation?

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.1, y == 0.5, z == 0.2, p == 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.

• If I change your code to sol=NDSolve[...{x[t],y[t],z[t],p}...][]; 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.
– Bill
May 9 '18 at 3:04
• 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? May 9 '18 at 14:12

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.1, y == 0.5, z == 0.2, p == 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]},
LabelStyle -> Directive[Black, Small]]

Phase portrait and time series 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.

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.1,
y == 0.5, z == 0.2, p == 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}] • What does dgls mean? May 9 '18 at 14:10
• Have a guess! (just: differential equations. My mixed german-english abreviation) May 10 '18 at 15:03