perturbsolve_ivp

probnum.diffeq.perturbsolve_ivp(f, t0, tmax, y0, rng, method='RK45', perturb='step-lognormal', noise_scale=10.0, adaptive=True, atol=1e-06, rtol=0.001, step=None, time_stops=None)[source]

Solve an initial value problem with a perturbation-based ODE solver.

Parameters

  • f (Callable) – ODE vector field.

  • t0 (FloatLike) – Initial time point.

  • tmax (FloatLike) – Final time point.

  • y0 (ArrayLike) – Initial value.

  • rng (np.random.Generator) – Random number generator.

  • method (str) –

    Integration method to use. The following are available (docs adapted from https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html):

    • RK45 (default): Explicit Runge-Kutta method of order 5(4) 2. The error is controlled assuming accuracy of the fourth-order method, but steps are taken using the fifth-order accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output 3. Can be applied in the complex domain.

    • RK23: Explicit Runge-Kutta method of order 3(2) 4. The error is controlled assuming accuracy of the second-order method, but steps are taken using the third-order accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.

    Other integrators are not supported currently.

  • perturb (str) –

    Which perturbation style to use. Currently, the following methods are supported:

    • step-lognormal: Perturbed-step (aka random time-step numerical integration) method with lognormally distributed perturbation of the step-size 1.

    • step-uniform: Perturbed-step (aka random time-step numerical integration) method with uniformly distributed perturbation of the step-size 1.

    Other integrators are not supported currently.

  • noise_scale (FloatLike) – Scale of the perturbation. Optional. Default is 10.0. The magnitude of this parameter significantly impacts the width of the posterior.

  • adaptive (bool) – Whether to use adaptive steps or not. Default is True.

  • atol (FloatLike) – Absolute tolerance of the adaptive step-size selection scheme. Optional. Default is 1e-6.

  • rtol (FloatLike) – Relative tolerance of the adaptive step-size selection scheme. Optional. Default is 1e-3.

  • step (Optional[FloatLike]) – Step size. If atol and rtol are not specified, this step-size is used for a fixed-step ODE solver. If they are specified, this only affects the first step. Optional. Default is None, in which case the first step is chosen as prescribed by propose_firststep().

  • time_stops (Optional[ArrayLike]) – Time-points through which the solver must step. Optional. Default is None.

Raises

  • ValueError – If the ‘method’ string does not correspond to a supported method.

  • ValueError – If the ‘perturb’ string does not correspond to a supported perturbation style.

Returns

ODE Solution of the perturbed-step-solver.

Return type

perturbed.step.PerturbedStepSolution

Examples

>>> from probnum.diffeq import perturbsolve_ivp
>>> import numpy as np

Solve a simple logistic ODE with fixed steps. Per default, perturbsolve_ivp uses a perturbed-step solver with lognormal perturbation.

>>> rng = np.random.default_rng(seed=2)
>>>
>>> def f(t, x):
...     return 4*x*(1-x)
>>>
>>> y0 = np.array([0.15])
>>> t0, tmax = 0., 1.5
>>> solution = perturbsolve_ivp(
...     f, t0, tmax, y0, rng=rng, step=0.25, method="RK23", adaptive=False
... )
>>> print(np.round(solution.states.mean, 3))
[[0.15 ]
 [0.325]
 [0.56 ]
 [0.772]
 [0.893]
 [0.964]
 [0.989]]

Each solution is the result of a randomly-perturbed call of an underlying Runge-Kutta solver. Therefore, if you call it again, the output will be different:

>>> other_solution = perturbsolve_ivp(
...     f, t0, tmax, y0, rng=rng, step=0.25, method="RK23", adaptive=False
... )
>>> print(np.round(other_solution.states.mean, 3))
[[0.15 ]
 [0.319]
 [0.57 ]
 [0.785]
 [0.908]
 [0.968]
 [0.989]]

Other methods, such as RK45 (as well as other perturbation styles) are easily accessible. Let us solve the same equation, with an adaptive RK45 solver that uses uniformly perturbed steps.

>>> solution = perturbsolve_ivp(
...     f, t0, tmax, y0, rng=rng, atol=1e-5, rtol=1e-6,
...     method="RK45", perturb="step-uniform", adaptive=True
... )
>>> print(np.round(solution.states.mean, 3))
[[0.15 ]
 [0.152]
 [0.167]
 [0.26 ]
 [0.431]
 [0.646]
 [0.849]
 [0.883]
 [0.915]
 [0.953]
 [0.976]
 [0.986]]

References

1(1,2)

Abdulle, A. and Garegnani, G.. Random time step probabilistic methods for uncertainty quantification in chaotic and geometric numerical integration. Statistics and Computing. 2020.

2

J. R. Dormand, P. J. Prince.. A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.

3

L. W. Shampine. Some Practical Runge-Kutta Formulas. Mathematics of Computation, Vol. 46, No. 173, pp. 135-150, 1986.

4

P. Bogacki, L.F. Shampine. A 3(2) Pair of Runge-Kutta Formulas. Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.