Skip to content

GlobalDiffEq: add GLEE23/GLEE24/GLEE35 global-error-estimating solvers#3954

Draft
ChrisRackauckas-Claude wants to merge 2 commits into
SciML:masterfrom
ChrisRackauckas-Claude:globaldiffeq-glee
Draft

GlobalDiffEq: add GLEE23/GLEE24/GLEE35 global-error-estimating solvers#3954
ChrisRackauckas-Claude wants to merge 2 commits into
SciML:masterfrom
ChrisRackauckas-Claude:globaldiffeq-glee

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Jul 18, 2026

Copy link
Copy Markdown
Member

Important

Ignore this PR until it has been reviewed by @ChrisRackauckas.

Stacked on #3953 (adjoint). #3929 (repo transfer) has merged and this branch is rebased onto master; only the top commit is new here. Merge #3953 first.

Summary

Adds GLEE23, GLEE24, GLEE35 — the explicit general linear methods with built-in global error estimation of Constantinescu, SIAM J. Numer. Anal. 54(6), 2016 (PETSc's TSGLEE23/24/35) — as proper OrdinaryDiffEqCore solver algorithms (alg_cache/initialize!/perform_step!). Addresses SciML/GlobalDiffEq.jl#6 and SciML/GlobalDiffEq.jl#22.

Design, following the plan in SciML/DifferentialEquations.jl#670 (comment):

  • The methods propagate the partitioned state (y, ε) in the paper's y-ε general linear form, where ε is an asymptotically correct estimate of the global error y_exact − y at every step. solve/init on a plain ODEProblem transparently extend it to ArrayPartition state (the "guiding layer" adjusts the user's input); global_error_estimate(sol) (new export) extracts the estimate trajectory.
  • The per-step increment of ε is an asymptotically correct local error estimate, which drives the standard step-size controllers, so the methods are fully adaptive with the usual abstol/reltol semantics while the global error estimate comes along for free.
  • Tableaus are transcribed exactly from the paper's LaTeX source and cross-checked against PETSc's glee.c (GLEE35's exact rationals exceed Int64 and are constructed as BigInt rationals). The y-ỹ→y-ε form conversion (paper Lemma 4.1) is done at construction. Consistency (ΣB₁ = 1, ΣB₂ = 0, U·[1,0]ᵀ = 1) is preserved exactly.
  • The paper's method A4 (eq. A.2) is defective as printed — it violates the order-2 condition b·c = 1/2 and converges at order 1 in a direct test; it is also absent from PETSc — so it was deliberately not implemented.
  • FSAL structure: stage 1 reuses fsalfirst when the tableau allows (GLEE23), and dense output gets true endpoint derivatives.

GLEE24 and GLEE35 satisfy both of the paper's decoupling conditions (B·U and B·A·U diagonal) and are the recommended choices; GLEE23 is the cheapest (3 stages).

Version 1.4.0.

Validation (run locally)

  • ODEDIFFEQ_TEST_GROUP=Core julia --project=lib/GlobalDiffEq -e 'using Pkg; Pkg.test()': pass — Basic 1/1, traits 3/3, Adjoint 14/14, GLEE solvers 30/30, BigFloat 2/2
  • ODEDIFFEQ_TEST_GROUP=QA julia --project=lib/GlobalDiffEq -e 'using Pkg; Pkg.test()': pass — 22 passed, 1 pre-existing declared broken
  • Convergence on the unstable Prince42 problem (y' = y − sin t + cos t, errors grow like eᵗ), in-place and out-of-place: solution orders 2/2/3 and estimate-accuracy orders 3/3/4 as claimed; est/true-error ratio ≈ 1.00 at the finest fixed step; adaptive solves stay faithful; the integrator interface (init/solve!) works on plain problems.

🤖 Generated with Claude Code

https://claude.ai/code/session_01MfUh8tgp3iBmMv9nkBqoyY

Carries SciML/GlobalDiffEq.jl#56 (GlobalAdjoint, adjoint_error_estimate;
Cao and Petzold 2004, closes SciML/GlobalDiffEq.jl#4) into the
monorepo, restructured as the package extension
GlobalDiffEqSciMLSensitivityExt on SciMLSensitivity + QuadGK weakdeps.
sensealg defaults to nothing and resolves to
QuadratureAdjoint(autojacvec = true) when the extension loads; solving
without the extension raises an instructive error. Public API and
semantics otherwise match the original PR.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Implements the explicit general linear methods with built-in global
error estimation from Constantinescu (2016), SIAM J. Numer. Anal. 54(6)
(arXiv:1503.05166; PETSc's TSGLEE23/TSGLEE24/TSGLEE35), addressing
SciML/GlobalDiffEq.jl#6 and SciML/GlobalDiffEq.jl#22, as proper
OrdinaryDiffEqCore algorithms. The methods propagate the partitioned
state (y, ε) where ε is an asymptotically correct global error
estimate; solve/init on plain ODEProblems transparently extend to the
ArrayPartition form and global_error_estimate(sol) extracts the
estimates. The per-step ε increment is an asymptotically correct local
error estimate driving standard step-size adaptivity. The paper's
method A4 is defective as printed (violates b·c = 1/2, converges at
order 1, absent from PETSc) and is deliberately not implemented.

Verified locally on the unstable Prince42 problem: solution orders
2/2/3, estimate accuracy orders 3/3/4, est/true-error ratio 1.00.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants