  In the last two decades, significant effort has been put in understanding and
designing so-called structure-preserving numerical methods for the simulation
of mechanical systems. Geometric integrators attempt to preserve the geometry
associated to the original system as much as possible, such as the structure of
the configuration space, the energy behaviour, preservation of constants of the
motion and of constraints or other structures associated to the continuous
system (symplecticity, Poisson structure...). In this article, we develop
high-order geometric (or pseudo-variational) integrators for nonholonomic
systems, i.e., mechanical systems subjected to constraint functions which are,
roughly speaking, functions on velocities that are not derivable from position
constraints. These systems realize rolling or certain kinds of sliding contact
and are important for describing different classes of vehicles.
