[ChrisRackauckas-Claude]

Draft for a decision by @ChrisRackauckas: please ignore until he has reviewed it. This PR changes default NLNewton behaviour for DIRK methods and needs a work-precision benchmark sweep before it can merge.

What and why

W is assembled as J - I/W_γdt, but the stage right-hand side is scaled by the current 1/γdt (invγdt is reset in initialize!). On an eigenmode λ, the stale iteration matrix W(W_γdt)⁻¹W(γdt) = (λ - 1/γdt)/(λ - 1/W_γdt) therefore runs from 1 (stiff) to 1/r (nonstiff), with r = γdt/W_γdt. The Kennedy–Carpenter factor 2/(1+r) (DIRK review, eq. 90) was derived for I - γhJ with an unscaled residual. Carried over to this formulation, the balanced factor is 2r/(1+r) = 2/(1 + W_γdt/γdt). For unchanged J, r > 0 and real λ ≤ 0, it bounds the contraction of a linear stage by |r-1|/(r+1). The code applied 2/(1 + γdt/W_γdt), which makes a linear stage contract worse than no correction, in both the stiff and the nonstiff limit.

The factor is now stale_W_scale and is used by in-place NLNewton and by NonlinearSolveAlg's rescale_stale_W_rhs! (the copy added in #4028). The DAE branch (W_γdt stores cj = α/γdt) and the coefficient-multistep branch (effective step γdt/α) are deliberately left on the old factor. The independent review below found the DAE orientation already balanced and the multistep one in need of a separate derivation. Out-of-place NLNewton applies no rescale at all, and this PR does not add one.

The ROBER testset in nsa_jacobian_reuse_tests.jl (red on master since the OrdinaryDiffEqCore 4.18.1 dt rounding, 45869 < 44710) compared accepted-step counts. Those counts are chaotic in the tolerances. KenCarp4 with nlsolve = NonlinearSolveAlg(NewtonRaphson()), ROBER, tspan 1e11, accepted steps:

reltol, abstol	NLNewton	master factor	no rescale	this PR
1e-9, 1e-12	35768	45869	51311	47504
1e-8, 1e-12	37114	45980	51776	10579
1e-9, 1e-11	11870	24116	19937	17927
1e-7, 1e-10	12292	13127	12630	12529

Raising the threshold would only hide this noise. The testset is replaced by a test of the property the rescale is supposed to provide. On the linear Prothero–Robinson problem u' = λ(u - cos t) - sin t, ImplicitEuler takes three steps, then dt is set to r·W_γdt so that W is reused (asserted). The Newton iterates are read off the points f is evaluated at and compared against the exact implicit-Euler stage root. The first contraction must equal (r-1)/(r+1)·(1+h₀λ)/(1-h₀λ) (rtol 1e-4), for r ∈ {0.9, 1.1}, λ ∈ {-1e6, -1e-2}, and both NonlinearSolveAlg and NLNewton. The problem is linear and the step sequence is fixed, so nothing in the test depends on step-count chaos. This replaces a chaotic step-count assertion with a stronger, exact property. It is not a loosening: the new test fails both on master and with the rescale removed.

Verification

Measured contraction e₁/e₀ (ImplicitEuler, r = 1.1):

λ	master factor	no rescale	this PR	expected
-1e6	+0.04762	+4.5e-6	-0.047614	-0.047614
-1e-2	+0.1342	+0.0909	+0.047600	+0.047600

New testset with the master factor restored (in-process @eval redefinition of stale_W_scale to 2/(1+γdt/W_γdt), which applies to both the NLNewton and NonlinearSolveAlg paths):

master                                   |   24     8     32  34.7s
  stale-W Newton directions are rescaled |   24     8     32  33.7s
  NonlinearSolveAlg λ=-1.0e6 r=1.1: Evaluated: 0.04762337640709256 ≈ -0.04761428595236919 (rtol=0.0001)

With the rescale removed (stale_W_scale ≡ 1):

off                                      |   24     8     32  33.4s
  NonlinearSolveAlg λ=-1.0e6 r=1.1: Evaluated: 4.545227488283e-6 ≈ -0.04761428595236919 (rtol=0.0001)

All 8 contraction assertions fail in each variant. The W-reuse, ratio and iterate-count assertions pass, so the test cannot pass vacuously when W is rebuilt.

With this PR, the whole nsa_jacobian_reuse_tests.jl file:

Test Summary: | Pass  Total     Time
file          |   50     50  4m54.5s

GROUP=Core for OrdinaryDiffEqNonlinearSolve (Pkg.test, Julia 1.12): Developer API Tests | 27 27, then Newton Tests errors with LineSearchException: Search direction is not a direction of descent (Trapezoid + NLNewton(relax = BackTracking()), LineSearches 7.8.2). That is the pre-existing master failure addressed by #4714. Because it is a top-level @safetestset error, it aborts the remaining files. A rerun of the rest of the group with Newton Tests excluded (local edit, not committed) was still running at push time: unverified.

OrdinaryDiffEqSDIRK GROUP=Core ([sources] pick up this branch's NonlinearSolve; Julia 1.12, --check-bounds=yes):

Tableau consistency |  103    103  39m57.0s
Stage predictors |  124    124  22.5s
Convergence   |  139       2    141  56m11.0s

The DAE tests file was cut off by my local 2-hour timeout and is unverified. The 2 broken results come from existing @test_broken assertions.

typos is clean on both changed files. Runic --check is clean on both changed files.

Validation still needed

From the independent review (Codex gpt-6.1-sol, round 1, CHANGES). These items are required before merging:

1. A representative stiff-ODE work-precision sweep, at matched accuracy, for the default NLNewton path and for NonlinearSolveAlg, comparing the old factor, this PR's factor, and no rescale. DIRK methods: at least TRBDF2, KenCarp4, Kvaerno5, and the ESDIRK family. Problems: the SciMLBenchmarks stiff ODE set (ROBER, HIRES, OREGO, POLLU, Van der Pol, E5). Measure error against a reference vs nf, Jacobian/W builds and wall time.
2. Chris's decision between inverting the factor (this PR) and removing the rescale. Removal is near-exact in the stiff limit.
3. The full GROUP=Core runs of OrdinaryDiffEqNonlinearSolve (blocked here by the pre-existing BackTracking error, fix NLNewton BackTracking when FD residual-norm slope is noisy #4714) and OrdinaryDiffEqSDIRK (DAE tests not completed locally), plus CI.

A small spot check (not a sweep): final-time error vs a Rodas5P reference (reltol 1e-13) for three tolerances (abstol = 1e-3·reltol). old swaps the previous factor in at runtime; new is this branch. All 72 solves returned Success. On ROBER (tspan 1e5) and HIRES, the new factor used fewer f evaluations in all 24 cases (often 5–20% fewer, up to 33%). It gave a lower final error in only 6 of the 24, so most points trade some accuracy for less work and these points cannot settle the question. On Van der Pol (μ = 1e6), the final-point error is phase-dominated and moves both ways. KenCarp4+nsa tracks KenCarp4 (NLNewton) closely. Script and raw output: wp.jl, old, new.

Spot-check table

problem	alg	tol	nf old→new	nreject old→new	final relerr old→new
ROBER	TRBDF2	1e-05	3135 → 2879	84 → 65	1.592e-05 → 1.551e-05
ROBER	TRBDF2	1e-07	10089 → 8097	9 → 73	1.519e-06 → 1.748e-06
ROBER	TRBDF2	1e-09	38355 → 34579	733 → 477	1.476e-07 → 1.509e-07
ROBER	KenCarp4	1e-05	2146 → 1941	0 → 1	3.926e-07 → 4.015e-07
ROBER	KenCarp4	1e-07	5128 → 4217	0 → 0	1.513e-08 → 2.932e-08
ROBER	KenCarp4	1e-09	12493 → 11068	0 → 0	4.836e-10 → 6.714e-10
ROBER	Kvaerno5	1e-05	3667 → 3415	1 → 3	4.104e-06 → 3.291e-06
ROBER	Kvaerno5	1e-07	7977 → 5311	2 → 7	1.666e-08 → 1.417e-08
ROBER	Kvaerno5	1e-09	11822 → 11650	4 → 4	1.106e-09 → 4.511e-09
ROBER	KenCarp4+nsa	1e-05	2149 → 1944	0 → 1	3.926e-07 → 4.015e-07
ROBER	KenCarp4+nsa	1e-07	5115 → 4204	0 → 0	1.513e-08 → 2.932e-08
ROBER	KenCarp4+nsa	1e-09	12481 → 11055	0 → 0	4.836e-10 → 6.714e-10
HIRES	TRBDF2	1e-05	1696 → 1667	37 → 36	1.821e-03 → 1.852e-03
HIRES	TRBDF2	1e-07	4537 → 4031	47 → 47	1.935e-04 → 1.934e-04
HIRES	TRBDF2	1e-09	14001 → 13358	324 → 152	1.399e-05 → 1.404e-05
HIRES	KenCarp4	1e-05	1646 → 1579	0 → 1	9.427e-05 → 3.761e-05
HIRES	KenCarp4	1e-07	4088 → 3874	18 → 18	2.741e-06 → 2.828e-06
HIRES	KenCarp4	1e-09	9961 → 9511	26 → 27	2.558e-08 → 3.081e-08
HIRES	Kvaerno5	1e-05	2314 → 2022	0 → 0	1.806e-05 → 8.770e-05
HIRES	Kvaerno5	1e-07	4544 → 4235	6 → 6	7.681e-07 → 1.317e-06
HIRES	Kvaerno5	1e-09	11175 → 11134	27 → 28	1.318e-08 → 1.617e-08
HIRES	KenCarp4+nsa	1e-05	1635 → 1574	0 → 1	9.427e-05 → 3.761e-05
HIRES	KenCarp4+nsa	1e-07	4065 → 3851	18 → 18	2.741e-06 → 2.828e-06
HIRES	KenCarp4+nsa	1e-09	9953 → 9494	26 → 27	2.558e-08 → 3.081e-08
VdP1e6	TRBDF2	1e-05	12811 → 12669	10 → 14	4.482e-01 → 5.049e-01
VdP1e6	TRBDF2	1e-07	16949 → 16921	186 → 184	1.944e-01 → 2.202e-01
VdP1e6	TRBDF2	1e-09	33738 → 33697	554 → 570	6.474e-03 → 5.210e-02
VdP1e6	KenCarp4	1e-05	19317 → 18862	0 → 0	2.636e-01 → 2.900e-01
VdP1e6	KenCarp4	1e-07	25472 → 26839	0 → 0	4.743e-03 → 2.936e-03
VdP1e6	KenCarp4	1e-09	33285 → 32899	73 → 73	3.014e-02 → 4.531e-02
VdP1e6	Kvaerno5	1e-05	25173 → 25104	0 → 0	1.279e-01 → 1.114e-01
VdP1e6	Kvaerno5	1e-07	32169 → 32947	0 → 0	1.115e-02 → 1.642e-02
VdP1e6	Kvaerno5	1e-09	38673 → 40010	0 → 0	5.454e-03 → 3.146e-03
VdP1e6	KenCarp4+nsa	1e-05	20118 → 19651	0 → 0	2.636e-01 → 2.900e-01
VdP1e6	KenCarp4+nsa	1e-07	26045 → 27454	0 → 0	4.743e-03 → 2.936e-03
VdP1e6	KenCarp4+nsa	1e-09	33665 → 33223	73 → 73	2.779e-02 → 4.532e-02

Not verified

• Work-precision or benchmark behaviour of NLNewton with DIRK methods on nonlinear problems. This changes the default NLNewton for every SDIRK/ESDIRK method and must be benchmarked (SciMLBenchmarks stiff ODE work-precision sweep) before merging. Dropping the rescale entirely is the other option: it is near-exact in the stiff limit, where W reuse matters most.
• The multistep and DAE stale-W factors (left unchanged, not tested).
• The new test fails against the master source only through the in-process redefinition shown above. I did not rerun it on a separate master checkout.
• The independent review: Codex gpt-6.1-sol verification (CONFIRMED). Notes and reproducer: NOTES.md, exp4.jl.

Unrelated, not fixed here: Trapezoid(nlsolve = NonlinearSolveAlg(NewtonRaphson(autodiff = AutoForwardDiff()))) returns Unstable at t = 0 on the stiff Prothero–Robinson problem above, where NLNewton succeeds.

Risk assessment

• Risk: medium
• Blast radius: the default NLNewton iteration for every DIRK/SDIRK/ESDIRK method on in-place problems, plus NonlinearSolveAlg on the same methods, wherever W is reused across a dt change. No public API change. BDF/multistep, DAE and out-of-place paths are unchanged.
• Evidence: a closed-form contraction test fails on the master factor and with no rescale, and passes with this PR. The CI-red ROBER step counts are shown above to be chaotic. Independent derivation by Codex gpt-6.1-sol (linked above) confirms the orientation and scope.
• Independent review: Codex CLI 0.163.0 + gpt-6.1-sol confirmed the diagnosis and the ODE-DIRK-only scope (verification). On the PR diff (round 1, head 107165c) it rated the PR medium, needs human review: CHANGES, blocked on the work-precision validation above; it found no wrong-answer defect, and its docstring finding is fixed in 2d660e2 (review).
• Merge: needs human review. It changes default solver behaviour and needs a benchmark sweep plus Chris's decision between inverting and removing the rescale.

🤖 Generated with Claude Code (model: claude-opus-5-5[1m])

https://claude.ai/code/session_011aK5Xvp9NpdkzKuEwtiMGd

🤖 Automated review section from an AI agent running as Chris Rackauckas's account — not written or reviewed by Chris.

Fix by Claude Code 2.1.325 claude-sonnet-5-5 (head f8d87c0)

Failing check: SciMLSensitivity.jl/Core4/1 / Downstream Tests - Core4. The failure is an Enzyme compile error, AssertionError: mi2 !== nothing in threadsfor_common (Enzyme parallelrules.jl), raised in ensembles.jl "EnsembleAlg = EnsembleThreads". It does not involve this PR's change.

Pre-existing on master at the time: the Downstream Core4 job failed with the identical assertion on master pushes between 2026-10-16 and 2026-10-20, for example run 38971870042. Core4 on master pushes from 2026-11-04 onward (runs 41772749111, 43401095529, 43443284988) is success. This PR had not merged master since 2026-10-16.

Change: a normal merge of upstream/master (66 commits) into the branch. There are no conflicts and no source changes of my own. The PR diff against master is unchanged: lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl and test/nsa_jacobian_reuse_tests.jl.

Validation, run locally on the merged head (Julia 1.13.3):

• GROUP=Core for OrdinaryDiffEqNonlinearSolve with only nsa_jacobian_reuse_tests.jl selected (runtests.jl trimmed temporarily, restored, not committed): nsa_jacobian_reuse | 50 50 44.1s, then Testing OrdinaryDiffEqNonlinearSolve tests passed.
• Full GROUP=Core, same Julia: Newton Tests errors with LineSearchException: Search direction is not a direction of descent (4 passed, 1 errored). A clean upstream/master worktree gives the identical error, so it is not caused by this PR. It aborts the rest of the group, so the remaining Core files were not run in the full-group run.
• Runic and typos: no files changed by this merge commit, so not rerun.

Not verified: the Downstream Core4 job itself (needs CI; not run locally). The OrdinaryDiffEqSDIRK Core group and the work-precision sweep listed under "Validation still needed" above are unchanged and still open.

Next counted review must come from a smart model in a fresh session, not this fix session (risk.md § Independent reviews).

---

🤖 Posted by an AI agent — harness: Claude Code 2.1.325 · model: claude-sonnet-5-5
Conversation: local fixer job SciML__OrdinaryDiffEq.jl__4882 (no session URL)