Plotting Limit Cycle

Question

I have Lotka Volterra Equations where the prey equation is modified. The Equations are

x'[t]=a*x[t]*(1-x[t]/k)-b*x[t]*y[t]
y'[t]=-c*y[t]+d*x[t]*y[t]

To plot the limit cycles, i used NDSolve to solve the coupled nonlinear differential equations for x[t] and y[t]. I then used ParametricPlot to plot for x[t] and y[t] with respect to t, but I'm not getting appropriate limit cycles.

I did a steady state analysis of the equations and got certain conditions on the value of a, b, c, d and k. The equation has a non-trivial steady state, at $(x,y)=(c/d,(a/b)-(ac/bdk))$. The commands are:

Equations = {x'[t] == a*x[t]*(1 - x[t]/k) - b*x[t]*y[t], x[0] == 2, y[0] == 1, y'[t] == -c*y[t] + d*x[t]*y[t]}; 
s = NDSolve[Equations, {x, y}, {t, 0, time}];
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 30}] 

For the value of a = 2; b = 1; c = 2.5; d = 1.2; k = 2.8, I'm supposed to get a unstable spiral. But instead, the trajectory approaches the steady state of $(x,y) = (2.08333 , 0.511905)$.

Any help will be greatly appreciated.

Answer 1

Hopf bifurcation analysis

The differential system:

f1[x_,y_]:=a x (1 - x/k) - b x y; 
f2[x_,y_]:=-c y + d x y; 
F[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@{f1[x,y],f2[x,y]};
X={x, y};
μ={a,b,c,d,k};

$$ \begin{align} &\dot{x}=a x\left(1-\frac{x}{k} \right)- bxy\\ &\dot{y}= dxy - cy \end{align} $$

The Jacobian matrix:

J[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@D[F[X,μ],{X}]

The non-trivial equilibrium point:

X0=Normal[Simplify[SolveValues[F[X,μ]==0&&Variables[F[X,μ]]>0,X]]][[1]];
MatrixForm@X0

$$ \begin{align} P_{0}(x,y)=\left(\frac{c}{d} ,\frac{a}{b}\left(1-\frac{c}{dk} \right)\right) \end{align} $$ The linear approximation at $P_{0}$ (coexistence equilibrium point):

J0=Simplify@J[X0,μ];
MatrixForm@J0

$$ \begin{align} J(P_{0})=\left( \begin{array}{cc} \hspace{-0.25cm}-\displaystyle\frac{a c}{d k} & -\displaystyle\frac{b c}{d}\hspace{0.3cm} \\\\ \hspace{0.2cm}\displaystyle\frac{a (d k-c)}{b k} &\hspace{0.2cm} 0 \\ \end{array} \right) \end{align} $$ Under the Hopf bifurcation conditions, $\text{tr}(J(P_{0},\mu_{0}))=0$ and $\text{det}(J(P_{0},\mu_{0}))>0$, where $\mu_{0}$ is the critical bifurcation value for some parameter of our system. In our case, the parameters are strictly positive and $\text{tr}(J(P_{0}))$ cannot be zero. Therefore, Hopf bifurcation not take place at $P_{0}$. The non-trivial equilibrium $P_{0}$ is always locally stable and the only condition that must be fulfilled is given by the following inequality $$ \frac{c}{d k}<1 $$ Code for time series and phase portrait:

Time series

s = ParametricNDSolve[{x'[t] == a x[t] (1 - x[t]/k) - bx[t]*y[t], 
y'[t] == -c y[t] + d x[t]*y[t], x[0] == 11/5, y[0] == 4/5}, {x, y}, {t, 0, 1000}, {a, b, c, d, k}];
Plot[Evaluate[x[1/4, 1/3, 1/2, 1/4, 10][t] /. s], {t, 0, 300}, PlotRange -> All, PlotPoints -> 500, 
PlotStyle -> {Blue, Thickness[0.003]},AxesStyle -> Directive[Black, Small], Background -> Lighter[Gray, 0.95]]

Phase portrait:

ParametricPlot[Evaluate[{x[1/4, 1/3, 1/2, 1/4, 10][t], y[1/4, 1/3, 1/2, 1/4, 10][t]} /. s], 
{t, 0, 300}, PlotRange -> All, PlotPoints -> 500, PlotStyle -> {Blue, Thickness[0.003]}, 
AxesStyle -> Directive[Black, Small],Background -> Lighter[Gray, 0.95]]

Example: Hopf bifurcation in the Brusselator system

Calculation of the first Lyapunov coefficient

The Brusselator system is given by: $$ \begin{align} &\dot{x}=\alpha-(\beta+1)x + x^2 y\\ &\dot{y}= \beta x - x^2 y \end{align} $$

Assuming $\alpha> 0$ fixed and taking $\beta$ as a bifurcation parameter, we show that at $\beta = 1 + \alpha^2$ the system exhibits a supercritical Hopf bifurcation.

The Brusselator system code:

f1[x_, y_] := α - (β + 1) x + x^2 y;
f2[x_, y_] := β x - x^2 y;
F[{x_, y_}, {α_, β_}] := Evaluate@{f1[x, y], f2[x, y]};
X = {x, y};
μ = {α, β};
U = {u, v};
R = {r, s};

The Jacobian matrix and its transpose:

J[{x_, y_}, {α_, β_}] = D[F[X, μ], {X}];
Jt[{x_, y_}, {α_, β_}] = Transpose[J[X, μ]];
MatrixForm[J[X, μ]]
MatrixForm[Jt[X, μ]]

Stability analysis (Routh-Hurwitz criterion):

X0[{α_, β_}] = SolveValues[F[X, μ] == 0, X][[1]]
polJX0 = Collect[CharacteristicPolynomial[J[X0[μ], μ], λ], λ,Simplify];
a0 = CoefficientList[polJX0, λ][[3]];
a1 = CoefficientList[polJX0, λ][[2]];
a2 = CoefficientList[polJX0, λ][[1]];
Reduce[a1 > 0 && a2 > 0 && α > 0 && β > 0, β]
(*α > 0 && 0 < β < 1 + α^2*)

Note that $a_{1}=0$ if and only if $β=β_{0}$, where $β_{0}=1+α^2$. Then, the Brusselator is locally asymptotically stable at $X_{0}(\mu)$ for $β<β_{0}$ and locally asymptotically unstable for $β>β_{0}$ (appears a stable limit cycle surrounded the unstable equilibrium point). We verify the previous conclusion (transversality condition) with the sign of the following derivative:

D[-a1, β]
(*1*)

The analysis at the critical bifurcation value $β_{0}$:

Solve the following system of equations: $$ \left\{\begin{align} F\left((x,y),(\alpha,\beta)\right) &=0, \\ \operatorname{tr}(J((x,y),(\alpha,\beta))) &=0, \end{align}\right. $$ for $(x,y,\beta)$ and we must check that det $\operatorname{det}J((x,y),(\alpha,\beta))>0$ when $\beta = \beta_{0}$ for the solution found, where $\beta_{0}$ is the Hopf critical bifurcation value.

The code for the above system of equations:

X0μ0 = Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, β}], 1], {3}]
μ0 = Prepend[Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, β}], 1], {{1}, {2}}], α]
Det[J[X0μ0, μ0]]

Here, the Hopf critical bifurcation value is $\beta_{0}=1+\alpha^2$ and $\operatorname{det}J((x,y),(\alpha,\beta_{0}))=\alpha^2>0$. Thus, the Brusselator at $\beta_{0}=1+\alpha^2$ has the equilibrium $$ \begin{align} X_{0}(\mu_{0})=\left(\alpha, \displaystyle\frac{1+\alpha^2}{\alpha} \right) \end{align} $$ and the linear approximation at $X_{0}(\mu_{0})$ has purely imaginary eigenvalues $\lambda_{1,2}=\pm \omega i$, $\omega=\alpha$.

The code for the linear approximation and its transpose at $X_{0}(\mu_{0})$:

α= ω;
JX0μ0 = Simplify@J[X0μ0, μ0];
JtX0μ0 = Simplify@Transpose@JX0μ0;
MatrixForm@JX0μ0
MatrixForm@JtX0μ0
Eigenvalues[JX0μ0]

The next step is to translate the equilibrium $X_{0}(\mu_{0})$ to the origin of coordinates:

bb = {0, 0};
F0[{x_, y_}, {α_, β_}] = Collect[Expand@F[X + X0μ0, μ0], {x, x^2, y, y^2, x y, x^2 y},Factor]
MatrixForm@F0[bb,μ]

Now, to obtain the normal form of the Hopf bifurcation, we need the Taylor expansion of the third order for $F_{0}((x,y),(\alpha,\beta))$:

(*Rank 3 tensor*)
D2[{x_, y_}, {α_, β_}] = Simplify@D[F0[X, μ], {X, 2}]
D2X0μ0 = Simplify@D2[bb, μ0]
(*Rank 4 tensor*)
D3[{x_, y_}, {α_, β_}] = Simplify@D[F0[X, μ], {X, 3}]
D3X0μ0 = Simplify@D3[bb, μ0]

Multilinear forms:

(*Bilinear form*)
BB[{x_, y_}, {u_, v_}] = Collect[Expand[D2X0μ0.X.U], {u x, v y, v x, v y}, FullSimplify];
MatrixForm[BB[{x, y}, {u, v}]]
(*Trilinear form*)
CC[{x_, y_}, {u_, v_}, {r_, s_}] = Collect[Expand[D3X0μ0.X.U.R], {r u x, r u y, r v x, r v y, s u x, s v x, s u y, s v y}, FullSimplify]
MatrixForm[CC[{x, y}, {u, v}, {r, s}]]

We verify that the first three terms of the Taylor series expansion of $F_{0}((x,y),(\alpha,\beta))$ are correct:

A=JX0μ0; (*linear approximation*)
MatrixForm@FullSimplify[F0[X, μ] - (A.X + 1/2! BB[X, X] + 1/3! CC[X, X, X]) /. {x -> t x, y -> t y}]
(*{0,0}*)

Now, we compute the critical eigenvectors of $J((0,0),\mu_{0})$ and its transpose:

(*Eigenvectors of A=J[X0,μ0]*)
vp = ComplexExpand@Eigenvectors[A]
q = vp[[2]];
qc = vp[[1]];
MatrixForm@q
MatrixForm@qc
MatrixForm@Simplify[A.q - I ω q]
(*Eigenvectors of Transpose[A]*)
At=JtX0μ0;
vpt = ComplexExpand[Eigenvectors[At]]
(*Normalization constant*)
cn = ComplexExpand[Conjugate[vp[[1]] . vpt[[1]]]]
p = Expand@Simplify[vpt[[2]]/cn];
pc = ComplexExpand[Conjugate[p]];
MatrixForm@p
MatrixForm@pc
Simplify[At.p-I ω p]

We verify the normalization condition $\langle p,q\rangle=1$

Simplify@(p.q)
(*1*)

Finally, we compute the first Lyapunov coefficient: $$ \begin{align} l_1(0,\mu_{0})= &\frac{1}{2\omega_0} {\rm Re}\left[\langle p,C(q,q,\bar{q}) \rangle - 2 \langle p, B(q,A_0^{-1}B(q,\bar{q}))\rangle +\\\hspace{0.5cm} \langle p, B(\bar{q},(2i\omega_0 I_n-A_0)^{-1}B(q,q))\rangle \right] \end{align} $$

Before to calculate $l_1(0,\mu_{0})$, we clean $\alpha$

Clear[α]
ω0=α;

The code for the first Lyapunov coefficient:

Factor@ComplexExpand[Re[1/(2 ω) (p.CC[q, q, qc] - 2 (p.BB[q, Inverse[A].BB[q, qc]]) + p.BB[qc, Inverse[2 I ω*IdentityMatrix[2] - A].BB[q, q]])] /. ω -> ω0]
(*-((2 + α^2)/(2 α (1 + α^2)))*)

$$ \begin{align} l_1(0,\mu_{0})=-\frac{\alpha^2+2}{2 \alpha \left(\alpha^2+1\right)} \end{align} $$ The first Lyapunov coefficient is clearly negative for all positive $\alpha$. Thus, the Hopf bifurcation is nondegenerate and always supercritical.

The above expression is the result that Kuznetsov arrives at on page 105 of his book (see Elements of Applied Bifurcation Theory).

Limit cycle:

For more details, see: Andronov-Hopf bifurcation.

Answer 2

I'm afraid that you have calculated the stability incorrectly. Here is the Jacobian of your system:

a = 2; b = 1; c = 2.5; d = 1.2; k = 2.8;
jac[x_, y_] := {{D[a x*(1 - x/k) - b x*y, x], D[a x*(1 - x/k) - b x*y, y]}, 
               {D[-c y + d x*y, x], D[-c y + d x*y, y]}};

At the equilibrium, this is:

jacEq = jac[x, y] //. {x -> 2.08333, y -> 0.511905}

The eigenvalues of this are:

Re[Eigenvalues[jacEq]]
{-0.744047, -0.744047}

so the system is stable about this equilibrium.