Always evaluate the Jacobian implicitly, and add an AD option

[oameye]

Closes #20.

The problem

Compiling eom.jacobian, the Jacobian of the rearranged system, was preferred whenever that field was filled. The rearrangement inverts the mass matrix symbolically, so the expression grows explosively with the number of harmonics. Staged timings for the system in #20 (two drive tones, η*x^2*d(x,t)):

stage	time
homotopy solve	0.07s
sort	0.003s
classify	0.005s
build the Jacobian	19s
first Jacobian call	583s
subsequent calls	0.6ms

The 583s is LLVM compiling a single 366,698-node expression. Nothing was wrong with the solve; it was all in a matrix that gets used only to check stability.

The fix

The implicit route reaches the same matrix by solving J1 x = -J0 once the parameters have values, and stays at ~1000 nodes on the same system. It is now used unconditionally, so the rearranged Jacobian is never compiled.

The two agree only at a steady state, which is the only place the Jacobian is ever evaluated: away from one they differ by ∂(M⁻¹)/∂u ⋅ f(u), which vanishes when f(u) = 0. Verified across five systems: max |ΔJ| and max eigenvalue difference both ≤ 2.5e-15 at every real steady state.

Node counts, rearranged vs implicit:

system	rearranged	implicit
Duffing	181	69
parametron	195	75
two-tone Duffing	1632	340
one-tone Duffing + η	1746	147
van der Pol, 2 harmonics	174458	661
#20	366698	1036

Issue #20 goes from 626s to 0.23s warm, same 9 branches, 3 real, 2 stable. The van der Pol first call goes from 865s to 0.2s. End-to-end on 100-point sweeps of the small systems, where the rearranged Jacobian was cheap, the difference is in the noise (-4.3% to +2.6%).

Also in here

ldiv!(lu!(J1; check=false), J0) instead of -inv(J1) * J0, a sevenfold cut in per-call cost (3.9µs to 0.56µs). check=false is load-bearing: lu! defaults to check=true and throws on NaN, whereas solution arrays are NaN-padded by pad_solutions and a Jacobian at a padded entry has to come back NaN. There is a regression test for this.

cse=true when building the compiled matrices. Jacobian entries share most of their subexpressions, so this cuts the first call by 2-3x on larger systems and about a tenth off per-call, bit-identical results:

variables	first call	with cse	per call	with cse
10	9.88s	4.62s	66.0µs	59.9µs
12	18.63s	6.74s	102.9µs	90.6µs
14	27.91s	10.75s	145.0µs	125.2µs

This is also what ModelingToolkit does by default when generating a Jacobian.

jacobian_backend, opt-in, differentiates the harmonic equations with any ADTypes backend rather than symbolically:

using DifferentiationInterface, ForwardDiff
get_steady_states(harmonic_eq, swept, fixed; jacobian_backend=AutoForwardDiff())

Defaults to nothing, so the implicit path is unchanged unless asked for. DifferentiationInterface is a dependency; the AD package itself is the caller's to load.

DI backends take real inputs, but the Jacobian must accept ComplexF64. The harmonic equations are polynomial, hence holomorphic, so a perturbation along the real axis already carries the full complex derivative: differentiating x -> [Re f; Im f] at fixed imaginary part and recombining gives the exact complex Jacobian. Agreement with the implicit route is 1.2e-16 including at genuinely complex solutions, so no complex-dual machinery is needed and stock backends work.

Worth knowing before reaching for it: on this workload AD costs far more to set up than it saves per call, and the gap widens with system size (at 12 variables, 77-205s of setup against 4.3s, to save tens of microseconds per call across ~10^3 calls). It is useful at small sizes and for differentiating through a sweep, not as a default. AutoEnzyme segfaults on RuntimeGeneratedFunctions and is not usable here.

A trap worth flagging for review

The first version of the AD path called prepare_jacobian once and reused the preparation object, which is the obvious performance move. It produced 25 stable solutions instead of 26 on a Duffing sweep. The Jacobians agreed to 1.2e-16 and the eigenvalues were solidly negative, so the matrix was never wrong: classify_solutions evaluates the Jacobian under Threads.@threads, and a DI preparation object holds buffers the backend writes into, so sharing one across threads races and silently misclassifies the odd solution.

The shared preparation is gone, and there is a test that runs the classification repeatedly and compares against the symbolic result, since this reproduces intermittently rather than always.

Testing

Full suite passes. New test/steady_states/jacobian.jl covers: the implicit and rearranged Jacobians agreeing at steady states on three systems where the rearranged one is small enough to compile; NaN padding returning NaN rather than throwing; #20 solving and classifying without the rearranged Jacobian; and the AD backend matching, including the threading case above.