Common Interface API

CommonSolve.init(args...; kwargs...) -> iter

Create an iterator or cache object that can be passed to CommonSolve.solve! or CommonSolve.step!. Generally, downstream packages extend:

iter = CommonSolve.init(prob::ProblemType, alg::SolverType; kwargs...)::IterType
CommonSolve.solve!(iter)::SolutionType

Arguments

• args...: Problem, algorithm, and implementation-specific positional arguments. The first positional argument must have a type owned by the package extending init.

Keywords

• kwargs...: Implementation-specific solver options.

Returns

An implementation-defined iterator or cache object that stores solver state.

Examples

struct MyProblem end
struct MyAlg end
struct MyIterator end

CommonSolve.init(::MyProblem, ::MyAlg; kwargs...) = MyIterator()

iter = CommonSolve.init(MyProblem(), MyAlg())

CommonSolve.solve!(iter) -> solution

Complete the solve using an iterator or cache object created by CommonSolve.init. Generally, downstream packages extend:

iter = CommonSolve.init(prob::ProblemType, alg::SolverType; kwargs...)::IterType
CommonSolve.solve!(iter)::SolutionType

Arguments

• iter: Solver state returned by CommonSolve.init. Its type must be owned by the package extending solve!.

Returns

The solution object defined by the downstream solver implementation.

Examples

struct MyIterator end

CommonSolve.solve!(::MyIterator) = :solution

CommonSolve.solve!(MyIterator())

CommonSolve.solve(args...; kwargs...) -> solution

Solve an equation or other mathematical problem using the algorithm specified in the arguments. Generally, downstream packages extend:

CommonSolve.solve(prob::ProblemType, alg::SolverType; kwargs...)::SolutionType

If a package only defines the iterator interface, solve falls back to:

solve(args...; kwargs...) = solve!(init(args...; kwargs...))

Arguments

• args...: Problem, algorithm, and implementation-specific positional arguments.

Keywords

• kwargs...: Implementation-specific solver options.

Interface

Extensions must dispatch the first positional argument on a type that they own. This prevents type piracy and ambiguities between independently developed solver packages.

Returns

The solution object defined by the downstream solver implementation.

Examples

struct MyProblem end
struct MyAlg end

CommonSolve.solve(::MyProblem, ::MyAlg; kwargs...) = :solution

CommonSolve.solve(MyProblem(), MyAlg())

CommonSolve.step!(iter, args...; kwargs...) -> step_result

Progress an iterator or cache object returned by CommonSolve.init. The additional arguments typically describe how far to advance the solve and are implementation-specific.

Arguments

• iter: Solver state returned by CommonSolve.init. Its type must be owned by the package extending step!.
• args...: Implementation-specific step controls.

Keywords

• kwargs...: Implementation-specific step options.

Returns

An implementation-defined value, commonly the updated iterator, a step result, or nothing.

Examples

mutable struct MyIterator
    steps::Int
end

function CommonSolve.step!(iter::MyIterator)
    iter.steps += 1
    return iter
end

iter = CommonSolve.step!(MyIterator(0))

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

Advance a differential equation integrator.

With one argument, perform one accepted solver step according to the concrete algorithm. With dt, repeatedly step until the signed time displacement from the starting time is at least dt. When stop_at_tdt is true, the generic fallback adds a temporary tstop so the integrator lands exactly at t + dt. Negative stepping relative to integ.tdir is rejected by the fallback.

struct DAEFunction{iip, specialize, F, Ta, Tt, TJ, TJU, TJD, JVP, VJP, JP, SP, TW, TWt, TPJ, O, TCV, SYS, ID, NLP<:Union{Nothing, SciMLBase.ODENLStepData}} <: SciMLBase.AbstractDAEFunction{iip}

A representation of an implicit DAE function f, defined by:

\[0 = f{\left(\frac{du}{dt},u,p,t\right)}\]

and all of its related functions, such as the Jacobian of f, its gradient with respect to time, and more. For all cases, u0 is the initial condition, p are the parameters, and t is the independent variable.

Constructor

DAEFunction{iip,specialize}(f;
                           analytic = __has_analytic(f) ? f.analytic : nothing,
                           jac = __has_jac(f) ? f.jac : nothing,
                           jac_u = __has_jac_u(f) ? f.jac_u : nothing,
                           jac_du = __has_jac_du(f) ? f.jac_du : nothing,
                           jvp = __has_jvp(f) ? f.jvp : nothing,
                           vjp = __has_vjp(f) ? f.vjp : nothing,
                           jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
                           sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
                           colorvec = __has_colorvec(f) ? f.colorvec : nothing,
                           sys = __has_sys(f) ? f.sys : nothing,
                           nlstep_data = __has_nlstep_data(f) ? f.nlstep_data : nothing)

Note that only the function f itself is required. This function should be given as f!(out,du,u,p,t) or out = f(du,u,p,t). See the section on iip for more details on in-place vs out-of-place handling.

All of the remaining functions are optional for improving or accelerating the usage of f. These include:

• analytic(u0,p,t): used to pass an analytical solution function for the analytical solution of the ODE. Generally only used for testing and development of the solvers.
• jac(J,du,u,p,gamma,t) or J=jac(du,u,p,gamma,t): returns the implicit DAE Jacobian defined as $γ \frac{dG}{d(du)} + \frac{dG}{du}$
• jac_u(J,du,u,p,t) or J=jac_u(du,u,p,t): returns the partial DAE Jacobian $\frac{dG}{du}$
• jac_du(J,du,u,p,t) or J=jac_du(du,u,p,t): returns the partial DAE Jacobian $\frac{dG}{d(du)}$ When jac_u and jac_du are provided, the solver can efficiently reuse them when only the coupling coefficient γ changes (e.g., step size or order changes), avoiding full Jacobian recomputation. If only jac is provided, the solver extracts the components automatically.
• jvp(Jv,v,du,u,p,gamma,t) or Jv=jvp(v,du,u,p,gamma,t): returns the directional derivative$\frac{df}{du} v$
• vjp(Jv,v,du,u,p,gamma,t) or Jv=vjp(v,du,u,p,gamma,t): returns the adjoint derivative$\frac{df}{du}^∗ v$
• jac_prototype: a prototype matrix matching the type that matches the Jacobian. For example, if the Jacobian is tridiagonal, then an appropriately sized Tridiagonal matrix can be used as the prototype and integrators will specialize on this structure where possible. Non-structured sparsity patterns should use a SparseMatrixCSC with a correct sparsity pattern for the Jacobian. The default is nothing, which means a dense Jacobian.
• colorvec: a color vector according to the SparseDiffTools.jl definition for the sparsity pattern of the jac_prototype. This specializes the Jacobian construction when using finite differences and automatic differentiation to be computed in an accelerated manner based on the sparsity pattern. Defaults to nothing, which means a color vector will be internally computed on demand when required. The cost of this operation is highly dependent on the sparsity pattern.
• nlstep_data: an ODENLStepData holding a structured nonlinear problem for the implicit stage solve, or nothing. Implicit DAE integrators which support it solve this problem in place of building a stage-equation closure. See the ODENLStepData documentation for the stage equation the nonlinear problem must represent in the fully implicit case.

iip: In-Place vs Out-Of-Place

For more details on this argument, see the ODEFunction documentation.

specialize: Controlling Compilation and Specialization

For more details on this argument, see the ODEFunction documentation.

Fields

The fields of the DAEFunction type directly match the names of the inputs.

Examples

Declaring Explicit Jacobians for DAEs

For fully implicit ODEs (DAEProblems), a slightly different Jacobian function is necessary. For the DAE

\[G(du,u,p,t) = \text{res}\]

The Jacobian should be given in the form gamma*dG/d(du) + dG/du where gamma is given by the solver. This means that the signature is:

f(J,du,u,p,gamma,t)

For example, for the equation

function testjac(res,du,u,p,t)
  res[1] = du[1] - 2.0 * u[1] + 1.2 * u[1]*u[2]
  res[2] = du[2] -3 * u[2] - u[1]*u[2]
end

we would define the Jacobian as:

function testjac(J,du,u,p,gamma,t)
  J[1,1] = gamma - 2.0 + 1.2 * u[2]
  J[1,2] = 1.2 * u[1]
  J[2,1] = - 1 * u[2]
  J[2,2] = gamma - 3 - u[1]
  nothing
end

Symbolically Generating the Functions

See the modelingtoolkitize function from ModelingToolkit.jl for automatically symbolically generating the Jacobian and more from the numerically-defined functions.

Defines an implicit ordinary differential equation (ODE) or differential-algebraic equation (DAE) problem. Documentation Page: https://docs.sciml.ai/DiffEqDocs/stable/types/dae_types/

Mathematical Specification of an DAE Problem

To define a DAE Problem, you simply need to give the function $f$ and the initial condition $u_0$ which define an ODE:

\[0 = f(du,u,p,t)\]

f should be specified as f(du,u,p,t) (or in-place as f(resid,du,u,p,t)). Note that we are not limited to numbers or vectors for u₀; one is allowed to provide u₀ as arbitrary matrices / higher dimension tensors as well.

Problem Type

Constructors

• DAEProblem(f::DAEFunction,du0,u0,tspan,p=NullParameters();kwargs...)
• DAEProblem{isinplace,specialize}(f,du0,u0,tspan,p=NullParameters();kwargs...) : Defines the DAE with the specified functions. isinplace optionally sets whether the function is inplace or not. This is determined automatically, but not inferred. specialize optionally controls the specialization level. See the specialization levels section of the documentation for more details. The default is AutoSpecialize.

For more details on the in-place and specialization controls, see the ODEFunction documentation.

Parameters are optional, and if not given, then a NullParameters() singleton will be used which will throw nice errors if you try to index non-existent parameters. Any extra keyword arguments are passed on to the solvers. For example, if you set a callback in the problem, then that callback will be added in every solve call.

For specifying Jacobians and mass matrices, see the SciMLFunctions interface page.

Fields

• f: The function in the ODE.
• du0: The initial condition for the derivative.
• u0: The initial condition.
• tspan: The timespan for the problem.
• differential_vars: A logical array which declares which variables are the differential (non-algebraic) vars (i.e. du' is in the equations for this variable). Defaults to nothing. Some solvers may require this be set if an initial condition needs to be determined.
• p: The parameters for the problem. Defaults to NullParameters
• kwargs: The keyword arguments passed onto the solves.

Example Problems

Examples problems can be found in DiffEqProblemLibrary.jl.

To use a sample problem, such as prob_dae_resrob, you can do something like:

#] add DAEProblemLibrary
using DAEProblemLibrary
prob = DAEProblemLibrary.prob_dae_resrob
sol = solve(prob,IDA())

struct DAESolution{T, N, uType, duType, uType2, DType, tType, P, A, ID, S, rateType, V} <: SciMLBase.AbstractDAESolution{T, N, uType}

Representation of the solution to an differential-algebraic equation defined by an DAEProblem.

DESolution Interface

For more information on interacting with DESolution types, check out the Solution Handling page of the DifferentialEquations.jl documentation.

https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/

Fields

• u: the representation of the DAE solution. Given as an array of solutions, where u[i] corresponds to the solution at time t[i]. It is recommended in most cases one does not access sol.u directly and instead use the array interface described in the Solution Handling page of the DifferentialEquations.jl documentation.
• du: the representation of the derivatives of the DAE solution.
• t: the time points corresponding to the saved values of the DAE solution.
• prob: the original DAEProblem that was solved.
• alg: the algorithm type used by the solver.
• stats: statistics of the solver, such as the number of function evaluations required, number of Jacobians computed, and more.
• retcode: the return code from the solver. Used to determine whether the solver solved successfully, whether it terminated early due to a user-defined callback, or whether it exited due to an error. For more details, see the return code documentation.

struct DynamicalODEFunction{iip, specialize, F1, F2, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, O, TCV, SYS, ID} <: SciMLBase.AbstractODEFunction{iip}

A representation of an ODE function f, defined by:

\[M \frac{du}{dt} = f(u,p,t)\]

as a partitioned ODE:

\[\begin{align*} M_1 \frac{du}{dt} = f_1(u,p,t) \\ M_2 \frac{du}{dt} = f_2(u,p,t) \end{align*}\]

and all of its related functions, such as the Jacobian of f, its gradient with respect to time, and more. For all cases, u0 is the initial condition, p are the parameters, and t is the independent variable.

Constructor

DynamicalODEFunction{iip,specialize}(f1,f2;
                                    mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
                                    analytic = __has_analytic(f) ? f.analytic : nothing,
                                    tgrad= __has_tgrad(f) ? f.tgrad : nothing,
                                    jac = __has_jac(f) ? f.jac : nothing,
                                    jvp = __has_jvp(f) ? f.jvp : nothing,
                                    vjp = __has_vjp(f) ? f.vjp : nothing,
                                    jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
                                    sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
                                    paramjac = __has_paramjac(f) ? f.paramjac : nothing,
                                    colorvec = __has_colorvec(f) ? f.colorvec : nothing,
                                    sys = __has_sys(f) ? f.sys : nothing)

Note that only the functions f_i themselves are required. These functions should be given as f_i!(du,u,p,t) or du = f_i(u,p,t). See the section on iip for more details on in-place vs out-of-place handling.

All of the remaining functions are optional for improving or accelerating the usage of f. These include:

• mass_matrix: the mass matrix M_i represented in the ODE function. Can be used to determine that the equation is actually a differential-algebraic equation (DAE) if M is singular. Note that in this case special solvers are required, see the DAE solver page for more details: https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/. Must be an AbstractArray or an AbstractSciMLOperator. Should be given as a tuple of mass matrices, i.e. (M_1, M_2) for the mass matrices of equations 1 and 2 respectively.
• analytic(u0,p,t): used to pass an analytical solution function for the analytical solution of the ODE. Generally only used for testing and development of the solvers.
• tgrad(dT,u,p,t) or dT=tgrad(u,p,t): returns $\frac{∂f(u,p,t)}{∂t}$
• jac(J,u,p,t) or J=jac(u,p,t): returns $\frac{df}{du}$
• jvp(Jv,v,u,p,t) or Jv=jvp(v,u,p,t): returns the directional derivative $\frac{df}{du} v$
• vjp(Jv,v,u,p,t) or Jv=vjp(v,u,p,t): returns the adjoint derivative $\frac{df}{du}^∗ v$
• jac_prototype: a prototype matrix matching the type that matches the Jacobian. For example, if the Jacobian is tridiagonal, then an appropriately sized Tridiagonal matrix can be used as the prototype and integrators will specialize on this structure where possible. Non-structured sparsity patterns should use a SparseMatrixCSC with a correct sparsity pattern for the Jacobian. The default is nothing, which means a dense Jacobian.
• paramjac(pJ,u,p,t): returns the parameter Jacobian $\frac{df}{dp}$.
• colorvec: a color vector according to the SparseDiffTools.jl definition for the sparsity pattern of the jac_prototype. This specializes the Jacobian construction when using finite differences and automatic differentiation to be computed in an accelerated manner based on the sparsity pattern. Defaults to nothing, which means a color vector will be internally computed on demand when required. The cost of this operation is highly dependent on the sparsity pattern.

iip: In-Place vs Out-Of-Place

For more details on this argument, see the ODEFunction documentation.

specialize: Controlling Compilation and Specialization

For more details on this argument, see the ODEFunction documentation.

Fields

The fields of the DynamicalODEFunction type directly match the names of the inputs.

DynamicalODEProblem(f::DynamicalODEFunction, v0, u0, tspan, p = NullParameters(), callback = CallbackSet())

Define a dynamical ODE function from a DynamicalODEFunction.

Defines a dynamical ordinary differential equation (ODE) problem. Documentation Page: https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/

Dynamical ordinary differential equations, such as those arising from the definition of a Hamiltonian system or a second order ODE, have a special structure that can be utilized in the solution of the differential equation. On this page, we describe how to define second order differential equations for their efficient numerical solution.

Mathematical Specification of a Dynamical ODE Problem

These algorithms require a Partitioned ODE of the form:

\[\begin{align*} \frac{dv}{dt} &= f_1(u,t) \\ \frac{du}{dt} &= f_2(v) \\ \end{align*}\]

This is a Partitioned ODE partitioned into two groups, so the functions should be specified as f1(dv,v,u,p,t) and f2(du,v,u,p,t) (in the inplace form), where f1 is independent of v (unless specified by the solver), and f2 is independent of u and t. This includes discretizations arising from SecondOrderODEProblems where the velocity is not used in the acceleration function, and Hamiltonians where the potential is (or can be) time-dependent, but the kinetic energy is only dependent on v.

Note that some methods assume that the integral of f2 is a quadratic form. That means that f2=v'*M*v, i.e. $\int f_2 = \frac{1}{2} m v^2$, giving du = v. This is equivalent to saying that the kinetic energy is related to $v^2$. The methods which require this assumption will lose accuracy if this assumption is violated. Methods listed make note of this requirement with "Requires quadratic kinetic energy".

Constructor

DynamicalODEProblem(f::DynamicalODEFunction,v0,u0,tspan,p=NullParameters();kwargs...)
DynamicalODEProblem{isinplace}(f1,f2,v0,u0,tspan,p=NullParameters();kwargs...)

Defines the ODE with the specified functions. isinplace optionally sets whether the function is inplace or not. This is determined automatically, but not inferred.

Parameters are optional, and if not given, then a NullParameters() singleton will be used which will throw nice errors if you try to index non-existent parameters. Any extra keyword arguments are passed on to the solvers. For example, if you set a callback in the problem, then that callback will be added in every solve call.

Fields

• f1 and f2: The functions in the ODE.
• v0 and u0: The initial conditions.
• tspan: The timespan for the problem.
• p: The parameters for the problem. Defaults to NullParameters
• kwargs: The keyword arguments passed onto the solves.

EnsembleContext{S, R, M}

Contextual information about the current trajectory within an ensemble solve.

An EnsembleContext is passed to rng_func, prob_func(prob, ctx), and output_func(sol, ctx). It is the stable interface for selecting trajectory-specific data without relying on global counters or backend-specific worker state.

Fields

• sim_id::Int: Unique trajectory index in 1:trajectories.
• repeat::Int: Rerun counter, starting at 1 and incremented when output_func requests a rerun.
• worker_id::Int: 0 for serial and threaded execution, or Distributed.myid() for distributed execution.
• sim_seed::S: Pre-generated seed for this trajectory, or nothing.
• rng::R: Per-trajectory RNG created by rng_func, or nothing while rng_func itself is running.
• master_rng::M: User-provided master RNG, or nothing. Distributed modes set this to nothing to avoid serializing mutable RNG state to workers.

struct EnsembleProblem{T, T2, T3, T4, T5} <: SciMLBase.AbstractEnsembleProblem

Container for a template problem and the user hooks used to run an ensemble of related SciML solves.

An EnsembleProblem is solved by repeatedly calling prob_func(prob, ctx) to construct the trajectory-specific problem, solving that problem with the requested numerical algorithm, passing the result through output_func(sol, ctx), and combining batches with reduction(u, data, I). The ctx argument is an EnsembleContext that identifies the trajectory and carries per-trajectory RNG state when rng or seed is supplied to solve.

Constructor

EnsembleProblem(prob::AbstractSciMLProblem;
    output_func = (sol, ctx) -> (sol, false),
    prob_func = (prob, ctx) -> prob,
    reduction = (u, data, I) -> (append!(u, data), false),
    u_init = [], safetycopy = prob_func !== DEFAULT_PROB_FUNC)

Positional Arguments

• prob: The canonical problem used as the template for each trajectory.

Keyword Arguments

• prob_func: A function (prob, ctx) that modifies the problem for each trajectory. ctx is an EnsembleContext providing ctx.sim_id (unique id 1:trajectories), ctx.repeat (rerun counter, starts at 1), ctx.rng (per-trajectory RNG or nothing), ctx.sim_seed, and ctx.master_rng. prob_func must preserve the problem type of prob; for example, a JumpProblem must remain a JumpProblem, an ODEProblem must remain an ODEProblem.
• output_func: A function (sol, ctx) that determines what is saved from each trajectory. It returns (out, rerun), where out is stored in the batch output and rerun requests that the same trajectory be rerun with ctx.repeat incremented.
• reduction: A function (u, data, I) that combines the current accumulator u with the outputs data from the trajectory index range I. It returns (new_u, converged), where converged=true stops the ensemble early.
• u_init: The initial accumulator passed to reduction. When nothing, the accumulator is initialized from the first batch output.
• safetycopy: Determines whether a safety deepcopy is called on the prob before the prob_func. By default, this is true for any user-given prob_func, as without this, modifying the arguments of something in the prob_func, such as parameters or caches stored within the user function, are not necessarily thread-safe. If you know that your function is thread-safe, then setting this to false can improve performance when used with threads. For nested problems, e.g., SDE problems with custom noise processes, deepcopy might be insufficient. In such cases, use a custom prob_func.

Example

function prob_func(prob, ctx)
    remake(prob, u0 = randn(ctx.rng, length(prob.u0)))
end

output_func(sol, ctx) = (sol[end, 2], false)
ensemble_prob = EnsembleProblem(prob; prob_func, output_func)

EnsembleProblem(
    prob;
    prob_func,
    output_func,
    reduction,
    u_init,
    safetycopy
)

Construct an EnsembleProblem from a template problem and optional trajectory, output, and reduction hooks.

User-supplied hooks are passed through prepare_function, and u_init is normalized through prepare_initial_state. The default safetycopy is true when a custom prob_func is supplied, since mutating shared objects inside prob_func is otherwise not thread-safe in threaded ensemble modes.

EnsembleProblem(
;
    prob,
    prob_func,
    output_func,
    reduction,
    u_init,
    p,
    safetycopy
)

Keyword-only constructor for EnsembleProblem. This is equivalent to EnsembleProblem(prob; kwargs...) with prob supplied as a keyword.

EnsembleProblem(
    prob;
    prob_func,
    output_func,
    reduction,
    u_init,
    safetycopy
)
EnsembleProblem(prob, u0s; kwargs...)

Deprecated constructor that builds an ensemble by selecting initial conditions from u0s by trajectory index.

struct ODEFunction{iip, specialize, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, WP, TPJ, VJP_P, O, TCV, SYS, ID<:Union{Nothing, SciMLBase.OverrideInitData}, NLP<:Union{Nothing, SciMLBase.ODENLStepData}} <: SciMLBase.AbstractODEFunction{iip}

A representation of an ODE function f, defined by:

\[M \frac{du}{dt} = f(u,p,t)\]

and all of its related functions, such as the Jacobian of f, its gradient with respect to time, and more. For all cases, u0 is the initial condition, p are the parameters, and t is the independent variable.

Constructor

ODEFunction{iip,specialize}(f;
                           mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
                           analytic = __has_analytic(f) ? f.analytic : nothing,
                           tgrad= __has_tgrad(f) ? f.tgrad : nothing,
                           jac = __has_jac(f) ? f.jac : nothing,
                           jvp = __has_jvp(f) ? f.jvp : nothing,
                           vjp = __has_vjp(f) ? f.vjp : nothing,
                           jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
                           sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
                           paramjac = __has_paramjac(f) ? f.paramjac : nothing,
                           vjp_p = __has_vjp_p(f) ? f.vjp_p : nothing,
                           colorvec = __has_colorvec(f) ? f.colorvec : nothing,
                           sys = __has_sys(f) ? f.sys : nothing)

Note that only the function f itself is required. This function should be given as f!(du,u,p,t) or du = f(u,p,t). See the section on iip for more details on in-place vs out-of-place handling.

All of the remaining functions are optional for improving or accelerating the usage of f. These include:

• mass_matrix: the mass matrix M represented in the ODE function. Can be used to determine that the equation is actually a differential-algebraic equation (DAE) if M is singular. Note that in this case special solvers are required, see the DAE solver page for more details: https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/. Must be an AbstractArray or an AbstractSciMLOperator.
• analytic(u0,p,t): used to pass an analytical solution function for the analytical solution of the ODE. Generally only used for testing and development of the solvers.
• tgrad(dT,u,p,t) or dT=tgrad(u,p,t): returns $\frac{∂f(u,p,t)}{∂t}$
• jac(J,u,p,t) or J=jac(u,p,t): returns $\frac{df}{du}$
• jvp(Jv,v,u,p,t) or Jv=jvp(v,u,p,t): returns the directional derivative$\frac{df}{du} v$
• vjp(Jv,v,u,p,t) or Jv=vjp(v,u,p,t): returns the adjoint derivative$\frac{df}{du}^∗ v$
• jac_prototype: a prototype matrix matching the type that matches the Jacobian. For example, if the Jacobian is tridiagonal, then an appropriately sized Tridiagonal matrix can be used as the prototype and integrators will specialize on this structure where possible. Non-structured sparsity patterns should use a SparseMatrixCSC with a correct sparsity pattern for the Jacobian. It must support the operations required by the selected differentiation and linear solver. The default is nothing, which means a dense Jacobian. Solvers may copy, allocate a similar object, or convert the prototype, so callers must not rely on object identity or aliasing.
• paramjac(pJ,u,p,t): returns the parameter Jacobian $\frac{df}{dp}$.
• vjp_p(Jpv,v,u,p,t) or Jpv=vjp_p(v,u,p,t): returns the parameter adjoint derivative $\frac{df}{dp}^∗ v$, i.e. the vector-Jacobian product with respect to parameters. This avoids materializing the full parameter Jacobian when only the VJP is needed (e.g. in adjoint sensitivity analysis). When not provided, falls back to paramjac or AD-based computation.
• colorvec: a column-color vector compatible with the selected sparse differentiation backend and the sparsity pattern of jac_prototype. This can accelerate Jacobian construction with finite differences or automatic differentiation. Defaults to nothing, which lets the selected backend compute coloring when required.

iip: In-Place vs Out-Of-Place

iip is the optional boolean for determining whether a given function is written to be used in-place or out-of-place. In-place functions are f!(du,u,p,t) where the return is ignored, and the result is expected to be mutated into the value of du. Out-of-place functions are du=f(u,p,t).

Normally, this is determined automatically by looking at the method table for f and seeing the maximum number of arguments in available dispatches. For this reason, the constructor ODEFunction(f) generally works (but is type-unstable). However, for type-stability or to enforce correctness, this option is passed via ODEFunction{true}(f).

specialize: Controlling Compilation and Specialization

The specialize parameter controls the specialization level of the ODEFunction on the function f. This allows for a trade-off between compile and run time performance. The available specialization levels are:

• SciMLBase.AutoSpecialize: this form performs a lazy function wrapping on the functions of the ODE in order to stop recompilation of the ODE solver, but allow for the prob.f to stay unwrapped for normal usage. This is the default specialization level and strikes a balance in compile time vs runtime performance.
• SciMLBase.FullSpecialize: this form fully specializes the ODEFunction on the constituent functions that make its fields. As such, each ODEFunction in this form is uniquely typed, requiring re-specialization and compilation for each new ODE definition. This form has the highest compile-time at the cost of being the most optimal in runtime. This form should be preferred for long-running calculations (such as within optimization loops) and for benchmarking.
• SciMLBase.NoSpecialize: this form fully unspecializes the function types in the ODEFunction definition by using an Any type declaration. As a result, it can result in reduced runtime performance, but is the form that induces the least compile-time.
• SciMLBase.FunctionWrapperSpecialize: this is an eager function wrapping form. It is unsafe with many solvers, and thus is mostly used for development testing.

For more details, see Specialization Levels.

Fields

The fields of the ODEFunction type directly match the names of the inputs.

More Details on Jacobians

The following example creates an inplace ODEFunction whose Jacobian is a Diagonal:

using LinearAlgebra
f = (du,u,p,t) -> du .= t .* u
jac = (J,u,p,t) -> (J[1,1] = t; J[2,2] = t; J)
jp = Diagonal(zeros(2))
fun = ODEFunction(f; jac=jac, jac_prototype=jp)

The prototype declares Jacobian shape, element type, and structure. It must support the operations required by the selected differentiation and linear solver, including writes for an in-place jac and any multiplication, diagonal-shift, or factorization operations that solver performs. Solvers may copy, allocate with similar, or convert the prototype, so code must not rely on its object identity or aliasing behavior.

When jac_prototype is an AbstractSciMLOperator and jac is omitted, the constructor creates a Jacobian update using update_coefficients! for the in-place form and update_coefficients for the out-of-place form. Refer to the AbstractSciMLOperators documentation for more information on setting up time/parameter dependent operators.

Examples

Declaring Explicit Jacobians for ODEs

The most standard case, declaring a function for a Jacobian is done by overloading the function f(du,u,p,t) with an in-place updating function for the Jacobian: f_jac(J,u,p,t) where the value type is used for dispatch. For example, take the Lotka-Volterra model:

function f(du,u,p,t)
  du[1] = 2.0 * u[1] - 1.2 * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

To declare the Jacobian, we simply add the dispatch:

function f_jac(J,u,p,t)
  J[1,1] = 2.0 - 1.2 * u[2]
  J[1,2] = -1.2 * u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

Then we can supply the Jacobian with our ODE as:

ff = ODEFunction(f;jac=f_jac)

and use this in an ODEProblem:

prob = ODEProblem(ff,ones(2),(0.0,10.0))

Symbolically Generating the Functions

See the modelingtoolkitize function from ModelingToolkit.jl for automatically symbolically generating the Jacobian and more from the numerically-defined functions.

ODEFunction(f)

Converts a NonlinearFunction into an ODEFunction.

ODENLStepData(nlprob, u0perm, set_gamma_c, set_outer_tmp, set_inner_tmp, nlprobmap)

A collection of hooks for custom nonlinear stage solves in implicit ODE and DAE algorithms.

ODENLStepData lets an ODEFunction, SplitFunction or DAEFunction provide a structured AbstractNonlinearProblem template for solver packages that form implicit stage equations. Before each nonlinear solve, the algorithm updates the stage guess, scaling factors, time information, and temporary vectors through the stored setter callables. After the nonlinear solve, nlprobmap converts the nonlinear unknown back to the state vector used by the original problem.

Mass-matrix form

For M * du/dt = f(u, p, t) the nonlinear problem should represent a stage equation of the form M * z = outer_tmp + gamma1 * f(gamma2 * z + inner_tmp, p, t_c), equivalently g(z, p') = gamma1 * f(gamma2 * z + inner_tmp, p, t_c) + outer_tmp - M * z. Here z is the nonlinear stage unknown, p is the ODE parameter object, t_c is the stage evaluation time, and gamma1, gamma2, outer_tmp, and inner_tmp are supplied by the ODE algorithm.

Fully implicit form

For 0 = F(du, u, p, t) (a DAEFunction) the stage equation has the same shape, with both arguments of F affine in the stage unknown: g(z, p') = F(gamma1 * z + outer_tmp, gamma2 * z + inner_tmp, p, t_c). gamma2 and inner_tmp build the state argument from the stage unknown exactly as in the mass-matrix form, while gamma1 and outer_tmp build the derivative argument. Taking the stage unknown to be the stage state (gamma2 = 1, inner_tmp = 0), a BDF-type step with du ≈ (u - tmp) / (γ * dt) gives gamma1 = inv(γ * dt) and outer_tmp = -tmp / (γ * dt). With that convention gamma1 is the gamma of the DAEFunction Jacobian signature jac(J, du, u, p, gamma, t): the Jacobian of the stage residual with respect to z is gamma1 * dF/d(du) + dF/du.

Fields

• nlprob::Any: The structured AbstractNonlinearProblem template solved for each implicit ODE stage.

• u0perm::Any: Callable used by the ODE algorithm to update the nonlinear problem's initial guess from the current stage data.

• set_γ_c::Any: Callable used by the ODE algorithm to update the stage scaling factors and stage time/abscissa data used by the nonlinear problem.

• set_outer_tmp::Any: Callable used by the ODE algorithm to update the outer_tmp vector in the nonlinear stage equation.

• set_inner_tmp::Any: Callable used by the ODE algorithm to update the inner_tmp vector in the nonlinear stage equation.

• nlprobmap::Any: Callable that maps the solution of nlprob back to the state vector or stage vector of the original ODE problem.

Extension Rules

Symbolic-system packages construct this value and store it as the nlstep_data of an ODEFunction, SplitFunction or DAEFunction. Solver packages may consume the six fields through their callable contracts, but must not assume concrete callable types or mutate the container. Each setter must update the object it closes over consistently with nlprob, and nlprobmap must map a completed nonlinear solution back to the stage representation of the original problem.

Defines an ordinary differential equation (ODE) problem. Documentation Page: https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/

Mathematical Specification of an ODE Problem

To define an ODE Problem, you simply need to give the function $f$ and the initial condition $u_0$ which define an ODE:

\[M \frac{du}{dt} = f(u,p,t)\]

There are two different ways of specifying f:

• f(du,u,p,t): in-place. Memory-efficient when avoiding allocations. Best option for most cases unless mutation is not allowed.
• f(u,p,t): returning du. Less memory-efficient way, particularly suitable when mutation is not allowed (e.g. with certain automatic differentiation packages such as Zygote).

u₀ should be an AbstractArray (or number) whose geometry matches the desired geometry of u. Note that we are not limited to numbers or vectors for u₀; one is allowed to provide u₀ as arbitrary matrices / higher dimension tensors as well.

For the mass matrix $M$, see the documentation of ODEFunction.

Problem Type

Constructors

ODEProblem can be constructed by first building an ODEFunction or by simply passing the ODE right-hand side to the constructor. The constructors are:

• ODEProblem(f::ODEFunction,u0,tspan,p=NullParameters();kwargs...)
• ODEProblem{isinplace,specialize}(f,u0,tspan,p=NullParameters();kwargs...) : Defines the ODE with the specified functions. isinplace optionally sets whether the function is inplace or not. This is determined automatically, but not inferred. specialize optionally controls the specialization level. See the Specialization Levels for more details. The default is AutoSpecialize.

For more details on the in-place and specialization controls, see the ODEFunction documentation.

Parameters are optional, and if not given, then a NullParameters() singleton will be used which will throw nice errors if you try to index non-existent parameters. Any extra keyword arguments are passed on to the solvers. For example, if you set a callback in the problem, then that callback will be added in every solve call.

For specifying Jacobians and mass matrices, see the ODEFunction documentation.

Fields

• f: The function in the ODE.
• u0: The initial condition.
• tspan: The timespan for the problem.
• p: The parameters.
• kwargs: The keyword arguments passed onto the solves.

Example Problem

using SciMLBase
function lorenz!(du,u,p,t)
  du[1] = 10.0(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)

# Test that it worked
using OrdinaryDiffEq
sol = solve(prob,Tsit5())
using Plots; plot(sol,vars=(1,2,3))

More Example Problems

Example problems can be found in DiffEqProblemLibrary.jl.

To use a sample problem, such as prob_ode_linear, you can do something like:

#] add ODEProblemLibrary
using ODEProblemLibrary
prob = ODEProblemLibrary.prob_ode_linear
sol = solve(prob)

ODEProblem(f::ODEFunction, u0, tspan, p = NullParameters(), callback = CallbackSet())

Define an ODE problem from an ODEFunction.

struct ODESolution{T, N, uType, uType2, DType, tType, rateType, discType, P, A, IType, S, AC<:Union{Nothing, Vector{Int64}}, R, O, V, GE} <: SciMLBase.AbstractODESolution{T, N, uType}

Representation of the solution to an ordinary differential equation defined by an ODEProblem.

DESolution Interface

For more information on interacting with DESolution types, check out the Solution Handling page of the DifferentialEquations.jl documentation.

https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/

Fields

• u: the representation of the ODE solution. Given as an array of solutions, where u[i] corresponds to the solution at time t[i]. It is recommended in most cases one does not access sol.u directly and instead use the array interface described in the Solution Handling page of the DifferentialEquations.jl documentation.
• t: the time points corresponding to the saved values of the ODE solution.
• prob: the original ODEProblem that was solved.
• alg: the algorithm type used by the solver.
• stats: statistics of the solver, such as the number of function evaluations required, number of Jacobians computed, and more.
• retcode: the return code from the solver. Used to determine whether the solver solved successfully, whether it terminated early due to a user-defined callback, or whether it exited due to an error. For more details, see the return code documentation.
• global_error: an estimate of the global (accumulated) error of the solution, or nothing when the algorithm does not compute one (see has_global_error). When present it is an array matching u, with global_error[i] the estimated global error of u[i] at t[i].

Defines a second order ordinary differential equation (ODE) problem. Documentation Page: https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/

Mathematical Specification of a 2nd Order ODE Problem

To define a 2nd Order ODE Problem, you simply need to give the function $f$ and the initial condition $u_0$ which define an ODE:

\[u'' = f(u',u,p,t)\]

f should be specified as f(du,u,p,t) (or in-place as f(ddu,du,u,p,t)), and u₀ should be an AbstractArray (or number) whose geometry matches the desired geometry of u. Note that we are not limited to numbers or vectors for u₀; one is allowed to provide u₀ as arbitrary matrices / higher dimension tensors as well.

From this form, a dynamical ODE:

\[\begin{align*} v' &= f(v,u,p,t) \\ u' &= v \end{align*}\]

is generated.

Constructors

SecondOrderODEProblem{isinplace}(f,du0,u0,tspan,callback=CallbackSet())

Defines the ODE with the specified functions.

Fields

• f: The function for the second derivative.
• du0: The initial derivative.
• u0: The initial condition.
• tspan: The timespan for the problem.
• callback: A callback to be applied to every solver which uses the problem. Defaults to nothing.

struct SplitFunction{iip, specialize, F1, F2, TMM, C, Ta, Tt, TJ, JVP, VJP, JP, WP, SP, TW, TWt, TPJ, O, TCV, SYS, ID<:Union{Nothing, SciMLBase.OverrideInitData}, NLP<:Union{Nothing, SciMLBase.ODENLStepData}} <: SciMLBase.AbstractODEFunction{iip}

A representation of a split ODE function f, defined by:

\[M \frac{du}{dt} = f_1(u,p,t) + f_2(u,p,t)\]

and all of its related functions, such as the Jacobian of f, its gradient with respect to time, and more. For all cases, u0 is the initial condition, p are the parameters, and t is the independent variable.

Generally, for ODE integrators the f_1 portion should be considered the "stiff portion of the model" with larger timescale separation, while the f_2 portion should be considered the "non-stiff portion". This interpretation is directly used in integrators like IMEX (implicit-explicit integrators) and exponential integrators.

Constructor

SplitFunction{iip,specialize}(f1,f2;
                             mass_matrix = __has_mass_matrix(f1) ? f1.mass_matrix : I,
                             analytic = __has_analytic(f1) ? f1.analytic : nothing,
                             tgrad= __has_tgrad(f1) ? f1.tgrad : nothing,
                             jac = __has_jac(f1) ? f1.jac : nothing,
                             jvp = __has_jvp(f1) ? f1.jvp : nothing,
                             vjp = __has_vjp(f1) ? f1.vjp : nothing,
                             jac_prototype = __has_jac_prototype(f1) ? f1.jac_prototype : nothing,
                             W_prototype = __has_W_prototype(f1) ? f1.W_prototype : nothing,
                             sparsity = __has_sparsity(f1) ? f1.sparsity : jac_prototype,
                             paramjac = __has_paramjac(f1) ? f1.paramjac : nothing,
                             colorvec = __has_colorvec(f1) ? f1.colorvec : nothing,
                             sys = __has_sys(f1) ? f1.sys : nothing)

Note that only the functions f_i themselves are required. These functions should be given as f_i!(du,u,p,t) or du = f_i(u,p,t). See the section on iip for more details on in-place vs out-of-place handling.

All of the remaining functions are optional for improving or accelerating the usage of the SplitFunction. These include:

• mass_matrix: the mass matrix M represented in the ODE function. Can be used to determine that the equation is actually a differential-algebraic equation (DAE) if M is singular. Note that in this case special solvers are required, see the DAE solver page for more details: https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/. Must be an AbstractArray or an AbstractSciMLOperator.
• analytic(u0,p,t): used to pass an analytical solution function for the analytical solution of the ODE. Generally only used for testing and development of the solvers.
• tgrad(dT,u,p,t) or dT=tgrad(u,p,t): returns $\frac{∂f_1(u,p,t)}{∂t}$
• jac(J,u,p,t) or J=jac(u,p,t): returns $\frac{df_1}{du}$
• jvp(Jv,v,u,p,t) or Jv=jvp(v,u,p,t): returns the directional derivative $\frac{df_1}{du} v$
• vjp(Jv,v,u,p,t) or Jv=vjp(v,u,p,t): returns the adjoint derivative $\frac{df_1}{du}^∗ v$
• jac_prototype: a prototype matrix matching the type that matches the Jacobian. For example, if the Jacobian is tridiagonal, then an appropriately sized Tridiagonal matrix can be used as the prototype and integrators will specialize on this structure where possible. Non-structured sparsity patterns should use a SparseMatrixCSC with a correct sparsity pattern for the Jacobian. The default is nothing, which means a dense Jacobian.
• W_prototype: a prototype matrix matching the type that matches the W matrix. For example, if the Jacobian is tridiagonal, and the mass_matrix is diagonal, then an appropriately sized Tridiagonal matrix can be used as the prototype and integrators will specialize on this structure where possible. Non-structured sparsity patterns should use a SparseMatrixCSC with a correct sparsity pattern for the W matrix. The default is nothing, which means a W of appropriate type for the jacobian and linear solver
• paramjac(pJ,u,p,t): returns the parameter Jacobian $\frac{df_1}{dp}$.
• colorvec: a color vector according to the SparseDiffTools.jl definition for the sparsity pattern of the jac_prototype. This specializes the Jacobian construction when using finite differences and automatic differentiation to be computed in an accelerated manner based on the sparsity pattern. Defaults to nothing, which means a color vector will be internally computed on demand when required. The cost of this operation is highly dependent on the sparsity pattern.

Note on the Derivative Definition

The derivatives, such as the Jacobian, are only defined on the f1 portion of the split ODE. This is used to treat the f1 implicit while keeping the f2 portion explicit.

iip: In-Place vs Out-Of-Place

For more details on this argument, see the ODEFunction documentation.

specialize: Controlling Compilation and Specialization

For more details on this argument, see the ODEFunction documentation.

Fields

The fields of the SplitFunction type directly match the names of the inputs.

Symbolically Generating the Functions

See the modelingtoolkitize function from ModelingToolkit.jl for automatically symbolically generating the Jacobian and more from the numerically-defined functions. See ModelingToolkit.SplitODEProblem for information on generating the SplitFunction from this symbolic engine.

Defines a split ordinary differential equation (ODE) problem. Documentation Page: https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/

Mathematical Specification of a Split ODE Problem

To define a SplitODEProblem, you simply need to give two functions $f_1$ and $f_2$ along with an initial condition $u_0$ which define an ODE:

\[\frac{du}{dt} = f_1(u,p,t) + f_2(u,p,t)\]

f should be specified as f(u,p,t) (or in-place as f(du,u,p,t)), and u₀ should be an AbstractArray (or number) whose geometry matches the desired geometry of u. Note that we are not limited to numbers or vectors for u₀; one is allowed to provide u₀ as arbitrary matrices / higher dimension tensors as well.

Many splits are at least partially linear. That is the equation:

\[\frac{du}{dt} = Au + f_2(u,p,t)\]

For how to define a linear function A, see the documentation for the AbstractSciMLOperator.

Constructors

SplitODEProblem(f::SplitFunction,u0,tspan,p=NullParameters();kwargs...)
SplitODEProblem{isinplace}(f1,f2,u0,tspan,p=NullParameters();kwargs...)

The isinplace parameter can be omitted and will be determined using the signature of f2. Note that both f1 and f2 should support the in-place style if isinplace is true or they should both support the out-of-place style if isinplace is false. You cannot mix up the two styles.

Parameters are optional, and if not given, then a NullParameters() singleton will be used which will throw nice errors if you try to index non-existent parameters. Any extra keyword arguments are passed on to the solvers. For example, if you set a callback in the problem, then that callback will be added in every solve call.

Under the hood, a SplitODEProblem is just a regular ODEProblem whose f is a SplitFunction. Therefore, you can solve a SplitODEProblem using the same solvers for ODEProblem. Solver packages document which methods specialize on split structure.

For specifying Jacobians and mass matrices, see the SciMLFunctions interface page.

Fields

• f1, f2: The functions in the ODE.
• u0: The initial condition.
• tspan: The timespan for the problem.
• p: The parameters for the problem. Defaults to NullParameters
• kwargs: The keyword arguments passed onto the solves.

SplitODEProblem(f, u0, tspan; ...)
SplitODEProblem(f, u0, tspan, p; kwargs...)

Define a split ODE problem from a SplitFunction.

struct CallbackSet{T1<:Union{Tuple, AbstractVector}, T2<:Union{Tuple, AbstractVector}} <: SciMLBase.DECallback

Container for the callbacks attached to a differential equation solve.

Multiple callbacks can be chained together to form a CallbackSet. A CallbackSet is constructed by passing ContinuousCallback, DiscreteCallback, VectorContinuousCallback, nothing, or other CallbackSet instances:

CallbackSet(cb1,cb2,cb3)

You can pass as many callbacks as needed. Nested callback sets are flattened into two ordered collections, continuous_callbacks and discrete_callbacks. Public constructors use tuples; solver paths may use vectors when callback types must be erased to reuse compilation.

When a solver encounters multiple callbacks, the following rules apply:

• ContinuousCallbacks and VectorContinuousCallbacks are applied before DiscreteCallbacks. (This is because they often implement event-finding that will backtrack the timestep to smaller than dt).
• For ContinuousCallbacks and VectorContinuousCallbacks, the event times are found by rootfinding and only the first ContinuousCallback or VectorContinuousCallback affect is applied.
• The DiscreteCallbacks are then applied in order. Note that the ordering only matters for the conditions: if a previous callback modifies u in such a way that the next callback no longer evaluates condition to true, its affect will not be applied.

ContinuousCallback(
    condition, affect!, affect_neg!;
    initialize = INITIALIZE_DEFAULT,
    finalize = FINALIZE_DEFAULT,
    idxs = nothing,
    rootfind = LeftRootFind,
    save_positions = (true, true),
    interp_points = 10,
    abstol = 10eps(), reltol = 0, repeat_nudge = 1 // 100,
    initializealg = nothing, maybe_discontinuity = true
)

ContinuousCallback(
    condition, affect!;
    initialize = INITIALIZE_DEFAULT,
    finalize = FINALIZE_DEFAULT,
    idxs = nothing,
    rootfind = LeftRootFind,
    save_positions = (true, true),
    affect_neg! = affect!,
    interp_points = 10,
    abstol = 10eps(), reltol = 0, repeat_nudge = 1 // 100,
    initializealg = nothing, maybe_discontinuity = true
)

Contains a single callback whose condition is a continuous function. The callback is triggered when this function evaluates to 0.

Arguments

• condition: This is a function condition(u,t,integrator) for declaring when the callback should be used. A callback is initiated if the condition hits 0 within the time interval. See the Integrator Interface documentation for information about integrator.
• affect!: This is the function affect!(integrator) where one is allowed to modify the current state of the integrator. If you do not pass an affect_neg! function, it is called when condition is found to be 0 (at a root) and the cross is either an upcrossing (from negative to positive) or a downcrossing (from positive to negative). You need to explicitly pass nothing as the affect_neg! argument if it should only be called at upcrossings, e.g. ContinuousCallback(condition, affect!, nothing). For more information on what can be done, see the Integrator Interface manual page. Modifications to u are safe in this function.
• affect_neg!=affect!: This is the function affect_neg!(integrator) where one is allowed to modify the current state of the integrator. This is called when condition is found to be 0 (at a root) and the cross is a downcrossing (from positive to negative). For more information on what can be done, see the Integrator Interface manual page. Modifications to u are safe in this function.
• rootfind=LeftRootFind: This is a flag to specify the type of rootfinding to do for finding event location. If this is set to LeftRootFind, the solution will be backtracked to the point where condition==0 and if the solution isn't exact, the left limit of root is used. If set to RightRootFind, the solution would be set to the right limit of the root. Otherwise, the systems and the affect! will occur at t+dt. Note that these enums are not exported, and thus one needs to reference them as SciMLBase.LeftRootFind, SciMLBase.RightRootFind, or SciMLBase.NoRootFind.
• interp_points=10: The number of interpolated points to check the condition. The condition is found by checking whether any interpolation point / endpoint has a different sign. If interp_points=0, then conditions will only be noticed if the sign of condition is different at t than at t+dt. This behavior is not robust when the solution is oscillatory, and thus it's recommended that one use some interpolation points (they're cheap to compute!). 0 within the time interval.
• save_positions=(true,true): Boolean tuple for whether to save before and after the affect!. This saving will occur just before and after the event, only at event times, and does not depend on options like saveat, save_everystep, etc. (i.e. if saveat=[1.0,2.0,3.0], this can still add a save point at 2.1 if true). For discontinuous changes like a modification to u to be handled correctly (without error), one should set save_positions=(true,true).
• idxs=nothing: The components which will be interpolated into the condition. Defaults to nothing which means u will be all components.
• initialize: This is a function (c,u,t,integrator) which can be used to initialize the state of the callback c. It should modify the argument c and the return is ignored.
• finalize: This is a function (c,u,t,integrator) which can be used to finalize the state of the callback c. It can modify the argument c and the return is ignored.
• abstol=10eps(): Tolerance for repeated event prevention. If the callback was just triggered and the new starting condition is less than the tolerance from its value at the root, then the next testing point will be nudged to avoid repeats. If the callback does not mutate the integrator in a way that affect the condition, this can be safely set to 0.0. reltol is deprecated.
• repeat_nudge = 1//100: This is used to set the next testing point after a previously found zero. Defaults to 1//100, which means after a callback, the next sign check will take place at t + dt*1//100 instead of at t to avoid repeats.
• initializealg = nothing: In the context of a DAE, this is the algorithm that is used to run initialization after the effect. The default of nothing defers to the initialization algorithm provided in the solve.
• maybe_discontinuity = true: Declares whether the condition time could have a discontinuity or the affect! could introduce a discontinuity. Defaults to true. This is only used if discontinuity detection is enabled in the controller (i.e. discontinuity_handling = true).

Extended help

• saved_clock_partitions: An iterable of clock partition indices to save after the callback triggers. MTK-only API.
• initialize_save_discretes = true: Whether callback initialization should save the discrete parameter partitions listed in saved_clock_partitions when save_positions[2] is true.

DiscreteCallback(
    condition, affect!;
    initialize = INITIALIZE_DEFAULT,
    finalize = FINALIZE_DEFAULT,
    save_positions = (true, true),
    initializealg = nothing
)

Arguments

• condition: This is a function condition(u,t,integrator) for declaring when the callback should be used. A callback is initiated if the condition evaluates to true. See the Integrator Interface documentation for information about integrator.
• affect!: This is the function affect!(integrator) where one is allowed to modify the current state of the integrator. For more information on what can be done, see the Integrator Interface manual page.
• save_positions: Boolean tuple for whether to save before and after the affect!. This saving will occur just before and after the event, only at event times, and does not depend on options like saveat, save_everystep, etc. (i.e. if saveat=[1.0,2.0,3.0], this can still add a save point at 2.1 if true). For discontinuous changes like a modification to u to be handled correctly (without error), one should set save_positions=(true,true).
• initialize: This is a function (c,u,t,integrator) which can be used to initialize the state of the callback c. It should modify the argument c and the return is ignored.
• finalize: This is a function (c,u,t,integrator) which can be used to finalize the state of the callback c. It can modify the argument c and the return is ignored.
• initializealg = nothing: In the context of a DAE, this is the algorithm that is used to run initialization after the effect. The default of nothing defers to the initialization algorithm provided in the solve.

Extended help

• saved_clock_partitions: An iterable of clock partition indices to save after the callback triggers. MTK-only API.
• initialize_save_discretes = true: Whether callback initialization should save the discrete parameter partitions listed in saved_clock_partitions when save_positions[2] is true.

VectorContinuousCallback(
    condition, affect!, len;
    initialize = INITIALIZE_DEFAULT,
    finalize = FINALIZE_DEFAULT,
    idxs = nothing,
    rootfind = LeftRootFind,
    save_positions = (true, true),
    interp_points = 10,
    abstol = 10eps(), reltol = 0, repeat_nudge = 1 // 100,
    initializealg = nothing, maybe_discontinuity = true
)

This is also a subtype of AbstractContinuousCallback. CallbackSet is not feasible when you have many callbacks, as it doesn't scale well. For this reason, we have VectorContinuousCallback - it allows you to have a single callback for multiple events.

VectorContinuousCallback intentionally does not have an affect_neg! callback. Its affect! receives the triggering event index and is responsible for applying the appropriate effect for that event.

Arguments

• condition: This is a function condition(out, u, t, integrator) which should save the condition value in the array out at the right index. Maximum index of out should be specified in the len property of callback. So, this way you can have a chain of len events, which would cause the ith event to trigger when out[i] = 0.
• affect!: This is a function affect!(integrator, event_index) which lets you modify integrator and it tells you about which event occurred using event_idx i.e. gives you index i for which out[i] came out to be zero.
• len: Number of callbacks chained. This is compulsory to be specified.

Rest of the arguments have the same meaning as in ContinuousCallback.

Extended help

• saved_clock_partitions: An iterable of len elements, where the ith element is an iterable of clock partition indices to save when the ith event triggers. MTK-only API.
• initialize_save_discretes = true: Whether callback initialization should save the discrete parameter partitions listed in saved_clock_partitions when save_positions[2] is true.

ODEAliasSpecifier(;alias_p = nothing, alias_f = nothing, alias_u0 = nothing, alias_du0 = nothing, alias_tstops = nothing, alias = nothing)

Control which ODE problem inputs and solver option arrays may be aliased.

alias_u0 controls the initial state, alias_du0 controls an initial derivative array when the problem representation has one, alias_p controls the parameter object, alias_f controls the function object, and alias_tstops controls the tstops vector passed to the solver. A value of nothing delegates to the solver default. Set alias = true or alias = false to apply the same policy to all fields.

Keywords

• alias_p::Union{Bool, Nothing}: alias the parameter object.
• alias_f::Union{Bool, Nothing}: alias the ODE function object.
• alias_u0::Union{Bool, Nothing}: alias the u0 array.
• alias_du0::Union{Bool, Nothing}: alias the du0 array, when present.
• alias_tstops::Union{Bool, Nothing}: alias the tstops array.
• alias::Union{Bool, Nothing}: set every field of the ODEAliasSpecifier.

add_saveat!(i::DEIntegrator, t)

Schedule solution output at the future physical time t.

An integrator must not accept a save point behind its current time in the direction of integration. saveat normally uses interpolation when t lies inside a step and therefore does not force the integrator to step exactly to t. Add a matching add_tstop! when an exact step endpoint is also required. Saving still follows the solver's save_on, save_idxs, and related output options.

add_tstop!(i::DEIntegrator, t)

Schedule a future stopping time at the physical time t.

An integrator must not accept a stop behind its current time in the direction of integration. A tstop constrains stepping so the integrator reaches t exactly when the method supports step-size changes or interpolation. It does not by itself request that the solution be saved there; use add_saveat! or the solver's saving options for output.

Implementations commonly store tstops as direction-normalized priority keys integrator.tdir * t. The companion queue accessors expose those keys so generic stepping code can compare them with integrator.tdir * integrator.t in both forward and reverse integration.

auto_dt_reset!(integrator::DEIntegrator)

Recompute the integrator's initial step size from its current state.

Concrete solvers should apply the same automatic step-size selection used during init, including the current state, time, parameters, tolerances, integration direction, and method-specific limits. They must update the active step size and any proposal state needed by the next step. This operation may evaluate the problem function and increment solver statistics. Its return value is not part of the interface.

derivative_discontinuity!(i::DEIntegrator, bool)

Record whether a callback or direct integrator mutation introduced a derivative discontinuity.

The flag describes whether f(u, p, t) may have changed discontinuously because u, p, t, or the definition of f changed. Solvers use this to decide whether to recompute derivatives, interpolation data, FSAL caches, or Jacobians before the next step. Callback code should leave the default discontinuity behavior in place after state-changing effects, and may call derivative_discontinuity!(integrator, false) only when it did not change the state, parameters, time, or dynamics.

has_global_error(alg::AbstractDEAlgorithm)

Trait declaring whether an algorithm computes an estimate of the global (accumulated) error of the solution.

Return true when the solver or solver wrapper populates the global_error field of the solution with a per-time-point estimate of the global error (the difference between the numerical and true solutions), rather than only controlling the local error of each step. When true, the solution's global_error is an array matching u; when false (the default) it is nothing.

The default is false.

reinit!(integrator::DEIntegrator, args...; kwargs...)

The reinit function lets you restart the integration at a new value.

Arguments

• u0: Value of u to start at. Default value is integrator.sol.prob.u0

Keyword Arguments

• t0: Starting timepoint. Default value is integrator.sol.prob.tspan[1]
• tf: Ending timepoint. Default value is integrator.sol.prob.tspan[2]
• erase_sol=true: Whether to start with no other values in the solution, or keep the previous solution.
• tstops, d_discontinuities, & saveat: Cache where these are stored. Default is the original cache.
• reset_dt: Set whether to reset the current value of dt using the automatic dt determination algorithm. Default is (integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive
• reinit_callbacks: Set whether to run the callback initializations again (and initialize_save is for that). Default is true.
• reinit_cache: Set whether to re-run the cache initialization function (i.e. resetting FSAL, not allocating vectors) which should usually be true for correctness. Default is true.

Additionally, once can access auto_dt_reset! which will run the auto dt initialization algorithm.

remake(thing; <keyword arguments>)

Re-construct thing with new field values specified by the keyword arguments.

remake(
    prob::DAEProblem; f = missing, du0 = missing, u0 = missing, tspan = missing,
    p = missing, differential_vars = missing, kwargs = missing, _kwargs...
)

Remake the given DAEProblem. If u0 or p are given as symbolic maps ModelingToolkit.jl has to be loaded.

remake(
    prob::NonlinearLeastSquaresProblem; f = missing, u0 = missing, p = missing,
    kwargs = missing, _kwargs...
)

Remake the given NonlinearLeastSquaresProblem.

remake(
    prob::NonlinearProblem; f = missing, u0 = missing, p = missing,
    problem_type = missing, kwargs = missing, _kwargs...
)

Remake the given NonlinearProblem. If u0 or p are given as symbolic maps ModelingToolkit.jl has to be loaded.

remake(
    prob::ODEProblem; f = missing, u0 = missing, tspan = missing,
    p = missing, kwargs = missing, _kwargs...
)

Remake the given ODEProblem. If u0 or p are given as symbolic maps ModelingToolkit.jl has to be loaded.

remake(
    prob::OptimizationProblem; f = missing, u0 = missing, p = missing,
    lb = missing, ub = missing, int = missing, lcons = missing, ucons = missing,
    sense = missing, kwargs = missing, _kwargs...
)

Remake the given OptimizationProblem. If u0 or p are given as symbolic maps ModelingToolkit.jl has to be loaded.

remake(
    prob::SCCNonlinearProblem; u0 = missing, p = missing, probs = missing,
    parameters_alias = prob.parameters_alias, sys = missing, explicitfuns! = missing
)

Remake the given SCCNonlinearProblem. u0 is the state vector for the entire problem, which will be chunked appropriately and used to remake the individual subproblems. p is the parameter object for prob. If parameters_alias, the same parameter object will be used to remake the individual subproblems. Otherwise if p !== missing, this function will error and require that probs be specified. probs is the collection of subproblems. Even if probs is explicitly specified, the value of u0 provided to remake will be used to override the values in probs. sys is the index provider for the full system.

remake(
    prob::SDEProblem; f = missing, g = missing, u0 = missing, tspan = missing,
    p = missing, noise = missing, noise_rate_prototype = missing,
    seed = missing, kwargs = missing, _kwargs...
)

Remake the given SDEProblem.

remake(func::AbstractSciMLFunction; f = missing, g = missing, f2 = missing, kwargs...)

remake the given func. Return an AbstractSciMLFunction of the same kind, isinplace and specialization as func. Retain the properties of func, except those that are overridden by keyword arguments. For stochastic functions (e.g. SDEFunction) the g keyword argument is used to override func.g. For split functions (e.g. SplitFunction) the f2 keyword argument is used to override func.f2, and f is used for func.f1. If f isa AbstractSciMLFunction and func is not a split function, properties of f will override those of func (but not ones provided via keyword arguments). Properties of f that are nothing will fall back to those in func (unless provided via keyword arguments). If f is a different type of AbstractSciMLFunction from func, the returned function will be of the kind of f unless func is a split function. If func is a split function, f and f2 will be wrapped in the appropriate AbstractSciMLFunction type with the same isinplace and specialization as func.

remake(
    prob::AbstractSciMLProblem; u0 = missing, p = missing,
    interpret_symbolicmap = true, use_defaults = false, kwargs...
)

Construct a problem of the same family as prob, replacing the supplied fields and preserving all other problem data. Extra keyword arguments are forwarded to the problem-family-specific remake implementation.

When u0 or p is a symbolic map and prob has an associated symbolic system, explicit entries take precedence. Missing entries with symbolic-expression defaults use those expressions so dependent values remain consistent. Other missing entries retain their values from prob unless use_defaults = true, in which case available numeric system defaults are preferred. use_defaults is meaningful only when an explicit symbolic map is supplied for a problem with an associated system.

Set interpret_symbolicmap = false to use a pair-valued p directly instead of interpreting it as a symbolic parameter map. Pair-valued u0 is still interpreted as a symbolic state map. Non-symbolic u0 and p values are used directly.

remake(
    prob::BVProblem; f = missing, u0 = missing, tspan = missing,
    p = missing, kwargs = missing, problem_type = missing, _kwargs...
)

Remake the given BVProblem.

set_proposed_dt!(i::DEIntegrator, dt)
set_proposed_dt!(i::DEIntegrator, i2::DEIntegrator)

Set the signed step-size proposal used for the next step.

The scalar form updates every step-size field that the concrete solver requires to honor a new proposal. It does not bypass error control, rejection, or tstop handling, and therefore does not guarantee that the next accepted step has exactly that size.

The two-integrator form synchronizes the first integrator's time-stepping state with the second. Adaptive implementations should copy the controller history or other state needed to reproduce the proposal, rather than only copying one dt field. This form is optional for integrators that cannot share compatible controller state.

successful_retcode(retcode::ReturnCode.T)::Bool
successful_retcode(sol::AbstractSciMLSolution)::Bool

Return whether a return code represents a successful solver outcome.

The solution form forwards to successful_retcode(sol.retcode). Use this predicate instead of comparing only with ReturnCode.Success, because requested termination and other informative terminal states can also be successful.

The successful codes are ReturnCode.Success, ReturnCode.Terminated, ReturnCode.ExactSolutionLeft, ReturnCode.ExactSolutionRight, ReturnCode.FloatingPointLimit, and ReturnCode.StalledSuccess. ReturnCode.Default is not successful because it means the solve is unfinished or its outcome is unknown. ReturnCode.Stalled and all failure codes are also unsuccessful.

AbstractADType

Abstract supertype for all AD choices.

Extension contract

External packages may subtype AbstractADType to describe an AD backend. They must also implement mode for their concrete subtype and return an instance of an AbstractMode subtype. Consumers should dispatch on mode(ad), rather than on a downstream concrete AD type.

AbstractColoringAlgorithm

Abstract supertype for Jacobian/Hessian coloring algorithms.

Extension contract

External algorithms implement the supported public coloring functions. Each result must be an AbstractVector of integers: column colorings have length size(M, 2), row colorings have length size(M, 1), and symmetric colorings require a square matrix and have length size(M, 1). The assigned colors must satisfy the structural orthogonality condition documented by each coloring function. Unsupported coloring forms must throw an error.

Note

The terminology and definitions are taken from the following paper:

What Color Is Your Jacobian? Graph Coloring for Computing Derivatives, Assefaw Hadish Gebremedhin, Fredrik Manne, and Alex Pothen (2005)

AbstractMode

Abstract supertype for the traits identifying differentiation modes.

Subtypes

• ForwardMode
• ReverseMode
• ForwardOrReverseMode
• SymbolicMode

AbstractSparsityDetector

Abstract supertype for sparsity pattern detectors.

Extension contract

External detectors implement the supported public forms of jacobian_sparsity and/or hessian_sparsity. A Jacobian pattern must be an AbstractMatrix{Bool} with shape (length(y), length(x)), where y is f(x) for the out-of-place form or the supplied output buffer for the in-place form. A Hessian pattern must be an AbstractMatrix{Bool} with shape (length(x), length(x)). Methods for unsupported operations must throw an error; they must not return an unrelated pattern.

New detectors should return Bool patterns, but consumers must not require Bool: they should treat every nonzero entry as a structural nonzero. Patterns given as integer or floating point matrices of ones and zeroes are common in the wild, are accepted unchanged by KnownJacobianSparsityDetector and KnownHessianSparsityDetector, and are handled correctly by the coloring and decompression implementations used downstream (SparseMatrixColorings.jl through DifferentiationInterface.jl). The element type of the pattern does not propagate to the differentiation result.

For a sparse matrix, the stored structure defines the pattern: an explicitly stored zero is treated as a structural nonzero, which yields a valid but more conservative coloring than the same matrix passed through dropzeros.

AutoFiniteDiff{T1,T2,T3}

Struct used to select the FiniteDiff.jl backend for automatic differentiation.

Defined by ADTypes.jl.

Constructors

AutoFiniteDiff(;
    fdtype=Val(:forward), fdjtype=fdtype, fdhtype=Val(:hcentral),
    relstep=nothing, absstep=nothing, dir=true
)

Fields

• fdtype::T1: finite difference type
• fdjtype::T2: finite difference type for the Jacobian
• fdhtype::T3: finite difference type for the Hessian
• relstep: relative finite difference step size
• absstep: absolute finite difference step size
• dir: direction of the finite difference step

AutoForwardDiff{chunksize,T}

Struct used to select the ForwardDiff.jl backend for automatic differentiation.

Defined by ADTypes.jl.

Constructors

AutoForwardDiff(; chunksize=nothing, tag=nothing)

Type parameters

• chunksize: the preferred chunk size to evaluate several derivatives at once

Fields

• tag::T: a custom tag to handle nested differentiation calls (usually not necessary)

AutoSparse{D,S,C}

Wraps an ADTypes.jl object to deal with sparse Jacobians and Hessians.

Fields

• dense_ad::D: the underlying AD package, subtyping AbstractADType
• sparsity_detector::S: the sparsity pattern detector, subtyping AbstractSparsityDetector
• coloring_algorithm::C: the coloring algorithm, subtyping AbstractColoringAlgorithm

Constructors

AutoSparse(
    dense_ad;
    sparsity_detector=ADTypes.NoSparsityDetector(),
    coloring_algorithm=ADTypes.NoColoringAlgorithm()
)

ForwardMode

Trait for AD choices that rely on forward mode algorithmic differentiation or finite differences.

These two paradigms are classified together because they can both efficiently compute Jacobian-vector products.

ForwardOrReverseMode

Trait for AD choices that can work either in ForwardMode or ReverseMode, depending on their configuration.

KnownHessianSparsityDetector(hessian_sparsity::AbstractMatrix) <: AbstractSparsityDetector

Trivial sparsity detector used to return a known Hessian sparsity pattern.

AbstractMatrix{Bool} is the canonical pattern type, but the element type is neither converted nor checked: hessian_sparsity hands the pattern back exactly as given, and consumers treat every nonzero entry as a structural nonzero, so integer or floating point matrices of ones and zeroes work as well (see the extension contract of AbstractSparsityDetector).

See also

• AbstractSparsityDetector
• KnownJacobianSparsityDetector

KnownJacobianSparsityDetector(jacobian_sparsity::AbstractMatrix) <: AbstractSparsityDetector

Trivial sparsity detector used to return a known Jacobian sparsity pattern.

AbstractMatrix{Bool} is the canonical pattern type, but the element type is neither converted nor checked: jacobian_sparsity hands the pattern back exactly as given, and consumers treat every nonzero entry as a structural nonzero, so integer or floating point matrices of ones and zeroes work as well (see the extension contract of AbstractSparsityDetector).

See also

• AbstractSparsityDetector
• KnownHessianSparsityDetector

ReverseMode

Trait for AD choices that rely on reverse mode algorithmic differentiation.

SymbolicMode

Trait for AD choices that rely on symbolic differentiation.

hessian_sparsity(f, x, sd::AbstractSparsityDetector)::AbstractMatrix{Bool}

Use detector sd to construct a (typically sparse) matrix S describing the pattern of nonzeroes in the Hessian of f applied at x.

jacobian_sparsity(f, x, sd::AbstractSparsityDetector)::AbstractMatrix{Bool}
jacobian_sparsity(f!, y, x, sd::AbstractSparsityDetector)::AbstractMatrix{Bool}

Use detector sd to construct a (typically sparse) matrix S describing the pattern of nonzeroes in the Jacobian of f (resp. f!) applied at x (resp. (y, x)).

mode(ad::AbstractADType)

Return the differentiation mode of ad, as a subtype of AbstractMode.

Extension contract

Every external concrete subtype of AbstractADType must provide this method. Return an instance of the most specific applicable mode trait; callers may use that trait for dispatch but must not require a particular backend implementation.

abstract type AbstractSciMLOperator{T}

AbstractSciMLOperator is the extension point for matrix-like and matrix-free operators. A subtype represents an operator $L(u,p,t)$ whose action on an array $v$ is written $L(u,p,t)v$. The subtype may be constant, state-dependent, or time-dependent, and may be composed with other SciML operators through the lazy algebra.

This is an interface type, not a constructor. The concrete type should be public only when users are expected to construct or extend it; otherwise use the qualified developer-facing API documented in this section.

Mathematical Notation

An AbstractSciMLOperator$L$ is an operator which is used to represent the following type of equation:

\[w = L(u,p,t)[v]\]

where L[v] is the operator application of $L$ on the vector $v$.

Construction and Extension Rules

AbstractSciMLOperator is an interface, not a concrete constructor. New operator types should subtype it with the scalar element type T and must preserve their mathematical action when they participate in lazy algebra.

Required Interface

A concrete subtype must implement size(L) -> (m, n), *(L, v), and mul!(w, L, v). The returned action has leading size m for an input whose leading size is n, and the in-place method must return w after writing the result. The scaling form mul!(w, L, v, α, β) is required when has_mul!(L) == true; it must compute $w \leftarrow α(Lv) + βw$.

has_mul(L), has_mul!(L), has_ldiv(L), and has_ldiv!(L) are promises, not capability probes: return true only when the corresponding operation is valid for all compatible inputs. convert(AbstractMatrix, L) is optional and should be defined only when isconvertible(L) == true; its result must have the same size and action as L in its current state.

An AbstractSciMLOperator can be called like a function in the following ways:

• L(v, u, p, t) - Out-of-place application where v is the action vector and u is the update vector
• L(w, v, u, p, t) - In-place application where w is the destination, v is the action vector, and u is the update vector
• L(w, v, u, p, t, α, β) - In-place application with scaling: w = α*(L*v) + β*w

Operator state can be updated separately from application:

• update_coefficients!(L, u, p, t) for in-place operator update
• L = update_coefficients(L, u, p, t) for out-of-place operator update

SciMLOperators also overloads Base.*, LinearAlgebra.mul!, LinearAlgebra.ldiv! for operator evaluation without updating operator state. An AbstractSciMLOperator behaves like a matrix in these methods. Allocation-free methods, suffixed with a ! often need cache arrays. To precache an AbstractSciMLOperator, call the function L = cache_operator(L, input_vector).

Required Interface For Subtypes

A concrete subtype must define Base.size(L) and one of the following application paths:

• Base.:*(L, v) for out-of-place matrix-like application.
• LinearAlgebra.mul!(w, L, v) and, when has_mul!(L) is true, LinearAlgebra.mul!(w, L, v, α, β) for in-place application.
• Base.convert(AbstractMatrix, L) when isconvertible(L) is true.

If the operator state depends on (u, p, t) or accepted keyword arguments, the subtype must implement update_coefficients for out-of-place state updates or update_coefficients! for in-place state updates. The out-of-place form returns a new operator and leaves L unchanged; the in-place form returns nothing. Composite operators assume these update methods may be called recursively on every operator returned by getops(L).

The positional arguments are forwarded unchanged through a composite operator. u is the state supplied by the caller and is not necessarily the same shape as the action vector v; an operator whose action is nonlinear in v should generally use FunctionOperator and report islinear(L) == false. For a constant leaf, the default update is a no-op. A stateful leaf must override isconstant rather than inheriting the empty-child default.

Subtypes that need preallocated work arrays for allocation-free application must implement cache_self(L, v) for their own caches, cache_internals(L, v) for child-operator caches, or both. cache_operator(L, v) calls these hooks and downstream solvers may call it before repeated mul! evaluations.

Caching Rules

cache_operator(L, v) may return either L or a cached replacement. A subtype that advertises has_mul!(L) == true must ensure the cached result is ready for repeated mul! calls with compatible vectors. Cache hooks may not change the mathematical action, dimensions, or trait values of the operator. A cached operator's scratch is mutable and is not safe to use concurrently unless the subtype explicitly provides that guarantee.

Composite types expose their children through the developer-facing getops method. A new composite must forward state updates, caching, and traits to every child that contributes to its action. The public action must remain unchanged by flattening, caching, or updating the composition.

Trait Rules

Trait functions such as isconstant, islinear, isconvertible, has_concretization, has_mul, has_mul!, has_ldiv, and has_ldiv! are part of the public operator interface. A trait returning true is a promise that the corresponding operation is valid for inputs with compatible sizes. For example, has_mul!(L) means mul!(w, L, v) is available, and has_concretization(L) means either convert(AbstractMatrix, L) or convert(Number, L) can materialize the operator state without changing its mathematical action.

isconstant(L) means repeated calls to update_coefficients[!] are not required to keep L current. islinear(L) means the action is linear in the vector being multiplied; state dependence on (u, p, t) is still allowed for a linear operator.

Keyword Arguments

When an operator accepts keywords during updates, its constructor must record the accepted names with accepted_kwargs, normally as Val((:name1, :name2)). Composite operators forward only those accepted keywords to each component. Extension authors must therefore accept (u, p, t; kwargs...) consistently in every update and application method they advertise. An unlisted keyword must not be silently passed to a leaf update function.

Standard Actions

The behavior of a SciMLOperator is indistinguishable from an AbstractMatrix. These operators can be passed to linear solver packages, and even to ordinary differential equation solvers. The list of overloads to the AbstractMatrix interface includes, but is not limited to, the following:

• Base: size, zero, one, +, -, *, /, \, ∘, inv, adjoint, transpose, convert
• LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize, issymmetric, ishermitian, isposdef
• SparseArrays: sparse, issparse

Multidimensional arrays and batching

SciMLOperator can also be applied to AbstractMatrix subtypes where operator-evaluation is done column-wise.

using LinearAlgebra, SciMLOperators

N = 4
K = 10
L = MatrixOperator(Matrix(I, N, N))
u_mat = rand(N, K)

v_mat = L(u_mat, nothing, nothing, 0.0)
size(v_mat) == (N, K) # true

L can also be applied to AbstractArrays that are not AbstractVecOrMats so long as their size in the first dimension is appropriate for matrix-multiplication. Internally, SciMLOperators reshapes an N-dimensional array to an AbstractMatrix, and applies the operator via matrix-multiplication.

Operator update

This package can also be used to write state-dependent, time-dependent, and parameter-dependent operators, whose state can be updated per a user-defined function. The updates can be done in-place, i.e. by mutating the object, or out-of-place, i.e. in a non-mutating, Zygote-compatible way.

For example,

using LinearAlgebra, SciMLOperators

n = 4
v = rand(n)
u = rand(n)
p = rand(n)
t = rand()

# out-of-place update
mat_update_func = (A, u, p, t) -> t * (p * u')
sca_update_func = (a, u, p, t) -> t * sum(p)

M = MatrixOperator(zeros(n, n); update_func = mat_update_func)
α = ScalarOperator(0.0; update_func = sca_update_func)

L = α * M
L = cache_operator(L, v)

# L is initialized with zero state
L * v == zeros(n) # true

# update operator state with `(u, p, t)`
L = update_coefficients(L, u, p, t)
# and multiply
L * v != zeros(n) # true

# updates state and evaluates L*v at (u, p, t)
L(v, u, p, t) != zeros(n) # true

The out-of-place evaluation function L(v, u, p, t) calls update_coefficients under the hood, which recursively calls the update_func for each component SciMLOperator. Therefore, the out-of-place evaluation function is equivalent to calling update_coefficients followed by Base.*. Notice that the out-of-place evaluation does not return the updated operator.

On the other hand, the in-place evaluation function, L(w, v, u, p, t), mutates L, and is equivalent to calling update_coefficients! followed by mul!. The in-place update behavior works the same way, with a few <!>s appended here and there. For example,

using LinearAlgebra, SciMLOperators

n = 4
w = rand(n)
v = rand(n)
u = rand(n)
p = rand(n)
t = rand()

# in-place update
_A = rand(n, n)
mat_update_func! = (A, u, p, t) -> (copy!(A, _A); lmul!(t, A); nothing)

M = MatrixOperator(zeros(n, n); update_func! = mat_update_func!)

L = M
L = cache_operator(L, v)

# L is initialized with zero state
L * v == zeros(n) # true

# update L in-place
update_coefficients!(L, v, p, t)
# and multiply
mul!(w, L, v) != zeros(n) # true

# updates L in-place, and evaluates w=L*v at (u, p, t)
L(w, v, u, p, t) != zeros(n) # true

The update behavior makes this package flexible enough to be used in OrdinaryDiffEq. As the parameter object p is often reserved for sensitivity computation via automatic-differentiation, a user may prefer to pass in state information via other arguments. For that reason, we allow update functions with arbitrary keyword arguments.

using SciMLOperators

n = 4
v = rand(n)
u = rand(n)
p = rand(n)
t = 0.0
mat_update_func = (A, u, p, t; scale = 0.0) -> scale * (p * u')

M = MatrixOperator(zeros(n, n); update_func = mat_update_func,
    accepted_kwargs = Val((:scale,)))

M(v, u, p, t) == zeros(n) # true
M(v, u, p, t; scale = 1.0) != zeros(n)

DefaultODEAlgorithm(; lazy = Val{true}(), stiffalgfirst = false, kwargs...)

Construct the automatic default ODE algorithm used when solve is called without an explicit algorithm for an ODE problem.

DefaultODEAlgorithm starts with explicit nonstiff methods and switches to stiff methods when stiffness is detected. It chooses among Tsit5, Vern7, Rosenbrock23, Rodas5P, and FBDF, using Krylov-based FBDF for larger stiff systems.

Keywords

• lazy: controls lazy tableau construction for Vern7.
• stiffalgfirst: start on the stiff solver branch when true.
• kwargs...: forwarded to the stiff solver constructors.

Examples

using OrdinaryDiffEqDefault
using SciMLBase: ODEProblem, solve

function f!(du, u, p, t)
    du[1] = -u[1]
    return
end

prob = ODEProblem(f!, [1.0], (0.0, 1.0))
sol = solve(prob, DefaultODEAlgorithm(); reltol = 1.0e-6, abstol = 1.0e-8)

DefaultImplicitODEAlgorithm(; lazy = Val{true}(), stol = 0, ntol = Inf, kwargs...)

Construct the default ODE algorithm with the stiff branch selected first.

This is useful when a problem is expected to be stiff but can still benefit from automatic switching. The nonstiff branch contains Tsit5 and Vern7; the stiff branch contains Rosenbrock23, Rodas5P, and FBDF variants.

Keywords

• lazy: controls lazy tableau construction for Vern7.
• stol: stiffness-detection tolerance passed as stifftol.
• ntol: nonstiff-detection tolerance passed as nonstifftol.
• kwargs...: forwarded to the stiff solver constructors.

Examples

using OrdinaryDiffEqDefault
using SciMLBase: ODEProblem, solve

function f!(du, u, p, t)
    du[1] = -1000u[1]
    return
end

prob = ODEProblem(f!, [1.0], (0.0, 1.0))
sol = solve(prob, DefaultImplicitODEAlgorithm(); reltol = 1.0e-8, abstol = 1.0e-10)