Higher-order linearly implicit full discretization of the Landau--Lifshitz--Gilbert equation

TL;DR

This paper develops and analyzes higher-order linearly implicit BDF time discretizations combined with advanced finite element methods for the Landau--Lifshitz--Gilbert equation, proving stability and error bounds.

Contribution

It introduces a novel combination of BDF methods up to order 5 with non-conforming finite element discretizations for LLG, providing rigorous stability and error analysis.

Findings

Stable and optimal error bounds for BDF methods of orders 3 to 5.

No time step restriction needed for BDF orders 1 and 2.

Discrete energy inequality established for certain methods.

Abstract

For the Landau--Lifshitz--Gilbert (LLG) equation of micromagnetics we study linearly implicit backward difference formula (BDF) time discretizations up to order $5$ combined with higher-order non-conforming finite element space discretizations, which are based on the weak formulation due to Alouges but use approximate tangent spaces that are defined by $L^2$-averaged instead of nodal orthogonality constraints. We prove stability and optimal-order error bounds in the situation of a sufficiently regular solution. For the BDF methods of orders $3$ to~$5$, this requires %a mild time step restriction $\tau \leqslant ch$ and that the damping parameter in the LLG equations be above a positive threshold; this condition is not needed for the A-stable methods of orders $1$ and $2$, for which furthermore a discrete energy inequality irrespective of solution regularity is proved.

Equations

$${\partial }_t\mathbf{m}=-\mathbf{m}\times {\mathbf{H}}_{\text{eff}}-\alpha \mathbf{m}\times (\mathbf{m}\times {\mathbf{H}}_{\text{eff}})$$

$$\alpha\,\partial_{t}\bm{m}+\bm{m}\times\partial_{t}\bm{m}=(1+\alpha^{2})\big{[}\bm{H}_{\text{eff}}-\big{(}\bm{m}\cdot\bm{H}_{\text{eff}}\big{)}\bm{m}\big{]}.$$

$$\mathbf{P}(\mathbf{m})=\mathbf{I}-\mathbf{m}{\mathbf{m}}^T$$

$$\bm{H}_{\text{eff}}=\frac{1}{1+\alpha^{2}}\bigl{(}\varDelta\bm{m}+\bm{H}\bigr{)},$$

$$\alpha {\partial }_t\mathbf{m}+\mathbf{m}\times {\partial }_t\mathbf{m}=\mathbf{P}(\mathbf{m})(\Delta \mathbf{m}+\mathbf{H}).$$

$$T(\bm{m}):=\big{\{}\bm{\varphi}\in L^{2}(\varOmega)^{3}\,:\,\bm{m}\cdot\bm{\varphi}=0\ \text{ a.e.}\big{\}}=\big{\{}\bm{\varphi}\in L^{2}(\varOmega)^{3}\,:\,{\mathbf{P}}(\bm{m})\bm{\varphi}=\bm{\varphi}\}$$

$$\displaystyle\alpha\bigl{(}\partial_{t}\bm{m},\bm{\varphi}\bigr{)}+\bigl{(}\bm{m}\times\partial_{t}\bm{m},\bm{\varphi}\bigr{)}+\bigl{(}\nabla\bm{m},\nabla\bm{\varphi}\bigr{)}=\bigl{(}\bm{H},\bm{\varphi}\bigr{)},$$

$$\delta(\zeta)=\sum_{\ell=1}^{k}\frac{1}{\ell}(1-\zeta)^{\ell}=\sum_{j=0}^{k}\delta_{j}\zeta^{j},\quad\gamma(\zeta)=\frac{1}{\zeta}\big{[}1-(1-\zeta)^{k}\big{]}=\sum_{i=0}^{k-1}\gamma_{i}\zeta^{i}.$$

$$\widehat{\bm{m}}^{n}:={\displaystyle{\sum_{j=0}^{k-1}\gamma_{j}\bm{m}^{n-j-1}}}\Big{/}{\Big{|}\displaystyle{\sum_{j=0}^{k-1}\gamma_{j}\bm{m}^{n-j-1}}\Big{|}}.$$

$$\dot{\bm{m}}^{n}=\frac{1}{\tau}\sum_{j=0}^{k}\delta_{j}\bm{m}^{n-j},\ \text{i.e., }\ {\bm{m}}^{n}=\Bigl{(}-\sum_{j=1}^{k}\delta_{j}\bm{m}^{n-j}+\tau\dot{\bm{m}}^{n}\Bigr{)}/\delta_{0}.$$

$$\displaystyle\alpha\bigl{(}\dot{\bm{m}}^{n},\bm{\varphi}\bigr{)}+\bigl{(}\widehat{\bm{m}}^{n}\times\dot{\bm{m}}^{n},\bm{\varphi}\bigr{)}+\bigl{(}\nabla{\bm{m}}^{n},\nabla\bm{\varphi}\bigr{)}=\bigl{(}\bm{H}(t_{n}),\bm{\varphi}\bigr{)}$$

$${\dot{\mathbf{m}}}^n\in T({\hat{\mathbf{m}}}^n),\text{ i.e., }{\hat{\mathbf{m}}}^n\cdot {\dot{\mathbf{m}}}^n=0.$$

$$\alpha\bigl{(}\dot{\bm{m}}^{n},\bm{\varphi}\bigr{)}+\bigl{(}\widehat{\bm{m}}^{n}\times\dot{\bm{m}}^{n},\bm{\varphi}\bigr{)}+\frac{\tau}{\delta_{0}}\bigl{(}\nabla\dot{\bm{m}}^{n},\nabla\bm{\varphi}\bigr{)}=\bigl{(}\bm{f}^{n},\bm{\varphi}\bigr{)}\quad\ \forall\bm{\varphi}\in T(\widehat{\bm{m}}^{n})\cap H^{1}(\varOmega)^{3},$$

$$T_h(\mathbf{m})$$

$$=\{{\mathbf{\phi }}_h\in V_h^3:{\Pi }_h(\mathbf{m}\cdot {\mathbf{\phi }}_h)=0\}.$$

$$\alpha\bigl{(}\partial_{t}\bm{m}_{h},\bm{\varphi}_{h}\bigr{)}+\bigl{(}\bm{m}_{h}\times\partial_{t}\bm{m}_{h},\bm{\varphi}_{h}\bigr{)}+\bigl{(}\nabla\bm{m}_{h},\nabla\bm{\varphi}_{h}\bigr{)}=\bigl{(}\bm{H},\bm{\varphi}_{h}\bigr{)}\quad\forall\bm{\varphi}_{h}\in T_{h}(\bm{m}_{h}),$$

$$\alpha\bigl{(}\dot{\bm{m}}_{h}^{n},\bm{\varphi}_{h}\bigr{)}+\bigl{(}\widehat{\bm{m}}_{h}^{n}\times\dot{\bm{m}}_{h}^{n},\bm{\varphi}_{h}\bigr{)}+\bigl{(}\nabla{\bm{m}}_{h}^{n},\nabla\bm{\varphi}_{h}\bigr{)}=\bigl{(}\bm{H}^{n},\bm{\varphi}_{h}\bigr{)}\quad\forall\bm{\varphi}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n}),$$

$$\dot{\bm{m}}^{n}_{h}=\frac{1}{\tau}\sum_{j=0}^{k}\delta_{j}\bm{m}^{n-j}_{h},\qquad\widehat{\bm{m}}^{n}_{h}={\displaystyle{\sum_{j=0}^{k-1}\gamma_{j}\bm{m}^{n-j-1}_{h}}}\Big{/}{\Big{|}\displaystyle{\sum_{j=0}^{k-1}\gamma_{j}\bm{m}^{n-j-1}_{h}}\Big{|}}.$$

$$K^n=\alpha \mathbf{I}\otimes M+\frac{\tau }{{\delta }_0}\mathbf{I}\otimes A+S^n.$$

$$K^n{\dot{m}}^n+(C^n{)}^T{\lambda }^n$$

$$C^n{\dot{m}}^n$$

$$(\mathbf{m}\cdot {\mathbf{\phi }}_h,{\psi }_h)=0\text{for all }{\psi }_h\in V_h\text{ with }\mathrm{supp}({\psi }_h)\subseteq \overline{\omega }.$$

$$(B^n{)}^T{K}^n{B}^n{x}^n=(B^n{)}^T{f}^n.$$

$$T_h^{\mathrm{pw}}(\mathbf{m})$$

$$=\{{\mathbf{\phi }}_h\in V_h^3:I_h(\mathbf{m}\cdot {\mathbf{\phi }}_h)=0\}={\mathbf{I}}_h\mathbf{P}(\mathbf{m})V_h^3,$$

$$\displaystyle\alpha\bigl{(}\dot{\bm{m}}_{h}^{n},\bm{\varphi}_{h}\bigr{)}$$

$$\displaystyle{}-\bigl{(}\nabla\widehat{\bm{m}}_{h}^{n},\nabla(\mathbf{I}-\mathbf{P}(\widehat{\bm{m}}_{h}^{n}))\bm{\varphi}_{h}\bigr{)}=\bigl{(}\mathbf{P}(\widehat{\bm{m}}_{h}^{n})\bm{H}(t_{n}),\bm{\varphi}_{h}\bigr{)},$$

$${\tau }^k⩽\bar{c}{h}^{1/2}$$

$$\Vert {\mathbf{m}}_h^n-\mathbf{m}(t_n){\Vert }_{H^1(\Omega {)}^3}⩽C({\tau }^k+h^r)\text{ for }{t}_n=n\tau ⩽\bar{t},$$

$$\mathbf{m}\in C^{k+1}([0,\bar{t}],L^{\infty }(\Omega {)}^3)\cap C^1([0,\bar{t}],W^{r+1,\infty }(\Omega {)}^3),$$

$$\Delta \mathbf{m}+\mathbf{H}\in C([0,\bar{t}],W^{r+1,\infty }(\Omega {)}^3).$$

$$\Vert 1-|{\mathbf{m}}_h^n|{\Vert }_{L^2(\Omega )}⩽C({\tau }^k+h^r)\text{ for }{t}_n=n\tau ⩽\bar{t},$$

$$\Bigl{\|}1-\bigl{|}\sum_{j=0}^{k-1}\gamma_{j}\bm{m}^{n-j-1}_{h}\bigr{|}\Bigr{\|}_{L^{\infty}(\varOmega)}\leqslant Ch^{-1/2}(\tau^{k}+h^{r})\leqslant\tfrac{1}{2}\quad\text{ for }\ t_{n}=n\tau\leqslant\bar{t},$$

Full text

1. Introduction

1.1. Scope

In this paper we study the convergence of higher-order time and space discretizations of the Landau–Lifshitz–Gilbert (LLG) equation, which is the basic model for phenomena in micromagnetism, such as in recording media [26, 36].

The main novelty of the paper lies in the construction and analysis of what is apparently the first numerical method for the LLG equation that is second-order convergent in both space and time to sufficiently regular solutions and that satisfies, as an important robustness property irrespective of regularity, a discrete energy inequality analogous to that of the continuous problem.

We study discretization in time by linearly implicit backward difference formulae (BDF) up to order $5$ and discretization in space by finite elements of arbitrary polynomial degree. For the BDF methods up to order $2$ we prove optimal-order error bounds in the situation of a sufficiently regular solution and a discrete energy inequality irrespective of solution regularity under very weak regularity assumptions on the data. For the BDF methods of orders $3$ to $5$, we prove optimal-order error bounds in the situation of a sufficiently regular solution under the additional condition that the damping parameter in the LLG equation be above a method-dependent positive threshold. However, no discrete energy inequality irrespective of solution regularity is obtained for the BDF methods of orders $3$ to $5$.

The discretization in space is done by a higher-order non-conforming finite element method based on the approach of Alouges [4, 5], which uses a projection to an approximate tangent space to the normality constraint. Contrary to the pointwise orthogonality constraints in the nodes, which define the approximate tangent space in those papers and yield only first-order convergence also for finite elements with higher-degree polynomials, we here enforce orthogonality averaged over the finite element basis functions. With these modified approximate tangent spaces we prove $H^{1}$-convergence of optimal order in space and time under the assumption of a sufficiently regular solution.

Key issues in the error analysis are the properties of the orthogonal projection onto the approximate tangent space, the higher-order consistency error analysis, and the proof of stable error propagation, which is based on non-standard energy estimates and uses both $L^{2}$ and maximum norm finite element analysis.

1.2. The Landau–Lifshitz–Gilbert equation

The standard phenomenological model for micromagnetism is provided by the Landau–Lifshitz (LL) equation

[TABLE]

where the unknown magnetization field $\bm{m}=\bm{m}(x,t)$ takes values on the unit sphere $\mathbb{S}^{2}$, $\alpha>0$ is a dimensionless damping parameter, and the effective magnetic field $\bm{H}_{\text{eff}}$ depends on the unknown $\bm{m}$. The Landau–Lifshitz equation (1.1) can be equivalently written in the Landau–Lifshitz–Gilbert form

[TABLE]

Indeed, in view of the vector identity $\bm{a}\times(\bm{b}\times\bm{c})=(\bm{a}\cdot\bm{c})\bm{b}-(\bm{a}\cdot\bm{b})\bm{c},$ for $\bm{a},\bm{b},\bm{c}\in{\mathbb{R}}^{3},$ we have -\bm{m}\times\big{(}\bm{m}\times\bm{H}_{\text{eff}}\big{)}=\bm{H}_{\text{eff}}-\big{(}\bm{m}\cdot\bm{H}_{\text{eff}}\big{)}\bm{m}, and taking the vector product of (1.1) with $\bm{m}$ and adding $\alpha$ times (1.1) then yields (1.2).

Since $\bm{m}\times\bm{a}$ is orthogonal to $\bm{m},$ for any $\bm{a}\in{\mathbb{R}}^{3},$ it is obvious from (1.1) that $\partial_{t}\bm{m}$ is orthogonal to $\bm{m}$: $\bm{m}\cdot\partial_{t}\bm{m}=0;$ we infer that the Euclidean norm satisfies $|\bm{m}(x,t)|=1$ for all $x$ and for all $t$, provided this is satisfied for the initial data.

The term in square brackets on the right-hand side in (1.2) can be rewritten as ${\mathbf{P}}(\bm{m})\bm{H}_{\text{eff}}$, where (with ${\mathbf{I}}$ the $3\times 3$ unit matrix)

[TABLE]

is the orthogonal projection onto the tangent plane to the unit sphere $\mathbb{S}^{2}$ at $\bm{m}$.

In this paper we consider the situation

[TABLE]

where $\bm{H}=\bm{H}(x,t)$ is a given external magnetic field. The factor $1/(1+\alpha^{2})$ is chosen for convenience of presentation, but is inessential for the theory; it can be replaced by any positive constant factor.

With this choice of $\bm{H}_{\text{eff}}$, we arrive at the Landau–Lifshitz–Gilbert (LLG) equation in the form

[TABLE]

We consider this equation as an initial-boundary value problem on a bounded domain $\varOmega\subset{\mathbb{R}}^{3}$ and a time interval $0\leqslant t\leqslant\bar{t}$, with homogeneous Neumann boundary conditions and initial data $\bm{m}_{0}$ taking values on the unit sphere, i.e., the Euclidean norm $|\bm{m}_{0}(x)|$ equals $1$ for all $x\in\varOmega$.

We consider the following weak formulation, first proposed by Alouges [4, 5]: Find the solution $\bm{m}:\varOmega\times[0,\bar{t}\,]\to\mathbb{S}^{2}$ with $\bm{m}(\cdot,0)=\bm{m}_{0}$ by determining, at $\bm{m}(t)\in H^{1}(\varOmega)^{3}$, the time derivative $\partial_{t}\bm{m}$ (omitting here and in the following the argument $t$) as that function in the tangent space

[TABLE]

that satisfies, for all $\bm{\varphi}\in T(\bm{m})\cap H^{1}(\varOmega)^{3}$,

[TABLE]

where the brackets $(\cdot,\cdot)$ denote the $L^{2}$ inner product over the domain $\varOmega$. The numerical methods studied in this paper are based on this weak formulation.

1.3. Previous work

There is a rich literature on numerical methods for Landau–Lifshitz(–Gilbert) equations; for the numerical literature up to $2007$ see the review by Cimrák [17].

Alouges & Jaisson [4, 5] propose linear finite element discretizations in space and linearly implicit backward Euler in time for the LLG equation in the weak formulation (1.5) and prove convergence without rates towards nonsmooth weak solutions, using a discrete energy inequality and compactness arguments. Convergence of this type was previously shown by Bartels & Prohl [11] for fully implicit methods that are based on a different formulation of the Landau–Lifshitz equation (1.1). In [6], convergence without rates towards weak solutions is shown for a method that is (formally) of “almost” order $2$ in time, based on the midpoint rule, for the LLG equation with an effective magnetic field of a more general type than (1.3).

In a complementary line of research, convergence with rates has been studied under sufficiently strong regularity assumptions, which can, however, not be guaranteed over a given time interval, since solutions of the LLG equation may develop singularities. A first-order error bound for a linearly implicit time discretization of the Landau–Lifshitz equation (1.1) was proved by Cimrák [16]. Optimal-order error bounds for linearly implicit time discretizations based on the backward Euler and Crank–Nicolson methods combined with finite element full discretizations for a different version of the Landau–Lifshitz equation (1.1) were obtained under sufficient regularity assumptions by Gao [23] and An [7], respectively. In contrast to [4, 5, 6, 11], these methods do not satisfy an energy inequality irrespective of the solution regularity.

Numerical discretizations for the coupled system of the LLG equation (1.5) with the eddy current approximation of the Maxwell equations are studied by Feischl & Tran [21], with first-order error bounds in space and time under sufficient regularity assumptions. This also yields the first result of first-order convergence of the method of Alouges & Jaisson [4, 5].

There are several methods for the LLG equations that are of formal order $2$ in time (though only of order $1$ in space), e.g., [35, 31, 19], but none of them comes with an error analysis. Fully implicit BDF time discretizations for LLG equations have been used successfully in the computational physics literature [37], though without giving any error analysis.

To the authors’ knowledge, the second-order linearly implicit method proposed and studied here is thus the first numerical method for the LLG (or LL) equation that has rigorous a priori error estimates of order $2$ in both space and time under high regularity assumptions and that satisfies a discrete energy inequality irrespective of regularity.

We conclude this brief survey of the literature with a remark: The existing convergence results either give convergence of a subsequence without rates to a weak solution (without imposing strong regularity assumptions), or they show convergence with rates towards sufficiently regular solutions (as we do here). Both approaches yield insight into the numerical methods and have their merits, and they complement each other. Clearly, neither approach is fully satisfactory, because convergence without rates of some subsequence is nothing to observe in actual computations, and on the other hand high regularity is at best provable for close to constant initial conditions [22] or over short time intervals. We regard the situation as analogous to the development of numerical methods and their analysis in other fields such as nonlinear hyperbolic conservation laws: second-order methods are highly popular in that field, even though they can only be shown to converge with very low order ($1/2$ or less or only without rates) for available regularity properties; see, e.g., [32, Chapter 3]. Nevertheless, second-order methods are favored over first-order methods in many applications, especially if they enjoy some qualitative properties that give them robustness in non-regular situations. A similar situation occurs with the LLG equation, where the most important qualitative property appears to be the energy inequality.

1.4. Outline

In Section 2 we describe the numerical methods studied in this paper. They use time discretization by linearly implicit BDF methods of orders up to $5$ and space discretization by finite elements of arbitrary polynomial degree in a numerical scheme that is based on the weak formulation (1.5), with an approximate tangent space that enforces the orthogonality constraint approximately in an $L^{2}$-projected sense.

In Section 3 we state our main results:

\bullet\ For the full discretization of (1.5) by linearly implicit BDF methods of orders $1$ and $2$ and finite element methods of arbitrary polynomial degree we give optimal-order error bounds in the $H^{1}$ norm, under very mild mesh conditions, in the case of sufficiently regular solutions (Theorem 3.1). For these methods we also show a discrete energy inequality that requires only very weak regularity assumptions on the data (Proposition 3.1). This discrete energy inequality is of the same type as the one used in [5, 11] for proving convergence without rates to a weak solution.

\bullet\ For the linearly implicit BDF methods of orders $3$ to $5$ and finite element methods with polynomial degree at least $2$, we have optimal-order error bounds in the $H^{1}$ norm only if the damping parameter $\alpha$ is larger than some positive threshold, which depends on the order of the BDF method (Theorem 3.2). Moreover, a stronger (but still mild) CFL condition $\tau\leqslant ch$ is required. A discrete energy inequality under very weak regularity conditions is not available for the BDF methods of orders $3$ to $5$, in contrast to the A-stable BDF methods of orders $1$ and $2$.

In Section 4 we prove a perturbation result for the continuous problem by energy techniques, as a preparation for the proofs of our error bounds for the discretization.

In Section 5 we study properties of the $L^{2}$-orthogonal projection onto the discrete tangent space, which are needed to ensure consistency of the full order and stability of the space discretization with the higher-order discrete tangent space.

In Section 6 we study consistency properties of the methods and present the error equation.

In Sections 7 and 8 we prove Theorems 3.1 and 3.2, respectively. The higher-order convergence proofs are separated into consistency (Section 6) and stability estimates. The stability proofs use the technique of energy estimates, in an unusual version where the error equation is tested with a projection of the discrete time derivative of the error onto the discrete tangent space. These proofs are different for the A-stable BDF methods of orders $1$ and $2$ and for the BDF methods of orders $3$ to $5$. For the control of nonlinearities, the stability proofs also require pointwise error bounds, which are obtained with the help of finite element inverse inequalities from the $H^{1}$ error bounds of previous time steps.

In Section 9 we illustrate our results by numerical experiments.

In an Appendix we collect basic results on energy techniques for BDF methods that are needed for our stability proofs.

2. Discretization of the LLG equation

We now describe the time and space discretization that is proposed and studied in this paper.

2.1. Time discretization by linearly implicit BDF methods

We shall discretize the LLG equation (1.5) in time by the linearly implicit $k$-step BDF methods, $1\leqslant k\leqslant 5$, described by the polynomials $\delta$ and $\gamma,$

[TABLE]

We let $t_{n}=n\tau,\ n=0,\dotsc,N,$ be a uniform partition of the interval $[0,\bar{t}\,]$ with time step $\tau=\bar{t}/N.$ For the $k$-step method we require $k$ starting values $\bm{m}^{i}$ for $i=0,\dotsc,k-1$. For $n\geqslant k$, we determine the approximation $\bm{m}^{n}$ to $\bm{m}(t_{n})$ as follows. We first extrapolate the known values $\bm{m}^{n-k},\dotsc,\bm{m}^{n-1}$ to a preliminary normalized approximation $\widehat{\bm{m}}^{n}$ at $t_{n}$,

[TABLE]

To avoid potentially undefined quantities, we define $\widehat{\bm{m}}^{n}$ to be an arbitrary fixed unit vector if the denominator in the above formula is zero.

The derivative approximation $\dot{\bm{m}}^{n}$ and the solution approximation ${\bm{m}}^{n}$ are related by the backward difference formula

[TABLE]

We determine ${\bm{m}}^{n}$ by requiring that for all $\bm{\varphi}\in T(\widehat{\bm{m}}^{n})\cap H^{1}(\varOmega)^{3}$,

[TABLE]

Here we note that on inserting the formula in (2.2) for ${\bm{m}}^{n}$ in the third term of (2.3), we obtain a linear constrained elliptic equation for $\dot{\bm{m}}^{n}\in T(\widehat{\bm{m}}^{n})\cap H^{1}(\varOmega)^{3}$ of the form

[TABLE]

where $\bm{f}^{n}$ consists of known terms. The bilinear form on the left-hand side is $H^{1}(\varOmega)^{3}$-coercive on $T(\widehat{\bm{m}}^{n})\cap H^{1}(\varOmega)^{3}$, and hence the above linear equation has a unique solution $\dot{\bm{m}}^{n}\in T(\widehat{\bm{m}}^{n})\cap H^{1}(\varOmega)^{3}$ by the Lax–Milgram lemma. Once this elliptic equation is solved for $\dot{\bm{m}}^{n}$, we obtain the approximation ${\bm{m}}^{n}\in H^{1}(\varOmega)^{3}$ to ${\bm{m}}(t_{n})$ from the second formula in (2.2).

2.2. Full discretization by BDF and higher-order finite elements

For a family of regular and quasi-uniform finite element triangulations of $\varOmega$ with maximum meshwidth $h>0$ we form the Lagrange finite element spaces $V_{h}\subset H^{1}(\varOmega)$ with piecewise polynomials of degree $r\geqslant 1$. We denote the $L^{2}$-orthogonal projections onto the finite element space by $\varPi_{h}\colon L^{2}(\varOmega)\to V_{h}$ and ${\mathbf{\Pi}}_{h}={\mathbf{I}}\otimes\varPi_{h}\colon L^{2}(\varOmega)^{3}\to V_{h}^{3}$. With a function $\bm{m}\in H^{1}(\varOmega)^{3}$ that vanishes nowhere on $\varOmega$, we associate the discrete tangent space

[TABLE]

This space is different from the discrete tangent space used in [4, 5], where the orthogonality constraint $\bm{m}\cdot\bm{\varphi}_{h}=0$ is required to hold pointwise at the finite element nodes. Here, the constraint is enforced weakly on the finite element space, as is done in various saddle point problems for partial differential equations, for example for the divergence-free constraint in the Stokes problem [14, 25]. In contrast to that example, here the bilinear form associated with the linear constraint, i.e., $b(\bm{m};\bm{\varphi}_{h},v_{h})=(\bm{m}\cdot\bm{\varphi}_{h},v_{h})$, depends on the state $\bm{m}$. This dependence substantially affects both the implementation and the error analysis.

Following the general approach of [4, 5] with this modified discrete tangent space, we discretize (1.5) in space by determining the time derivative $\partial_{t}\bm{m}_{h}(t)\in T_{h}(\bm{m}_{h}(t))$ such that (omitting the argument $t$)

[TABLE]

where the brackets $(\cdot,\cdot)$ denote again the $L^{2}$ inner product over the domain $\varOmega$.

The full discretization with the linearly implicit BDF method is then readily obtained from (2.3): determine $\dot{\bm{m}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}^{n}_{h})$ such that

[TABLE]

where $\widehat{\bm{m}}_{h}^{n}$ and $\dot{\bm{m}}_{h}^{n}$ are related to $\bm{m}_{h}^{n-j}$ for $j=0,\dotsc,k$ in the same way as in (2.1) and (2.2) above with ${\bm{m}}_{h}^{n-j}$ in place of $\bm{m}^{n-j}$, viz.,

[TABLE]

To avoid potentially undefined quantities, we define $\widehat{\bm{m}}^{n}_{h}$ to be an arbitrary fixed unit vector if the denominator in the above formula is zero. (We will, however, show that this does not occur in the situation of sufficient regularity.)

To implement the discrete tangent space $T_{h}(\widehat{\bm{m}}^{n}_{h})$, there are at least two options: using the constraints $\varPi_{h}(\bm{m}\cdot\bm{\varphi}_{h})=0$ or constructing a local basis of $T_{h}(\bm{m})$.

(a) Constraints: Let $\phi_{i}$ for $i=1,\dotsc,N:=\text{dim}V_{h}$ denote the nodal basis of $V_{h}$ and denote the basis functions of $V_{h}^{3}$ by $\mathbf{\phi}_{\bm{i}}=\bm{e}_{k}\otimes\phi_{i}$ for $\bm{i}=(i,k)$, where $\bm{e}_{k}$ for $k=1,2,3$ are the standard unit vectors of ${\mathbb{R}}^{3}$. We denote by $M$ and $A$ the usual mass and stiffness matrices, respectively, with entries $m_{ij}=(\phi_{i},\phi_{j})_{L^{2}(\varOmega)}$ and $a_{ij}=(\nabla\phi_{i},\nabla\phi_{j})_{L^{2}(\varOmega)^{3}}$. We further introduce the sparse skew-symmetric matrix $S^{n}=(s_{\bm{i},\bm{j}}^{n})\in{\mathbb{R}}^{3N\times 3N}$ with entries $s_{\bm{i},\bm{j}}^{n}=(\widehat{\bm{m}}^{n}_{h}\times\phi_{\bm{i}},\phi_{\bm{j}})_{L^{2}(\varOmega)^{3}}$ and the sparse constraint matrix $C^{n}=(c_{\bm{i},j}^{n})\in{\mathbb{R}}^{3N\times N}$ by $c_{\bm{i},j}^{n}=(\widehat{\bm{m}}^{n}_{h}\cdot\phi_{\bm{i}},\phi_{j})_{L^{2}(\varOmega)}$. Finally, we denote the matrix of the unconstrained time-discrete problem as

[TABLE]

Let $\dot{m}^{n}\in{\mathbb{R}}^{3N}$ denote the nodal vector of $\dot{\bm{m}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$. In this setting, (2.6) yields a system of linear equations of saddle point type

[TABLE]

where $\lambda^{n}\in{\mathbb{R}}^{N}$ is the unknown vector of Lagrange multipliers and $f^{n}\in{\mathbb{R}}^{3N}$ is a known right-hand side.

(b) Local basis: It is possible to compute a local basis of $T_{h}(\bm{m})$ by solving small local problems. To see that, let $\omega\subset\varOmega$ denote a collection of elements of the mesh and let $\overline{\omega}\supset\omega$ denote the same set plus the layer of elements touching $\omega$ (the patch of $\omega$). A sufficient (and necessary) condition for $\bm{\varphi}_{h}\in V_{h}^{3}$ with $\operatorname{supp}(\bm{\varphi}_{h})\subseteq\omega$ to belong to $T_{h}(\bm{m})$ is

[TABLE]

If we denote by $\#\omega$ the number of generalized hat functions of $V_{h}$ supported in $\omega$, the space of functions in $V_{h}^{3}$ with support in $\omega$ is $3\#\omega$-dimensional. On the other hand, the space of test functions in (2.8) is $\#\overline{\omega}$-dimensional. We may choose $\omega$ sufficiently large (depending only on shape regularity) such that $3\#\omega>\#\overline{\omega}$ and hence (2.8) has at least one solution which is then a local basis function of $T_{h}(\bm{m})$. Choosing different $\omega$ to cover $\varOmega$ yields a full basis of $T_{h}(\bm{m})$.

Let us denote the so obtained basis of $T_{h}(\widehat{\bm{m}}_{h}^{n})$ by $(\psi_{\ell}^{n})$, given via $\psi_{\ell}^{n}=\sum_{\bm{i}}\phi_{\bm{i}}b_{\bm{i}\ell}^{n}$, and the sparse basis matrix by $B^{n}=(b_{\bm{i}\ell}^{n})$. Then, the nodal vector $\dot{m}^{n}=B^{n}x^{n}$ is obtained by solving the linear system

[TABLE]

An advantage of this approach is that the dimension is roughly halved compared to the formulation with constraints. However, the efficiency of one approach versus the other depends heavily on the numerical linear algebra used. Such comparisons are outside the scope of this paper.

Remark 2.1**.**

The algorithm described above does not enforce the norm constraint $|\bm{m}|=1$ at the nodes. The user might add a normalization step in the definition of $\bm{m}^{n}$ in (2.2). However, here we do not consider this normalized variant of the method, whose convergence properties are not obvious to derive.

Remark 2.2**.**

Differently to [4], we do not use the pointwise discrete tangent space

[TABLE]

where $I_{h}:C(\bar{\varOmega})\to V_{h}$ denotes finite element interpolation and $\mathbf{I}_{h}={\mathbf{I}}\otimes I_{h}:C(\bar{\varOmega})^{3}\to V_{h}^{3}$. It is already reported in [4, Section 4] that an improvement of the order with higher-degree finite elements could not be observed in numerical experiments when using the pointwise tangent spaces in the discretization (2.5). Our analysis shows a lack of consistency of optimal order in the discretization with $T_{h}^{\rm pw}(\bm{m})$, which originates from the fact that $\mathbf{I}_{h}{\mathbf{P}}(\bm{m})$ is not self-adjoint. The order reduction can, however, be cured by adding a correction term: in the $n$th time step, determine $\dot{\bm{m}}_{h}^{n}\in T_{h}^{\rm pw}(\widehat{\bm{m}}_{h}^{n})$ such that for all $\bm{\varphi}_{h}\in T_{h}^{\rm pw}(\widehat{\bm{m}}_{h}^{n})$,

[TABLE]

with notation $\widehat{\bm{m}}_{h}^{n}$ and $\dot{\bm{m}}_{h}^{n}$ as in (2.7). With the techniques of the present paper, it can be shown that like (2.6), also this discretization converges with optimal order in the $H^{1}$ norm under sufficient regularity conditions. Since this paper is already rather long, we do not include the proof of this result. In contrast to (2.6) for the first- and second-order BDF methods, the method (2.9) does not admit an $h$- and $\tau$-independent bound of the energy that is irrespective of the smoothness of the solution.

3. Main results

3.1. Error bound and energy inequality for BDF of orders 1 and 2

For the full discretization with first- and second-order BDF methods and finite elements of arbitrary polynomial degree $r\geqslant 1$ we will prove the following optimal-order error bound in Sections 5 to 7.

Theorem 3.1** (Error bound for orders $k=1,2$).**

Consider the full discretization (2.6) of the LLG equation (1.4) by the linearly implicit $k$-step BDF time discretization for $k\leqslant 2$ and finite elements of polynomial degree $r\geqslant 1$ from a family of regular and quasi-uniform triangulations of $\varOmega$. Suppose that the solution $\bm{m}$ of the LLG equation is sufficiently regular. Then, there exist $\bar{\tau}>0$ and $\bar{h}>0$ such that for numerical solutions obtained with step sizes $\tau\leqslant\bar{\tau}$ and meshwidths $h\leqslant\bar{h}$, which are restricted by the very mild CFL-type condition

[TABLE]

with a sufficiently small constant $\bar{c}$ $($independent of $h$ and \tau$$), the errors are bounded by

[TABLE]

where $C$ is independent of $h,\tau$ and $n$ $($but depends on $\alpha$ and exponentially on $\bar{t}$ $)$, provided that the errors of the starting values also satisfy such a bound.

The precise regularity requirements are as follows:

[TABLE]

Remark 3.1** (Discrepancy from normality).**

Since $\bm{m}(x,t_{n})$ are unit vectors, an immediate consequence of the error estimate (3.1) is that

[TABLE]

with a constant $C$ independent of $n,\tau$ and $h$. The proof of Theorem 3.1 also shows that the denominator in the definition of the normalized extrapolated value ${\widehat{\bm{m}}}^{n}_{h}$ satisfies

[TABLE]

which in particular ensures that ${\widehat{\bm{m}}}^{n}_{h}$ is unambiguously defined.

Testing with $\bm{\varphi}=\partial_{t}\bm{m}\in T(\bm{m})$ in (1.5), we obtain (only formally, if $\partial_{t}\bm{m}$ is not in $H^{1}(\varOmega)^{3}$)

[TABLE]

which, by integration in time and the Cauchy–Schwarz and Young inequalities, implies the energy inequality

[TABLE]

Similarly, we test with $\bm{\varphi}_{h}={\dot{\bm{m}}}^{n}_{h}\in T_{h}({\widehat{\bm{m}}}^{n}_{h})$ in (2.6). Then we can prove the following discrete energy inequality, which holds under very weak regularity assumptions on the data.

Proposition 3.1** (Energy inequality for orders $k=1,2$).**

Consider the full discretization (2.6) of the LLG equation (1.4) by the linearly implicit $k$-step BDF time discretization for $k\leqslant 2$ and finite elements of polynomial degree $r\geqslant 1$. Then, the numerical solution satisfies the following discrete energy inequality:* for $n\geqslant k$ with $n\tau\leqslant\bar{t}$,*

[TABLE]

where $\gamma_{1}^{\pm}=1$ and $\gamma_{2}^{\pm}=(3\pm 2\sqrt{2})/4$.

This energy inequality is an important robustness indicator of the numerical method. In [5, 11], such energy inequalitys are used to prove convergence without rates (for a subsequence $\tau_{n}\to 0$ and $h_{n}\to 0$) to a weak solution of the LLG equation for the numerical schemes considered there (which have $\gamma^{\pm}=1$, but this is inessential in the proofs).

As the proof of Proposition 3.1 is short, we give it here.

Proof.

The proof relies on the A-stability of the first- and second-order BDF methods via Dahlquist’s G-stability theory as expressed in Lemma 10.1 of the Appendix, used with $\delta(\zeta)=\sum_{\ell=1}^{k}(1-\zeta)^{\ell}/\ell$ and $\mu(\zeta)=1$. The positive definite symmetric matrices $G=(g_{ij})_{i,j=1}^{k}$ are known to be $G=1$ for $k=1$ and (see [27, p. 309])

[TABLE]

which has the eigenvalues $\gamma^{\pm}=(3\pm 2\sqrt{2})/4$.

We test with $\bm{\varphi}_{h}={\dot{\bm{m}}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$ in (2.6) and note \bigl{(}\widehat{\bm{m}}_{h}^{n}\times\dot{\bm{m}}_{h}^{n},{\dot{\bm{m}}}_{h}^{n}\bigr{)}=0, so that

[TABLE]

The right-hand side is bounded by

[TABLE]

Recalling the definition of ${\dot{\bm{m}}}_{h}^{n}$, we have by Lemma 10.1

[TABLE]

We fix $\bar{n}$ with $k\leqslant\bar{n}\leqslant\bar{t}/\tau$ and sum from $n=k$ to $\bar{n}$ to obtain

[TABLE]

Noting that

[TABLE]

we obtain the stated result. ∎

3.2. Error bound for BDF of orders $3$ to $5$

For the BDF methods of orders $3$ to $5$ we prove the following result in Section 8. Here we require a stronger, but still moderate stepsize restriction in terms of the meshwidth. More importantly, we must impose a positive lower bound on the damping parameter $\alpha$ of (1.1).

Theorem 3.2** (Error bound for orders $k=3,4,5$).**

Consider the full discretization (2.6) of the LLG equation (1.4) by the linearly implicit $k$-step BDF time discretization for $3\leqslant k\leqslant 5$ and finite elements of polynomial degree $r\geqslant 2$ from a family of regular and quasi-uniform triangulations of $\varOmega$. Suppose that the solution $\bm{m}$ of the LLG equation has the regularity (3.2), and that the damping parameter $\alpha$ satisfies

[TABLE]

Then, for an arbitrary constant $\bar{C}>0$, there exist $\bar{\tau}>0$ and $\bar{h}>0$ such that for numerical solutions obtained with step sizes $\tau\leqslant\bar{\tau}$ and meshwidths $h\leqslant\bar{h}$ that are restricted by

[TABLE]

*the errors are bounded by *

[TABLE]

where $C$ is independent of $h,\tau$ and $n$ $($but depends on $\alpha$ and exponentially on \bar{C}\bar{t}$$), provided that the errors of the starting values also satisfy such a bound.

Theorem 3.2 limits the use of the BDF methods of orders higher than $2$ (and more severely for orders higher than $3$) to applications with a large damping parameter $\alpha$, such as cases described in [24, 39]. We remark, however, that in many situations $\alpha$ is of magnitude $10^{-2}$ or even smaller [10]. A very small damping parameter $\alpha$ affects not only the methods considered here. To our knowledge, the error analysis of any numerical method proposed in the literature breaks down as $\alpha\to 0$, as does the energy inequality.

It is not surprising that a positive lower bound on $\alpha$ arises for the methods of orders $k\geqslant 3$, since they are not A-stable and a lower bound on $\alpha$ is required also for the simplified linear problem $(\alpha+{\rm i})\partial_{t}u=\varDelta u$, which arises from (1.4) by freezing $\bm{m}$ in the term $\bm{m}\times\partial_{t}\bm{m}$ and diagonalizing this skew-symmetric linear operator (with eigenvalues $\pm{\rm i}$ and [math]) and by omitting the projection ${\mathbf{P}}(\bm{m})$ on the right-hand side of (1.4).

The proof of Theorem 3.2 uses a variant of the Nevanlinna–Odeh multiplier technique [34], which is described in the Appendix for the convenience of the reader. While for sufficiently large $\alpha$ we have an optimal-order error bound in the case of a smooth solution, there is apparently no discrete energy inequality under weak regularity assumptions similar to Proposition 3.1 for the BDF methods of orders $3$ to $5$.

As in Remark 3.1, the error bounds also allow us to bound the discrepancy from normality.

4. A continuous perturbation result

In this section we present a perturbation result for the continuous problem, because we will later transfer the arguments of its proof to the discretizations to prove stability and convergence of the numerical methods.

Let $\bm{m}(t)$ be a solution of (1.4) for $0\leqslant t\leqslant\bar{t}$, and let $\bm{m}_{\star}(t)$, also of unit length, solve the same equation up to a defect $\bm{d}(t)$ for $0\leqslant t\leqslant\bar{t}$:

[TABLE]

with

[TABLE]

Then, $\bm{m}_{\star}$ also solves the perturbed weak formulation

[TABLE]

and the error $\bm{e}=\bm{m}-\bm{m}_{\star}$ satisfies the error equation

[TABLE]

Before we turn to the perturbation result, we need Lipschitz-type bounds for the orthogonal projection ${\mathbf{P}}(\bm{m})={\mathbf{I}}-\bm{m}\bm{m}^{T}$ applied to sufficiently regular functions.

Lemma 4.1**.**

The projection ${\mathbf{P}}(\cdot)$ satisfies the following estimates, for functions $\bm{m},\bm{m}_{\star},\bm{v}:\varOmega\to{\mathbb{R}}^{3}$, where $\bm{m}$ and $\bm{m}_{\star}$ take values on the unit sphere and $\bm{m}_{\star}\in W^{1,\infty}(\varOmega)^{3}$:**

[TABLE]

Proof.

Setting $\bm{e}=\bm{m}-\bm{m}_{\star}$, we start by rewriting

[TABLE]

The first inequality then follows immediately by taking the $L^{2}$ norm of both sides of the above equality, using the fact that $\bm{m}$ and $\bm{m}_{\star}$ are of unit length. The second inequality is proved similarly, using the product rule

[TABLE]

the $L^{\infty}$ bound of $\partial_{i}\bm{m}_{\star}$, and the fact that $\|\bm{e}\|_{L^{\infty}}\leqslant\|\bm{m}\|_{L^{\infty}}+\|\bm{m}_{\star}\|_{L^{\infty}}\leqslant 2$. ∎

We have the following perturbation result.

Lemma 4.2**.**

Let $\bm{m}(t)$ and $\bm{m}_{\star}(t)$ be solutions of unit length of (1.5) and (4.1), respectively, and suppose that, for $0\leqslant t\leqslant\bar{t}$, we have

[TABLE]

Then, the error $\bm{e}(t)=\bm{m}(t)-\bm{m}_{\star}(t)$ satisfies, for $0\leqslant t\leqslant\bar{t}$,

[TABLE]

where the constant $C$ depends only on $\alpha,R,K$, and $\bar{t}$.

Proof.

Let us first assume that $\partial_{t}\bm{m}(t)\in H^{1}(\varOmega)^{3}$ for all $t$. Following [21], we test in the error equation (4.2) with $\bm{\varphi}={\mathbf{P}}(\bm{m})\partial_{t}\bm{e}\in T(\bm{m})$. By the following argument, this test function is then indeed in $H^{1}(\varOmega)^{3}$ and can be viewed as a perturbation of $\partial_{t}\bm{e}$:

[TABLE]

and so we have

[TABLE]

By Lemma 4.1 and using (4.3) we have

[TABLE]

Testing the error equation (4.2) with $\bm{\varphi}=\partial_{t}\bm{e}+\bm{q}$, we obtain

[TABLE]

where, by (4.1) and Lemma 4.1 with (4.3), $\bm{r}$ is bounded as

[TABLE]

By collecting terms, and using the fact that $(\bm{m}\times\partial_{t}\bm{e},\partial_{t}\bm{e})$ vanishes, we altogether obtain

[TABLE]

For the right-hand side, the Cauchy–Schwarz inequality and $\|\bm{m}\|_{L^{\infty}}=1$ yield

[TABLE]

Young’s inequality and absorptions, together with the bounds in (4.6) and (4.7), yield

[TABLE]

Here, we note that

[TABLE]

Combining these inequalities and integrating in time, we obtain

[TABLE]

By Gronwall’s inequality, we then obtain the stated error bound.

Finally, if $\partial_{t}\bm{m}(t)$ is not in $H^{1}(\varOmega)^{3}$ for some $t$, then a regularization and density argument, which we do not present here, yields the result, since the error bound does not depend on the $H^{1}$ norm of $\partial_{t}\bm{m}$. ∎

5. Orthogonal projection onto the discrete tangent space

For consistency and stability of the full discretization, we need to study properties of the $L^{2}(\varOmega)$-orthogonal projection onto the discrete tangent space $T_{h}(\bm{m})$, which we denote by

[TABLE]

We do not have an explicit expression for this projection, but the properties stated in Lemmas 5.1 to 5.3 will be used for proving consistency and stability. We recall that we consider a quasi-uniform, shape-regular family $\mathcal{T}_{h}$ of triangulations with Lagrange finite elements of polynomial degree $r$.

The first lemma states that the projection ${\mathbf{P}}_{h}(\bm{m})$ approximates the orthogonal projection ${\mathbf{P}}(\bm{m})=\mathbf{I}-\bm{m}\bm{m}^{T}$ onto the tangent space $T(\bm{m})$ with optimal order. It will be used in the consistency error analysis of Section 6.

Lemma 5.1**.**

For $\bm{m}\in W^{r+1,\infty}(\varOmega)^{3}$ with $|\bm{m}|=1$ almost everywhere we have

[TABLE]

for all $\bm{v}\in H^{r+1}(\varOmega)^{3}$, where $C$ depends on a bound of $\|\bm{m}\|_{W^{r+1,\infty}(\varOmega)^{3}}$.

The second lemma states that the projection ${\mathbf{P}}_{h}(\bm{m})$ has Lipschitz bounds of the same type as those of the orthogonal projection ${\mathbf{P}}(\bm{m})$ given in Lemma 4.1. It will be used in the stability analysis of Sections 7 and 8.

Lemma 5.2**.**

Let $\bm{m}\in W^{1,\infty}(\varOmega)^{3}$ and $\widetilde{\bm{m}}\in H^{1}(\varOmega)^{3}$ with $|\bm{m}|=|\widetilde{\bm{m}}|=1$ almost everywhere and $\|\bm{m}\|_{W^{1,\infty}}\leqslant R$. There exist $C_{R}>0$ and $h_{R}>0$ such that for $h\leqslant h_{R}$, for all ${\bm{v}_{h}}\in V_{h}^{3}$,

[TABLE]

The next lemma shows the $W^{s,p}$-stability of the projection. It is actually used for $p=2$ in the proof of Lemmas 5.1 and 5.2 and will be used for $p=2$ in Section 6 and for $p=\infty$ in Sections 7 and 8.

Lemma 5.3**.**

There exists a constant depending only on $p\in[1,\infty]$ and the shape regularity of the mesh such that for all $\bm{m}\in W^{1,\infty}(\varOmega)^{3}$ with $|\bm{m}|=1$ almost everywhere,

[TABLE]

for all $\bm{v}_{h}\in V_{h}^{3}$ and $s\in\{-1,0,1\}$.

These three lemmas will be proved in the course of this section, in which we formulate also three more lemmas that are of independent interest but will not be used in the following sections.

In the following, we use the dual norms

[TABLE]

The space $W^{-1,1}(\varOmega)$ is not the dual space of $W^{1,\infty}(\varOmega)$ but rather defined as the closure of $L^{2}(\varOmega)$ with respect to the norm $\|\cdot\|_{W^{-1,1}}$. We also recall that $\varPi_{h}\colon W^{s,p}(\varOmega)$ $\to W^{s,p}(\varOmega)$ is uniformly bounded for $s\in\{0,1\}$ and $p\in[1,\infty]$ (see, e.g., [20] for proofs in a much more general setting). By duality, we also obtain uniform boundedness for $s=-1$ and $p\in[1,\infty]$. A useful consequence is that for $v_{h}\in V_{h}$,

[TABLE]

Lemma 5.4**.**

There holds $\|v\|_{W^{s,p}(\varOmega)}\simeq\sup_{w\in W^{-s,q}(\varOmega)}\frac{(v,w)}{\|w\|_{W^{-s,q}(\varOmega)}}$ with $1/p+1/q=1$ for $p\in[1,\infty]$ and $s\in\{-1,0,1\}$.

Proof.

The interesting case is $(s,p)=(1,\infty)$ since all other cases follow by duality. For $v\in W^{1,\infty}(\varOmega)$, there exists a sequence of functions $\bm{q}_{n}\in C^{\infty}_{0}(\varOmega)^{3}$ with $\|\bm{q}_{n}\|_{L^{1}}=1$ such that

[TABLE]

Moreover, there holds

[TABLE]

Combining the last two estimates shows

[TABLE]

Since

[TABLE]

we conclude the proof. ∎

Let the discrete normal space $N_{h}(\bm{m}):=V_{h}^{3}\ominus T_{h}(\bm{m})$ be given as the $L^{2}$-orthogonal complement of $T_{h}(\bm{m})$ in $V_{h}^{3}$. We note that

[TABLE]

by the definition of $T_{h}(\bm{m})$. The functions in the discrete normal space are bounded from below as follows.

Lemma 5.5**.**

For every $R>0$, there exist $h_{R}>0$ and $c>0$ such that for all $\bm{m}\in W^{1,\infty}(\varOmega)^{3}$ with $|\bm{m}|=1$ almost everywhere and $\|\bm{m}\|_{W^{1,\infty}(\varOmega)}\leqslant R$ and for all $h\leqslant h_{R}$,

[TABLE]

for all $\psi_{h}\in V_{h}$ and $(s,p)\in\{-1,0,1\}\times[1,\infty]$.

Proof.

(a) We first prove the result for $s\in\{-1,0\}$. Let $\mathbf{I}_{h}\colon C(\varOmega)\to V_{h}^{3}$ denote the nodal interpolation operator and define $\bm{m}_{h}:=\mathbf{I}_{h}\bm{m}\in V_{h}^{3}$.

There holds

[TABLE]

Moreover, stability of ${\mathbf{\Pi}}_{h}$ in $L^{p}(\varOmega)^{3}$, for $1\leqslant p\leqslant\infty$, see [20], implies the estimate

[TABLE]

In turn, this implies

[TABLE]

For each element, the approximation properties of $\mathbf{I}_{h}$ show

[TABLE]

Thus, multiple inverse estimates yield

[TABLE]

Moreover, we have

[TABLE]

provided that $\|\bm{m}-\bm{m}_{h}\|_{L^{\infty}}\leqslant\tfrac{1}{2}$, which in view of

[TABLE]

is satisfied for $h\leqslant h_{R}$ with a sufficiently small $h_{R}>0$ that depends only on $R$. Altogether, this shows

[TABLE]

for $h\leqslant h_{R}$. Similarly we estimate

[TABLE]

Altogether, we obtain

[TABLE]

for $h\leqslant h_{R}$. This concludes the proof for $s=0$. Finally, for $s=-1$ we note that by using the result for $s=0$ and an inverse inequality,

[TABLE]

Since $\|\bm{m}\psi_{h}\|_{W^{-1,p}}\gtrsim\|\bm{m}\|_{W^{1,\infty}}^{-1}\|\psi_{h}\|_{W^{-1,p}}$, this concludes the proof for $s\in\{-1,0\}$.

(b) It remains to prove the result for $s=1$. Note that the result follows from duality if we show

[TABLE]

for all $\bm{w}_{h}\in N_{h}(\bm{m})$. To see this, note that (5.2) implies

[TABLE]

where we used in the second to last equality that part (a) for $s=0$ already shows that $\dim(N_{h}(\bm{m}))=\dim(V_{h})$ and since (5.2) implies that the map $N_{h}(\bm{m})\to V_{h},\,\bm{w}_{h}\mapsto\varPi_{h}(\bm{m}\cdot\bm{w}_{h})$ is injective, it is already bijective. It remains to prove (5.2). To that end, we first show for $\bm{w}_{h}={\mathbf{\Pi}}_{h}(\bm{m}\omega_{h})\in N_{h}(\bm{m})$ for some $\omega_{h}\in V_{h}$, using the reverse triangle inequality, that

[TABLE]

With $\bm{m}_{h}:=\mathbf{I}_{h}(\bm{m})\in V_{h}^{3}$, the last term satisfies

[TABLE]

where we used the same arguments as in the proof of part (a) to get the estimate $\|(\mathbf{I}-{\mathbf{\Pi}}_{h})(\bm{m}_{h}\omega_{h})\|_{L^{q}}\lesssim h\|\bm{m}_{h}\|_{W^{1,\infty}}\|\omega_{h}\|_{L^{q}}$. The fact $\|\bm{m}_{h}\|_{W^{1,\infty}}\lesssim\|\bm{m}\|_{W^{1,\infty}}$, the approximation property $\|\bm{m}-\bm{m}_{h}\|_{L^{\infty}}\lesssim h\|\bm{m}\|_{W^{1,\infty}}$, and an inverse inequality conclude

[TABLE]

with (hidden) constants depending only on $\|\bm{m}\|_{W^{1,\infty}}$ and shape regularity of the mesh.

To prove (5.2), it remains to bound the left-hand side above by $\|\varPi_{h}(\bm{m}\cdot\bm{w}_{h})\|_{W^{-1,q}}$. To that end, we note

[TABLE]

where we used part (a) for $s=0$ for the last inequality. An inverse inequality and the combination with (5.3) imply (5.2) for $h>0$ sufficiently small in terms of $\|\bm{m}\|_{W^{1,\infty}}^{-1}$. This concludes the proof. ∎

Lemma 5.6**.**

Define the matrix $M\in{\mathbb{R}}^{N\times N}$, where $N$ denotes the dimension of $V_{h}$, by $M_{ij}:=h^{-3}({\mathbf{\Pi}}_{h}(\bm{m}\phi_{j}),{\mathbf{\Pi}}_{h}(\bm{m}\phi_{i}))$. Under the assumptions of Lemma 5.5, there exists $C>0$ such that for $h\leqslant h_{R}$,

[TABLE]

where $C$ depends only on the shape regularity.

Proof.

Lemma 5.5 shows for $x\in{\mathbb{R}}^{N}$

[TABLE]

where $|\cdot|$ denotes the Euclidean norm on ${\mathbb{R}}^{N}$. Let $d(i,j):={\rm dist}(z_{i},z_{j})h^{-3}$ denote the metric which (approximately) measures the number of elements between the supports of $\phi_{i}$ and $\phi_{j}$, corresponding to the nodes $z_{i}$ and $z_{j}$, and let $B_{d}(z)$ denote the corresponding ball. In the following, we use a localization property of the $L^{2}$-projection, i.e., there exist $a,b>0$ such that for all $\ell\in{\mathbb{N}}$,

[TABLE]

The proof of this bound is essentially contained in the proof of [9, Lemma 3.1]. Since we use the very same arguments below, we briefly recall the strategy: First, one observes that the mass matrix $\widetilde{M}\in{\mathbb{R}}^{N\times N}$ with entries $\widetilde{M}_{ij}:=h^{-3}(\phi_{j},\phi_{i})$ is banded in the sense that $d(i,j)\gtrsim 1$ implies $\widetilde{M}_{ij}=0$, and it satisfies $\widetilde{M}x\cdot x\gtrsim|x|^{2}$. As shown below, this implies that the inverse matrix $\widetilde{M}^{-1}$ satisfies $|(\widetilde{M}^{-1})_{ij}|\lesssim{\rm e}^{-bd(i,j)}$ for some $b>0$ independent of $h>0$. Note that each entry of the vector field ${\mathbf{\Pi}}_{h}(\bm{m}\phi_{i})\in V_{h}^{3}$ can be represented by $\sum_{j=1}^{N}x_{k,j}\phi_{j},k=1,2,3,$ and is computed by solving $\widetilde{M}x_{k}=g_{k}\in{\mathbb{R}}^{N}$ with $\bm{m}=(m_{1},m_{2},m_{3})^{T}$ and $g_{k,j}:=(m_{k}\phi_{i},\phi_{j})$. Hence, the exponential decay of $\widetilde{M}^{-1}$ directly implies (5.5).

From the decay property (5.5), we immediately obtain

[TABLE]

for all $1\leqslant i,j\leqslant N$ and some $\widetilde{a},\widetilde{b}>0$. This already proves $\|M\|_{p}\leqslant C$. We follow the arguments from [28] to show that also $M^{-1}$ decays exponentially. To that end, note that (5.4) implies the existence of $c>0$ such that $\|I-cM\|_{2}=:q<1$ and hence

[TABLE]

Clearly, $I-cM$ inherits the decay properties from $M$ and therefore

[TABLE]

The value of $\max_{s=1,\dotsc,N}\sum_{r=1}^{N}{\rm e}^{-\widetilde{b}d(s,r)/2}$ depends only on the shape regularity of the triangulation and on $\widetilde{b}$, but is independent of $h$ (it just depends on the number of elements contained in an annulus of thickness $\approx h$). This implies the existence of $\widetilde{c}\geqslant 1$ such that

[TABLE]

Thus, for $\widetilde{c}^{k+1}\leqslant{\rm e}^{\widetilde{b}d(i,j)/4}$, we have $|((I-cM)^{k+1})_{ij}|\leqslant{\rm e}^{-\widetilde{b}d(i,j)/4}$, whereas for $\widetilde{c}^{k+1}>{\rm e}^{\widetilde{b}d(i,j)/4}$, we have $|((I-cM)^{k+1})_{ij}|\leqslant q^{k+1}<q^{\widetilde{b}d(i,j)/(4\log(\widetilde{c}))}$. Altogether, we find some $\widetilde{b}>0$ (we reuse the symbol), independent of $h$ such that

[TABLE]

Plugging this into (5.6), we obtain

[TABLE]

This yields the stated result. ∎

We are now in a position to prove Lemma 5.3.

Proof of Lemma 5.3.

(a) We first consider the case $s=0$. In view of (5.1), we write $(\mathbf{I}-{\mathbf{P}}_{h}(\bm{m}))\bm{v}_{h}\in N_{h}(\bm{m})$ as

[TABLE]

for some coefficient vector $x\in{\mathbb{R}}^{N}$ and let $b_{i}:=h^{-3/2}(\bm{v}_{h},\bm{m}\phi_{i})$ for $i=1,\dotsc,N$. Then, there holds $Mx=b$ with the matrix $M$ from Lemma 5.6. This lemma and the $L^{p}$-stability of the $L^{2}$-orthogonal projection $\Pi_{h}$ [20] imply that for $p\in[1,\infty]$,

[TABLE]

With $|b_{i}|\leqslant h^{-3/2}\|\bm{v}_{h}\|_{L^{p}(\operatorname{supp}(\phi_{i}))^{3}}h^{3(1-1/p)}=\|\bm{v}_{h}\|_{L^{p}(\operatorname{supp}(\phi_{i}))^{3}}h^{3/2-3/p}$, this shows

[TABLE]

(b) We now turn to the cases $s=\pm 1$. Define the operator

[TABLE]

and note that $\widetilde{\mathbf{P}}_{h}^{\perp}(\bm{m})\bm{v}_{h}\in N_{h}(\bm{m})$ as well as $\operatorname{ker}\widetilde{\mathbf{P}}_{h}^{\perp}(\bm{m})=T_{h}(\bm{m})$ (due to Lemma 5.5). However, $\widetilde{\mathbf{P}}_{h}^{\perp}(\bm{m})$ is no projection. We observe for $\bm{v}_{h}={\mathbf{\Pi}}_{h}(\bm{m}\psi_{h})\in N_{h}(\bm{m})$ that

[TABLE]

With Lemma 5.5 we conclude

[TABLE]

Since $\widetilde{\mathbf{P}}_{h}^{\perp}(\bm{m}){\mathbf{P}}_{h}(\bm{m})=0$ by definition of $T_{h}(\bm{m})$, we obtain with part (a) and an inverse inequality that for all $\bm{v}_{h}\in V_{h}^{3}$,

[TABLE]

The $W^{-1,p}(\varOmega)$-stability of $\varPi_{h}$ implies $\|\widetilde{\mathbf{P}}_{h}^{\perp}(\bm{m})\bm{v}_{h}\|_{W^{-1,p}}\lesssim\|\bm{m}\|_{W^{1,\infty}}^{2}\|\bm{v}_{h}\|_{W^{-1,p}}$ and the triangle inequality concludes the proof for $s=-1$. The case $s=1$ follows by duality. ∎

Proof of Lemma 5.2.

(a) ($s=0$) The projection $\bm{v}_{h}:={\mathbf{P}}_{h}(\bm{m})\bm{v}$ is given by the equation

[TABLE]

which in view of the definition of $T_{h}(\bm{m})$ is equivalent to the solution of the saddle point problem (with the Lagrange multiplier $\lambda_{h}\in V_{h}$)

[TABLE]

By the first equation, we also obtain the identity ${\mathbf{\Pi}}_{h}(\bm{m}\lambda_{h})=(\mathbf{I}-{\mathbf{P}}_{h}(\bm{m}))\bm{v}_{h}$, which will be used below. Furthermore, $\widetilde{\bm{v}}_{h}:={\mathbf{P}}_{h}(\widetilde{\bm{m}})\bm{v}$ is given by the same system with $\widetilde{\bm{m}}$ in place of $\bm{m}$, yielding a corresponding Lagrange multiplier $\widetilde{\lambda}_{h}$. Hence, the differences $\bm{e}_{h}:=\bm{v}_{h}-\widetilde{\bm{v}}_{h}$ and $\delta_{h}:=\lambda_{h}-\widetilde{\lambda}_{h}$ satisfy

[TABLE]

The classical results on saddle-point problems (see [13, Proposition 2.1]) require two inf-sup conditions to be satisfied. First,

[TABLE]

holds uniformly in $h$ due to Lemma 5.5. Second,

[TABLE]

holds uniformly in $h$ due to the stability estimates from Lemma 5.3 (noting that $\bm{v}_{h}={\mathbf{P}}_{h}(\bm{m})\bm{v}_{h}$ and $\bm{w}_{h}={\mathbf{P}}_{h}(\bm{m})\bm{w}_{h}$ for $\bm{v}_{h},\bm{w}_{h}\in T_{h}(\bm{m})$). For the above saddle-point problems, these bounds for $s=0$ give us an $L^{2}$ bound for $\bm{e}_{h}={\mathbf{P}}_{h}(\bm{m})\bm{v}-{\mathbf{P}}_{h}(\widetilde{\bm{m}})\bm{v}$: From [13] we obtain

[TABLE]

and

[TABLE]

With the stability from Lemma 5.3 and Lemma 5.5, we also obtain

[TABLE]

Altogether, this implies

[TABLE]

for $(p,q)\in\{(2,\infty),(\infty,2)\}$.

(b) ($s=1$) For the $H^{1}(\varOmega)$-estimate, we introduce the Riesz mapping $J_{h}$ between $V_{h}\subset H^{1}(\varOmega)$ and its dual $V_{h}\subset H^{1}(\varOmega)^{\prime}$, i.e., the isometry defined by

[TABLE]

By $\mathbf{J}_{h}:=\mathbf{I}\otimes J_{h}$ we denote the corresponding vector-valued mapping on $V_{h}^{3}$. We consider the bilinear form on ${V_{h}^{3}}\times{V_{h}^{3}}$ defined by

[TABLE]

and reformulate the saddle-point problem for $(\bm{v}_{h},\lambda_{h})\in V_{h}^{3}\times V_{h}\subset H^{1}(\varOmega)^{3}\times H^{1}(\varOmega)^{\prime}$ as

[TABLE]

As in the case $s=0$ (algebraically it is the same system), we have $\bm{v}_{h}={\mathbf{P}}_{h}(\bm{m})\bm{v}$ and ${\mathbf{\Pi}}_{h}(\bm{m}\lambda_{h})=(\mathbf{I}-{\mathbf{P}}_{h}(\bm{m}))\bm{v}$. The system for $\bm{e}_{h}=\bm{v}_{h}-\widetilde{\bm{v}}_{h}$ and $\delta_{h}=\lambda_{h}-\widetilde{\lambda}_{h}$ reads

[TABLE]

The above inf-sup bounds for $s=1$ and $s=-1$ are precisely the inf-sup conditions that need to be satisfied for these generalized saddle-point problems (see [15, Theorem 2.1]), whose right-hand sides are bounded by

[TABLE]

and

[TABLE]

As in the case $s=0$, we obtain from Lemma 5.3 and Lemma 5.5 that

[TABLE]

Hence, we obtain from [15, Theorem 2.1], for $(p,q)\in\{(2,\infty),(\infty,2)\}$,

[TABLE]

This implies the $H^{1}(\varOmega)^{3}$ estimate and hence concludes the proof. ∎

Proof of Lemma 5.1.

Since ${\mathbf{P}}_{h}(\bm{m})\bm{v}$ is the Galerkin approximation of the saddle point problem for ${\mathbf{P}}(\bm{m})\bm{v}$ (as in the previous proof), the Céa lemma for saddle-point problems (see [13, Theorem 2.1]) shows in $L^{2}$

[TABLE]

and similarly in $H^{1}$, using [15, Theorem 2.1],

[TABLE]

This concludes the proof. ∎

6. Consistency error and error equation

To study the consistency errors, we find it instructive to separate the issues of consistency for the time and space discretizations. Therefore, we first show defect estimates for the semidiscretization in time, and then turn to the full discretization.

6.1. Consistency error of the semi-discretization in time

The order of both the fully implicit $k$-step BDF method, described by the coefficients $\delta_{0},\dotsc,\delta_{k}$ and $1,$ and the explicit $k$-step BDF method, that is the method described by the coefficients $\delta_{0},\dotsc,\delta_{k}$ and $\gamma_{0},\dotsc,\gamma_{k-1},$ is $k,$ i.e.,

[TABLE]

We first rewrite the linearly implicit $k$-step BDF method (2.3) in strong form,

[TABLE]

with Neumann boundary conditions.

The consistency error $\bm{d}^{n}$ of the linearly implicit $k$-step BDF method (6.2) for the solution $\bm{m}$ is the defect by which the exact solution misses satisfying (6.2), and is given by

[TABLE]

for $n=k,\dotsc,N$, where we use the notation $\bm{m}^{n}_{\star}=\bm{m}(t_{n})$ and

[TABLE]

Note that the definition of $\dot{\bm{m}}_{\star}^{n}$ contains the projection ${\mathbf{P}}(\widehat{\bm{m}}_{\star}^{n})$, while $\dot{\bm{m}}^{n}$ was defined without a projection (see the first formula in (2.2)), since $\dot{\bm{m}}^{n}={\mathbf{P}}(\widehat{\bm{m}}^{n})\dot{\bm{m}}^{n}$ is automatically satisfied due to the constraint in (2.3).

The consistency error is bounded as follows.

Lemma 6.1**.**

If the solution of the LLG equation (1.4) has the regularity

[TABLE]

then the consistency error (6.3) is bounded by

[TABLE]

for $n=k,\dotsc,N$.

Proof.

We begin by rewriting the equation for the defect as

[TABLE]

In view of (1.4), we have

[TABLE]

and can rewrite (6.5) as

[TABLE]

i.e.,

[TABLE]

Therefore,

[TABLE]

with

[TABLE]

Now, in view of the first estimate in Lemma 4.1, we have

[TABLE]

i.e.,

[TABLE]

Therefore, it suffices to estimate $\dot{\bm{d}}^{n}$ and $\widehat{\bm{d}}^{n}$.

To estimate $\widehat{\bm{d}}^{n}$, we shall proceed in two steps. First we shall estimate the extrapolation error

[TABLE]

and then $\widehat{\bm{d}}^{n}.$

By Taylor expanding about $t_{n-k},$ the leading terms of order up to $k-1$ cancel, due to the second equality in (6.1), and we obtain

[TABLE]

with $\bm{m}^{(\ell)}:=\frac{\partial^{\ell}\bm{m}}{\partial t^{\ell}},$ whence

[TABLE]

Now, for a normalized vector $\bm{a}$ and a non-zero vector $\bm{b},$ we have

[TABLE]

whence

[TABLE]

Therefore, (6.11) yields

[TABLE]

To bound $\dot{\bm{d}}^{n},$ we use the fact that ${\mathbf{P}}(\bm{m}(t_{n}))\partial_{t}\bm{m}(t_{n})=\partial_{t}\bm{m}(t_{n})\in T(\bm{m}(t_{n}))$, so that we have

[TABLE]

By Lemma 4.1 and (6.12), we have for the last term

[TABLE]

By Taylor expanding the first term about $t_{n-k},$ we see that, due to the order conditions of the implicit BDF method, i.e., the first equality in (6.1), the leading terms of order up to $k-1$ cancel, and we obtain

[TABLE]

whence

[TABLE]

provided the solution $\bm{m}$ is sufficiently regular. Now, (6.6), (6.8), (6.14), and (6.12) yield

[TABLE]

This is the desired consistency estimate, which is valid for BDF methods of arbitrary order $k$. ∎

6.2. Consistency error of the full discretization

We define the Ritz projection $R_{h}\colon H^{1}(\varOmega)\to V_{h}$ corresponding to the Poisson–Neumann problem via

[TABLE]

for all $\psi\in V_{h}$, and we denote ${\mathbf{R}}_{h}={\mathbf{I}}\otimes R_{h}\colon H^{1}(\varOmega)^{3}\to V_{h}^{3}$. We denote again the $L^{2}$-orthogonal projections onto the finite element space by $\varPi_{h}\colon L^{2}(\varOmega)\to V_{h}$ and ${\mathbf{\Pi}}_{h}={\mathbf{I}}\otimes\varPi_{h}\colon L^{2}(\varOmega)^{3}\to V_{h}^{3}$. As in the previous section, we write ${\mathbf{P}}_{h}({\bm{m}})$ for the $L^{2}$-orthogonal projection onto the discrete tangent space at $\bm{m}$. We insert the following quantities, which are related to the exact solution,

[TABLE]

into the linearly implicit $k$-step BDF method (2.6) and obtain a defect $\bm{d}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{\star,h}^{n})$ from

[TABLE]

for all $\bm{\varphi}_{h}\in T_{h}(\widehat{\bm{m}}_{\star,h}^{n})$. By definition, there holds $(R_{h}\varphi,1)=(\varphi,1)$ (this can be seen by testing with $\psi=1$) and hence

[TABLE]

Thus, we obtain the consistency error for the full discretization by

[TABLE]

for $n=k,\dotsc,N$. The consistency error is bounded as follows.

Lemma 6.2**.**

If the solution of the LLG equation (1.4) has the regularity

[TABLE]

then the consistency error (6.18) is bounded by

[TABLE]

for $n$ with $k\tau\leqslant n\tau\leqslant\bar{t}$.

Proof.

We begin by defining

[TABLE]

and note that ${\mathbf{P}}(\bm{m}_{\star}^{n})\bm{D}^{n}=0$. Here we denote again $\bm{m}_{\star}^{n}=\bm{m}(t_{n})$ and in the following we use also the notations $\dot{\bm{m}}_{\star}^{n}$ and $\widehat{\bm{m}}_{\star}^{n}$ as defined in (6.4). With this, we rewrite the equation for the defect as

[TABLE]

For the term IV we have by Lemma 4.1

[TABLE]

where the last term $\widehat{\bm{m}}_{\star}^{n}-{\bm{m}}_{\star}^{n}$ has been bounded in the $L^{2}$ norm by $C\tau^{k}$ in the proof of Lemma 6.1.

The term III is estimated using the first bound from Lemma 5.1, under our regularity assumptions, as

[TABLE]

For the bound on II we use Lemma 5.2 ($i$) (with $p=2$ and $q=\infty$), to obtain

[TABLE]

where, using (7.11), we obtain

[TABLE]

The denominator is bounded from below by $1-C\tau^{k}$, because $|\bm{m}_{*}^{n}|=1$ and $|\sum_{i=1}^{k}\gamma_{i}\bm{m}_{*}^{n-i}-\bm{m}_{*}^{n}|\leqslant C\tau^{k}$. For the first term we have

[TABLE]

The terms $\|\partial_{t}\bm{m}(t_{n})-\dot{\bm{m}}_{\star}^{n}\|_{L^{2}}$ and $\|\bm{m}_{\star}^{n}\times\partial_{t}\bm{m}(t_{n})-\widehat{\bm{m}}_{\star}^{n}\times\dot{\bm{m}}_{\star}^{n}\|_{L^{2}}$ can be handled as in the proof of Lemma 6.1. Standard error estimates for the Ritz projection $R_{h}$ (we do not exploit the Aubin–Nitsche duality here) imply

[TABLE]

Together this yields, under the stated regularity assumption,

[TABLE]

and the result follows. ∎

6.3. Error equation

We recall, from (2.6), the fully discrete problem with the linearly implicit BDF method: find $\dot{\bm{m}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$ such that for all $\bm{\varphi}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$,

[TABLE]

Then, similarly as we have done in Section 4, we first rewrite (6.17): for all $\bm{\varphi}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$,

[TABLE]

with

[TABLE]

The error $\bm{e}_{h}^{n}=\bm{m}_{h}^{n}-\bm{m}_{\star,h}^{n}$ satisfies the error equation that is obtained by subtracting (6.20) from (6.19). We use the notations

[TABLE]

We have the following bound for $\bm{s}_{h}^{n}$.

Lemma 6.3**.**

Under the regularity assumptions of Lemma 6.2, we have

[TABLE]

Proof.

We use Lemmas 5.1 and 5.3, and the bounds in the proof of Lemma 6.2. We start by subtracting $(\mathbf{I}-{\mathbf{P}}(\widehat{\bm{m}}_{\star,h}^{n}))\partial_{t}\bm{m}_{\star}^{n}=0$, and obtain (with $\partial^{\tau}\bm{m}_{\star,h}^{n}:=\frac{1}{\tau}\sum_{j=0}^{k}\delta_{j}\bm{m}_{\star,h}^{n-j}$)

[TABLE]

The first term above is bounded as $O(\tau^{k}+h^{r})$ via the techniques of the consistency proofs, Lemma 6.1 and 6.2. For the second term we have

[TABLE]

where the first term is bounded as $O(\tau^{k}+h^{r})$, using Lemma 5.3 and the previous estimate, while the second term is bounded as $O(h^{r})$ by the $H^{1}$ estimate from Lemma 5.1. Altogether, we obtain the stated $H^{1}$ bound for $\bm{s}_{h}^{n}$. ∎

We then have the error equation

[TABLE]

for all $\bm{\varphi}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$, which is to be taken together with (6.21)–(6.23).

7. Stability of the full discretization for BDF of orders 1 and 2

For the A-stable BDF methods (those of orders 1 and 2) we obtain the following stability estimate, which is analogous to the continuous perturbation result Lemma 4.2.

Lemma 7.1** (Stability for orders $k=1,2$).**

Consider the linearly implicit $k$-step BDF discretization (2.6) for $k\leqslant 2$ with finite elements of polynomial degree $r\geqslant 1$. Let $\bm{m}_{h}^{n}$ and $\bm{m}_{\star,h}^{n}={\mathbf{R}}_{h}\bm{m}(t_{n})$ satisfy equations (2.6) and (6.17), respectively, and suppose that the exact solution $\bm{m}(t)$ is bounded by (4.3) and $\|\mathbf{H}(t)\|_{L^{\infty}}\leqslant M$ for $0\leqslant t\leqslant\bar{t}$. Then, for sufficiently small $h\leqslant\bar{h}$ and $\tau\leqslant\bar{\tau}$, the error $\bm{e}_{h}^{n}=\bm{m}_{h}^{n}-\bm{m}_{\star,h}^{n}$ satisfies the following bound, for $k\tau\leqslant n\tau\leqslant\bar{t}$,

[TABLE]

*where the constant $C$ is independent of $h,\tau$ and $n$, but depends on $\alpha,R,K,M$, and $\bar{t}$. This estimate holds under the smallness condition that the right-hand side is bounded by $\hat{c}h$ with a sufficiently small constant $\hat{c}$ $($note that the right-hand side is of size $O((\tau^{k}+h^{r})^{2})$ in the case of a sufficiently regular solution$)$. *

Combining Lemmas 7.1, 6.2 and 6.3 yields the proof of Theorem 3.1: These lemmas imply the estimate

[TABLE]

in the case of a sufficiently regular solution. Since then $\|{\mathbf{R}}_{h}\bm{m}(t_{n})-\bm{m}(t_{n})\|_{H^{1}(\varOmega)^{3}}\leqslant Ch^{r}$ and because of $\bm{m}_{h}^{n}-\bm{m}(t_{n})=\bm{e}_{h}^{n}+({\mathbf{R}}_{h}\bm{m}(t_{n})-\bm{m}(t_{n}))$, this implies the error bound (3.1).

The smallness condition imposed in Lemma 7.1 is satisfied under the very mild CFL condition, for a sufficiently small $\bar{c}>0$ (independent of $h,\tau$ and $n$),

[TABLE]

Taken together, this proves Theorem 3.1.

Proof.

(a) Preparations. The proof of this lemma transfers the arguments of the proof of Lemma 4.2 to the fully discrete situation, using energy estimates obtained by testing with (essentially) the discrete time derivative of the error, as presented in the Appendix, which is based on Dahlquist’s $G$-stability theory.

However, testing the error equation (6.25) directly with $\dot{\bm{e}}_{h}^{n}$ is not possible, since $\dot{\bm{e}}_{h}^{n}$ is not in the tangent space $T_{h}(\widehat{\bm{m}}_{h}^{n})$. Therefore, as in the proof of Lemma 4.2, we again start by showing that the test function $\bm{\varphi}_{h}={\mathbf{P}}_{h}(\widehat{\bm{m}}_{h}^{n})\dot{\bm{e}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n})\cap H^{1}(\varOmega)^{3}$ is a perturbation of $\dot{\bm{e}}_{h}^{n}$ itself:

[TABLE]

Here we note that ${\mathbf{P}}_{h}(\widehat{\bm{m}}_{h}^{n})\dot{\bm{m}}_{h}^{n}=\dot{\bm{m}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$ by construction of the method (2.6), and ${\mathbf{P}}_{h}(\widehat{\bm{m}}_{\star,h}^{n})\dot{\bm{m}}_{\star,h}^{n}=\dot{\bm{m}}_{\star,h}^{n}\in T_{h}(\widehat{\bm{m}}_{\star,h}^{n})$ by the definition of $\dot{\bm{m}}_{\star,h}^{n}$ in (6.4). So we have

[TABLE]

and hence

[TABLE]

The proof now transfers the proof of the continuous perturbation result Lemma 4.2 to the discrete situation with some notable differences, which are emphasized here:

(i) Instead of using the continuous quantities it uses their spatially discrete counterparts, in particular the discrete projections ${\mathbf{P}}_{h}(\widehat{\bm{m}}_{h}^{n})$ and ${\mathbf{P}}_{h}(\widehat{\bm{m}}_{\star,h}^{n})$, defined and studied in Section 5. In view of the definition (2.1) and (6.16) of $\widehat{\bm{m}}_{h}^{n}$ and $\widehat{\bm{m}}_{\star,h}^{n}$, respectively, this requires that $\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{h}^{n-j-1}(x)$ and $\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{\star,h}^{n-j-1}(x)$ are bounded away from zero uniformly for all $x\in\varOmega$.

(ii) Instead of Lemma 4.1 we use Lemma 5.2 (with $\widehat{\bm{m}}_{h}^{n}$ and $\widehat{\bm{m}}_{\star,h}^{n}$ in the role of $\widetilde{\bm{m}}$ and $\bm{m}$, respectively) to bound the quantity $\bm{q}_{h}^{n}$. This requires that $\widehat{\bm{m}}_{\star,h}^{n}$ and $\dot{\bm{m}}_{\star,h}^{n}$ are bounded in $W^{1,\infty}$ independently of $h$.

Ad (i): In order to show that $|\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{h}^{n-j-1}(x)|$ stays close to $1$ for all $x\in\varOmega$, we need to establish an $L^{\infty}$ bound for the errors $\bm{e}_{h}^{n-j-1}=\bm{m}_{h}^{n-j-1}-\bm{m}_{\star,h}^{n-j-1}$.

We use an induction argument and assume that for some time step number $\bar{n}$ with $\bar{n}\tau\leqslant\bar{t}$ we have

[TABLE]

where we choose $\rho$ sufficiently small independent of $h$ and $\tau$. (In this proof it suffices to choose $\rho\leqslant 1/(4C_{\gamma})$, where $C_{\gamma}=\sum_{j=0}^{k-1}|\gamma_{j}|=2^{k}-1$.)

Note that the smallness condition of the lemma implies that (7.3) is satisfied for $\bar{n}=k$, because for the $L^{\infty}$ errors of the starting values we have by an inverse inequality, for $i=0,\dotsc,k-1$,

[TABLE]

provided that $\hat{c}$ is sufficiently small (independent of $\tau$ and $h$), as is assumed.

We will show in part (b) of the proof that with the induction hypothesis (7.3) we obtain also $\|\bm{e}_{h}^{\bar{n}}\|_{L^{\infty}}\leqslant\rho$ so that finally we obtain (7.3) for all $\bar{n}$ with $\bar{n}\tau\leqslant\bar{t}$.

Using reverse and ordinary triangle inequalities, the error bound of [12, Corollary 8.1.12] (noting that $\bm{m}(t)\in W^{2,\infty}(\varOmega)$ under our assumptions) and the $L^{\infty}$ boundedness of $\partial_{t}\bm{m}$, and the bound (7.3), we estimate

[TABLE]

provided that $h$ and $\tau$ are sufficiently small. The same argument also yields that \bigl{\|}|\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{\star,h}^{n-j-1}|-1\bigr{\|}_{L^{\infty}}\leqslant\frac{1}{2}, and so we have

[TABLE]

for all $x\in\varOmega$. In particular, it follows that $\widehat{\bm{m}}_{h}^{n}$ and $\widehat{\bm{m}}_{\star,h}^{n}$ are unambiguously defined.

Ad (ii): The required $W^{1,\infty}$ bound for $\bm{m}_{\star,h}^{n}={\mathbf{R}}_{h}\bm{m}(t_{n})$ follows from the $W^{1,\infty}$-stability of the Ritz projection: by [12, Theorem 8.1.11] and by the assumed $W^{1,\infty}$ bound (4.3) for $\bm{m}(t)$,

[TABLE]

The bounds (7.5) and (7.6) for $n\leqslant\bar{n}$ imply that also

[TABLE]

for $n\leqslant\bar{n}$ (with a different constant $C$). Using this bound in Lemma 5.3 and the assumed $W^{1,\infty}$ bound (4.3) for $\partial_{t}\bm{m}(t)$, we obtain with $\delta(\zeta)/(1-\zeta)=\sum_{\ell=1}^{k}(1-\zeta)^{\ell-1}/\ell=:\sum_{j=0}^{k-1}\mu_{j}\zeta^{j}$ that

[TABLE]

We can now establish a bound for $\bm{q}_{h}^{n}$ as defined in (7.2), using Lemma 5.2 together with the above $W^{1,\infty}$ bounds for $\widehat{\bm{m}}_{\star,h}^{n}$ and $\dot{\bm{m}}_{\star,h}^{n}$ to obtain

[TABLE]

With the $W^{1,\infty}$ bound of $\widehat{\bm{m}}_{\star,h}^{n}$ we also obtain a bound of $\bm{r}_{h}^{n}$ defined in (6.21). Using Lemma 5.2 ($i$) and recalling the $L^{\infty}$ bound of $\varDelta\bm{m}+\bm{H}$ of (4.3), we find that $\bm{r}_{h}^{n}$ is bounded by

[TABLE]

(b) Energy estimates. For $n\leqslant\bar{n}$ with $\bar{n}$ of (7.3), we test the error equation (6.25) with $\bm{\varphi}_{h}=\dot{\bm{e}}_{h}^{n}+\bm{q}_{h}^{n}$ and obtain

[TABLE]

By collecting the terms, and using the fact that $(\widehat{\bm{m}}_{h}^{n}\times\dot{\bm{e}}_{h}^{n},\dot{\bm{e}}_{h}^{n})=0$, we altogether obtain

[TABLE]

We now estimate the term $(\nabla\bm{e}_{h}^{n},\nabla\dot{\bm{e}}_{h}^{n})$ on the left-hand side from below using Dahlquist’s Lemma 10.1, so that the ensuing relation (8.7b) yields

[TABLE]

where ${\mathbf{E}}_{h}^{n}=(\bm{e}_{h}^{n-k+1},\dotsc,\bm{e}_{h}^{n})$ and the $G$-weighted semi-norm is given by

[TABLE]

This semi-norm satisfies the relation

[TABLE]

where $\gamma^{-}$ and $\gamma^{+}$ are the smallest and largest eigenvalues of the positive definite symmetric matrix $G=(g_{ij})$ from Lemma 10.1.

The remaining terms are estimated using the Cauchy–Schwarz inequality and $\|\widehat{\bm{m}}_{h}^{n}\|_{L^{\infty}}=1$; we altogether obtain

[TABLE]

We now show an $L^{2}$ error bound for $\widehat{\bm{e}}_{h}^{n}$ in terms of $(\bm{e}_{h}^{n-j-1})_{j=0}^{k-1}$. Using the fact that for $\bm{a},\bm{b}\in{\mathbb{R}}^{3}\setminus\{0\}$,

[TABLE]

and the lower bounds in (7.5) for both $|\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{h}^{n-j-1}|$ and $|\sum_{j=0}^{k-1}\gamma_{j}\bm{m}_{\star,h}^{n-j-1}|$, we can estimate

[TABLE]

To show a similar bound for $\|\nabla\widehat{\bm{e}}_{h}^{n}\|_{L^{2}}$ we need the following two observations: First, the $W^{1,\infty}$ bounds for $\bm{m}_{\star,h}^{n-j-1}$ from (7.6) imply $W^{1,\infty}$ boundedness for $\widehat{\bm{m}}_{\star,h}^{n}$ by

[TABLE]

Second, similarly we have

[TABLE]

Combining these two observations, again with $\bm{m}_{h}$ and $\bm{m}_{\star,h}$ in the role of $\bm{a}$ and $\bm{b}$, respectively, and the upper and lower bounds from (7.5) altogether yield

[TABLE]

We estimate further using Young’s inequality and absorptions into the term $\|\dot{\bm{e}}^{n}\|_{L^{2}}^{2}$, together with the bounds in (7.8) and (7.9), to obtain

[TABLE]

Multiplying both sides by $\tau$, summing up from $k$ to $n\leqslant\bar{n}$, and using an absorption yield

[TABLE]

We then arrive, using (7.10), at

[TABLE]

with $c$ depending on $\alpha$.

Similarly as in the time continuous case in the proof of Lemma 4.2, we connect $\|\bm{e}_{h}^{n}\|_{L^{2}}^{2}$ and $\tau\sum_{j=k}^{n}\|\dot{\bm{e}}_{h}^{j}\|_{L^{2}}^{2}$. We rewrite the identity

[TABLE]

as

[TABLE]

with $\delta_{\ell}=0$ for $\ell>k$ and where

[TABLE]

depends only on the starting errors and satisfies $\bm{g}_{h}^{n}=0$ for $n\geqslant 2k$. With the inverse power series of $\delta(\zeta)$,

[TABLE]

we then have, for $n\geqslant k$,

[TABLE]

By the zero-stability of the BDF method of order $k\leqslant 6$, the coefficients $\kappa_{n}$ are uniformly bounded: $|\kappa_{n}|\leqslant c$ for all $n\geqslant 0$. Therefore we obtain via the Cauchy–Schwarz inequality

[TABLE]

Inserting this bound into (7.14) then yields

[TABLE]

and a discrete Gronwall inequality implies the stated stability result for $n\leqslant\bar{n}$. It then follows from this stability bound, the smallness condition of the lemma and the inverse estimate from $H^{1}$ to $L^{\infty}$ that (7.3) is satisfied also for $\bar{n}+1$. This completes the induction step for (7.3) and proves the stated error bound. ∎

8. Stability of the full discretization for BDF of orders 3 to 5

Stability for full discretizations using the BDF methods of orders $3$ to $5$ can be shown under additional conditions on the damping parameter $\alpha$ and the stepsize $\tau$.

Lemma 8.1** (Stability for orders $k=3,4,5$).**

Consider the linearly implicit $k$-step BDF discretization (2.6) for $3\leqslant k\leqslant 5$ with finite elements of polynomial degree $r\geqslant 2$. Let $\bm{m}_{h}^{n}$ and $\bm{m}_{\star,h}^{n}$ satisfy (2.6) and (6.17), respectively, and suppose that the regularity assumptions of Lemma 7.1 hold. Furthermore, assume that the damping parameter $\alpha$ satisfies

[TABLE]

with the multiplier $\eta_{k}$ of Lemma 10.2, and that $\tau$ and $h$ satisfy the mild CFL-type condition, for some $\bar{c}>0$,

[TABLE]

Then, for sufficiently small $h\leqslant\bar{h}$ and $\tau\leqslant\bar{\tau}$, the error $\bm{e}_{h}^{n}=\bm{m}_{h}^{n}-\bm{m}_{\star,h}^{n}$ satisfies the following bound, for $k\tau\leqslant n\tau\leqslant\bar{t}$,

[TABLE]

where the constant $C$ is independent of $\tau,h$ and $n$, but depends on $\alpha,R,K,M$, and exponentially on $\bar{c}\bar{t}$. This estimate holds under the smallness condition that the right-hand side is bounded by $\hat{c}h^{3}$ with a constant $\hat{c}$ $($note that the right-hand side is of size $O((\tau^{k}+h^{r})^{2})$ in the case of a sufficiently regular solution$)$.

Together with the defect bounds of Section 6, this stability lemma proves Theorem 3.2. We remark that the thresholds $\alpha_{k}>0$ defined here are the same as those appearing in Theorem 3.2.

Proof.

The proof of this lemma combines the arguments of the proof of Lemma 7.1 with a nonstandard variant of the multiplier technique of Nevanlinna and Odeh, as outlined in the Appendix. Since the size of the parameter $\alpha$ determines which BDF methods satisfy the stability estimate, the dependence on $\alpha$ will be carefully traced all along the proof.

(a) Preparations. As in the previous proof, we make again the induction hypothesis (7.3) for some $\bar{n}$ with $\bar{n}\tau\leqslant\bar{t}$, but this time with $\rho=c_{0}h$ for some positive constant $c_{0}$:

[TABLE]

By an inverse inequality, this implies that $\|\bm{e}_{h}^{n}\|_{W^{1,\infty}}$ has an $h$- and $\tau$-independent bound, and hence also $\|\bm{m}_{h}^{n}\|_{W^{1,\infty}}$ for $n<\bar{n}$. Together with (7.5), this implies

[TABLE]

and further

[TABLE]

As in the Appendix, we aim to subtract $\eta_{k}$ times the error equation for time step $n-1$ from the error equation for time step $n$, and then to test with $\bm{\varphi}_{h}={\mathbf{P}}_{h}(\widehat{\bm{m}}_{h}^{n})\dot{\bm{e}}^{n}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n})$ (similarly as in the proof of Lemma 7.1). However, this is not possible directly due to the different test spaces at different time steps:

[TABLE]

for all $\bm{\psi}_{h}\in T_{h}(\widehat{\bm{m}}_{h}^{n-1})$.

As in (7.2), we have

[TABLE]

where $\bm{q}_{h}^{n}$ is bounded by (7.8).

In turn, the test function $\bm{\psi}_{h}={\mathbf{P}}_{h}(\widehat{\bm{m}}_{h}^{n-1})\dot{\bm{e}}_{h}^{n}\in T_{h}(\widehat{\bm{m}}_{h}^{n-1})$ is a perturbation of $\bm{\varphi}_{h}=\dot{\bm{e}}_{h}^{n}+\bm{q}_{h}^{n}$, since using (8.7h) we obtain

[TABLE]

The perturbation $\bm{p}_{h}^{n}$ is estimated using the second bound in Lemma 5.2 ($i$) with $p=\infty$, $q=2$, and noting (8.5). We obtain

[TABLE]

We have $\|\mathbf{R}_{h}\partial_{t}\bm{m}(t)\|_{L^{\infty}}\leqslant c\|\partial_{t}\bm{m}(t)\|_{W^{1,\infty}}$ by [12, Theorem 8.1.11]. In view of (8.6) we obtain, for $\tau\leqslant\bar{C}h$,

[TABLE]

and by an inverse estimate,

[TABLE]

We also recall the bound (7.9) for $\|\bm{r}_{h}^{n}\|_{L^{2}}$.

(b) Energy estimates. By subtracting (8.7)$-\eta_{k}$(8.7) with the above choice of test functions, we obtain

[TABLE]

We estimate the terms of the error equation (8.7k) separately and track carefully the dependence on $\eta_{k}$ and $\alpha$.

The term $\alpha(\dot{\bm{e}}_{h}^{n}-\eta_{k}\dot{\bm{e}}_{h}^{n-1},\dot{\bm{e}}_{h}^{n})$ is bounded from below, using Young’s inequality and absorptions, by

[TABLE]

while the term $(\nabla\bm{e}_{h}^{n}-\eta_{k}\nabla\bm{e}_{h}^{n-1},\nabla\dot{\bm{e}}_{h}^{n})$ is bounded from below, via the relation (8.7b) and (6.23), by

[TABLE]

with ${\mathbf{E}}_{h}^{n}=(\bm{e}_{h}^{n-k+1},\dotsc,\bm{e}_{h}^{n})$, and where the $G$-weighted semi-norm is generated by the matrix $G=(g_{ij})$ from Lemma 10.1 for the rational function $\delta(\zeta)/(1-\eta_{k}\zeta)$.

The remaining terms outside the rectangular bracket are estimated using the Cauchy–Schwarz and Young inequalities (the latter often with a sufficiently small but fixed $h$- and $\tau$-independent weighting factor $\mu>0$) and $\|\widehat{\bm{m}}_{h}^{n}\|_{L^{\infty}}=1$ and orthogonality. We obtain, with varying constants $c$ (which depend on $\alpha$ and are inversely proportional to $\mu$)

[TABLE]

where in the last inequality we used (7.12) and (7.13) to estimate $\widehat{\bm{e}}_{h}^{n}$.

The terms inside the rectangular bracket are bounded similarly, using (8.7i) and (8.7j) and the condition $\tau\leqslant\bar{C}h$, by

[TABLE]

Here $\mu$ is an arbitrarily small positive constant (independent of $\tau$ and $h$), and $c$ depends on the choice of $\mu$.

In view of (7.9), the terms with the defects $\bm{r}_{h}^{n}$ are bounded by

[TABLE]

Combination of these inequalities yields

[TABLE]

Under condition (8.1) we have

[TABLE]

Multiplying both sides by $\tau$ and summing up from $k$ to $n$ with $n\leqslant\bar{n}$ yields, for sufficiently small $\mu$,

[TABLE]

The proof is then completed using exactly the same arguments as in the last part of the proof of Lemma 7.1, by establishing an estimate between $\|\bm{e}_{h}^{n}\|_{L^{2}}^{2}$ and $\tau\sum_{j=k}^{n}\|\dot{\bm{e}}_{h}^{j}\|_{L^{2}}^{2}$ and using a discrete Gronwall inequality, and completing the induction step for (8.4). ∎

9. Numerical experiments

To obtain significant numerical results, we prescribe the exact solution $\bm{m}$ on given three-dimensional domains $\varOmega:=[0,1]\times[0,1]\times[0,L]$ with $L\in\{1/100,1/4\}$. The discretizations of these domains will consist of a few layers of elements in $z$-direction (one layer for $L=1/100$ and ten layers for $L=1/4$) and a later specified number of elements in $x$ and $y$ directions. This mimics the common case of thin film alloys as for example in the standard problems of the Micromagnetic Modeling Activity Group at NIST Center for Theoretical and Computational Materials Science (ctcms.nist.gov). Moreover, this mesh structure helps to keep the computational requirements reasonable and allow us to compute the experiments on a desktop PC. We are aware that these experiments are only of preliminary nature and are just supposed to confirm the theoretical results. A more thorough investigation of the numerical properties of the developed method is needed. This will require us to incorporate preconditioning, parallelization of the computations, as well as lower order energy contributions in the effective field (1.3) to be able to compare to benchmark results from computational physics. This, however, is beyond the scope of this paper, and will be the topic of a subsequent work.

We consider the time interval $[0,\bar{t}\,]$ with $\bar{t}=0.2$ and define two different exact solutions. Since within our computational budget either the time discretization error or the space discretization error dominates, we construct the solutions such that the first one is harder to approximate in space, while the second one is harder to approximate in time. Both solutions are constant in $z$-direction as is often observed in thin-film applications.

9.1. Implementation

The numerical experiments were conducted using the finite element package FEniCS (www.fenicsproject.org) on a desktop computer. As already discussed in Section 2.2, there are several ways to implement the tangent space restriction. We decided to solve a saddle point problem (variant (a) in Section 2.2) for simplicity of implementation. For preconditioning, we used the black-box AMG preconditioner that comes with FEniCS. Although this might not be the optimal solution, it keeps the number of necessary iterative solver steps within reasonable bounds. Assuming perfect preconditioning, the cost per time-step is then proportional to the number of mesh-elements. We observed this behavior approximately, although further research beyond the scope of this work is required to give a definite conclusion.

9.2. Exact solutions

We choose the damping parameter $\alpha=0.2$ and define $g(t):=(\bar{t}+0.1)/(\bar{t}+0.1-t)$ as well as $d(x):=(x_{1}-1/2)^{2}+(x_{2}-1/2)^{2}$, which is the squared distance of the projection of $x$ to $[0,1]\times[0,1]$ and the point $(1/2,1/2)$. For some constant $C=400$ (a choice made to have pronounced effects), define

[TABLE]

It is easy to check that $|\bm{m}(x,t)|=1$ for all $(x,t)\in\varOmega\times[0,\bar{t}\,]$. Moreover, $\partial_{n}\bm{m}(x,t)=0$ for all $x\in\partial\varOmega$. We may calculate the time derivative of $\bm{m}$ in a straightforward fashion, i.e., $\partial_{t}\bm{m}(x,t)=0$ for $d(x)>1/4$ and

[TABLE]

Here, $\bm{m}_{3}$ denotes the third component of $\bm{m}$ as defined above.

The second exact solution is defined via

[TABLE]

Due to the polynomial nature in the first and the third component, and the well-behaved square-root, the space approximation error does not dominate the time approximation.

9.3. The experiments

We now may compute the corresponding forcings $\bm{H}$ resp. $\widetilde{\bm{H}}$ to obtain the prescribed solutions by inserting into (1.4), i.e.,

[TABLE]

(Note that we may disregard the projection $\mathbf{P}(\bm{m})$ from (1.4) since we solve in the tangent space anyway.) We compute $\bm{H}$ numerically by first interpolating $\bm{m}$ and $\partial_{t}\bm{m}$ and then computing the derivatives. This introduces an additional error which is not accounted for in the theoretical analysis. However, the examples below confirm the expected convergence rates and hence conclude that this additional perturbation is negligible. Figure 9.3 shows slices of the exact solution at different time steps. Figure 9.2 shows the convergence with respect to the time step size $\tau$, while Figure 9.3 shows convergence with respect to the spatial mesh size $h$. All the experiments confirm the expected rates for smooth solutions.

Figure 9.1. The first row shows the exact solution $\bm{m}(x,t)$ from (8.7a) for $x\in[0,1]\times[0,1]\times\{0\}$ and $t\in\{0,0.05,\bar{t}\}$ (from left to right), whereas the second row shows the exact solution $\widetilde{\bm{m}}(x,t)$ from (8.7b) for $x\in[0,1]\times[0,1]\times\{0\}$ and $t\in\{0,0.2/6,0.2/3\}$ (from left to right). While the problems are three-dimensional, the solutions are constant in $z$-direction and we only show one slice of the solution.

Finally, we consider an example with nonsmooth initial data and constant right-hand side. The initial data are given by

[TABLE]

With the constant forcing field $\bm{H}:=(0,1,1)^{T}$ we compute a numerical approximation to the unknown exact solution. Note that we do not expect any smoothness of the solution (even the initial data is not smooth). Figure 9.4 nevertheless shows a physically consistent decay of the energy $\|\nabla\bm{m}(t)\|_{L^{2}(\varOmega)^{3}}$ over time as well as a good agreement between different orders of approximation. Moreover, the computed approximation shows little deviation from unit length as would be expected for smooth solutions.

10. Appendix: Energy estimates for backward difference formulae

The stability proofs of this paper rely on energy estimates, that is, on the use of positive definite bilinear forms to bound the error $e$ in terms of the defect $d$. This is, of course, a basic technique for studying the time-continuous problem and also for backward Euler and Crank–Nicolson time discretizations (see, e.g., Thomée [38]), but energy estimates still appear to be not well known for backward difference formula (BDF) time discretizations of order up to $5$, which are widely used for solving stiff ordinary differential equations. To illustrate the basic mechanism, we here just consider the prototypical linear parabolic evolution equation in its weak formulation, given by two positive definite symmetric bilinear forms $(\cdot,\cdot)$ and $a(\cdot,\cdot)$ on Hilbert spaces $H$ and $V$ with induced norms $|\cdot|$ and $\|\cdot\|$, respectively, and with $V$ densely and continuously embedded in $H$. The problem then is to find $u(t)\in V$ such that

[TABLE]

with initial condition $u(0)=u_{0}$. If $u^{\star}$ is a function that satisfies the equation up to a defect $d$, that is,

[TABLE]

then the error $e=u-u^{\star}$ satisfies, in this linear case, an equation of the same form,

[TABLE]

with initial value $e_{0}=u_{0}-u_{0}^{\star}$. Testing with $v=e$ yields

[TABLE]

Estimating the right-hand side by $(d,e)\leqslant\|d\|_{\star}\,\|e\|\leqslant\tfrac{1}{2}\|d\|_{\star}^{2}+\tfrac{1}{2}\|e\|^{2}$, with the dual norm $\|\cdot\|_{\star}$, and integrating from time [math] to $t$ results in the error bound

[TABLE]

On the other hand, testing with $v=\partial_{t}e$ yields

[TABLE]

which leads similarly to the error bound

[TABLE]

This procedure is all-familiar, but it is not obvious how to extend it to time discretizations beyond the backward Euler and Crank–Nicolson methods. The use of energy estimates for BDF methods relies on the following remarkable results.

Lemma 10.1**.**

(Dahlquist [18]; see also [8] and [27, Section V.6])*

Let $\delta(\zeta)=\delta_{k}\zeta^{k}+\dotsb+\delta_{0}$ and $\mu(\zeta)=\mu_{k}\zeta^{k}+\dotsb+\mu_{0}$ be polynomials of degree at most $k\ ($and at least one of them of degree $k)$ that have no common divisor. Let $(\cdot,\cdot)$ be an inner product with associated norm $|\cdot|.$ If*

[TABLE]

then there exists a positive definite symmetric matrix $G=(g_{ij})\in{\mathbb{R}}^{k\times k}$ such that for $v_{0},\dotsc,v_{k}$ in the real inner product space,

[TABLE]

In combination with the preceding result for the multiplier $\mu(\zeta)=1-\eta_{k}\zeta,$ the following property of BDF methods up to order $5$ becomes important.

Lemma 10.2**.**

(Nevanlinna & Odeh [34])* For $k\leqslant 5,$ there exists $0\leqslant\eta_{k}<1$ such that for $\delta(\zeta)=\sum_{\ell=1}^{k}\frac{1}{\ell}(1-\zeta)^{\ell}$,*

[TABLE]

The smallest possible values of $\eta_{k}$ are

[TABLE]

Precise expressions for the optimal multipliers for the BDF methods of orders $3,4$ and $5$ are given by Akrivis & Katsoprinakis [1].

An immediate consequence of Lemma 10.2 and Lemma 10.1 is the relation

[TABLE]

with a positive definite symmetric matrix $G=(g_{ij})\in{\mathbb{R}}^{k\times k}$; it is this inequality that plays a crucial role in our energy estimates, and the same inequality for the inner product $a(\cdot,\cdot)$.

The error equation for the BDF time discretization of the linear parabolic problem (8.7a) reads

[TABLE]

with starting errors $e^{0},\dotsc,e^{k-1}$. When we test with $v=e^{n}-\eta_{k}e^{n-1}$, the first term can be estimated from below by (8.7b), the second term is bounded from below by $(1-\tfrac{1}{2}\eta_{k})\|e^{n}\|^{2}-\tfrac{1}{2}\eta_{k}\|e^{n-1}\|^{2}$, and the right-hand term is estimated from above by the Cauchy-Schwarz inequality. Summing up from $k$ to $n$ then yields the error bound

[TABLE]

where $C_{k}$ depends only on the order $k$ of the method. This kind of estimate for the BDF error has recently been used for a variety of linear and nonlinear parabolic problems [33, 3, 2, 30].

On the other hand, when we first subtract $\eta_{k}$ times the error equation for $n-1$ from the error equation with $n$ and then test with $\dot{e}^{n}$, we obtain

[TABLE]

Here, the second term is bounded from below by (8.7b) with the $a(\cdot,\cdot)$ inner product, the first term is bounded from below by $(1-\tfrac{1}{2}\eta_{k})|\dot{e}^{n}|^{2}-\tfrac{1}{2}\eta_{k}|\dot{e}^{n-1}|^{2}$, and the right-hand term is estimated from above by the Cauchy–Schwarz inequality. Summing up from $k$ to $n$ then yields the error bound

[TABLE]

It is this type of estimate that we use in the present paper for the nonlinear problem considered here. It has previously been used in [29].

Acknowledgment

The work of Michael Feischl, Balázs Kovács and Christian Lubich is supported by Deutsche Forschungsgemeinschaft – Project-ID 258734477 – SFB 1173.