A note on statistical consistency of numerical integrators for multi-scale dynamics
For researchers simulating multi-scale dynamical systems, this work quantifies statistical consistency of numerical integrators and provides a practical correction for bias.
This note shows that numerical simulations of multi-scale systems converge to the probability density function of homogenized modified equations rather than the exact homogenized system, and that statistical bias is exacerbated by insufficiently rapid decorrelation of fast chaotic dynamics. A second-order Taylor method is proposed to correct the lowest-order bias in Euler's method.
A minimal requirement for simulating multi-scale systems is to reproduce the statistical behavior of the slow variables. In particular, a good numerical method should accurately aproximate the probability density function of the continuous-time slow variables. In this note we use results from homogenization and from backward error analysis to quantify how errors of time integrators affect the mean behavior of trajectories. We show that numerical simulations converge, not to the exact probability density function (pdf) of the homogenized multi-scale system, but rather to that of the homogenized modified equations following from backward error analysis. Using homogenization theory we find that the observed statistical bias is exacerbated for multi-scale systems driven by fast chaotic dynamics that decorrelate insufficiently rapidly. This suggests that to resolve the statistical behavior of trajectories in certain multi-scale systems solvers of sufficiently high order are necessary. Alternatively, backward error analysis suggests the form of an amended vector field that corrects the lowest order bias in Euler's method. The resulting scheme, a second order Taylor method, avoids any statistical drift bias. We corroborate our analysis with a numerical example.