  This paper develops a weak Galerkin (WG) finite element method of arbitrary
order for the steady incompressible Magnetohydrodynamics equations. The WG
scheme uses piecewise polynomials of degrees $k(k\geq 1),k,k-1$, and $k-1$
respectively for the approximations of the velocity, the magnetic field, the
pressure, and the magnetic pseudo-pressure in the interior of elements, and
uses piecewise polynomials of degree $k$ for their numerical traces on the
interfaces of elements. The method is shown to yield globally divergence-free
approximations of the velocity and magnetic fields. We give existence and
uniqueness results for the discrete scheme and derive optimal a priori error
estimates. We also present a convergent linearized iterative algorithm.
Numerical experiments are provided to verify the obtained theoretical results.
