  A novel probabilistic numerical method for quantifying the uncertainty
induced by the time integration of ordinary differential equations (ODEs) is
introduced. Departing from the classical strategy to randomize ODE solvers by
adding a random forcing term, we show that a probability measure over the
numerical solution of ODEs can be obtained by introducing suitable random
time-steps in a classical time integrator. This intrinsic randomization allows
for the conservation of geometric properties of the underlying deterministic
integrator such as mass conservation, symplecticity or conservation of first
integrals. Weak and mean-square convergence analysis are derived. We also
analyse the convergence of the Monte Carlo estimator for the proposed random
time step method and show that the measure obtained with repeated sampling
converges in the mean-square sense independently of the number of samples.
Numerical examples including chaotic Hamiltonian systems, chemical reactions
and Bayesian inferential problems illustrate the accuracy, robustness and
versatility of our probabilistic numerical method.
