The Landau–Lifshitz equation is a well-known mathematical model for describing the evolution of magnetization in ferromagnets. It is a partial differential equation that, due to its nonlinear structure, poses a number of challenges when one tries to construct a numerical solution. Assuming a homogeneous Neumann boundary condition and an initial condition of Dirichlet type, the paper begins with a review of already known numerical methods and their shortcomings. The most important problem in those methods is usually the high computational cost incurred when attempting to obtain a reasonably accurate approximate solution. As an alternative to these known approaches, An suggests a linearization of the problem, followed by a Galerkin approximation based on piecewise polynomials of degree 1 or 2 in space and a Crank–Nicolson scheme with respect to the time variable. Under a certain assumption on the smoothness of the exact solution, a lengthy and very technical analysis of the error term shows that a reasonably fast convergence can be expected. A sufficient condition for the smoothness assumption to be valid is given; however, this condition appears to be very restrictive, so it remains unclear whether the theory developed in the paper can be used in a wide class of realistic applications.