Integrator Interface

Integrator Interface

The integrator interface gives one the ability to interactively step through the numerical solving of a differential equation. Through this interface, one can easily monitor results, modify the problem during a run, and dynamically continue solving as one sees fit.

Initialization and Stepping

To initialize an integrator, use the syntax:

integrator = init(prob,alg;kwargs...)

The keyword args which are accepted are the same Common Solver Options used by solve and the returned value is an integrator which satisfies typeof(integrator)<:DEIntegrator. One can manually choose to step via the step! command:

step!(integrator)

which will take one successful step. Additonally:

step!(integrator,dt[,stop_at_tdt=false])

passing a dt will make the integrator keep stepping until integrator.t+dt, and setting stop_at_tdt=true will add a tstop to force it to step to integrator.t+dt

To check whether or not the integration step was successful, you can call check_error(integrator) which returns one of the Return Codes (RetCodes).

This type also implements an iterator interface, so one can step n times (or to the last tstop) using the take iterator:

for i in take(integrator,n) end

One can loop to the end by using solve!(integrator) or using the iterator interface:

for i in integrator end

In addition, some helper iterators are provided to help monitor the solution. For example, the tuples iterator lets you view the values:

for (u,t) in tuples(integrator)
  @show u,t
end

and the intervals iterator lets you view the full interval:

for (tprev,uprev,u,t) in intervals(integrator)
  @show tprev,t
end

Additionally, you can make the iterator return specific time points via the TimeChoiceIterator:

ts = linspace(0,1,11)
for (u,t) in TimeChoiceIterator(integrator,ts)
  @show u,t
end

Lastly, one can dynamically control the "endpoint". The initialization simply makes prob.tspan[2] the last value of tstop, and many of the iterators are made to stop at the final tstop value. However, step! will always take a step, and one can dynamically add new values of tstops by modifiying the variable in the options field: add_tstop!(integrator,new_t).

Finally, to solve to the last tstop, call solve!(integrator). Doing init and then solve! is equivalent to solve.

DiffEqBase.step!Function.
step!(integ::DEIntegrator [, dt [, stop_at_tdt]])

Perform one (successful) step on the integrator.

Alternative, if a dt is given, then step! the integrator until there is a temporal difference ≥ dt in integ.t. When true is passed to the optional third argument, the integrator advances exactly dt.

check_error(integrator)

Check state of integrator and return one of the Return Codes

check_error!(integrator)

Same as check_error but also set solution's return code (integrator.sol.retcode) and run postamble!.

Handing Integrators

The integrator<:DEIntegrator type holds all of the information for the intermediate solution of the differential equation. Useful fields are:

The p is the data which is provided by the user as a keyword arg in init. opts holds all of the common solver options, and can be mutated to change the solver characteristics. For example, to modify the absolute tolerance for the future timesteps, one can do:

integrator.opts.abstol = 1e-9

The sol field holds the current solution. This current solution includes the interpolation function if available, and thus integrator.sol(t) lets one interpolate efficiently over the whole current solution. Additionally, a a "current interval interpolation function" is provided on the integrator type via integrator(t). This uses only the solver information from the interval [tprev,t] to compute the interpolation, and is allowed to extrapolate beyond that interval.

Note about mutating

Be cautious: one should not directly mutate the t and u fields of the integrator. Doing so will destroy the accuracy of the interpolator and can harm certain algorithms. Instead if one wants to introduce discontinuous changes, one should use the Event Handling and Callback Functions. Modifications within a callback affect! surrounded by saves provides an error-free handling of the discontinuity.

As low-level alternative to the callbacks, one can use set_t!, set_u! and set_ut! to mutate integrator states. Note that certain integrators may not have efficient ways to modify u and t. In such case, set_*! are as inefficient as reinit!.

DiffEqBase.set_t!Function.
set_t!(integrator::DEIntegrator, t::Real)

Set current time point of the integrator to t.

DiffEqBase.set_u!Function.
set_u!(integrator::DEIntegrator, u)

Set current state of the integrator to u.

DiffEqBase.set_ut!Function.
set_ut!(integrator::DEIntegrator, u, t)

Set current state of the integrator to u and t

Integrator vs Solution

The integrator and the solution have very different actions because they have very different meanings. The typeof(sol) <: DESolution type is a type with history: it stores all of the (requested) timepoints and interpolates/acts using the values closest in time. On the other hand, the typeof(integrator)<:DEIntegrator type is a local object. It only knows the times of the interval it currently spans, the current caches and values, and the current state of the solver (the current options, tolerances, etc.). These serve very different purposes:

If one wants the solution object, then one can find it in integrator.sol.

Function Interface

In addition to the type interface, a function interface is provided which allows for safe modifications of the integrator type, and allows for uniform usage throughout the ecosystem (for packages/algorithms which implement the functions). The following functions make up the interface:

Saving Controls

Caches

Stepping Controls

Resizing

Reinit

The reinit function lets you restart the integration at a new value. The full function is of the form:

reinit!(integrator::ODEIntegrator,u0 = integrator.sol.prob.u0;
  t0 = integrator.sol.prob.tspan[1], tf = integrator.sol.prob.tspan[2],
  erase_sol = true,
  tstops = integrator.opts.tstops_cache,
  saveat = integrator.opts.saveat_cache,
  d_discontinuities = integrator.opts.d_discontinuities_cache,
  reset_dt = (integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive,
  reinit_callbacks = true, initialize_save = true,
  reinit_cache = true)

u0 is the value to start at. The starting time point and end point can be changed via t0 and tf. erase_sol allows one to start with no other values in the solution, or keep the previous solution. tstops, d_discontinuities, and saveat are reset as well, but can be ignored. reset_dt is a boolean for whether to reset the current value of dt using the automatic dt determination algorithm. reinit_callbacks is whether to run the callback initializations again (and initialize_save is for that). reinit_cache is whether to re-run the cache initialization function (i.e. resetting FSAL, not allocating vectors) which should usually be true for correctness.

Additionally, once can access auto_dt_reset!(integrator::ODEIntegrator) which will run the auto dt initialization algorithm.

Misc

Note

Note that not all of these functions will be implemented for every algorithm. Some have hard limitations. For example, Sundials.jl cannot resize problems. When a function is not limited, an error will be thrown.

Additional Options

The following options can additionally be specified in init (or be mutated in the opts) for further control of the integrator:

For example, if one wants to iterate but only stop at specific values, one can choose:

integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5],advance_to_tstop=true)
for (u,t) in tuples(integrator)
  @test t ∈ [0.5,1.0]
end

which will only enter the loop body at the values in tstops (here, prob.tspan[2]==1.0 and thus there are two values of tstops which are hit). Addtionally, one can solve! only to 0.5 via:

integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5])
integrator.opts.stop_at_next_tstop = true
solve!(integrator)

Plot Recipe

Like the DESolution type, a plot recipe is provided for the DEIntegrator type. Since the DEIntegrator type is a local state type on the current interval, plot(integrator) returns the solution on the current interval. The same options for the plot recipe are provided as for sol, meaning one can choose variables via the vars keyword argument, or change the plotdensity / turn on/off denseplot.

Additionally, since the integrator is an iterator, this can be used in the Plots.jl animate command to iteratively build an animation of the solution while solving the differential equation.

For an example of manually chaining together the iterator interface and plotting, one should try the following:

using DifferentialEquations, DiffEqProblemLibrary, Plots

# Linear ODE which starts at 0.5 and solves from t=0.0 to t=1.0
prob = ODEProblem((u,p,t)->1.01u,0.5,(0.0,1.0))

using Plots
integrator = init(prob,Tsit5();dt=1//2^(4),tstops=[0.5])
pyplot(show=true)
plot(integrator)
for i in integrator
  display(plot!(integrator,vars=(0,1),legend=false))
end
step!(integrator); plot!(integrator,vars=(0,1),legend=false)
savefig("iteratorplot.png")

Iterator Plot