Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion src/DiffEqCallbacks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ using DataStructures: DataStructures, BinaryMaxHeap, BinaryMinHeap
using DiffEqBase: get_tstops, get_tstops_array, get_tstops_max
using DifferentiationInterface: DifferentiationInterface, Constant
using LinearAlgebra: LinearAlgebra, adjoint, axpy!, copyto!, diagind, mul!
using Markdown: @doc_str
using PrecompileTools: PrecompileTools
using RecipesBase: @recipe
using RecursiveArrayTools: RecursiveArrayTools, DiffEqArray, copyat_or_push!
Expand Down
17 changes: 17 additions & 0 deletions src/autoabstol.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,23 @@ is set to the maximum value that the state has thus far reached times the relati

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.

## Returns

A `DiscreteCallback` that updates the integrator absolute tolerance from the largest
state magnitude observed so far.

## Examples

```julia
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)
```
"""
function AutoAbstol(save = true; init_curmax = 0.0)
affect! = AutoAbstolAffect(abs.(init_curmax))
Expand Down
36 changes: 30 additions & 6 deletions src/domain.jl
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,20 @@ inside the domain. Thus, a `PositiveDomain` callback should generally be preferr
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

```julia
using 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)
sol = solve(prob, Tsit5(); callback = cb)
```
"""
function GeneralDomain(
g, u = nothing; save = true, abstol = nothing, scalefactor = nothing,
Expand Down Expand Up @@ -310,15 +324,13 @@ function GeneralDomain(
return CallbackSet(manifold_projection, domain_cb)
end

@doc doc"""
```julia
PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
```
"""
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
``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](https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html)
show.
Expand Down Expand Up @@ -346,7 +358,7 @@ 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)\}``,
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)``.

Expand All @@ -370,6 +382,18 @@ is the ``i``-th component of function ``f`` in an ODE system ``x' = f(x, t)``.
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

```julia
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)
```
"""
function PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
if u isa Nothing
Expand Down
36 changes: 28 additions & 8 deletions src/function_caller.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,23 +60,43 @@ function functioncalling_initialize(cb, u, t, integrator)
end

"""
```julia
FunctionCallingCallback(func;
funcat = Vector{Float64}(),
func_everystep = isempty(funcat),
func_start = true,
tdir = 1)
```
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:
which gets called at the time points of interest.

## Arguments

- `func(u, t, integrator)` is the function to be called.

## Keyword Arguments

- `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]`.

## Returns

A `DiscreteCallback` that calls `func` without modifying the integrator state.

## Examples

```julia
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)
```
"""
function FunctionCallingCallback(
func;
Expand Down
17 changes: 17 additions & 0 deletions src/integrating.jl
Original file line number Diff line number Diff line change
Expand Up @@ -235,6 +235,23 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
The outputted values are saved into `integrand_values`. The values are found
via `integrand_values.integrand`.

## Returns

A `DiscreteCallback` that saves one quadrature estimate per accepted step.

## Examples

```julia
using DiffEqCallbacks, OrdinaryDiffEq

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
values = IntegrandValues(Float64, Float64)
cb = IntegratingCallback((u, t, integrator) -> u^2, values, 0.0)

sol = solve(prob, Tsit5(); callback = cb)
integrals_by_step = values.integrand
```

!!! note

This method is currently limited to ODE solvers of order 10 or lower. Open an issue if other
Expand Down
18 changes: 18 additions & 0 deletions src/integrating_GK_affect.jl
Original file line number Diff line number Diff line change
Expand Up @@ -269,6 +269,24 @@ returns Integral(integrand_func(u(t),t)dt) over the problem tspan.
The outputted values are saved into `integrand_values`. The values are found
via `integrand_values.integrand`.

## Returns

A `DiscreteCallback` that saves one Gauss-Kronrod quadrature estimate per
accepted step.

## Examples

```julia
using DiffEqCallbacks, OrdinaryDiffEq

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
values = IntegrandValues(Float64, Float64)
cb = IntegratingGKCallback((u, t, integrator) -> u^2, values, 0.0)

sol = solve(prob, Tsit5(); callback = cb)
integrals_by_step = values.integrand
```

!!! note

Method has automatic error control (h-adaptive quadrature).
Expand Down
26 changes: 23 additions & 3 deletions src/integrating_GK_sum.jl
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,11 @@ end

"""
```julia
IntegratingCallback(integrand_func,
integrand_values::IntegrandValues,
cache = nothing)
IntegratingGKSumCallback(integrand_func,
integrand_values::IntegrandValuesSum,
integrand_prototype,
tol = 1.0e-7;
integrand_inplace = nothing)
```

Lets one define a function `integrand_func(u, t, integrator)` which
Expand Down Expand Up @@ -119,6 +121,24 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
The outputted values are saved into `integrand_values`. The values are found
via `integrand_values.integrand`.

## Returns

A `DiscreteCallback` that accumulates Gauss-Kronrod quadrature estimates into
`integrand_values.integrand`.

## Examples

```julia
using DiffEqCallbacks, OrdinaryDiffEq

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
values = IntegrandValuesSum(0.0)
cb = IntegratingGKSumCallback((u, t, integrator) -> u^2, values, 0.0)

sol = solve(prob, Tsit5(); callback = cb)
total = values.integrand
```

!!! note

This method uses Gauss-Kronrod quadrature rule to allow for error control.
Expand Down
25 changes: 22 additions & 3 deletions src/integrating_sum.jl
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,10 @@ end

"""
```julia
IntegratingCallback(integrand_func,
integrand_values::IntegrandValues,
cache = nothing)
IntegratingSumCallback(integrand_func,
integrand_values::IntegrandValuesSum,
integrand_prototype;
integrand_inplace = nothing)
```

Lets one define a function `integrand_func(u, t, integrator)` which
Expand Down Expand Up @@ -153,6 +154,24 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
The outputted values are saved into `integrand_values`. The values are found
via `integrand_values.integrand`.

## Returns

A `DiscreteCallback` that accumulates the quadrature estimate into
`integrand_values.integrand`.

## Examples

```julia
using DiffEqCallbacks, OrdinaryDiffEq

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
values = IntegrandValuesSum(0.0)
cb = IntegratingSumCallback((u, t, integrator) -> u^2, values, 0.0)

sol = solve(prob, Tsit5(); callback = cb)
total = values.integrand
```

!!! note

This method is currently limited to ODE solvers of order 10 or lower. Open an issue if other
Expand Down
53 changes: 44 additions & 9 deletions src/iterative_and_periodic.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
"""
```julia
IterativeCallback(time_choice, user_affect!, tType = Float64;
initial_affect = false, kwargs...)
```
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.
Expand All @@ -16,6 +14,28 @@ effect at `t₁`, you can define `t₂` to apply the next effect.
## Keyword Arguments

- `initial_affect` is whether to apply the affect at `t=0` which defaults to `false`
- `initialize` is the callback initialization function.
- `kwargs` are keyword arguments accepted by `DiscreteCallback`.

## Returns

A `DiscreteCallback` that repeatedly schedules the next time returned by
`time_choice`.

## Examples

```julia
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)
```
"""
function IterativeCallback(
time_choice, user_affect!, tType = Float64;
Expand Down Expand Up @@ -120,11 +140,8 @@ function add_next_tstop!(integrator, S::PeriodicCallbackAffect)
end

"""
```julia
PeriodicCallback(f, Δt::Number; phase = 0, initial_affect = false,
final_affect = false,
kwargs...)
```
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`,
Expand All @@ -149,6 +166,24 @@ discrete-time controller for a continuous-time system, running at a fixed rate.
- `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.

## Returns

A `DiscreteCallback` that schedules an `affect!` call every `Δt` units of
integration time.

## Examples

```julia
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)
```
"""
function PeriodicCallback(
f, Δt::Number;
Expand Down
20 changes: 20 additions & 0 deletions src/manifold.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,26 @@ order of the integrator even with the ManifoldProjection.
[1] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration:
Structure-Preserving Algorithms for Ordinary Differential Equations. Berlin ;
New York :Springer, 2002.

## Examples

```julia
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; residual_prototype = [0.0], autodiff = AutoFiniteDiff())

sol = solve(prob, Tsit5(); callback = cb)
```
"""
@concrete mutable struct ManifoldProjection
manifold
Expand Down
Loading
Loading