Add getStability() for steady-state dynamic stability analysis

[gustavdelius]

Summary

Adds two new experimental features for analysing the dynamic stability of a mizer steady state.

New function: getStability()

Analyses the dynamic stability of the steady state stored in params@initial_n by numerically differentiating the one-step-ahead map G(N) = A(N)⁻¹ S(N) at the fixed point N*. Uses the same project_n_loop() C++ Thomas solver as the regular dynamics (centred finite differences with relative step h).

Returns a named list with:

• eigenvalues: complex vector of eigenvalues of L = dG/dN|_{N*}, sorted by decreasing modulus
• spectral_radius: max|λᵢ|; the steady state is stable iff < 1
• stable: logical
• dominant_period: period of the dominant eigenvalue (in time steps)
• hopf_period: period 2π/|Arg(λ)| of the complex eigenvalue closest to the unit circle — the expected limit-cycle period near a Hopf bifurcation; NULL if no complex eigenvalues exist
• n_active: dimension of the Jacobian

Two analysis modes via include_resource = FALSE (default):

• FALSE: resource treated as a quasi-static fast variable (valid for semichemostat). For each perturbed N, n_pp is set to its analytic steady-state value. Jacobian dimension = number of active fish cells. This is the slow-manifold reduction.
• TRUE: both fish and resource cells perturbed independently; resource stepped with the full resource_semichemostat() exact solver. Jacobian dimension = fish cells + all resource cells. Useful for verifying the quasi-static approximation makes little difference.

Updated: steadyNewton(..., stability = FALSE)

When stability = TRUE, calls getStability() after finding the steady state and attaches the result as attr(params, "stability").

Usage

params <- steadyNewton(NS_params, stability = TRUE)
stab <- attr(params, "stability")
stab$stable          # TRUE/FALSE
stab$spectral_radius # dominant eigenvalue modulus
stab$hopf_period     # period in years near Hopf bifurcation, or NULL

# Or call separately:
stab_red  <- getStability(params)                        # quasi-static resource
stab_full <- getStability(params, include_resource = TRUE) # full coupled system

Implementation notes

The key mathematical point is that the Newton residual Jacobian J_F is not equal to dG/dN − I when A depends on N (through predation/growth rates). An earlier implementation using J_F + I gave a spurious spectral radius of ~147. Direct numerical differentiation of G(N) is required.

Tests

7 new tests in test-steadyNewton.R:

• Well-formed list structure
• Stable = TRUE for the NS model at its steady state
• Eigenvalue/spectral_radius consistency
• stability = TRUE argument attaches result correctly
• hopf_period validity
• include_resource = TRUE structural checks
• Full and reduced analyses agree on spectral radius (tolerance 5%)

All 64 tests pass (57 pre-existing + 7 new).

[gustavdelius]

Pushed an additional commit (01ee290e) that extends this work to steady()/projectToSteady():

Detect and report limit cycles in steady()/projectToSteady()

These functions previously had a single convergence test — whether two states t_per years apart are close — which never fires on a limit cycle, so a run that settled onto a cycle went to t_max and returned an arbitrary-phase snapshot as if it were a steady state.

• Fine sampling: each t_per block is now sub-stepped at a new t_save resolution (default dt), recording a per-species biomass summary. This is numerically bit-identical to the old single-block step (a test asserts max abs diff == 0).
• Cycle detection: from that series a limit cycle is identified via the first autocorrelation peak (period) plus a three-period amplitude-stability test with a no-net-decay condition (rejects decaying spirals). Detection no longer requires the cycle period to be commensurate with t_per.
• Structured reporting: both functions attach attr(result, "convergence") = list(type, converged, distance, years, period, amplitude), type ∈ {steady, cycle, not_converged, extinction} — mirroring steadyNewton()'s "stability" attribute. steady() re-stamps it after setBevertonHolt()/setResource(), which would otherwise drop it.

Verification: full NS model at effort 2 → type="cycle", period 5.1 yr, amplitude 0.436 (matching a direct project() measurement); synthetic unit tests confirm correct period/amplitude and no false positives on decaying/flat series. No regressions in test-steady.R/test-steadyNewton.R/test-getLimitCycleSim.R/test-steadySingleSpecies.R; new test-cycle_detection.R added (slow full-NS cycle test is skip_on_cran()).