# Lyapunov exponent and stability of limit cycles

I have a four dimensional dynamical system where I have a Hopf bifurcation. My aim is to figure out if limit cycles are stable or not.

The calibration in the system is

paramFinal = {ρ -> 0.85, u -> 0.95, h -> 60, a -> 1, b -> 0.1, δ -> 0.001, γ -> 0.005, ψ -> 60, θ -> 0.04586401037966653, β ->  0.8772372953238555};

KK4 = 33.5081
SS4 = 0.212842
lambdaSS4 = 0.19493
muSS4 = 4.02085
Subscript[a, 1] = 151.467


I have the following system :

kdot = (a - δ) k - (h (2 - u) - λ Subscript[a, 1]);
sdot = s (1 - s) - γ a k;
lambdadot = ((ρ + θ) - (a - δ)) λ + μ γ a;
mudot = (ρ + θ) μ - μ (1 - 2 s) - θ ψ b s^-1 ;


To be able to have Jacobian matrix, I write

Subscript[f, 11] = D[kdot, k];
Subscript[f, 12] = D[kdot, s];
Subscript[f, 13] = D[kdot, λ];
Subscript[f, 14] = D[kdot, μ];
Subscript[f, 21] = D[sdot, k];
Subscript[f, 22] = D[sdot, s];
Subscript[f, 23] = D[sdot, λ];
Subscript[f, 24] = D[sdot, μ];
Subscript[f, 41] = D[mudot, k];
Subscript[f, 42] = D[mudot, s];
Subscript[f, 43] = D[mudot, λ];
Subscript[f, 44] = D[mudot, μ];


To write the fourth order characteristic polynominal equation, I write down the the sum of minors of matrix and determinant.

Subscript[b, 1] = (Subscript[f, 11] + Subscript[f, 22] + Subscript[f, 33] + Subscript[f, 44]) /. paramFinal /. k -> KK4 /. s ->  SS4 /. μ ->  muSS4 /. λ ->   lambdaSS4

Subscript[b, 2] = Det[{{Subscript[f, 11], Subscript[f, 12]}, {Subscript[f, 21], Subscript[f, 22]}}] + Det[{{Subscript[f, 11], Subscript[f, 13]}, {Subscript[f, 31], Subscript[f, 33]}}] + Det[{{Subscript[f, 11], Subscript[f, 14]}, {Subscript[f, 41], Subscript[f, 44]}}] + Det[{{Subscript[f, 22], Subscript[f, 23]}, {Subscript[f, 32], Subscript[f, 33]}}] + Det[{{Subscript[f, 22], Subscript[f, 24]}, {Subscript[f, 42], Subscript[f, 44]}}] + Det[{{Subscript[f, 33], Subscript[f, 34]}, {Subscript[f, 43], Subscript[f,44]}}]/. paramFinal /. k -> KK4 /. s ->  SS4 /. μ ->  muSS4 /. [Lambda] ->   lambdaSS4;

Subscript[b, 3] = (Det[{{Subscript[f, 22], Subscript[f, 23], Subscript[f, 24]}, {Subscript[f, 32], Subscript[f, 33], Subscript[f, 34]}, {Subscript[f, 42], Subscript[f, 43], Subscript[f, 44]}}] + Det[{{Subscript[f, 11], Subscript[f, 13], Subscript[f, 14]}, {Subscript[f, 31], Subscript[f, 33], Subscript[f, 34]}, {Subscript[f, 41], Subscript[f, 43], Subscript[f, 44]}}] + Det[{{Subscript[f, 11], Subscript[f, 12], Subscript[f, 14]}, {Subscript[f, 21], Subscript[f, 22], Subscript[f, 24]}, {Subscript[f, 41], Subscript[f, 42], Subscript[f, 44]}}] +Det[{{Subscript[f, 11], Subscript[f, 12], Subscript[f, 13]}, {Subscript[f, 21],Subscript[f, 22], Subscript[f, 23]}, {Subscript[f, 31], Subscript[f, 32],Subscript[f, 33]}}]) /. paramFinal /. k -> KK4 /. s ->  SS4 /. μ ->muSS4/. λ -> lambdaSS4

Subscript[b, 4] = Det[{{Subscript[f, 11], Subscript[f, 12], Subscript[f, 13], Subscript[f, 14]}, {Subscript[f, 21], Subscript[f, 22], Subscript[f, 23], Subscript[f, 24]}, {Subscript[f, 31], Subscript[f, 32], Subscript[f, 33], Subscript[f, 34]}, {Subscript[f, 41], Subscript[f, 42], Subscript[f, 43],Subscript[f, 44]}}] /. paramFinal /. k -> KK4 /. s ->  SS4 /. μ -> muSS4 /. λ ->   lambdaSS4


Then, I compute the eigenvalues.

Solve[λ^4 + Subscript[b, 1] λ^3 + Subscript[b, 2] λ^2 + Subscript[b, 3] λ + Subscript[b, 4] == 0, {λ}


There are four roots, where a pair of eigenvalues are purely imaginary and two others are conjugate with real and complex parts, which shows the existence of limit cycles. My problem is that I can not figure out how to find numerically Lyapunov number to know if cycles are stable or not. Any hints or suggestions ? or another way to know if cycles are stable or not ?

• May I suggest that you use f[n] instead of Subscript[f, n] in the question. It is more compact and occasionally more reliable. Also, consider including a plot of the bifurcation diagram. May 18, 2017 at 20:05

I recommend that you review the Liu criterion for the Hopf bifurcation, since it is the most suitable for detecting limit cycles in systems with three or more equations.

It is more convenient to write your system and the Jacobian matrix as follows:

f1 = (a - δ) k - (h (2 - u) - λ a1);
f2 = s (1 - s) - γ a k;
f3 = ((ρ + θ) - (a - δ)) λ + μ γ a;
f4 = (ρ + θ) μ - μ (1 - 2 s) - θ ψ b s^-1;
F = {f1, f2, f3, f4};
V = {k, s, λ, μ};

J = FullSimplify[D[F, {V}]]