GitHub

Numerical Integration Callbacks

Sometimes one may want to solve an integral simultaneously to the solution of a differential equation. For example, assume we want to solve:

\[u^\prime = f(u,p,t) h = \int_{t_0}^{t_f} g(u,p,t) dt\]

While one can use the ODE solver's dense solution to call an integration scheme on sol(t) after the solve, this can be memory intensive. Another way one can solve this problem is by extending the system, i.e.:

\[u^\prime = f(u,p,t) h^\prime = g(u,p,t)\]

with $h(t_0) = 0$, so then $h(t_f)$ would be the solution to the integral. However, many differential equation solvers scale superlinearly with the equation size and thus this could add an extra cost to the solver process.

The IntegratingCallback allows one to be able to solve such definite integrals in a way that is both memory and compute efficient. It uses the free local interpolation of a given step in order to approximate the Gaussian quadrature for a given step to the order of the numerical differential equation solve, thus achieving accuracy while not requiring the post-solution dense interpolation to be saved. By doing this via a callback, this method is able to easily integrate with functionality that introduces discontinuities, like other callbacks, in a way that is more accurate than a direct integration post solve.

The IntegratingSumCallback is the same, but instead of returning the timeseries of the interval results of the integration, it simply returns the final integral value.

The IntegratingGKCallback uses Gauss-Kronrod quadrature method in order to allow for error control.

The IntegratingGKSumCallback is the same, but instead of returning the timeseries of the interval results of the integration, it simply returns the final integral value.

DiffEqCallbacks.IntegratingCallback — Function
IntegratingCallback(integrand_func, integrand_values::IntegrandValues,
    integrand_prototype) -> DiscreteCallback

Construct a callback that uses fixed-order Gaussian quadrature to save the integral of integrand_func over each accepted solver step. The callback appends the step end time to integrand_values.ts and the corresponding integral estimate to integrand_values.integrand.

Arguments

  • integrand_func: for an out-of-place problem, define integrand_func(u, t, integrator) to return the integrand. For an in-place problem, define integrand_func(out, u, t, integrator) to write the integrand into out; when integrand_prototype === nothing, the allocating three-argument form is used instead. A returned value must not alias u.
  • integrand_values::IntegrandValues: storage for accepted-step end times and quadrature estimates. Its element types must accept the integration times and integrand outputs.
  • integrand_prototype: representative integrand output used to allocate the in-place output and accumulation buffers. Pass nothing to use the allocating form.

Returns

  • DiscreteCallback: a callback that saves one Gaussian quadrature estimate per accepted step.

Examples

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.

For the four-argument in-place form of integrand_func, pass an integrand_prototype that can be used as its output buffer.

source
DiffEqCallbacks.IntegrandValues — Type
IntegrandValues{tType, integrandType}

Storage used by IntegratingCallback and IntegratingGKCallback for one quadrature estimate per accepted solver step. Construct it before solving and pass the same instance to the callback; the callback appends to it in place.

Fields

  • ts::Vector{tType}: accepted-step end times associated with the stored estimates.
  • integrand::Vector{integrandType}: quadrature estimates over the corresponding solver steps. Entry integrand[i] approximates the integral from the previous saved solver time to ts[i].

Do not mutate either vector while a solve using the storage is active.

Constructors

  • IntegrandValues(tType, integrandType): create empty storage with time element type tType and estimate element type integrandType.

Examples

values = IntegrandValues(Float64, Float64)
source
DiffEqCallbacks.IntegratingSumCallback — Function
IntegratingSumCallback(integrand_func, integrand_values::IntegrandValuesSum,
    integrand_prototype; integrand_inplace = nothing) -> DiscreteCallback

Construct a callback that uses fixed-order Gaussian quadrature to accumulate the integral of integrand_func over accepted solver steps in integrand_values.integrand.

Arguments

  • integrand_func: define either integrand_func(u, t, integrator) to return the integrand or integrand_func(out, u, t, integrator) to write it into out. Returned or written values must be compatible with integrand_values.integrand.
  • integrand_values::IntegrandValuesSum: storage for the running integral. Construct it as IntegrandValuesSum(initial_value) with an initial value compatible with the integrand.
  • integrand_prototype: representative integrand output used as an in-place output buffer.

Keywords

  • integrand_inplace::Union{Nothing, Bool} = nothing: select the integrand calling form. With nothing, use the in-place form for an in-place problem when a cache can be allocated, and otherwise use the allocating form. Set this to true to force the in-place form or false to force the allocating form. true requires a non-nothingintegrand_prototype.

Returns

  • DiscreteCallback: a callback that adds each Gaussian quadrature estimate to integrand_values.integrand.

Throws

  • ArgumentError: if integrand_inplace = true and integrand_prototype === nothing.

Examples

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 solvers are required.

source
DiffEqCallbacks.IntegrandValuesSum — Type
IntegrandValuesSum{integrandType}

A struct used to save the accumulated integrand value in integrand::integrandType.

Fields

  • integrand::integrandType: the running quadrature total. The callback mutates this value in place for mutable values such as arrays; scalar totals are replaced with the updated scalar value.

Do not mutate integrand while a solve using this storage is active.

Constructors

  • IntegrandValuesSum(initial_value): Create with an initial value (recommended). The type is inferred from the value. Use zeros(n) for arrays or zero(T) for scalars.

Returns

An IntegrandValuesSum container for IntegratingSumCallback or IntegratingGKSumCallback.

Examples

# For array-valued integrands
integrated = IntegrandValuesSum(zeros(3))

# For scalar-valued integrands
integrated = IntegrandValuesSum(0.0)
source
DiffEqCallbacks.IntegratingGKCallback — Function
IntegratingGKCallback(integrand_func, integrand_values::IntegrandValues,
    integrand_prototype, tol = 1.0e-7) -> DiscreteCallback

Construct a callback that uses adaptive Gauss-Kronrod quadrature to save one integral estimate over each accepted solver step in integrand_values.integrand.

Arguments

  • integrand_func: for out-of-place problems, define integrand_func(u, t, integrator) to return the integrand. For in-place problems, define integrand_func(out, u, t, integrator) to write the integrand into out. Returned values must be compatible with integrand_values.
  • integrand_values::IntegrandValues: storage for accepted-step end times and quadrature estimates. Construct it as IntegrandValues(time_type, integrand_type) with an integrand_type that accepts the returned integrand values.
  • integrand_prototype: representative integrand output used to allocate the in-place cache.
  • tol::Real = 1.0e-7: absolute error tolerance for adaptive quadrature on each accepted solver step.

Returns

  • DiscreteCallback: a callback that appends one Gauss-Kronrod estimate after each accepted step.

Throws

  • An exception when used with an SDEProblem or RODEProblem, for which this Gauss-Kronrod algorithm is not guaranteed to converge.

Examples

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

This method is currently limited to ODE solvers of order 10 or lower.

For the four-argument in-place form of integrand_func, pass an integrand_prototype that can be used as its output buffer.

source
DiffEqCallbacks.IntegratingGKSumCallback — Function
IntegratingGKSumCallback(integrand_func, integrand_values::IntegrandValuesSum,
    integrand_prototype, tol = 1.0e-7; integrand_inplace = nothing) -> DiscreteCallback

Construct a callback that uses adaptive Gauss-Kronrod quadrature to accumulate the integral of integrand_func over accepted solver steps in integrand_values.integrand.

Arguments

  • integrand_func: define either integrand_func(u, t, integrator) to return the integrand or integrand_func(out, u, t, integrator) to write it into out. Returned or written values must be compatible with integrand_values.integrand.
  • integrand_values::IntegrandValuesSum: storage for the running integral. Construct it as IntegrandValuesSum(initial_value) with an initial value compatible with the integrand.
  • integrand_prototype: representative integrand output used as an in-place output buffer.
  • tol::Real = 1.0e-7: absolute error tolerance for adaptive quadrature on each accepted solver step.

Keywords

  • integrand_inplace::Union{Nothing, Bool} = nothing: select the integrand calling form. With nothing, use the in-place form for an in-place problem when a cache can be allocated, and otherwise use the allocating form. Set this to true to force the in-place form or false to force the allocating form. true requires a non-nothingintegrand_prototype.

Returns

  • DiscreteCallback: a callback that adds each Gauss-Kronrod estimate to integrand_values.integrand.

Throws

  • ArgumentError: if integrand_inplace = true and integrand_prototype === nothing.
  • An exception when used with an SDEProblem or RODEProblem, for which this Gauss-Kronrod algorithm is not guaranteed to converge.

Examples

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.

This method is currently limited to ODE solvers of order 10 or lower.

source

Example

using OrdinaryDiffEq, OrdinaryDiffEqLowOrderRK, DiffEqCallbacks, Test
prob = ODEProblem((u, p, t) -> [1.0], [0.0], (0.0, 1.0))
integrated = IntegrandValues(Float64, Vector{Float64})
sol = solve(prob, Euler(),
    callback = IntegratingCallback(
        (u, t, integrator) -> [1.0], integrated, Float64[0.0]),
    dt = 0.1)
@test all(integrated.integrand .≈ [[0.1] for i in 1:10])

integrated = IntegrandValues(Float64, Vector{Float64})
sol = solve(prob, Euler(),
    callback = IntegratingCallback(
        (u, t, integrator) -> [u[1]], integrated, Float64[0.0]),
    dt = 0.1)
@test all(integrated.integrand .≈ [[((n * 0.1)^2 - ((n - 1) * (0.1))^2) / 2] for n in 1:10])
@test sum(integrated.integrand)[1] ≈ 0.5

integrated = IntegrandValuesSum(zeros(1))
sol = solve(prob, Euler(),
    callback = IntegratingSumCallback(
        (u, t, integrator) -> [1.0], integrated, Float64[0.0]),
    dt = 0.1)
@test integrated.integrand[1] == 1
integrated = IntegrandValuesSum(zeros(1))
sol = solve(prob, Euler(),
    callback = IntegratingSumCallback(
        (u, t, integrator) -> [u[1]], integrated, Float64[0.0]),
    dt = 0.1)
@test integrated.integrand[1] == 0.5

integrated = IntegrandValues(Float64, Vector{Float64})
sol = solve(prob, Euler(),
    callback = IntegratingGKCallback(
        (u, t, integrator) -> [cos.(1000*u[1])], integrated, Float64[0.0], 1e-7),
    dt = 0.1)
@test sum(integrated.integrand)[1] .≈ sin(1000)/1000

integrated = IntegrandValuesSum(zeros(1))
sol = solve(prob, Euler(),
    callback = IntegratingGKSumCallback(
        (u, t, integrator) -> [cos.(1000*u[1])], integrated, Float64[0.0], 1e-7),
    dt = 0.1)
@test integrated.integrand[1] ≈ sin(1000)/1000
Test Passed