A note on statistical consistency of numerical integrators for multiscale dynamics

A minimal requirement for simulating multiscale systems is to reproduce the statistical behavior of the slow variables. In particular, a good numerical method should accurately appropximate the probability density function (pdf) 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 of the homogenized multiscale 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 multiscale systems driven by fast chaotic dynamics that decorrelate insufficiently rapidly. This suggests that to resolve the statistical behavior of trajectories in certain multiscale systems solvers of suffciently 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.