Make the rearranged Jacobian opt-in

[oameye]

get_harmonic_equations derived the Jacobian of the rearranged system on every call. The rearrangement is a symbolic mass-matrix inversion whose cost grows combinatorially with the number of harmonics, and on larger systems it dominates the runtime of the whole function. A van der Pol ansatz in three harmonics is already prohibitive, and a two-tone Duffing with nonlinear damping produces a 366,698-node expression that LLVM then needs about ten minutes to compile on first use.

Nothing in the stability analysis reads the field. HarmonicSteadyState evaluates the Jacobian implicitly, solving for the derivatives numerically once the parameters have values. That agrees with the rearranged matrix at every steady state, which is the only place the Jacobian is ever evaluated: away from one the two differ by ∂(M⁻¹)/∂u ⋅ f(u), which vanishes when f(u) = 0.

So deriving the matrix becomes something you ask for when you want to look at it:

• jacobian is renamed to explicit_jacobian and now defaults to false
• get_Jacobian(eom) remains available for computing it on demand

Breaking

This is a breaking change on two levels.

The kwarg rename is loud: get_harmonic_equations(deq; jacobian=true) now raises a MethodError rather than silently doing something different. That is deliberate, and it covers everyone who passed the kwarg explicitly.

The default flip is quiet, and worth stating plainly: code that called get_harmonic_equations(deq) with no kwarg and then read harmonic_eq.jacobian used to get a filled symbolic matrix and now gets the NaN placeholder, with no error. There is no shim for that case; the migration is explicit_jacobian=true or get_Jacobian(harmonic_eq).

Note that the placeholder was already the outcome whenever the symbolic solve exceeded the Bareiss budget, so callers reading that field already had to handle it.

Testing

Full suite passes. The four call sites in the test suite that passed jacobian= were updated.