Abstract The paper focuses on the analysis of the Euler projection Galerkin finite element method (FEM) for the dynamics of magnetization in ferromagnetic materials, described by the Landau–Lifshitz equation with the point-wise constraint $|{\textbf{m}}|=1$. The method is based on a simple sphere projection that projects the numerical solution onto a unit sphere at each time step, and the method has been used in many areas in the past several decades. However, error analysis for the commonly used method has not been done since the classical energy approach cannot be applied directly. In this paper we present an optimal $\textbf{L}^2$ error analysis of the backward Euler sphere projection method by using quadratic or higher order finite elements under a time step condition $\tau =O(\epsilon _0 h)$ with some small $\epsilon _0>0$. The analysis is based on more precise estimates of the extra error caused by the sphere projection in both $\textbf{L}^2$ and $\textbf{H}^1$ norms, and the classical estimate of dual norm. Numerical experiment is provided to confirm our theoretical analysis.