Probabilistic ODE Solvers with Runge-Kutta Means

Runge-Kutta methods are the classic family of solvers for ordinary differential equations (ODEs), and the basis for the state of the art. Like most numerical methods, they return point estimates. We construct a family of probabilistic numerical methods that instead return a Gauss-Markov process defining a probability distribution over the ODE solution. In contrast to prior work, we construct this family such that posterior means match the outputs of the Runge-Kutta family exactly, thus inheriting their proven good properties. Remaining degrees of freedom not identified by the match to Runge-Kutta are chosen such that the posterior probability measure fits the observed structure of the ODE. Our results shed light on the structure of Runge-Kutta solvers from a new direction, provide a richer, probabilistic output, have low computational cost, and raise new research questions.

💡 Research Summary

This paper introduces a novel family of probabilistic ordinary differential equation (ODE) solvers that retain the proven high‑order convergence properties of classical Runge‑Kutta (RK) methods while delivering a full probability distribution over the solution. Traditional RK methods are deterministic point‑estimate algorithms whose coefficients (the Butcher tableau) are carefully chosen so that the numerical solution matches the Taylor expansion of the true solution up to a prescribed order p. Consequently, they enjoy strong stability and error‑control guarantees, but they provide no quantification of uncertainty.

Recent attempts to endow ODE solvers with probabilistic outputs have employed Gaussian process (GP) regression: the ODE’s derivative is treated as noisy observations of the latent solution, and a GP prior on the solution yields a posterior mean and covariance. However, earlier GP‑based solvers selected kernels and evaluation points ad‑hoc, resulting in no guarantee of the desired convergence order; in fact, some configurations fail to even reproduce first‑order (Euler) accuracy.

The authors resolve this gap by constructing Gauss‑Markov‑Runge‑Kutta (GMRK) methods. The key idea is to use integrated Wiener processes of increasing order as GP priors. By choosing the evaluation (or “node”) points appropriately, the posterior mean after s = p observations coincides exactly with the RK estimate of order p. Three theorems are proved:

1. First order (Euler) – A once‑integrated Wiener process (kernel k₁) together with nodes placed at the current posterior mean yields a posterior mean μ(t₀ + h) = x₀ + h y₁, which is precisely Euler’s update.

2. Second order – A twice‑integrated Wiener process (kernel k₂) is employed, but the process must be shifted by a parameter τ that is taken to infinity. In the limit τ → ∞ the posterior mean becomes μ(t₀ + h) = x₀ + h