Convergence of a mass-lumped finite element method for the Landau-Lifshitz equation
aa r X i v : . [ m a t h . NA ] J u l CONVERGENCE OF A MASS-LUMPED FINITE ELEMENTMETHOD FOR THE LANDAU-LIFSHITZ EQUATION
EUGENIA KIM AND JON WILKENING
Abstract.
The dynamics of the magnetic distribution in a ferromagnetic ma-terial is governed by the Landau-Lifshitz equation, which is a nonlinear geo-metric dispersive equation with a nonconvex constraint that requires the mag-netization to remain of unit length throughout the domain. In this article, wepresent a mass-lumped finite element method for the Landau-Lifshitz equation.This method preserves the nonconvex constraint at each node of the finite ele-ment mesh, and is energy nonincreasing. We show that the numerical solutionof our method for the Landau-Lifshitz equation converges to a weak solutionof the Landau-Lifshitz-Gilbert equation using a simple proof technique thatcancels out the product of weakly convergent sequences. Numerical tests forboth explicit and implicit versions of the method on a unit square with periodicboundary conditions are provided for structured and unstructured meshes.
Micromagnetics is the study of the behavior of ferromagnetic materials at sub-micron length scales, including magnetization reversal and hysteresis effects [22].The dynamics of the magnetic distribution of a ferromagnetic material occupying aregion Ω ⊂ R or R is governed by the Landau-Lifshitz (LL) equation [22, 25, 36].The magnetization m ( x, t ) : Ω × [0 , T ] → S ⊂ R satisfies(0.1) ∂ t m = − m × h − αm × ( m × h ) in Ω ∂m∂ν = 0 on ∂ Ω m ( x,
0) = m ( x )where α is a dimensionless damping parameter and h is an effective field given by(0.2) h ( m ) := − δ E δm ( m ) = η ∆ m − Q ( m e + m e ) + h s ( m ) + h e . Here η is the exchange constant, Q is an anisotropy constant, h s is the stray field, h e is an external field, and δ E δm is the functional derivative of the Landau-Lifshitzenergy, defined by(0.3) E ( m ) = η Z Ω |∇ m | + Q Z Ω m + m − Z Ω h s ( m ) · m − Z Ω h e · m. Mathematics Subject Classification.
Primary 65M60 35Q60; Secondary 78M10.
Key words and phrases.
Landau-Lifshitz equation, Landau-Lifshitz-Gilbert equation, micro-magnetics, finite element methods, mass-lumped method, convergence, weak solutions.The first author was supported in part by the U.S. Department of Energy, Office of Science,Office of Workforce Development for Teachers and Scientists, Office of Science Graduate StudentResearch (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institutefor Science and Education for the DOE under contract number DE-AC05-06OR23100.The second author was supported in part by the US Department of Energy, Office of Science,Applied Scientific Computing Research, under award number DE-AC02-05CH11231.
The first term is the exchange energy, which tries to align the magnetization locally;the second term is the anisotropy energy, which tries to orient the magnetization incertain easy direction taken to be e ; the third term is the stray field energy, whichis induced by the magnetization distribution inside the material; and the last termis the external field energy, which tries to align the magnetization with an externalfield. We denote the lower order terms in (0.2) by(0.4) ¯ h ( m ) := − Q ( m e + m e ) + h s ( m ) + h e . When considering mathematical properties such as existence and regularity of thesolution, these terms can be considered lower order compared to the exchange term[6]. They also have fewer derivatives than the exchange term, and can be treatedas lower order when developing numerical methods.The stray field h s depends on m via h s = −∇ φ , where the potential φ satisfies(0.5) ∆ φ = ( ∇ · m in Ω0 on ∂ Ω[ φ ] ∂ Ω = 0 , (cid:20) ∂φ∂ν (cid:21) ∂ Ω = − m · ν. Here [ v ] ∂ Ω ( x ) = v ( x + ) − v ( x − ) is the jump in v across the boundary ∂ Ω from inside( − ) to outside (+); see [25].There are several equivalent forms of the Landau-Lifshitz (LL) equation. Thefollowing is the Landau-Lifshitz-Gilbert (LLG) equation :(0.6) ∂ t m − αm × ∂ t m = − (1 + α )( m × h ) . Also, the equation(0.7) α∂ t m + m × ∂ t m = (1 + α )( h − ( h · m ) m )is equivalent to LL and LLG; see [6]. If only the gyromagnetic term is present inequation (0.1), i.e. if ∂ t m = − m × ∆ m , it is called a Schr¨odinger map into S [29].This is a geometric generalization of the linear Schr¨odinger equation. If only thedamping term is present, i.e. if ∂ t m = − m × ( m × ∆ m ), it is called a harmonic mapheat flow into S [29].In 1935, Landau and Lifshitz [37] calculated the structure of the domain wallsbetween antiparallel domains, which started the theory of micromagnetics. Thetheory was further developed by W. F. Brown Jr in [15]. Applications of micro-magnetics include magnetic sensor technology, magnetic recording, and magneticstorage devices such as hard drives and magnetic memory (MRAM) [22].Finite difference methods for micromagnetics can be derived in two different ways[41]. The first is a field-based approach in which the effective field h is discretizeddirectly, and the other is an energy-based approach in which the effective field isderived from the discretized energy. Finite element methods can also be derivedin a number of ways. In [44], the Landau-Lifshitz-Gilbert equation (0.6) is used toobtain a discrete system by approximating the magnetization by piecewise linearfunction on a finite element mesh and then applying time integration in the resultingsystem of ODEs. In [22], the effective field h is calculated by taking a functionalderivative of the discretized energy, where the magnetization in the energy formulais approximated by piecewise linear functions. Extensive work has also been donedeveloping time stepping schemes for micromagnetics [18,25]. In [31], semi-analyticintegration in time was introduced, which is explicit and first order in time and ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 3 allows stepsize control. In [34, 38], geometric integration methods were appliedto the Landau-Lifshitz equation. In [26, 47], a Gauss-Seidel projection methodwas developed that treats the gyromagnetic term and damping term separately toovercome the difficulty of the stiffness and the nonlinearity of the equation.However, relatively little work has been done deriving error estimates or es-tablishing rigorous convergence results for weak solutions. In a series of papers[4, 6, 8], Alouges and various co-authors introduced a convergent finite elementmethod based on equation (0.7), an equivalent form of the Landau-Lifshitz equa-tion, that is first order in time. The idea is to use a tangent plane formulation ateach timestep, where the velocity vector lies in the finite element space perpendicu-lar to the magnetization at each node. One advantage of this method is that at eachstep, only a linear system has to be solved, although the Landau-Lifshitz-Gilbertequation is nonlinear. More recently, they developed a formally second order intime scheme [7, 35] that performs better than first order in practice, though notfully at second order. Another finite element scheme was introduced by Bartelsand Prohl in [11] based on the Landau-Lifshitz-Gilbert equation, which is an im-plicit, unconditionally stable method, but involves solving nonlinear equations ateach timestep. This method is second order in time; however, there is still a timestep constraint, namely that k/h remain bounded, to guarantee the existence of thesolutions of the nonlinear systems. Cimr´ak [19] introduced a finite element methodbased on the Landau-Lifshitz equation, which is an implicit, unconditionally stablemethod, but also has nonlinear inner iterations. We note that the Backward Eulermethod and higher-order diagonally implicit Runge-Kutta (DIRK) methods [30]generally involve solving nonlinear equations at each internal Runge-Kutta stagewhen applied to nonlinear PDEs.In this article, we introduce a family of mass-lumped finite element methodsfor the Landau-Lifshitz equation. The implicit version is similar in computationalcomplexity to the algorithms in [4,6,8] in that each timestep involves solving a largesparse linear system. The explicit method is more efficient than the explicit versionof [4, 6, 8] as it is completely explicit — it does not even require that a linear systeminvolving a mass matrix be solved as the effective mass matrix is diagonal. Themethod involves finding the velocity vector in the tangent plane of the magnetiza-tion by discretizing the Landau-Lifshitz equation instead of the Landau-Lifshitz-Gilbert equation, as was done in [4, 6–8, 35]. By building a numerical scheme basedon the Landau-Lifshitz equation instead of the Landau-Lifshitz-Gilbert equation,we can naturally apply the scheme to limiting cases such as the Schr¨odinger mapor harmonic map heat flow [12, 20, 29, 32, 33]. The main result of the paper is aproof that the numerical solution of our scheme for the Landau-Lifshitz equationconverges to a weak solution of the Landau-Lifshitz-Gilbert equation, using a sim-ple technique that cancels out the product of two weakly convergent sequences.Our proof builds on tools developed in [8]. For simplicity, we defer the treat-ment of the stray field to future work. This term poses computational challenges[1, 2, 13, 16, 21, 23, 27, 39, 40, 43, 46, 48], but does not affect the convergence resultssince it is a lower order term in comparison to the exchange term; see [8].The paper is organized as follows. In section 1, we introduce a finite elementmesh and review the weak formulation of the Landau-Lifshitz-Gilbert equation. Insection 2, the main algorithm and the main theorem will be introduced. In section3, we conduct a numerical test for the equation h = ∆ m on the unit square with EUGENIA KIM AND JON WILKENING periodic boundary conditions, where an exact analytical solution is known from[24]. In section 4, the proof of the main theorem will be presented.1.
Weak solutions, meshes and the finite element space
Let us denote Ω T = Ω × (0 , T ). The definition of a weak solution of the Landau-Lifshitz-Gilbert equation is given by Definition 1.1.
Let m ( x ) ∈ H (Ω) with | m ( x ) | = 1 a.e. Then m is a weaksolution of (0.6) if for all T > m ( x, t ) ∈ H (Ω T ) , | m ( x, t ) | = 1 a.e. ,(ii) m ( x,
0) = m ( x ) in the trace sense,(iii) m satisfies Z Ω T ∂ t m · φ − α Z Ω T ( m × ∂ t m ) · φ (1.1) = (1 + α ) η d X l =1 Z Ω T ( m × ∂ x l m ) · ∂ x l φ − (1 + α ) d X l =1 Z Ω T ( m × ¯ h ( m )) · φ. for all φ ∈ H (Ω T ) .(iv) m satisfies an energy inequality(1.2) C Z Ω T | ∂ t m | + E ( m ( x, T )) ≤ E ( m ( x, . for some constant C >
0, where the energy E ( m ) is defined in equation(0.3).In [6, 11], the value of C in ( iv ) is taken to be C = α α . The existence of globalweak solution of the Landau-Lifshitz equation in Ω ⊂ R into S was proved in[9, 28]. The nonuniqueness of weak solution was proved in [9].Let the domain Ω ⊂ R d where d = 2 or 3 be discretized into triangular ortetrahedral elements {T h } h of mesh size at most h with vertices ( x i ) Ni =1 . Let thefamily of partitions T = {T h } h be admissible, shape regular and uniform [14]. Let { φ i } ≤ i ≤ N be piecewise linear nodal basis functions for T , such that φ i ( x j ) = δ ij ,where δ ij is a Kronecker delta function. We will consider a vector-valued finiteelement space F h defined by F h = { w h | w h ( x ) = P Ni =1 w hi φ i ( x ) , w hi ∈ R } . The discrete magnetization m h is required to belong to the submanifold M h of F h defined by M h = { m h ∈ F h | m h ( x ) = P Ni =1 m hi φ i ( x ) , | m hi | = 1 } . We define theinterpolation operator I h : C (Ω , R ) → F h by(1.3) I h ( m ) = N X i =1 m ( x i ) φ i ( x ) . ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 5
We need the following additional conditions for our finite element method: Thereexist some constants C , C , C , C > C h d ≤ b i = Z Ω φ i ≤ C h d , | M ij | = (cid:12)(cid:12)(cid:12)(cid:12)Z Ω φ i φ j (cid:12)(cid:12)(cid:12)(cid:12) ≤ C h d , | ∂ x l φ i | ≤ C h Z Ω ∇ φ i · ∇ φ j ≤ , for i = j, for all h > i, j = 1 , . . . , N and l = 1 , . . . , d .2. The finite element scheme, the algorithm, and the main theorem
To illustrate how we obtain Algorithm 1 below, consider the simple case inwhich only the exchange energy term is present in the effective field, i.e. h = η △ m from (0.2). Let’s first consider the weak form of the Landau-Lifshitz equation with h = η △ m ,(2.1) Z Ω T ∂ t m · w = η d X l =1 Z Ω T ( m × ∂ x l m ) · ∂ x l w − αη d X l =1 Z Ω T ∂ x l m · ∂ x l w + αη d X l =1 Z Ω T ( ∂ x l m · ∂ x l m )( m · w ) . Taking this weak form as a hint, we would like to find v = P Nj =1 v j φ j ∈ F h suchthat(2.2) Z Ω N X j =1 v j φ j · w i φ i = η d X l =1 N X j =1 Z Ω ( m i × m j ∂ x l φ j ) · ∂ x l φ i w i − αη d X l =1 N X j =1 Z Ω m j ∂ x l φ j · ∂ x l φ i w i + αη d X l =1 N X j =1 Z Ω ( ∂ x l φ j m j · m i )( m i · w i ∂ x l φ i )for i = 1 , . . . , N , where m = P Nj =1 m j φ j ( x ) ∈ M h , w ∈ ( C ∞ (Ω)) and w i = I h ( w )( x i ) = w ( x i ). Then, with w i as (1 , , , ,
0) or (0 , ,
1) in equation (2.2),we obtain(2.3) ( M v ) i = η m i × ( A m ) i + αη m i × ( m i × ( A m ) i )for i = 1 , . . . , N , where M = M M
00 0 M and A = A A
00 0 A are 3 N × N block diagonal matrices with each block M and A a mass or stiffness matrix,i.e. M ij = R Ω φ i φ j , and A ij = P dl =1 R Ω ∂ x l φ i ∂ x l φ j . Note that m i · ( M v ) i = 0, soapproximating v by ˆ v = M vb yields a tangent vector to the constraint manifold M h ,where b i = R Ω φ i . The left hand side of (2.3) is then b i ˆ v i which is a mass-lumpingapproximation. This suggests the following algorithm. Algorithm 1.
For a given time ¯
T >
0, set J = [ ¯ Tk ]. EUGENIA KIM AND JON WILKENING (1) Set an initial discrete magnetization m at the nodes of the finite elementmesh described in § j = 0 , . . . , J ,a. compute a velocity vector ˆ v ji at each node by(2.4) ˆ v ji = ( M v j ) i b i = η m ji × ( A m + θk A ˆ v ) ji + αη m ji × ( m ji × ( A m + θk A ˆ v ) ji ) b i − m ji × ( M ¯ h ( m ) + θk M ¯ h (ˆ v )) ji + αm ji × ( m ji × ( M ¯ h ( m ) + θk M ¯ h (ˆ v )) ji ) b i . for θ ∈ [0 ,
1] and for i = 1 , . . . , N .b. Compute m j +1 i = m ji + k ˆ v ji | m ji + k ˆ v ji | for i = 1 , . . . , N .We define the time-interpolated magnetization and velocity as in [8]: Definition 2.1.
For ( x, t ) ∈ Ω × [ jk, ( j + 1) k ) ⊂ Ω × [0 , T ), where T = Jk , define m h,k ( x, t ) = m j ( x ) , ¯ m h,k ( x, t ) = t − jkk m j +1 ( x ) + ( j + 1) k − tk m j ( x ) , ˆ v h,k ( x, t ) = ˆ v j ( x ) ,v h,k ( x, t ) = v j ( x ) . The main theorem in this article is the following theorem, which is proved insection 4.
Theorem 2.2.
Let m ∈ H (Ω , S ) and suppose m h → m in H (Ω) as h → .Let θ ∈ [0 , , and for ≤ θ < , assume that kh ≤ C , for some C > . Ifthe triangulation T = {T h } h satisfies condition (1.4), then the sequence { m h,k } ,constructed by Algorithm 1 and definition 2.1, has a subsequence that convergesweakly to a weak solution of the Landau-Lifshitz equation. Numerical Results
Before proving Theorem 2.2, we demonstrate the effectiveness of the schemeon a test problem. We conduct a numerical experiment for the Landau-Lifshitzequation (0.1) with effective field involving only the exchange energy term, with h = ∆ m in equation (0.2), on the unit square with periodic boundary conditions.This corresponds to setting η = 1 and ¯ h = 0 in equation (2.4) in Algorithm 1. Forthe convergence study, we used an explicit method ( θ = 0) and an implicit method( θ = 0 .
5) on a structured and unstructured mesh. The unstructured mesh withpoint and line sources, which is an arbitrary mesh, was generated using DistMesh[42], with an example shown in figure 1.The L ∞ and L errors were measured relative to an exact solution for theLandau-Lifshitz equation with h = ∆ m from [24], namely(3.1) m x ( x , x , t ) = 1 d ( t ) sin β cos( k ( x + x ) + g ( t )) ,m y ( x , x , t ) = 1 d ( t ) sin β sin( k ( x + x ) + g ( t )) ,m z ( x , x , t ) = 1 d ( t ) e k αt cos β. ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 7
Figure 1.
Unstructured mesh with point and line sources, with h = 1 / β = π , k = 2 π , d ( t ) = p sin β + e k αt cos β and g ( t ) = α log( d ( t )+ e k αt cos β β ).The numerical results are summarized in the tables 1 and 2. Figure 2 shows theconvergence rate of the methods, which is first order in the time step k and secondorder in the mesh size h . -3 -2 -1 spatial step h10 -6 -5 -4 -3 -2 -1 E rr o r At time 0.001 1 2structured - L errorstructured - L errorunstructured - L errorunstructured - L error 10 -3 -2 -1 spatial step h10 -6 -5 -4 -3 -2 -1 E rr o r At time 0.001 1 2structured - L errorstructured - L errorunstructured - L errorunstructured - L error Figure 2.
Convergence plot, Left : Explicit method, Right: Im-plicit method.3.1.
Going beyond first order in time.
In this section, we propose a methodwhich is second order in time, by replacing the nonlinear projection step 2 (b)in Algorithm 1 by a linear projection step, and test the convergence order. InAlgorithm 1, step 2 (a) can be viewed as the predictor step and 2 (b) as the correctorstep. The corrector step was used to conserve the length of the magnetization ateach node. By replacing this nonlinear projection by a linear projection step, it notonly preserves the length of the magnetization, but also makes the method higher
EUGENIA KIM AND JON WILKENING
Structured mesh Unstructured mesh h || m − m h || L ∞ rate || m − m h || L rate || m − m h || L ∞ rate || m − m h || L rate32 8.22e-05 2.00 7.40e-04 2.00 5.61e-03 1.28 3.83e-03 1.6564 2.06e-05 2.00 1.85e-04 2.00 2.32e-03 1.56 1.22e-03 1.97128 5.15e-06 2.00 4.63e-05 2.00 7.87e-04 1.81 3.13e-04 2.01256 1.29e-06 1.16e-05 2.25e-04 7.77e-05 Table 1.
Explicit method ( θ = 0 ) : L ∞ and L error and con-vergence rates on a structured and unstructured mesh with spatialstep h , time step k = 8 · − h and time 0 . h || m − m h || L ∞ rate || m − m h || L rate || m − m h || L ∞ rate || m − m h || L rate32 8.26e-05 2.00 7.40e-04 2.00 5.61e-03 1.28 3.83e-03 1.6564 2.07e-05 2.00 1.85e-04 2.00 2.32e-03 1.56 1.22e-03 1.97128 5.17e-06 2.00 4.63e-05 2.00 7.87e-04 1.81 3.13e-04 2.01256 1.29e-06 1.16e-05 2.25e-04 7.77e-05 Table 2.
Implicit method ( θ = ) : L ∞ and L error andconvergence rates on a structured and unstructured mesh , withspatial step h , time step k = 0 . h and time 0 . × Algorithm 2.
For a given time ¯
T >
0, set J = [ ¯ Tk ].(1) Set an initial discrete magnetization m at the nodes of the finite elementmesh described in § j = 0 , . . . , J ,a. compute an intermediate magnetization vector m ∗ i at each node by m ∗ i − m ji k = ˆ v ji = ( M v j ) i b i = η m ji × ( A m + θk A ˆ v ) ji + αη m ji × ( m ji × ( A m + θk A ˆ v ) i ) j b i − m ji × ( M ¯ h ( m ) + θk M ¯ h (ˆ v )) ji + αm ji × ( m ji × ( M ¯ h ( m ) + θk M ¯ h (ˆ v )) i ) j b i . for θ ∈ [0 ,
1] and for i = 1 , . . . , N . ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 9 b. Compute m j +1 i for i = 1 , . . . , N . m j +1 i − m ji k == η m j +1 i + m ji × ( A m j +1 / i ) + αη m j +1 i + m ji × ( m j +1 / i × ( A m j +1 / ) b i − m j +1 i + m ji × M ¯ h ( m j +1 / ) i + α m j +1 i + m ji × ( m j +1 / i × ( M ¯ h ( m j +1 / i )) b i . where m j +1 / = m j + m ∗ for i = 1 , . . . , N .As before, we conduct a numerical test for the Landau-Lifshitz equation (0.1)with effective field involving only the exchange energy term, with h = ∆ m inequation (0.2), on the unit square with periodic boundary conditions, to comparethe two algorithms. For the convergence study, we used an implicit method ( θ =0 .
5) on a structured and unstructured mesh. One of the unstructured meshes wasshown in figure 1. The L errors were measured relative to an analytical solution(3.1) for the Landau-Lifshitz equation with h = ∆ m . The numerical results aresummarized in the tables 3. Figure 3 shows the convergence rates of the methods,which shows first order in k for Algorithm 1 and second order convergence forAlgorithm 2. -2 spatial step h10 -6 -5 -4 -3 E rr o r At time 0.01 1 21Algorithm 1Algorithm 2 10 -2 spatial step h10 -5 -4 -3 -2 E rr o r At time 0.01 1 21Algorithm 1Algorithm 2
Figure 3.
Convergence plot, Left : Structured mesh, Right : Un-structured mesh. 4.
Proof of Theorem 2.2
In this section, we present the proof of the theorem, which states that the se-quence { m h,k } , constructed by Algorithm 1 and Definition 2.1, has a subsequencethat converges weakly to a weak solution m of the Landau-Lifshitz-Gilbert equationunder some conditions. That is, we show that the limit m satisfies Definition 1.1.In section 4.1, we derive a discretization of the weak form of the Landau-Lifshitz-Gilbert equation satisfied by the { m h,k } , namely (4.3). In section 4.2, we deriveenergy estimates to show that the sequences m h,k , ¯ m h,k and ˆ v h,k converge to m in acertain sense made precise in section 4.3. In section 4.4, we show that each term ofthe discretization of the weak form converges to the appropriate limit, so that the Structured mesh Unstructured mesh h Alg. 1 rate Alg. 2 rate Alg. 1 rate Alg. 2 rate32 5.14e-04 1.32 1.87e-04 2.01 2.10e-03 1.61 1.91e-03 1.7264 2.06e-04 1.18 4.66e-05 2.00 6.87e-04 1.75 5.81e-04 2.01128 9.12e-05 1.10 1.16e-05 2.00 2.04e-04 1.57 1.44e-04 2.03256 4.27e-05 2.91e-06 6.84e-05 3.53e-05
Table 3.
Implicit method ( θ = ) : L error and convergencerates on a structured and unstructured mesh , with spatial step h ,time step k = 0 . h and time 0 . m satisfies the weak form of the Landau-Lifshitz-Gilbert equation. In section4.5, we show that the limit m satisfies the energy inequality (1.2) in Definition 1.1(iv). Finally, in section 4.6, we establish that the magnitude of m is 1 a.e. in Ω T .4.1. Equations that m h,k and v h,k satisfy. In this section, we derive a discretiza-tion of the weak form of the Landau-Lifshitz-Gilbert equation. This form is easierto use for the proof of Theorem 2.2, since it does not involve the product of theweakly convergent sequences. In general, a product of weakly convergent sequencesis not weakly convergent. It is convergent only in some certain cases, such as whenthe sequences satisfy the hypothesis of the div-curl lemma [17, 45].The generalized version of equation (2.2) including all the terms in the effectivefield h (0.2) and with 0 ≤ θ ≤ Z Ω v h,k · w h = η X l,i Z Ω ( m h,ki × ∂ x l ( m h,k + θk ˆ v h,k )) · ∂ x l φ i w hi − αη X l,i Z Ω ∂ x l ( m h,k + θk ˆ v h,k ) · ∂ x l φ i w hi + αη X l,i Z Ω ( ∂ x l ( m h,k + θk ˆ v h,k ) · m h,ki )( m h,ki · w hi ) ∂ x l φ i − X i Z Ω ( m h,ki × ¯ h ( m h,k + θk ˆ v h,k )) · φ i w hi + α X i Z Ω ¯ h ( m h,k + θk ˆ v h,k ) · φ i w hi − α X i Z Ω (¯ h ( m h,k + θk ˆ v h,k ) · m i )( m i · w hi φ i ) . ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 11
In fact, by taking w hi as (1 , , , (0 , ,
0) and (0 , ,
1) in (4.1), we get (2.4) in Algo-rithm 1. Setting w h = P Nj =1 ( m h,kj × u hj ) φ j in (4.1), we have(4.2) − X i Z Ω ( m h,ki × v h,k ) · u hi φ i = η X l Z Ω ∂ x l ( m h,k + θk ˆ v h,k ) · ∂ x l u h − η X l,i Z Ω ( ∂ x l ( m h,k + θk ˆ v h,k ) · m h,ki ) ( m h,ki · u hi ) ∂ x l φ i + αη X l,i Z Ω ( m h,ki × ∂ x l ( m h,k + θk ˆ v h,k )) · ∂ x l φ i u hi − X i Z Ω (¯ h ( m h,k + θk ˆ v h,k )) · φ i u hi + X i Z Ω (¯ h ( m h,k + θk ˆ v h,k ) · m h,ki ) ( m h,ki · u hi ) φ i − α X i Z Ω ( m h,ki × ¯ h ( m h,k + θk ˆ v h,k )) · u hi φ i . Equations (4.1) and (4.2) have terms that contain the product of weakly convergentsequences, namely the third term of the right hand side of (4.1), and the second termof the right hand side of (4.2), αη P l,i R Ω ( ∂ x l ( m h,k + θk ˆ v h,k ) · m h,ki )( m h,ki · w hi ) ∂ x l φ i . By adding α times equation (4.2) to equation (4.1), we eliminate the terms thatcontain the product of weakly convergent sequences :(4.3) Z Ω (cid:2) v h,k · w h − α X i ( m h,ki × v h,k ) · w hi φ i (cid:3) = (1 + α ) (cid:20) η X l,i Z Ω (cid:2) ( m h,ki × ∂ x l m h,k ) · ∂ x l φ i w hi + θk ( m h,ki × ∂ x l ˆ v h,k ) · ∂ x l φ i w hi (cid:3) − X i Z Ω (cid:2) ( m h,ki × ¯ h ( m h,k )) · w hi φ i + θk ( m h,ki × ¯ h (ˆ v h,k )) · w hi φ i (cid:3)(cid:21) . This is a similar procedure to subtracting α times the following equation(4.4) m × ∂ t m = − m × ( m × h ) + αm × h from the Landau-Lifshitz equation (0.1) to get the Landau-Lifshitz-Gilbert equation(0.6). Here, equation (4.4) is obtained by taking m × the Landau-Lifshitz equation(0.1).4.2. Energy inequality.
In this section, we derive the energy inequalities we willneed to prove Theorem 2.2, namely (4.17) for 0 ≤ θ < and (4.18) for ≤ θ ≤ Theorem 4.1.
For the P approximation in Ω ⊂ R , if (4.5) Z Ω ∇ φ i · ∇ φ j ≤ , for i = j, then for all w = P Ni =1 w i φ i ∈ F h such that | w i | ≥ for i = 1 , . . . , N , we have (4.6) Z Ω (cid:12)(cid:12)(cid:12)(cid:12) ∇ I h ( w | w | ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z Ω |∇ w | . In 3D, we have (4.6), if an additional condition that all dihedral angles of thetetrahedra of the mesh are smaller than π is satisfied, along with (4.5). Also, wewill use inequality (14) of [8],(4.7) (cid:13)(cid:13) ¯ h ( m ) (cid:13)(cid:13) L ≤ C k m k L + C and equation (25) from [7],(4.8) k h s ( m ) k L ≤ C k m k L where C are positive constants, depending only on Ω. Furthermore, we will usean inequality (20) of [8] in the proof, which states there exists C > ≤ p < ∞ and all φ h ∈ F h , we have(4.9) 1 C k φ h k pL p ≤ h d N X i =1 | φ h ( x i ) | p ≤ C k φ h k pL p . Moreover, we will assume that there exists C > Z Ω |∇ v h | ≤ C h Z Ω | v h | for all v h ∈ F h .Taking w h = P Nj =1 ( m h,kj × u hj ) φ j in (4.3), and setting u h = ˆ v h,k , we have(4.11) − α X i Z Ω v h,ki · ˆ v i φ i = (1 + α ) (cid:20) η X l,i Z Ω (cid:2) ( ∂ x l m h,k · ∂ x l φ i ˆ v i ) + θk ( ∂ x l ˆ v h,k · ∂ x l φ i ˆ v i ) (cid:3) − X i Z Ω (cid:2) (¯ h ( m h,k ) · ˆ v i ) φ i + θk (¯ h (ˆ v h,k ) · ˆ v i ) φ i (cid:3)(cid:21) where we have used the fact m h,ki · ˆ v h,ki = 0 for i = 1 , . . . , N . This equation can bewritten as(4.12)( ∇ m, ∇ ˆ v ) = − θk k∇ ˆ v k L − α α η X i | ( M v ) j | b j + 1 η (¯ h ( m ) , ˆ v ) + θkη (¯ h (ˆ v ) , ˆ v ) . We now derive an energy estimate. We have(4.13)12 (cid:13)(cid:13) ∇ m j +1 (cid:13)(cid:13) L ≤ (cid:13)(cid:13) ∇ m j + k ∇ ˆ v j (cid:13)(cid:13) L = 12 (cid:13)(cid:13) ∇ m j (cid:13)(cid:13) L + k ( ∇ m j , ∇ ˆ v j ) + 12 k (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L ≤ (cid:13)(cid:13) ∇ m j (cid:13)(cid:13) L − k ( α α ) 1 η X i | ( M v ) ji | b i b i + 12 k (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L − θk (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L + kη (¯ h ( m j ) , ˆ v j ) + θ k η (¯ h (ˆ v j ) , ˆ v j ) ≤ (cid:13)(cid:13) ∇ m j (cid:13)(cid:13) L − k ( α α ) 1 η C C (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L − ( θ −
12 ) k (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L + kη (¯ h ( m j ) , ˆ v j )+ θ k η (¯ h e , ˆ v j ) ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 13 where the first inequality is obtained by Theorem 4.1, the second inequality byequation (4.12), and the last inequality by the fact ( h s (ˆ v j ) , ˆ v j ) <
0. We have theestimate for the last two terms of the above inequality :(4.14) (cid:12)(cid:12) (¯ h ( m j ) + θk ¯ h e , ˆ v j ) (cid:12)(cid:12) ≤ (cid:13)(cid:13) ¯ h ( m j ) + θk ¯ h e (cid:13)(cid:13) L (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L ≤ C (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L ≤ ǫ (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L + 14 ǫ C for some C >
0, where the second inequality is obtained by equation (4.7) and thelast inequality by Young’s inequality with ǫ = α α C C . Summing the inequality(4.13) from j = 0 , . . . , J − (cid:13)(cid:13) ∇ m J (cid:13)(cid:13) L + k ( 12 η ( α α ) C C − C ( 12 − θ ) kh ) J − X j =0 (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L ≤ (cid:13)(cid:13) ∇ m (cid:13)(cid:13) L + C T with kh ≤ C < α α C C C η , for 0 ≤ θ < , and(4.16)12 (cid:13)(cid:13) ∇ m J (cid:13)(cid:13) L + k ( 12 η α α ) C C J − X j =0 (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L + ( θ −
12 ) k J − X j =0 (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L ≤ (cid:13)(cid:13) ∇ m (cid:13)(cid:13) L + C T for ≤ θ ≤
1, and for some C > Z Ω |∇ m h,k ( · , T ) | + ( 12 η ( α α ) C C − C C ) Z Ω T | ˆ v h,k | ≤ Z Ω |∇ m h,k ( · , | + C T with C < α α C C C η , for 0 ≤ θ < and(4.18) 12 Z Ω (cid:12)(cid:12) ∇ m h,k ( · , T ) (cid:12)(cid:12) + ( 12 η α α ) C C Z Ω T (cid:12)(cid:12) ˆ v h,k (cid:12)(cid:12) + ( θ −
12 ) k Z Ω T (cid:12)(cid:12) ∇ ˆ v h,k (cid:12)(cid:12) ≤ Z Ω (cid:12)(cid:12) ∇ m h,k ( · , (cid:12)(cid:12) + C T. for ≤ θ ≤ Weak convergence of m h,k , ¯ m h,k and ˆ v h,k . In this section, we show theweak convergence of ¯ m h,k and ˆ v h,k and strong convergence of m h,k in some sense,based on the energy estimates (4.17) and (4.18). We follow similar arguments fromsection 6 of [7].Since we have(4.19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m j +1 i − m ji k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) ˆ v ji (cid:12)(cid:12)(cid:12) for i = 1 , . . . , N and j = 0 , . . . , J −
1, we have(4.20) (cid:13)(cid:13) ∂ t ¯ m h,k (cid:13)(cid:13) L (Ω) = (cid:13)(cid:13)(cid:13)(cid:13) m j +1 − m j k (cid:13)(cid:13)(cid:13)(cid:13) L (Ω) ≤ C (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω) . Thus, we have(4.21) (cid:13)(cid:13) ∂ t ¯ m h,k (cid:13)(cid:13) L (Ω T ) = (cid:13)(cid:13)(cid:13)(cid:13) m j +1 − m j k (cid:13)(cid:13)(cid:13)(cid:13) L (Ω T ) ≤ C (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω T )4 EUGENIA KIM AND JON WILKENING which is bounded by the energy inequalities, (4.17) for 0 ≤ θ < and (4.18) for ≤ θ ≤
1. Hence, ¯ m h,k is bounded in H (Ω T ) and ˆ v h,k is bounded in L (Ω T ) by(4.21) and by the energy inequalities, (4.17) for 0 ≤ θ < and (4.18) for ≤ θ ≤ m ∈ H (Ω T ) and ˆ v ∈ L (Ω T ) suchthat(4.22) ¯ m h,k → m weakly in H (Ω T ) , ¯ m h,k → m strongly in L (Ω T ) , ˆ v h,k → ˆ v weakly in L (Ω T ) . Moreover, we have(4.23) (cid:12)(cid:12)(cid:12) m j +1 i − m ji − k ˆ v ji (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m ji + k ˆ v ji | m ji + k ˆ v ji | − m ji − k ˆ v ji (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) − | m ji + k ˆ v ji | (cid:12)(cid:12)(cid:12) ≤ k (cid:12)(cid:12)(cid:12) ˆ v ji (cid:12)(cid:12)(cid:12) , since | m ji + k ˆ v ji | = q k | ˆ v ji | ≤ k | ˆ v ji | , for i = 1 , . . . , N and j = 0 , . . . , J − (cid:13)(cid:13) ∂ t ¯ m h,k − ˆ v h,k (cid:13)(cid:13) L (Ω T ) ≤ kC C (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω T ) which converges to 0 as h, k →
0, so(4.25) ∂ t m = ˆ v. Furthermore, since(4.26) (cid:13)(cid:13) m h,k − ¯ m h,k (cid:13)(cid:13) L (Ω T ) = (cid:13)(cid:13)(cid:13)(cid:13) ( t − jk ) m j +1 − m j k (cid:13)(cid:13)(cid:13)(cid:13) L (Ω T ) ≤ k (cid:13)(cid:13) ∂ t ¯ m h,k (cid:13)(cid:13) L (Ω T ) and the right hand side goes to 0 as h, k →
0, we have(4.27) m h,k → m strongly in L (Ω T ) . In summary, we have shown that there exist a subsequence of { ¯ m h,k } that con-verges weakly in H (Ω × (0 , T )), a subsequence of { ¯ v h,k } that converges weakly in L (Ω × (0 , T )), and a subsequence of { m h,k } converges strongly in L (Ω T ) based onthe energy estimates (4.17) and (4.18). However, in our numerical tests in section3, it was not necessary to take subsequences and the method was in fact secondorder in space and first order in time. Thus, there is still a gap in what we are ableto prove and the practical performance of the algorithm in cases where the weaksolution is unique and sufficiently smooth.4.4. The proof that the limit m actually satisfies Landau-Lifshitz-Gilbertequation. In this section, we show that each term of equation (4.3) converges tothe appropriate limit, so that the limit m of the sequences { ¯ m h,k } and { m h,k } satisfies the weak form of the Landau-Lifshitz-Gilbert equation (1.1) in Definition1.1. Lemma 4.2.
Let the sequences { m h,k } , { ¯ m h,k } , { ˆ v h,k } , and { v h,k } be definedby Definition 2.1. Also, let m ∈ H (Ω T ) be the limit as in (4.22) and (4.27). ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 15
Moreover, let’s assume w ∈ ( C ∞ (Ω T ) ∩ ( H (Ω T )) , and w h = I h ( w ) ∈ F h as inequation (1.3). Then we have (4.28) lim h,k → Z Ω T v h,k · w h = lim h,k → Z T N X j =1 ˆ v h,kj · w hj Z Ω φ j = Z Ω T ∂ t m · w. Proof.
The difference between the last two terms is bounded by(4.29) (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T I h (ˆ v h,k · w h ) − ˆ v h,k · w h (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T ˆ v h,k · w h − ∂ t m · w (cid:12)(cid:12)(cid:12)(cid:12) . The first term of (4.29) has the following estimate. For each element L , we haveˆ v h,k ( · , t ) · w h ( · , t ) ∈ C ∞ ( L ) and(4.30) (cid:13)(cid:13) I h (ˆ v h,k · w h ) − ˆ v h,k · w h (cid:13)(cid:13) L ( L ) ≤ C h (cid:13)(cid:13) ∆(ˆ v h,k · w h ) (cid:13)(cid:13) L ( L ) ≤ C h ( (cid:13)(cid:13) ∆ˆ v h,k · w h (cid:13)(cid:13) L ( L ) + (cid:13)(cid:13) ∇ ˆ v h,k · ∇ w h (cid:13)(cid:13) L ( L ) + (cid:13)(cid:13) ˆ v h,k · ∆ w h (cid:13)(cid:13) L ( L ) ) ≤ C h (( (cid:13)(cid:13) ∇ ˆ v h,k · ∇ w h (cid:13)(cid:13) L ( L ) )for some C >
0, where the first inequality is obtained by the Bramble-Hilbertlemma [14], and in the last inequality we have used ∆ˆ v h,k = 0 and ∆ w h = 0 in L ,since ˆ v h,k and w h are the sum of piecewise linear functions. We have the estimate(4.31) (cid:13)(cid:13) ∇ ˆ v h,k (cid:13)(cid:13) L (Ω) ≤ X L Z L | X i (ˆ v h,k ) i ∇ φ i | ≤ C h X L | X i ∈ I L (ˆ v h,k ) i | | L |≤ C h d − N X i =1 | (ˆ v h,k ) i | ≤ C h (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω) . for some constants C , C , C > I L is the index of nodes of L , where thesecond inequality is obtained by (4.10), and the last inequality by (4.9). Hence,(4.32) (cid:13)(cid:13) I h (ˆ v h,k · w h ) − ˆ v h,k · w h (cid:13)(cid:13) L (Ω T ) ≤ C h (cid:13)(cid:13) ∇ ˆ v h,k · ∇ w h (cid:13)(cid:13) L (Ω T ) ≤ C h (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω T ) for some constant C >
0. Therefore, the first term of (4.29) goes to 0 as h, k → v h,k to ∂ t m which are equations (4.22) and (4.25) . (cid:3) Lemma 4.3.
Under the same assumptions of Lemma 4.2, we have (4.33) lim h,k → Z Ω T X i ( m h,ki × v h,k ) · w hi φ i = lim h,k → Z Ω T X i ( m h,ki × ˆ v h,ki ) · w hi φ i = Z Ω T ( m × ∂ t m ) · w Proof.
The difference between the last two terms is bounded by(4.34) (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T I h (( m h,k ) a (ˆ v h,k ) b ( w h ) c ) − ( m h,k ) a (ˆ v h,k ) b ( w h ) c (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T ( m h,k ) a (ˆ v h,k ) b ( w h ) c − m a ( ∂ t m ) b w c (cid:12)(cid:12)(cid:12)(cid:12) for some a, b, c ∈ { , , } . The first term of (4.34), has the following estimate. Foreach element L , we have ( m h,k ) a (ˆ v h,k ) b ( w h ) c ∈ C ∞ ( L ) and(4.35) (cid:13)(cid:13) I h (( m h,k ) a (ˆ v h,k ) b ( w h ) c ) − ( m h,k ) a (ˆ v h,k ) b ( w h ) c (cid:13)(cid:13) L ( L ) ≤ C h ( (cid:13)(cid:13) ∆(( m h,k ) a (ˆ v h,k ) b ( w h ) c ) (cid:13)(cid:13) L ( L ) ) ≤ C h ( (cid:13)(cid:13) ∇ ( m h,k ) a ∇ (ˆ v h,k ) b ( w h ) c (cid:13)(cid:13) L ( L ) + (cid:13)(cid:13) ∇ ( m h,k ) a (ˆ v h,k ) b ∇ ( w h ) c (cid:13)(cid:13) L ( L ) + (cid:13)(cid:13) ( m h,k ) a ∇ (ˆ v h,k ) b ∇ ( w h ) c (cid:13)(cid:13) L ( L ) )for some constant C >
0, where the first inequality is obtained by Bramble-Hilbertlemma, and in the last inequality we have used ∆ ˆ m h,k = 0, ∆ˆ v h,k = 0 and ∆ w h = 0in L , since m h,k ˆ v h,k and w h are the sum of piecewise linear functions. Hence, wehave the estimate(4.36) (cid:13)(cid:13) I h (( m h,k ) a (ˆ v h,k ) b ( w h ) c ) − ( m h,k ) a (ˆ v h,k ) b ( w h ) c (cid:13)(cid:13) L (Ω T ) ≤ C h (cid:13)(cid:13) (ˆ v h,k ) b (cid:13)(cid:13) L (Ω T ) ( (cid:13)(cid:13) ∇ ( m h,k ) a (cid:13)(cid:13) L (Ω T ) + h (cid:13)(cid:13) ∇ ( m h,k ) a (cid:13)(cid:13) L (Ω T ) + (cid:13)(cid:13) ( m h,k ) a (cid:13)(cid:13) L (Ω T ) ) . for some constant C >
0, where we have used H¨older’s inequality for all theterms and used (4.31) for the first and the third terms. Therefore, the first term of(4.34) goes to 0 as h, k →
0. Moreover, the second term of (4.34) goes to 0 by theweak convergence of (ˆ v h,k ) b to ( ∂ t m ) b established in (4.22) and (4.25), and strongconvergence of ( m h,k ) a to m a . (cid:3) Lemma 4.4.
Under the same assumptions of Lemma 4.2, we have (4.37) lim h,k → X l,i Z Ω T ( m h,ki × ∂ x l m h,k ) · w hi ∂ x l φ i = X l Z Ω T ( m × ∂ x l m ) · ∂ x l w . Proof.
The difference between the last two terms is bounded by(4.38) (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T ( ∂ x l m h,k ) b ∂ x l I h (( m h,k ) c ( w h ) a ) − Z Ω T ( ∂ x l m h,k ) b (( m h,k ) c ( ∂ x l w h ) a ) (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T ( m h,k ) c ( ∂ x l m h,k ) b ( ∂ x l w h ) a − Z Ω T m c ( ∂ x l m ) b ( ∂ x l w ) a ) (cid:12)(cid:12)(cid:12)(cid:12) , for some a, b, c ∈ { , , } . The first term is bounded by(4.39) (cid:13)(cid:13) ( ∂ x l m h,k ) b (cid:13)(cid:13) L (Ω T ) (cid:13)(cid:13) ∂ x l I h (( m h,k ) c ( w h ) a ) − ∂ x l (( m h,k ) c ( w h ) a ) (cid:13)(cid:13) L (Ω T ) . For each element L , we have m h,k ( · , t ) w ( · , t ) ∈ C ∞ ( L ), and we have the estimate,(4.40) (cid:13)(cid:13) ∂ x l I h (( m h,k ) c ( w h ) a ) − ∂ x l (( m h,k ) c ( w h ) a ) (cid:13)(cid:13) L ( L ) ≤ C h | ( m h,k ) c ( w h ) a | H ( L ) for some constant C >
0, by the Bramble-Hilbert lemma. Moreover, we have theestimate,(4.41) | ( m h,k ) c ( w h ) a | H ( L ) = Z L | ∆(( m h,k ) c ( w h ) a ) | ≤ C Z L |∇ ( m h,k ) c | |∇ ( w h ) a | ≤ C (cid:13)(cid:13) ( m h,k ) c (cid:13)(cid:13) H ( L )ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 17 for some constants C , C >
0, since ∆ m h,k = 0 and ∆ w h = 0 in L , since m h,k and w h are the sum of piecewise linear functions. We get the estimate(4.42) (cid:13)(cid:13) ∂ x l I h (( m h,k ) c ( w h ) a ) − ∂ x l (( m h,k ) c ( w h ) a ) (cid:13)(cid:13) L (Ω T ) ≤ C C h (cid:13)(cid:13) ( m h,k ) c (cid:13)(cid:13) H (Ω T ) . Therefore, we may conclude that the first term of (4.38) goes to 0 as h, k → ∂ x l m h,k ) b to ( ∂ x l m ) b and strong convergence of ( m h,k ) c to m c , which gives (4.22) and (4.27). (cid:3) Lemma 4.5.
Under the same assumptions of Lemma 4.2, we have (4.43) lim h,k → (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k X i Z Ω T ( m h,ki × ∂ x l ˆ v h,k ) a ( ∂ x l w hi ) a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 . for ≤ θ ≤ .Proof. An upper bound for the sequence above is(4.44) √ k (cid:13)(cid:13)(cid:13) √ k ∂ x l (ˆ v h,k ) c (cid:13)(cid:13)(cid:13) L (Ω T ) (cid:13)(cid:13) ∇ ( I h ( m h,k ) b ( w h ) a ) (cid:13)(cid:13) L (Ω T ) . for some a, b, c ∈ { , , } . The term (cid:13)(cid:13)(cid:13) √ k ∂ x l (ˆ v h,k ) c (cid:13)(cid:13)(cid:13) L (Ω T ) in (4.44) is uniformlybounded, since (cid:13)(cid:13)(cid:13) √ k∂ x l (ˆ v h,k ) c (cid:13)(cid:13)(cid:13) L (Ω) ≤ C √ kh (cid:13)(cid:13) (ˆ v h,k ) c (cid:13)(cid:13) L (Ω) is uniformly boundedby (4.17) for 0 ≤ θ < , which is obtained by (4.10), and (cid:13)(cid:13)(cid:13) √ k∂ x l (ˆ v h,k ) c (cid:13)(cid:13)(cid:13) L (Ω) isuniformly bounded by equation (4.18) for ≤ θ ≤
1. For each element L , we have m h,k ( · , t ) w ( · , t ) ∈ C ∞ ( L ), so(4.45) (cid:13)(cid:13) ∇ I h (( m h,k ) b ( w h ) a ) − ∇ (( m h,k ) b ( w h ) a ) (cid:13)(cid:13) L ( L ) ≤ C h ( (cid:13)(cid:13) ( ∇ ( m h,k ) b ) (cid:13)(cid:13) L ( L ) ) , for some constant C >
0, by the Bramble-Hilbert lemma, and using ∆ m h,k = 0and ∆ w h = 0 in L , since m h,k and w h are the sum of piecewise linear functions.Thus, we have(4.46) (cid:13)(cid:13) ∇ ( I h ( m h,k ) b ( w h ) a ) (cid:13)(cid:13) L (Ω T ) ≤ (cid:13)(cid:13) ∇ ( m h,k ) b (cid:13)(cid:13) L (Ω T ) + C h ( (cid:13)(cid:13) ( ∇ ( m h,k ) b ) (cid:13)(cid:13) L (Ω T ) ) , which is uniformly bounded. Hence, (4.44) goes to 0 as h, k → (cid:3) Lemma 4.6.
Under the same assumptions of Lemma 4.2, we have (4.47) lim h,k → X i Z Ω T ( m h,ki × ¯ h ( m h,k )) · w hi φ i = Z Ω T ( m × ¯ h ( m )) · w. Proof.
An upper bound for the difference between the sequence and the limit isgiven by(4.48) (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T (¯ h ( m h,k )) a I h (( m h,k ) b ( w h ) c ) − Z Ω T (¯ h ( m h,k )) a ( m h,k ) b ( w h ) c ) (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)Z Ω T (¯ h ( m h,k )) a ( m h,k ) b ( w h ) c − Z Ω T (¯ h ( m )) a m b w c ) (cid:12)(cid:12)(cid:12)(cid:12) for some a, b, c ∈ { , , } . The first term of (4.48) is bounded by(4.49) (cid:13)(cid:13) ¯ h ( m h,k ) a (cid:13)(cid:13) L (Ω T ) (cid:13)(cid:13) I h (( m h,k ) b ( w h ) c ) − ( m h,k ) b ( w h ) c (cid:13)(cid:13) L (Ω T ) For each element L , we have m h,k ( · , t ) w ( · , t ) ∈ C ∞ ( L ), and we get the estimate,(4.50) (cid:13)(cid:13) I h (( m h,k ) b w c ) − (( m h,k ) b w c ) (cid:13)(cid:13) L ( L ) ≤ C h | ( m h,k ) b w c | H ( L ) for some constant C >
0, by the Bramble-Hilbert lemma. Moreover,(4.51) | ( m h,k ) b w c | H ( L ) ≤ C Z L |∇ ( m h,k ) b | |∇ ( w h ) c | + | ( m h,k ) b | | ∆( w h ) c | ≤ C (cid:13)(cid:13) ( m h,k ) b (cid:13)(cid:13) H ( L ) for some constant C >
0, and using the fact ∆ m h,k = ∆ w h = 0 in L , since m h,k and w h are the sum of piecewise linear functions. We get the estimate(4.52) (cid:13)(cid:13) I h (( m h,k ) b w c ) − (( m h,k ) b w c ) (cid:13)(cid:13) L (Ω T ) ≤ C h (cid:13)(cid:13) ( m h,k ) b (cid:13)(cid:13) H (Ω T ) . for some constant C >
0. Thus, the first term of (4.48) goes to 0 as h, k → h, k →
0, because of the strongconvergence of (¯ h ( m h,k )) a and ( m h,k ) b . (cid:3) Lemma 4.7.
Under the same assumptions of Lemma 4.2, we have (4.53) lim h,k → (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k X i Z Ω T ( m h,ki × ¯ h (ˆ v h,k )) · w hi φ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 . Proof.
An upper bound for the sequence above is(4.54) k (cid:13)(cid:13) ¯ h (ˆ v h,k )) (cid:13)(cid:13) L (Ω T ) (cid:13)(cid:13) w h (cid:13)(cid:13) L (Ω T ) . Since, (cid:13)(cid:13) ¯ h (ˆ v h,k )) (cid:13)(cid:13) L (Ω T ) ≤ ( C (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω T ) + C ) by (4.7), the term (cid:13)(cid:13) ¯ h (ˆ v h,k )) (cid:13)(cid:13) L (Ω T ) in (4.54) is uniformly bounded. Therefore, (4.54) goes to 0 as h, k → (cid:3) Energy of m . Recall the definition of the energy E ( m ) in (0.3). We followthe same arguments in section 6 of [7]. We have an energy estimate of m h,k as(4.55) E ( m j +1 ) − E ( m j ) ≤ − k ( α α ) C C || ˆ v j || L − ( θ −
12 ) k η ||∇ ˆ v j || L + k (¯ h ( m j ) , ˆ v j )+ θk (¯ h e , ˆ v j ) − Z Ω (¯ h ( m j +1 ) + ¯ h ( m j )) · ( m j +1 − m j ) . by (4.13) from section 4.2. For 0 ≤ θ < , the second term on the right has anupper bound(4.56)( θ −
12 ) k η (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L (Ω) ≤ k η (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L (Ω) ≤ C kη kh (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) ≤ C C ηk (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) and by choosing C ≤ α α C C C η , this term and the first term on the right handside of (4.55) can be combined to be less than equal to(4.57) − k α α ) C C || ˆ v j || L (Ω)ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 19 The second term on the right of equation(4.55) can be disregarded for ≤ θ ≤ (cid:12)(cid:12)(cid:12)(cid:12) k (¯ h ( m j ) , ˆ v j ) − Z Ω (¯ h ( m j +1 ) + ¯ h ( m j )) · ( m j +1 − m j ) (cid:12)(cid:12)(cid:12)(cid:12) . and has an upper bound(4.59) (cid:12)(cid:12)(cid:12)(cid:12)Z Ω ¯ h ( m j ) · ( m j +1 − m j − k ˆ v j ) (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) Z Ω (¯ h ( m j +1 ) − ¯ h ( m j )) · ( m j +1 − m j ) (cid:12)(cid:12)(cid:12)(cid:12) . The first term of (4.59) is bounded by(4.60) C k (cid:16)(cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) (cid:17) ≤ C k (cid:16)(cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) (cid:17) for some constant C >
0, by (4.23), and (4.9). The second term of (4.59) isbounded by C k (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) for some constant C >
0, by (4.19) and (4.9).The fourth term on the right has the upper bound | θk (¯ h e , ˆ v j ) | ≤ C k (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) for some constant C >
0. Then (4.55) has an upper bound(4.61) E ( m j +1 ) − E ( m j ) + k (cid:18) α α (cid:19) C C || ˆ v j || L (Ω) ≤ C k (cid:16)(cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) (cid:17) ≤ C k (cid:16)(cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) (cid:17) for some constants C , C >
0, by using Sobolev embedding theorem [3], (cid:13)(cid:13) ˆ v j (cid:13)(cid:13) L (Ω) ≤ C (cid:13)(cid:13) ∇ ˆ v j (cid:13)(cid:13) L (Ω) for some constant C >
0. Summing from j = 0 , . . . , J −
1, we get(4.62) E ( m J ) − E ( m ) + 12 (cid:18) α α (cid:19) C C Z Ω T | ˆ v h,k | ≤ C k (cid:16)(cid:13)(cid:13) ∇ ˆ v h,k (cid:13)(cid:13) L (Ω T ) + (cid:13)(cid:13) ˆ v h,k (cid:13)(cid:13) L (Ω T ) (cid:17) Therefore, taking h, k →
0, we get the energy inequality (1.2).4.6.
Magnitude of m . By the same argument in [8], we have | m ( x, t ) | = 1 a.e. for ( x, t ) ∈ Ω T ( See equation (28) and (29) on page 1347 of [8] ).5. Conclusion
We have presented a mass-lumped finite element method for the Landau-Lifshitzequation. We showed that the numerical solution of our method has a subsequencethat converges weakly to a weak solution of the Landau-Lifshitz-Gilbert equation.Numerical tests show that the method is second order accurate in space and firstorder accurate in time when the underlying solution is smooth. A second-order intime variant was also presented and tested numerically, but not analyzed rigorouslyin the present work.
References [1] Claas Abert, Lukas Exl, Gunnar Selke, Andr´e Drews, and Thomas Schrefl,
Numerical methodsfor the stray-field calculation: A comparison of recently developed algorithms , Journal ofMagnetism and Magnetic Materials (2013), 176–185. [2] Claas Abert, Gunnar Selke, Benjamin Kr¨uger, and Andr´e Drews,
A fast finite-differencemethod for micromagnetics using the magnetic scalar potential , IEEE Transactions on Mag-netics (2012), no. 3, 1105–1109.[3] Robert A Adams and John JF Fournier, Sobolev spaces , Vol. 140, Academic press, 2003.[4] Fran¸cois Alouges,
A new finite element scheme for Landau-Lifshitz equations , Discrete Con-tin. Dyn. Syst. Ser. S (2008), no. 2, 187–196.[5] Fran¸cois Alouges, Sergio Conti, Antonio DeSimone, and Yvo Pokern, Energetics and switchingof quasi-uniform states in small ferromagnetic particles , ESAIM: Mathematical Modellingand Numerical Analysis-Mod´elisation Math´ematique et Analyse Num´erique (2004), no. 2,235–248.[6] Fran¸cois Alouges and Pascal Jaisson, Convergence of a finite element discretization for theLandau-Lifshitz equations in micromagnetism , Mathematical Models and Methods in AppliedSciences (2006), no. 02, 299–316.[7] Fran¸cois Alouges, Evaggelos Kritsikis, Jutta Steiner, and Jean-Christophe Toussaint, A con-vergent and precise finite element scheme for Landau-Lifschitz-Gilbert equation , NumerischeMathematik (2014), no. 3, 407–430.[8] Fran¸cois Alouges, Evaggelos Kritsikis, and Jean-Christophe Toussaint,
A convergent finiteelement approximation for Landau-Lifschitz-Gilbert equation , Physica B: Condensed Matter (2012), no. 9, 1345–1349.[9] Fran¸cois Alouges and Alain Soyeur,
On global weak solutions for Landau-Lifshitz equations:existence and nonuniqueness , Nonlinear Analysis: Theory, Methods & Applications (1992),no. 11, 1071–1084.[10] S¨oren Bartels, Stability and convergence of finite-element approximation schemes for har-monic maps , SIAM journal on numerical analysis (2005), no. 1, 220–238.[11] S¨oren Bartels and Andreas Prohl, Convergence of an implicit finite element method for theLandau-Lifshitz-Gilbert equation , SIAM journal on numerical analysis (2006), no. 4, 1405–1419.[12] , Constraint preserving implicit finite element discretization of harmonic map flowinto spheres , Mathematics of Computation (2007), no. 260, 1847–1859.[13] James L Blue and MR Scheinfein, Using multipoles decreases computation time for magne-tostatic self-energy , IEEE Transactions on Magnetics (1991), no. 6, 4778–4780.[14] Dietrich Braess, Finite elements: Theory, fast solvers, and applications in solid mechanics ,Cambridge University Press, 2007.[15] W. F. Brown,
Magnetostatic Principles , North-Holland Amsterdam, 1962.[16] Xavier Brunotte, G´erard Meunier, and Jean-Fran¸cois Imhoff,
Finite element modeling ofunbounded problems using transformations: a rigorous, powerful and easy solution , IEEETransactions on Magnetics (1992), no. 2, 1663–1666.[17] Snorre H Christiansen, A div-curl lemma for edge elements , SIAM Journal on NumericalAnalysis (2005), no. 1, 116–126.[18] Ivan Cimr´ak, A survey on the numerics and computations for the Landau-Lifshitz equationof micromagnetism , Archives of Computational Methods in Engineering (2007), no. 3,1–37.[19] , Convergence result for the constraint preserving mid-point scheme for micromag-netism , Journal of Computational and Applied Mathematics (2009), no. 1, 238–246.[20] Francisco de la Hoz, Carlos J Garc´ıa-Cervera, and Luis Vega,
A numerical study of the self-similar solutions of the Schr¨odinger map , SIAM Journal on Applied Mathematics (2009),no. 4, 1047–1077.[21] Lukas Exl, Winfried Auzinger, Simon Bance, Markus Gusenbauer, Franz Reichel, and ThomasSchrefl, Fast stray field computation on tensor grids , Journal of Computational Physics (2012), no. 7, 2840–2850.[22] Josef Fidler and Thomas Schrefl,
Micromagnetic modelling-the current state of the art , Jour-nal of Physics D: Applied Physics (2000), no. 15, R135.[23] DR Fredkin and TR Koehler, Hybrid method for computing demagnetizing fields , IEEE Trans-actions on Magnetics (1990), no. 2, 415–417.[24] Atsushi Fuwa, Tetsuya Ishiwata, and Masayoshi Tsutsumi, Finite difference scheme for theLandau–Lifshitz equation , Japan journal of industrial and applied mathematics (2012),no. 1, 83–110.[25] Carlos J Garc´ıa-Cervera, Numerical micromagnetics: A review (2007).
ONVERGENCE OF A MASS-LUMPED FINITE ELEMENT METHOD 21 [26] Carlos J Garc´ıa-Cervera and Weinan E,
Improved Gauss-Seidel projection method for micro-magnetics simulations , IEEE transactions on magnetics (2003), no. 3, 1766–1770.[27] Carlos J Garc´ıa-Cervera and Alexandre M Roma, Adaptive mesh refinement for micromag-netics simulations , IEEE Transactions on Magnetics (2006), no. 6, 1648–1654.[28] Boling Guo and Min-Chun Hong, The Landau-Lifshitz equation of the ferromagnetic spinchain and harmonic maps , Calculus of Variations and Partial Differential Equations (1993),no. 3, 311–334.[29] Stephen Gustafson, Kenji Nakanishi, and Tai-Peng Tsai, Asymptotic Stability, Concentration,and Oscillation in Harmonic Map Heat-Flow, Landau-Lifshitz, and Schr¨odinger Maps on R ,Communications in Mathematical Physics (2010), no. 1, 205–242.[30] E. Hairer and G. Wanner, Solving ordinary differential equations II, stiff and differential-algebraic problems , 2nd ed., Springer, New York, 1996.[31] J Samuel Jiang, Hans G Kaper, and Gary K Leaf,
Hysteresis in layered spring magnets ,Discrete and Continuous Dynamical Systems Series B (2001), 219–323.[32] Stavros Komineas, Rotating vortex dipoles in ferromagnets , Physical Review Letters (2007), no. 11, 117202.[33] Stavros Komineas and Nikos Papanicolaou, Skyrmion dynamics in chiral ferromagnets , Phys-ical Review B (2015), no. 6, 064412.[34] Perinkulam S Krishnaprasad and Xiaobo Tan, Cayley transforms in micromagnetics , PhysicaB: Condensed Matter (2001), no. 1, 195–199.[35] Evaggelos Kritsikis, A Vaysset, LD Buda-Prejbeanu, F Alouges, and J-C Toussaint,
Beyondfirst-order finite element schemes in micromagnetics , Journal of Computational Physics (2014), 357–366.[36] Martin Kruzik and Andreas Prohl,
Recent developments in the modeling, analysis, and nu-merics of ferromagnetism , SIAM review (2006), no. 3, 439–483.[37] Lev D Landau and E Lifshitz, On the theory of the dispersion of magnetic permeability inferromagnetic bodies , Phys. Z. Sowjetunion (1935), no. 153, 101–114.[38] Debra Lewis and Nilima Nigam, Geometric integration on spheres and some interestingapplications , Journal of Computational and Applied Mathematics (2003), no. 1, 141–170.[39] Boris Livshitz, Amir Boag, H Neal Bertram, and Vitaliy Lomakin,
Nonuniform grid algo-rithm for fast calculation of magnetostatic interactions in micromagnetics , Journal of AppliedPhysics (2009), no. 7, 07D541.[40] HH Long, ET Ong, ZJ Liu, and EP Li,
Fast fourier transform on multipoles for rapid calcu-lation of magnetostatic fields , IEEE Transactions on Magnetics (2006), no. 2, 295–300.[41] Jacques E Miltat and Michael J Donahue, Numerical micromagnetics: Finite difference meth-ods , Handbook of magnetism and advanced magnetic materials (2007).[42] Per-Olof Persson and Gilbert Strang,
A simple mesh generator in MATLAB , SIAM review (2004), no. 2, 329–345.[43] N Popovi´c and Dirk Praetorius, Applications of H -matrix techniques in micromagnetics ,Computing (2005), no. 3, 177–204.[44] Thomas Schrefl, Finite elements in numerical micromagnetics: Part i: granular hard mag-nets , Journal of magnetism and magnetic materials (1999), no. 1, 45–65.[45] Luc Tartar,
The general theory of homogenization: a personalized introduction , Vol. 7,Springer Science & Business Media, 2009.[46] Igor Tsukerman, Alexander Plaks, and H Neal Bertram,
Multigrid methods for computationof magnetostatic fields in magnetic recording problems , Journal of Applied Physics (1998),no. 11, 6344–6346.[47] Xiao-Ping Wang, Carlos J Garcıa-Cervera, and E Weinan, A Gauss–Seidel projection methodfor micromagnetics simulations , Journal of Computational Physics (2001), no. 1, 357–372.[48] Samuel W Yuan and H Neal Bertram,
Fast adaptive algorithms for micromagnetics , IEEETransactions on Magnetics (1992), no. 5, 2031–2036. Department of Mathematics, University of California, Berkeley, CA 94720 and MS-B284, Los Alamos National Laboratory, Los Alamos, NM 87544
E-mail address : [email protected] Department of Mathematics and Lawrence Berkeley National Laboratory, Univer-sity of California, Berkeley, CA 94720
E-mail address ::