Add tests, docs, and monorepo wiring for GlobalDiffEq

ChrisRackauckas · ChrisRackauckas · commit 2c286784bd33 · 2026-07-17T05:47:36.000-04:00
Adds GLEE and companion-estimator test suites (convergence orders,
estimate fidelity, adaptive stepping, integrator interface, function
evaluation counts, argument validation), a Global Error Control docs
section with all public names rendered, and the docs env/make wiring.
The GLEE __init override extends plain problems to the partitioned
state via an exact-arity invoke into OrdinaryDiffEqCore so the
integrator interface works directly. QA's rendered-docs check points at
the monorepo docs and excludes reexported DiffEqBase names, which are
documented at their owners.

Co-Authored-By: Chris Rackauckas &lt;accounts@chrisrackauckas.com&gt;
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -2,6 +2,7 @@
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+GlobalDiffEq = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c"
 ImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqAMF = "08082164-f6a1-4363-b3df-3aa6fcf571ad"
@@ -39,6 +40,7 @@ OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 [sources]
 DiffEqBase = {path = "../lib/DiffEqBase"}
 DiffEqDevTools = {path = "../lib/DiffEqDevTools"}
+GlobalDiffEq = {path = "../lib/GlobalDiffEq"}
 ImplicitDiscreteSolve = {path = "../lib/ImplicitDiscreteSolve"}
 OrdinaryDiffEqAMF = {path = "../lib/OrdinaryDiffEqAMF"}
 OrdinaryDiffEqAdamsBashforthMoulton = {path = "../lib/OrdinaryDiffEqAdamsBashforthMoulton"}
@@ -76,6 +78,7 @@ OrdinaryDiffEqVerner = {path = "../lib/OrdinaryDiffEqVerner"}
 DiffEqBase = "7.6.1"
 DiffEqDevTools = "3, 4.0"
 Documenter = "0.27, 1"
+GlobalDiffEq = "1.3"
 ImplicitDiscreteSolve = "2"
 OrdinaryDiffEq = "7"
 OrdinaryDiffEqAMF = "2"
diff --git a/docs/make.jl b/docs/make.jl
@@ -36,6 +36,7 @@ using OrdinaryDiffEqStabilizedRK
 using OrdinaryDiffEqSymplecticRK
 using OrdinaryDiffEqTsit5
 using OrdinaryDiffEqVerner
+using GlobalDiffEq
 
 cp(joinpath(@__DIR__, "Manifest.toml"), joinpath(@__DIR__, "src", "assets", "Manifest.toml"), force = true)
 cp(joinpath(@__DIR__, "Project.toml"), joinpath(@__DIR__, "src", "assets", "Project.toml"), force = true)
@@ -83,6 +84,7 @@ makedocs(
         OrdinaryDiffEqAMF,
         ImplicitDiscreteSolve,
         DiffEqDevTools,
+        GlobalDiffEq,
     ],
     linkcheck_ignore = [r"https://github.com/JuliaDiff/ForwardDiff.jl"],
     warnonly = [:docs_block, :missing_docs, :eval_block],
diff --git a/docs/pages.jl b/docs/pages.jl
@@ -53,6 +53,9 @@ pages = [
     "Misc Solvers" => [
         "misc.md",
     ],
+    "Global Error Control" => [
+        "globalerrorcontrol/GlobalDiffEq.md",
+    ],
     "APIs" => [
         "api/diffeqbase.md",
         "api/ordinarydiffeqcore.md",
diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md
@@ -0,0 +1,76 @@
+# GlobalDiffEq: Global Error Estimation and Control
+
+Standard adaptive ODE solvers control the *local* error of each step. The
+local tolerances only indirectly control the *global* (accumulated) error of
+the solution, which can grow arbitrarily large over long integrations or on
+unstable problems even when every step satisfies its local tolerance. The
+GlobalDiffEq sublibrary provides solvers and solver wrappers that estimate
+the global error, and in several cases control it to a requested global
+tolerance.
+
+To use these methods:
+
+```julia
+using GlobalDiffEq
+```
+
+## Choosing a method
+
+  - For a solution accompanied by a running, asymptotically correct estimate of
+    its global error at every time point, use the global-error-estimating
+    solvers [`GLEE24`](@ref), [`GLEE35`](@ref) (Constantinescu 2016), or the
+    Dormand-Prince-based [`MM5GEE`](@ref) (Makazaga and Murua 2003). These cost
+    only a few extra stages per step over a plain method of the same order and
+    require nothing beyond the right-hand side `f`.
+  - To *control* the endpoint global error to a tolerance `gtol`, wrap any
+    adaptive solver in [`GlobalErrorTransport`](@ref) (linearized
+    error-transport equation, Jacobian-vector products via automatic
+    differentiation), [`GlobalDefectCorrection`](@ref) (solving for the
+    correction; no Jacobian needed), or [`GlobalAdjoint`](@ref) (adjoint-based,
+    for endpoint functionals; requires SciMLSensitivity and QuadGK to be
+    loaded). Each solves the problem, estimates the endpoint global error, and
+    tightens the local tolerances until the requested global tolerance is met.
+  - [`GlobalRichardson`](@ref) wraps any fixed-step method in global Richardson
+    extrapolation over whole solves, interpreting `abstol` and `reltol` as
+    global tolerances. It is the most robust and most expensive option.
+
+For example, solving with a controlled endpoint global error of `1e-8`:
+
+```julia
+using GlobalDiffEq, OrdinaryDiffEqTsit5
+
+function lorenz!(du, u, p, t)
+    du[1] = 10.0(u[2] - u[1])
+    du[2] = u[1] * (28.0 - u[3]) - u[2]
+    du[3] = u[1] * u[2] - (8 / 3) * u[3]
+end
+prob = ODEProblem(lorenz!, [1.0; 0.0; 0.0], (0.0, 10.0))
+sol = solve(prob, GlobalDefectCorrection(Tsit5(); gtol = 1.0e-8))
+```
+
+or solving while tracking the global error along the trajectory:
+
+```julia
+sol = solve(prob, GLEE35(); abstol = 1.0e-8, reltol = 1.0e-8)
+errs = global_error_estimate(sol)  # global error estimate at every sol.t
+```
+
+## Global-error-estimating solvers
+
+```@docs
+GLEE23
+GLEE24
+GLEE35
+MM5GEE
+global_error_estimate
+```
+
+## Global error controlling wrappers
+
+```@docs
+GlobalRichardson
+GlobalErrorTransport
+GlobalDefectCorrection
+GlobalAdjoint
+adjoint_error_estimate
+```
diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml
@@ -64,10 +64,11 @@ OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"]
+test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "RecursiveArrayTools", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"]
diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md
@@ -1,6 +1,31 @@
 # GlobalDiffEq
 
-[![Build Status](https://github.com/SciML/GlobalDiffEq.jl/workflows/CI/badge.svg)](https://github.com/SciML/GlobalDiffEq.jl/actions?query=workflow%3ACI)
+Differential equation solvers and solver wrappers with global (accumulated)
+error estimates and global error control for the
+[DifferentialEquations.jl common solver interface](https://docs.sciml.ai/DiffEqDocs/stable/).
+Formerly the standalone [GlobalDiffEq.jl](https://github.com/SciML/GlobalDiffEq.jl)
+repository, now maintained as a sublibrary of OrdinaryDiffEq.jl.
 
-This is a repository for differential equation solvers with global error estimates
-for the [DifferentialEquations.jl common solver interface](diffeq.sciml.ai).
+Standard adaptive ODE solvers control the *local* error of each step, which
+does not bound the error accumulated over the whole integration. This package
+provides:
+
+  - `GlobalRichardson`: global Richardson extrapolation of whole solves of any
+    fixed-step method, interpreting `abstol`/`reltol` as global tolerances.
+  - `GLEE23`, `GLEE24`, `GLEE35`: explicit general linear methods with built-in
+    global error estimation (Constantinescu 2016), which propagate the solution
+    together with an asymptotically correct estimate of its global error.
+  - `MM5GEE`: the Makazaga-Murua (2003) Dormand-Prince-based order-5 scheme
+    with a cheap built-in global error estimate.
+  - `GlobalErrorTransport`: global error estimation and `gtol`-based control by
+    integrating the linearized error-transport equation driven by the
+    dense-output defect (Shampine 1986, Berzins 1988, Lang-Verwer 2007).
+  - `GlobalDefectCorrection`: global error estimation and control by solving
+    for the correction (Zadunaisky 1976, Dormand-Duckers-Prince 1984/1989),
+    requiring no Jacobian.
+  - `GlobalAdjoint`: adjoint-based a posteriori endpoint error estimation and
+    control (Cao and Petzold 2004), available as a package extension when
+    SciMLSensitivity and QuadGK are loaded.
+
+See the [OrdinaryDiffEq.jl documentation](https://docs.sciml.ai/OrdinaryDiffEq/stable/)
+for the full API reference.
diff --git a/lib/GlobalDiffEq/src/correction.jl b/lib/GlobalDiffEq/src/correction.jl
@@ -126,4 +126,4 @@ function SciMLBase.__solve(
         estimator, prob, alg.alg, alg.gtol, alg.options, args...;
         abstol, reltol, kwargs...
     )
-end
+end
diff --git a/lib/GlobalDiffEq/src/glee/algorithms.jl b/lib/GlobalDiffEq/src/glee/algorithms.jl
@@ -1,7 +1,7 @@
 # Explicit general linear methods with built-in global error estimation
 # (GLEE methods) from Constantinescu (2016), doi:10.1137/15M1014633.
 abstract type AbstractGLEEAlgorithm <:
-              OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end
+OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end
 
 const _GLEE_DOCS_SHARED = """
 GLEE methods propagate the solution `y` together with an asymptotically
diff --git a/lib/GlobalDiffEq/src/glee/caches.jl b/lib/GlobalDiffEq/src/glee/caches.jl
@@ -1,5 +1,5 @@
 @cache struct GLEECache{uType, yType, rateType, yRateType, uNoUnitsType, TabType} <:
-              OrdinaryDiffEqCore.OrdinaryDiffEqMutableCache
+    OrdinaryDiffEqCore.OrdinaryDiffEqMutableCache
     u::uType
     uprev::uType
     tmp::uType
diff --git a/lib/GlobalDiffEq/src/glee/solve.jl b/lib/GlobalDiffEq/src/glee/solve.jl
@@ -46,12 +46,34 @@ function _glee_extended_problem(prob)
     return SciMLBase.remake(prob; f = extended_f, u0 = u0)
 end
 
+_is_glee_extended(prob) = prob.f.f isa GLEEExtendedRHS
+
+# init/solve on a plain ODEProblem transparently extend it to the partitioned
+# (y, ε) state; the invoke dispatches into OrdinaryDiffEqCore's generic __init
+# (its exact five-argument form, so the DiffEqBase default-algorithm catch-all
+# cannot be selected).
+function SciMLBase.__init(
+        prob::SciMLBase.AbstractODEProblem, alg::AbstractGLEEAlgorithm,
+        timeseries_init = (), ts_init = (), ks_init = ();
+        kwargs...
+    )
+    extended_prob = _is_glee_extended(prob) ? prob : _glee_extended_problem(prob)
+    return invoke(
+        SciMLBase.__init,
+        Tuple{
+            SciMLBase.AbstractODEProblem,
+            OrdinaryDiffEqCore.OrdinaryDiffEqAlgorithm,
+            Any, Any, Any,
+        },
+        extended_prob, alg, timeseries_init, ts_init, ks_init; kwargs...
+    )
+end
+
 function SciMLBase.__solve(
         prob::SciMLBase.AbstractODEProblem, alg::AbstractGLEEAlgorithm, args...;
         kwargs...
     )
-    extended_prob = _glee_extended_problem(prob)
-    integrator = SciMLBase.__init(extended_prob, alg, args...; kwargs...)
+    integrator = SciMLBase.__init(prob, alg, args...; kwargs...)
     SciMLBase.solve!(integrator)
     return integrator.sol
 end
diff --git a/lib/GlobalDiffEq/src/glee/tableaus.jl b/lib/GlobalDiffEq/src/glee/tableaus.jl
@@ -65,16 +65,16 @@ function GLEE23Tableau(::Type{T}, ::Type{T2}) where {T, T2}
     A = [
         0 0 0
         1 0 0
-        1//4 1//4 0
+        1 // 4 1 // 4 0
     ]
     U = [
         1 0
         1 10
         1 -1
     ]
     B = [
-        1//12 1//12 5//6
-        1//12 1//12 -1//6
+        1 // 12 1 // 12 5 // 6
+        1 // 12 1 // 12 -1 // 6
     ]
     return _glee_tableau(T, T2, A, U, B, 0, 2)
 end
@@ -90,22 +90,22 @@ integration.
 function GLEE24Tableau(::Type{T}, ::Type{T2}) where {T, T2}
     A = [
         0 0 0 0
-        3//4 0 0 0
-        1//4 29//60 0 0
-        -21//44 145//44 -20//11 0
+        3 // 4 0 0 0
+        1 // 4 29 // 60 0 0
+        -21 // 44 145 // 44 -20 // 11 0
     ]
     # y-ỹ form as printed in the paper
     U = [
         0 1
-        75//58 -17//58
+        75 // 58 -17 // 58
         0 1
         0 1
     ]
     B = [
-        109//275 58//75 -37//110 1//6
-        3//11 0 75//88 -1//8
+        109 // 275 58 // 75 -37 // 110 1 // 6
+        3 // 11 0 75 // 88 -1 // 8
     ]
-    Uyε, Byε = _yytilde_to_yeps(U, B, 0//1)
+    Uyε, Byε = _yytilde_to_yeps(U, B, 0 // 1)
     return _glee_tableau(T, T2, A, Uyε, Byε, 0, 2)
 end
 
diff --git a/lib/GlobalDiffEq/test/companion_tests.jl b/lib/GlobalDiffEq/test/companion_tests.jl
@@ -0,0 +1,68 @@
+using GlobalDiffEq, OrdinaryDiffEqTsit5, LinearAlgebra
+using Test
+import SciMLBase
+
+lv!(du, u, p, t) = (du[1] = 1.5u[1] - u[1] * u[2]; du[2] = -3.0u[2] + u[1] * u[2]; nothing)
+lv(u, p, t) = [1.5u[1] - u[1] * u[2], -3.0u[2] + u[1] * u[2]]
+const lv_tspan = (0.0, 10.0)
+
+@testset "Companion global error estimators" begin
+    prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan)
+    ref = solve(prob, Tsit5(); abstol = 1.0e-13, reltol = 1.0e-13)
+
+    for alg in (GlobalErrorTransport(Tsit5()), GlobalDefectCorrection(Tsit5()))
+        for tol in (1.0e-4, 1.0e-6)
+            est = global_error_estimate(prob, alg; abstol = tol, reltol = tol)
+            sol = solve(prob, Tsit5(); abstol = tol, reltol = tol)
+            true_err = norm(sol.u[end] - ref.u[end])
+            @test est / true_err ≈ 1 rtol = 0.1
+        end
+    end
+
+    # out-of-place problems
+    prob_oop = ODEProblem(lv, [1.0, 1.0], lv_tspan)
+    for alg in (GlobalErrorTransport(Tsit5()), GlobalDefectCorrection(Tsit5()))
+        est = global_error_estimate(prob_oop, alg; abstol = 1.0e-6, reltol = 1.0e-6)
+        sol = solve(prob_oop, Tsit5(); abstol = 1.0e-6, reltol = 1.0e-6)
+        true_err = norm(sol.u[end] - ref.u[end])
+        @test est / true_err ≈ 1 rtol = 0.1
+    end
+end
+
+@testset "Companion global error control" begin
+    prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan)
+    ref = solve(prob, Tsit5(); abstol = 1.0e-13, reltol = 1.0e-13)
+    gtol = 1.0e-7
+
+    for alg in (
+            GlobalErrorTransport(Tsit5(); gtol),
+            GlobalDefectCorrection(Tsit5(); gtol),
+        )
+        sol = solve(prob, alg; abstol = 1.0e-3, reltol = 1.0e-3)
+        @test SciMLBase.successful_retcode(sol)
+        @test norm(sol.u[end] - ref.u[end]) <= gtol
+    end
+end
+
+@testset "Companion estimator argument validation" begin
+    prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan)
+
+    for Alg in (GlobalErrorTransport, GlobalDefectCorrection)
+        alg = Alg(Tsit5())
+        @test !SciMLBase.allows_arbitrary_number_types(alg)
+        @test !SciMLBase.allowscomplex(alg)
+        @test !SciMLBase.isautodifferentiable(alg)
+        @test_throws ArgumentError Alg(Tsit5(); gtol = -1.0)
+        @test_throws ArgumentError Alg(Tsit5(); maxiters = 0)
+        @test_throws ArgumentError Alg(Tsit5(); safety = 2.0)
+        # gtol is required to solve with the wrapper
+        @test_throws ArgumentError solve(prob, alg)
+
+        callback = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing)
+        callback_prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan; callback)
+        @test_throws ArgumentError global_error_estimate(callback_prob, alg)
+
+        backwards_prob = ODEProblem(lv!, [1.0, 1.0], (10.0, 0.0))
+        @test_throws ArgumentError global_error_estimate(backwards_prob, alg)
+    end
+end
diff --git a/lib/GlobalDiffEq/test/glee_tests.jl b/lib/GlobalDiffEq/test/glee_tests.jl
diff --git a/lib/GlobalDiffEq/test/qa/qa.jl b/lib/GlobalDiffEq/test/qa/qa.jl
diff --git a/lib/GlobalDiffEq/test/runtests.jl b/lib/GlobalDiffEq/test/runtests.jl

Original file line number	Diff line number	Diff line change
`@@ -126,4 +126,4 @@ function SciMLBase.__solve(`
`126`	`126`	`estimator, prob, alg.alg, alg.gtol, alg.options, args...;`
`127`	`127`	`abstol, reltol, kwargs...`
`128`	`128`	`)`
`129`		`-end`
	`129`	`+end`