Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation
Panchi Li Changjian Xie Rui Du Jingrun Chen Xiao-Ping Wang
math.NA cs.NA 35Q99 65Z05 65M06








1]Panchi Li
LiPanchi1994@163.com
1]Changjian Xie

20184007005@stu.suda.edu.cn
1,2]Rui Ducor1
[cor1]Corresponding authors
durui@suda.edu.cn
1,2]Jingrun Chencor1
jingrunchen@suda.edu.cn
3]Xiao-Ping Wangcor1
mawang@ust.hk

[1]School of Mathematical Sciences, Soochow University, Suzhou, 215006, China.
[2]Mathematical Center for Interdisciplinary Research, Soochow University, Suzhou, 215006, China.
[3]Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China



Micromagnetic simulation is an important tool to study various dynamic behaviors of magnetic order in ferromagnetic materials.
The underlying model is the Landau-Lifshitz-Gilbert equation, where the magnetization dynamics is driven by
the gyromagnetic torque term and the Gilbert damping term. Numerically, considerable progress has been made in the past decades.
One of the most popular methods is the Gauss-Seidel projection method developed by Xiao-Ping Wang, Carlos García-Cervera, and Weinan E in 2001.
It first solves a set of heat equations with constant coefficients and updates the gyromagnetic term in the Gauss-Seidel manner,
and then solves another set of heat equations with constant coefficients for the damping term.
Afterwards, a projection step is applied to preserve the length constraint in the pointwise sense.
This method has been verified to be unconditionally stable numerically and successfully applied to study magnetization dynamics
under various controls.

In this paper, we present two improved Gauss-Seidel projection methods with unconditional stability.
The first method updates the gyromagnetic term and the damping term simultaneously and follows
by a projection step. The second method introduces two sets of approximate solutions, where
we update the gyromagnetic term and the damping term simultaneously for one set of approximate
solutions and apply the projection step to the other set of approximate solutions in an alternating manner.
Compared to the original Gauss-Seidel projection method which has to solve heat equations 7 times at each time step,
the improved methods solve heat equations 5 times and 3 times, respectively.
First-order accuracy in time and second-order accuracy in space are verified by examples in both 1D and 3D.
In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed
numerically. Application of both methods to a realistic material is also presented with hysteresis loops
and magnetization profiles. Compared with the original method, the recorded running times suggest that savings
of both methods are about 2/7 and 4/7 for the same accuracy requirement, respectively.









Landau-Lifshitz-Gilbert equationGauss-Seidel projection methodunconditional stabilitymicromagnetic simulation

[2000] 35Q99 65Z05 65M06









§ INTRODUCTION


In ferromagnetic materials, the intrinsic magnetic order, known as magnetization = (M_1, M_2, M_3)^T,  is modeled by
the following Landau-Lifshitz-Gilbert (LLG) equation <cit.>

    ∂𝐌/∂ t = -γ𝐌×ℋ -
    γα/M_s𝐌×(𝐌×ℋ)

with γ the gyromagnetic ratio and |𝐌| = M_s the saturation magnetization.
On the right-hand side of (<ref>), the first term is the gyromagnetic
term and the second term is the Gilbert damping term with α the dimensionless
damping coefficient <cit.>. Note that the gyromagnetic term is a conservative
term, whereas the damping term is a dissipative term.
The local field ℋ = -δ F/δ𝐌 is computed from the Landau-Lifshitz energy functional

    F[𝐌] = 1/2∫_Ω{A/M_s^2|∇𝐌|^2 +
    Φ(𝐌/M_s)
     -2μ_0𝐇_e·𝐌}d
    + μ_0/2∫_ℝ^3|∇ U|^2 d,

where A is the exchange constant, A/M_s^2|∇𝐌| is the
exchange interaction energy; Φ(𝐌/M_s)
is the anisotropy energy, and for simplicity the material is assumed to be uniaxial with
Φ(𝐌/M_s) = K_u/M_s^2(M_2^2+M_3^2) with
K_u the anisotropy constant; -2μ_0𝐇_e·𝐌 is the Zeeman energy due to the
external field with μ_0 the permeability of vacuum. Ω is the volume occupied by the material.
The last term in (<ref>) is the energy resulting from the field
induced by the magnetization distribution inside the material.
This stray field 𝐇_s = -∇ U where U() satisfies

    U() = ∫_Ω∇ N(-)·𝐌()d,

where N(-) = -1/4π1/|-| is the Newtonian potential.

For convenience, we rescale the original LLG equation (<ref>) by changes of variables
t →(μ_0γ M_s)^-1t and  x→ Lx with L the diameter of Ω.
Define 𝐦=/M_s and 𝐡 = M_sℋ. The dimensionless LLG equation reads as

    ∂𝐦/∂ t = -𝐦×𝐡 -
    α𝐦×(𝐦×𝐡),

where

    𝐡 = -Q(m_2𝐞_2 + m_3𝐞_3) + ϵΔ𝐦
    + 𝐡_e + 𝐡_s

with dimensionless parameters Q=K_u/(μ_0M_s^2) and ϵ=A/(μ_0M_s^2L^2).
Here 𝐞_2=(0,1,0), 𝐞_3=(0,0,1).
Neumann boundary condition is used

    ∂𝐦/∂ν_|∂Ω = 0,

where ν is the outward unit normal vector on ∂Ω.

The LLG equation is a weakly nonlinear equation. In the absence of Gilbert damping,
α = 0, equation (<ref>) is a degenerate equation of parabolic type and is
related to the sympletic flow of harmonic maps <cit.>. In the large
damping limit, α→∞, equation (<ref>) is related to
the heat flow for harmonic maps <cit.>. It is easy to check that |𝐦| = 1
in the pointwise sense in the evolution. All these properties possesses interesting
challenges for designing numerical methods to solve the LLG equation.
Meanwhile, micromagnetic simulation is an important tool to study magnetization
dynamics of magnetic materials <cit.>.
Over the past decades, there has been increasing progress on numerical methods
for the LLG equation; see <cit.> for reviews and
references therein. Finite difference method and finite element method have been
used for the spatial discretization.

For the temporal discretization, there are explicit schemes such as Runge-Kutta
methods <cit.>. Their stepsizes are subject to strong stability constraint.
Another issue is that the length of magnetization cannot be preserved and thus a projection
step is needed. Implicit schemes <cit.>
are unconditionally stable and usually can preserve the length of magnetization automatically.
The difficulty of implicit schemes is how to solve a nonlinear system of equations at each step.
Therefore, semi-implicit methods <cit.>
provide a compromise between stability and the difficult for solving the equation at each step.
A projection step is also needed to preserve the length of magnetization.

Among the semi-implicit schemes, the most popular one is the Gauss-Seidel projection method (GSPM) proposed
by Wang, García-Cervera, and E <cit.>.
GSPM first solves a set of heat equations with constant coefficients and updates the gyromagnetic term in the Gauss-Seidel manner,
and then solves another set of heat equations with constant coefficients for the damping term.
Afterwards, a projection step is applied to preserve the length of magnetization.
GSPM is first-order accurate in time and has been verified to be unconditionally stable numerically.


In this paper, we present two improved Gauss-Seidel projection methods with unconditional stability.
The first method updates the gyromagnetic term and the damping term simultaneously and follows
by a projection step. The second method introduces two sets of approximate solutions, where
we update the gyromagnetic term and the damping term simultaneously for one set of approximate
solutions and apply the projection step to the other set of approximate solutions in an alternating manner.
Compared to the original Gauss-Seidel projection method, which solves heat equations 7 times at each time step,
the improved methods solve heat equations 5 times and 3 times, respectively.
First-order accuracy in time and second-order accuracy in space are verified by examples in both 1D and 3D.
In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed
numerically. Application of both methods to a realistic material is also presented with hysteresis loops
and magnetization profiles. Compared with the original method, the recorded running times suggest that savings
of both methods are about 2/7 and 4/7 for the same accuracy requirement, respectively.

The rest of the paper is organized as follows. For completeness and comparison, we first introduce
GSPM in <Ref>. Two improved GSPMs are presented in <Ref>.
Detailed numerical tests are given in <Ref>, including accuracy check and efficiency check
in both 1D and 3D, unconditional stability with respect to both the grid size and the damping parameter,
hysteresis loops, and magnetization profiles. Conclusions are drawn in <Ref>.




§ GAUSS-SEIDEL PROJECTION METHOD FOR LANDAU-LIFSHITZ-GILBERT EQUATION


Before the introduction of the GSPM <cit.>,
we first use the finite difference method for spatial discretization. Figure <ref> shows
a schematic picture of spatial grids in 1D. Let i=0,1,⋯,M,M+1,
j=0,1,⋯,N,N+1, and k=0,1,⋯,K,K+1 be the indices of grid points in 3D.


Second-order centered difference
for Δ𝐦 reads as

    Δ_h𝐦_i,j,k   =𝐦_i+1,j,k-
    	2𝐦_i,j,k+𝐦_i-1,j,k/Δ x^2
       +𝐦_i,j+1,k-2𝐦_i,j,k+
    	𝐦_i,j-1,k/Δ y^2
       +𝐦_i,j,k+1-2𝐦_i,j,k+
    	𝐦_i,j,k-1/Δ z^2,

where 𝐦_i,j,k=𝐦((i - 1/2)Δ x,
(j - 1/2)Δ y, (k - 1/2)Δ z). For the Neumann boundary condition, a second-order approximation yields

    𝐦_0,j,k    =𝐦_1,j,k,  𝐦_M,j,k  =𝐦_M+1,j,k,   j = 1,⋯,N,k=1,⋯,K, 
    𝐦_i,0,k    =𝐦_i,1,k,  𝐦_i,N,k  =𝐦_i,N+1,k,   i = 1,⋯,M,k=1,⋯,K, 
    𝐦_i,j,0    =𝐦_i,j,1,  𝐦_i,j,K  =𝐦_i,j,K+1,   i = 1,⋯,M,j=1,⋯,N.


To illustrate the main ideas, we first consider the following simplified equation

    𝐦_t = - 𝐦×Δ𝐦 -
    α𝐦×(𝐦×Δ𝐦),

which can be rewritten as

    𝐦_t = - 𝐦×Δ𝐦 - α𝐦
    (𝐦·Δ𝐦) + αΔ𝐦.


We split (<ref>) into two equations

    𝐦_t     = - 𝐦×Δ𝐦, 
    𝐦_t     =  αΔ𝐦.


However, (<ref>) is still nonlinear. Therefore, we consider a fractional step scheme
to solve (<ref>)

    𝐦^* - 𝐦^n/Δ t   = Δ_h𝐦^*
    𝐦^n+1   = 𝐦^n - 𝐦^n×𝐦^*

or

    𝐦^n+1 = 𝐦^n - 𝐦^n×(I-Δ t
    Δ_h)^-1𝐦^n,

where I is the identity matrix. This scheme is subject to strong stability constraint, and thus the implicit Gauss-Seidel scheme is introduced to overcome this issue.
Let

    g_i^n = (I-Δ tΔ_h)^-1m_i^n,   i = 1,2,3.

We then have

    [ m_1^n+1; m_2^n+1; m_3^n+1 ]=
    [         m_1^n + (g_2^nm_3^n - g_3^nm_2^n);     m_2^n + (g_3^nm_1^n+1 - g_1^n+1m_3^n); m_3^n + (g_1^n+1m_2^n+1 - g_2^n+1m_1^n+1) ].

This scheme solve (<ref>) with unconditional stability.
(<ref>) is linear heat equation which can be solved easily.
However, the splitting scheme (<ref>) - (<ref>)
cannot preserve || = 1, and thus a projection step needs to
be added.

For the full LLG equation (<ref>), the GSPM works as follows.
Define

    𝐡 = ϵΔ𝐦 + 𝐟̂,

where 𝐟̂= -Q(m_2𝐞_2 + m_3𝐞_3) + 𝐡_e + 𝐡_s.

The original GSPM <cit.> solves the equation (<ref>) in three steps:


  * Implicit Gauss-Seidel

    g_i^n    = (I - Δ tϵΔ_h)^-1(m_i^n + Δ tf̂_i^n),  i = 2,3, 
    g_i^*    = (I - Δ tϵΔ_h)^-1(m_i^* + Δ tf̂_i^n),  i = 1,2,


    [ m_1^*; m_2^*; m_3^* ] =
    [ m_1^n + (g_2^nm_3^n - g_3^nm_2^n); m_2^n + (g_3^nm_1^* - g_1^*m_3^n); m_3^n + (g_1^*m_2^* - g_2^*m_1^*) ].


  * Heat flow without constraints

    𝐟̂^* = -Q(m_2^*𝐞_2+m_3^*𝐞_3) + 𝐡_e + 𝐡_s^n,


    [ m_1^**; m_2^**; m_3^** ] =
    [ m_1^* + αΔ t(ϵΔ_hm_1^** + f̂_1^*); m_2^* + αΔ t(ϵΔ_hm_2^** + f̂_2^*); m_3^* + αΔ t(ϵΔ_hm_3^** + f̂_3^*) ].


  * Projection onto S^2

    [ m_1^n+1; m_2^n+1; m_3^n+1 ] = 1/|𝐦^**|[ m_1^**; m_2^**; m_3^** ].


Here the numerical stability of the original GSPM <cit.> was founded to be independent of gridsizes
but depend on the damping parameter α. This issue was solved in <cit.> by replacing
(<ref>) and (<ref>) with

    g_i^* = (I - Δ tϵΔ_h)^-1(m_i^* + Δ tf̂_i^*),  i = 1,2,

and

    𝐟̂^* = -Q(m_2^*𝐞_2+m_3^*𝐞_3) + 𝐡_e + 𝐡_s^*,

respectively.
Update of the stray field is done using fast Fourier transform <cit.>.
It is easy to see that the GSPM solves 7 linear systems of equations with constant coefficients
and updates the stray field using FFT 6 times at each step.




§ TWO IMPROVED GAUSS-SEIDEL PROJECTION METHODS FOR LANDAU-LIFSHITZ-GILBERT EQUATION


Based on the description of the original GSPM in <Ref>, we introduce two improved GSPMs
for LLG equation. The first improvement updates both the gyromagnetic term and the damping term
simultaneously, termed as Scheme A. The second improvement introduces two sets of approximate solution
with one set for implicit Gauss-Seidel step and the other set for projection in an alternating manner,
termed as Scheme B. Details are given in below.




 §.§ Scheme A

The main improvement of Scheme A over the original GSPM is the combination of (<ref>) - (<ref>),
or (<ref>) - (<ref>).


  * Implicit-Gauss-Seidel

    g_i^n    = (I - Δ tϵΔ_h)^-1(m_i^n+Δ tf̂_i^n),   i = 1,2,3, 
    g_i^*   = (I - Δ tϵΔ_h)^-1(m_i^*+Δ tf̂_i^*),   i = 1,2,


    [ m_1^*; m_2^*; m_3^* ] =
    [ m_1^n - (m_2^ng_3^n - m_3^ng_2^n) -
    α(m_1^ng_1^n+m_2^ng_2^n + m_3^ng_3^n)m_1^n + α g_1^n; m_2^n - (m_3^ng_1^* - m_1^*g_3^n) -
    α(m_1^*g_1^*+m_2^ng_2^n + m_3^ng_3^n)m_2^n + α g_2^n; m_3^n - (m_1^*g_2^* - m_2^*g_1^*) -
    α(m_1^*g_1^*+m_2^*g_2^* + m_3^ng_3^n)m_3^n + α g_3^n ].


  * Projection onto S^2

    [ m_1^n+1; m_2^n+1; m_3^n+1 ] = 1/|𝐦^*|[ m_1^*; m_2^*; m_3^* ].


It is easy to see that Scheme A solves 5 linear systems of equations with constant coefficients
and uses FFT 5 times at each step.




 §.§ Scheme B


The main improvement of Scheme B over Scheme A is the introduction of two sets of approximate solutions,
one for (<ref>) - (<ref>) and the other for (<ref>) and the update
of these two sets of solutions in an alternating manner.

Given the initialized 𝐠^0

    g_i^0 = (I - Δ tϵΔ_h)^-1(m_i^0+Δ tf̂_i^0),  i = 1,2,3,

Scheme B works as follows


  * Implicit Gauss-Seidel

    g_i^n+1   = (I - Δ tϵΔ_h)^-1(m_i^*+Δ tf̂_i^*),  i = 1,2,3


    m_1^*   = m_1^n - (m_2^ng_3^n - m_3^ng_2^n) -
    α(m_1^ng_1^n+m_2^ng_2^n + m_3^ng_3^n)m_1^n + 
          α((m_1^n)^2 + (m_2^n)^2 + (m_3^n)^2)g_1^n
    
    m_2^*   = m_2^n - (m_3^ng_1^n+1 - m_1^*g_3^n) -
    α(m_1^*g_1^n+1+m_2^ng_2^n + m_3^ng_3^n)m_2^n + 
          α((m_1^*)^2 + (m_2^n)^2 + (m_3^n)^2)g_2^n 
    
    m_3^*   = m_3^n - (m_1^*g_2^n+1 - m_2^*g_1^n+1) -
    α(m_1^*g_1^n+1+m_2^*g_2^n+1 + m_3^ng_3^n)m_3^n + 
          α((m_1^*)^2 + (m_2^*)^2 + (m_3^n)^2)g_3^n


  * Projection onto S^2

    [ m_1^n+1; m_2^n+1; m_3^n+1 ] = 1/|𝐦^*|[ m_1^*; m_2^*; m_3^* ].


Here one set of approximate solution {^*} is updated in the implicit Gauss-Seidel step
and the other set of approximate solution {^n+1} is updated in the projection step.
Note that (<ref>) is defined only for {^*} which can be used in two successive temporal steps, and thus only 3 linear systems of equations
with constant coefficients are solved at each step and 3 FFT executions are used for
the stray field. The length of magnetization can be preserved in the time evolution.

The computational cost of GSPM and its improvements comes from solving
the linear systems of equations with constant coefficients.
To summarize, we list the number of linear systems of equations to be solved and
the number of FFT executions to be used at each step for the original GSPM <cit.>,
Scheme A, and Scheme B in  <ref>.
The savings represent the ratio between costs of two improved schemes over that of the original GSPM.









§ NUMERICAL EXPERIMENTS


In this section, we compare the original GSPM <cit.>, Scheme A, and Scheme B
via a series of examples in both 1D and 3D, including accuracy check and efficiency check, unconditional stability
with respect to both the grid size and the damping parameter, hysteresis loops, and magnetization profiles.
For convenience, we define

    ratio-i = Time(GSPM) - Time(Scheme  i)/Time(GSPM),

for i = A and B, which quantifies the improved efficiency of Scheme A and Scheme B over the original GSPM <cit.>.




 §.§ Accuracy Test


[1D case]

In 1D, we choose the exact solution over the unit interval Ω=[0,1]

    𝐦_e = (cos(x̅)sin(t), sin(x̅)sin(t), cos(t)),

which satisfies

    _t = - ×_xx - α×(×_xx) + 𝐟

with x̅=x^2(1-x)^2, and 𝐟=𝐦_et +
𝐦_e×𝐦_exx + α𝐦_e×(𝐦_e×𝐦_exx).
Parameters are α=0.00001 and T=5.0e-2.

We first show the error 𝐦_e - 𝐦_h_∞ with 𝐦_h being the numerical solution
with respect to the temporal stepsize Δ t and the spatial stepsize Δ x.
As shown in <ref>(a) and <ref>(c), suggested by the least squares fitting,
both first-order accuracy in time and second-order accuracy in space
are observed. Meanwhile, we record the CPU time as a function of accuracy (error) by varying
the temporal stepsize and the spatial stepsize in <ref>(b) and <ref>(d),
<ref> and <ref>, respectively. In addition, from <ref>
and <ref>, the saving of Scheme A over GSPM is about 2/7, which equals 1-5/7,
and the saving of Scheme B over GSPM is about 4/7, respectively. This observation is in good agreement
with the number of linear systems being solved at each step for these three methods, as shown in <ref>.









[3D case]

In 3D, we choose the exact solution over Ω=[0,2]×[0,1]×[0,0.2]

    𝐦_e = (cos(x̅y̅z̅)sin(t),sin(x̅y̅z̅)sin(t), cos(t)),

which satisfies

    _t=-×Δ - α×(×Δ) + 𝐟

with x̅=x^2(1-x)^2, y̅=y^2(1-y)^2, z̅=z^2(1-z)^2 and 𝐟=𝐦_et +
𝐦_e×Δ𝐦_e + α𝐦_e×(𝐦_e×Δ𝐦_e).
Parameters are T = 1.0e-05 and α = 0.01.

Like in the 1D case, we first show the error 𝐦_e - 𝐦_h_∞ with 𝐦_h being the numerical solution
with respect to the temporal stepsize Δ t and the spatial stepsize Δ x.
As shown in <ref>(a) and <ref>(c), suggested by the least squares fitting,
both first-order accuracy in time and second-order accuracy in space
are observed. Meanwhile, we record the CPU time as a function of accuracy (error) by varying
the temporal stepsize and the spatial stepsize in <ref>(b) and <ref>(d),
<ref> and <ref>, respectively. In addition, from <ref>
and <ref>, the saving of Scheme A over GSPM is about 2/7,
and the saving of Scheme B over GSPM is about 4/7, respectively. This observation is in good agreement
with the number of linear systems being solved at each step for these three methods, as shown in <ref>.









It worths mentioning that all these three methods are tested to be unconditionally stable with respect to the spatial gridsize and the temporal stepsize.




 §.§ Micromagnetic Simulations


To compare the performance of Scheme A and Scheme B with GSPM, we have
carried out micromagnetic simulations of the full LLG equation with realistic material parameters.
In all our following simulations, we consider a thin film ferromagnet of size Ω=1 μm×1 μm×0.02 μm
with the spatial gridsize 4 nm×4 nm×4 nm and the temporal stepsize =1 picosecond.
The demagnetization field (stray field) is calculated via FFT <cit.>.



  §.§.§ Comparison of hysteresis loops

The hysteresis loop is calculated in the following way. First, a positive external field
H_0=μ_0 H is applied and the system is allowed to reach a stable state. Afterwards,
the external field is reduced by a certain amount and the system is relaxed to a stable state again.
The process continues until the external field attains a negative field of strength H_0.
Then the external field starts to increase and the system relaxes until the initial applied external
field H_0 is approached. In the hysteresis loop, we can monitor the magnetization dynamics and
plot the average magnetization at the stable state as a function of the strength of the external field.
The stopping criterion for a steady state is that the relative change of the total energy is less than 10^-7.
The applied field is parallel to the x axis. The initial state we take is the uniform state and the damping parameter α=0.1.

In <ref>, we compare the average magnetization in the hysteresis loop simulated by GSPM, Scheme A and
Scheme B. Profiles of the average magnetization of these three methods are in quantitative agreements with
approximately the same switch field 9 (± 0.4) mT.






  §.§.§ Comparison of magnetization profiles


It is tested that GSPM in <cit.> was unstable with a very small damping parameter α
and was resolved in <cit.>.
This section is devoted to the unconditional stability of Scheme A and Scheme B with respect to α. We consider a thin film ferromagnet of size 1 μm×1 μm×0.02 μm
with the spatial gridsize  4 nm×4 nm× 4 nm and the temporal stepsize is 1 picosecond.
Following <cit.>, we consider the full LLG equation with α=0.1 and α=0.01 and without the external field.
The initial state is _0=(0,1,0) if x∈[0,L_x/5]∪[4L_x/5,L_x] and _0=(1,0,0) otherwise. The final time is 10 ns.
In <Ref>, we present a color plot of the angle between the in-plane magnetization and the x axis,
and an arrow plot of the in-plane magnetization for the original GSPM <cit.>, Scheme A, and Scheme B, respectively.
In these figures, α=0.1 is presented in the top row and α=0.01 is presented in the bottom row;
a color plot of the angle between the in-palne magnetization and the x axis is presented in the left column
and  an arrow plot of the in-plane magnetization is presented in the right column.











§ CONCLUSION


In this paper, based on the original Gauss-Seidel projection methods,
we present two improved Gauss-Seidel projection methods with the first-order accuracy in time
and the second-order accuracy in space.
The first method updates the gyromagnetic term and the damping term simultaneously and follows
by a projection step, which requires to solve heat equations 5 times at each time step.
The second method introduces two sets of approximate solutions, where
we update the gyromagnetic term and the damping term simultaneously for one set of approximate
solutions and apply the projection step to the other set of approximate solutions in an alternating manner.
Therefore, only 3 heat equations are needed to be solved at each step.
Compared to the original Gauss-Seidel projection method, which solves heat equations 7 times at each step,
savings of these two improved methods are about 2/7 and 4/7, which is verified by
both 1D and 3D examples for the same accuracy requirement.
In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed
numerically. Application of both methods to a realistic material is also presented with hysteresis loops
and magnetization profiles.



§ ACKNOWLEDGMENTS

This work is supported in part by the grants NSFC 21602149 (J. Chen), NSFC 11501399 (R. Du), the Hong Kong Research Grants Council (GRF grants 16302715,
16324416, 16303318 and NSFC-RGC joint research grant N-HKUST620/15) (X.-P. Wang), and the Innovation Program for postgraduates in Jiangsu province via grant KYCX19_1947 (C. Xie).


model1-num-names


