Choose the right solver

Good solvers are problem-dependent. However, some guidelines exist:

Problem characteristics

Choosing the right approach matters because problem size and behaviour directly impact solver efficiency, stability, and the accuracy of the uncertainty quantification.

Dimensionality: For low-dimensional problems, use dense covariances, which track full correlations between state variables and offer the best stability and uncertainty quantification. For larger problems, use blockdiagonal or isotropic state-space models, which are more efficient by tracking only partial uncertainty correlations. However, their uncertainty quantification is typically worse. The general trade-off is between accuracy and speed: dense models scale cubically in the dimension but provide the best accuracy; the other two models scale linearly in the dimension.

Stiffness: Stiff problems have rapid changes or very different timescales. For these, use dense state-space models with first-order linearization. See also the prior recommendations below. Avoid zeroth-order methods and isotropic state-space models for stiff problems. Block-diagonal state-space models with first-order linearization may suffice for moderately stiff cases, but expect that all solvers except first-order linearisation in dense state-space models have worse stability than, for example, implicit Runge-Kutta methods.

Filters vs smoothers

Choosing between filters and smoothers matters because it balances computational cost with the accuracy of uncertainty estimates across the trajectory. Use fixed-point smoothing for adaptive timestepping and fixed-interval smoothing for fixed timestepping. When only computing the terminal value of a differential equation, choose a filter.

Calibration

Calibration matters not just because it ensures uncertainty estimates reflect the real error, but also because it can considerably affect adaptive time-stepping. Use dynamic calibration when output scales vary significantly, for example, in \(u'(t) = 10u(t)\), \(0 \leq t \leq 10\). Use maximum-likelihood calibration for other cases. Remove automatic calibration for parameter estimation. When solving multidimensional problems where each dimension has a different magnitude, adjust the output scale of the prior manually before the simulation.

Prior distributions

Prior distributions encode assumptions about solution dynamics. The default prior is the integrated Wiener process; however, integrated Ornstein-Uhlenbeck processes work well for discretised semilinear partial differential equations (and other semilinear problems), especially in fixed-step simulations. For other needs, consult:

Bosch, Nathanael, Philipp Hennig, and Filip Tronarp. "Probabilistic exponential integrators." Advances in Neural Information Processing Systems 36 (2023): 40450-40467.

Adjust the base-output-scales of the Wiener process if state variables have vastly different magnitudes, like in the Robertson problem, where one dimension is \(10^{5}\) times smaller than the other.

Number of Taylor coefficients (''order'')

The number of Taylor coefficients trades off accuracy against computational cost. Use 4-5 for most problems, 7-8 when simulating with tolerances close to machine precision, and 2-3 in low-precision arithmetic (for instance, on a GPU).

Error estimation

In adaptive time-stepping, there exist different error estimates. The default for solving (explicit) ODEs is the residual-based one, because it has proven effective over many years. When solving implicit differential equations (like DAEs), use the state-based estimate instead because the constraints may live on arbitrary scales, which the residual-based method struggles with. For error-normalisation, use scale-then-RMS when applicable, and RMS-then-scale only if necessary.

For controllers, the choice does not matter much. Integral controllers seem to be slightly more effective for most problems except for stiff ODEs, where proportional-integral controllers work better. Try all and report back with any insights.

Summary: Choosing a solver

For beginners: Start with integrated Wiener processes and four Taylor coefficients, fixed-point smoothing, first-order linearization, and dense state-space models. For high-dimensional problems, use zeroth-order linearization and block-diagonal state-space models. For parameter estimation, use fixed steps with a fixed-interval smoother.

For advanced users: Use the guidelines above based on your problem's dimensionality, stiffness, and requirements. Consult the examples in the documentation, and reach out with any questions.