SDIRK smooth_est: W⁻¹-filtered error estimate underestimates VDP local error by ~3 orders; silently disabled under NonlinearSolveAlg (isnewton gate)

[ChrisRackauckas-Claude]

Investigation follow-up from the NonlinearSolveAlg benchmarking work (#3819/#3820 series). Two intertwined findings, both empirically verified on VdP μ=1e5, tspan=(0,6.3), rtol=1e-8/atol=1e-11, reference Rodas5P@1e-13.

Finding 1: the smoothed estimate can catastrophically underestimate

TRBDF2() (default smooth_est=true) returns final relative error 4.3e-2 at rtol=1e-8. A per-step audit of the accepted trajectory shows:

• Replaying every accepted step with fully-converged nonlinear solves (NLNewton(κ=1e-8, max_iter=200, always_new=true)) changes results by ≤ 0.042 tolerance units (median 0.0016) → the nonlinear solves and all Newton acceptance heuristics (η·ndz<κ, iter-1 shortcut, θ≈1 branch, relaxation) are sound; this is not a nonlinear-solver problem.
• True local error of accepted steps (vs Rodas5P@1e-13 restarts): median 3.3e3× tol, p90 6.1e5×, max 3.4e6×; 927/1133 accepted steps exceed 100× tol. The W⁻¹-filtered (Hosea–Shampine smoothed) embedded estimate is blind to the dominant error mode here.
• κ-sweep control: tightening NLNewton's κ from 1e-2 to 1e-6 leaves the error at ~2–4e-2 — confirming the acceptance threshold is not the lever.

One-variable isolation (only smooth_est changed):

config	naccept	nreject	final rel err
TRBDF2 smooth_est=true (default)	1,133	443	4.3e-2
TRBDF2 smooth_est=false	45,261	3,839	9.6e-7
KenCarp4 smooth_est=true	498	1	3.2e-2
KenCarp4 smooth_est=false	5,825	674	6.9e-8

The smoothing is the literature-standard stiff filter and clearly intentional — but a 1e-8 request returning 4e-2 with Success is an accuracy trap on relaxation-oscillator problems. (FBDF/QNDF at the same tolerances deliver ~1e-6 — this is SDIRK-estimator-specific.) Whether the default or the filter formulation should change is a design call; filing the quantified evidence.

Finding 2: NonlinearSolveAlg silently loses the smoothing (isnewton gate)

generic_imex_perform_step.jl (~L1255, ~L2299):

if isnewton(nlsolver) && _esdirk_smooth_est(alg)
    est = W \ tmp          # smoothed
else
    est = tmp              # raw
end

isnewton is false for NonlinearSolveAlg, so NSA always error-controls on the raw estimate regardless of smooth_est. Consequences:

• Silent behavioral divergence: the same TRBDF2() request produces a 40× different step count and 4 orders different accuracy depending on the inner nonlinear solver choice, which should be an implementation detail.
• This is the entire remaining TRBDF2/KenCarp4 performance gap measured in the NonlinearSolveAlg: give inner NonlinearFunction W's structure as jac_prototype #3819/NonlinearSolveAlg: reuse stored J on dt-only W updates (do_newJW parity) #3820 benchmark series (NSA at 1.5–2.2× NLNewton wall time from taking 4–6× more steps): NLNewton(smooth_est=false) reproduces NSA's step counts and accuracy nearly exactly (table above).

For NSA parity, the gate needs (isnewton(nlsolver) || nsa_with_W(nlsolver)) && smooth_est, with the smoothing solve done against the NSA-held W (nlsolver.cache.W); the wrinkle is that NSA has no ODE-side linsolve cache — the factorization lives inside the inner NonlinearSolve cache, so either a dedicated linsolve for the estimator or reuse of the inner cache's factorization is needed.

Repro scripts available (per-step audit + isolation experiment) — happy to attach.

🤖 Filed by an agent after local reproduction; benchmarks and audit run on the #3819+#3820+NonlinearSolve#1003 stack, Julia 1.12.6.

[ChrisRackauckas-Claude]

Finding 2 (the silent isnewton gating under NonlinearSolveAlg) is fixed in #3823 — smooth_est now acts identically for NLNewton and NSA (VdP: NSA 1341 steps/3.2e-2 with smoothing vs 45237/6.8e-7 without, matching NLNewton's 1133/4.3e-2 and 45261/9.6e-7; HIRES@1e-8 step counts 715 vs 690). Finding 1 (the smoothed estimator underestimating true local error by ~3 orders on relaxation oscillators — now affecting both solvers equally) remains open as the design question here.