DiffEqCallbacks.jl

ManifoldProjection(
    manifold; nlsolve = missing, save = true, autonomous = nothing,
    manifold_jacobian = nothing, autodiff = nothing,
    resid_prototype = nothing, kwargs...) -> DiscreteCallback

In many cases, you may want to declare a manifold on which a solution lives. Mathematically, a manifold M is defined by a function g as the set of points where g(u) = 0. An embedded manifold can be a lower dimensional object which constrains the solution. For example, g(u) = E(u) - C where E is the energy of the system in state u, meaning that the energy must be constant (energy preservation). Thus by defining the manifold the solution should live on, you can retain desired properties of the solution.

ManifoldProjection projects the solution of the differential equation to the chosen manifold g, conserving a property while conserving the order. It is a consequence of convergence proofs both in the deterministic and stochastic cases that post-step projection to manifolds keep the same convergence rate, thus any algorithm can be easily extended to conserve properties. If the solution is supposed to live on a specific manifold or conserve such property, this guarantees the conservation law without modifying the convergence properties.

Arguments

• manifold: residual function defining the manifold. For an in-place problem, define manifold(resid, u, p) or manifold(resid, u, p, t). For an out-of-place problem, define manifold(u, p) or manifold(u, p, t). The residual is zero on the manifold.

Keywords

• nlsolve = missing: use the built-in single-factorization projection algorithm. Pass a nonlinear solver in the NonlinearSolve.jl format to select another algorithm. Pass nothing to use the NonlinearSolve polyalgorithm.
• save::Bool = true: save immediately after the projection.
• autonomous = nothing: whether manifold omits the time argument. Pass Val(true) or Val(false) to avoid runtime branching; nothing infers the form during initialization.
• resid_prototype = nothing: prototype defining the residual shape for an in-place problem. Without one, the residual is assumed to have the same shape as u.
• autodiff = nothing: DifferentiationInterface automatic-differentiation backend used when manifold_jacobian is not supplied.
• manifold_jacobian = nothing: analytic Jacobian of manifold with respect to u. Use the same calling form as manifold, with the Jacobian output as the first argument for an in-place problem.
• kwargs...: additional keywords passed to the built-in projection algorithm or to NonlinearSolve.jl if nlsolve is not missing.

Returns

• DiscreteCallback: a callback that projects the state after each accepted step. If the nonlinear projection does not converge, the callback terminates the integrator with the projection solver's unsuccessful return code.

Throws

• ErrorException: during callback initialization if both manifold_jacobian and autodiff are nothing.

Saveat Warning

Note that the ManifoldProjection callback modifies the endpoints of the integration intervals and thus breaks assumptions of internal interpolations. Because of this, the values for given by saveat will not be order-matching. However, the interpolation error can be proportional to the change by the projection, so if the projection is making small changes then one is still safe. However, if there are large changes from each projection, you should consider only saving at stopping/projection times. To do this, set tstops to the same values as saveat. There is a performance hit by doing so because now the integrator is forced to stop at every saving point, but this is guaranteed to match the order of the integrator even with the ManifoldProjection.

References

[1] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Berlin ; New York :Springer, 2002.

Examples

using ADTypes, DiffEqCallbacks, OrdinaryDiffEq

function unit_circle(resid, u, p, t)
    resid[1] = sum(abs2, u) - 1
end

function rotation!(du, u, p, t)
    du[1] = -u[2]
    du[2] = u[1]
end

prob = ODEProblem(rotation!, [1.0, 0.0], (0.0, 10.0))
cb = ManifoldProjection(unit_circle; resid_prototype = [0.0], autodiff = AutoFiniteDiff())

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

SavingCallback(save_func, saved_values::SavedValues;
    saveat = Vector{eltype(saved_values.t)}(),
    save_everystep = isempty(saveat),
    save_start = save_everystep || isempty(saveat) || saveat isa Number,
    save_end = save_everystep || isempty(saveat) || saveat isa Number,
    tdir = 1) -> DiscreteCallback

The saving callback lets you define a function save_func(u, t, integrator) which returns quantities of interest that shall be saved.

Arguments

• save_func: function called as save_func(u, t, integrator). It must return a value compatible with eltype(saved_values.saveval) and must not return a view of u.
• saved_values::SavedValues: storage whose time and value element types match the integration time and the output of save_func.

Keywords

• saveat = Vector{eltype(saved_values.t)}(): selected integration times, or a scalar interval at which to evaluate save_func throughout the problem time span.
• save_everystep::Bool = isempty(saveat): whether to save after every accepted step.
• save_start::Bool = ...: whether to save at the initial condition.
• save_end::Bool = ...: whether to save at the final time.
• tdir = 1: integration direction used to order saveat. Set this to sign(tspan[end] - tspan[1]) for reverse-time problems.

The output values are saved into saved_values. Time points are found via saved_values.t and the values are saved_values.saveval.

Returns

• DiscreteCallback: a callback that evaluates save_func at the requested times and appends the results to saved_values.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
saved_values = SavedValues(Float64, Float64)
cb = SavingCallback((u, t, integrator) -> u^2, saved_values; saveat = 0.0:0.25:1.0)

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

SavedValues{tType<:Real, savevalType}

Container used by SavingCallback to store saved time points and user-defined values.

Fields

• t::Vector{tType}: saved time points.
• saveval::Vector{savevalType}: values returned by the saving function.

Construct empty storage with SavedValues(tType, savevalType). The callback appends to both vectors in place; do not mutate them while a solve is active.

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)

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)

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)

FunctionCallingCallback(func;
    funcat = Vector{Float64}(),
    func_everystep = isempty(funcat),
    func_start = true,
    tdir = 1) -> DiscreteCallback

The function calling callback lets you define a function func(u,t,integrator) which gets called at the time points of interest.

Arguments

• func: function called as func(u, t, integrator) at each selected time. Its return value is ignored, and it should not modify u or the integrator.

Keywords

• funcat = Vector{Float64}(): selected integration times, or a scalar interval at which to call func throughout the problem time span.
• func_everystep::Bool = isempty(funcat): whether to call func after every accepted step.
• func_start::Bool = true: whether to call func at the initial condition.
• tdir = 1: integration direction used to order funcat. Set this to sign(tspan[end] - tspan[1]) for reverse-time problems.

Returns

• DiscreteCallback: a callback that calls func without modifying the integrator state.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

seen = Float64[]
func = (u, t, integrator) -> push!(seen, t)
cb = FunctionCallingCallback(func; funcat = 0.0:0.25:1.0)

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
sol = solve(prob, Tsit5(); callback = cb)

TerminateSteadyState(abstol = 1.0e-8, reltol = 1.0e-6, test = allDerivPass;
    min_t = nothing, wrap_test::Val = Val(true)) -> DiscreteCallback

TerminateSteadyState can be used to solve the problem for the steady-state by running the solver until the derivatives of the problem converge to 0 or tspan[2] is reached. This is an alternative approach to root finding; see the Steady State Solvers documentation.

Arguments

• abstol = 1.0e-8: absolute termination tolerance. It may be a scalar or an array with the same length as the state.
• reltol = 1.0e-6: relative termination tolerance. It may be a scalar or an array with the same length as the state.
• test = allDerivPass: function that evaluates the termination condition. By default, every derivative must be smaller than abstol or the corresponding state magnitude times reltol. A custom wrapped test must accept integrator, abstol, reltol, and min_t.

Keywords

• min_t = nothing: optional minimum integration time before termination is allowed.
• wrap_test::Val = Val(true): with Val(true), call test as test(integrator, abstol, reltol, min_t). With Val(false), use test directly as the callback condition test(u, t, integrator).

Returns

• DiscreteCallback: a callback that terminates the integrator when test returns true.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

f(u, p, t) = 1 - u
prob = ODEProblem(f, 0.0, (0.0, 100.0))
cb = TerminateSteadyState(1.0e-8, 1.0e-8)

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

IterativeCallback(time_choice, user_affect!, tType = Float64;
    initial_affect = false, initialize = ..., kwargs...) -> DiscreteCallback

Construct a callback that applies user_affect! at the sequence of integration times returned by time_choice.

Arguments

• time_choice: a function time_choice(integrator) that returns the next callback time or nothing to stop scheduling further affects.
• user_affect!: a function user_affect!(integrator) applied at each scheduled time.
• tType::Type = Float64: type used to store the next callback time. Set this to the problem time type when it is not Float64.

Keywords

• initial_affect::Bool = false: apply user_affect! during callback initialization at the initial integration time before asking time_choice for the next time.
• initialize = ...: callback initialization function called as initialize(callback, u, t, integrator) before the initial affect or first scheduled stop. By default, it marks a derivative discontinuity according to initial_affect.
• kwargs...: keyword arguments forwarded to DiscreteCallback.

Returns

• DiscreteCallback: a callback that schedules each time returned by time_choice until it returns nothing.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

count = Ref(0)
hits = Float64[]
time_choice = integrator -> (count[] += 1; count[] <= 3 ? integrator.t + 0.1 : nothing)
affect! = integrator -> push!(hits, integrator.t)
cb = IterativeCallback(time_choice, affect!)

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
sol = solve(prob, Tsit5(); callback = cb)

PeriodicCallback(f, Δt::Number; phase = 0, initial_affect = false,
    final_affect = false, initialize = ..., kwargs...) -> DiscreteCallback

Construct a callback that applies f at regular intervals of integration time. Scheduled stops are separated by Δt and are offset from the initial time by phase. When initial_affect = true, f is also applied during callback initialization.

Arguments

• f: a function f(integrator) applied at each periodic stop.
• Δt::Number: signed integration-time period. Its sign must match the integration direction.

Keywords

• phase = 0: nonnegative offset of scheduled periodic stops from the initial integration time. A negative phase throws an ArgumentError.
• initial_affect::Bool = false: apply f during callback initialization at the initial integration time.
• final_affect::Bool = false: apply f when the integrator finishes, even when the final time is not a periodic stop.
• initialize = ...: callback initialization function called as initialize(callback, u, t, integrator) before periodic stops are scheduled. By default, it marks a derivative discontinuity according to initial_affect.
• kwargs...: keyword arguments forwarded to DiscreteCallback.

Returns

• DiscreteCallback: a callback that schedules f at periodic integration-time stops.

Throws

• ArgumentError: if phase < 0.
• AssertionError: during callback initialization if the sign of Δt does not match the integration direction.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

samples = Float64[]
affect! = integrator -> push!(samples, integrator.u)
cb = PeriodicCallback(affect!, 0.1; initial_affect = true)

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
sol = solve(prob, Tsit5(); callback = cb)

PresetTimeCallback(tstops, user_affect!; initialize = INITIALIZE_DEFAULT,
    filter_tstops = true, sort_inplace = false, kwargs...) -> DiscreteCallback

Construct a callback that schedules user_affect! at the supplied integration-time stops.

Arguments

• tstops::Union{Number, AbstractVector}: one time or a collection of callback times. Vector inputs are sorted before use.
• user_affect!: a function user_affect!(integrator) applied at each scheduled stop.

Keywords

• initialize = INITIALIZE_DEFAULT: callback initialization function called as initialize(callback, u, t, integrator) before the stops are scheduled.
• filter_tstops::Bool = true: schedule only stops strictly inside the integration interval. Set this to false to schedule all supplied stops, including values outside that interval.
• sort_inplace::Bool = false: sort a vector tstops in place. By default, sort a copy and leave the supplied vector unchanged.
• kwargs...: keyword arguments forwarded to DiscreteCallback.

Returns

• DiscreteCallback: a callback that schedules the requested stops and calls user_affect!.

Throws

• ArgumentError: if tstops is neither a number nor a vector.

Examples

using DiffEqCallbacks, OrdinaryDiffEq

hits = Float64[]
affect! = integrator -> push!(hits, integrator.t)
cb = PresetTimeCallback([0.25, 0.5, 0.75], affect!)

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
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)