  Two numerical schemes are proposed and investigated for the Yang--Mills
equations, which can be seen as a nonlinear generalisation of the Maxwell
equations set on Lie algebra-valued functions, with similarities to certain
formulations of General Relativity. Both schemes are built on the Discrete de
Rham (DDR) method, and inherit from its main features: an arbitrary order of
accuracy, and applicability to generic polyhedral meshes. They make use of the
complex property of the DDR, together with a Lagrange-multiplier approach, to
preserve, at the discrete level, a nonlinear constraint associated with the
Yang--Mills equations. We also show that the schemes satisfy a discrete energy
dissipation (the dissipation coming solely from the implicit time stepping).
Issues around the practical implementations of the schemes are discussed; in
particular, the assembly of the local contributions in a way that minimises the
price we pay in dealing with nonlinear terms, in conjunction with the
tensorisation coming from the Lie algebra. Numerical tests are provided using a
manufactured solution, and show that both schemes display a convergence in
$L^2$-norm of the potential and electrical fields in $\mathcal O(h^{k+1})$
(provided that the time step is of that order), where $k$ is the polynomial
degree chosen for the DDR complex. We also numerically demonstrate the
preservation of the constraint.
