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 GaussMarkov 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 RungeKutta type.
This function turns a priorstring into an
ODEPrior
, a methodstring into a filter/smoother of classGaussFiltSmooth
, creates aGaussianIVPFilter
object and calls thesolve()
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 zeroth order Taylor approximation (EKF0), extended Kalman filtering (EKF1), unscented Kalman filtering (UKF), extended Kalman smoothing based on a zeroth 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
andtol
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
andtol
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 areExtended 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 RHSvector 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:
 t
np.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\).
 y
list
ofRandomVariable
, length=N Discretetime solution at times \(t_1, ..., t_N\), as a list of random variables. The means and covariances can be accessed with
solution.y.mean
andsolution.y.cov
.
 t
 Return type
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 nonlinear 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]]