# Estimation of parameters of limit cycles for systems of high-order differential equations (n> = 3)

There is a system of differential equations:

Then, call the limit cycle the projection of the phase trajectory onto the plane in a pairwise combination of state variables ($$x-y,y-z,x-z$$).

where $$x,y,z$$ - state variables, $$a,b,c$$ - constants.

Is it possible to use Mathematica to estimate the amplitude and frequency of the limit cycle? (it is possible by approximate numerical methods, most importantly, not graphical).

I did like this: 1. Using NDSolve, I solve the system of differential equations numerically.

s = NDSolve[{x'[t] == -3 (x[t] - y[t]),
y'[t] == -x[t] z[t] + 26.5 x[t] - y[t], z'[t] == x[t] y[t] - z[t],
x[0] == z[0] == 0.1, y[0] == 0.25}, {x, y, z}, {t, 0, 400}]
1. Using ParametricPlot, I build a phase plane for a pairwise combination of state variables (see Figure 1 for an $$x-y$$ pair).

ParametricPlot[Evaluate[{x[t], y[t]} /. First[%]], {t, 0, 100}]

2. Using the Plot command, I build a graph for the state variable in time and try to estimate the frequency of the alternating signal from the graph. (see Figure 1 for an $$x$$ variable).

Plot[Evaluate[x[t] /. s], {t, 0, 100}]

EDIT:

After several hours of calculations and on the advice of one of the users, I applied data sampling and Fourier expansion with the construction of a frequency spectrum.

xsol[t_] := x[t] /. s[[1]]

xdis = Table[xsol[i], {i, 0, 100, 0.1}];

ListPlot[xdis]

fft = Fourier[xdis, FourierParameters -> {1, -1}];

ListLinePlot[shortFFT = Abs[fft[[5 ;; 400]]], PlotRange -> All]

f = Abs[Fourier[xdis]];

peaksize = Last[TakeLargest[f, 2]];

peaks = Flatten[Position[f, i_ /; i >= peaksize]];

pos = First[peaks];

Show[ListPlot[f], Graphics[{Red, Point[{pos, f[[pos]]}]}],
PlotRange -> All]

n = 100/0.1 + 1;

fr = Abs[Fourier[xdis Exp[2 Pi I (pos - 2) N[Range[0, n - 1]]/n],
FourierParameters -> {0, 2/n}]];

frpos = Position[fr, Max[fr]][[1, 1]]

Show[ListPlot[fr], Graphics[{Red, Point[{frpos, fr[[frpos]]}]}],
PlotRange -> All]

N[n/(pos - 2 + 2 (frpos - 1)/n)]

Fourier - > Applications - > Frequency Identification

This code gives an estimate of a period of ~ 564 sec and a frequency of 1 / T ~ 0.002 Hz. Which, of course, does not look like the results of NDSolve.

EDIT №2:

There is my code for Lorenz System. Nothing unusual, only classical continuous Fourier series.

In[49]:= pars = {n = 15, T = 20, \[Omega] = 2 Pi/T}

Out[49]= {15, 20, \[Pi]/10}

In[61]:= s =
NDSolve[{x'[t] == -3 (x[t] - y[t]),
y'[t] == -x[t] z[t] + 26.5 x[t] - y[t], z'[t] == x[t] y[t] - z[t],
x[0] == z[0] == 0.1, y[0] == 0.25}, {x, y, z}, {t, 0, 20}]

In[66]:= Plot[Evaluate[x[t] /. s], {t, 0, T}, PlotRange -> Full]

In[67]:= ifun = First[x /. s]

In[68]:= a0 = 2/T NIntegrate[ifun[t], {t, 0, T}]

Out[68]= -4.74859

In[69]:= f =
a0/2 + Sum[
2/T NIntegrate[
ifun[t] Cos[\[Omega] k t], {t, 0, T}] Cos[\[Omega] k t] +
2/T NIntegrate[
ifun[t] Sin[\[Omega] k t], {t, 0, T}] Sin[\[Omega] k t], {k, 1,
n}];

In[70]:= Plot[{ifun[t], f}, {t, 0, T}, PlotRange -> Full]

QUESTION: Is it possible to speed up this code, for example, apply a faster algorithm of numerical integration?

• Probably. What have you tried? Commented May 26, 2020 at 16:40
• Include the code you used with NDSolve etc. If you give people something to start from, they are much more likely to help. If you have a starting point, even just the equations formatted in MMA format rather than having to type them back in, please share it by editing the question. In this case, please share the NDSolve and the ParametricPlot, together with any definition of constants you used, etc. Commented May 26, 2020 at 18:03
• As you can see, I did it "manually." This is not very convenient and most likely not correct. And I do not know how to evaluate the parameters of limit cycles in Mathematica.
– ayr
Commented May 27, 2020 at 4:25
• Interesting question. Perhaps Fourier on a (large) discrete equi-spaced sample set might be useful here? Commented May 27, 2020 at 16:03
• Do you mean to expand the data from NDSolve in a discrete Fourier series and analyze them?
– ayr
Commented May 27, 2020 at 16:11

These aren't actual limit cycles and what you're looking for has a fuzzy definition (notice how the amplitude increases in each pass). However this is still fun to play with, so let's see what we can find. I learned a lot of theory and practical tips from reading The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors by Colin Sparrow.

Although this isn't an ecological model, my EcoEvo package has some functions that are useful, so I'll make use of it. To install it, use:

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]

Load the package and define the model:

<<EcoEvo;

SetModel[{
Aux[x] -> {Equation :> σ (y - x)},
Aux[y] -> {Equation :> r x - y - x z},
Aux[z] -> {Equation :> x y - b z},
Parameters :> {σ, r}
}]

σ = 3; r = 26.5; b = 1;

First, simulate for 400 time steps to get on the attractor:

s = EcoSim[{x -> 0.1, y -> 0.25, z -> 0.1}, 400];
PlotDynamics[FinalSlice[s, 100], x]

There are three equilibria. Find them and plot with the attractor in phase space:

eq = SolveEcoEq[]
(* {{x -> 0, y -> 0, z -> 0}, {x -> -5.04975, y -> -5.04975, z -> 25.5},
{x -> 5.04975, y -> 5.04975, z -> 25.5}} *)

Show[
RuleListPlot[eq[[2 ;; 3]]],
RuleListPlot[FinalSlice[s, 100], PlotPoints -> 200]
]

To get an approximate idea of the "period" in each wing, we could find the eigenvalues of the equilibria and calculate from their imaginary parts.

EcoEigenvalues[eq[[2]]]
(* {0.0495488 + 5.47749 I, 0.0495488 - 5.47749 I, -5.0991} *)

2 π/5.477486069462778
(* 1.14709 *)

Looks reasonable at least!

Now for the fun part. There are a bunch of periodic orbits in there, but they're all unstable. To find them, we'll first make a Poincaré section at z = r - 1 using WhenEvent.

ps = Reap[
EcoSim[{x -> 0.1, y -> 0.25, z -> 0.1}, 10000,
WhenEvents -> {WhenEvent[z[t] < r - 1, Sow[{t, x[t], y[t], z[t]}]]}]
][[2, 1]];
ListPlot[ps[[All, 2 ;; 3]], PlotStyle -> PointSize[0.001], AxesLabel -> {x, y}]

Plotting the return map of x[t]:

Show[
ListPlot[Partition[ps[[All, 2]], 2, 1], PlotStyle -> PointSize[0.002]],
Plot[x, {x, -4, 4}], AxesLabel -> {x[t], x[t + 1]}
]

Now, to find an unstable limit cycle we will use Newton's method, which requires a good initial guess. To get one, we'll scan through the Poincaré section for near misses to an n-peak cycle. First, a 2-peak cycle:

n = 2;
ics = Table[
If[Abs[ps[[i, 2]] - ps[[i + n, 2]]] < 10^-3,
{Thread[{x, y, z} -> Mean[{ps[[i + n, 2 ;; 4]], ps[[i + n, 2 ;; 4]]}]],
Period -> ps[[i + n, 1]] - ps[[i, 1]]},
Nothing
]
, {i, Length[ps] - n}]
(* {{{x -> 2.66067, y -> -3.507, z -> 25.5}, Period -> 2.78729},
{{x -> 2.6621, y -> -3.50162, z -> 25.5}, Period -> 2.78703}} *)

Then use that initial guess in my FindEcoCycle:

lc2 = FindEcoCycle[Sequence @@ ics[[1]], Method -> "FindRoot"];
RuleListPlot[lc2]

We can verify that it is an unstable cycle by calculating its Floquet multipliers:

EcoEigenvalues[lc2, Multipliers -> True]
(* {4.9097, 0.999996, 1.80293*10^-7} *)

Greater than 1 means unstable.

We can do the same for 3-peak cycles and 4-peak cycles (there are two kinds: LLRR and LLLR).

lc3 = FindEcoCycle[{x -> -2.783111360797086, y -> 3.054488804793205, z -> 25.5},
Period -> 4.12801717847924, Method -> "FindRoot"];
lc4 = FindEcoCycle[{x -> -3.256897972998302, y -> 1.2980664828293624, z -> 25.5},
Period -> 5.419166840014896, Method -> "FindRoot"];
lc4b = FindEcoCycle[{x -> 3.0826624273683545, y -> -1.9196610967506293, z -> 25.5},
Period -> 5.517358070569571`, Method -> "FindRoot"];
RuleListPlot[{lc3, lc4, lc4b}]

Notice that the periods are close to n multiples of ~1.4, another way to estimate the time spend in each wing.

Finally, put these unstable orbits together with the attractor:

RuleListPlot[{FinalSlice[s, 100], lc2, lc3, lc4, lc4b},
PlotPoints -> 200, PlotStyle -> {{Thin, Gray}, Red, Orange, Green, Blue}]

• This, of course, is a very advanced answer. I have not studied it completely yet, but it seems to be what we need. After a detailed study I will write questions if they arise.I can’t immediately formulate them yet, until I try this calculation myself.
– ayr
Commented May 31, 2020 at 6:25
• @chris. This is really beautiful, and exactly what I'm looking for! The code no longer seems to work with the current version of the packlet, is that correct, or am I doing something wrong? Commented Oct 4, 2022 at 15:13
• @JonathanShock Yes I made a change to the EcoEvo package to remove the [t]s in the SetModel function. After updating that part, it seems to still work. Sorry about the confusion! Commented Oct 4, 2022 at 15:33
• Perfect, cheers! Commented Oct 4, 2022 at 20:33
• @ChrisK. Thank you for great sharing! And I have a question here, how to find an unstable limit cycle in the fixed point region? In other words, are there any efficient ways to choose an appropriate initial state in non-chaos region? Commented Nov 29, 2022 at 16:34