Lyapunov exponent and stability of limit cycles

Question

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, 31] = D[lambdadot, k];
Subscript[f, 32] = D[lambdadot, s];
Subscript[f, 33] = D[lambdadot, λ];
Subscript[f, 34] = D[lambdadot, μ];
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 ?

Answer 1

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}]]