ODE-Solvers from ScratchΒΆ

All the other tutorials show how to use the ODE-solver with the probsolve_ivp function. This is great, though probnum has more customisation to offer.

import probnum.diffeq as pnd
import probnum.filtsmooth as pnfs
import probnum.random_variables as pnrv
import numpy as np

import matplotlib.pyplot as plt

First we define the ODE problem. As always, we use Lotka-Volterra. It is important that initrv is an actual random variable. Arrays can be wrapped with probnum.random_variables.Constant.

ivp = pnd.lotkavolterra([0.0, 20.0], initrv=pnrv.Constant(np.array([20.0, 20.0])))

Next, we define a prior distribution and a measurement model. The former allows any Integrator, which currently restricts the choice to IBM, IOUP, and Matern. The measurement model requires a choice between EK0, EK1 (extended Kalman filters of order 0 or 1, respectively) and UK (unscented Kalman filter). After this choice is made, each of those classes offer a constructor that defines the model directly from the initial value problem, e.g. DiscreteUKFComponent.from_ode(ivp, prior, evlvar). evlvar is usually zero. The prior is the Integrator chosen above.

prior = pnfs.statespace.IBM(ordint=4, spatialdim=ivp.dimension, diffconst=1.0)
ekf = pnfs.DiscreteEKFComponent.from_ode(ivp, prior, np.zeros((2, 2)), ek0_or_ek1=0)

Next, we turn this into a Kalman object which will do all the heavy lifting for the ODE solver. To this end, we need to define an initial distribution. In probsolve_ivp, this is automatised by using as many derivatives as the IVP object allows (max. 2).

This time, assume that some bird told us about the real derivatives of the initial values (stay tuned for more on this). This can be used to define a more robust ODE solver.

initial_values = np.array(
initrv = pnrv.Normal(
    initial_values, np.zeros((len(initial_values), len(initial_values)))
kalman = pnfs.Kalman(prior, ekf, initrv)
solver = pnd.GaussianIVPFilter(ivp, kalman, with_smoothing=True)

Now we can solve the ODE. To this end, define a StepRule, e.g. ConstantSteps or AdaptiveSteps.

steprule = pnd.ConstantSteps(0.1)
odesol = solver.solve(steprule=steprule)

GaussianIVPFilter.solve returns an ODESolution object, which is sliceable and callable. The latter can be used to plot the solution on a uniform grid, even though the solution was computed on an adaptive grid. Be careful: the return values of __call__, etc., are always random variable-like objects. We decide to plot the mean.

evalgrid = np.arange(ivp.t0, ivp.tmax, step=0.1)
y = odesol(evalgrid).mean

Et voila: this is the solution to the Lotka-Volterra model.

plt.plot(evalgrid, y, "o-", linewidth=1)