DiffEqCallbacks.jl

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

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: The residual function for the manifold. If the ODEProblem is an inplace problem, then g must be an inplace function of form g(resid, u, p) or g(resid, u, p, t) and similarly if the ODEProblem is an out-of-place problem then g must be a function of form g(u, p) or g(u, p, t).

Keyword Arguments

• nlsolve: Defaults to a special single-factorize algorithm (denoted by missing) that would work in most cases (See [1] for details). Alternatively, a nonlinear solver as defined in the NonlinearSolve.jl format can be specified. Additionally if NonlinearSolve.jl is loaded and nothing is specified a polyalgorithm is used.
• save: Whether to do the standard saving (applied after the callback)
• autonomous: Whether g is an autonomous function of the form g(resid, u, p) or g(u, p). Specify it as Val(::Bool) to disable runtime branching. If nothing, we attempt to infer it from the signature of g.
• residual_prototype: The size of the manifold residual. If it is not specified, we assume it to be same as u in the inplace case. Else we run a single evaluation of manifold to determine the size.
• kwargs: All other keyword arguments are passed to NonlinearSolve.jl if nlsolve is not missing.
• autodiff: The autodifferentiation algorithm to use to compute the Jacobian if manifold_jacobian is not specified. This must be specified if manifold_jacobian is not specified.
• manifold_jacobian: The Jacobian of the manifold (wrt the state). This has the same signature as manifold and the first argument is the Jacobian if inplace.

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 guerenteed 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.

SavingCallback(save_func, saved_values::SavedValues;
    saveat = Vector{eltype(saved_values.t)}(),
    save_everystep = isempty(saveat),
    save_start = true,
    tdir = 1)

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(u, t, integrator) returns the quantities which shall be saved. Note that this should allocate the output (not as a view to u).
• saved_values::SavedValues is the types that save_func will return, i.e. save_func(t, u, integrator)::savevalType. It's specified via SavedValues(typeof(t),savevalType), i.e. give the type for time and the type that save_func will output (or higher compatible type).

Keyword Arguments

• saveat mimics saveat in solve from solve.
• save_everystep mimics save_everystep from solve.
• save_start mimics save_start from solve.
• save_end mimics save_end from solve.
• tdir should be sign(tspan[end]-tspan[1]). It defaults to 1 and should be adapted if tspan[1] > tspan[end].

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

SavedValues{tType<:Real, savevalType}

A struct used to save values of the time in t::Vector{tType} and additional values in saveval::Vector{savevalType}.

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

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 additional solver option, which will reject any step that leads to non-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's et al. 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: A prototype of the state vector of the integrator. A copy of it 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.

Keyword Arguments

• save: Whether to do the standard saving (applied after the callback).
• abstol: Tolerance up to, which negative 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: Factor by which an unaccepted time step is reduced. If it is not specified, time steps are halved.

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.

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

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: A prototype of the state vector of the integrator. A copy of it 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.

Keyword Arguments

• save: Whether to do the standard saving (applied after the callback).
• abstol: Tolerance up to, 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: Factor by which an unaccepted time step is reduced. If it is not specified, time steps are halved.
• autonomous: Whether g is an autonomous function of the form g(resid, u, p). If it is not specified, it is determined automatically.
• kwargs: All other keyword arguments are passed to ManifoldProjection.
• nlsolve_kwargs: All keyword arguments are passed to the nonlinear solver in ManifoldProjection. The default is (; abstol = 10 * eps()).
• domain_jacobian: The Jacobian of the domain (wrt the state). This has the same signature as g and the first argument is the Jacobian if inplace. This corresponds to the manifold_jacobian argument of ManifoldProjection. Note that passing a manifold_jacobian is not supported for GeneralDomain and results in an error.

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.

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

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 is the maximal timestep and is calculated using the previous t and u.

Keyword Arguments

• safety_factor is the factor below the true maximum that will be stepped to which defaults to 9//10.
• max_step=true makes every step equal to safety_factor*dtFE(u,p,t) when the solver is set to adaptive=false.
• cached_dtcache should be set to match the type for time when not using Float64 values.

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

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

• func(u, t, integrator) is the function to be called.
• funcat values or interval that the function is sure to be evaluated at.
• func_everystep whether to call the function after each integrator step.
• func_start whether the function is called at the initial condition.
• tdir should be sign(tspan[end]-tspan[1]). It defaults to 1 and should be adapted if tspan[1] > tspan[end].

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

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 section).

Arguments

• abstol and reltol are the absolute and relative tolerance, respectively. These tolerances may be specified as scalars or as arrays of the same length as the states of the problem.
• test represents the function that evaluates the condition for termination. The default condition is that all derivatives should become smaller than abstol or the states times reltol. The user can pass any other function to implement a different termination condition. Such function should take four arguments: integrator, abstol, reltol, and min_t.
• wrap_test can be set to Val(false), in which case test must have the definition test(u, t, integrator). Otherwise, test must have the definition test(integrator, abstol, reltol, min_t).

Keyword Arguments

• min_t specifies an optional minimum t before the steady state calculations are allowed to terminate.

IterativeCallback(time_choice, user_affect!, tType = Float64;
    initial_affect = false, kwargs...)

A callback to be used to iteratively apply some affect. For example, if given the first effect at t₁, you can define t₂ to apply the next effect.

Arguments

• time_choice(integrator) determines the time of the next callback. If nothing is returned for the time choice, then the iterator ends.
• user_affect! is the effect applied to the integrator at the stopping points.

Keyword Arguments

• initial_affect is whether to apply the affect at t=0 which defaults to false

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

PeriodicCallback can be used when a function should be called periodically in terms of integration time (as opposed to wall time), i.e. at t = tspan[1], t = tspan[1] + Δt, t = tspan[1] + 2Δt, and so on.

If a non-zero phase is provided, the invocations of the callback will be shifted by phase time units, i.e., the calls will occur at t = tspan[1] + phase, t = tspan[1] + phase + Δt, t = tspan[1] + phase + 2Δt, and so on.

This callback can, for example, be used to model a discrete-time controller for a continuous-time system, running at a fixed rate.

Arguments

• f the affect!(integrator) function to be called periodically
• Δt is the period

Keyword Arguments

• phase is a phase offset
• initial_affect is whether to apply the affect at the initial time, which defaults to false
• final_affect is whether to apply the affect at the final time, which defaults to false
• kwargs are keyword arguments accepted by the DiscreteCallback constructor.

PresetTimeCallback(tstops, user_affect!;
    initialize = DiffEqBase.INITIALIZE_DEFAULT,
    filter_tstops = true,
    kwargs...)

A callback that adds callback affect! calls at preset times. No playing around with tstops or anything is required: this callback adds the triggers for you to make it automatic.

Arguments

• tstops: the times for the affect! to trigger at.
• user_affect!: an affect!(integrator) function to use at the time points.

Keyword Arguments

• filter_tstops: Whether to filter out tstops beyond the end of the integration timespan. Defaults to true. If false, then tstops can extend the interval of integration.

AutoAbstol(save = true; init_curmax = 1e-6)

Provides a way to automatically adapt the absolute tolerance to the problem. This helps the solvers automatically “learn” what appropriate limits are. This callback set starts the absolute tolerance at init_curmax (default 1e-6), and at each iteration it is set to the maximum value that the state has thus far reached times the relative tolerance.

Keyword Arguments

• save determines whether this callback has saving enabled
• init_curmax is the initial abstol.

If this callback is used in isolation, save=true is required for normal saving behavior. Otherwise, save=false should be set to ensure extra saves do not occur.