# Adaptive step-size selection for ODE filters¶

In ProbNum the step size adaptation scheme from Schober et al., 2019 respectively Bosch et al., 2021 is implemented.

[1]:

# Make inline plots vector graphics instead of raster graphics
%matplotlib inline
from IPython.display import set_matplotlib_formats

set_matplotlib_formats("pdf", "svg")

# Plotting
import matplotlib.pyplot as plt

plt.style.use("../../probnum.mplstyle")

[2]:

import numpy as np
from probnum import randvars
from probnum.diffeq import probsolve_ivp


We begin by defining an ODE: the usual Lotka-Volterra problem.

[3]:

def f(t, y):
y1, y2 = y
return np.array([0.5 * y1 - 0.05 * y1 * y2, -0.5 * y2 + 0.05 * y1 * y2])

def df(t, y):
y1, y2 = y
return np.array([[0.5 - 0.05 * y2, -0.05 * y1], [0.05 * y2, -0.5 + 0.05 * y1]])

t0 = 0.0
tmax = 20.0
y0 = np.array([20, 20])


## IBM(1) with EK0¶

We start by showing how solving an IVP with adaptive steps works.

[4]:

sol = probsolve_ivp(f, t0, tmax, y0, df=df, algo_order=4, method="EK1")

evalgrid = np.linspace(t0, tmax, 200)
means = sol(evalgrid).mean

plt.plot(evalgrid, means[:, 0], label="Prey")
plt.plot(evalgrid, means[:, 1], label="Predators")
plt.legend(loc="upper right")
plt.show()


Lets visualise the individual steps.

[5]:

plt.plot(evalgrid, means[:, 0], label="Prey")
plt.plot(evalgrid, means[:, 1], label="Predators")
for t in sol.locations:
plt.axvline(t, linewidth=0.5, color="gray")
plt.legend(loc="upper right")
plt.show()


Note how more steps are taken near the peaks.

### The same, for a few priors¶

Let’s consider how other priors fare in this setting.

[6]:

algo_orders = [2, 3]
filters = ["EK0", "EK1"]

fig, ax = plt.subplots(
nrows=len(algo_orders),
ncols=len(filters),
sharex=True,
sharey=True,
figsize=(12, 12),
)

evalgrid = np.linspace(t0, tmax, 100)

for idx in range(len(algo_orders)):  # algo_orders in rows
for jdx in range(len(filters)):  # filters in cols
sol = probsolve_ivp(
f,
t0,
tmax,
y0,
df=df,
algo_order=algo_orders[idx],
method=filters[jdx],
)

solution = sol(evalgrid)
ts, means, stds = evalgrid, solution.mean, solution.std
ax[idx][jdx].plot(ts, means)

for t in sol.locations:
ax[idx][jdx].axvline(t, linewidth=0.2, color="black")
ax[idx][jdx].set_title(f"Order {algo_orders[idx]} w/{filters[jdx]}")
ax[idx][jdx].set_xlim((t0, tmax))
ax[idx][jdx].set_ylim((-10, 35))

plt.show()


They all seem to capture the true solution fairly well. The higher order takes significantly fewer steps than the lower order.