initialize_odefilter_with_rk

probnum.diffeq.initialize_odefilter_with_rk(f, y0, t0, prior_process, df=None, h0=0.01, method='DOP853')[source]

Initialize an ODE filter by fitting the prior process to a few steps of an approximate ODE solution computed with Scipy’s RK.

It goes as follows:

1. The ODE integration problem is set up on the interval [t0, t0 + (2*order+1)*h0] and solved with a call to scipy.integrate.solve_ivp. The solver is uses adaptive steps with atol=rtol=1e-12, but is forced to pass through the events (t0, t0+h0, t0 + 2*h0, ..., t0 + (2*order+1)*h0). The result is a vector of time points and states, with at least (2*order+1). Potentially, the adaptive steps selected many more steps, but because of the events, fewer steps cannot have happened.

  1. A prescribed prior is fitted to the first (2*order+1) (t, y) pairs of the solution. order is the order of the prior.

3. The value of the resulting posterior at time t=t0 is an estimate of the state and all its derivatives. The resulting marginal standard deviations estimate the error. This random variable is returned.

Parameters

  • f – ODE vector field.

  • y0 – Initial value.

  • t0 – Initial time point.

  • prior_process – Prior Gauss-Markov process used for the ODE solver. For instance an integrated Brownian motion prior (IBM).

  • df – Jacobian of the ODE vector field. Optional. If specified, more components of the result will be exact.

  • h0 – Maximum step-size to use for computing the approximate ODE solution. The smaller, the more accurate, but also, the smaller, the less stable. The best value here depends on the ODE problem, and probably the chosen method. Optional. Default is 1e-2.

  • method – Which solver to use. This is communicated as a string that is compatible with scipy.integrate.solve_ivp(..., method=method). Optional. Default is DOP853.

Returns

Estimated (improved) initial random variable. Compatible with the specified prior.

Return type

Normal

Examples

>>> from dataclasses import astuple
>>> from probnum.randvars import Normal
>>> from probnum.statespace import IBM
>>> from probnum.problems.zoo.diffeq import vanderpol
>>> from probnum.randprocs import MarkovProcess

Compute the initial values of the van-der-Pol problem as follows

>>> f, t0, tmax, y0, df, *_ = astuple(vanderpol())
>>> print(y0)
[2. 0.]
>>> prior = IBM(ordint=3, spatialdim=2)
>>> initrv = Normal(mean=np.zeros(prior.dimension), cov=np.eye(prior.dimension))
>>> prior_process = MarkovProcess(transition=prior, initrv=initrv, initarg=t0)
>>> improved_initrv = initialize_odefilter_with_rk(f, y0, t0, prior_process=prior_process, df=df)
>>> print(prior_process.transition.proj2coord(0) @ improved_initrv.mean)
[2. 0.]
>>> print(np.round(improved_initrv.mean, 1))
[    2.      0.     -2.     58.2     0.     -2.     60.  -1745.7]
>>> print(np.round(np.log10(improved_initrv.std), 1))
[-13.8 -11.3  -9.   -1.5 -13.8 -11.3  -9.   -1.5]