Solvers

ProbNumDiffEq.jl provides two solvers: the EK1 and the EK0. Both based on extended Kalman filtering and smoothing, but the latter relies on evaluating the Jacobian of the vector field.

Which solver should I use?

Use the EK1 to get the best uncertainty quantification and to solve stiff problems.
Use the EK0 to get the fastest runtimes and to solve high-dimensional problems.

All solvers are compatible with DAEs in mass-matrix ODE form. They also specialize on second-order ODEs: If the problem is of type SecondOrderODEProblem, it solves the second-order problem directly; this is more efficient than solving the transformed first-order problem and provides more meaningful posteriors [1].

API

ProbNumDiffEq.EK1 — Type

EK1(; order=3,
      smooth=true,
      prior=IWP(order),
      diffusionmodel=DynamicDiffusion(),
      initialization=TaylorModeInit(num_derivatives(prior)),
      kwargs...)

Gaussian ODE filter with first-order vector field linearization.

This is a semi-implicit, L-stable ODE solver so it can handle stiffness quite well [4], and it generally produces more expressive posterior covariances than the EK0. However, as typical implicit ODE solvers it scales cubically with the ODE dimension [5], so if you're solving a high-dimensional non-stiff problem you might want to give the EK0 a try.

Arguments

order::Integer: Order of the integrated Wiener process (IWP) prior.
smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
prior::AbstractGaussMarkovProcess: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Some additional kwargs relating to implicit solvers are supported; check out DifferentialEquations.jl's Extra Options page. Right now, we support autodiff, chunk_size, and diff_type. In particular, autodiff=false can come in handy to use finite differences instead of ForwardDiff.jl to compute Jacobians.

Examples

julia> solve(prob, EK1())

References

source

ProbNumDiffEq.EK0 — Type

EK0(; order=3,
      smooth=true,
      prior=IWP(order),
      diffusionmodel=DynamicDiffusion(),
      initialization=TaylorModeInit(num_derivatives(prior)))

Gaussian ODE filter with zeroth-order vector field linearization.

This is an explicit ODE solver. It is fast and scales well to high-dimensional problems [5], but it is not L-stable [4]. So for stiff problems, use the EK1.

Whenever possible this solver will use a Kronecker-factored implementation to achieve its linear scaling and to get the best runtimes. This can currently be done only with an IWP prior (default), with a scalar diffusion model (either DynamicDiffusion or FixedDiffusion). For other configurations the solver falls back to a dense implementation which scales cubically with the problem size.

Arguments

order::Integer: Order of the integrated Wiener process (IWP) prior.
smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
prior::AbstractGaussMarkovProcess: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Examples

julia> solve(prob, EK0())

References

source

Probabilistic Exponential Integrators

ProbNumDiffEq.ExpEK — Function

ExpEK(; L, order=3, kwargs...)

Probabilistic exponential integrator

Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs that provide improved stability by essentially solving the linear part of the ODE exactly. In probabilistic numerics, this amounts to including the linear part into the prior model of the solver.

ExpEK is therefore just a short-hand for EK0 with IOUP prior:

ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)

See also RosenbrockExpEK, EK0, EK1.

Arguments

See EK0 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
julia> solve(prob, ExpEK(L=-1))

Reference

[2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)

source

ProbNumDiffEq.RosenbrockExpEK — Function

RosenbrockExpEK(; order=3, kwargs...)

Probabilistic Rosenbrock-type exponential integrator

A probabilistic exponential integrator similar to ExpEK, but with automatic linearization along the mean numerical solution. This brings the advantage that the linearity does not need to be specified manually, and the more accurate local linearization can sometimes also improve stability; but since the "prior" is adjusted at each step the probabilistic interpretation becomes more complicated.

RosenbrockExpEK is just a short-hand for EK1 with locally-updated IOUP prior:

RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)

References

[4]: F. Tronarp, H. Kersting, S. Särkkä and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: a new perspective. Statistics and Computing 29, 1297–1315 (2019).
[5]: N. Krämer, N. Bosch, J. Schmidt and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 11634–11649.
[2]: N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).