Abstract. Often, Lagrange's formalism based on the Lagrangian is the preferred way in order to derive equations of motion for a mechanical system. Therefore, the first step is the formulation of the total kinetic and total potential energy. In this formalism, holonomic constraints can be taken into account by using the Lagrange multiplier method, and external and dissipative forces can be included variationally consistently by means of the Lagrange-D'Alembert principle. In this way, we are directly arrive at a weak formulation of the equations of motion without using a strong form and integration by parts. The variational time integration method is now based on the discretization of the action functional and the discrete variational calculus. This method leads to time stepping schemes called variational integrators.In this paper, we show that variational integrators based on higher-order shape functions and quadrature rules leads to a higher-order time approximation, which preserve the total linear and total angular momentum of the mechanical system. A further topic of the paper is a new implementation of the Lagrange multiplier method for holonomic constraints in the higher-order variational integrator, in order to compute bearing forces by means of Lagrange multipliers. Usually, variational integrators show an excellent long time behavior and bind the total energy error per time step, but are not able to preserve the total energy exactly. Therefore, we introduce a discrete gradient of the total potential energy in the variational integrator to preserve the total energy. We show different discrete gradients with different results for numerical stability. We can show that these time stepping schemes preserve the total linear momentum, total angular momentum as well as the total energy for every time step size. As numerical examples, we show motions of three-dimensional continua, which are discretized in space by the finite element method. We start with a free flying hyper-elastic rotor to show the preservation of total linear and total angular momentum in combination with higher-order accuracy in time. In the second example, we consider a hyper-elastic beam and apply the Lagrange multiplier method in order to calculate the bearing forces at fixed nodes. In the last example, the energy conservation is shown for different discrete gradients.