  We study numerical methods for the generalized Langevin equation (GLE) with a
positive Prony series memory kernel, in which case the GLE can be written in an
extended variable Markovian formalism. We propose a new splitting method that
is easy to implement and is able to substantially improve the accuracy and
robustness of GLE simulations in a wide range of the parameters. An error
analysis is performed in the case of a one-dimensional harmonic oscillator,
revealing that several averages are exact for the newly proposed method.
Various numerical experiments in both equilibrium and nonequilibrium
simulations are also conducted to compare the method with popular alternatives
in interacting multi-particle systems.
