probsolve_ivp

probnum.diffeq.probsolve_ivp(ivp, method='ekf0', which_prior='ibm1', atol=None, rtol=None, step=None, firststep=None, **kwargs)[source]

Solve initial value problem with Gaussian filtering and smoothing.

Numerically computes a Gauss-Markov process which solves numerically the initial value problem (IVP) based on a system of first order ordinary differential equations (ODEs)

\[\dot x(t) = f(t, x(t)), \quad x(t_0) = x_0, \quad t \in [t_0, T]\]

by regarding it as a (nonlinear) Gaussian filtering (and smoothing) problem 3. For some configurations it recovers certain multistep methods 1. Convergence rates of filtering 2 and smoothing 4 are comparable to those of methods of Runge-Kutta type.

This function turns a prior-string into an ODEPrior, a method-string into a filter/smoother of class GaussFiltSmooth, creates a GaussianIVPFilter object and calls the solve() method. For advanced usage we recommend to do this process manually which enables advanced methods of tuning the algorithm.

This function supports the methods: extended Kalman filtering based on a zero-th order Taylor approximation (EKF0), extended Kalman filtering (EKF1), unscented Kalman filtering (UKF), extended Kalman smoothing based on a zero-th order Taylor approximation (EKS0), extended Kalman smoothing (EKS1), and unscented Kalman smoothing (UKS).

Parameters
  • ivp (IVP) – Initial value problem to be solved.

  • step (float) – Step size \(h\) of the solver. This defines the discretisation mesh as each proposed step is equal to \(h\) and all proposed steps are accepted. Only one of out of step and tol is set.

  • tol (float) – Tolerance \(\varepsilon\) of the adaptive step scheme. We implement the scheme proposed by Schober et al., accepting a step if the absolute as well as the relative error estimate are smaller than the tolerance, \(\max\{e, e / |y|\} \leq \varepsilon\). Only one of out of step and tol is set.

  • which_prior (str, optional) –

    Which prior is to be used. Default is an IBM(1), further support for IBM(\(q\)), IOUP(\(q\)), Matern(\(q+1/2\)), \(q\in\{1, 2, 3, 4\}\) is provided. The available options are

    IBM(\(q\))

    'ibm1', 'ibm2', 'ibm3', 'ibm4'

    IOUP(\(q\))

    'ioup1', 'ioup2', 'ioup3', 'ioup4'

    Matern(\(q+0.5\))

    'matern32', 'matern52', 'matern72', 'matern92'

    The type of prior relates to prior assumptions about the derivative of the solution. The IBM(\(q\)) prior leads to a \(q\)-th order method that is recommended if little to no prior information about the solution is available. On the other hand, if the \(q\)-th derivative is expected to regress to zero, an IOUP(\(q\)) prior might be suitable.

  • method (str, optional) –

    Which method is to be used. Default is ekf0 which is the method proposed by Schober et al.. The available options are

    Extended Kalman filtering/smoothing (0th order)

    'ekf0', 'eks0'

    Extended Kalman filtering/smoothing (1st order)

    'ekf1', 'eks1'

    Unscented Kalman filtering/smoothing

    'ukf', 'uks'

    First order extended Kalman filtering and smoothing methods require Jacobians of the RHS-vector field of the IVP. The uncertainty estimates as returned by EKF1/S1 and UKF/S appear to be more reliable than those of EKF0/S0. The latter is more stable when it comes to very small steps.

  • firststep (float, optional) – First suggested step \(h_0\) for adaptive step size scheme. Default is None which lets the solver start with the suggestion \(h_0 = T - t_0\). For low accuracy it might be more efficient to start out with smaller \(h_0\) so that the first acceptance occurs earlier.

Returns

solution – Solution of the ODE problem.

Contains fields:

tnp.ndarray, shape=(N,)

Mesh used by the solver to compute the solution. It includes the initial time \(t_0\) but not necessarily the final time \(T\).

ylist of RandomVariable, length=N

Discrete-time solution at times \(t_1, ..., t_N\), as a list of random variables. The means and covariances can be accessed with solution.y.mean and solution.y.cov.

Return type

ODESolution

See also

GaussianIVPFilter()

Solve IVPs with Gaussian filtering and smoothing

ODESolution()

Solution of ODE problems

References

1

Schober, M., Särkkä, S. and Hennig, P.. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing, 2019.

2

Kersting, H., Sullivan, T.J., and Hennig, P.. Convergence rates of Gaussian ODE filters. 2019.

3

Tronarp, F., Kersting, H., Särkkä, S., and Hennig, P.. Probabilistic solutions to ordinary differential equations as non-linear Bayesian filtering: a new perspective. Statistics and Computing, 2019.

4

Tronarp, F., Särkkä, S., and Hennig, P.. Bayesian ODE solvers: the maximum a posteriori estimate. 2019.

Examples

>>> from probnum.diffeq import logistic, probsolve_ivp
>>> from probnum import random_variables as rvs
>>> import numpy as np
>>> initrv = rvs.Constant(0.15)
>>> ivp = logistic(timespan=[0., 1.5], initrv=initrv, params=(4, 1))
>>> solution = probsolve_ivp(ivp, method="ekf0", step=0.1)
>>> print(np.round(solution.y.mean, 2))
[[0.15]
 [0.21]
 [0.28]
 [0.36]
 [0.46]
 [0.56]
 [0.65]
 [0.74]
 [0.81]
 [0.86]
 [0.9 ]
 [0.93]
 [0.95]
 [0.97]
 [0.98]
 [0.98]]
>>> initrv = rvs.Constant(0.15)
>>> ivp = logistic(timespan=[0., 1.5], initrv=initrv, params=(4, 1))
>>> solution = probsolve_ivp(ivp, method="eks1", which_prior="ioup3", step=0.1)
>>> print(np.round(solution.y.mean, 2))
[[0.15]
 [0.21]
 [0.28]
 [0.37]
 [0.47]
 [0.57]
 [0.66]
 [0.74]
 [0.81]
 [0.87]
 [0.91]
 [0.93]
 [0.96]
 [0.97]
 [0.98]
 [0.99]]