Step Control Callbacks

StepsizeLimiter(dtFE; safety_factor = 9 // 10, max_step = false,
    cached_dtcache = 0.0) -> DiscreteCallback

In many cases, there is a known maximal stepsize for which the computation is stable and produces correct results. For example, in hyperbolic PDEs one normally needs to ensure that the stepsize stays below some $\Delta t_{FE}$ determined by the CFL condition. For nonlinear hyperbolic PDEs this limit can be a function dtFE(u,p,t) which changes throughout the computation. The stepsize limiter lets you pass a function which will adaptively limit the stepsizes to match these constraints.

Arguments

• dtFE: function called as dtFE(u, p, t) to compute the current maximum stable step.

Keywords

• safety_factor = 9 // 10: factor applied to the maximum returned by dtFE.
• max_step::Bool = false: when true, set every proposed step to safety_factor * dtFE(u, p, t), including for a non-adaptive solver.
• cached_dtcache = 0.0: initial cache for the unconstrained step. Set it to a value with the problem time type when that type is not Float64.

Returns

• DiscreteCallback: a callback that updates integrator.opts.dtmax before every step.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

f(u, p, t) = -u
prob = ODEProblem(f, 1.0, (0.0, 1.0))
dtFE(u, p, t) = 0.05
cb = StepsizeLimiter(dtFE; safety_factor = 0.8)

sol = solve(prob, Tsit5(); callback = cb)

GeneralDomain(
    g, u = nothing; save = true, abstol = nothing, scalefactor = nothing,
    autonomous = nothing, domain_jacobian = nothing, manifold_jacobian = missing,
    nlsolve_kwargs = (; abstol = 10 * eps()), kwargs...) -> CallbackSet

A GeneralDomain callback in DiffEqCallbacks.jl generalizes the concept of a PositiveDomain callback to arbitrary domains.

Domains are specified by

• in-place functions g(resid, u, p) or g(resid, u, p, t) if the corresponding ODEProblem is an inplace problem, or
• out-of-place functions g(u, p) or g(u, p, t) if the corresponding ODEProblem is an out-of-place problem.

The function calculates residuals of a state vector u at time t relative to that domain, with p the parameters of the corresponding integrator.

As for PositiveDomain, steps are accepted if residuals of the extrapolated values at the next time step are below a certain tolerance. Moreover, this callback is automatically coupled with a ManifoldProjection that keeps all calculated state vectors close to the desired domain, but in contrast to a PositiveDomain callback the nonlinear solver in a ManifoldProjection cannot guarantee that all state vectors of the solution are actually inside the domain. Thus, a PositiveDomain callback should generally be preferred.

Arguments

• g: the implicit definition of the domain as a function as described above which is zero when the value is in the domain.
• u = nothing: a prototype of the state vector of the integrator. A copy is saved and extrapolated values are written to it. If it is not specified, every application of the callback allocates a new copy of the state vector.

Keywords

• save::Bool = true: whether to save immediately after applying the domain callback.
• abstol = nothing: tolerance below which residuals are accepted. Element-wise tolerances are allowed. If it is not specified, every application of the callback uses the current absolute tolerances of the integrator.
• scalefactor = nothing: factor by which an unaccepted time step is reduced. If it is not specified, time steps are halved.
• autonomous = nothing: whether g is an autonomous function of the form g(resid, u, p) or g(u, p). If it is not specified, it is determined automatically.
• domain_jacobian = nothing: analytic Jacobian of g with respect to the state, using the same calling form as g and a leading Jacobian output for an in-place problem.
• manifold_jacobian = missing: unsupported compatibility keyword. Supplying any value throws an ArgumentError; use domain_jacobian instead.
• nlsolve_kwargs = (; abstol = 10 * eps()): keywords passed to the nonlinear solver in ManifoldProjection. The default is (; abstol = 10 * eps()).
• kwargs...: additional keywords passed to ManifoldProjection, including autodiff, nlsolve, and resid_prototype. Either domain_jacobian or autodiff must be provided.

Returns

• CallbackSet: a manifold projection followed by a discrete callback that restricts the proposed step to the requested domain.

Throws

• ArgumentError: if manifold_jacobian is supplied, or if the callback is applied with a non-adaptive integrator.
• DimensionMismatch: if an element-wise abstol does not match the residual length.
• ErrorException: during callback initialization if both domain_jacobian and the forwarded autodiff keyword are nothing.

References

Shampine, Lawrence F., Skip Thompson, Jacek Kierzenka and G. D. Byrne. Non-negative solutions of ODEs. Applied Mathematics and Computation 170 (2005): 556-569.

Examples

using ADTypes, DiffEqCallbacks, OrdinaryDiffEq

function nonnegative_residual(resid, u, p, t)
    @. resid = max(-u, 0)
end

prob = ODEProblem((du, u, p, t) -> (du .= -u), [1.0, 2.0], (0.0, 2.0))
cb = GeneralDomain(
    nonnegative_residual, [1.0, 2.0]; abstol = 1.0e-8,
    autodiff = AutoForwardDiff()
)
sol = solve(prob, Tsit5(); callback = cb)

PositiveDomain(u = nothing; save = true, abstol = nothing,
    scalefactor = nothing) -> DiscreteCallback

Especially in biology and other natural sciences, a desired property of dynamical systems is the positive invariance of the positive cone, i.e. non-negativity of variables at time $t_0$ ensures their non-negativity at times $t \geq t_0$ for which the solution is defined. However, even if a system satisfies this property mathematically it can be difficult for ODE solvers to ensure it numerically, as these MATLAB examples show.

To deal with this problem, one can specify isoutofdomain=(u,p,t) -> any(x -> x < 0, u) as an additional solver option, which will reject any step that leads to negative values and reduce the next time step. However, since this approach only rejects steps and hence calculations might be repeated multiple times until a step is accepted, it can be computationally expensive.

Another approach is taken by a PositiveDomain callback in DiffEqCallbacks.jl, which is inspired by Shampine et al.'s paper about non-negative ODE solutions. It reduces the next step by a certain scale factor until the extrapolated value at the next time point is non-negative with a certain tolerance. Extrapolations are cheap to compute but might be inaccurate, so if a time step is changed it is additionally reduced by a safety factor of 0.9. Since extrapolated values are only non-negative up to a certain tolerance and in addition actual calculations might lead to negative values, also any negative values at the current time point are set to 0. Hence, by this callback non-negative values at any time point are ensured in a computationally cheap way, but the quality of the solution depends on how accurately extrapolations approximate next time steps.

Please note, that the system should be defined also outside the positive domain, since even with these approaches, negative variables might occur during the calculations. Moreover, one should follow Shampine's et al. advice and set the derivative $x'_i$ of a negative component $x_i$ to $\max \{0, f_i(x, t)\}$, where $t$ denotes the current time point with state vector $x$ and $f_i$ is the $i$-th component of function $f$ in an ODE system $x' = f(x, t)$.

Arguments

• u = nothing: a prototype of the state vector of the integrator. A copy is saved and extrapolated values are written to it. If it is not specified, every application of the callback allocates a new copy of the state vector.

Keywords

• save::Bool = true: whether to save immediately after applying the domain callback.
• abstol = nothing: tolerance above the negative of which extrapolated values are accepted. Element-wise tolerances are allowed. If it is not specified, every application of the callback uses the current absolute tolerances of the integrator.
• scalefactor = nothing: factor by which an unaccepted time step is reduced. If it is not specified, time steps are halved.

Returns

• DiscreteCallback: a callback that restricts proposed steps to the positive domain and replaces negative entries in each accepted state with zero.

Throws

• ArgumentError: if the callback is applied with a non-adaptive integrator.
• DimensionMismatch: if an element-wise abstol does not match the state length.

References

Shampine, Lawrence F., Skip Thompson, Jacek Kierzenka and G. D. Byrne. Non-negative solutions of ODEs. Applied Mathematics and Computation 170 (2005): 556-569.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

f(u, p, t) = -u
prob = ODEProblem(f, [1.0], (0.0, 2.0))
cb = PositiveDomain()

sol = solve(prob, Tsit5(); callback = cb)

AutoAbstol(save = true; init_curmax = 0.0) -> DiscreteCallback

Construct a callback that updates integrator.opts.abstol after every accepted step to the largest magnitude observed in the state so far, multiplied by integrator.opts.reltol.

Arguments

• save::Bool = true: save the solution immediately before the callback affect. Set this to false when another callback controls saving.

Keywords

• init_curmax = 0.0: initial maximum state magnitude. A zero value is replaced during initialization with the integrator's configured abstol; arrays update elementwise.

Returns

• DiscreteCallback: a callback that updates the absolute tolerance after each accepted step without marking the state as modified.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

f(u, p, t) = 0.5u
prob = ODEProblem(f, 1.0, (0.0, 2.0))
cb = AutoAbstol(; init_curmax = 1.0e-8)

sol = solve(prob, Tsit5(); callback = cb, reltol = 1.0e-6)