  In this paper we construct higher-order variational integrators for a class
of degenerate systems described by Lagrangians that are linear in velocities.
We analyze the geometry underlying such systems and develop the appropriate
theory for variational integration. Our main observation is that the evolution
takes place on the primary constraint and the 'Hamiltonian' equations of motion
can be formulated as an index-1 differential-algebraic system. We also
construct variational Runge-Kutta methods and analyze their properties. The
general properties of Runge-Kutta methods depend on the 'velocity' part of the
Lagrangian. If the 'velocity' part is also linear in the position coordinate,
then we show that non-partitioned variational Runge-Kutta methods are
equivalent to integration of the corresponding first-order Euler-Lagrange
equations, which have the form of a Poisson system with a constant structure
matrix, and the classical properties of the Runge-Kutta method are retained. If
the 'velocity' part is nonlinear in the position coordinate, we observe a
reduction of the order of convergence, which is typical of numerical
integration of DAEs. We verify our results through numerical experiments for
various dynamical systems.
