Add NordsieckBDF and DNordsieckBDF: Newton-based Nordsieck BDF solvers

[ChrisRackauckas-Claude]

Please ignore this PR until reviewed by @ChrisRackauckas. Opened as a draft.

Adds a variable-order variable-step BDF pair built on a propagated Nordsieck
history array zn[j] = h^j/j! * y^(j)(t_n), ported from SUNDIALS CVODE:

• NordsieckBDF — ODE / mass-matrix form
• DNordsieckBDF — fully implicit DAE form, f(du, u, p, t) = 0

Both support in-place and out-of-place problems.

Why

FBDF stores raw (t_i, u_i) history and reconstructs the predictor (Lagrange)
and the error/order estimates (finite-difference stencils) every step. A residual
left in that stored history by the nonlinear solve gets amplified by the
extrapolating predictor, which inflates EEst and shrinks the step — so FBDF has
to over-solve the corrector to keep its own step size. Measured on the stiff suite,
FBDF runs ~2.5 Newton iterations/step against CVODE's ~1.5, and that accounts for
essentially the whole FBDF↔CVODE gap.

The Nordsieck array propagates rather than reconstructs: predicting is a
Pascal-triangle shift, a step-size change is zn[j] *= eta^j, and accepting is a
rank-1 update. Nothing is rebuilt from stored points, so a loose corrector cannot
collapse the step size.

I arrived at this after implementing and benchmarking seven localized fixes to
FBDF (cheaper nonlinear tolerance, EEst-driven step sizing, Gustafsson
predictive control, a true predictor-polynomial FLC corrector, Jacobian/W reuse
policy, CVODE's crate convergence test, CVODE's never-shrink step policy) — every
one was neutral or negative, because each policy works and the history
representation then cancels it.

Results

Matched accuracy over ROBER/HIRES/OREGO/POLLU/VDPOL/BRUSS1D (min-of-5 timing, L2
error vs Rodas5P at 1e-14, log-log interpolation onto common accuracy targets):

vs FBDF	time	nf	nW	steps
NordsieckBDF	0.65x	0.98x	0.68x	1.19x
CVODE_BDF	0.52x	0.60x	0.82x	1.31x

So ~35% faster than FBDF, with a third fewer W factorizations. It does not yet
reach CVODE — the residual gap is the nonlinear convergence test's iteration
economy (the package's Hairer–Wanner η·‖dz‖ < κ test with max_iter = 10
against CVODE's crate test with MAXCOR = 3), which is shared infrastructure
rather than anything about the history representation. I left the default
NLNewton() alone rather than ship a tuned default I can't fully justify.

On DAEs at a fixed tolerance, DNordsieckBDF vs DFBDF: index-1 test 92 steps /
116 nf vs 150 / 174; Robertson DAE 459 / 962 vs 616 / 1408.

Relation to the existing JVODE_BDF

OrdinaryDiffEqNordsieck already has JVODE_BDF, but it uses functional
iteration (nlsolve_functional!), so it is not usable on stiff problems — on
ROBER it hits MaxIters at 2e5 steps, on HIRES it returns 55% error, and on OREGO
it throws StackOverflowError. It is also marked experimental and still carries a
# TODO: BDF in its order-increase path. These are the first Newton-based
Nordsieck BDF solvers in the package. Happy to move them into
OrdinaryDiffEqNordsieck instead if you'd prefer them to live alongside JVODE.

Implementation notes

• Reuses the package nonlinear solver rather than rolling its own. ODE:
  COEFFICIENT_MULTISTEP with gamma = 1/l1. DAE: correction relative to
  cache.u₀, giving cj = l1/h — the same role IDA's cj plays. Both give
  gamma*W = dt/l1, so the existing W and Jacobian-reuse machinery is untouched.
• NLNewton(κ = …) acts as CVODE's NLSCOEF (the fraction of the local error
  budget the corrector may consume), because the increment norm is scaled by the
  test quantity tq[2] through the existing has_special_newton_error hook.
• Dense output is the Nordsieck polynomial itself, so it is free; fsallast is
  recovered algebraically with no extra f evaluation.
• The diff is additive: two new files plus 5 lines in OrdinaryDiffEqBDF.jl and
  runtests.jl. No existing solver's code path is modified.

Testing

New nordsieck_tests.jl (9 testsets, all passing locally against master
d8cd163): adaptive accuracy and tolerance scaling for both iip and oop,
max_order 1–5, five stiff problems, dense output, mass matrices, callback
discontinuity restart, index-1 DAEs (iip + oop), Robertson DAE, and the
loose-corrector step-count property. Existing bdf_regression_tests and
dae_convergence_tests pass unchanged. Runic clean.

Correctness was additionally checked by verifying the Nordsieck invariant
zn[j] = h^j/j! * y^(j) holds on every column against an analytic solution.

Two notes for reviewers, neither introduced here: prob_ode_filament throws a
DimensionMismatch under both FBDF and NordsieckBDF on current master; and
FBDF returns Unstable on prob_ode_vanderpol_stiff at abstol=1e-12/reltol=1e-10.

🤖 Generated with Claude Code

https://claude.ai/code/session_016bHrkkvJzfYPMMJMRJ5GqX

[ChrisRackauckas-Claude]

NonlinearSolveAlg backend

Checked both solvers against the pluggable NonlinearSolve.jl corrector; added a
testset for the cases that work (9ce7cd1).

ODE — works, identical results. NordsieckBDF(nlsolve = NonlinearSolveAlg())
matches the built-in Newton to the digit on both paths (linear oop and 2Dlinear
iip: same step count, same error), and tracks it on stiff problems:

problem	NLNewton steps/nf	NonlinearSolveAlg steps/nf
ROBER	349 / 700	369 / 793
HIRES	312 / 900	315 / 824
OREGO	770 / 2316	853 / 2273

The spread is the same shape as FBDF vs FBDF(nlsolve = NonlinearSolveAlg()).

Worth noting: has_special_newton_error is applied in the NonlinearSolveAlg
compute_step! methods as well as the NLNewton ones, so the tq[2] scaling
holds on both backends and κ keeps its NLSCOEF meaning either way.

DAE — does not work, but not because of this PR. NonlinearSolveAlg currently
fails on fully implicit DAEs for the existing DFBDF with exactly the same
errors it gives for DNordsieckBDF:

• in-place: type DAEFunction has no field nlstep_data
• out-of-place: MethodError: no method matching _compute_rhs(::Vector{Float64}, …)

So it looks like a pre-existing gap in the DAE ↔ NonlinearSolveAlg plumbing rather
than anything specific to the new solver — DNordsieckBDF fails identically to
DFBDF, which is the consistent behavior. I have not tried to fix that here; happy
to open it as a separate issue if useful.

[ChrisRackauckas-Claude]

Follow-ups: STALD, failure-driven order reduction, precompile (0d74c18)

Plus one negative result worth recording.

Crate-based convergence test — implemented, measured, reverted. My hypothesis
was that the remaining gap to CVODE was the corrector's stopping rule (the package's
Hairer–Wanner η·‖δ‖ < κ with η = θ/(1-θ) versus CVODE's min(1,crate)·‖δ‖ < tol).
I wired it in as an opt-in trait and A/B'd it at identical κ/max_iter:

κ=1/10, max_iter=10	time	nf	steps
Hairer–Wanner	1.15x	1.60x	2.52x
crate	1.14x	1.60x	2.52x

Identical, and identical at the default κ too. The two rules coincide in practice
because θ is small on converging steps, so θ/(1-θ) ≈ min(1,θ). I reverted it rather
than land an unused trait and an extra branch in a shared hot loop. So the residual
gap to CVODE is not the convergence test — I no longer have a specific explanation
for it, and would rather say that than keep the earlier guess in the PR description.

STALD (stald = true, default false as in CVODE). Wires the existing
StabilityLimitDetectionState with cvBDFStab's quantities; a violation caps the
newly chosen order and forbids growth.

Caveat I'd rather state than bury: I could not get the detector to fire — not on
the stiff suite, not on purely imaginary spectra (a rotation block at ω = 5/20/50),
and it returns false on synthetic geometrically-growing ratio data 200/200 times.
FBDF's STALD behaves identically (stald=true and stald=false give bit-identical
step counts on the same problems). Since my wiring mirrors FBDF's, this looks like a
property of the shared stald.jl rather than of either caller — possibly worth a
separate look at whether #3088's detector can trigger at all. Either way the path is
unexercised here, so treat it as untested infrastructure rather than a working feature.

Order reduction on repeated corrector failures. post_newton_controller! now
drops the order after three consecutive failures instead of only shrinking dt
(matching FBDF). On a discontinuous-RHS test: 446 steps / 4 convergence failures
versus FBDF's 670 / 8.

Precompile. NordsieckBDF added to the workload.

Cleanup. Removed an etamax ternary whose two branches were identical — harmless
today since CVODE's ETA_MAX_ES and ETA_MAX_GS are both 10, but it read as a bug.

Tests are now 11 testsets, all passing against master. Runic clean.

[ChrisRackauckas-Claude]

Instrumenting the remaining gap to CVODE — and a correction

Correction first. In my previous comment I said the crate-based convergence test
made "no measurable difference". That was wrong: I compared coarse matched-accuracy
aggregates from two separate runs. Re-running it with a counter on the branch to
confirm it actually fires, it does reduce Newton iterations — HIRES 2.56 → 2.33
per step (~9%). It just doesn't convert into time: at matched accuracy it gives
0.66x versus 0.65x without. The saved iterations come back as extra steps. Same
conclusion (reverted, not worth a shared-infrastructure change), corrected reason.

Per-step wall time (ns per accepted+rejected step, abstol 1e-8 / reltol 1e-6).
"alone" is my standalone CVODE port running the identical algorithm outside the
package, which is the useful control:

problem	N	FBDF	NordsieckBDF	alone	CVODE	pkg/alone	pkg/CVODE
ROBER	3	2938	1643	1521	1349	1.08x	1.22x
HIRES	8	3963	2592	1893	1637	1.37x	1.58x
OREGO	3	3472	2040	1367	1257	1.49x	1.62x
POLLU	20	6180	4103	3500	2648	1.17x	1.55x
VDPOL	2	2950	1702	1299	1061	1.31x	1.60x

So per step this solver is already 1.5–1.8x cheaper than FBDF, and the residual
gap to CVODE splits into two roughly equal halves: ~1.1–1.5x of package machinery
(the same algorithm is that much faster outside the package), and ~1.0–1.3x of
Julia-vs-C on top of that.

Where the time goes (profile, self-time). getrs! — the LAPACK triangular solve
inside the Newton iteration — is the single dominant cost for both solvers: 38.8%
on HIRES and 33.7% on OREGO for NordsieckBDF, against 24.5%/23.8% for FBDF. The share
is higher precisely because everything else got cheaper; absolute LAPACK time per step
is about the same. After that it is a long tail with no second hotspot.

Components of the package-vs-standalone gap, measured rather than guessed:

• LinearSolve.jl per-call overhead: ~30–40 ns for N = 2…63 (raw ldiv! vs a
  primed LinearProblem cache). At ~2.5 solves/step that is well under a third of
  the gap.
• FunctionWrappers / specialization: ~3% (FullSpecialize versus the default).
• The remainder is diffuse integrator machinery — residual weighting, stats,
  controller dispatch, buffer copies — and it is shared with FBDF, not specific to
  this solver.

Conclusion. There is no single remaining hotspot to fix in NordsieckBDF. The one
real lever is Newton iterations per step (CVODE 1.5, this 2.3–2.8), and it is already
at its optimum: loosening κ trades iterations for steps at roughly 1:1 in wall time,
with the best measured setting NLNewton(κ = 1//30, max_iter = 4) at 0.61–0.64x versus
the default's 0.65x — not enough to justify overriding the conventional default on six
problems. Getting closer to CVODE would mean reducing the shared per-step machinery
cost, which would speed up FBDF equally and is a separate piece of work.