Stiff-aware nested solver for multi-trait adaptive dynamics (canonical equation + differentiable demographic equilibrium)

[dfalster]

Context

Multi-species / multi-trait adaptive dynamics can be very stiff — there is a
strong timescale separation between trait evolution, demography, and physiology.
The current explicit RK solver in odelia degrades badly on stiff systems (step
size collapses for stability, not accuracy). With AD now providing exact
Jacobians, we can build a stiff-aware, differentiable solver stack for
eco-evolutionary dynamics.

This issue captures the intended architecture for multi-trait attractor solving.
It depends on odelia gaining an implicit stepper + steady-state/sensitivity
support: traitecoevo/odelia#35 (cross-package, blocked).

Proposed architecture: three nested odelia systems

Level	System	Timescale	Solver
Outer	Trait dynamics — the canonical equation dx/dτ = M(x)·g(x)	slow	RODAS4(3) (implicit, adaptive)
Middle	Demographic model, solved to fixed-point equilibrium f(n, x) = 0	fast	Newton (reusing RODAS Jacobian); RODAS4(3) as fallback / stability check
Inner	plant physiology	fastest	explicit RK45 (unchanged)

Scope assumption: the demographic attractor is a fixed point. Limit-cycle /
non-equilibrium attractors (needing monodromy-matrix / periodic-orbit
sensitivity) are explicitly out of scope for now.

Key AD design — differentiate the equilibrium implicitly

The selection gradient g(x) = ∂s/∂x_m (invasion fitness s of a rare mutant,
evaluated at mutant = resident in the environment set by the resident
equilibrium) is the RHS of the canonical equation. To get it we need the
derivative of the demographic fixed point w.r.t. the trait.

Do not tape through the Newton iteration. Instead give the fixed-point solve
a hand-written derivative via the implicit function theorem:

dn*/dx = − f_n⁻¹ · f_x        (reuse the LU factorization of f_n)

registered as a custom derivative (XAD supports custom derivatives /
checkpointing). Payoffs:

• Tape stays bounded regardless of how many stiff Newton/RODAS steps reaching
  equilibrium took — it never records the iteration.
• Exact, not iteration-count-dependent.

This is what makes the fixed-point case tractable: the whole environmental
feedback x → n*(x) → E(x) reduces to one linear solve we already have factored.

AD type stack

• T0 = double — inner physiology.
• T1 = active<double> — first order: f_n, f_x, ∂s/∂x_m, and the IFT
  sensitivity dn*/dx.
• T2 = active<active<double>> — second order: the outer canonical-equation
  Jacobian ∂g/∂x (needed for RODAS on the trait dynamics) and singular-
  strategy classification we want anyway — ESS (∂²s/∂x_m² < 0) and convergence
  stability (∂g/∂x < 0).

Prefer forward-over-forward (tangent-over-tangent): trait dimension is
typically small relative to demographic state, and it avoids a live inner reverse
tape crossing the odelia↔consumer DLL boundary — keeping closest to odelia's
single-active_tape_ design. Each of T1, T2 is a distinct Tape<…>
requiring explicit instantiation in odelia (src/Tape.cpp) — tracked in
odelia#35.

Fit with the regnans roadmap

Extends current capabilities (1D selection gradients, viable bounds,
community_demography equilibria, 1-species/1-trait attractor) toward the
envisioned 2D/3D multi-trait attractor solving — with a solver that stays
stable and efficient under stiffness, and gives singular-point classification
(CSS/ESS) as a by-product.

Tasks

• Depends on odelia#35 (implicit stepper + steady-state/IFT sensitivity)
• Custom-derivative registration for the demographic fixed-point solve (IFT)
• Selection gradient g(x) at the resident equilibrium environment via AD
• Outer canonical-equation integration with RODAS4(3) + its ∂g/∂x Jacobian
• Singular-strategy classification (CSS / ESS) from T2 derivatives
• Validate against existing 1D selection-gradient / attractor results