GitHub

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.

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

source
GlobalDiffEq.GLEE24 — Type
GLEE24()

4-stage, second-order explicit general linear method with global error estimation (Constantinescu 2016, eq. (A.3); PETSc's TSGLEE24). Satisfies both decoupling conditions (B·U and B·A·U diagonal), which keeps the error estimate faithful in long-time integration; this is the recommended second-order GLEE method.

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

source
GlobalDiffEq.GLEE35 — Type
GLEE35()

5-stage, third-order explicit general linear method with global error estimation (Constantinescu 2016, eq. (4.9); PETSc's TSGLEE35). Satisfies both decoupling conditions and has a large negative-real-axis stability region; this is the recommended third-order GLEE method.

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

source
GlobalDiffEq.MM5GEE — Type
MM5GEE()

Fifth-order Dormand-Prince-based scheme with cheap global error estimation (Makazaga and Murua, BIT Numerical Mathematics 43, 2003). Propagates the standard DOPRI5 solution together with an order-6 companion solution through three extra stages (9 function evaluations per step in total), whose difference is an asymptotically correct global error estimate. This is the recommended method of the Runge-Kutta triple family of Dormand, Duckers and Prince, as it is the coefficient-complete published scheme of that type.

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

  • J. Makazaga and A. Murua, New Runge-Kutta based schemes for ODEs with cheap global error estimation, BIT Numerical Mathematics 43 (2003).
source
GlobalDiffEq.global_error_estimate — Function
global_error_estimate(sol)
global_error_estimate(sol, i)

Extract the global error estimate from a solution computed with a GLEE method (GLEE23, GLEE24, GLEE35, MM5GEE).

These solvers report the global error through the standard SciMLBase interface: sol.global_error is a vector matching sol.u, with sol.global_error[i] the estimated global error of sol.u[i] at sol.t[i] (see SciMLBase.has_global_error). global_error_estimate(sol) returns that vector and global_error_estimate(sol, i) its i-th element.

source

Global error controlling wrappers