Convergence Analysis of A Second-order Semi-implicit Projection Method for Landau-Lifshitz EquationSubmitted to the editors DATE:July 5, 2019.

Convergence Analysis of A Second-order Semi-implicit Projection Method for Landau-Lifshitz Equationthanks: Submitted to the editors DATE:July 5, 2019.

Jingrun Chen School of Mathematical Sciences, Soochow University, Suzhou, China. (E-mail addresses: jingrunchen@suda.edu.cn, 20164207002@stu.suda.edu.cn)    Cheng Wang Mathematics Department; University of Massachusetts; North Dartmouth, MA 02747, USA. (E-mail address: cwang1@umassd.edu)    Changjian Xie22footnotemark: 2
Abstract

The numerical approximation for the Landau-Lifshitz equation, the dynamics of magnetization in a ferromagnetic material, is taken into consideration. This highly nonlinear equation, with a non-convex constraint, has several equivalent forms, and involves solving an auxiliary problem in the infinite domain. All these features have posed interesting challenges in developing numerical methods. In this paper, we first present a fully discrete semi-implicit method for solving the Landau-Lifshitz equation based on the second-order backward differentiation formula and the one-sided extrapolation (using previous time-step numerical values). A projection step is further used to preserve the length of the magnetization. Subsequently, we provide a rigorous convergence analysis for the fully discrete numerical solution, with second-order accuracy in both time and space, provided that the spatial step-size is the same order as the temporal step-size. And also, the unique solvability of the numerical solution is theoretically justified, which turns out to be the first such result for the micromagnetics model. All these theoretical properties are verified by numerical examples in both one- and three- dimensional spaces.

Key words. Landau-Lifshitz equation, backward differentiation formula, semi-implicit scheme, second-order accuracy

AMS subject classifications. 35K61, 65N06, 65N12

1 Introduction

Micromagnetics is a continuum theory describing magnetization patterns inside ferromagnetic media. The dynamics of magnetization is governed by the Landau-Lifshitz (LL) equation [30]. This highly nonlinear equation indicates a non-convex constraint, which has always been a well-known difficulty in the numerical analysis. And also, this equation has several equivalent forms, and an auxiliary problem in the infinite domain has to be involved. All these features have posed interesting challenges in developing numerical methods. In the past several decades, many works have focused on the mathematical theory and numerical analysis of the LL equation [29, 32, 38]. The well-posedness of LL-type equations can be found in [21, 35, 39]; two structures of the solution regularity have been investigated. In the framework of weak solution, the existence of global weak solution in was proved in [5] and in [22] on a bounded domain ; the nonuniqueness of weak solutions was demonstrated in [5] as well. In the framework of strong solution, local existence and uniqueness, and global existence and uniqueness with small-energy initial data for strong solutions to the LL equation in was shown in [10]. Local existence and uniqueness of strong solutions on a bounded domain was proved in [9]; global existence and uniqueness of strong solutions for small-energy initial data on a 2-D bounded domain was established, provided that is small enough for the bounded domain . A similar uniqueness analysis was provided in [33] as well. We may refer to [22, 45] for the existence of unique local strong solution.

Accordingly, numerous numerical approaches have been proposed to demonstrate the mathematical theory; review articles could be found in [13, 29]. The first finite element work was introduced by Alouges and his collaborators [1, 2, 3, 4], in which rigorous convergence proof was included with first-order accuracy in time and second-order accuracy in space. This method was further developed to reach almost the second-order temporal accuracy [3, 28]. In another finite element work by Bartels and Prohl [6], they presented an implicit time integration method with second-order accuracy and unconditional stability. However, a nonlinear solver is needed at each time step, and a theoretical justification of the unique solvability of the numerical solution has not been available. And also, a step-size condition is needed to guarantee the existence of the solution for the fixed point iteration (with the temporal step-size and the spatial mesh-size), which is highly restrictive. A similar finite element scheme was reported by Cimrák [14]. Again, a nonlinear solver is necessary at each time step, and the same step-size condition has to be imposed. The existing works of finite difference method to the LL equation may be referred to [18, 19, 23, 26, 44]. In [18], a time stepping method in the form of a projection method was proposed; this method is implicit and unconditionally stable, and the rigorous proof was provided with the first-order accuracy in time and second-order accuracy in space. In [23], an updated source term was used, and an iteration algorithm was repeatedly performed until the numerical solution converges. In [26], the explicit and implicit mimetic finite difference algorithm was developed.

Regarding to the temporal discretization, the first kind of time-stepping scheme is the Gauss-Seidel projection method proposed by Wang, García-Cervera, and E [42], in which was treated as the Lagrange multiplier for the non-convex constraint in the point-wise sense with the magnetization vector field. The resulting method is first-order accurate in time and is unconditionally stable. The second kind of time-stepping scheme is called geometric integration method. In [25], Jiang, Kaper, and Leaf developed the semi-analytic integration method by analytically integrating the system of ODEs, obtained after a spatial discretization of the LL equation. This is an explicit method with first-order accuracy, hence is subject to the CFL constraint. Such an approach has been applied in [27] (which yields the same numerical solution as the mid-point method, with second-order accuracy in time), and in a more general setting in [31] using the Cayley transform to lift the LL equation to the Lie algebra of the three-dimensional rotation group. In addition, the first, second and fourth-order accurate temporal approximations were examined in [31], which is more amenable for building numerical schemes with the high-order accuracy. The third kind of time-stepping scheme is called the mid-point method [7, 16], which is second-order accurate, unconditionally stable, and preserves the Lyapunov and Hamiltonian structures of the LL equation. Moreover, the fourth kind of time-stepping method is the high-order Runge-Kutta algorithms [36]. Also see other related works [15, 17, 24, 28], etc.

Based on the linearity of the discrete system, we can also classify numerical methods into the explicit scheme [2, 25], the fully implicit scheme [6, 19, 35] and the semi-implicit scheme [12, 18, 20, 31, 42]. In particular, the semi-discrete schemes are introduced in [35] for 2-D and in [12] for 3-D formulation of the LL equation. Error estimates are derived under the existence assumption for the strong solution.

From the perspective of convergence analysis, it is worthy of mentioning [15], in which the fixed point iteration technique was used for handling the nonlinearities; the second-order convergence in time was proved, and was confirmed by numerical examples. It is noticed that, for all above-mentioned works with the established convergence analysis, a nonlinear solver has to be used at each time step, for the sake of numerical stability. However, the unique solvability analysis for these nonlinear numerical schemes has been a very challenging issue at the theoretical level, due to the highly complicated form in the nonlinear term. The only relevant analysis was reported in [19], in which the unique solvability was proved under a very restrictive condition, . And also, a projection step has been used in many existing works, to preserve the length of the magnetization. Its nonlinear nature makes a theoretical analysis highly non-trivial. In turn, a derivation of the following numerical scheme is greatly desired: second-order accuracy in time and linearity of the scheme at each time step, so that the length of magnetization is preserved in the point-wise sense, and an optimal rate error estimate and unconditionally unique solvability analysis could be established at a theoretical level.

In this work, we propose and analyze a second-order accurate scheme that satisfies these desired properties. The second-order backward differentiation formula (BDF) approximation is applied to obtain an intermediate magnetization , and the right-hand-side nonlinear terms are treated in a semi-implicit style with a second-order extrapolation applied to the explicit coefficients. Such a numerical algorithm leads to a linear system of equations with variable coefficients to solve at each time step. Its unconditionally unique solvability (no condition is needed for the temporal step-size in terms of spatial step-size) is guaranteed by a careful application of the monotonicity analysis, the so-called Browder-Minty lemma. A projection step is further used to preserve the unit length of magnetization at each time step, which poses a non-convex constraint. More importantly, we provide a rigorous convergence and error estimate, by the usage of the linearized stability analysis for the numerical error functions. In particular, we notice that, an a priori bound assumption for the numerical solution at the previous time steps has to be imposed to pass through the convergence analysis. As a consequence, the standard error estimate is insufficient to recover such a bound for the numerical solution. Instead, we have to perform the error estimate, and such a bound could be obtained at the next time step as a consequence of the estimate, via the help of the inverse inequality combined with a mild time step-size condition . Careful error estimates for both the original magnetization and the intermediate magnetization have to be taken into consideration at the projection step (a highly nonlinear operation). To the best of our knowledge, it is the first such result to report an optimal convergence analysis with second order accuracy in both time and space.

The rest of this paper is organized as follows. In LABEL:sec:main_theory, we introduce the fully discrete numerical scheme and state the main theoretical results: unique solvability analysis and optimal rate convergence analysis. Detailed proofs are also provided in this section. Numerical results are presented in LABEL:sec:experiments, including both the 1-D and 3-D examples to confirm the theoretical analysis. Conclusions are drawn in LABEL:sec:conclusions.

2 Main theoretical results

The LL equation reads as

(2.1)

with

(2.2)

where and is the unit outward normal vector along . Here represents the magnetization vector field with , is the spatial dimension, and is the damping parameter. The first term on the right hand side of LABEL:c1 is the gyromagnetic term, and the second term is the damping term. Compared to the original LL equation [30], LABEL:c1 only includes the exchange term which poses the main difficulty in numerical analysis, as done in the literature [18, 15, 6, 20]. Application of the scheme LABEL:eqn:bdf2 to the original LL equation under external fields will be presented in another publication [43]. To ease the presentation, we set when and when .

2.1 Finite difference discretization and the fully discrete scheme

The finite difference method is used to approximate LABEL:c1 and LABEL:boundary. Denote the spatial step-szie by in the 1-D case and divide into equal segments; see the schematic mesh in LABEL:spatial_mesh. Define , , with , . and , . Denote the magnetization obtained by the numerical scheme at by . To approximate the boundary condition LABEL:boundary, we introduce ghost points and apply Taylor expansions for , at , and , at , respectively. We then obtain a third order extrapolation formula:

In the 3-D case, we have spatial step-sizes , , and grid points , with , and ( , , ). The extrapolation formula along the direction near and is

(2.3)

Extrapolation formulas for the boundary condition along other directions can be derived similarly.

ghost point

ghost point

Figure 1: Illustration of the 1-D spatial mesh.

The standard second-order centered difference applied to results in

and the discrete gradient operator with reads as

Denote the temporal step-size by , and define , with the final time. The second-order BDF approximation is applied to the temporal derivative:

Note that the right hand side of the above equation is evaluated at , a direct application of the BDF method leads to a fully nonlinear scheme. To overcome this difficulty, we come up with a semi-implicit scheme, in which the nonlinear coefficient is approximated by the second-order extrapolation formula:

(2.4)

A projection step is then added to preserve the length of magnetization. This scheme has been used to study domain wall dynamics under external magnetic fields [43]. However, this scheme is difficult to conduct the convergence analysis due to the lack of numerical stability of Lax-Richtmyer type. To overcome this difficulty, we separate the time-marching step and the projection step in the following way:

(2.5)
(2.6)
(2.7)

2.2 Some notations and a few preliminary estimates

For simplicity of presentation, we assume that so that . An extension to the general case is straightforward.

First, we introduce the discrete inner product and discrete norm.

Definition 2.1 (Inner product and norm).

For grid functions and over the uniform numerical grid, we define

where is the index set and is the index which closely depends on . In turn, the discrete norm is given by

In addition, the discrete -norm is given by .

Definition 2.2 (Discrete norm).

For the grid function over the uniform numerical grid, we define

Definition 2.3.

For the grid function , we define the average of summation as

Definition 2.4.

For the grid function , we define the discrete -norm as

The proof of inverse inequality, discrete Gronwall inequality, and summation by parts formula could be obtained in many existing textbooks; we just cite the results here.

Lemma 2.1.

(Inverse inequality). The classical inverse inequality implies that

Lemma 2.2.

(Discrete Gronwall inequality). Let , and be sequences of real numbers such that

Then it holds that

Lemma 2.3 (Summation by parts).

For any grid functions and , with satisfying the discrete boundary condition LABEL:BC-1, the following identity is valid:

(2.8)

The following estimate will be utilized in the convergence analysis. In the sequel, for simplicity of our notation, we will use the uniform constant to denote all the controllable constants in this paper.

Lemma 2.4 (Discrete gradient acting on cross product).

For grid functions and over the uniform numerical grid, we have

(2.9)
(2.10)
(2.11)

Proof.

Without loss of generality, we only look at the 1-D case; an extension to the 3-D case is straightforward. We begin with the following expansion

(2.12)

In turn, an application of the discrete Hölder inequality to LABEL:expansion-1 yields  LABEL:lem27_1. Also note that

and

    

The following estimate will be used in the error estimate at the projection step.

Lemma 2.5.

Consider with the exact solution to LABEL:c1 and at a point-wise level, and . For any numerical solution , we define . Suppose both numerical profiles satisfy the following bounds

(2.13)
(2.14)

and we denote the numerical error functions as , . Then the following estimate is valid

(2.15)

Proof.

A direct calculation shows that

(2.16)

Since , we get

(2.17)

For the last term on the right hand side of LABEL:lem_6-3, we observe that

(2.18)
(2.19)

As a result, a substitution of LABEL:lem_6-4 and LABEL:lem_6-5-2 into LABEL:lem_6-3 leads to the first estimate in LABEL:lem_6-2.

For the second inequality, we notice that

(2.20)

The analysis for the first part is straightforward:

(2.21)

For the second part, we rewrite it as

based on the fact . In turn, the following two bounds could be derived:

and

Therefore, we obtain

(2.22)

Finally, a substitution of LABEL:lem_6-7 and LABEL:lem_6-8-4 into LABEL:lem_6-6 yields the second inequality in LABEL:lem_6-2. This completes the proof of LABEL:lem_6-0.     

2.3 The main theoretical results

The first theoretical result is the unique solvability analysis of scheme LABEL:scheme-1-1-LABEL:scheme-1-2. We observe that the unique solvability for LABEL:scheme-1-1 could be simplified as the analysis for

(2.23)

with , given.

Theorem 2.1.

Given , , the numerical scheme LABEL:scheme-alt-1 is uniquely solvable.

To facilitate the unique solvability analysis for LABEL:scheme-alt-1, we denote . Note that , due to the Neumann boundary condition for . Meanwhile, we observe that in general, since . Instead, could be represented as follows:

and given by LABEL:cc. LABEL:scheme-alt-1 is then rewritten as

(2.24)
Lemma 2.2 (Browder-Minty lemma [8, 34]).

Let X be a real, reflexive Banach space and let (the dual space of ) be bounded, continuous, coercive (i.e., , as ) and monotone. Then for any there exists a solution of the equation .

Furthermore, if the operator is strictly monotone, then the solution is unique.

Then we proceed into the proof of LABEL:thm:solvability.

Proof.

Recall that LABEL:scheme-alt-1 is equivalent to LABEL:scheme-alt-2. For any , with , we denote and derive the following monotonicity estimate:

Note that the following equality and inequality have been applied in the second step:

The third step is based on the fact that both and are constants, and