  I propose a method, based on a set of Langevin equations, for bringing
classical gauge theories to thermal equilibrium while respecting the set of
Gauss' constraints exactly. The algorithm is described in detail for the SU(2)
gauge theory with or without the Higgs doublet. As an example of application,
canonical average of the maximal Lyapunov exponent is computed for the SU(2)
Yang-Mills theory.
