# Limit cycles and phase portraits of a two-dimensional polynomial quadratic systems using Mathematica

Cross-posted: https://community.wolfram.com/groups/-/m/t/2589679

I have elementary questions about the construction of phase portraits and limit cycles using Mathematica.

It should be remarked that I have also searched and read previous posts from this Mathematica Forum. In addition, I have read Mathematica's documentation. However, these did not help me with the questions to be presented here. Owing to the fact that phase portrait and phase spaces are quite interesting topics on nonlinear dynamics, I think that these questions, and their possible answers, may be of great interest to those that subscribe to the Mathematica forum.

With that said, suppose that we have the following nonlinear dynamical system

$$\displaystyle \frac{dx}{dt}=x^{2}+xy+y,$$

$$\displaystyle \frac{dy}{dt}=Ax^{2}+Bxy+Cy^{2}+Dx+Ey$$,

in which $$A$$, $$B$$, $$C$$, $$D$$, and $$E$$ are real constant parameters. It is well-known that the mathematical model above has four isolated periodic trajectories, that is to say, four limit cycles when we take representative values for the parameters, notably $$A=-10$$, $$B=2.2$$, $$C=0.7$$, $$D=-72.7778$$, and $$E=0.0015$$.

I have tried to construct the phase portrait of the model above by employing the StreamPlot command below:

curve = StreamPlot[{x^2 + y + x y, (-72.7778* x)-10 x^2 + (0.0015* y) + (2.2* x y) + (0.7*y^2)}, {x, -5.0,20.0}, {y, -60.0, 30.0},FrameLabel -> (Style[#, Black, 20, Bold] & /@ {"x","y"}),TicksStyle -> Directive[Bold, Black]];


Despite being correct, the code does not allow one to clearly see isolated periodic trajectories. Based on the above, I ask:

How may I construct the phase portrait for those differential equations so that the limit cycles become visible and highlighted? That is to say, I want a construct a phase portrait like the one presented below:

• Try: LineIntegralConvolutionPlot[{{x^2 + y + x y, (-72.7778*x) - 10 x^2 + (0.0015*y) + (2.2*x y) + (0.7*y^2)}, {"noise", 500, 500}}, {x, -1, 1}, {y, -3, 3}, ColorFunction -> "Rainbow", LineIntegralConvolutionScale -> 3, Frame -> False]
– Moo
Aug 1 at 14:34
• You can also use StreamPlot and then use NDSolve with various initial conditions and then overlay the two.
– Moo
Aug 1 at 14:43
• "It is well-known that the mathematical model above has four isolated periodic trajectories" -- reference? Aug 1 at 20:52
• Thanks. The reference was helpful. Aug 2 at 7:30
• Sorry, there was a typo in the code. (Edit mistake: After changing the name of a function, I missed one place it occurred and didn't update it to the new name.) Let me know if it doesn't work. Aug 2 at 14:43

We can examine the Poincaré return map to locate the limit cycles:

ode = {x'[t], y'[t]} ==
({x^2 + y + x y, -((727778 x)/10000) -
10 x^2 + (15 y)/10000 + (22 x y)/10 + (7 y^2)/10} /.
v : x | y :> v[t]);

loop // ClearAll;
loop // Attributes = {HoldRest};
loop[x0_?NumericQ, action_ : Null] := Module[{wp},
wp = Precision[x0];
wp = Replace[
wp, {Infinity | MachinePrecision -> MachinePrecision,
p_ :> p + 8}];
NDSolveValue[
{ode,
{x[0], y[0]} == SetPrecision[{x0, 0}, wp],
With[{evt = If[x0 < 0, y[t] > 0, y[t] < 0]},
WhenEvent[evt, action; "StopIntegration"]
]},
{x, y},
{t, 0, Sign[x0] 100000},
WorkingPrecision -> wp]
];
return // ClearAll;
return[x0_?NumericQ] := Module[{dx},
loop[x0, dx = x[t] - x0];
dx /. Except[_?NumericQ] -> 0.];

GraphicsRow[{
Table[
{x0, return[x0]},
{x0, Range[-4000, -100, 100]~Join~Range[-90, 0, 10]}] //
ListLinePlot,
Table[
{x0, return[x0]},
{x0, (Range[0, 229]/100)(*~Join~Range[17]*)}] // ListLinePlot,
Table[
{x0, return[x0]},
{x0, Range[0, 17]}] // ListLinePlot
}]


Then we can use FindRoot to find the precise initial condition for each:

cycles =
FindRoot[return[x1], {x1, ##}, Method -> "Brent",
WorkingPrecision -> 24] & @@@ {{-3800, 3700}, {1/2, 1}, {2,
5/2}, {15, 17}}
(*
{{x1 -> -3711.56080639431786733064},
{x1 -> 0.683210217398156351542212},
{x1 -> 2.18369982492124843474047},
{x1 -> 15.9627839815812680483170}}
*)


And plot them:

plot = ListLinePlot[
Transpose@Through[loop[#]["ValuesOnGrid"]] & /@ (x1 /. cycles),
PlotRange -> All];

GraphicsRow[{
Show[
VectorPlot[
Evaluate@SolveValues[ode, {x'[t], y'[t]}],
{x[t], -4000, 20}, {y[t], -6000, 15000},
VectorStyle -> Opacity[0.3]],
plot],
Show[
VectorPlot[
Evaluate@SolveValues[ode, {x'[t], y'[t]}],
{x[t], -10, 20}, {y[t], -70, 40}, VectorStyle -> Opacity[0.3]],
plot]
}]


Things can get a little iffy near the origin, so I left the high working precision in. WorkingPrecision -> MachinePrecision yields accurate-looking graphs, despite a few error messages from FindRoot. At MachinePrecision, the solution computed by NDSolve in return[]/loop[] is not accurate enough for FindRoot. But it's accurate enough for a nice phase portrait.

## Update: A more general solution

The initial answer took the x-axis for the section in the return map, and I had to carefully avoid stopping NDSolve when it crossed it at the wrong place. Below we modify the code of return[]/loop[] to allow the specification of the return event (that specifies the Poincaré section). In the present case, taking the line connecting the two equilibria makes a better section than in the initial answer.

We also use the action argument in loop[] to return the period, which was asked for in a comment below. The basic idea is to sow the stopping time with loop[x1, event, Sow[t, "Period"]]. Reap[] then returns both the cycle and the period.

I borrow some code/ideas: For xysample, compare FindAllCrossings[]. The function zc[] ("zero crossings") below is copied from Find zero crossing in a list — it is the same as davidZC2[].

ode = {x'[t], y'[t]} ==
({x^2 + y + x y, -((727778 x)/10000) -
10 x^2 + (15 y)/10000 + (22 x y)/10 + (7 y^2)/10} /.
v : x | y :> v[t]);

equil = Solve[ode /. x'[t] | y'[t] -> 0, Reals];

section =
Det[Join[{{x[t], y[t], 1}}, {x[t], y[t], 1} /. equil]] == 0;
interval = Less @@ Insert[x[t] /. equil // NumericalSort, x[t], 2];
event = section && interval; (* = line segment *)

loop // ClearAll;
loop // Attributes = {HoldRest};
loop[x0_?NumericQ, returnEvent_, action_ : Null] := Module[{wp},
wp = Precision[x0];
wp = Replace[
wp, {Infinity | MachinePrecision -> MachinePrecision,
p_ :> p + 8}];
NDSolveValue[
{ode,
{x[0], y[0]} ==
SetPrecision[{x0,
y[t] /. First@Solve[returnEvent /. x[t] -> x0, y[t]]}, wp],
With[{yp = y'[t] /. First@Solve[ode, {x'[t], y'[t]}]},
WhenEvent[returnEvent,
action;
Sow[x[t], loop];
"StopIntegration"]
]
},
{x, y},
{t, 0, 100000},
WorkingPrecision -> wp]
];
return // ClearAll;
return[x0_?NumericQ, returnEvent_] :=
Reap[loop[x0, returnEvent], loop][[2, 1, 1]] - x0;

zc[l_] :=
SApos_ :>
With[{c = SApos[l]}, {c[[#]], c[[# + 1]]}\[Transpose] &@
SApos@Differences@Sign@l[[c]]]

xysample = First@Cases[
Replace[interval,
Less[a_, __, b_] :> (* save plot in case (to debug) *)
(foo = Plot[return[x1, event], {x1, a, b}])],
Line[p_] :> p,
Infinity];
rootIntervals = xysample[[#, 1]] & /@ zc@ xysample[[All, 2]]
cycles = FindRoot[
return[x1, event], {x1, ##},
Method -> "Brent"(*,WorkingPrecision->24*)
] & @@@ rootIntervals
(*
{{-1.17469,  -1.15925},  {-0.696595, -0.694518},
{-0.514333, -0.512432}, {-0.343169, -0.326669}}

{{x1 -> -1.16671},  {x1 -> -0.6962},
{x1 -> -0.512579}, {x1 -> -0.33355}}
*)

plotCycles = Replace[
Reap[
Transpose@ (* {{x1,..}, {y1,..} --> {{x1,y1},..} *)
N@Through[
loop[#, Evaluate@event, Sow[t, "Period"]][
"ValuesOnGrid"]], "Period"] & /@ (x1 /. cycles) //
Transpose, (* {{c, p},...} -> {cycles, periods} *)
{cyc_, per_} :>
ListLinePlot[cyc,
PlotLegends -> (Row[{"Period = ", #[[1, 1]]}] & /@ N@per),
PlotRange -> All]
] /. Graphics[g_, opts___] :>
Graphics[
{{Magenta, Dashing[{0.02, 0.01}],
Line[{x[t], y[t]} /. equil]},
g},
opts] /.
LineLegend[g_, lbl_, opts___] :>
LineLegend[
Append[g,
Directive[PointSize[1/120], Magenta, Dashing[{0.02, 0.01}],
AbsoluteThickness[1.6]]
],
Append[lbl,
"Section"
],
{opts} /. Verbatim[Rule][LegendMarkers, lm_] :>
LegendMarkers -> Append[lm, {False, Automatic}]];

finalPlot = Row[{
Show[
VectorPlot[
Evaluate@SolveValues[ode, {x'[t], y'[t]}],
{x[t], -4000, 50}, {y[t], -6000, 15000},
VectorStyle -> Opacity[0.3]],
plotCycles,
ImageSize -> {Automatic, 240}],
Spacer[20],
Show[
VectorPlot[
Evaluate@SolveValues[ode, {x'[t], y'[t]}],
{x[t], -10, 20}, {y[t], -70, 40},
VectorStyle -> Opacity[0.3]],
plotCycles,
ImageSize -> {Automatic, 240}]
}];
finalPlot =
With[{leg = First@Cases[finalPlot,
Legended[_, l_] :> l, Infinity]},
Replace[finalPlot,
gr_ :> Legended[gr /. Legended[g_, _] :> g, leg]]]


• Although it was not a question of OP but what is the period of each cycle, how to compute it? I mean value of parameter t for which we return to the exact same point we started. Aug 2 at 15:40
• @azerbajdzan Map[With[{y = #[[2]]}, FindRoot[y[t], {t, y["Domain"][[1, 2]], Sequence @@ y["Domain"][[1]]}]] &, loop[x1] /. cycles] is one way to get the periods. Another would be to save the time of the WhenEvent in loop[]. While the event action is to stop integration, the solution returned will include the last full step, which usually goes a little beyond the event. Refactoring loop this way would make it somewhat complicated. Since WhenEvent already does the FindRoot work, it also seems a slightly more efficient way. Aug 2 at 16:39
• @azerbajdzan See updated answer. Aug 2 at 18:59
• Great job! And now find all rational points on these curves, derive a group law and use it for a new cryptography method ;-) Aug 2 at 21:11
• Hi @MichaelE2, I hope you are doing well. This year, I plan to write a scientific paper in which I will employ your code. Based on that, it would be a pleasure to acknowledge you. Could you tell me your real name and institution? Thanks once more and God bless you. Aug 3 at 16:40