  A variational integrator of arbitrarily high-order on the special orthogonal
group $SO(n)$ is constructed using the polar decomposition and the constrained
Galerkin method. It has the advantage of avoiding the second-order derivative
of the exponential map that arises in traditional Lie group variational
methods. In addition, a reduced Lie--Poisson integrator is constructed and the
resulting algorithms can naturally be implemented by fixed-point iteration. The
proposed methods are validated by numerical simulations on $SO(3)$ which
demonstrate that they are comparable to variational Runge--Kutta--Munthe-Kaas
methods in terms of computational efficiency. However, the methods we have
proposed preserve the Lie group structure much more accurately and and exhibit
better near energy preservation.
