# How to determine the frequency of oscillations in system of three ODEs?

I want to determine the frequency of oscillations in a system of three stiff ODEs (Oregonator model). That model describes a chemical oscillator.

I have a slightly more advanced model of the default or regular Oregonator. It consists of three ODEs:

ode1=ε*x'[t]==x[t](1-x[t])-2*(x[t]-μ)/(x[t]+μ)*(r*α1*y[t]+q*((α2*z[t])/(1-z[t])));
ode2=ξ1*y'[t]==x[t](1+β2*z[t])-α1*y[t]-((x[t]*(1+β1*y[t]+β2*z[t]))/((1-z[t])+η*(1-y[t])))*(1-z[t]);
ode3=ξ2*z'[t]==x[t]*(1+β1*y[t])-(α2*z[t])/(1-z[t])-((x[t]*(1+β1*y[t]+β2*z[t]))/((1-z[t])+η*(1-y[t])))*η*(1-y[t]);

with the initial (example) conditions ic

ic2 = {x[0] == .5, y[0] == 0.1, z[0] == 0.3};

I use NDSolveValue for this:

{xSol, ySol, zSol} = NDSolveValue[{ode1, ode2, ode3, ic2}, {x, y, z}, {t, 0, 200}]

This looks like this:

So far, so fine. I now need to determine the frequency of the oscillations in this model with three ODEs.

I found this related question, but that only features a single ODE. And as I'm really a Mathematica novice, I also didn't understand how the Reap and Sow worked.

The suggested solution was as follows:

pts =
Reap[s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1,
WhenEvent[y'[x] == 0, Sow[x]]}, {y, y'}, {x, 0, 30}]][[2, 1]]

(* Out[290]= {0.448211158984, 4.6399193764, 7.44068279785, 10.953122261, \
13.8722260952, 17.2486864443, 20.2244048853, 23.5386505821, \
26.5478466115, 29.8261176372} *)

Plot[{Evaluate[y[x] /. s], Evaluate[y'[x] /. s]}, {x, 0, 30},
PlotRange -> All]

and then finding the differences:

diffs = Differences[pts, 1, 2]

(* Out[288]= {6.99247163887, 6.31320288463, 6.43154329733, \
6.29556418327, 6.35217879014, 6.28996413777, 6.32344172616, \
6.28746705515} *)

Mean[diffs]

(* Out[289]= 6.41072921417 *)

This looks exactly what I need, but I don't know how to apply this to my three ODEs? I preferably want to keep the initial conditions, ic, in a separate variable like I now have.

Can anyone show me how to modify the above solution so that it works with my system? I want to determine the frequency separately for x[t], y[t] and z[t]. If people have a different solution than proposed in the related question, you're of course very welcome!

Many thanks in advance!

Current file with all the needed variables, ODEs and my use of NDSolveValue. You can copy and modify this.

• Another related Q&A: mathematica.stackexchange.com/questions/129993/… – Michael E2 Oct 14 '20 at 18:34
• To determine the frequency of oscillatory data the FFT: Fourier is commonly used. – Daniel Huber Oct 14 '20 at 19:46
• @MichaelE2 thanks! Will look into that one. – ralphjsmit Oct 15 '20 at 5:41
• Hi @DanielHuber thanks for mentioning that! Is there a very simple beginner example in which the frequency is calculated? And some explanation about FFT? – ralphjsmit Oct 15 '20 at 5:52
• It is not too easy if you never heard about it, but it is a standard tool in data analysis. There are zillion of articles on the web,.You may e.g. look at: earlevel.com/main/2002/08/31/a-gentle-introduction-to-the-fft. You should play around with signals of the form Sin[f1 t]+ Sin[f2 t]+.. and see what the output of Fourier is. Note that in general the output is complex and that Fourier returns a list: first element: zero frequency component, second: base frequency, first overtone,.. up to the highest frequency, then highest negative frequency, then the second highest,.. – Daniel Huber Oct 15 '20 at 8:47

Using my EcoEvo package as in my answer here.

First, you need to install it with

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

Then, load the package and define your model:

<< EcoEvo`;

SetModel[{
Aux[x] -> {Equation :> (x[t] (1 - x[t]) - 2*(x[t] - μ)/(x[t] + μ)*(r*α1*y[t] + q*((α2*z[t])/(1 - z[t]))))/ε},
Aux[y] -> {Equation :> (x[t] (1 + β2*z[t]) - α1*y[t] - ((x[t]*(1 + β1*y[t] + β2*z[t]))/((1 - z[t]) + η*(1 - y[t])))*(1 - z[t]))/ξ1},
Aux[z] -> {Equation :> (x[t]*(1 + β1*y[t]) - (α2*z[t])/(1 - z[t]) - ((x[t]*(1 + β1*y[t] + β2*z[t]))/((1 - z[t]) + η*(1 - y[t])))*η*(1 - y[t]))/ξ2}
}];

Simulate for 40 time steps to get on the limit cycle:

sol = EcoSim[{x -> 0.5, y -> 0.1, z -> 0.3}, 40];
GraphicsRow[{PlotDynamics[sol, x], PlotDynamics[sol, y], PlotDynamics[sol, z]}]

Find the limit cycle using FindEcoCycle:

ec = FindEcoCycle[FinalSlice[sol]];
GraphicsRow[{PlotDynamics[ec, x], PlotDynamics[ec, y], PlotDynamics[ec, z]}]

Make sure the initial and final values are the same:

InitialSlice[ec]
FinalSlice[ec]

(* {x -> 0.617907, y -> 0.522312, z -> 0.989451} *)
(* {x -> 0.617907, y -> 0.522312, z -> 0.989451} *)

Finally, the period can be extracted as the final time of ec:

FinalTime[ec]
(* 2.71597 *)

To answer E. Chan-López's question, to find equilibria use SolveEcoEq.

eq = SolveEcoEq[]
(* {{x -> 0, y -> 0, z -> 0}, {x -> 0.000395349, y -> 1., z -> 1.}, {x -> 0.0279224, y -> 0.06391, z -> 0.881397}, {x -> 0.0279224, y -> 0.124848, z -> 1.08789}, {x -> 0.0279224, y -> 0.988952, z -> 1.00342}} *)

Apparently they're all unstable:

EcoStableQ[eq]
(* {False, False, False, False, False} *)
• Hey @ChrisK thanks for your answer! I’ll try it out (at last Saturday) and let you know! – ralphjsmit Oct 15 '20 at 5:53
• Hi Chris, thanks, it works nicely👍 – ralphjsmit Oct 15 '20 at 15:35
• Hi! How many equilibrium points can you calculate in your model? – E. Chan-López Feb 8 at 5:05
• @E.Chan-López I added the equilibria to my answer - looks like there are five, all of which are unstable. – Chris K Feb 8 at 15:59
• @ChrisK I added a calculation to show the existence of a Hopf bifurcation at a non-trivial equilibrium point. – E. Chan-López Feb 9 at 2:09

Contributing to Chris's analysis:

Taking $$q=\displaystyle\frac{33}{4}$$, $$r=\displaystyle\frac{1}{90}$$, $$\alpha_{1}=\displaystyle\frac{64}{33}$$, $$\alpha_{1}=\displaystyle\frac{1}{22}$$, $$\beta_{1}=24$$, $$\beta_{2}=3/4$$, $$\xi_{1}=1$$, $$\xi_{2}=1$$, $$\eta=32$$, $$\mu=\displaystyle\frac{16}{1517}$$, and $$\varepsilon$$ as a free parameter, we obtain five non-trivial equilibrium points.

The system:

f1 = 1/ε (x (1 - x) - 2*(x - μ)/(x + μ)*(r*α1*y + q*((α2*z)/(1 - z))));
f2 = 1/ξ1 (x (1 + β2*z) - α1*y - ((x*(1 + β1*y + β2*z))/((1 - z) + η*(1 - y)))*(1- z));
f3 = 1/ξ2 (x*(1 + β1*y) - (α2*z)/(1 - z) - ((x*(1 + β1*y + β2*z))/((1 - z) + η*(1 - y)))*η*(1 - y));
F = {f1, f2, f3};
X = {x, y, z};

The Jacobian matrix:

J = FullSimplify[D[F, {X}]];

The parameter values and non-trivial equilibrium points:

q = 33/4;
r = 1/90;
α1 = 64/33;
α2 = 1/22;
β1 = 24;
β2 = 3/4;
ξ1 = 1;
ξ2 = 1;
η = 32;
μ = 16/1517;

X0 = N[Normal[Solve[F == 0 && Variables[F] > 0, X]]]
(*{{x -> 1/2, y -> 1/4, z -> 1/4},
{x -> 13.3111, y -> 12.0089, z -> 1.00458},
{x -> 0.0105402, y -> 0.999138, z -> 1.02416},
{x -> 2.38478, y -> 1.33557, z -> 1.28417},
{x -> 0.0103128, y -> 0.0337578, z -> 5.68495}}*)

The linear approximation at the point $$X_{01}=(1/2, 1/4, 1/4)$$ and its characteristic polynomial:

J01 = FullSimplify[J /. X01];
polJ01 = -Collect[CharacteristicPolynomial[J01, λ], λ, FullSimplify]
(*polJ0 = A0 λ^3 + A1 λ^2 + A2 λ^1 + A3;
A0 = 1; A1 = (8 (66748 + 9580565 ε))/(25302915 ε);
A2 = (4 (268141961 + 2208997920 ε))/(7514965755 ε);
A3 = 571980408512/(743981609745 ε) *)

According to the Routh-Hurwitz criterion, $$X_{01}$$ is stable if and only if

Reduce[A1 > 0 && A2 > 0 && A3 > 0 && A1 A2 - A3 > 0 && ε > 0]
(*0 < ε < (239125592953 - Sqrt[31922603976781931264305])/5465060854080 ||
ε > (239125592953 + Sqrt[31922603976781931264305])/5465060854080*)

The last result means that we have two critical Hopf values for $$\varepsilon$$ at point $$X_{01}$$. According to the Liu criterion for the Hopf bifurcation, these critical values are calculated as follows

ε0 = Solve[A1 A2 - A3 == 0 && ε > 0, ε]];
(*{{ε -> (239125592953 - Sqrt[31922603976781931264305])/5465060854080},
{ε -> (239125592953 + Sqrt[31922603976781931264305])/5465060854080}}*)

At the other equilibrium points, the analysis is similar, except for the second equilibrium, since it is unstable.