Micromagnetics is a continuum theory describing magnetization patterns inside ferromagnetic media. The dynamics of magnetization is governed by the Landau-Lifshitz (LL) equation [35]. This highly nonlinear equation indicates a nonconvex 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 [34,37,43]. The solvability of LL-type equations can be found in [26,40,44]; two structures of the solution regularity have been investigated. In the framework of weak solution, the existence of global weak solution in R 3 was proved in [6]a n d in [27]o n a bounded domain ⊂ R 2 ; the nonuniqueness of weak solutions was demonstrated in [6]a s 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 R 3 was shown in [12]. Local existence and uniqueness of strong solutions on a bounded domain was proved in [11]; global existence and uniqueness of strong solutions for small-energy initial data on a 2-D bounded domain was established, provided that ∇m 0 H 1 ( ) is small enough for the bounded domain ⊂ R 2 . A similar uniqueness analysis was provided in [38]a s well.
We may refer to [27,52]f o r 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 [16,34]. The first finite element work was introduced by Alouges and his collaborators [2,3,5], 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 [4,33]. In another finite element work by Bartels and Prohl [7], 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 step-size condition k = O(h 2 ) is needed to guarantee the existence of the solution for the fixed point iteration (with k the temporal step-size and h the spatial mesh-size), which is highly restrictive. In fact, a theoretical justification of the unique solvability of the numerical solution is expected to come from the convergence of the fixed-point iteration by applying Banach fixed-point theorem, under the CFL-like condition k = O(h 2 ). A similar finite element scheme was reported by Cimrák [17]. Again, a nonlinear solver is necessary at each time step, and the same step-size condition has to be imposed. Also see the related finite element works [1,21], as well as the corresponding analysis. Meanwhile, the existing works of finite difference method to the LL equation may be referred to [20,23,28,31,51]. In [20], 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 secondorder accuracy in space. In [28], an updated source term was used, and an iteration algorithm was repeatedly performed until the numerical solution converges. In [31], 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 [48], in which |∇m| 2 was treated as the Lagrange multiplier for the non-convex constraint |m| = 1i n the pointwise sense with m 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 [30], 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 a CFL constraint, k = O(h 2 ). Such an approach has been applied in [32] (which yields the same numerical solution as the mid-point method, with second-order accuracy in time), and in a more general setting in [36]u s i n g the Cayley transform to lift the LL equation to the Lie algebra of the 3D rotation group. In addition, the first, second and fourth-order accurate temporal approximations were examined in [36], 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 [9,19], 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 [41]. Also see other related works [18,22,29,33], etc. In fact, there have been many existing numerical works, in which a time step constraint k = O(h 2 ) has to be imposed, such as explicit numerical schemes [3,30]( s o that the constraint k = O(h 2 ) has played a role of CFL condition), the fully implicit schemes [7,23,40]( s o that a nonlinear solver is needed), and a fixed point iteration approach [18], etc.
Based on the linearity of the discrete system, we can also classify numerical methods into the explicit scheme [3,30], the fully implicit scheme [7,23,40] and the semi-implicit scheme [15,20,24,36,48]. In particular, the semi-discrete schemes are introduced in [40]f o r 2-D and in [15]f o r 3D 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 [18], 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 [23], in which the unique solvability was proved under a very restrictive condition, k ≤ Ch 2 . 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 pointwise 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 secondorder backward differentiation formula (BDF) approximation is applied to obtain an intermediate magnetization m, 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 stepsize in terms of spatial stepsize) 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 W
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 L 2 error estimate is insufficient to recover such a bound for the numerical solution. Instead, we have to perform the H 1 error estimate, and such a W 1,∞ h bound could be obtained at the next time step as a consequence of the H 1 estimate, via the help of the inverse inequality combined with a mild time stepsize condition k = O(h). Careful error estimates for both the original magnetization m and the intermediate magnetization m 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 Section 2, 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 Section 3, including both the 1D and 3D examples to confirm the theoretical analysis. Conclusions are drawn in Section 4.
The LL equation reads as
where = ∂ and ν is the unit outward normal vector along . Here m : ⊂ R d → S 2 represents the magnetization vector field with |m| = 1, ∀x ∈ , d = 1, 2, 3i s the spatial dimension, and α > 0i s the damping parameter. The first term on the right hand side of (2.1)i s the gyromagnetic term, and the second term is the damping term. Compared to the original LL equation [35], (2.1)o n l y includes the exchange term which poses the main difficulty in numerical analysis, as done in the literature [7,18,20,24]. Application of the scheme (2.5)t o the original LL equation under external fields will be presented in another publication [50]. To ease the presentation, we set =[0, 1] d , in which d is the dimension.
The finite difference method is used to approximate (2.1) and (2.2). Denote the spatial step-size by h in the 1D case and divide [0, 1] into N x equal segments; see the schematic mesh in Fig. 1. Define
In the 3D case, we have spatial stepsizes
. The extrapolation formula along the z direction near z = 0 and z = 1i s m i, j,1 = m i, j,0 , m i, j,N z +1 = m i, j,N z .
(2.3) Extrapolation formulas for the boundary condition along other directions can be derived similarly.
In addition, given m e (•, t = 0) as the exact initial data at t = 0( w i t h m e the exact solution), the numerical initial data for m is set as m 0 i, j,k = P h m e (x i , ŷ j , ẑk , t = 0), P h is the point-wise interpolation.
(2.4)
The standard second-order centered difference applied to m results in
and the discrete gradient operator ∇ h m with m = (u, v, w) T reads as
Denote the temporal stepsize by k, and define t n = nk, n ≤ T k with T 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 t n+2 , 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:
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 [50]. 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, namely Algorithm 2.1:
(2.8)
The discrete boundary condition (2.3)i s imposed for mn+2 h in (2.7). In fact, this boundary condition could be rewritten as
Remark 2.1. To kick start the iteration of our method, we can use the first-order semi-implicit projection scheme using the first-order BDF and the first-order one-sided interpolation and the method is still second-order accurate.
Remark 2.2. We use the sparse LU factorization solver and also the Generalized Minimum Residual Method to solve the linear system in (2.7).
For simplicity of presentation, we assume that
An extension to the general case is straightforward. In the finite difference approximation, all the numerical values are assigned on the numerical grid points. As a result, the discrete grid functions (with notations f h , g h ), which are only defined over the corresponding numerical grid points, are introduced.
First, we introduce the discrete 2 inner product and discrete • 2 norm.
where d is the index set and I is the index which closely depends on d. In turn, the discrete • 2 norm is given by
(2.10)
Definition 2.3. For the grid function f h , we define the average of summation as
Definition 2.4. For any grid function f h with f h = 0, a discrete inverse Laplacian operator is defined as:
It is noticed that the zero-average constraint, ψ h = 0, makes the operator (h ) -1 uniquely defined. In turn, a discrete H -1
The first theoretical result is the unique solvability analysis of scheme (2.6)-(2.8). We observe that the unique solvability for (2.7) could be simplified as the analysis for
with p h , mh given. Theorem 2.1. Given p h , mh , the numerical scheme (2.11) is uniquely solvable.
The second theoretical result is the optimal rate convergence analysis.
(2.12)
in which the constant C > 0 is independent of kand h.
The proof of the standard inverse inequality and discrete Gronwall inequality could be obtained in existing textbooks; we just cite the results here. The inverse inequality presented in [14]i s in the finite element version; its extension to the finite difference version is straightforward. Lemma 2.1. (Inverse inequality) [14]. The inverse inequality implies that
, in which constant γ depends on the form of the discrete • 2 norm. Under the definition (2.9) and (2.10) for the cell-centered grid function, such a constant could be taken as γ = 1. [25]. Let {α j } j≥0 , {β j } j≥0 and {ω j } j≥0 be sequences of real numbers such that α j ≤ α j+1 ,β j ≥ 0, and ω j ≤ α j + j-1 i=0
Then it holds that
Proof. For simplicity of presentation, we only focus on the 1D case; an extension to the 2D and 3D formulas will be straightforward. A careful calculation reveals that
This identity is exactly the summation by parts formula (2.13). The proof of Lemma 2.3 is complete.
The following estimate will be utilized in the convergence analysis. In the sequel, for simplicity of our notation, we will use the uniform constant C to denote all the controllable constants in this paper.
)
.
(2.16)
Proof. Without loss of generality, we only look at the 1D case; an extension to the 3D case is straightforward. We begin with the following expansion
In turn, an application of the discrete Hölder inequality to (2.17)y i e l d s ( 2.14). Also note that
and
.
The following estimate will be used in the error estimate at the projection step.
Lemma 2.5. Consider m h = P h m e + h 2 m (1) , in which P h stands for the point-wise interpolation of a continuous function over the numerical grid points, the continuous function m e satisfies a regularity requirement m e W 1,∞ ≤ C, |m e | = 1 at a point-wise level, and the grid function m (1) satisfies m (1) ∞ + ∇ h m (1)
and we denote the numerical error functions as e h = m hm h , ẽh = mhm h . Then the following estimate is valid
(2.20)
(2.21)
(2.22)
For the last term on the right hand side of (2.21), we observe that
which in turn yields
(2.24)
As a result, a substitution of (2.22) and (2.24)i n t o( 2.21)l e a d s to the first estimate in (2.20).
For the second inequality, we notice that
The analysis for the first part is straightforward:
(2.26)
For the second part, we rewrite it as
based on the fact |m e | ≡ 1. For the nonlinear coefficient terms, we observe the following • ∞ bounds
, since m h = P h m e + h 2 m (1) , |m e |≡1, m (1)
)
As a result of these inequalities, we arrive at
2
in which C has a different value in the last step, which stands for another controllable constant. In turn, the following two estimates could be derived for the two expansion terms:
and
Therefore, we obtain
(2.27)
Finally, a substitution of (2.26) and (2.27)i n t o( 2.25)y i e l d s the second inequality in (2.20). This completes the proof of Lemma 2.5.
To facilitate the unique solvability analysis for (2.11), we denote q h =h mh . Note that q h = 0, due to the Neumann boundary condition for mh . Meanwhile, we observe that mh = (h ) -1 q h in general, since mh = 0. Instead, mh could be represented as follows:
and mh given by (2.6). (2.11)i s then rewritten as
(2.28)
Lemma 2.6 (Browder-Minty lemma [10,39]). Let X be a real, reflexive Banach space and let T : X → X (the dual space of X) be bounded, continuous, coercive (i.e.,
(T (u),u)
→+∞, as u X →+∞) and monotone. Then for any g ∈ X there exists a solution u ∈ Xo f the equation T (u) = g. Furthermore, if the operator Ti s strictly monotone, then the solution uis unique.
Then we proceed into the proof of Theorem 2.1.
Proof. Recall that (2.11)i s equivalent to (2.28). For any q 1,h , q 2,h with q 1,h = q 2,h = 0, we denote qh = q 1,hq 2,h and derive the following monotonicity estimate:
Note that the following equality and inequality have been applied in the second step:
mh × qh , qh =0, mh × ( mh × qh ), qh ≤0.
The third step is based on the fact that both C * q 1,h and C * q 2,h are constants, and q 1,h = q 2,h = 0, so that C * q 1,h -C * q 2,h , qh = 0. Moreover, for any q 1,h , q 2,h with q 1,h = q 2,h = 0, we get
and the equality only holds when q 1,h = q 2,h .
Therefore, an application of Lemma 2.6 implies a unique solution of both (2.28)a n d (2.11), which completes the proof of Theorem 2.1.
Proof. First, we construct an approximate solution m: m = m e + h 2 m (1) , (2.29) in which the auxiliary field m (1) satisfies the following Poisson equation
∂ z m (1) | z=0 =-
with boundary conditions along x and y directions defined in a similar way.
The purpose of such a construction will be illustrated later. Then we extend the approximate profile m to the numerical "ghost" points, according to the extrapolation formula (2.3): m i, j,0 = m i, j,1 , m i, j,N z +1 = m i, j,N z , (2.31) and the extrapolation for other boundaries can be formulated in the same manner. Subsequently, we prove that such an extrapolation yields a higher order O(h 5 ) approximation, instead of the standard O(h 3 ) accuracy. Also see the related works [42,45,46]i n the existing literature.
Performing a careful Taylor expansion for the exact solution around the boundary section z = 0, combined with the mesh point values: ẑ0 =-1 2 h, ẑ1 = 1 2 h, we get m e (x i , ŷ j , ẑ0 ) = m e (x i , ŷ j , ẑ1 ) -h∂ z m e (x i , ŷ j , 0) -
= m e (x i , ŷ j , ẑ1 ) -
in which the homogenous boundary condition has been applied in the second step. A similar Taylor expansion for the constructed profile m (1) reveals that m (1) (x i , ŷ j , ẑ0 ) = m (1) (x i , ŷ j , ẑ1 ) -h∂ z m (1)
In other words, the extrapolation formula (2.31)i s indeed O(h 5 ) accurate.
As a result of the boundary extrapolation estimate (2.34), we see that the discrete Laplacian of m yields the second-order accuracy, even at the mesh points around the boundary sections:
Moreover, a detailed calculation of Taylor expansion, in both time and space, leads to the following truncation error estimate:
In addition, a higher order Taylor expansion in space and time reveals the following estimate for the discrete gradient of the truncation error:
(2.37)
In fact, such a discrete • H 1 h bound for the truncation comes from the regularity assumption for the exact solution,
), as stated in Theorem 2.2, as well as the fact that m (1)
In turn, we introduce the numerical error functions ẽn
h , at a point-wise level. In other words, instead of a direct comparison between the numerical solution and the exact solution, we analyze the error function between the numerical solution and the constructed solution m h , due to its higher order consistency estimate (2.34)a r o u n d the boundary. A subtraction of (2.7)-(2.8)f r o m the consistency estimate (2.36)l e a d s to the error function evolution system:
Before we proceed into the formal error estimate, we establish the bound for the constructed approximate solution m and the numerical solution m h . For the approximate profile m ∈ L ∞ ([0, T ], C 5 ), which turns out to be the exact solution and an O(h 2 ) correction term, we still use C to denote its bound:
in which m h = P h m, the point-wise interpolation of the constructed continuous function m. In addition, we make the following a priori assumption for the numerical error function:
, for k = , + 1.
(2.40)
Such an assumption will be recovered by the convergence analysis at time step t +2 . In turn, an application of triangle inequality yields the desired W 1,∞ h bound for the numerical solutions m h and mh :
,
• Estimate of Ĩ1 : A combination of the summation by parts formula (2.13)( n o t i c e that the numerical error function ẽ satisfies the homogeneous Neumann boundary condition (2.3)) and inequality (2.14)r e s u l t s in
).
in which the bound for h m +2 h ∞ is given by the preliminary estimate (2.39), with r = 2.
• Estimate of the truncation error term Ĩ3 : An application of Cauchy inequality gives
(2.46)
• Estimate of Ĩ4 : It follows from (2.15)i n Lemma 2.4 that
).
• Estimates of Ĩ5 and Ĩ6 :
).
h ∞ is given by the preliminary estimate (2.39), with r = 2.
Meanwhile, the inner product of the left hand side of (2.38)w i t h ẽ +2 h turns out to be
Its combination with eqs. (2.44)t o (2.49) and (2.43)l e a d s to
(2.50)
However, the standard L 2 error estimate (2.50)d o e s not allow one to apply discrete Gronwall inequality, due to the H 1 h norms of the error function involved on the right hand side. To overcome this difficulty, we take a discrete inner product with the numerical error equation (2.38)b yh ẽ +2 h and see that
• Estimate of I 1 :
(2.52)
• Estimate of I 2 :
.
Similarly, the bound for ∇ h h m +2 h ∞ comes from the preliminary estimate (2.39), by taking r = 3. • Estimate of the truncation error term I 3 : (2.54) in which the discrete H 1 h estimate (2.37)f o r the local truncation error has been recalled.
• Estimate of I 4 : It follows from Lemma 2.16 in Lemma 2.4 that
• Estimates of I 5 and I 6 :
.
.
Again, the bound for ∇ h h m +2 h ∞ comes from the preliminary estimate (2.39), by taking r = 3.
And also, the inner product on the left hand side becomes
(2.58)
Substituting
(2.59)
As a consequence, a combination of (2.50) and (2.59)y i e l d s
(2.60)
+ Ck(k 4 + h 4 ).
At this point, recalling the W
Its substitution into (2.60)l e a d s to
+ Ck(k 4 + h 4 ).
In turn, an application of discrete Gronwall inequality (in Lemma 2.2) yields the desired convergence estimate for ẽh :
). An application of Lemma 2.1, as well as the time step constraint k ≤ Ch, leads to
, so that the second part of the a priori assumption (2.40)h a s been recovered at time step k = n. In turn, the W
Similar to the derivation of (2.61), we also get
, so that the first part of the a priori assumption (2.40) has been recovered at time step k = + 2. This completes the proof of Theorem 2.2.
Remark 2.3. The regularity assumption for the exact solution, namely m e ∈ C 3 ([0, T ]; In other words, the optimal rate error estimate (2.12) stands for a local-in-time theoretical result. In addition, since the finite difference numerical method is evaluated at the collocation grid points, instead of the ones based on a weak formulation, it usually requires higher order regularity requirement for the exact solution in the optimal rate convergence estimate than that of the finite element approach; see the related finite difference analysis for various gradient flows [8,47,49], etc.
In this section, we perform 1D and 3D numerical experiments for the final time T = 1t o verify the theoretical analysis in Section 2. Rate of convergence is obtained via the least-squares fitting for a sequence of error data recorded with successive stepsize refinements.
In details, we test four examples: 1D example with a forcing term and the given exact solution, 1D example without the exact solution, 3D example with a forcing term and the given exact solution, and 3D example with the full Landau-Lifshitz equation for simulating domain wall dynamics. Solutions in these four cases satisfy the homogenous Neumann boundary condition (2.2). In the presence of an forcing term, the LL equation reads as
and m e the exact solution. In more details, the forcing term f is evaluated at t n+2 in the numerical scheme (2.7). Only one linear system of equations needs to solve at each time step. In all examples, we find that the scheme is unconditionally stable.
. The exact solution is given by m e = (cos(x 2 (1x) 2 ) sin t, sin(x 2 (1x) 2 ) sin t, cos t) T , which satisfies the homogeneous Neumann boundary condition. A point-wise interpolation is applied to m e (•, t = 0) to obtain the initial data at the grid points, following formula (2.4). Results in Table 1 and Fig. 2 suggest the second-order accuracy in both time and space of the proposed numerical method in the discrete H 1 -norm; Results in Table 2 indicate the unconditional stability in the 1D case.
Example 3.2 (1D example without the exact solution). For this example, in the absence of the forcing term, we do not have the exact solution. For comparison, we first set h and k small enough to obtain a numerical solution which will be used as the reference solution. The initial condition is chosen as m 0 (x, 0) = (0, 0, 1) T for x ∈ . To get the temporal accuracy, we set h = 1D -4 and k = 1D -4t o get the reference solution and then record the temporal error with varying k in Table 3 and Fig. 3a. To get the spatial accuracy, we set h = 1/3 8 and k = 1D -4t o get the reference solution and record the error in Table 4 and Fig. 3b. Again, the second-order accuracy in both time and space in the discrete H 1 -norm have been confirmed.
The given exact solution is given by m e = (cos(XY Z) sin t, sin(XY Z) sin t, cos t) T , 2 . Again, a point-wise interpolation is applied to m e (•, t = 0) to obtain the initial data at the grid points, following formula (2.4).
Table 5, Table 6 and Fig. 4 shows the second-order convergence in both time and space in the 3D case. Results in Table 7 indicate the unconditional stability of the proposed numerical method in the 3D case. We visualize the magnetization in
where the effective field h eff includes not only the exchange field m, but also the anisotropy field, the external field h e and the demagnetization (stray) field h s . The uniaxial material with the x axis as the easy direction is considered, thus
where and Q are both dimensionless constants, e 2 = (0, 1, 0) and e 3 = (0, 0, 1) are unit vectors, h s takes the nondimensionalized form Consider a ferromagnetic nanostrip of size 0.8 × 0.1 × 0.004 µm 3 with 128 × 32 × 2g r i d points. In our simulations, the material parameters used for Permalloy, which is an alloy of Nickel (80%) and Iron (20%), are the exchange constant 1.3 × 10 -11 J/m, the anisotropy constant 1.0 × 10 2 J/m 3 , the saturation magnetization constant 8.0 × 10 5 A/m, the damping coefficient α = 0.1 and the temporal step-size k = 1 ps. A steady state is reached if the relative change in the total energy is less than 10 -7 . A transverse domain wall is formed with the in plane head-to-head Néel wall as the initial state (Fig. 6a), specified by m 0 h = (tanh(s), sech(s), 0) T , and6c. Details of the implementation can be found in [50].
Fig. 4. Accuracy of the proposed numerical method when α = 0.01 in the 3D case. (a) Temporal accuracy of our method on the uniform mesh when h x = h y = h z = 1/32 and α = 0.01; (b) Spatial accuracy of the proposed numerical method on the uniform mesh when k = 1/2048 and α = 0.01.
Unconditional stability of the proposed numerical method in the 3D case when α = 0.01.
In this paper, we have proposed and analyzed a second-order time stepping scheme to solve the LL equation. The second-order BDF is applied for temporal discretization and a linearized multistep approximation is used for the nonlinear coefficients on the right hand side of the equation. The resulting scheme avoids a well-known difficulty associated with the nonlinearity of the system, and its unique solvability is established via the monotonicity analysis of the system. In addition, an optimal rate convergence analysis is provided, by making use of a linearized stability analysis for the numerical error functions, in which the W 1,∞ h error estimate at the projection step has played an important role. Numerical experiments in both 1D and 3D cases are presented to verify the unconditional stability and the second-order accuracy in both space and time, and applied to the domain wall dynamics driven by an external field. The technique presented here may be applicable to the model for current-driven domain wall dynamics [13], which shall be explored as a future project.