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:

using GlobalDiffEq

The global-error-estimating solvers GLEE24, GLEE35 (Constantinescu 2016), and the Dormand-Prince-based MM5GEE (Makazaga and Murua 2003) carry a running, asymptotically correct estimate of the solution's global error at every time point, at the cost of a few extra stages per step. They report it through the standard SciMLBase global-error interface: solving returns an ordinary solution (sol.u and sol(t) are the solution), with the estimate in sol.global_error — sol.global_error[i] is the estimated global error of sol.u[i] at sol.t[i], and SciMLBase.has_global_error is true for these algorithms. global_error_estimate returns sol.global_error.

GlobalErrorEstimation wraps any adaptive solver with global error estimation and gtol-based control by integrating a defect-driven companion equation. Two orthogonal choices configure it: GlobalErrorEquation selects the companion ODE — the nonlinear DefectCorrection (no Jacobian) or the linearized ErrorTransport (matrix-free Jacobian-vector products) — and GlobalErrorMode selects how it is integrated — InterpolatingMode (store the dense solution, integrate the companion in a second pass) or SimultaneousMode (advance the companion inline over each step for O(1) extra memory and no second solve). global_error_estimate(prob, alg) exposes the standalone estimator.

To control the endpoint global error to a tolerance gtol, wrap any adaptive solver in GlobalAdjoint (adjoint-based, for endpoint functionals; requires SciMLSensitivity to be loaded); adjoint_error_estimate exposes the standalone estimator. Its GlobalAdjointScope selects whether the estimate is reported only at the endpoint (EndpointError) or along the whole trajectory (TrajectoryError, filling sol.global_error at every saved time), and its GlobalAdjointControl selects whether gtol is met by tolerance tightening (ToleranceRefinement) or by re-solving non-adaptively on the accepted step grid (StepGridRefinement).

GlobalRichardson 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 while tracking the global error along the trajectory:

using GlobalDiffEq

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]
    return
end
prob = ODEProblem(lorenz!, [1.0; 0.0; 0.0], (0.0, 10.0))
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

GlobalDiffEq.GLEE23 — Type

GLEE23()

3-stage, second-order explicit general linear method with global error estimation (Constantinescu 2016, eq. (4.6); PETSc's TSGLEE23). The cheapest GLEE method: only the first decoupling condition holds, so prefer GLEE24 for long-time integration or mildly stiff problems.

GLEE methods propagate the solution y together with an asymptotically correct estimate ε of its global (accumulated) error, at the cost of a few extra stages per step. Solving any ODEProblem with a GLEE method returns an ordinary solution — sol.u[i] is the solution at sol.t[i] and sol(t) interpolates it at full order — with the global error estimate reported through the standard SciMLBase interface: SciMLBase.has_global_error is true, and sol.global_error[i] is the estimated global error of sol.u[i] (also available as global_error_estimate(sol)). The per-step increment of ε is an asymptotically correct local error estimate, which drives standard step-size adaptivity, so local tolerances behave exactly as for ordinary adaptive Runge-Kutta methods while the global error is estimated for free.

Only explicit, mass-matrix-free ODEs are supported. The reference for the methods and their theory is:

Emil M. Constantinescu, Estimating Global Errors in Time Stepping, SIAM Journal on Numerical Analysis 54(6), 2016. arXiv:1503.05166

GlobalDiffEq: Global Error Estimation and Control

Global-error-estimating solvers

Global error controlling wrappers