Add Rodas3d: L-stable stiffly accurate Rosenbrock method for steady-state DAE integration (#3931)

ChrisRackauckas-Claude · ChrisRackauckas · claude · web-flow · commit 754b504c71d4 · 2026-07-23T09:43:12.000Z
Adds the 4-stage 3rd-order (2nd-order embedded) stiffly accurate Rosenbrock method from Yu et al., arXiv:2312.02809, constructed for the semi-implicit Continuous Newton method (SciML/NonlinearSolve.jl#330). The damping parameter γ = 0.57281606 gives R(±∞) = 0 so both stable and unstable modes are damped, and is a root of the 4th-order linear order condition, so the method is additionally 4th order on linear problems. Coefficients come from the authors' reference implementation (github.com/rzyu45/MATPOWER-SICNM) and were verified against all order conditions stated in the paper; the transformation to the RodasTableau convention was computed in BigFloat and cross-checked against an independent implementation of the paper's untransformed scheme. The companion SciML/SteadyStateDiffEq.jl PR adds the SICNM nonlinear solver that uses this method. Claude-Session: https://claude.ai/code/session_01NNigabD3sZaUKvDZ1NVQgv Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
diff --git a/docs/src/massmatrixdae/Rosenbrock.md b/docs/src/massmatrixdae/Rosenbrock.md
@@ -129,6 +129,7 @@ Rosenbrock23
 Rosenbrock32
 ROS3P
 Rodas3
+Rodas3d
 Rodas23W
 Rodas3P
 Rodas4
diff --git a/docs/src/semiimplicit/Rosenbrock.md b/docs/src/semiimplicit/Rosenbrock.md
@@ -72,6 +72,7 @@ Rosenbrock23
 Rosenbrock32
 ROS3P
 Rodas3
+Rodas3d
 Rodas23W
 Rodas3P
 Rodas4
diff --git a/lib/OrdinaryDiffEqRosenbrock/Project.toml b/lib/OrdinaryDiffEqRosenbrock/Project.toml
@@ -1,6 +1,6 @@
 name = "OrdinaryDiffEqRosenbrock"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "2.4.2"
+version = "2.5.0"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
 
 [deps]
diff --git a/lib/OrdinaryDiffEqRosenbrock/README.md b/lib/OrdinaryDiffEqRosenbrock/README.md
@@ -15,6 +15,7 @@ While completely independent and usable on its own, users wanting the full ODE s
 - `Rosenbrock32`
 - `RosShamp4`
 - `Rodas3`
+- `Rodas3d`
 - `Rodas4`
 - `Rodas4P`
 - `Rodas5`
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl
@@ -40,7 +40,7 @@ using OrdinaryDiffEqRosenbrockTableaus: OrdinaryDiffEqRosenbrockTableaus,
     ROS2SRodasTableau, ROS34PRwRodasTableau, ROS34PW1aRodasTableau,
     ROS34PW1bRodasTableau, ROS34PW2RodasTableau, ROS3PRL2RodasTableau,
     ROS3PRLRodasTableau, ROS3PRRodasTableau, ROS3PRodasTableau, ROS3RodasTableau,
-    Rodas3PRodasTableau, Rodas3RodasTableau, Rodas42Tableau, Rodas4P2Tableau,
+    Rodas3PRodasTableau, Rodas3RodasTableau, Rodas3dRodasTableau, Rodas42Tableau, Rodas4P2Tableau,
     Rodas4PTableau, Rodas4PWTableau, Rodas4Tableau, Rodas5Tableau, RodasTableau,
     Ros4LStabRodasTableau, RosShamp4RodasTableau, RosenbrockW6S4OSRodasTableau,
     Scholz4_7RodasTableau, Veldd4RodasTableau, Velds4RodasTableau
@@ -209,7 +209,7 @@ PrecompileTools.@compile_workload begin
 end
 
 export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
-    Ros4LStab, ROS3P, Rodas3, Rodas23W, Rodas3P, Rodas4, Rodas42, Rodas4P, Rodas4P2,
+    Ros4LStab, ROS3P, Rodas3, Rodas3d, Rodas23W, Rodas3P, Rodas4, Rodas42, Rodas4P, Rodas4P2,
     Rodas4PW, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P, HybridExplicitImplicitRK,
     Tsit5DA, RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3, ROS34PRw,
     ROS3PRL, ROS3PRL2, ROK4a,
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/alg_utils.jl b/lib/OrdinaryDiffEqRosenbrock/src/alg_utils.jl
@@ -6,6 +6,7 @@ alg_order(alg::ROS3PRL) = 3
 alg_order(alg::ROS3PRL2) = 3
 alg_order(alg::ROS3P) = 3
 alg_order(alg::Rodas3) = 3
+alg_order(alg::Rodas3d) = 3
 alg_order(alg::Rodas23W) = 3
 alg_order(alg::ROS2) = 2
 alg_order(alg::ROS2PR) = 2
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl b/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl
@@ -26,6 +26,12 @@ for (Alg, desc, refs, is_W) in [
             "- Sandu, Verwer, Van Loon, Carmichael, Potra, Dabdub, Seinfeld, Benchmarking stiff ode solvers for atmospheric chemistry problems-I. \n  implicit vs explicit, Atmospheric Environment, 31(19), 3151-3166, 1997.",
             false,
         ),
+        (
+            :Rodas3d,
+            "3rd order L-stable stiffly accurate Rosenbrock method with a 2nd order embedded method,\nconstructed with damping parameter γ = 0.57281606 so that R(±∞) = 0 (damps both stable and\nunstable modes). Since this γ is a root of the 4th-order linear order condition, the method\nis 4th order on linear problems. Designed for steady-state focused DAE integration such as\nthe semi-implicit continuous Newton method, where fast damping toward the equilibrium\nmatters more than trajectory accuracy.",
+            "- Yu, R., Gu, W., Xu, Y., Lu, S. (2024). Semi-implicit Continuous Newton Method\n  for Power Flow Analysis. arXiv:2312.02809. https://arxiv.org/abs/2312.02809",
+            false,
+        ),
         (
             :Rodas23W,
             "An Order 2/3 L-Stable Rosenbrock-W method for stiff ODEs and DAEs in mass matrix form. 2nd order stiff-aware interpolation and additional error test for interpolation.",
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
@@ -408,7 +408,7 @@ end
 _get_step_limiter(alg) = trivial_limiter!
 _get_stage_limiter(alg) = trivial_limiter!
 for Alg in (
-        :Rosenbrock23, :Rosenbrock32, :ROS3P, :Rodas3, :Rodas23W, :Rodas3P,
+        :Rosenbrock23, :Rosenbrock32, :ROS3P, :Rodas3, :Rodas3d, :Rodas23W, :Rodas3P,
         :Rodas4, :Rodas42, :Rodas4P, :Rodas4P2, :Rodas4PW, :Rodas5,
         :Rodas5P, :Rodas5Pe, :Rodas5Pr, :Rodas6P,
     )
@@ -431,6 +431,7 @@ tabtype(::Rodas6P) = Rodas6PTableau
 # Consolidated methods: tableau type dispatch
 tabtype(::ROS3P) = ROS3PRodasTableau
 tabtype(::Rodas3) = Rodas3RodasTableau
+tabtype(::Rodas3d) = Rodas3dRodasTableau
 tabtype(::Rodas3P) = Rodas3PRodasTableau
 tabtype(::Rodas23W) = Rodas23WRodasTableau
 tabtype(::ROS2) = ROS2RodasTableau
@@ -459,7 +460,7 @@ tabtype(::RosenbrockW6S4OS) = RosenbrockW6S4OSRodasTableau
 const RodasTableauAlgorithms = Union{
     Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas4PW,
     Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P,
-    ROS3P, Rodas3, Rodas3P, Rodas23W,
+    ROS3P, Rodas3, Rodas3d, Rodas3P, Rodas23W,
     ROS2, ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7,
     ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3,
     ROS34PRw, ROS3PRL, ROS3PRL2, ROK4a,
diff --git a/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl b/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl
@@ -159,6 +159,38 @@ end
     end
 
 
+    println("Rodas3d")
+
+    # Rodas3d's damping parameter γ = 0.57281606 is a root of the 4th-order
+    # linear order condition (an endpoint of the L-stability interval in
+    # Hairer & Wanner Table 6.4), so the method superconverges at order 4 on
+    # linear problems; the design order 3 is checked on a nonlinear problem below.
+    prob = prob_ode_linear
+
+    sim = test_convergence(dts, prob, Rodas3d())
+    @test sim.𝒪est[:final] ≈ 4 atol = testTol
+
+    sol = solve(prob, Rodas3d())
+    @test length(sol.t) < 20
+    @test SciMLBase.successful_retcode(sol)
+
+    prob = prob_ode_2Dlinear
+
+    sim = test_convergence(dts, prob, Rodas3d())
+    @test sim.𝒪est[:final] ≈ 4 atol = testTol
+
+    sol = solve(prob, Rodas3d())
+    @test length(sol.t) < 20
+    @test SciMLBase.successful_retcode(sol)
+
+    prob = ODEProblem(
+        (u, p, t) -> [-2u[1] + u[2]^2, -u[2] + sin(u[1]) + 0.1t],
+        [1.0, 0.5], (0.0, 1.0)
+    )
+    test_setup = Dict(:alg => Rodas4(), :reltol => 1.0e-13, :abstol => 1.0e-13)
+    sim = analyticless_test_convergence((1 / 2) .^ (7:-1:4), prob, Rodas3d(), test_setup)
+    @test sim.𝒪est[:final] ≈ 3 atol = testTol
+
     println("Rodas5P")
 
     prob = prob_ode_linear
@@ -418,7 +450,7 @@ end
     # interpolation errors orders of magnitude larger than the knot errors.
     for (method, expected_interp) in (
             (Rodas4P, 1.0e-10), (ROS34PW2, 5.0e-2),
-            (ROS34PW3, 5.0e-2), (Rodas3, 5.0e-2),
+            (ROS34PW3, 5.0e-2), (Rodas3, 5.0e-2), (Rodas3d, 5.0e-2),
         )
         sol = solve(prob, method(); dense = true, reltol = 1.0e-4, abstol = 1.0e-4)
         err_knot = maximum(abs.(sol[1, :] .- anasol(sol.t, p)))
diff --git a/lib/OrdinaryDiffEqRosenbrock/test/qa/Project.toml b/lib/OrdinaryDiffEqRosenbrock/test/qa/Project.toml
@@ -2,14 +2,15 @@
 AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[sources.OrdinaryDiffEqCore]
-path = "../../../OrdinaryDiffEqCore"
+[sources]
+OrdinaryDiffEqCore = {path = "../../../OrdinaryDiffEqCore"}
 
 [compat]
 AllocCheck = "0.2"
diff --git a/lib/OrdinaryDiffEqRosenbrock/test/qa/allocation_tests.jl b/lib/OrdinaryDiffEqRosenbrock/test/qa/allocation_tests.jl
@@ -15,7 +15,7 @@ using Test
 
     rosenbrock_solvers = [
         Rosenbrock23(), Rosenbrock32(), RosShamp4(), Veldd4(), Velds4(), GRK4T(), GRK4A(),
-        Rodas3(), Rodas23W(), Rodas3P(), Rodas4(), Rodas42(), Rodas4P(), Rodas4P2(), Rodas5(),
+        Rodas3(), Rodas3d(), Rodas23W(), Rodas3P(), Rodas4(), Rodas42(), Rodas4P(), Rodas4P2(), Rodas5(),
         Rodas5P(), Rodas5Pe(), Rodas5Pr(), Rodas6P(),
     ]
 
diff --git a/lib/OrdinaryDiffEqRosenbrock/test/qa/inference_tests.jl b/lib/OrdinaryDiffEqRosenbrock/test/qa/inference_tests.jl
@@ -22,7 +22,7 @@ using Test
 @testset "Rosenbrock perform_step! Inference Tests" begin
     rosenbrock_solvers = [
         Rosenbrock23(), Rosenbrock32(), RosShamp4(), Veldd4(), Velds4(), GRK4T(), GRK4A(),
-        Rodas3(), Rodas23W(), Rodas3P(), Rodas4(), Rodas42(), Rodas4P(), Rodas4P2(), Rodas5(),
+        Rodas3(), Rodas3d(), Rodas23W(), Rodas3P(), Rodas4(), Rodas42(), Rodas4P(), Rodas4P2(), Rodas5(),
         Rodas5P(), Rodas5Pe(), Rodas5Pr(), Rodas6P(),
     ]
 
diff --git a/lib/OrdinaryDiffEqRosenbrock/test/rosenbrock_ad_tests.jl b/lib/OrdinaryDiffEqRosenbrock/test/rosenbrock_ad_tests.jl
@@ -71,6 +71,7 @@ const ROSENBROCK_ALGS_TO_TEST = (
     # several of them hit dtmin on this simple non-stiff problem under
     # nested ForwardDiff and that failure mode is not what we're testing.
     Rodas3(),
+    Rodas3d(),
     Rodas23W(),
     ROS3P(),
     Rodas4(),
diff --git a/lib/OrdinaryDiffEqRosenbrockTableaus/Project.toml b/lib/OrdinaryDiffEqRosenbrockTableaus/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqRosenbrockTableaus"
 uuid = "b4bd8bb3-f80f-41d2-9b21-73a655b304b9"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "2.2.2"
+version = "2.3.0"
 
 [extras]
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
diff --git a/lib/OrdinaryDiffEqRosenbrockTableaus/src/OrdinaryDiffEqRosenbrockTableaus.jl b/lib/OrdinaryDiffEqRosenbrockTableaus/src/OrdinaryDiffEqRosenbrockTableaus.jl
@@ -4,7 +4,7 @@ include("rosenbrock_tableaus.jl")
 
 export RodasTableau, Rodas5Tableau,
     Rodas4Tableau, Rodas42Tableau, Rodas4PTableau, Rodas4P2Tableau,
-    ROS3PRodasTableau, Rodas3RodasTableau, Rodas3PRodasTableau,
+    ROS3PRodasTableau, Rodas3RodasTableau, Rodas3dRodasTableau, Rodas3PRodasTableau,
     RosShamp4RodasTableau, Veldd4RodasTableau, Velds4RodasTableau,
     GRK4TRodasTableau, GRK4ARodasTableau, Ros4LStabRodasTableau,
     ROS2RodasTableau, ROS2PRRodasTableau, ROS2SRodasTableau,
diff --git a/lib/OrdinaryDiffEqRosenbrockTableaus/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrockTableaus/src/rosenbrock_tableaus.jl
@@ -340,6 +340,48 @@ function Rodas3RodasTableau(::Type{T}, ::Type{T2}) where {T, T2}
     return RodasTableau(A, C, gamma, c, d, H, b, btilde)
 end
 
+################################################################################
+# Rodas3d (4-stage, hand-written)
+################################################################################
+
+"""
+    Rodas3dRodasTableau(T, T2)
+
+A 3rd order stiffly accurate Rosenbrock method with 4 stages and an embedded
+2nd order method for adaptivity. Constructed with the damping parameter
+γ = 0.57281606 so that R(±∞) = 0, which damps both stable and unstable modes.
+Designed for steady-state focused DAE integration such as the semi-implicit
+continuous Newton method, where damping speed matters more than trajectory
+accuracy.
+Reference: Yu, R., Gu, W., Xu, Y., Lu, S. (2024). Semi-implicit Continuous
+Newton Method for Power Flow Analysis. arXiv:2312.02809.
+Coefficients: https://github.com/rzyu45/MATPOWER-SICNM (coeff_rodas3d.m).
+"""
+function Rodas3dRodasTableau(::Type{T}, ::Type{T2}) where {T, T2}
+    A = zeros(T, 4, 4)
+    A[2, 1] = convert(T, 2.1736562342774159)
+    A[3, 1] = convert(T, 1.745761108723104)
+    A[4, 1] = convert(T, 1.745761108723104)
+    A[4, 3] = convert(T, 1)
+
+    C = zeros(T, 4, 4)
+    C[2, 1] = convert(T, -13.387001858207178)
+    C[3, 1] = convert(T, 0.30442314006596932)
+    C[3, 2] = convert(T, 0.30745278826153299)
+    C[4, 1] = convert(T, 0.57287646414081528)
+    C[4, 2] = convert(T, 0.34771098605699358)
+    C[4, 3] = convert(T, -2.7425340696473901)
+
+    gamma = convert(T2, 0.57281606)
+    c = T2[convert(T2, 0), convert(T2, 1.2451051999132263), convert(T2, 1), convert(T2, 1)]
+    d = T[convert(T, 0.57281606), convert(T, -3.819703409768521), zero(T), zero(T)]
+    b = T[convert(T, 1.745761108723104), zero(T), convert(T, 1), convert(T, 1)]
+    btilde = T[zero(T), zero(T), zero(T), one(T)]
+    H = zeros(T, 0, 4)
+
+    return RodasTableau(A, C, gamma, c, d, H, b, btilde)
+end
+
 ################################################################################
 # Rodas3P (5-stage, with H matrix for dense output)
 ################################################################################
diff --git a/lib/OrdinaryDiffEqRosenbrockTableaus/test/qa/Project.toml b/lib/OrdinaryDiffEqRosenbrockTableaus/test/qa/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqRosenbrockTableaus = "b4bd8bb3-f80f-41d2-9b21-73a655b304b9"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
