OrdinaryDiffEqCore API

ODEIntegrator

Fundamental struct allowing interactively stepping through the numerical solving of a differential equation. The full documentation is hosted here: https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/. This docstring describes basic functionality only!

Initialize using integrator = init(prob::ODEProblem, alg; kwargs...). The keyword args which are accepted are the same common solver options used by solve.

For reference, relevant fields of the ODEIntegrator are:

• t - time of the proposed step
• u - value at the proposed step
• opts - common solver options
• alg - the algorithm associated with the solution
• f - the function being solved
• sol - the current state of the solution
• tprev - the last timepoint
• uprev - the value at the last timepoint

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 = 1.0e-9

For more info see the linked documentation page.

AbstractThreadingOption

Abstract supertype for the threading = … option controlling how solvers that expose independent internal work (extrapolation columns, parallel Runge-Kutta stages) execute it. The concrete choices are Sequential, BaseThreads, and PolyesterThreads; isthreaded reports whether a given option enables multithreading.

This is distinct from the thread = … option of the FastBroadcast-based solvers, which takes FastBroadcast.Serial() or FastBroadcast.Threaded().

Sequential() <: AbstractThreadingOption

Threading option that disables internal multithreading — all per-element work runs on a single thread. isthreaded(Sequential()) is false.

BaseThreads() <: AbstractThreadingOption

Threading option that parallelizes internal broadcasting with Julia's built-in Threads.@threads. isthreaded(BaseThreads()) is true.

PolyesterThreads() <: AbstractThreadingOption

Threading option that parallelizes internal broadcasting with Polyester.jl's low-overhead @batch. isthreaded(PolyesterThreads()) is true.

isthreaded(opt) -> Bool

Return whether the threading option opt enables multithreaded internal broadcasting. true for BaseThreads/PolyesterThreads (and true for a Boolopt equal to true), false for Sequential.

Sequential() <: OrdinaryDiffEqCore.AbstractThreadingOption

Use one thread for solver work that can otherwise be executed in parallel. This is the default threading option for deterministic execution.

BaseThreads() <: OrdinaryDiffEqCore.AbstractThreadingOption

Use Julia's built-in Threads.@threads for solver work that can be executed in parallel. The active Julia process must have more than one thread for this to provide parallelism.

PolyesterThreads() <: OrdinaryDiffEqCore.AbstractThreadingOption

Use Polyester.jl's low-overhead threaded execution for solver work that can be executed in parallel. Load Polyester before selecting this option.

AutoAlgSwitch(nonstiffalg, stiffalg; kwargs...)
AutoAlgSwitch(nonstiffalgs::Tuple, stiffalgs::Tuple; kwargs...)

Build a CompositeAlgorithm whose choice function is an AutoSwitch, i.e. an algorithm that automatically switches between a nonstiff and a stiff solver based on runtime stiffness detection. Keyword arguments are forwarded to AutoSwitch.

AutoSwitch(nonstiffalg, stiffalg; kwargs...)

An automatic algorithm switching method that dynamically chooses between a nonstiff and stiff solver based on the problem's stiffness detection. This provides robust performance across a wide range of problems without requiring the user to know the problem's stiffness characteristics a priori.

Arguments

• nonstiffalg: Algorithm to use for nonstiff regions (default: Tsit5())
• stiffalg: Algorithm to use for stiff regions (default: Rodas5P())

Keywords

• maxstiffstep: Maximum number of consecutive steps before switching from nonstiff to stiff (default: 10)
• maxnonstiffstep: Maximum number of consecutive steps before switching from stiff to nonstiff (default: 3)
• nonstifftol: Tolerance for detecting nonstiff behavior (default: 3//4)
• stifftol: Tolerance for detecting stiff behavior (default: 9//10)
• dtfac: Factor for step size adjustment during switches (default: 2.0)
• stiffalgfirst: Whether to start with the stiff algorithm (default: false)
• switch_max: Maximum number of algorithm switches allowed (default: 10)

The switching decision is based on step size rejections and stability estimates.

ssp_coefficient(alg)

Return the SSP coefficient of the ODE algorithm alg. If one time step of size dt with alg can be written as a convex combination of explicit Euler steps with step sizes cᵢ * dt, the SSP coefficient is the minimal value of 1/cᵢ.

Examples

julia> ssp_coefficient(SSPRK104())
6

OrdinaryDiffEqCore.has_stage_limiter(alg)

Trait declaring whether alg's perform_step! applies a stage limiter, i.e. whether it calls integrator.opts.stage_limiter! on each stage value. Defaults to false; a method that implements stage limiting must opt in with

OrdinaryDiffEqCore.has_stage_limiter(::MyAlg) = true

Passing a non-trivial stage_limiter keyword to solve/init for an algorithm whose has_stage_limiter returns false triggers the stage_limiter_unused verbosity toggle, which defaults to ErrorLevel (it errors, rather than silently ignoring the limiter) and can be lowered to WarnLevel/Silent via verbose. (The step_limiter keyword is applied centrally on every accepted step for every method and needs no such opt-in.)

Predictor

Enumeration of per-stage initial-guess strategies for implicit Runge-Kutta methods.

Values

• Predictor.Trivial: Use a zero increment.
• Predictor.Linear: Use a linear extrapolation from the first-same-as-last derivative.
• Predictor.MaxOrder: Use the full previous-step interpolation order.
• Predictor.VariableOrder: Reduce interpolation order as the stage extrapolates farther from the previous step.
• Predictor.CutoffOrder: Use full interpolation order below a cutoff and order one above it.
• Predictor.CopyPrev: Reuse the previous stage derivative.
• Predictor.StageExtrap: Extrapolate recent stage derivatives.
• Predictor.Tableau: Use the tableau-derived stage guess.

NLAnderson(;
    κ = 1 // 100, max_iter = 10, max_history = 5, aa_start = 1,
    droptol = nothing, fast_convergence_cutoff = 1 // 5
)

Anderson-accelerated fixed-point iteration for the implicit stage equations. Like NLFunctional but mixes in max_history previous residuals via a least-squares update to accelerate convergence.

Keywords

• κ, max_iter, fast_convergence_cutoff: as in NLFunctional.
• max_history: number of past iterates kept for the acceleration.
• aa_start: iteration at which acceleration starts.
• droptol: optional condition-number threshold for dropping history columns.

Examples

using OrdinaryDiffEq, OrdinaryDiffEqNonlinearSolve

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
alg = ImplicitEuler(nlsolve = OrdinaryDiffEqNonlinearSolve.NLAnderson())
sol = solve(prob, alg)

NLFunctional(; κ = 1 // 100, max_iter = 10, fast_convergence_cutoff = 1 // 5)

Functional (fixed-point) iteration solver for the implicit stage equations, z ← g(z). No Jacobian/W is formed, so it is cheap per iteration but only converges for mildly stiff problems.

Keywords

• κ: relative tolerance on the increment used in the convergence test.
• max_iter: maximum number of fixed-point iterations per solve.
• fast_convergence_cutoff: convergence-rate threshold for fast convergence.

Examples

using OrdinaryDiffEq, OrdinaryDiffEqNonlinearSolve

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
alg = ImplicitEuler(nlsolve = OrdinaryDiffEqNonlinearSolve.NLFunctional())
sol = solve(prob, alg)

NLNewton(;
    κ = 1 // 100, max_iter = 10, fast_convergence_cutoff = 1 // 5,
    new_W_dt_cutoff = 1 // 5, always_new = false, check_div = true, relax = nothing
)

Quasi-Newton nonlinear solver for the implicit stage equations. Uses the W = M/(γΔt) - J matrix (reused/refactorized across steps and stages when possible) to solve g(z) = 0.

Keywords

• κ: relative tolerance on the Newton increment used in the convergence test.
• max_iter: maximum number of Newton iterations per solve.
• fast_convergence_cutoff: convergence-rate threshold below which convergence is deemed fast (allowing W reuse).
• new_W_dt_cutoff: relative change in γΔt above which W is refactorized.
• always_new: force recomputation of W on every solve.
• check_div: enable early divergence detection.
• relax: optional relaxation parameter in [0, 1) damping the Newton update.

Examples

using OrdinaryDiffEq, OrdinaryDiffEqNonlinearSolve

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
alg = ImplicitEuler(nlsolve = OrdinaryDiffEqNonlinearSolve.NLNewton())
sol = solve(prob, alg)

NonlinearSolveAlg(
    alg = NewtonRaphson(autodiff = AutoFiniteDiff());
    κ = 1 // 100, max_iter = 10, fast_convergence_cutoff = 1 // 5,
    new_W_dt_cutoff = 1 // 5, always_new = false, check_div = true,
    precondition = nothing, postcondition = nothing
)

Use a NonlinearSolve.jl algorithm for the nonlinear stage equations of an implicit OrdinaryDiffEq method. Pass this algorithm as the nlsolve keyword to an implicit solver constructor.

Arguments

• alg: a NonlinearSolve.jl algorithm. The default is finite-difference Newton.

Keywords

• κ: relative tolerance on the Newton increment used in the convergence test.
• max_iter: maximum number of nonlinear iterations per stage solve.
• fast_convergence_cutoff: convergence-rate threshold below which convergence is considered fast and the W matrix may be reused.
• new_W_dt_cutoff: relative change in γΔt above which W is refactorized.
• always_new: force recomputation of W on every nonlinear solve.
• check_div: enable early divergence detection.
• precondition, postcondition: NonlinearSolve.jl's nonlinear preconditioning options, applied to the stage solve. See the section below.

Nonlinear preconditioning of the stage solve

precondition (a left preconditioner G on the residual) and postcondition (an iterate corrector H, the corrector phase of the PCNR method that replaces SPICE-style limiting) are forwarded to the NonlinearSolve.jl solve of each implicit stage. Both are stated in terms of the ODE state at the stage and the ODE parameters, not in terms of the raw unknown of the stage system:

• postcondition(u_stage, u_stage_prev, p, cache). u_stage is the state the stage equations are being solved for — uₙ₊₁ for ImplicitEuler, the stage value tmp + γ⋅z at time t + cΔt for a DIRK stage — so a limiter written for a physical ODE variable applies unchanged. A multi-stage method calls the corrector once per implicit stage, at that stage's own time point, not only at the end of the step. p is the ODE's parameter object and cache is the inner NonlinearSolve cache (nothing for the once-per-stage correction of the predictor). In-place problems overwrite the first argument.
• precondition(fu, u_stage, p). Only u_stage and p are remapped: fu is the stage residual(Δt⋅f(u_stage, p, t + cΔt) - z)/(γΔt), not f itself, and there is no state-space reading of it to map onto.

Internally the stage unknown is the increment z, and the corrector is conjugated with the affine map z ↦ u_stage so that H never sees z. The correction is applied to the stage predictor and then at every iterate the inner solver commits, before its residual is evaluated there. u_stage_prev is the previous iterate of the same stage — the corrected predictor on the stage's first correction — so a limiter that clips relative to the last iterate restarts at each stage rather than reaching back into the previous one.

Two consequences are worth knowing:

• A postcondition that actively corrects makes the stage iteration's displacement larger than the raw Newton step, so a limiter that is still clamping at the end of the iteration shows up as a non-converged stage and the integrator rejects the step and retries with a smaller Δt. That is the intended response.
• precondition changes the residual the inner solver differentiates, so the ODE's W matrix is no longer its Jacobian. W reuse is therefore disabled for the stage solve and the inner solver builds its own Jacobian of the composed residual; this also disables the W-based smoothed error estimate.

On an in-place problem that last point carries a restriction. The in-place stage residual writes the stage state and f's output through preallocated Float64 buffers, and an AutoSpecializef is a FunctionWrapper accepting only Float64 and OrdinaryDiffEq's own one-chunk duals, so it cannot be evaluated at ForwardDiff.Dual. Without W to reuse, the inner solver has to differentiate it, so a precondition on an in-place problem requires an inner algorithm with a non-dual AD backend:

NonlinearSolveAlg(NewtonRaphson(autodiff = AutoFiniteDiff()); precondition = G!)

which is what the default alg already is. Anything ForwardDiff-based is rejected with an ArgumentError rather than allowed to fail inside ForwardDiff. Out-of-place problems have an allocating residual and are unrestricted, and postcondition is unaffected either way — it does not change the residual, so W reuse stays on.

Neither option is supported when the stage system comes from ModelingToolkit's nlstep_data or from a DAEProblem; both throw rather than apply a corrector to an unknown whose relation to the ODE state is not available here. The other nonlinear solvers in this package (NLNewton, NLFunctional, NLAnderson, HomotopyNonlinearSolveAlg) run their own iteration and reject both options rather than ignore them, so nlsolve = NonlinearSolveAlg(...) is how an implicit method gets a corrector.

Examples

using OrdinaryDiffEq, OrdinaryDiffEqNonlinearSolve, NonlinearSolve

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
alg = ImplicitEuler(
    nlsolve = OrdinaryDiffEqNonlinearSolve.NonlinearSolveAlg(NewtonRaphson())
)
sol = solve(prob, alg)

Keeping a positive concentration positive throughout the stage iteration:

H(u, uprev, p, cache) = max(u, zero(u))
alg = ImplicitEuler(
    nlsolve = OrdinaryDiffEqNonlinearSolve.NonlinearSolveAlg(
        NewtonRaphson(); postcondition = H
    )
)

HomotopyNonlinearSolveAlg(
    alg = HomotopySweep(inner = NewtonRaphson(autodiff = AutoFiniteDiff()));
    κ = 1 // 100, max_iter = 10, fast_convergence_cutoff = 1 // 5,
    abstol = nothing, reltol = nothing
)

Solve the implicit stage equations by homotopy continuation in the step size instead of a plain Newton iteration. The stage equation

\[0 = dt⋅f(\mathrm{tmp} + γ⋅z, p, t + c⋅dt) - M z\]

is embedded into the one-parameter family

\[H(z, λ) = λ⋅dt⋅f(\mathrm{tmp} + γ⋅z, p, t + c⋅dt) - M z, \quad λ ∈ [0, 1]\]

(with the analogous embedding tmp + λ⋅f(z) - α/(γ dt)⋅M z for multistep-form methods), built as a SciMLBase.HomotopyProblem and handed to a NonlinearSolve.jl continuation solver. At λ = 0 the solution is known (z = 0 for DIRK-form methods, a single linear solve for multistep form), and the continuation tracks the solution branch of the implicit method as the effective step size grows from 0 to dt — the "principal branch" of Green, Patrick & Spiteri (On theoretical upper limits for valid timesteps of implicit ODE methods, AIMS Mathematics 4(6), 2019). This is much more expensive than NLNewton per stage, but it converges from the exact anchor at λ = 0 regardless of predictor quality, which makes it useful on problems where the Newton iteration fails even after step-size reduction.

For continuation algorithms that implement init/reinit!/solve! (currently HomotopySweep and KantorovichHomotopy), the homotopy problem and solver cache are initialized with the ODE cache and reused across implicit stages. Other homotopy algorithms retain the one-shot solve path.

When the continuation cannot reach λ = 1 (for example at a fold, where the connected solution branch of the stage equation turns back and no consistent solution at the full dt exists), the solve reports divergence and the step is rejected, so the integrator retries with a smaller dt — exactly the semantically correct response to a fold in the step-size homotopy.

Arguments

• alg: the continuation algorithm used for the HomotopyProblem. Defaults to HomotopySweep(inner = NewtonRaphson(autodiff = AutoFiniteDiff())) (natural-parameter continuation). Any NonlinearSolve.jl homotopy solver works, e.g. ArcLengthContinuation(inner = NewtonRaphson(autodiff = AutoFiniteDiff())) to track the branch around folds. Note that the inner solver must not use ForwardDiff-based autodiff: the stage residual closes over preallocated Float64 buffers (the same restriction as NonlinearSolveAlg).

Keywords

• κ, max_iter, fast_convergence_cutoff: kept for interface compatibility with the other nonlinear-solver algorithms; the continuation solve does not run the shared Newton convergence loop, so only max_iter (as a bound on inner iterations reported to the integrator statistics) has an effect.
• abstol, reltol: tolerances passed to the continuation solve. The residual is kept in u units, and abstol = nothing (the default) resolves at solve time to κ * abstol_integrator, mirroring the κ ⋅ tol convergence criterion of the Newton solvers. reltol = nothing uses the NonlinearSolve.jl default.

Examples

using OrdinaryDiffEq, OrdinaryDiffEqNonlinearSolve

prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
alg = ImplicitEuler(
    nlsolve = OrdinaryDiffEqNonlinearSolve.HomotopyNonlinearSolveAlg()
)
sol = solve(prob, alg)