  Discrete variational methods show excellent performance in numerical
simulations of different mechanical systems. In this paper, we introduce an
iterative procedure for the solution of discrete variational equations for
boundary value problems. More concretely, we explore a parallelization strategy
that leverages the capabilities of multicore CPUs and GPUs (graphics cards). We
study this parallel method for higher-order Lagrangian systems, which appear in
fully-actuated problems and beyond. The most important part of the paper is
devoted to a precise study of different convergence conditions for these
methods. We illustrate their excellent behavior in some interesting examples,
namely Zermelo's navigation problem, a fuel-optimal navigation problem,
interpolation problems or in a fuel optimization problem for a controlled
4-body problem in astrodynamics showing the potential of our method.
