Step Control Callbacks

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.

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.

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.

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.