Add Rodas3d: L-stable stiffly accurate Rosenbrock method for SICNM (NonlinearSolve.jl#330)

[ChrisRackauckas-Claude]

Note

This PR is part of an agent-assisted workflow. Please ignore until it has been reviewed by @ChrisRackauckas.

Add Rodas3d: 4-stage 3rd-order L-stable stiffly accurate Rosenbrock method for steady-state focused DAE integration

Part 1 of 2 for SciML/NonlinearSolve.jl#330 (semi-implicit Continuous Newton method, SICNM). As discussed in that issue, the special DAE solver goes into OrdinaryDiffEq while the nonlinear-solve wrapper goes into SteadyStateDiffEq.jl (companion PR: SciML/SteadyStateDiffEq.jl).

What is Rodas3d

Rodas3d is the 4-stage, 3rd-order (2nd-order embedded), stiffly accurate Rosenbrock method constructed in Yu et al., Semi-implicit Continuous Newton Method for Power Flow Analysis, arXiv:2312.02809. Its distinguishing feature is the damping parameter γ = 0.57281606, chosen (per Hairer & Wanner Table 6.4) at the upper end of the range where R(±∞) = 0, so the method damps both stable and unstable modes as fast as possible ("hyper-stability" in the paper's terminology). It is designed for integrating DAEs to steady state — where damping speed matters more than trajectory accuracy — specifically the continuous-Newton DAE

$$\dot y = z,\qquad 0 = J(y)z + g(y)$$

whose equilibrium is the solution of the nonlinear system $0 = g(y)$. In the paper's ill-conditioned power flow benchmarks (9241pegase, ACTIVSg25k) the Rodas3d-based SICNM beat both Rodas4-based SICNM and implicit continuous Newton by significant margins.

Coefficient provenance and verification

The paper states the order conditions but not the coefficients; the coefficients come from the authors' reference implementation rzyu45/MATPOWER-SICNM coeff_rodas3d.m. Verification performed while developing this PR:

• All order conditions stated in the paper (3rd order, 2nd-order embedded, stiff accuracy, consistency of the stage-4 argument with the embedded solution) verified to machine precision from the raw (α, β, γ) tableau. (Note the paper's β′ᵢ = Σⱼβᵢⱼ excludes the diagonal γ.)
• The stability function computed from the tableau matches the paper's closed-form R(z), with R(±∞) = 0.
• The transformation to the Hairer/OrdinaryDiffEq RodasTableau form (A = αΓ⁻¹, C = diag(1/γ) − Γ⁻¹, transformed b/btilde, c, d) was computed in BigFloat and cross-verified: an independent implementation of the untransformed scheme (6) from the paper and the OrdinaryDiffEq-convention loop produce identical steps to ~1e-16 on random singular-mass-matrix problems.
• Observed convergence order ≈ 3.0–3.16 on a nonautonomous nonlinear stiff ODE (exercising the d coefficients) and ≈ 2.9–3.0 on the index-1 continuous-Newton DAE.
• An interesting property found during testing: γ = 0.57281606 is (to ~1e-9) a root of the 4th-order linear order condition (it is an endpoint of the L-stability γ-interval in Hairer & Wanner Table 6.4), so the method superconverges at order 4 on linear problems (verified: observed order 4.09 on prob_ode_linear, exact Taylor coefficient of R(z) − e^z at z⁴ is ≈ −1.0e-9). The convergence tests assert order 4 on the linear test problems (with an explanatory comment) and order 3 via analyticless_test_convergence on a nonlinear problem.
• Robertson mass-matrix DAE solves and matches a tight-tolerance Rodas4 reference to 1e-4.
• Full SICNM pipeline (extended DAE + termination on ‖g(y)‖) solves classic Newton-divergent problems (atan from y₀=3, Powell badly scaled) through the public ODEProblem + mass-matrix API with this method.

Changes

• OrdinaryDiffEqRosenbrockTableaus: new Rodas3dRodasTableau (+ export), version 2.2.1 → 2.3.0.
• OrdinaryDiffEqRosenbrock: new Rodas3d algorithm via the standard consolidated RodasTableau machinery (tabtype, RodasTableauAlgorithms union, alg_order = 3, step/stage limiter support, export, docstring with references), version 2.4.1 → 2.5.0.
• Tests: convergence tests (order 4 on the linear test problems per the superconvergence note above, order 3 on a nonlinear problem), added to the empty-H DAE dense-output regression test (Wrong interpolation in DAE case for Rosenbrock methods with no H matrix #3631 pattern), the nested-ForwardDiff AD tests, and the QA inference/allocation solver lists.
• Docs: added Rodas3d to the Rosenbrock solver lists (semiimplicit + massmatrixdae pages) and the sublibrary README.

Like Rodas3, the method has no stiff-aware dense output (empty H); it falls back to the guarded Hermite path. A stiff-aware interpolant could be derived as follow-up work, but for the steady-state use case dense output is not on the critical path.

Test results

• ODEDIFFEQ_TEST_GROUP=Core for lib/OrdinaryDiffEqRosenbrock: PASS (local, 178 tests incl. the new ones)
• lib/OrdinaryDiffEqRosenbrockTableaus test suite: only pre-existing failures (see note below)
• Runic formatting: clean on all changed files

🤖 Generated with Claude Code

Note on pre-existing QA failures

While running the QA group locally, two failures were observed that also reproduce identically on unmodified master (aa895b5) and are unrelated to this PR: public API has docstrings (HybridExplicitImplicitRK in OrdinaryDiffEqRosenbrock; 8 tableau names in OrdinaryDiffEqRosenbrockTableaus) and all_explicit_imports_are_public (pre-existing OrdinaryDiffEqDifferentiation internal imports). A separate investigation of those master failures is being run in parallel.

[ChrisRackauckas-Claude]

Companion PR (SICNM wrapper): SciML/SteadyStateDiffEq.jl#142

[ChrisRackauckas-Claude]

CI triage — none of the failures are introduced by this PR:

• tests / *(julia lts) failures (AlgConvergence_I/III, Integrators_I, InterfaceI, Regression_I/II): all share one cause — the in-repo lib/DiffEqDevTools fails to precompile on Julia 1.10 with UndefVarError: RungeKuttaMethod not defined before any tests run. Julia-1 lanes of the same groups pass. This hits every PR's LTS lane right now (recently merged PRs were green here a week ago, so a recent dep release is suspected); a separate investigation/fix is in progress.
• sublibrary-ci [QA] failures (Rosenbrock, RosenbrockTableaus, DelayDiffEq): pre-existing on unmodified master (verified locally at aa895b5) — undocumented HybridExplicitImplicitRK/tableau exports and non-public OrdinaryDiffEqDifferentiation imports, latent until SciMLTesting 2.1's new checks (Update SciMLTesting QA docs coverage #3867). A fix PR restoring the public declarations (removed in Document public developer APIs #3832) and adding the missing docstrings is being prepared separately.
• downgrade-sublibraries failures (BDF, Core, Differentiation, FIRK, NonlinearSolve, SDIRK, StochasticDiffEq): identical list failed on merged PR Fix reinit! crash on null-u0 integrators (uType Nothing) #3921; none of these sublibraries are touched here.
• Rosenbrock [GPU]: also failed on merged PR Fix reinit! crash on null-u0 integrators (uType Nothing) #3921; not related to this change.

Locally, ODEDIFFEQ_TEST_GROUP=Core for lib/OrdinaryDiffEqRosenbrock passes with the new tests, and the QA group fails only with the pre-existing master failures above.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Rebased onto current master (38f01c5) — clean rebase, no conflicts. Since #3935 rewrote the RosenbrockCache stage-accumulation kernels, the full ODEDIFFEQ_TEST_GROUP=Core suite was re-run locally on the rebased branch: all tests pass, and Rodas3d's observed convergence orders and Robertson-DAE results are unchanged. The LTS CI lanes should now be green here too since #3939/#3941 are in the base.

[ChrisRackauckas-Claude]

Post-rebase CI triage update — the rebase cleared all LTS lanes (thanks to #3939/#3941). The two checks that newly appeared were investigated and neither is introduced by this PR:

• tests / QA (julia 1) ("public API is rendered in docs", 33 umbrella re-exports): reproduces identically on clean master 38f01c5f. Root cause: SciMLTesting v2.3.0 (released 2026-07-20) flipped run_api_docs(rendered=...) on by default. Already fixed on master by Document umbrella re-exported API #3970 (merged a few hours after this PR's CI ran, adding docs/src/api/common_interface.md), so this lane clears on the next CI run against updated master.
• tests / Integrators_II (julia 1) (resize! with SplitODEProblem, @test_nowarn step! seeing a dt-below-epsilon warning): passes locally on this branch and on clean master; the same group also failed on the previously merged PR Fix reinit! crash on null-u0 integrators (uType Nothing) #3921. The test steps the integrator after resize!(integrator, 2), i.e. with an uninitialized state entry, so whether the step warns depends on memory garbage — inherently flaky, and this PR touches nothing in the Tsit5/SplitODEProblem path.

Remaining failures are the previously-triaged set: Rosenbrock/Tableaus [QA] (fixed by #3940), and the chronic downgrade/GPU/DelayDiffEq-Default-QA failures that also fail on merged PRs.

🤖 Generated with Claude Code