Add SICNM: semi-implicit Continuous Newton method (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 SICNM: semi-implicit Continuous Newton method

Part 2 of 2 for SciML/NonlinearSolve.jl#330, implementing the method of Yu et al., Semi-implicit Continuous Newton Method for Power Flow Analysis, arXiv:2312.02809. Companion PR: SciML/OrdinaryDiffEq.jl (adds the Rodas3d Rosenbrock method the paper constructs for this purpose). As discussed on the issue, the DAE solver lives in OrdinaryDiffEq and the wrapper lives here to avoid a circular dependency.

What it does

SICNM(odealg) solves a SteadyStateProblem or NonlinearProblem 0 = g(y) by integrating the continuous Newton flow written as a DAE:

$$\dot y = z,\qquad 0 = J(y),z + g(y),\qquad z(0) = -J(y_0)^{-1} g(y_0)$$

Along the exact flow g(y(t)) = g(y₀)e^{−t}, so integrating to steady state solves the system; Newton's method is exactly explicit Euler on this flow, so integrating it with an L-stable, stiffly accurate, step-adaptive Rosenbrock method is dramatically more robust (this is the "implicitly stable but explicitly solvable" tradeoff the paper targets — no inner nonlinear iterations, one LU per step).

Implementation notes:

• The Jacobian-vector product J(y)z and the residual g(y) for the DAE right-hand side are computed simultaneously from a single ForwardDiff dual evaluation of g at y + εz — no full Jacobian is materialized in the RHS, and no manual Hessian is needed: the ODE solver's own AD of the extended RHS automatically produces the extended Jacobian [0 I; H(y)⊗z + J(y) J(y)] from the paper (nested duals, works with the default AutoForwardDiff of the Rosenbrock solvers). One dense Jacobian is used only for the consistent initialization of z(0).
• Termination is on the nonlinear residual g(y) (default AbsNormTerminationMode(infnorm), matching the paper's ‖g‖∞ < ε), not on du like DynamicSS, and the standard NonlinearSolveBase termination modes are supported.
• iip and oop problems, Float32, save_idxs, custom callbacks, and odesolve_kwargs are supported, mirroring DynamicSS.
• New hard dependency: ForwardDiff (was already a test dependency); version 2.11.3 → 2.12.0.

The recommended ODE solver is Rodas3d() from the companion OrdinaryDiffEq PR (constructed by the paper specifically to damp all modes of this DAE as fast as possible), but any stiffly accurate mass-matrix Rosenbrock method works (Rodas4(), Rodas5P(), ...). The tests here use Rodas4/Rodas5P so they run against currently-released OrdinaryDiffEq; a follow-up can switch/extend to Rodas3d once released.

Verification highlights (run locally)

• Solves problems where Newton diverges: atan(y)=0 from y₀=3, Powell's badly scaled problem from (0,1), and a near-singular Chandrasekhar H-equation (c = 0.9999, N=32) from a far initial guess — all to ‖g‖∞ ≤ 1e-8.
• With the dev'd companion branch, SICNM(Rodas3d()) was verified end-to-end on all of the above.
• Full test suite passes locally including new test/sicnm.jl (iip/oop, both problem types, termination modes incl. safe/safe-best, deliberate short-tspan failure reporting, Float32, save_idxs).

One behavior note: a continuous Newton flow cannot cross det(J(y)) = 0 manifolds — if the initial guess is separated from every root by such a manifold, the integrator's step size collapses and the solve reports failure (this is a property of the method, shared by every CNM variant, not an implementation limitation; the paper's convergence guarantee is for initial values inside the region of attraction).

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Companion PR (Rodas3d solver): SciML/OrdinaryDiffEq.jl#3931