RungeKuttaInitialization¶

class probnum.diffeq.odefilter.initialization_routines.RungeKuttaInitialization(dt=0.01, method='DOP853')

Initialize a probabilistic ODE solver by fitting the prior process to a few steps of an approximate ODE solution computed with Scipy’s Runge-Kutta methods.

Parameters

Examples

```>>> import numpy as np
>>> from probnum.randvars import Normal
>>> from probnum.problems.zoo.diffeq import vanderpol
>>> from probnum.randprocs.markov.integrator import IntegratedWienerProcess
```

Compute the initial values of the van-der-Pol problem as follows. First, we set up the ODE problem and the prior process.

```>>> ivp = vanderpol()
>>> print(ivp.y0)
[2. 0.]
>>> prior_process = IntegratedWienerProcess(initarg=ivp.t0, num_derivatives=3, wiener_process_dimension=2)
```

Next, we call the initialization routine.

```>>> rk_init = RungeKuttaInitialization()
>>> improved_initrv = rk_init(ivp=ivp, prior_process=prior_process)
>>> 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]
```

Attributes Summary

 `is_exact` Exactness of the computed initial values. `requires_jax` Whether the implementation of the routine relies on JAX.

Methods Summary

 `__call__`(ivp, prior_process) Compute the initial distribution.

Attributes Documentation

is_exact

Exactness of the computed initial values.

Some initialization routines yield the exact initial derivatives, some others only yield approximations.

Return type

`bool`

requires_jax

Whether the implementation of the routine relies on JAX.

Return type

`bool`

Methods Documentation

__call__(ivp, prior_process)[source]

Compute the initial distribution.

For Runge-Kutta initialization, 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
Returns

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

Return type

Normal