Abstract. The magnetization of a ferromagnetic sample solves a non-convex variational problem, where its relaxation by convexifying the energy density resolves relevant macroscopic information. The numerical analysis of the relaxed model has to deal with a constrained convex but degenerated, nonlocal energy functional in mixed formulation for magnetic potential u and magnetization m. In [C. Carstensen and A. Prohl, Numer. Math. 90 (2001) 65-99], the conforming P 1 − (P 0) d -element in d = 2, 3 spatial dimensions is shown to lead to an ill-posed discrete problem in relaxed micromagnetism, and suboptimal convergence. This observation motivated a non-conforming finite element method which leads to a well-posed discrete problem, with solutions converging at optimal rate. In this work, we provide both an a priori and a posteriori error analysis for two stabilized conforming methods which account for inter-element jumps of the piecewise constant magnetization. Both methods converge at optimal rate; the new approach is applied to a macroscopic nonstationary ferromagnetic model [M. Kružík and A.