1 Introduction
Electrons in a material pose local magnetic orders, but typically do not exhibit a spontaneous macroscopic magnetic ordering unless a collective motion of these local magnetic orders is present. This results in a net magnetization even in the absence of an external magnetic field. Such a material is called a ferromagnet. It has binary stable configurations, which makes it an ideal material for data recording and storage. Recent advances in experiment and theory Žutić et al. (2004) have demonstrated effective and precise control of ferromagnetic configurations by means of external fields.
A very common phenomenological model for magnetization dynamics is the LandauLifshitzGilbert (LLG) equation Landau and Lifshitz (1935); Gilbert (1955). This model has been successfully used to interpret various experimental observations. The LLG equation is technically quasilinear, nonlocal and has a nonconvex constraint, which poses interesting challenges in designing efficient numerical methodologies. In addition, the magnetization reversal process requires numerical methods to resolve different length and temporal scales in the presence of domain walls and vortices, due to their important roles in the switching process (Shi et al., 1999, 2000). Therefore, numerical methods for the LLG equation with high accuracy and efficiency are highly demanding.
There has been a continuous progress of developing numerical algorithms in the past few decades; see for example (Kruzík and Prohl, 2006; Cimrák, 2008) and references therein. The spatial derivative is typically approximated by the finite element method (FEM) (Abert et al., 2014; Bruckner et al., 2014; Baňas et al., 2014; Aurada et al., 2015; Le et al., 2015; Baňas et al., 2015; Praetorius et al., 2018; Alouges and Jaisson, 2006; Alouges, 2008; Alouges et al., 2014a, b) and the finite difference method (E and Wang, 2000; Fuwa et al., 2012; Jeong and Kim, 2010; Kim and Lipnikov, 2017; Romeo et al., 2008). As for the temporal discretization, explicit schemes (Alouges and Jaisson, 2006; Jiang et al., 2001), fully implicit schemes (Prohl, 2001; Bartels and Prohl, 2006; Fuwa et al., 2012), and semiimplicit schemes (E and Wang, 2000; Wang et al., 2001; GarcíaCervera and E, 2003; Lewis and Nigam, 2003; Cimrák, 2005; Gao, 2014; Boscarino et al., 2016) have been extensively studied. Explicit schemes suffer from severe stability constraints. Fully implicit schemes can overcome this severe stability constraints. However, a (nonsymmetric) nonlinear system of equations needs to be solved at each time step, which is timeconsuming. A nonlinear multigrid method is used to handle the nonlinearity at each time step in (Jeong and Kim, 2014), and the fixed point iteration technique is used to deal with the nonlinearity in (Cimrák and Slodička, 2004). In (Fuwa et al., 2012), the existence and uniqueness of a solution to the nonlinear system is proved under the condition that the temporal stepsize be proportional to the square of the spatial gridsize. This, however, is contrary to the unconditional stability of implicit schemes.
Semiimplicit schemes achieve a desired balance between stability and efficiency. One of the most popular methods is the GaussSeidel projection method (GSPM) developed by Wang, GarcíaCervera, and E (E and Wang, 2000; Wang et al., 2001; GarcíaCervera and E, 2003). This method is based on a combination of a GaussSeidel implementation of a fractional step implicit solver for the gyromagnetic term, and the projection method for the heat flow of harmonic maps to overcome the difficulties associated with the stiffness and nonlinearity. Only several linear systems of equations need to be solved at each step, whose complexity is comparable to solving the scalar heat equation implicitly. It is tested that GSPM is unconditionally stable with firstorder accuracy in time. In order to get secondorder accuracy in time, two nonsymmetric linear systems of equations need to be solved at each step. Note that a projection step is needed to preserve the pointwise length constraint.
In this work, we propose two secondorder semiimplicit projection methods for LLG equation based on the secondorder backward differentiation formula and the secondorder interpolation formula using the information at previous two temporal steps. The unconditional unique solvability of both methods is proved, with their secondorder accuracy verified through numerical examples in both 1D and 3D. The efficiency of both methods is compared to that of another two popular schemes in the literature. In addition, we test the robustness of both methods using the first benchmark problem for a ferromagnetic thin film material developed by the micromagnetic modeling activity group from National Institute of Standards and Technology (NIST). It is worth mentioning that a modification of the proposed method has been proved to be secondorder accurate in time under mild conditions Chen et al. (2019).
The rest of the paper is organized as follows. In Section 2, we first introduce the micromagnetics model based on the LLG equation. The secondorder semiimplicit projection methods are described in Section 2.2 with their unique solvability given in Section 2.3. The calculation of the demagnetization field (stray filed) is given in Section 2.4. Numerical results in Section 3 are used to test the accuracy and the efficiency of the proposed methods in both 1D and 3D. Moreover, the first benchmark for a ferromagnetic thin film material developed by the micromagnetics modeling activity group from NIST is used to check the applicability of the proposed methods in Section 4. Conclusions are drawn in Section 5.
2 Secondorder semiimplicit methods
2.1 LandauLifshitzGilbert equation
The magnetization dynamics in a ferromagnetic material are described by the LLG equation Landau and Lifshitz (1935); Jr. (1963), which take the following nondimensionalized form:
(1) 
with
(2) 
on . Here the magnetization
is a threedimensional vector field with a pointwise constraint
and is the unit outward normal vector. is a bounded domain occupied by the ferromagnetic material. The first term on the right hand side in eq. 1 is the gyromagnetic term and the second term is the damping term with being the dimensionless damping coefficient.The effective field consists of the exchange field, the anisotropy field, the external field and the demagnetization or stray field . For a uniaxial material,
(3) 
Here, the dimensionless parameters are and with the diameter of the ferromagnetic body, the anisotropy constant, the exchange constant, the permeability of vacuum, and the saturation magnetization. , and denotes the standard Laplacian operator. is the applied (external) magnetic field and the detailed description of will be given in §2.4. Typical values of the physical parameters for Permalloy are included in Table 1.
Physical Parameters for Permalloy  

2.2 Secondorder semiimplicit projection methods
Denote the temporal stepsize by , the spatial mesh size by , the standard secondorder centered difference for Laplacian operator by , and , with the final time. For convenience, we use the finite difference method to approximate the spatial derivatives in eqs. 6 and 5. For the temporal discretization, we employ the secondorder backward differentiation formulas (BDFs) to approximate the temporal derivative
(7) 
Such a discretization results the following fully implicit scheme for eq. 5:
(8)  
As expected, at each time step, a nonlinear system of equations needs to be solved in eq. 8. Moreover, the nonsymmetric structure of the system introduces additional difficulties.
To overcome this severe difficulty, we approximate the nonlinear prefactors in front of the discrete Laplacian term using the information from previous time steps (onesided interpolation) with its accuracy the same as the corresponding BDF scheme. For (8), we have
(9) 
where
(10)  
(11) 
and . However, such a scheme cannot preserve the magnitude of magnetization, we therefore add a projection step and obtain the following scheme for eq. 5
where is the intermediate magnetization. For eq. 6, using the same idea, we have
Remark 1.
Using the same idea, we can construct highorder semiimplicit projection schemes for eq. 5 and eq. 6 using highorder BDFs and highorder onesided interpolations. However, if the order is greater than , then such a scheme cannot be Astable. Numerical tests for LLG equation also indicates the conditional stability of higher order semiimplicit projection schemes. Firstorder semiimplicit projection schemes are Astable, but they do not have obvious advantages over the GaussSeidel projection method Wang et al. (2001).
Remark 2.
To kick start Scheme A and Scheme B in implementation, we use the firstorder semiimplicit projection scheme with the firstorder BDF and the firstorder onesided interpolation for the first temporal step, and thus the whole method is still secondorder accurate.
2.3 Unconditionally unique solvability
For simplicity of illustration, we assume that the spatial mesh size . An extension to the general case is straightforward.
We firstly introduce the discrete inner product.
Definition 1 (Discrete inner product).
For grid functions and over the uniform numerical grid, we define
where is the index set and is the index which closely depends on the dimension .
Definition 2.
For the grid function , we define the average of summation as
Definition 3.
For the grid function with , we define the discrete norm as
For ease of notation, we drop the temporal indices and rewrite Scheme A and Scheme B in a more compact form
(12)  
(13) 
where , , and are given.
Theorem 2.1 (Solvability for Scheme A).
Given and , the numerical scheme eq. 12 admits a unique solution.
For the unique solvability analysis for eq. 12, we denote . Note that , due to the Neumann boundary condition for . Meanwhile, we observe that in general, since . Instead, we rewrite eq. 12 into
and take the average on both sides. Therefore, can be represented as follows:
and is given by eq. 10. In turn, eq. 12 is then rewritten as
(14) 
To proceed, we need the following lemma.
Lemma 2.0 (BrowderMinty lemma Browder (1963); Minty (1963)).
Let 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.
Proof of Theorem 2.1.
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 , so that . Moreover, for any , with , we get
and the equality only holds when .
Therefore, an application of Lemma 2.2 implies a unique solution of Scheme A.
∎
Theorem 2.3 (Solvability for Scheme B).
Given and , the numerical scheme eq. 13 admits a unique solution.
Proof of Theorem 2.3.
We first rewrite eq. 13 in a compact form , where
Here
is the identity matrix,
is the antisymmetric matrix corresponding to the discrete operator , and is the symmetric positive definite matrix corresponding to which admits a decomposition with being nonsingular.Thus, we have
Due to the antisymmetry of matrix , we have
It follows from the spectral lemma for antisymmetric matrices (Prasolov, 1994)
that the eigenvalues of
are either or purely imaginary, thus the eigenvalues of are either or purely imaginary as well. This unique solvability comes as a consequence of the fact that all eigenvalues of have as real parts and . ∎Remark 3.
Note that the unique solvability of Scheme A and Scheme B does not impose any condition on and . This is in contrast with earlier results for the fully implicit schemes where is needed for the unique solution of the nonlinear system of equations at each time step; see (Fuwa et al., 2012) for example.
Remark 4.
2.4 Computation of the stray field
The stray field with the scalar function in which satisfies
(15) 
together with jump conditions
(16) 
Here denotes the jump of at the material boundary as
(17) 
and is defined similarly. The solution to eq. 15  section 2.4 is
and thus the stray field
(18) 
where is the Newtonian potential.
It follows from eq. 18 combined with the divergence theorem that
(19) 
The evaluation of the stray field can be carried out by performing an integration over the entire material for every point . In terms of the computational complexity, a direct evaluation of (18) requires with
the degree of freedoms. Moreover, we need to evaluate (
18) at each time step. Therefore, the direct evaluation is computationally expensive and thus a fast solver is highly desirable. The complexity for solving stray field using FFT is Wang et al. (2001); GarcíaCervera and E (2003).3 Accuracy and efficiency test
We use examples in both 1D and 3D to show the secondorder accuracy of Scheme A and Scheme B.
In order to have the exact magnetization profile, we consider the simplified LLG equation with only the exchange term and set with the forcing term
(20) 
where with the exact solution. The ferromagetic body in 1D and in 3D. The final time . Since the exchange term is the stiffest term in the original LLG equation, it is adequate to use (20) to test accuracy and efficiency of the proposed methods.
The exact solution in 1D is
(21) 
which satisfies the homogeneous Neumann boundary condition.
The exact solution in 3D is
(22) 
where , , and .
Since both schemes are semiimplicit, we compare their efficiency with another two semiimplicit methods: GaussSeidel projection method Wang et al. (2001) and the secondorder implicitexplicit method (Boscarino et al., 2016). For completeness, we first state these two methods.
3.1 The GaussSeidel projection method
It follows from eqs. 5 and 4 that the GaussSeidel projection method (GSPM) Wang et al. (2001) is given as the following three steps:

Implicit GaussSeidel:
(23) (24) 
Heat flow without constraints :
(25) (26) 
Projection onto :
(27)
where denotes the intermediate values of . As in GarcíaCervera and E (2003), the stray field is computed using the intermediate values in eq. 23 and eq. 25.
3.2 The secondorder implicitexplicit method
For simplicity, we only formulate the secondorder implicitexplicit (IMEX2) scheme (Boscarino et al., 2016) without Gilbert damping. Define
and
The secondorder IMEX scheme (IMEX2) reads as
(28) 
where . Note that IMEX2 solves two linear systems of equations at each step.
Remark 5.
The computational cost for GSPM, IMEX2, and BDF2 at each temporal step is as follows. Seven symmetric linear systems of equations with constant coefficients and dimension need to be solved for GSPM, and two nonsymmetric linear systems of equations with variable coefficients and dimension need to be solved for IMEX2. BDF2 only needs to solve one nonsymmetric linear systems of equations with variable coefficients and dimension for both Scheme A and Scheme B. Here is the number of unknowns in each spatial dimension.
3.3 Accuracy test
Since Scheme A and Scheme B are comparable numerically, we only show results of Scheme A, termed as BDF2. For comparison, we also list results of the other two semiimplicit methods: GSPM and IMEX2.
In 1D, for (21), we fix and record the temporal error in terms of the temporal stepsize in Table 2 and Figure (a)a to get the temporal accuracy. Both BDF2 and IMEX2 are secondorder accurate, while GSPM is firstorder accurate. To get the spatial accuracy, we fix and record the spatial error in terms of in Table 3 and Figure (b)b. BDF2, IMEX2, and GSPM are all secondorder accurate.
BDF2  GSPM  IMEX2  

2.0D2  1.0753D4  3.2024D2  4.1119D5 
1.0D2  2.7384D5  1.6289D2  1.3107D5 
5.0D3  6.8538D6  7.9890D3  3.7623D6 
2.5D3  1.6513D6  4.1923D3  7.2361D7 
1.25D3  3.4152D7  2.0662D3  1.2711D7 
order  2.065  0.987  2.085 
BDF2  GSPM  IMEX2  

4.0D2  6.2092D4  6.2094D4  6.2092D4 
2.0D2  1.5516D4  1.5517D4  1.5516D4 
1.0D2  3.8789D5  3.8805D5  3.8789D5 
5.0D3  9.6973D6  9.7145D6  9.6972D6 
2.5D3  2.4243D6  2.4439D6  2.4243e6 
order  2.000  1.998  2.000 
In 3D, for (22), we fix and record the temporal error in terms of in table 4 and fig. (a)a to get the temporal accuracy. Both BDF2 and IMEX2 are secondorder accurate, while GSPM is firstorder accurate. To get the spatial accuracy, we fix the temporal stepsize and record the spatial error in terms of in table 5 and fig. (b)b. BDF2, IMEX2, and GSPM are all secondorder accurate.
BDF2  GSPM  IMEX2  

1/8  3.1853D3  1/32  3.0024D4  1/4  4.1759D3 
1/16  9.3167D4  1/64  1.5214D4  1/8  1.0659D3 
1/32  2.4929D4  1/128  7.7368D5  1/16  2.6611D4 
1/64  6.1027D5  1/256  3.8171D5  1/32  6.3198D5 
1/128  1.1785D5  1/512  2.0295D5  1/64  1.2107D5 
order  2.009    0.977    2.094 
BDF2  GSPM  IMEX2  

1/2  3.5860D4  3.5856D4  3.5860D4 
1/4  8.1251D5  8.1245D5  8.1252D5 
1/8  2.0187D5  2.0187D5  2.0189D5 
1/16  5.0532D6  5.0519D6  5.0552D6 
1/32  1.2635D6  1.3278D6  1.2656D6 
order  2.030  2.016  2.030 
3.4 Efficiency comparison
To compare the efficiency, we plot the CPU time (in seconds) of BDF2, GSPM and IMEX2 in terms of the error in fig. (a)a for the 1D case and in fig. (b)b for the 3D case.
In 1D, for the same tolerance, costs of BDF2, GSPM, and IMEX2 in fig. (a)a satisfy: BDF2IMEX2GSPM. In 3D, for the same tolerance, costs of BDF2, GSPM, and IMEX2 in fig. (b)b satisfy: BDF2IMEX2GSPM. For both cases, BDF2 is slightly better than IMEX2 since two linear systems of equations need to be solved in IMEX2 while only one linear system needs to solved in BDF2. Both BDF2 and IMEX2 are better than GSPM.
Remark 6.
When the spatial mesh is very fine, we observe that to achieve the same tolerance, costs of BDF2, GSPM, and IMEX2 satisfy: GSPMBDF2IMEX2. The reason is that fast solvers can solve symmetric linear systems with constant coefficients in GSPM, while nonsymmetric linear systems with variable coefficients are involved in both BDF2 and IMEX2. It becomes increasingly difficult to solve such systems using the Generalized Minimum Residual Method (GMRES), for example. This issue will be further explored in a subsequent work.
4 Benchmark problem from NIST
To examine our methods in the realistic case, we simulate the first standard problem established by the micromagnetic modelling activity group at National Institute of Standards and Technology (NIST) Micromagnetic Modeling Activity Group (2000). This problem asks for simulating the hysteresis loop of a thinfilm element with material parameters that are not too different from Permalloy. The hysteresis loop is obtained in the following way: A positive external field of strength , in the unit of is applied. The magnetization is able to reach a steady state. Once this steady state is approached, the applied external field is reduced by a certain amount, and the material sample is again allowed to reach a steady state. The process continues until the hysteresis system attains a negative field of strength . The process then is repeated, increasing the field in small steps until it reaches the initial applied external field. As a consequence, we are able to plot the average magnetization at the steady state as a function of the external filed strength during the hysteresis loop. For BDF2, the unsymmetric linear systems of equations are solved by GMRES from library called high performance preconditioners Center for Applied Scientific Computing, Lawrence Livermore National Laboratory (1952) which was developed by Center for Applied Scientific Computing, Lawrence Livermore National Laboratory.
4.1 Magnetization profile
For comparison, we use the same setup of the available code mo96a of the first standard problem from NIST. Its setup is grid points and the canting angle of applied field from nominal axis. The calculation of the demagnetization field is done by FFT. The initial state is uniform. In the loop, successive steps between and are adopted for both xloop and yloop.
Due to the presence of metastable symmetric states, the applied fields should be rotated one degree counterclockwise off the nominal axis. The damping coefficient , the temporal stepsize and the cell size is . A stopping criterion is used to determine that a steady state is reached when the relative change in the total energy is less than . For mo96a, LABEL:NIST_long_magnetization and LABEL:NIST_short_magnetization plot the average remanent magnetization on the bottom surface of the sample in the xy plane when when the applied fields are approximately parallel (canting angle ) to the y (long) axis and the x (short) axis, respectively. For BDF2, the corresponding results are shown in LABEL:BDF2_long_magnetization and LABEL:BDF2_short_magnetization, respectively. The inplane magnetization components are represented by arrows in LABEL:NIST_BDF2_magnetization.
Furthermore, magnetization components are visualized by the grayscale value in fig. 5. For mo96a, fig. (a)a  fig. (d)d plot the x component and the y component when the applied field is along the y axis, the x component and the y component when the applied field is along the x axis, respectively. Results of BDF2 are shown in Figures (h)h, (g)g, (f)f and (e)e.
From LABEL:NIST_BDF2_magnetization and fig. 5, we observe that results of BDF2 are in qualitative agreements with those of mo96a.
4.2 Hysteresis loop
Hysteresis loops generated by the code mo96a are shown in fig. (a)a when the applied field is approximately parallel to the long axis and in fig. (b)b when the applied field is approximately parallel to the short axis, respectively. The average remanent magnetization in reduced units is for the yloop and for the xloop. The coercive fields are in fig. (a)a and in fig. (b)b.
For BDF2, hysteresis loops are presented in fig. (c)c when the applied field is approximately parallel to the long axis and in fig. (d)d when the applied field is approximately parallel to the short axis, respectively. The average remanent magnetization in reduced units is for the yloop and for the xloop. The coercive fields are in fig. (c)c and in fig. (d)d.
Based on these results, we conclude that results of BDF2 agree well with those of NIST, both qualitatively and quantitatively.
5 Conclusions
In this paper, we proposed two secondorder semiimplicit projection methods to solve the LandauLifshitzGilbert equation, which possess unconditionally unique solvability at each time step. Examples in both 1D and 3D are used to verify the accuracy and the efficiency of both schemes. Additionally, the first benchmark problem from NIST is used to check the applicability of our methods in the realistic situation using the full LandauLifshitizGilbert equation. Results show that our schemes can produce the correct hysteresis loop with quantities of interest agreeing with other methods in a quantitative manner.
One issue associated to the proposed methods is that nonsymmetric linear systems with variable coefficients are involved at each time step. It becomes increasingly difficult to solve such systems using the Generalized Minimum Residual Method. However, such linear systems have some unique features, as a consequence of the unique structure at the continuous level. This shall be used to develop more efficient linear solvers. Meanwhile, the technique presented here may be applicable to the model for currentdriven domain wall dynamics Chen et al. (2015) and the SchrödingerLandauLifshitz system Chen et al. (2016), which shall be explored later.
Acknowledgments
This work is supported in part by the grants NSFC 21602149 (J. Chen), AFORS grant FA95501810095 (C. J. GarcíaCervera), NSF DMS1418689 (C. Wang), the Postgraduate Research and Practice Innovation Program of Jiangsu Province via grant KYCX19_1947 (C. Xie), and NSFC 11801016 (Z. Zhou).
References
 Žutić et al. (2004) I. Žutić, J. Fabian, S. Das Sarma, Spintronics: Fundamentals and applications, Rev. Mod. Phys. 76 (2004) 323–410.
 Landau and Lifshitz (1935) L. Landau, E. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies, Phys. Z. Sowjet. 8 (1935) 153–169.
 Gilbert (1955) T. Gilbert, Phys. Rev. 100 (1955) 1243. [Abstract only; full report, Armor Research Foundation Project No. A059, Supplementary Report, May 1, 1956 (unpublished)].
 Shi et al. (1999) J. Shi, S. Tehrani, M. R. Scheinfein, Magnetization vortices and anomalous switching in patterned NiFeCo submicron arrays, Appl. Phys. Lett. 74 (1999) 2525–2527.
 Shi et al. (2000) J. Shi, T. Zhu, S. Tehrani, Y. Zheng, J.G. Zhu, Geometry dependence of magnetization vortices in patterned submicron NiFe elements, Appl. Phys. Lett. 76 (2000) 2588–2590.
 Kruzík and Prohl (2006) M. Kruzík, A. Prohl, Recent developments in the modeling, analysis, and numerics of ferromagnetism, SIAM Rev. 48 (2006) 439–483.
 Cimrák (2008) I. Cimrák, A survey on the numerics and computations for the LandauLifshitz equation of micromagnetism, Arch. Comput. Method Eng. 15 (2008) 277–309.
 Abert et al. (2014) C. Abert, G. Hrkac, M. Page, D. Praetorius, M. Ruggeri, D. Suess, Spinpolarized transport in ferromagnetic multilayers: An unconditionally convergent FEM integrator, Comput. Math. Appl. 68 (2014) 639–654.
 Bruckner et al. (2014) F. Bruckner, D. Suess, M. Feischl, T. Führer, P. Goldenits, M. Page, D. Praetorius, M. Ruggeri, Multiscale modeling in micromagnetics: Existence of solutions and numerical integration, Math. Models Methods Appl. Sci. 24 (2014) 2627–2662.
 Baňas et al. (2014) L. Baňas, M. Page, D. Praetorius, J. Rochat, A decoupled and unconditionally convergent linear FEM integrator for the Landau–Lifshitz–Gilbert equation with magnetostriction, IMA J. Numer. Anal. 34 (2014) 1361–1385.
 Aurada et al. (2015) M. Aurada, J. Melenk, D. Praetorius, FEM–BEM coupling for the largebody limit in micromagnetics, J. Comput. Appl. Math. 281 (2015) 10–31.
 Le et al. (2015) K.N. Le, M. Page, D. Praetorius, T. Tran, On a decoupled linear FEM integrator for eddycurrentLLG, Appl. Anal. 94 (2015) 1051–1067.
 Baňas et al. (2015) L. Baňas, M. Page, D. Praetorius, A convergent linear finite element scheme for the MaxwellLandauLifshitzGilbert equation, Electron. Trans. Numer. Anal. 44 (2015) 250–270.
 Praetorius et al. (2018) D. Praetorius, M. Ruggeri, B. Stiftner, Convergence of an implicit–explicit midpoint scheme for computational micromagnetics, Comput. Math. Appl. 75 (2018) 1719–1738.
 Alouges and Jaisson (2006) F. Alouges, P. Jaisson, Convergence of a finite element discretization for the LandauLifshitz equations in micromagnetism, Math. Models Methods Appl. Sci. 16 (2006) 299–316.
 Alouges (2008) F. Alouges, A new finite element scheme for LandauLifchitz equations, Discrete Contin. Dyn. Syst. Ser. S 1 (2008) 187–196.
 Alouges et al. (2014a) F. Alouges, E. Kritsikis, J. Steiner, J.C. Toussaint, A convergent and precise finite element scheme for Landau–Lifschitz–Gilbert equation, Numer. Math. 128 (2014a) 407–430.
 Alouges et al. (2014b) F. Alouges, A. de Bouard, A. Hocquet, A semidiscrete scheme for the stochastic Landau–Lifshitz equation, Stoch. Partial Differ. Equ. Anal. Comput. 2 (2014b) 281–315.
 E and Wang (2000) W. E, X. Wang, Numerical methods for the LandauLifshitz equation, SIAM J. Numer. Anal. 38 (2000) 1647–1665.
 Fuwa et al. (2012) A. Fuwa, T. Ishiwata, M. Tsutsumi, Finite difference scheme for the LandauLifshitz equation, Japan J. Indust. Appl. Math. 29 (2012) 83–110.
 Jeong and Kim (2010) D. Jeong, J. Kim, A Crank–Nicolson scheme for the Landau–Lifshitz equation without damping, J. Comput. Appl. Math. 234 (2010) 613–623.
 Kim and Lipnikov (2017) E. Kim, K. Lipnikov, The mimetic finite difference method for the Landau–Lifshitz equation, J. Comput. Phys. 328 (2017) 109–130.
 Romeo et al. (2008) A. Romeo, G. Finocchio, M. Carpentieri, L. Torres, G. Consolo, B. Azzerboni, A numerical solution of the magnetization reversal modeling in a permalloy thin film using fifth order Runge–Kutta method with adaptive step size control, Phys. B: Condens. Matter 403 (2008) 464 – 468.
 Jiang et al. (2001) J. S. Jiang, H. G. Kaper, G. K. Leaf, Hysteresis in layered spring magnets, Discrete Continuous Dyn. Syst. Ser. B. 1 (2001) 219–232.
 Prohl (2001) A. Prohl, Computational micromagnetism, volume 1, Springer, 2001.
 Bartels and Prohl (2006) S. Bartels, A. Prohl, Convergence of an implicit finite element method for the LandauLifshitzGilbert equation, SIAM J. Numer. Anal. 44 (2006) 1405–1419.
 Wang et al. (2001) X. Wang, C. J. GarcíaCervera, W. E, A GaussSeidel projection method for micromagnetics simulations, J. Comput. Phys. 171 (2001) 357–372.
Comments
There are no comments yet.