Los integradores variacionales para sistemas dinámicos con restricciones holonómicas se proponen basados en el principio de Hamilton. El principio variacional se discretiza aproximando las coordenadas generalizadas y los multiplicadores de Lagrange mediante polinomios de Lagrange, aproximando las integrales mediante reglas de cuadratura. Mientras tanto, se definen puntos de restricción para discretizar las restricciones holonómicas. La función del principio variacional se divide en dos partes, es decir, la acción del término no restringido y el término restringido y las acciones del término no restringido y el término restringido se integran por separado utilizando diferentes reglas de cuadratura numérica. La influencia de los puntos de interpolación, la regla de cuadratura y los puntos de restricción en la precisión de los algoritmos se analiza exhaustivamente. Las propiedades de los algoritmos propuestos se investigan utilizando ejemplos. Los resultados numéricos muestran que los algoritmos propuestos tienen un orden arbitrariamente alto, satisfacen las restricciones holonómicas con alta precisión y proporcionan un buen rendimiento para la integración a largo plazo.