Stability and error analysis of fully discrete original energy-dissipative and length-preserving scheme for the Landau-Lifshitz-Gilbert equation

The Landau-Lifshitz-Gilbert (LLG) equation, regarded as a gradient flow with manifold constraint, is the fundamental model describing magnetization dynamics in ferromagnetic materials. It is well known that the normalized tangent plane method is able to simultaneously achieve the non-convex manifold constraint and original energy dissipation. However, the associated computational cost of this numerical approach is exceedingly high. By contrast, the projection method is more straightforward to implement, while it often compromises the inherent energy dissipative property of the continuous model, and the error analysis turns out to be even more challenging. In this work, we first construct a linear and fully discrete finite difference numerical scheme, based on the projection method for the LLG equation, which is capable of simultaneously preserving the non-convex manifold constraint (|\mathbf{m}| = 1) and an unconditional original energy dissipation. In the error analysis, the classical theoretical technique becomes ineffective, due to the presence of the nonlinear Laplacian term, which in turn poses a significant challenge. To overcome this subtle difficulty, we carefully rewrite the numerical method in an equivalent weak form, in which a point-wise length preserving feature of the numerical solution plays an essential role. As a result of these estimates in the reformulated weak form, an optimal convergence rate could be theoretically established. In our knowledge, this numerical method is the first linear algorithm that preserves the following combined theoretical properties: (i) point-wise length preservation, (ii) unconditional original energy dissipation, (iii) a theoretical justification of convergence analysis and optimal rate error estimate.

💡 Research Summary

**
The paper addresses the long‑standing challenge of designing a numerical scheme for the Landau‑Lifshitz‑Gilbert (LLG) equation that simultaneously preserves the point‑wise unit‑length constraint (|m| = 1) and the original energy‑dissipation law of the continuous model. While the normalized tangent‑plane method achieves both properties, its computational cost is prohibitive because each time step requires solving a nonlinear system. The projection method is computationally cheap but typically destroys the original energy‑dissipation property, and rigorous error analysis for such schemes has been elusive.

The authors propose a linear, fully discrete finite‑difference scheme based on a semi‑implicit projection approach. The algorithm proceeds in two stages at each time step: (i) solve a linear system for an intermediate magnetization (\tilde{m}^{n+1}) using the discrete projection operator (P(m^n)) applied to the discrete Laplacian, and (ii) renormalize (\tilde{m}^{n+1}) point‑wise to obtain the next solution (m^{n+1} = \tilde{m}^{n+1}/|\tilde{m}^{n+1}|). Because the renormalization is performed after the linear solve, the scheme remains linear and straightforward to implement, yet the intermediate solution already satisfies the unit‑length condition, which is crucial for the subsequent analysis.

A key technical contribution is the reformulation of the scheme in an equivalent weak form. By applying discrete summation‑by‑parts and exploiting the point‑wise length preservation, the authors eliminate the troublesome nonlinear Laplacian term that would otherwise prevent a clean error estimate. Instead, they derive a “law‑of‑cosine” type bound for the error introduced by the renormalization step, showing that the constant in this bound scales like the time step (\Delta t). This scaling is essential to close the (H^1) error analysis.

The stability analysis proves an unconditional discrete energy dissipation: the discrete energy (E_h^{n}= \frac12|\nabla_h m^{n}|^2) satisfies (E_h^{n+1}\le E_h^{n}) for any choice of mesh size (h) and time step (\Delta t). No restrictive CFL‑type condition is required.

For convergence, the authors combine the weak‑form reformulation with discrete Sobolev interpolation inequalities and inverse estimates. They obtain optimal error bounds \