# 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