IMEX-RK methods for Landau-Lifshitz equation with arbitrary damping^†^†thanks: Received date, and accepted date (The correct dates will be entered by the editor).

Yan Gui School of Mathematical Sciences, Soochow University, Suzhou,China, ([email protected]). Cheng Wang Mathematics Department, University of Massachusetts, North Dartmouth, MA 02747, USA, ([email protected]). Jingrun Chen School of Mathematical Sciences, University of Science and Technology of China, Hefei, China and Suzhou Institute for Advanced Research, University of Science and Technology of China, Suzhou, China, ([email protected]).

Abstract

Magnetization dynamics in ferromagnetic materials is modeled by the Landau-Lifshitz (LL) equation, a nonlinear system of partial differential equations. Among the numerical approaches, semi-implicit schemes are widely used in the micromagnetics simulation, due to a nice compromise between accuracy and efficiency. At each time step, only a linear system needs to be solved and a projection is then applied to preserve the length of magnetization. However, this linear system contains variable coefficients and a non-symmetric structure, and thus an efficient linear solver is highly desired. If the damping parameter becomes large, it has been realized that efficient solvers are only available to a linear system with constant, symmetric, and positive definite (SPD) structure. In this work, based on the implicit-explicit Runge-Kutta (IMEX-RK) time discretization, we introduce an artificial damping term, which is treated implicitly. The remaining terms are treated explicitly. This strategy leads to a semi-implicit scheme with the following properties: (1) only a few linear system with constant and SPD structure needs to be solved at each time step; (2) it works for the LL equation with arbitrary damping parameter; (3) high-order accuracy can be obtained with high-order IMEX-RK time discretization. Numerically, second-order and third-order IMEX-RK methods are designed in both the 1-D and 3-D domains. A comparison with the backward differentiation formula scheme is undertaken, in terms of accuracy and efficiency. The robustness of both numerical methods is tested on the first benchmark problem from National Institute of Standards and Technology. The linearized stability estimate and optimal rate convergence analysis are provided for an alternate IMEX-RK2 numerical scheme as well.

keywords:

Micromagnetics simulation; Landau-Lifshitz equation; Implicit-explicitly Runge-Kutta scheme; Second-order accuracy; Third-order accuracy; Hysteresis loop.

{AMS}

35K61; 65N06; 65N12

1 Introduction

Ferromagnetic materials have the intrinsic magnetic order (magnetization), whose dynamics is modeled by the Landau-Lifshitz (LL) equation [1, 2]. The equation stands for a non-local nonlinear problem with non-convex constraint and possible degeneracy. The past few decades have witnessed the progress of numerical methods for the LL equation; see [4, 44, 46, 45] for reviews and references therein. There are explicit algorithms (e.g. [35, 6]), fully implicit ones (e.g. [5, 16]), and semi-implicit schemes (e.g. [15, 26, 37, 38, 39, 3, 43]). An explicit method updates the magnetization without any need to solver linear or nonlinear system of equations, while it suffers from the severe stability constraint on the temporal step-size. Implicit schemes are unconditionally stable and preserve the physical properties of the LL equation, such as energy dissipation and length conservation, while a nonlinear system of equation needs to be solved at each time step. Semi-implicit schemes only solve a linear system of equations at each time step and are unconditionally stable in micromagnetics simulations. Therefore, a semi-implicit method provides a nice balance between accuracy and efficiency.

One typical semi-implicit method is the Gauss-Seidel projection method (GSPM) [15, 26, 36, 47]. Using the vectorial structure of LL equation, GSPM achieves unconditional stability in micromagnetics simulations; and only a few linear systems need to be solved, with constant, symmetric, and positive definite coefficients. However, the temporal accuracy of GSPM is limited to the first-order. Another semi-implicit approach is based on the backward differentiation formula (BDF) for temporal derivative and the one-sided extrapolation for nonlinear terms [27, 40, 48]. High-order accuracy can be obtained in time using the BDF approach. However, only first-order and second-order BDFs are unconditionally stable. Meanwhile, the linear system of equations has variable coefficients and non-symmetric structure, thus no fast solver is available. Meanwhile, for time-dependent nonlinear partial differential equations in general, implicit-explicit (IMEX) schemes have been extensively applied [39]. The basic idea is to treat dominant linear term implicitly and the remaining nonlinear terms explicitly. For the LL equation, the second-order IMEX has been studied in [27]. Two linear systems, with variable coefficients and non-symmetric structure, need to be solved. Thus IMEX2 can hardly compete with BDF2 in terms of accuracy and efficiency.

In this work, we propose an implicit-explicit Runge-Kutta (IMEX-RK) scheme for solving the LL equation, based on the recent development of IMEX-RK method for the nonlinear diffusion equation [30]. The basic idea is to introduce an artificial linear diffusion term and treat it implicitly. All the remaining terms are treated explicitly. RK methods are employed for the time discretization. Only a few linear systems, with constant coefficients and SPD structure, need to be solved. In the existing literature, such linear numerical schemes have only be reported for the large damping parameter [41]. Instead, the IMEX-RK method works for the LL equation general damping parameters, which is very important since the damping parameter may be small in most magnetic materials [49]. Moreover, higher-order numerical schemes could be constructed using RK approaches. Numerical results have demonstrated an advantage of IMEX-RK schemes over the BDF2 approach in terms of accuracy and efficiency. The performance of IMEX-RK schemes has also been verified over a large range of artificial damping parameters and the first benchmark problem for a ferromagnetic thin film material from National Institute of Standards and Technology (NIST).

Because of the robust numerical results, a theoretical analysis of the proposed IMEX-RK numerical schemes is highly desired. However, such an analysis turns out to be very challenging, due to the multi-stage nature, as well as the highly complicated nonlinear terms in the vector form. Some convergence analysis works have been reported for the IMEX-RK numerical methods to various nonlinear PDEs in the existing literature, such as the convection-diffusion equation [29], the porous media equation [28, 31], etc. Meanwhile, the degree of nonlinearity of the LL equation is even higher than these reported PDE models, and the associated theoretical analysis is expected to be more challenging. In this article, we choose an alternate IMEX-RK2 numerical algorithm, in which the explicit part satisfies the strong stability-preserving (SSP) property [11, 17] (denoted as the SSP-IMEX-RK2 scheme), and provide a linearized stability estimate and convergence analysis for the SSP-IMEX-RK2 method to a simplified version of the LL equation. Many state-of-arts techniques, such as a rough error estimate to control the discrete $\ell^{\infty}$ and $W_{h}^{1,8}$ norms of the numerical solution, combined with a refined error estimate, have to be applied to obtain an optimal rate convergence analysis for the SSP-IMEX-RK2 scheme to the nonlinear LL equation.

The rest of paper is organized as follows. In Section 2, we introduce the LL model and then propose the IMEX-RK schemes, including the SSP-IMEX-RK2 algorithm. For the convenience of comparison, we also briefly review the BDF2 method. Accuracy and efficiency tests of the IMEX-RK schemes are provided in Section 3 with a detailed check for the dependence on the artificial damping parameter. Section 4 is devoted to the micromagnetics simulations of the IMEX-RK methods, including equilibrium structures and the first benchmark problem from NIST. The linearized stability estimate and the convergence analysis for the SSP-IMEX-RK2 scheme to a simplified version of nonlinear LL equation are provided in Section 5. Finally, some concluding remarks are made in Section 6.

2 The model and the proposed methods

2.1 Landau-Lifshitz equation

The LL equation is a phenomenological model for magnetization dynamics in a ferromagnetic material, which takes the following form

\displaystyle\boldsymbol{m}_{t}=-\boldsymbol{m}\times\boldsymbol{h}_{\mathrm{eff}}-\alpha\boldsymbol{m}\times\left(\boldsymbol{m}\times\boldsymbol{h}_{\mathrm{eff}}\right),

(2.1)

with the homogeneous Neumann boundary condition

\displaystyle{{\left.\frac{\partial\boldsymbol{m}}{\partial\nu}\right|}_{\partial\Omega}}=0.

(2.2)

Here $\Omega$ is a bounded domain occupied by the ferromagnetic material, and the magnetization $\boldsymbol{m}:\Omega\subset{\mathbb{R}^{d}}\to{\mathbb{S}^{2}},d=1,2,3$ is a 3-D vector field with the length constraint $\left|\boldsymbol{m}\right|=1$ , and $\nu$ is the unit outward normal vector along $\partial\Omega$ . The first term of the right hand side in (2.1) is the gyromagnetic term, while the second term represents the damping term with a dimensionless parameter $\alpha>0$ .

For a uniaxial material, the free energy is given by

\displaystyle F[\boldsymbol{m}]=\frac{\mu_{0}M_{s}^{2}}{2}\int_{\Omega}\left(\epsilon|\nabla\boldsymbol{m}|^{2}+Q\left(m_{2}^{2}+m_{3}^{2}\right)-\boldsymbol{h}_{s}\cdot\boldsymbol{m}-2\boldsymbol{h}_{e}\cdot\boldsymbol{m}\right)\mathrm{d}\boldsymbol{x},

in which the listed terms correspond to the exchange energy, the anisotropy energy, the magnetostatic energy, and the Zeeman energy, respectively. The effective field ${{\boldsymbol{h}}_{\text{eff}}}$ can be obtained by taking the variation of $F[\boldsymbol{m}]$ with respect to $\boldsymbol{m}$ , and consists of the exchange field, the anisotropy field, the stray field ${{\boldsymbol{h}}_{\text{s}}}$ , and the external field ${{\boldsymbol{h}}_{\text{e}}}$ of the following form

\displaystyle{{\boldsymbol{h}}_{\text{eff}}}=\epsilon\Delta\boldsymbol{m}-Q({{m}_{2}}{{\boldsymbol{e}}_{2}}+{{m}_{3}}{{\boldsymbol{e}}_{3}})+{{\boldsymbol{h}}_{s}}+{{\boldsymbol{h}}_{e}}.

Table 1: Typical values of the physical parameters for Permalloy, which is an alloy of Nickel (80%) and Iron (20%) frequently used in magnetic storage devices.

Physical Parameters for Permalloy
$K_{u}$	$1.0\times 10^{2}\;\mathrm{J/m^{3}}$
$C_{ex}$	$1.3\times 10^{-11}\;\mathrm{J/m}$
$M_{s}$	$8\times 10^{5}\;\mathrm{A/m}$
$\mu_{0}$	$4\pi\times 10^{-7}\;\mathrm{N/A^{2}}$
$\alpha$	$0.01$

In the above representation, $Q={{K}_{u}}/({{\mu}_{0}}M_{s}^{2})$ and $\epsilon={{C}_{ex}}/({{\mu}_{0}}M_{s}^{2}{{L}^{2}})$ are the dimensionless parameters with ${{C}_{ex}}$ the exchange constant, ${{K}_{u}}$ the anisotropy constant, $L$ the diameter of ferromagnetic body, ${{\mu}_{0}}$ the permeability of vacuum, and $M_{s}$ the saturation magnetization, respectively.. Typical values of the physical parameters for Permalloy are included as shown in Table 1.

The two unit vectors ${{\boldsymbol{e}}_{2}}={{(0,1,0)}^{T}}$ , ${{\boldsymbol{e}}_{3}}={{(0,0,1)}^{T}}$ , and $\Delta$ stands for the standard Laplacian operator. The stray field ${{\boldsymbol{h}}_{\text{s}}}$ takes the form

\displaystyle{{\boldsymbol{h}}_{s}}=-\nabla\int_{\Omega}{\nabla N(\boldsymbol{x}-\boldsymbol{y})\cdot\boldsymbol{m}(\boldsymbol{y})d\boldsymbol{y}},

where $N(\boldsymbol{x})=-\frac{1}{4\pi\left|\boldsymbol{x}\right|}$ is the Newtonian potential.

For simplicity, we denote

\displaystyle\emph{ {f}}=-Q({{m}_{2}}{{\boldsymbol{e}}_{2}}+{{m}_{3}}{{\boldsymbol{e}}_{3}})+{{\boldsymbol{h}}_{s}}+{{\boldsymbol{h}}_{e}}.

(2.3)

Accordingly, the LL equation can be rewritten as

\displaystyle{{\boldsymbol{m}}_{t}}=-\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\emph{ {f}})-\alpha\boldsymbol{m}\times\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\emph{ {f}}).

(2.4)

Thanks to the constraint $|\boldsymbol{m}|=1$ , we obtain an equivalent form

\displaystyle\boldsymbol{m}_{t}=\alpha(\epsilon\Delta\boldsymbol{m}+\boldsymbol{f})+\alpha\left(\epsilon|\nabla\boldsymbol{m}|^{2}-\boldsymbol{m}\cdot\boldsymbol{f}\right)\boldsymbol{m}-\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\boldsymbol{f}).

(2.5)

Some notations are introduced for discretization and numerical approximation. The 1-D domain is set as $\Omega=(0,1)$ , and the 3-D version becomes $\Omega=(0,1)^{3}$ , and the final time is denoted as $T$ . In the 1-D domain, we divide $\Omega$ into $N$ equal parts with $h=1/N$ . Fig. 1 displays a schematic picture of 1-D spatial grids, with ${{x}_{i-\frac{1}{2}}}=(i-\frac{1}{2})h,i=1,2,\cdots,N$ .

Refer to caption — Figure 1: Spatial grids in 1D where ${{x}_{-\frac{1}{2}}}$ and ${{x}_{N+\frac{1}{2}}}$ are two ghost points.

The 3-D grid points in can be similarly constructed. For convenience, we set $h_{x}=h_{y}=h_{z}=h$ , and ${{\mathbf{m}}_{i,j,k}}=\mathbf{m}((i-\frac{1}{2})h,(j-\frac{1}{2})h,(k-\frac{1}{2})h),0\leq i,j,k\leq N+1$ .

The second-order centered difference for $\Delta\mathbf{m}$ in the 3-D domain is formulated as

\begin{split}{{\Delta}_{h}}{{\mathbf{m}}_{i,j,k}}=&\frac{{{\mathbf{m}}_{i+1,j,k}}-2{{\mathbf{m}}_{i,j,k}}+{{\mathbf{m}}_{i-1,j,k}}}{{h^{2}}}\\ &+\frac{{{\mathbf{m}}_{i,j+1,k}}-2{{\mathbf{m}}_{i,j,k}}+{{\mathbf{m}}_{i,j-1,k}}}{{h^{2}}}\\ &+\frac{{{\mathbf{m}}_{i,j,k+1}}-2{{\mathbf{m}}_{i,j,k}}+{{\mathbf{m}}_{i,j,k-1}}}{{h^{2}}}.\end{split}

The second-order approximation of Neumann boundary condition in (2.2) gives

\displaystyle\begin{aligned} &{{\mathbf{m}}_{0,j,k}}={{\mathbf{m}}_{1,j,k}},\,\,\,{{\mathbf{m}}_{N,j,k}}={{\mathbf{m}}_{N+1,j,k}},\,\,\,j,k=1,\cdots,N,\\ &{{\mathbf{m}}_{i,0,k}}={{\mathbf{m}}_{i,1,k}},\,\,\,{{\mathbf{m}}_{i,N,k}}={{\mathbf{m}}_{i,N+1,k}},\,\,\,i,k=1,\cdots,N,\\ &{{\mathbf{m}}_{i,j,0}}={{\mathbf{m}}_{i,j,1}},\,\,\,{{\mathbf{m}}_{i,j,N}}={{\mathbf{m}}_{i,j,N+1}},\,\,\,i,j=1,\cdots,N.\end{aligned}

(2.6)

The discrete gradient operator $\nabla_{h}\boldsymbol{m}$ with $\boldsymbol{m}=(u,v,w)^{T}$ becomes

\nabla_{h}\boldsymbol{m}_{i,j,k}=\left[\begin{array}[]{lll}\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2h}&\frac{v_{i+1,j,k}-v_{i-1,j,k}}{2h}&\frac{w_{i+1,j,k}-w_{i-1,j,k}}{2h}\\ \frac{u_{i,j+1,k}-u_{i,j-1,k}}{2h}&\frac{v_{i,j+1,k}-v_{i,j-1,k}}{2h}&\frac{w_{i,j+1,k}-w_{i,j-1,k}}{2h}\\ \frac{u_{i,j,k+1}-u_{i,j,k-1}}{2h}&\frac{v_{i,j,k+1}-v_{i,j,k-1}}{2h}&\frac{w_{i,j,k+1}-w_{i,j,k-1}}{2h}\end{array}\right].

For temporal discretization, we set ${{t}^{n}}=nk$ with $k$ the step-size and $n\leq\left[\frac{T}{k}\right]$ .

2.2 Second-order backward differentiation formula method

The BDF numerical method [27, 40] is based on the BDF temporal discretization and the one-sided extrapolation. For comparison, we recall the BDF2 algorithm as

\displaystyle\left\{\begin{array}[]{c}\dfrac{\frac{3}{2}\tilde{\boldsymbol{m}}_{h}^{n+2}-2\boldsymbol{m}_{h}^{n+1}+\frac{1}{2}\boldsymbol{m}_{h}^{n}}{k}=-\hat{\boldsymbol{m}}_{h}^{n+2}\times\left(\epsilon\Delta_{h}\tilde{\boldsymbol{m}}_{h}^{n+2}+\hat{\boldsymbol{f}}_{h}^{n+2}\right)\\ -\alpha\hat{\boldsymbol{m}}_{h}^{n+2}\times\left(\hat{\boldsymbol{m}}_{h}^{n+2}\times\left(\epsilon\Delta_{h}\tilde{\boldsymbol{m}}_{h}^{n+2}+\hat{\boldsymbol{f}}_{h}^{n+2}\right)\right)\\ \boldsymbol{m}_{h}^{n+2}=\frac{\tilde{\boldsymbol{m}}_{h}^{n+2}}{\left|\tilde{\boldsymbol{m}}_{h}^{n+2}\right|},\end{array}\right.

(2.10)

where $\tilde{\boldsymbol{m}}_{h}^{n+2}$ is an intermediate magnetization, and $\hat{\boldsymbol{m}}_{h}^{n+2},\hat{\boldsymbol{f}}_{h}^{n+2}$ are given by the following extrapolation formula:

\displaystyle\hat{\boldsymbol{m}}_{h}^{n+2}

\displaystyle=2\boldsymbol{m}_{h}^{n+1}-\boldsymbol{m}_{h}^{n},\quad\hat{\boldsymbol{f}}_{h}^{n+2}

\displaystyle=2\boldsymbol{f}_{h}^{n+1}-\boldsymbol{f}_{h}^{n},

with $\boldsymbol{f}_{h}^{n}=-Q\left(m_{2}^{n}\boldsymbol{e}_{2}+m_{3}^{n}\boldsymbol{e}_{3}\right)+\boldsymbol{h}_{s}^{n}+\boldsymbol{h}_{e}^{n}$ , corresponding to (2.3). The presence of cross product in (2.10) yields a linear system of equations with variable coefficients and non-symmetric structure. Often, GMRES is employed as the numerical solver. The BDF2 algorithm (2.10) treats both the gyromagentic and damping terms semi-implicitly, i.e., $\Delta m$ is treated implicitly while the prefactors are treated explicitly.

Since $|\boldsymbol{m}|=1$ , the strength of gyromagnetic term is controlled by $\epsilon\Delta\boldsymbol{m}+\boldsymbol{f}$ . Meanwhile, the strength of damping term is controlled by $\alpha(\epsilon\Delta\boldsymbol{m}+\boldsymbol{f})$ . Since $\alpha<1$ for many magnetic materials, it is reasonable to treat both the gyromagentic and damping terms semi-implicitly. However, for large $\alpha$ , it is possible to treat $\alpha\epsilon\Delta\boldsymbol{m}$ part in the damping term implicitly, and the gyromagnetic term and all remaining terms explicitly, as demonstrated in [41]. The stability and convergence of the scheme is proved under the condition $\alpha>3$ [50]. Starting with (2.5), the BDF2 algorithm in [41] becomes

\displaystyle\left\{\begin{array}[]{l}\dfrac{\frac{3}{2}\tilde{\boldsymbol{m}}_{h}^{n+2}-2\boldsymbol{m}_{h}^{n+1}+\frac{1}{2}\boldsymbol{m}_{h}^{n}}{k}=-\hat{\boldsymbol{m}}_{h}^{n+2}\times\left(\epsilon\Delta_{h}\hat{\boldsymbol{m}}_{h}^{n+2}+\hat{\boldsymbol{f}}_{h}^{n+2}\right)\\ \quad+\alpha\left(\epsilon\Delta_{h}\tilde{\boldsymbol{m}}_{h}^{n+2}+\hat{\boldsymbol{f}}_{h}^{n+2}\right)\\ \quad+\alpha\left(\epsilon\left|\nabla_{h}\hat{\boldsymbol{m}}_{h}^{n+2}\right|^{2}-\hat{\boldsymbol{m}}_{h}^{n+2}\cdot\hat{\boldsymbol{f}}_{h}^{n+2}\right)\hat{\boldsymbol{m}}_{h}^{n+2}\\ \boldsymbol{m}_{h}^{n+2}=\frac{\tilde{\boldsymbol{m}}_{h}^{n+2}}{\left|\tilde{\boldsymbol{m}}_{h}^{n+2}\right|}\end{array}\right.

(2.15)

where

\displaystyle\hat{\boldsymbol{m}}_{h}^{n+2}

\displaystyle=2\boldsymbol{m}_{h}^{n+1}-\boldsymbol{m}_{h}^{n},\quad\hat{\boldsymbol{f}}_{h}^{n+2}

\displaystyle=2\boldsymbol{f}_{h}^{n+1}-\boldsymbol{f}_{h}^{n}.

At each time step, only one linear system with constant coefficients and SPD structure needs to be solved with fast solvers, such as fast Fourier transform (FFT).

2.3 Implicit-explicit Runge-Kutta methods

For a given time-dependent nonlinear problem, the basic idea of implicit-explicit methods is to treat a dominant linear term implicitly and the remaining terms explicitly. The success of IMEX methods relies on the dominance of linear term that is implicitly treated [32]. While it is not the case for all problems, the introduction of an artificial term may help, such as the linear diffusion term introduced for the nonlinear diffusion equation [30]. For the LL equation, the linear diffusion term does not dominate the magnetization dynamics, and thus a direct application of IMEX method does not work. Motivated by this observation and the work in [30], we propose to introduce an artificial diffusion term in the LL equation which will be implicitly treated while all other terms are explicitly treated explicitly. RK methods are employed for the time discretization.

For ease of description, we first list the Butcher tableau for second-order and third-order RK methods in (2.16) and (2.17).

\begin{array}[]{c|ccc|ccc}0&0&0&0&0&0&0\\ 1/2&0&1/2&0&1/2&0&0\\ 1&1/2&0&1/2&0&1&0\\ \hline\cr&1/2&0&1/2&0&1&0\end{array}

(2.16)

\begin{array}[]{c|ccccc|ccccc}\rule{0.0pt}{20.0pt}0&0&0&0&0&0&0&0&0&0&0\\ 1/2&0&1/2&0&0&0&1/2&0&0&0&0\\ 2/3&0&1/6&1/2&0&0&11/18&1/18&0&0&0\\ 1/2&0&-1/2&1/2&1/2&0&5/6&-5/6&1/2&0&0\\ 1&0&3/2&-3/2&1/2&1/2&1/4&7/4&3/4&-7/4&0\\ \hline\cr&0&3/2&-3/2&1/2&1/2&1/4&7/4&3/4&-7/4&0\end{array}

(2.17)

We add an artificial damping term $\beta\Delta\boldsymbol{m}$ into (2.4) and rewrite the LL equation as

\displaystyle{{\boldsymbol{m}}_{t}}=\underbrace{-\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\emph{ {f}})-\alpha\boldsymbol{m}\times\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\emph{ {f}})-\beta{{\Delta}}\boldsymbol{m}}_{{N(t,\boldsymbol{m})}}+\underbrace{\beta{{\Delta}}\boldsymbol{m}}_{{L(t,\boldsymbol{m})}},

(2.18)

where the artificial term is denoted as $L(t,\boldsymbol{m})$ , and all the remaining terms are included in $N(t,\boldsymbol{m})$ .

Therefore, at time step $t_{n}$ , the corresponding marching algorithms in IMEX-RK2 and IMEX-RK3 methods are

\left\{\begin{array}[]{l}\boldsymbol{\tilde{m}}_{1}=\boldsymbol{{m}}_{n}\\ {\boldsymbol{\tilde{m}}_{2}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(L(t_{n+\frac{1}{2}},{\boldsymbol{\tilde{m}}_{2}})+N\left(t_{n},\boldsymbol{{\tilde{m}}}_{1}\right)\right)\\ {\boldsymbol{\tilde{m}}_{3}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(L\left(t_{n},\boldsymbol{{\tilde{m}}}_{1}\right)+L\left(t_{n+1},{\boldsymbol{\tilde{m}}_{3}}\right)+2N(t_{n+\frac{1}{2}},\boldsymbol{\tilde{m}}_{2})\right)\\ \boldsymbol{{m}}_{n+1}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(L\left(t_{n},\boldsymbol{{\tilde{m}}}_{1}\right)+L\left(t_{n+1},\boldsymbol{\tilde{m}}_{3}\right)+2N(t_{n+\frac{1}{2}},\boldsymbol{\tilde{m}}_{2})\right)\end{array}\right.,

(2.19)

and

\left\{\begin{array}[]{l}\boldsymbol{{\tilde{m}}}_{1}=\boldsymbol{{m}}_{n}\\ {\boldsymbol{\tilde{m}}_{2}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(L(t_{n+\frac{1}{2}},{\boldsymbol{\tilde{m}}_{2}})+N\left(t_{n},\boldsymbol{{\tilde{m}}}_{1}\right)\right)\\ {\boldsymbol{\tilde{m}}_{3}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{6}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})+\frac{k}{2}L({{t}_{n+\frac{2}{3}}},{\boldsymbol{\tilde{m}}_{3}})+\frac{11k}{18}N({{t}_{n}},\boldsymbol{{\tilde{m}}}_{1})+\frac{k}{18}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})\\ {\boldsymbol{\tilde{m}}_{4}}=\boldsymbol{{\tilde{m}}}_{1}-\frac{k}{2}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})+\frac{k}{2}L({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})+\frac{k}{2}L({{t}_{n+\frac{1}{2}}},{\boldsymbol{\tilde{m}}_{4}})\\ +\frac{5k}{6}N({{t}_{n}},\boldsymbol{{\tilde{m}}}_{1})-\frac{5k}{6}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})+\frac{k}{2}N({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})\\ {\boldsymbol{\tilde{m}}_{5}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{3k}{2}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})-\frac{3k}{2}L({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})+\frac{k}{2}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{4})+\frac{k}{2}L({{t}_{n+1}},{\boldsymbol{\tilde{m}}_{5}})\\ +\frac{k}{4}N({{t}_{n}},\boldsymbol{{\tilde{m}}}_{1})+\frac{7k}{4}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})+\frac{3k}{4}N({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})-\frac{7k}{4}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{4})\\ \boldsymbol{{m}}_{n+1}=\boldsymbol{{\tilde{m}}}_{1}+\frac{3k}{2}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})-\frac{3k}{2}L({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})+\frac{k}{2}L({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{4})+\frac{k}{2}L({{t}_{n+1}},\boldsymbol{\tilde{m}}_{5})\\ +\frac{k}{4}N({{t}_{n}},\boldsymbol{{\tilde{m}}}_{1})+\frac{7k}{4}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{2})+\frac{3k}{4}N({{t}_{n+\frac{2}{3}}},\boldsymbol{\tilde{m}}_{3})-\frac{7k}{4}N({{t}_{n+\frac{1}{2}}},\boldsymbol{\tilde{m}}_{4})\end{array}\right..

(2.20)

In addition, if ${\boldsymbol{f}}$ is time-independent, so that the rewritten LL equation (2.18) is autonomous, an alternate IMEX-RK2 numerical algorithm could be chosen:

\left\{\begin{array}[]{l}\boldsymbol{\tilde{m}}_{1}=\boldsymbol{{m}}_{n}\\ {\boldsymbol{\tilde{m}}_{2}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{4}L({\boldsymbol{\tilde{m}}_{2}})\\ {\boldsymbol{\tilde{m}}_{3}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}N(\boldsymbol{{\tilde{m}}}_{2})+\frac{k}{4}L({\boldsymbol{\tilde{m}}_{3}})\\ {\boldsymbol{\tilde{m}}_{4}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(N(\boldsymbol{{\tilde{m}}}_{2})+N(\boldsymbol{{\tilde{m}}}_{3})\right)+\frac{k}{3}\left(L({\boldsymbol{\tilde{m}}_{2}})+L({\boldsymbol{\tilde{m}}_{3}})+L({\boldsymbol{\tilde{m}}_{4}})\right)\\ \boldsymbol{{m}}_{n+1}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{3}\left(N(\boldsymbol{{\tilde{m}}}_{2})+N(\boldsymbol{{\tilde{m}}}_{3})+N(\boldsymbol{{\tilde{m}}}_{4})\right)\\ \qquad\qquad\qquad+\frac{k}{3}\left(L({\boldsymbol{\tilde{m}}_{2}})+L({\boldsymbol{\tilde{m}}_{3}})+L({\boldsymbol{\tilde{m}}_{4}})\right).\end{array}\right.,

(2.21)

This IMEX-RK2 algorithm contains four stages, with three intermediate numerical solutions at each time step. On the other hand, the explicit part satisfies the strong stability-preserving property [11, 17]; as a result, we denote it as the SSP-IMEX-RK2 scheme. In addition, this numerical method contains stronger diffusion coefficients, in comparison with the standard IMEX-RK2 algorithm (2.19), and this feature will greatly facilitate a theoretical justification of numerical stability and convergence analysis, as will be presented in Section 5.

Remark 2.1.

In the current work, the finite difference method is employed for the spatial discretization. It is worth mentioning that other spatial discretization, such as finite element method and discontinuous Galerkin method, can also be employed.

Remark 2.2.

Two linear systems with constant coefficients and SPD structure are solved in IMEX-RK2 (2.19). Although only one linear system solver is needed, with constant coefficients and SPD structure, the method in [41] only works in the large damping parameter case. Similarly, three linear system solvers are needed in the SSP-IMEX-RK2 method (2.21), and four linear system solvers are needed in the IMEX-RK3 method (2.20), with constant coefficients and SPD structure. It will be verified in Section 3 that IMEX-RK methods work for any damping parameter.

3 Numerical results

In this section, we provide a series of numerical experiments for IMEX-RK methods, including accuracy check, efficiency comparison, and stability test in terms of different $\beta$ values. Denote $\boldsymbol{m}_{e}$ the exact solution and $\boldsymbol{m}_{h}$ the numerical solution. To measure the error, we introduce the following notations in the discrete case.

Definition 3.1 ( $\ell^{2}$ inner product, ${{\left\|\cdot\right\|}_{2}}$ norm).

For grid functions ${{\boldsymbol{f}}_{h}}$ and ${{\boldsymbol{g}}_{h}}$ that take values on a uniform numerical grid, we define

\left\langle\boldsymbol{f}_{h},\boldsymbol{g}_{h}\right\rangle=h^{d}\sum_{\mathcal{I}\in\Lambda_{d}}\boldsymbol{f}_{\mathcal{I}}\cdot\boldsymbol{g}_{\mathcal{I}},

where $\Lambda_{d}$ is the set of grid point, and $\mathcal{I}$ is an index.

In turn, the ${{\left\|\cdot\right\|}_{2}}$ norm is defined as

\left\|\boldsymbol{f}_{h}\right\|_{2}=\left(\left\langle\boldsymbol{f}_{h},\boldsymbol{f}_{h}\right\rangle\right)^{1/2}.

In addition, the discrete ${H^{1}}$ norm is defined as

\left\|\boldsymbol{f}_{h}\right\|_{H^{1}}^{2}:=\left\|\boldsymbol{f}_{h}\right\|_{2}^{2}+\left\|\nabla_{h}\boldsymbol{f}_{h}\right\|_{2}^{2}.

Definition 3.2 ( ${{\left\|\cdot\right\|}_{\infty}}$ and ${{\left\|\cdot\right\|}_{p}}$ norms in the discrete sense).

For grid functions ${{\boldsymbol{f}}_{h}}$ that take values on a uniform numerical grid, we define

\left\|\boldsymbol{f}_{h}\right\|_{\infty}=\max_{\mathcal{I}\in\Lambda_{d}}\left\|\boldsymbol{f}_{\mathcal{I}}\right\|_{\infty},\quad\left\|\boldsymbol{f}_{h}\right\|_{p}=\left(h^{d}\sum_{I\in\Lambda_{d}}\left|\boldsymbol{f}_{\mathcal{I}}\right|^{p}\right)^{\frac{1}{p}},1\leq p<+\infty.

Lemma 1 (Summation by parts).

For any grid functions ${{\boldsymbol{f}}_{h}}$ and ${{\boldsymbol{g}}_{h}}$ , with ${{\boldsymbol{f}}_{h}}$ satisfying the discrete boundary condition (2.6), the following identity is valid:

\displaystyle\left\langle-\Delta_{h}\boldsymbol{f}_{h},\boldsymbol{g}_{h}\right\rangle=\left\langle\nabla_{h}\boldsymbol{f}_{h},\nabla_{h}\boldsymbol{g}_{h}\right\rangle.

(3.22)

The proof of the standard inverse inequality can be obtained in existing textbooks and references; we just cite the results here. In the sequel, for simplicity of notation, we will use the uniform constant $\mathcal{C}$ to denote all the controllable constants.

Lemma 2 (Inverse inequality).

[7, 8, 10] For each vector-valued grid function $\boldsymbol{f}_{h}\in X$ , we have

	$\displaystyle\\|{\boldsymbol{f}}_{h}\\|_{\infty}\leq\gamma{h}^{-1/2}(\\|{\boldsymbol{f}}_{h}\\|_{2}+\\|\nabla_{h}\boldsymbol{f}_{h}\\|_{2}),$		(3.23)
	$\displaystyle\\|\boldsymbol{f}_{h}\\|_{q}\leq\gamma{h}^{-(\frac{3}{2}-\frac{3}{q})}\\|\boldsymbol{f}_{h}\\|_{2},\,\,\,\forall 2<q\leq+\infty,$		(3.24)

in which constant $\gamma$ depends on $\Omega$ , as well as the form of the discrete $\|\cdot\|_{2}$ norm.

The following discrete Sobolev inequality has been derived in the existing works [20, 19], for the discrete grid function with periodic boundary condition; an extension to the discrete homogeneous Neumann boundary condition can be made in a similar fashion.

Lemma 3 (Discrete Sobolev inequality).

[20, 19] For a grid function $\boldsymbol{f}_{h}\in\boldsymbol{X}$ , we have the following discrete Sobolev inequality:

\displaystyle\|\boldsymbol{f}_{h}\|_{4}\leq\mathcal{C}\|\boldsymbol{f}_{h}\|_{2}^{\frac{1}{4}}\cdot\|\boldsymbol{f}_{h}\|_{H_{h}^{1}}^{\frac{3}{4}}\leq\mathcal{C}(\|\boldsymbol{f}_{h}\|_{2}+\|\boldsymbol{f}_{h}\|_{2}^{\frac{1}{4}}\cdot\|\nabla_{h}\boldsymbol{f}_{h}\|_{2}^{\frac{3}{4}}),

(3.25)

in which the positive constant $\mathcal{C}$ only depends on the domain $\Omega$ .

In the numerical simulation, we set $\epsilon=1$ and $\boldsymbol{f}=0$ in (2.5) for convenience. The 1-D exact solution is chosen to be

\boldsymbol{m}_{e}=(\cos(X)\sin t,\sin(X)\sin t,\cos t)^{T},

and the 3-D exact solution is set as

\boldsymbol{m}_{e}=(\cos(XYZ)\sin t,\sin(XYZ)\sin t,\cos t)^{T},

where $X=x^{2}(1-x)^{2},Y=y^{2}(1-y)^{2},Z=z^{2}(1-z)^{2}$ . It is easy to check that the homogeneous Neumann boundary condition (2.2) is satisfied. A forcing term $\boldsymbol{f}_{e}=\partial_{t}\boldsymbol{m}_{e}-\alpha\Delta\boldsymbol{m}_{e}-\alpha\left|\nabla\boldsymbol{m}_{e}\right|^{2}+\boldsymbol{m}_{e}\times$ $\Delta\boldsymbol{m}_{e}$ is included into the nonlinear part $N(t,\boldsymbol{m})$ .

3.1 Accuracy and efficiency test of IMEX-RK2

In the 1-D computation, we fix $h=1/5000$ and record the error in terms of $k$ in Table 2, and fix $k=1e-3/10000$ and record the error in terms of $h$ in Table 3. The second-order accuracy of IMEX-RK2 is observed in both space and time.

Table 2: Temporal accuracy of IMEX-RK2 in the 1-D computation (

\tau=1e-4,h=1/5000,\alpha=0.01

, and

\beta=5

$k$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$\tau/5$	$1.7011e-10$	$2.7861e-11$	$1.7582e-09$
$\tau/10$	$4.4130e-11$	$8.1126e-12$	$6.3812e-10$
$\tau/15$	$2.1077e-11$	$3.8709e-12$	$3.3541e-10$
$\tau/20$	$1.2028e-11$	$2.2727e-12$	$2.1471e-10$
$\tau/25$	$8.1095e-12$	$1.4985e-12$	$1.6165e-10$
order	$1.8930$	$1.8152$	$1.5001$

Table 3: Spatial accuracy of IMEX-RK2 in the 1-D computation (

k=1e-7,\alpha=0.01

, and

\beta=5

$h$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$1/50$	$2.7361e-09$	$8.3040e-10$	$4.5150e-07$
$1/100$	$7.3177e-10$	$2.1625e-10$	$1.1307e-07$
$1/150$	$3.2921e-10$	$9.6819e-11$	$5.0272e-08$
$1/200$	$1.8597e-10$	$5.4601e-11$	$2.8281e-08$
$1/250$	$1.1926e-10$	$3.4986e-11$	$1.8101e-08$
order	$1.9472$	$1.9681$	$1.9986$

In the 3-D computation, we fix $h=1/16$ and record the error in terms of $k$ in Table 4, and fix $k=1/10000$ and record the error in terms of $h$ in Table 5. Again, the second-order accuracy of IMEX-RK2 is observed in both space and time.

Table 4: Temporal accuracy of IMEX-RK2 in the 3-D computation (

h=1/16,\alpha=0.01

, and

\beta=5

$k$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$1/4$	$0.0022$	$0.0025$	$0.0025$
$1/6$	$0.0010$	$0.0011$	$0.0011$
$1/8$	$5.5504e-04$	$6.3857e-04$	$6.6930e-04$
$1/10$	$3.6163e-04$	$4.1354e-04$	$4.6041e-04$
order	$1.9773$	$1.9600$	$1.8594$

Table 5: Spatial accuracy of IMEX-RK2 in the 3-D computation (

k=1/10000,\alpha=0.01

, and

\beta=5

$h$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$1/3$	$1.4756e-04$	$1.7447e-04$	$2.5609e-04$
$1/5$	$5.2131e-05$	$6.1669e-05$	$9.3673e-05$
$1/7$	$2.6372e-05$	$3.1402e-05$	$4.7958e-05$
$1/9$	$1.5935e-05$	$1.8985e-05$	$2.9167e-05$
$1/11$	$1.0670e-05$	$1.2707e-05$	$1.9579e-05$
order	$2.0223$	$2.0155$	$1.9794$

In terms of the efficiency comparison, we plot the CPU time (in seconds) vs. the error $\left\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\|_{\infty}$ . Results of BDF2 and IMEX-RK2 are visualized in Fig. LABEL:time_1D and Fig. LABEL:space_1D for the 1-D case, and in Fig. LABEL:time_3D and Fig. LABEL:space_3D for the 3-D case.

By Fig. LABEL:rk2_cpu_time, IMEX-RK2 is superior to BDF2 in both the 1-D and 3-D computations.

3.2 Accuracy and efficiency test of IMEX-RK3

Using the same spatial resolution, we only test the temporal accuracy of IMEX-RK3 here. In the 1-D computation, we fix $k=0.01\times{{h}^{\frac{2}{3}}}$ and record the error in terms of $k$ in Table 6.

Table 6: Temporal accuracy of IMEX-RK3 in the 1-D computation (

k=0.01\times{{h}^{\frac{2}{3}}}

\alpha=0.01

, and

\beta=5

$k$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$1/208$	$0.0435$	$0.0508$	$0.0921$
$1/252$	$0.0244$	$0.0287$	$0.0527$
$1/292$	$0.0155$	$0.0184$	$0.0341$
$1/330$	$0.0109$	$0.0128$	$0.0236$
order	$3.0235$	$2.9930$	$2.8900$

In the 3-D computation, we fix $k=0.001\times{{h}^{\frac{2}{3}}}$ and record the error in terms of $k$ in Table 7. The third-order temporal accuracy is observed for IMEX-RK3 in time, in both the 1-D and 3-D compuations.

Table 7: Temporal accuracy of IMEX-RK3 in the 3-D compuation (

k=0.001\times{{h}^{\frac{2}{3}}}

\alpha=0.01

, and

\beta=5

$k$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$1/2080$	$1.4755e-04$	$1.7446e-04$	$2.5608e-04$
$1/2520$	$8.1182e-05$	$9.6730e-05$	$1.4581e-04$
$1/2924$	$5.2128e-05$	$6.1666e-05$	$9.3671e-05$
$1/3302$	$3.6028e-05$	$4.2767e-05$	$6.5140e-05$
order	$3.0460$	$3.0422$	$2.9621$

Since IMEX-RK2 and IMEX-RK3 only differ in the temporal discretization, we further plot the CPU time (in seconds) of IMEX-RK2 and IMEX-RK3 in terms of the temporal error in the 1-D and 3-D cases; see Fig. 3. These results indicate that IMEX-RK3 is more efficient than IMEX-RK2, and thus is more efficient than BDF2.

3.3 Accuracy test of SSP-IMEX-RK2

In the 1-D test, we keep $h=5\times{{10}^{1}}k$ and record the error in terms of $k$ in Table 8. In the 3-D computation, we fix $h=1\times{{10}^{3}}k$ and record the error in terms of $k$ in Table 9. Second-order accuracy of SSP-IMEX-RK2 has been observed in both the 1-D and 3-D cases. The accuracy curves are displayed in Fig. 4.

Table 8: Temporal accuracy of SSP-IMEX-RK2 in the 1-D computation (

\tau=2e-2,\alpha=0.01

, and

\beta=5

$k$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$\tau/3$	$5.4161e-05$	$1.1077e-05$	$0.0019$
$\tau/4$	$3.1570e-05$	$5.9625e-06$	$0.0012$
$\tau/5$	$2.0783e-05$	$3.6472e-06$	$7.7420e-04$
$\tau/6$	$1.4901e-05$	$2.4208e-06$	$5.4551e-04$
order	$1.8640$	$2.1927$	$1.8076$

Table 9: Temporal accuracy of SSP-IMEX-RK2 in the 3-D computation (

\tau=1e-3,\alpha=0.01

, and

\beta=5

$h$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{\infty}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{2}$	$\left\\|\boldsymbol{m}_{h}-\boldsymbol{m}_{e}\right\\|_{H^{1}}$
$\tau/8$	$1.7697e-10$	$6.6171e-11$	$4.5228e-08$
$\tau/10$	$1.1599e-10$	$4.5129e-11$	$2.9500e-08$
$\tau/12$	$8.1594e-11$	$3.3029e-11$	$2.0741e-08$
$\tau/14$	$6.0414e-11$	$2.5397e-11$	$1.5370e-08$
order	$1.9203$	$1.7115$	$1.9284$

3.4 Dependence on the damping parameter

There are a few numerical methods that only several linear systems with constant coefficients and SPD structure need to be solved at each time step, including the first-order-in-time GSPM [26, 47], the second-order-in-time method [41], and the current method. Numerically, the method in [41] only works when $\alpha>1$ , most magnetic materials correspond to $\alpha\ll 1$ . If $\alpha>1$ , we can set $\beta=\alpha$ and then apply the idea of IMEX-RK. Therefore, the current method works for a general damping parameter.

Next, we examine the performance of IMEX-RK on the choice of artificial damping parameter $\beta$ . The 3-D results are recorded in Table 10 and Table 11. On the basis of these results, it is clear that IMEX-RK methods work for general artificial damping parameters. The 1-D results are similar and are not listed here. Therefore, we can set $\beta>1$ if $\alpha\ll 1$ and $\beta=\alpha$ if $\alpha\geq 1$ numerically.

Table 10: 3-D errors of IMEX-RK2, with

h=500k

$\beta$	$k$	$\alpha=0.001$			$\alpha=0.01$
$\beta$	$k$	$L^{\infty}$	$L^{2}$	$H^{1}$	$L^{\infty}$	$L^{2}$	$H^{1}$
	$1/1000$	$3.59e-04$	$4.26e-04$	$5.44e-04$	$3.59e-04$	$4.26e-04$	$5.44e-04$
$1$	$1/2000$	$8.13e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/4000$	$2.02e-05$	$2.40e-05$	$3.68e-05$	$2.02e-05$	$2.40e-05$	$3.68e-05$
order		$2.0758$	$2.0749$	$1.9429$	$2.0758$	$2.0749$	$1.9429$
	$1/1000$	$3.59e-04$	$4.26e-04$	$5.44e-04$	$3.59e-04$	$4.26e-04$	$5.44e-04$
$3$	$1/2000$	$8.13e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/4000$	$2.02e-05$	$2.40e-05$	$3.68e-05$	$2.02e-05$	$2.40e-05$	$3.68e-05$
order		$2.0758$	$2.0749$	$1.9429$	$2.0758$	$2.0749$	$1.9429$
	$1/1000$	$3.59e-04$	$4.26e-04$	$5.44e-04$	$3.59e-04$	$4.26e-04$	$5.44e-04$
$4$	$1/2000$	$8.13e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/4000$	$2.03e-05$	$2.40e-05$	$3.68e-05$	$2.02e-05$	$2.40e-05$	$3.68e-05$
order		$2.0722$	$2.0749$	$1.9429$	$2.0758$	$2.0749$	$1.9429$

Table 11: 3-D error of IMEX-RK3, with

k=0.001{{h}^{\frac{2}{3}}}

$\beta$	$k$	$\alpha=0.001$			$\alpha=0.01$
$\beta$	$k$	$L^{\infty}$	$L^{2}$	$H^{1}$	$L^{\infty}$	$L^{2}$	$H^{1}$
	$1/2080$	$1.48e-04$	$1.74e-04$	$2.56e-04$	$1.48e-04$	$1.74e-04$	$2.56e-04$
$1$	$1/2520$	$8.12e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/3302$	$3.61e-05$	$4.28e-05$	$6.51e-05$	$3.60e-05$	$4.28e-05$	$6.51e-05$
order		$3.0494$	$3.0335$	$2.9644$	$3.0557$	$3.0335$	$2.9644$
	$1/2080$	$1.48e-04$	$1.74e-04$	$2.56e-04$	$1.48e-04$	$1.74e-04$	$2.56e-04$
$3$	$1/2520$	$8.12e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/3302$	$3.61e-05$	$4.28e-05$	$6.51e-05$	$3.60e-05$	$4.28e-05$	$6.51e-05$
order		$3.0494$	$3.0335$	$2.9644$	$3.0557$	$3.0335$	$2.9644$
	$1/2080$	$1.48e-04$	$1.74e-04$	$2.56e-04$	$1.48e-04$	$1.74e-04$	$2.56e-04$
$4$	$1/2520$	$8.12e-05$	$9.67e-05$	$1.46e-04$	$8.12e-05$	$9.67e-05$	$1.46e-04$
	$1/3302$	$3.61e-05$	$4.28e-05$	$6.51e-05$	$3.60e-05$	$4.28e-05$	$6.51e-05$
order		$3.0494$	$3.0335$	$2.9644$	$3.0557$	$3.0335$	$2.9644$

Remark 3.3.

We have demonstrated the high-order-in-time accuracy of IMEX-RK, including second-order and third-order ones. It is worth noting that even higher-order-in-time IMEX-RK methods for the LL equation can be designed; see the $s$ -level method in [51] for a general purpose. This method does not necessarily have the $s$ -order accuracy when $s>4$ in [51]. For instance, the $s$ -level method has at most $(s-1)$ -order accuracy when $s=5,6,7$ .

4 Micromagnetics simulations

In this section, we apply IMEX-RK2 and IMEX-RK3 to conduct micromagnetics simulations, including different equilibrium structures and a benchmark problem from NIST [42] to examine the performance of our proposed methods in the real-world applications.

4.1 Equilibrium states

We use a spatial resolution $64\times 128\times 1$ on a $1\times 2\times 0.02$ $\mu{{m}^{3}}$ thin-film element with $\alpha=0.1$ , a temporal step size $k=1$ picosecond (ps), and set $\beta=3$ in IMEX-RK methods. In the absence of an external field, multiple metastable states are often observed in ferromagnets, both experimentally and numerically [33, 34].

Starting with three different initial magnetization distributions, we obtain three equilibrium states in Fig. LABEL:metastable_states, including Landau state, C-state, and S-state. The arrow denotes the first two components of the magnetization vector and the color denotes the angle between them. It is clearly that IMEX-RK2,3 produce consistent results.

4.2 Benchmark peoblem from NIST

To investigate the dynamical performance of IMEX-RK methods, we simulate a benchmark problem proposed by the Micromagnetic Modeling Activity Group (muMag) from NIST. A positive external field of strength ${{H}_{0}}={{\mu}_{0}}{{H}_{e}}$ in the unit of $mT$ is applied. The magnetization is able to reach a steady state. Once this state is reached, the applied external field is reduced by a certain amount and the material sample is allowed to reach another steady state again. Repeat the process until the hysteresis system attains a negative external field of strength ${{H}_{0}}$ . This process is then implemented in reverse, increasing the field in small steps until the initial applied external field is reached. Afterward, we can plot the average magnetization at steady states as a function of the external field strength during the hysteresis loop. The stopping criterion for a steady state is that the relative change of the total energy is less than ${{10}^{-9}}$ . For comparison with the available code $mo96a$ of the first standard problem from NIST, we set $100\times 50\times 1$ spatial resolution with mesh size $20\times 20\times 20$ $\rm n{{m}^{3}}$ and the canting angle $+{{1}^{{}^{\circ}}}$ of applied field from the nominal axis. The initial state is uniform and $\left[\rm-50mT,+50mT\right]$ is split into 200 steps for both $x$ -loop and $y$ -loop.

Hysteresis loops generated by $mo96a$ are displayed in Fig. 6a and Fig. 6b when the applied field is approximately parallel to the $y$ -(long) axis and the $x$ -(short) axis. The average remanent magnetization in reduced units is given by $(-1.5120\times{{10}^{-1}},8.6964\times{{10}^{-1}},0)$ for the $y$ -loop and $(1.5257\times{{10}^{-1}},8.6870\times{{10}^{-1}},0)$ for the $x$ -loop. The coercive fields are 4.8871 $\rm mT$ in Fig. 6a and 2.5253 $\rm mT$ in Fig. 6b, respectively. Hysteresis loops generated by IMEX-RK2 method are presented in Fig. 6c and Fig. 6d when the applied field is approximately parallel to the long axis and the short axis, respectively. The average remanent magnetization in reduced units is $(-1.6099\times{{10}^{-1}},8.6096\times{{10}^{-1}},4.2423\times{{10}^{-7}})$ for the $y$ -loop and $(-1.4274\times{{10}^{-1}},8.6656\times{{10}^{-1}},1.5753\times{{10}^{-7}})$ for the $x$ -loop. The coercive fields are $5.4688~{}(\pm 0.7)~{}\rm mT$ in Fig. 6c and $2.7188~{}(\pm 0.2)~{}\rm mT$ in Fig. 6d. Similarly, the IMEX-RK3 results are presented in the bottom row of Fig. 6. The average remanent magnetization in reduced units is $(-1.6107\times{{10}^{-1}},8.6102\times{{10}^{-1}},9.5492\times{{10}^{-8}})$ for the $y$ -loop and $(-1.4292\times{{10}^{-1}},8.6659\times{{10}^{-1}},4.5836\times{{10}^{-9}})$ for the $x$ -loop. The coercive fields are $5.4688~{}(\pm 0.7)~{}\rm mT$ in Fig. 6e and $2.7188~{}(\pm 0.2)~{}\rm mT$ in Fig. 6f. Based on these results, we conclude that IMEX-RK methods work well for the benchmark problems from NIST, both qualitatively and quantitatively.

5 Numerical stability and convergence analysis for the proposed SSP-IMEX-RK2 scheme

A theoretical analysis for the proposed IMEX-RK numerical schemes turns out to be highly challenging, due to the multi-stage nature, as well as the highly complicated nonlinear terms in the vector form. For simplicity, we focus on the SSP-IMEX-RK2 numerical algorithm (2.21). In the first step, a numerical stability is stablished for the linear part, i.e., in the simple case with only linear diffusion part (the term $\beta\Delta\boldsymbol{m}$ ) taken into consideration. Afterward, we provide a convergence analysis of the SSP-IMEX-RK2 scheme (2.21) for a simplified nonlinear model of LL equation, in which only the damping term is considered, while the gyromagnetic term is skipped.

5.1 Linear stability estimate for the SSP-IMEX-RK2 scheme (2.21)

In the simple case with only linear diffusion term is considered, we denote $L_{h}=\beta\Delta_{h}$ . The SSP-IMEX-RK2 scheme (2.21) is simplified as

\left\{\begin{array}[]{l}\boldsymbol{\tilde{m}}_{1}=\boldsymbol{{m}}_{n}\\ {\boldsymbol{\tilde{m}}_{2}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{4}L_{h}({\boldsymbol{\tilde{m}}_{2}})\\ {\boldsymbol{\tilde{m}}_{3}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{4}L_{h}({\boldsymbol{\tilde{m}}_{3}})\\ {\boldsymbol{\tilde{m}}_{4}}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{3}\left(L_{h}({\boldsymbol{\tilde{m}}_{2}})+L_{h}({\boldsymbol{\tilde{m}}_{3}})+L_{h}({\boldsymbol{\tilde{m}}_{4}})\right)\\ \boldsymbol{{m}}_{n+1}=\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{3}\left(L_{h}({\boldsymbol{\tilde{m}}_{2}})+L_{h}({\boldsymbol{\tilde{m}}_{3}})+L_{h}({\boldsymbol{\tilde{m}}_{4}})\right)\end{array}\right..

(5.26)

For the convenience of the stability analysis, numerical system could be rewritten as

	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\boldsymbol{m}}_{2},$		(5.27)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{3}-\boldsymbol{m}_{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\boldsymbol{m}}_{3},$		(5.28)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{4}-\boldsymbol{m}_{n}}{k}=\frac{\beta}{3}\Delta_{h}(\tilde{\boldsymbol{m}}_{2}+\tilde{\boldsymbol{m}}_{3}+\tilde{\boldsymbol{m}}_{4}),$		(5.29)
	$\displaystyle\boldsymbol{m}_{n+1}=\tilde{\boldsymbol{m}}_{4}.$		(5.30)

Moreover, to reveal the numerical stability of this Runge-Kutta style algorithm, we subtract (5.27) from (5.28), (5.28) from (5.29), and arrive at the following equivalent numerical system:

	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\boldsymbol{m}}_{2},$		(5.31)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}}{k}=\frac{\beta}{4}\Delta_{h}(\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}),$		(5.32)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{4}-\tilde{\boldsymbol{m}}_{3}}{k}=\beta\Delta_{h}(\frac{1}{3}\tilde{\boldsymbol{m}}_{2}+\frac{1}{12}\tilde{\boldsymbol{m}}_{3}+\frac{1}{3}\tilde{\boldsymbol{m}}_{4}),$		(5.33)
	$\displaystyle\boldsymbol{m}_{n+1}=\tilde{\boldsymbol{m}}_{4}.$		(5.34)

Taking a discrete inner product with (5.31) by $2\tilde{\boldsymbol{m}}_{2}$ gives

\|\tilde{\boldsymbol{m}}_{2}\|_{2}^{2}-\|\boldsymbol{m}_{n}\|_{2}^{2}+\|\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}\|_{2}^{2}+\frac{\beta}{2}k\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\|_{2}^{2}=0,

(5.35)

in which the summation-by-parts formula (3.22) has been applied. Similarly, taking a discrete inner product with (5.32) by $2\tilde{\boldsymbol{m}}_{3}$ , with (5.33) by $2\tilde{\boldsymbol{m}}_{4}$ , leads to

	$\displaystyle\\|\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}-\\|\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\frac{\beta}{2}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}=\frac{\beta}{2}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{3}\rangle,$		(5.36)
	$\displaystyle\\|\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}-\\|\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{4}-\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\frac{2\beta}{3}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}$
	$\displaystyle\qquad\qquad=-\frac{2\beta}{3}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle-\frac{\beta}{6}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{3},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle.$		(5.37)

In turn, a combination of (5.35) and (5.37) indicates that

	$\displaystyle\\|\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}-\\|\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{4}-\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}$
	$\displaystyle\qquad+\frac{\beta}{2}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\frac{\beta}{2}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\frac{2\beta}{3}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}$
	$\displaystyle\qquad=\frac{\beta}{2}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{3}\rangle-\frac{2\beta}{3}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle-\frac{\beta}{6}k\langle\nabla_{h}\tilde{\boldsymbol{m}}_{3},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle.$		(5.38)

Meanwhile, a careful application of Cauchy inequality implies that

	$\displaystyle\frac{1}{2}\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{3}\rangle\leq\frac{1}{4}(\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}),$		(5.39)
	$\displaystyle-\frac{2}{3}\langle\nabla_{h}\tilde{\boldsymbol{m}}_{2},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle\leq\frac{2}{9}\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\frac{1}{2}\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2},$		(5.40)
	$\displaystyle-\frac{1}{6}\langle\nabla_{h}\tilde{\boldsymbol{m}}_{3},\nabla_{h}\tilde{\boldsymbol{m}}_{4}\rangle\leq\frac{1}{12}\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\frac{1}{12}\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}.$		(5.41)

Therefore, a substitution of these estimates into (5.38) yields

	$\displaystyle\\|\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}-\\|\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{4}-\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}$
	$\displaystyle\qquad+\frac{\beta}{36}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\frac{\beta}{6}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\frac{\beta}{12}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{2}^{2}\leq 0.$		(5.42)

By the fact that $\boldsymbol{m}_{n+1}=\tilde{\boldsymbol{m}}_{4}$ , we obtain the linear stability estimate for the SSP-IMEX-RK2 scheme:

	$\displaystyle\\|\boldsymbol{m}_{n+1}\\|_{2}^{2}-\\|\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\\|\boldsymbol{m}_{n+1}-\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}$
	$\displaystyle\qquad+\frac{\beta}{36}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{2}^{2}+\frac{\beta}{6}k\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{2}^{2}+\frac{\beta}{12}k\\|\nabla_{h}\boldsymbol{m}_{n+1}\\|_{2}^{2}\leq 0,$		(5.43)

which in turn gives the $\ell^{\infty}(0,T:\ell^{2})\cap\ell^{2}(0,T;H_{h}^{1})$ bound of the numerical solution, if only is the linear diffusion part is considered:

\|\boldsymbol{m}_{n+1}\|_{2}+\Big{(}\frac{\beta}{12}k\sum_{j=1}^{n+1}\|\nabla_{h}\boldsymbol{m}_{j}\|_{2}^{2}\Big{)}^{\frac{1}{2}}\leq\|\boldsymbol{m}_{0}\|_{2}.

(5.44)

Remark 5.1.

In comparison with the standard three-stage IMEX-RK2 method (2.19), the SSP-IMEX-RK2 algorithm (2.21) contains four stages, with three intermediate numerical solutions, so that more computations are needed at each time step. Meanwhile, the stability analysis in this section reveals that, this numerical algorithm contains stronger diffusion coefficients than the standard IMEX-RK2 algorithm. In more details, the diffusion part in the standard IMEX-RK2 method (2.19) essentially correspond to the Crank-Nicolson approximation, which may face a serious theoretical difficulty in the nonlinear analysis, while additional diffusion terms appear in the stability estimate (5.42) for the SSP-IMEX-RK2 algorithm (2.21). This subtle fact will greatly facilitate the convergence analysis in the next subsection.

5.2 Convergence analysis of the SSP-IMEX-RK2 scheme (2.21) for a simplified nonlinear LL equation

We consider a simplified nonlinear LL equation (2.4) in this subsection, in which only the damping term is included, while the gyromagnetic term is skipped for simplicity:

\displaystyle{{\boldsymbol{m}}_{t}}=-\alpha\boldsymbol{m}\times(\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\emph{ {f}})).

(5.45)

For a vector function $\boldsymbol{m}$ with $|\boldsymbol{m}|\equiv 1$ , the following identity is recalled

-\boldsymbol{m}\times(\boldsymbol{m}\times\Delta\boldsymbol{m})=\Delta\boldsymbol{m}+|\nabla\boldsymbol{m}|^{2}\boldsymbol{m}.

(5.46)

In turn, by taking $\beta=\alpha\epsilon$ , the nonlinear term $N(\boldsymbol{m})$ could be rewritten as

$\displaystyle N(\boldsymbol{m})=$	$\displaystyle-\alpha\boldsymbol{m}\times(\boldsymbol{m}\times(\epsilon\Delta\boldsymbol{m}+\boldsymbol{f}))-\beta\Delta\boldsymbol{m}$	(5.47)
$\displaystyle=$	$\displaystyle\beta(\Delta\boldsymbol{m}+\|\nabla\boldsymbol{m}\|^{2}\boldsymbol{m})-\alpha\boldsymbol{m}\times(\boldsymbol{m}\times\boldsymbol{f})-\beta\Delta\boldsymbol{m}$
$\displaystyle=$	$\displaystyle\beta\|\nabla\boldsymbol{m}\|^{2}\boldsymbol{m}-\alpha\boldsymbol{m}\times(\boldsymbol{m}\times\boldsymbol{f}).$

Subsequently, the discrete form of the nonlinear term becomes

N_{h}(\boldsymbol{m})=\beta|{\mathcal{A}}_{h}\nabla_{h}\boldsymbol{m}|^{2}\boldsymbol{m}-\alpha\boldsymbol{m}\times(\boldsymbol{m}\times\boldsymbol{f}),

(5.48)

in which ${\mathcal{A}}_{h}\nabla_{h}$ (second approximation to the gradient operator) is an average gradient operator defined for the gird function $\boldsymbol{m}=(u_{h},v_{h},w_{h})^{T}\in\boldsymbol{X}$ as $\mathcal{A}_{h}\nabla_{h}\boldsymbol{m}_{h}=\nabla_{h}\mathcal{A}_{h}\boldsymbol{m}_{h}$ and $\mathcal{A}_{h}\boldsymbol{m}=(\mathcal{A}_{x}u_{h},\mathcal{A}_{y}v_{h},\mathcal{A}_{z}w_{h})$ :

\mathcal{A}_{x}u_{i,j,\ell}=\frac{u_{i,j,\ell}+u_{i-1,j,\ell}}{2},\,\mathcal{A}_{y}v_{i,j,\ell}=\frac{v_{i,j,\ell}+v_{i,j-1,\ell}}{2},\,\mathcal{A}_{z}w_{i,j,\ell}=\frac{w_{i,j,\ell}+w_{i,j,\ell-1}}{2}.

As a result, the SSP-IMEX-RK2 numerical algorithm is formulated as

$\displaystyle\boldsymbol{\tilde{m}}_{1}=$	$\displaystyle\boldsymbol{{m}}_{n},$	(5.49)
$\displaystyle{\boldsymbol{\tilde{m}}_{2}}=$	$\displaystyle\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{4}L_{h}({\boldsymbol{\tilde{m}}_{2}}),$
$\displaystyle{\boldsymbol{\tilde{m}}_{3}}=$	$\displaystyle\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}N_{h}(\boldsymbol{{\tilde{m}}}_{2})+\frac{k}{4}L_{h}({\boldsymbol{\tilde{m}}_{3}}),$
$\displaystyle{\boldsymbol{\tilde{m}}_{4}}=$	$\displaystyle\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{2}\left(N_{h}(\boldsymbol{{\tilde{m}}}_{2})+N_{h}(\boldsymbol{{\tilde{m}}}_{3})\right)+\frac{k}{3}\left(L_{h}({\boldsymbol{\tilde{m}}_{2}})+L_{h}({\boldsymbol{\tilde{m}}_{3}})+L_{h}({\boldsymbol{\tilde{m}}_{4}})\right),$
$\displaystyle\boldsymbol{{m}}_{n+1}=$	$\displaystyle\boldsymbol{{\tilde{m}}}_{1}+\frac{k}{3}\left(N_{h}(\boldsymbol{{\tilde{m}}}_{2})+N_{h}(\boldsymbol{{\tilde{m}}}_{3})+N_{h}(\boldsymbol{{\tilde{m}}}_{4})\right)$
	$\displaystyle+\frac{k}{3}\left(L_{h}({\boldsymbol{\tilde{m}}_{2}})+L_{h}({\boldsymbol{\tilde{m}}_{3}})+L_{h}({\boldsymbol{\tilde{m}}_{4}})\right).$

For the convenience of the Runge-Kutta analysis, this numerical system could be equivalently rewritten as

	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{2}-\boldsymbol{m}_{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\boldsymbol{m}}_{2},$		(5.50)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}}{k}=\frac{1}{2}N_{h}(\tilde{\boldsymbol{m}}_{2})+\frac{\beta}{4}\Delta_{h}(\tilde{\boldsymbol{m}}_{3}-\tilde{\boldsymbol{m}}_{2}),$		(5.51)
	$\displaystyle\frac{\tilde{\boldsymbol{m}}_{4}-\tilde{\boldsymbol{m}}_{3}}{k}=\frac{1}{2}N_{h}(\tilde{\boldsymbol{m}}_{3})+\beta\Delta_{h}(\frac{1}{3}\tilde{\boldsymbol{m}}_{2}+\frac{1}{12}\tilde{\boldsymbol{m}}_{3}+\frac{1}{3}\tilde{\boldsymbol{m}}_{4}),$		(5.52)
	$\displaystyle\frac{\boldsymbol{m}_{n+1}-\tilde{\boldsymbol{m}}_{4}}{k}=-\frac{1}{6}(N_{h}(\tilde{\boldsymbol{m}}_{2})+N_{h}(\tilde{\boldsymbol{m}}_{3}))+\frac{1}{3}N_{h}(\tilde{\boldsymbol{m}}_{4}).$		(5.53)

Denote $\Phi$ as the exact solution to the LL equation (5.45), with the regularity

\Phi\in{\mathcal{R}}=C^{3}([0,T];[C^{0}(\bar{\Omega})]^{3})\cap C^{2}([0,T];[C^{2}(\bar{\Omega})]^{3})\cap L^{\infty}([0,T];[C^{4}(\bar{\Omega})]^{3}).

(5.54)

The main theoretical result is stated in the following theorem.

Theorem 4.

Assume that the exact solution $\Phi$ of (5.45) has the regularity ${\mathcal{R}}$ . Denote ${\boldsymbol{m}}^{n}$ ( $n\geq 0$ ) as the numerical solution obtained from (5.49), or equivalently (5.50)-(5.53), with the initial error satisfying $\|\mathcal{P}_{h}\Phi(\cdot,t_{0})-\boldsymbol{m}_{0}\|_{2}+\|\nabla_{h}(\mathcal{P}_{h}\Phi(\cdot,t_{0})-\boldsymbol{m}_{0})\|_{2}=\mathcal{O}(h^{2})$ . In addition, a linear refinement assumption is made for the time step size: $C_{1}h\leq k\leq C_{2}h$ . Then the following convergence result holds for $1\leq n\leq\left\lfloor\frac{T}{k}\right\rfloor$ as $k,h\to 0^{+}$ :

\displaystyle\|\Phi(\cdot,t_{n})-\boldsymbol{m}^{n}\|_{2}

\displaystyle\leq\mathcal{C}(k^{2}+h^{2}),

(5.55)

in which the constant $\mathcal{C}>0$ is independent of $k$ and $h$ .

Proof 5.2.

(Proof of Theorem 4.) Around the boundary section $z=0$ , we set $\hat{z}_{0}=-\frac{1}{2}h$ , $\hat{z}_{1}=\frac{1}{2}h$ , and we can extend the profile $\Phi$ to the numerical “ghost” points, according to the extrapolation formula (2.6):

\Phi_{i,j,0}=\Phi_{i,j,1},\quad\Phi_{i,j,N+1}=\Phi_{i,j,N},

(5.56)

and the extrapolation for other boundaries can be formulated in the same manner. The proof of such an extrapolation yields a higher order $\mathcal{O}(h^{5})$ approximation, instead of the standard $\mathcal{O}(h^{3})$ accuracy. See the related derivation in [50], as well as the related consistency analysis works [40, 24, 23, 25].

Given the exact solution $\Phi$ , we denote $\Phi^{n}=\Phi(\cdot,t^{n})$ . To facilitate the Runge-Kutta analysis, three more intermediate approximate solutions are constructed at each time step, following the same algorithm as in (5.49):

$\displaystyle\tilde{\Phi}^{n,(2)}=$	$\displaystyle\Phi^{n}+\frac{\beta k}{4}\Delta_{h}\tilde{\Phi}^{n,(2)},$	(5.57)
$\displaystyle\tilde{\Phi}^{n,(3)}=$	$\displaystyle\tilde{\Phi}^{n}+\frac{k}{2}N_{h}(\tilde{\Phi}^{n,(2)})+\frac{\beta k}{4}\Delta_{h}\tilde{\Phi}^{n,(3)},$	(5.58)
$\displaystyle\tilde{\Phi}^{n,(4)}=$	$\displaystyle\tilde{\Phi}^{n}+\frac{k}{2}(N_{h}(\tilde{\Phi}^{n,(2)})+N_{h}(\tilde{\Phi}^{n,(3)}))$
	$\displaystyle+\frac{\beta k}{3}\Delta_{h}(\tilde{\Phi}^{n,(2)}+\tilde{\Phi}^{n,(3)}+\tilde{\Phi}^{n,(4)}),$	(5.59)

in which the homogeneous discrete Neumann boundary condition (similar to (2.6), (5.56)) is imposed for $\tilde{\Phi}^{n,(j)}$ , $j=2,3,4$ . Subsequently, a careful Taylor expansion (associated with the SSP-IMEX-RK schemes) reveals the following consistency estimate, for the exact solution at the next time step:

	$\displaystyle\Phi^{n+1}=$	$\displaystyle\tilde{\Phi}^{n}+\frac{k}{3}(N_{h}(\tilde{\Phi}^{n,(2)})+N_{h}(\tilde{\Phi}^{n,(3)})+N_{h}(\tilde{\Phi}^{n,(4)}))$
		$\displaystyle+\frac{\beta k}{3}\Delta_{h}(\tilde{\Phi}^{n,(2)}+\tilde{\Phi}^{n,(3)}+\tilde{\Phi}^{n,(4)})+k\tau_{0}^{n},\quad\\|\tau_{0}^{n}\\|_{2}\leq{\mathcal{C}}(k^{2}+h^{2}).$		(5.60)

Of course, with a similar transformation as in (5.50)-(5.53), the exact solution $\Phi^{n}$ , $\Phi^{n+1}$ and the constructed profiles $\tilde{\Phi}^{n,(j)}$ ( $j=2,3,4$ ) satisfy the following numerical system:

	$\displaystyle\frac{\tilde{\Phi}^{n,(2)}-\Phi^{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\Phi}^{n,(2)},$		(5.61)
	$\displaystyle\frac{\tilde{\Phi}^{n,(3)}-\tilde{\Phi}^{n,(2)}}{k}=\frac{1}{2}N_{h}(\tilde{\Phi}^{n,(2)})+\frac{\beta}{4}\Delta_{h}(\tilde{\Phi}^{n,(3)}-\tilde{\Phi}^{n,(2)}),$		(5.62)
	$\displaystyle\frac{\tilde{\Phi}^{n,(4)}-\tilde{\Phi}^{n,(3)}}{k}=\frac{1}{2}N_{h}(\tilde{\Phi}^{n,(3)})+\beta\Delta_{h}(\frac{1}{3}\tilde{\Phi}^{n,(2)}+\frac{1}{12}\tilde{\Phi}^{n,(3)}+\frac{1}{3}\tilde{\Phi}^{n,(4)}),$		(5.63)
	$\displaystyle\frac{\Phi^{n+1}-\tilde{\Phi}^{n,(4)}}{k}=-\frac{1}{6}(N_{h}(\tilde{\Phi}^{n,(2)})+N_{h}(\tilde{\Phi}^{n,(3)}))+\frac{1}{3}N_{h}(\tilde{\Phi}^{n,(4)})+\tau_{0}^{n}.$		(5.64)

It is clear that the constructed profiles $\tilde{\Phi}^{n,(j)}$ ( $j=2,3,4$ ) only depend on the exact solution $\Phi^{n}$ , and the consistency estimate indicates that

\|\tilde{\Phi}^{n,(j)}\|_{\infty}\leq\frac{9}{8},\quad\|\nabla_{h}\tilde{\Phi}^{n,(j)}\|_{\infty}\leq{\mathcal{C}}^{*},\quad j=2,3,4.

(5.65)

The following numerical error functions are defined:

		$\displaystyle\boldsymbol{e}^{k}=\Phi^{k}-\boldsymbol{m}_{k},\,\,\,k=n,n+1,$		(5.66)
		$\displaystyle\tilde{\boldsymbol{e}}^{n,(j)}=\tilde{\Phi}^{n,(j)}-\tilde{\boldsymbol{m}}_{j},\,\,\,j=2,3,4,$		(5.66)

at a point-wise level. In addition, the following nonlinear error terms are introduced:

{\mathcal{N}LE}^{n,(j)}=N_{h}(\tilde{\Phi}^{n,(j)})-N_{h}(\tilde{\boldsymbol{m}}_{j}),\quad j=2,3,4.

(5.67)

Therefore, subtracting the numerical scheme (5.50)-(5.53) from the consistency estimate (5.61)-(5.64) yields

	$\displaystyle\frac{\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}}{k}=\frac{\beta}{4}\Delta_{h}\tilde{\boldsymbol{e}}^{n,(2)},$		(5.68)
	$\displaystyle\frac{\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}}{k}=\frac{1}{2}{\mathcal{N}LE}^{n,(2)}+\frac{\beta}{4}\Delta_{h}(\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}),$		(5.69)
	$\displaystyle\frac{\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}}{k}=\frac{1}{2}{\mathcal{N}LE}^{n,(3)}+\beta\Delta_{h}(\frac{1}{3}\tilde{\boldsymbol{e}}^{n,(2)}+\frac{1}{12}\tilde{\boldsymbol{e}}^{n,(3)}+\frac{1}{3}\tilde{\boldsymbol{e}}^{n,(4)}),$		(5.70)
	$\displaystyle\frac{\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}}{k}=-\frac{1}{6}({\mathcal{N}LE}^{n,(2)}+{\mathcal{N}LE}^{n,(3)})+\frac{1}{3}{\mathcal{N}LE}^{n,(4)}+\tau_{0}^{n}.$		(5.71)

In addition, the discrete homogeneous Neumann boundary condition (2.6) is satisfied for both $\boldsymbol{e}^{n+1}$ , $\boldsymbol{e}^{n}$ , as well as the intermediate error functions $\tilde{\boldsymbol{e}}^{n,(j)}$ , $j=2,3,4$ .

To facilitate the convergence proof, the following functional bound of the nonlinear error terms is needed.

Lemma 5.

Under the regularity estimate (5.65) for the constructed profiles, and the following bound for the numerical solution in the IMEX-RK stages

\|\tilde{\boldsymbol{m}}_{j}\|_{\infty}\leq\frac{5}{4},\quad\|\nabla_{h}\tilde{\boldsymbol{m}}_{j}\|_{8}\leq\tilde{\mathcal{C}}:={\mathcal{C}}^{*}+1,\quad j=2,3,4.

(5.72)

we have an $\|\cdot\|_{\frac{8}{5}}$ estimate for the nonlinear error terms:

\|{\mathcal{N}LE}^{n,(j)}\|_{\frac{8}{5}}\leq\tilde{M}(\|\tilde{\boldsymbol{e}}^{n,(j)}\|_{2}+\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\|_{2}),\quad j=2,3,4,

(5.73)

in which $\tilde{M}$ only depends on $\alpha$ , $\beta$ , ${\mathcal{C}}^{*}$ , $\tilde{\mathcal{C}}$ , and the external force term $\boldsymbol{f}$ .

Proof 5.3.

(Proof of Lemma 5.) A careful expansion of the nonlinear terms $N_{h}(\tilde{\Phi}^{n,(j)})$ and $N_{h}(\boldsymbol{m}_{j})$ gives

$\displaystyle{\mathcal{N}LE}^{n,(j)}=$	$\displaystyle N_{h}(\tilde{\Phi}^{n,(j)})-N_{h}(\tilde{\boldsymbol{m}}_{j})$
$\displaystyle=$	$\displaystyle\beta\|{\mathcal{A}}_{h}\nabla_{h}\tilde{\Phi}^{n,(j)}\|^{2}\tilde{\boldsymbol{e}}^{n,(j)}+\beta\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\tilde{\Phi}^{n,(j)}+\tilde{\boldsymbol{m}}_{j})\cdot{\mathcal{A}}_{h}\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\Big{)}\tilde{\boldsymbol{m}}_{j}$
	$\displaystyle-\alpha\tilde{\boldsymbol{m}}_{j}\times(\tilde{\boldsymbol{e}}^{n,(j)}\times\boldsymbol{f})-\alpha\tilde{\boldsymbol{e}}^{n,(j)}\times(\tilde{\Phi}^{n,(j)}\times\boldsymbol{f}).$	(5.74)

In turn, an application of discrete Hölder inequality leads to

	$\displaystyle\Big{\\|}\beta\|{\mathcal{A}}_{h}\nabla_{h}\tilde{\Phi}^{n,(j)}\|^{2}\tilde{\boldsymbol{e}}^{n,(j)}\Big{\\|}_{\frac{8}{5}}\leq\beta\\|\nabla_{h}\tilde{\Phi}^{n,(j)}\\|_{\infty}^{2}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\leq\beta({\mathcal{C}}^{})^{2}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\beta({\mathcal{C}}^{})^{2}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.75)
	$\displaystyle\Big{\\|}\beta\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\tilde{\Phi}^{n,(j)}+\tilde{\boldsymbol{m}}_{j})\cdot{\mathcal{A}}_{h}\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\Big{)}\tilde{\boldsymbol{m}}_{j}\Big{\\|}_{\frac{8}{5}}$
$\displaystyle\leq$	$\displaystyle\beta(\\|\nabla_{h}\tilde{\Phi}^{n,(j)}\\|_{8}+\\|\nabla_{h}\tilde{\boldsymbol{m}}_{j}\\|_{8})\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2}\cdot\\|\tilde{\boldsymbol{m}}_{j}\\|_{\infty}$
$\displaystyle\leq$	$\displaystyle C\beta({\mathcal{C}}^{*}+\tilde{\mathcal{C}})\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.76)
	$\displaystyle\Big{\\|}\alpha\tilde{\boldsymbol{m}}_{j}\times(\tilde{\boldsymbol{e}}^{n,(j)}\times\boldsymbol{f})\Big{\\|}_{\frac{8}{5}}\leq\alpha\\|\tilde{\boldsymbol{m}}_{j}\\|_{\infty}\cdot\\|\boldsymbol{f}\\|_{\infty}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\leq\frac{5\alpha}{4}C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\alpha C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.77)
	$\displaystyle\Big{\\|}\alpha\tilde{\boldsymbol{e}}^{n,(j)}\times(\tilde{\Phi}^{n,(j)}\times\boldsymbol{f})\Big{\\|}_{\frac{8}{5}}\leq\alpha\\|\tilde{\Phi}^{n,(j)}\\|_{\infty}\cdot\\|\boldsymbol{f}\\|_{\infty}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\leq\frac{9\alpha}{8}C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\alpha C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.78)

in which the regularity estimate (5.65) and the functional bound (5.72) have been repeatedly applied, along with the fact that $\|g\|_{8}\leq C\|g\|_{\infty}$ , $\|g\|_{\frac{8}{5}}\leq C\|g\|_{2}$ (for any grid function $g$ ). Also notice an $\|\cdot\|_{\infty}$ bound for the external force term: $\|\boldsymbol{f}\|_{\infty}\leq C_{0}$ . As a result, a substitution of (5.75)-(5.78) into (5.74) yields the nonlinear error estimate (5.73), by taking $\tilde{M}=C(\beta(({\mathcal{C}}^{*})^{2}+{\mathcal{C}}^{*}+\tilde{\mathcal{C}})+\alpha C_{0})$ . This completes the proof of Lemma 5.

Before proceeding into the formal error estimate, we make the following a-priori assumption for the numerical error function at the previous time step:

\|\boldsymbol{e}^{n}\|_{2}\leq k^{\frac{15}{8}}+h^{\frac{15}{8}}.

(5.79)

Such an assumption will be recovered by the convergence analysis at the next time step $t^{n+1}$ .

Error estimate at Stage 1 Taking a discrete inner product with (5.68) by $2\tilde{\boldsymbol{e}}^{n,(2)}$ gives

\|\tilde{\boldsymbol{e}}^{n,(2)}\|_{2}^{2}-\|\boldsymbol{e}^{n}\|_{2}^{2}+\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\|_{2}^{2}+\frac{\beta}{2}k\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\|_{2}^{2}=0,

(5.80)

with an application of the summation-by-parts formula (3.22), due to the discrete homogeneous Neumann boundary condition for $\tilde{\boldsymbol{e}}^{n,(2)}$ . As a result, the following estimates available:

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq$	$\displaystyle\\|\boldsymbol{e}^{n}\\|_{2}\leq k^{\frac{15}{8}}+h^{\frac{15}{8}},$	(5.81)
$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq$	$\displaystyle\sqrt{2}\beta^{-\frac{1}{2}}k^{-\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq\sqrt{2}\beta^{-\frac{1}{2}}\Big{(}k^{\frac{11}{8}}+\frac{h^{\frac{15}{8}}}{k^{\frac{1}{2}}}\Big{)}$
$\displaystyle\leq$	$\displaystyle C(k^{\frac{11}{8}}+h^{\frac{11}{8}})\leq k^{\frac{5}{4}}+h^{\frac{5}{4}},$

in which the linear refinement requirement, $C_{1}h\leq k\leq C_{2}h$ , has been applied. Subsequently, the $\|\cdot\|_{\infty}$ and $\|\cdot\|_{W_{h}^{1,8}}$ bound for the numerical error function $\tilde{\boldsymbol{e}}^{n,(2)}$ could be derived, with the help of inverse inequalities (3.23), (3.24) (by taking $q=8$ ) in Lemma 2:

	$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{\infty}\leq\gamma{h}^{-1/2}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2})\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{1}{2}}}+h^{\frac{3}{4}}\Big{)}\leq\frac{1}{8},$		(5.82)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{8}\leq\gamma{h}^{-\frac{9}{8}}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{9}{8}}}+h^{\frac{1}{8}}\Big{)}\leq 1,$		(5.83)

provided that $k$ and $h$ are sufficiently small, the linear refinement requirement, $C_{1}h\leq k\leq C_{2}h$ . In turn, we obtain the following functional bound for the numerical solution $\tilde{\boldsymbol{m}}_{2}$ at the first Runge-Kutta stage, which will be useful in the later analysis:

	$\displaystyle\\|\tilde{\boldsymbol{m}}_{2}\\|_{\infty}\leq\\|\tilde{\Phi}^{n,(2)}\\|_{\infty}+\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{\infty}\leq\frac{9}{8}+\frac{1}{8}=\frac{5}{4},$		(5.84)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{m}}_{2}\\|_{8}\leq\\|\nabla_{h}\tilde{\Phi}^{n,(2)}\\|_{8}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{8}\leq{\mathcal{C}}^{*}+1=\tilde{\mathcal{C}},$		(5.85)

with an application of triangular inequality, along with the regularity estimate (5.65).

Error estimate at Stage 2 Taking a discrete inner product with (5.69) by $2\tilde{\boldsymbol{e}}^{n,(3)}$ gives

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{2}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}$		(5.86)
		$\displaystyle\quad=\frac{\beta}{2}k\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\rangle+k\langle{\mathcal{N}LE}^{n,(2)},\tilde{\boldsymbol{e}}^{n,(3)}\rangle,$		(5.86)

with an application of the summation-by-parts formula (3.22). The first term on the right hand side could be controlled by the Cauchy inequality:

\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\rangle\leq\frac{1}{2}(\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\|_{2}^{2}+\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\|_{2}^{2}).

(5.87)

Regarding the inner product term associated with the nonlinear error, we observe that an application of Lemma 5 implies the following estimates:

		$\displaystyle\\|{\mathcal{N}LE}^{n,(2)}\\|_{\frac{8}{5}}\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}),\quad\mbox{so that}$		(5.88)
		$\displaystyle\langle{\mathcal{N}LE}^{n,(2)},\tilde{\boldsymbol{e}}^{n,(3)}\rangle\leq\\|{\mathcal{N}LE}^{n,(2)}\\|_{\frac{8}{5}}\cdot\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2})\cdot\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{2}\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{144}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+(\frac{\tilde{M}}{2}+\frac{36\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}^{2},$

based on the regularity estimate (5.65) and the functional bound (5.84)-(5.85). Moreover, an application of the discrete Sobolev inequality (3.25) (in Lemma 3) indicates that

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}\leq C\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{4}\leq$	$\displaystyle\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{1}{4}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{3}{4}}),\quad\mbox{so that}$
$\displaystyle(\frac{\tilde{M}}{2}+\frac{36\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}^{2}\leq$	$\displaystyle\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{1}{2}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{3}{2}})$
$\displaystyle\leq$	$\displaystyle\mathcal{C}\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{24}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2},$	(5.89)

in which the Young’s inequality has been applied in the last step. Therefore, a substitution of (5.87)-(5.89) into (5.86) yields

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{5\beta}{24}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}$		(5.90)
		$\displaystyle-\frac{37\beta}{144}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}\leq\frac{\tilde{M}k}{2}\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\mathcal{C}k\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}.$		(5.90)

Furthermore, its combination with (5.80) gives

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}-\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}$		(5.91)
		$\displaystyle+\frac{35\beta}{144}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{5\beta}{24}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}\leq\frac{\tilde{M}k}{2}\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\mathcal{C}k\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}.$		(5.91)

Consequently, with an application of the a-priori estimate (5.81), we obtain

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq$	$\displaystyle\Big{(}\frac{1+\frac{\tilde{M}k}{2}}{1-\mathcal{C}k}\Big{)}^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq 2(k^{\frac{15}{8}}+h^{\frac{15}{8}}),$	(5.92)
$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq$	$\displaystyle\sqrt{4.8}\beta^{-\frac{1}{2}}k^{-\frac{1}{2}}(1+\frac{\tilde{M}k}{2})^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq\sqrt{5}\beta^{-\frac{1}{2}}\Big{(}k^{\frac{11}{8}}+\frac{h^{\frac{15}{8}}}{k^{\frac{1}{2}}}\Big{)}$
$\displaystyle\leq$	$\displaystyle C(k^{\frac{11}{8}}+h^{\frac{11}{8}})\leq k^{\frac{5}{4}}+h^{\frac{5}{4}},$

under the linear refinement requirement, $C_{1}h\leq k\leq C_{2}h$ . Similarly, the $\|\cdot\|_{\infty}$ and $\|\cdot\|_{W_{h}^{1,8}}$ bound for both the numerical error function $\tilde{\boldsymbol{e}}^{n,(3)}$ and the numerical solution $\tilde{\boldsymbol{m}}_{3}$ could be derived as follows

	$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\infty}\leq\gamma{h}^{-1/2}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2})\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{1}{2}}}+h^{\frac{3}{4}}\Big{)}\leq\frac{1}{8},$		(5.93)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{8}\leq\gamma{h}^{-\frac{9}{8}}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{9}{8}}}+h^{\frac{1}{8}}\Big{)}\leq 1,$		(5.94)
	$\displaystyle\\|\tilde{\boldsymbol{m}}_{3}\\|_{\infty}\leq\\|\tilde{\Phi}^{n,(3)}\\|_{\infty}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\infty}\leq\frac{9}{8}+\frac{1}{8}=\frac{5}{4},$		(5.95)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{m}}_{3}\\|_{8}\leq\\|\nabla_{h}\tilde{\Phi}^{n,(3)}\\|_{8}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{8}\leq{\mathcal{C}}^{*}+1=\tilde{\mathcal{C}}.$		(5.96)

Error estimate at Stage 3 Taking a discrete inner product with (5.70) by $2\tilde{\boldsymbol{e}}^{n,(4)}$ gives

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{2\beta}{3}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$		(5.97)
		$\displaystyle=-\frac{2\beta}{3}k\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\rangle-\frac{\beta}{6}k\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\rangle+k\langle{\mathcal{N}LE}^{n,(3)},\tilde{\boldsymbol{e}}^{n,(4)}\rangle,$		(5.97)

with an application of the summation-by-parts formula (3.22). The first two terms on the right hand side could be analyzed in the same way as in (5.40)-(5.41)

	$\displaystyle-\frac{2}{3}\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\rangle\leq\frac{2}{9}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{1}{2}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2},$		(5.98)
	$\displaystyle-\frac{1}{6}\langle\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)},\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\rangle\leq\frac{1}{12}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{1}{12}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2},$		(5.99)

Again, the nonlinear error term, as well as the corresponding inner product, could be analyzed in a similar fashion:

		$\displaystyle\\|{\mathcal{N}LE}^{n,(3)}\\|_{\frac{8}{5}}\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}),$		(5.100)
		$\displaystyle\langle{\mathcal{N}LE}^{n,(3)},\tilde{\boldsymbol{e}}^{n,(4)}\rangle\leq\\|{\mathcal{N}LE}^{n,(3)}\\|_{\frac{8}{5}}\cdot\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2})\cdot\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{2}\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{24}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+(\frac{\tilde{M}}{2}+\frac{6\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2},$
		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}\leq C\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{4}\leq\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{1}{4}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{3}{4}}),$
		$\displaystyle(\frac{\tilde{M}}{2}+\frac{6\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2}\leq\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{1}{2}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{3}{2}})$
		$\displaystyle\qquad\qquad\qquad\leq\mathcal{C}\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\frac{\beta}{48}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2},$

based on the a-priori bound estimate (5.95)-(5.96) in the second RK stage and the regularity estimate (5.65). Subsequently, a substitution of (5.98)-(5.100) into (5.97) gives

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{16}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$		(5.101)
		$\displaystyle-\frac{2\beta}{9}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}-\frac{\beta}{8}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}\leq\frac{\tilde{M}k}{2}\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\mathcal{C}k\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2},$		(5.101)

and its combination with (5.91) yields

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}-\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}$		(5.102)
		$\displaystyle+\frac{\beta}{48}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{12}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{16}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$
		$\displaystyle\leq\frac{\tilde{M}k}{2}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2})+\mathcal{C}k(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}).$

Similarly, with the help of the a-priori estimates (5.81), (5.92), in the first and second RK stages, respectively, the following rough error estimates could be derived:

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}\leq$	$\displaystyle\Big{(}\frac{1+{\mathcal{C}}k}{1-\mathcal{C}k}\Big{)}^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq 2(k^{\frac{15}{8}}+h^{\frac{15}{8}}),$	(5.103)
$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}\leq$	$\displaystyle 4\beta^{-\frac{1}{2}}k^{-\frac{1}{2}}(1+{\mathcal{C}}k)^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq 5\beta^{-\frac{1}{2}}\Big{(}k^{\frac{11}{8}}+\frac{h^{\frac{15}{8}}}{k^{\frac{1}{2}}}\Big{)}$
$\displaystyle\leq$	$\displaystyle C(k^{\frac{11}{8}}+h^{\frac{11}{8}})\leq k^{\frac{5}{4}}+h^{\frac{5}{4}},$

and the $\|\cdot\|_{\infty}$ and $\|\cdot\|_{W_{h}^{1,8}}$ bound for both the numerical error function $\tilde{\boldsymbol{e}}^{n,(4)}$ and the numerical solution $\tilde{\boldsymbol{m}}_{4}$ also becomes available:

	$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\infty}\leq\gamma{h}^{-1/2}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2})\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{1}{2}}}+h^{\frac{3}{4}}\Big{)}\leq\frac{1}{8},$		(5.104)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{8}\leq\gamma{h}^{-\frac{9}{8}}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}\leq\gamma\Big{(}\frac{k^{\frac{5}{4}}}{h^{\frac{9}{8}}}+h^{\frac{1}{8}}\Big{)}\leq 1,$		(5.105)
	$\displaystyle\\|\tilde{\boldsymbol{m}}_{4}\\|_{\infty}\leq\\|\tilde{\Phi}^{n,(4)}\\|_{\infty}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\infty}\leq\frac{9}{8}+\frac{1}{8}=\frac{5}{4},$		(5.106)
	$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{m}}_{4}\\|_{8}\leq\\|\nabla_{h}\tilde{\Phi}^{n,(4)}\\|_{8}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{8}\leq{\mathcal{C}}^{*}+1=\tilde{\mathcal{C}}.$		(5.107)

Error estimate at Stage 4 Taking a discrete inner product with (5.71) by $2\\ e^{n+1}$ gives

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}-k\langle\tau_{0}^{n},\boldsymbol{e}^{n+1}\rangle$		(5.108)
		$\displaystyle=-\frac{1}{3}k\langle{\mathcal{N}LE}^{n,(2)},\boldsymbol{e}^{n+1}\rangle-\frac{1}{3}k\langle{\mathcal{N}LE}^{n,(3)},\boldsymbol{e}^{n+1}\rangle+\frac{2}{3}k\langle{\mathcal{N}LE}^{n,(4)},\boldsymbol{e}^{n+1}\rangle.$		(5.108)

The local truncation error inner product term could be controlled in a straightforward way:

\langle\tau_{0}^{n},\boldsymbol{e}^{n+1}\rangle\leq\frac{1}{2}(\|\tau_{0}^{n}\|_{2}^{2}+\|\boldsymbol{e}^{n+1}\|_{2}^{2}).

(5.109)

The nonlinear error terms could be analyzed in a similar manner:

		$\displaystyle\\|{\mathcal{N}LE}^{n,(4)}\\|_{\frac{8}{5}}\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}),$		(5.110)
		$\displaystyle-\frac{1}{3}\langle{\mathcal{N}LE}^{n,(2)},\boldsymbol{e}^{n+1}\rangle\leq\frac{1}{3}\\|{\mathcal{N}LE}^{n,(2)}\\|_{\frac{8}{5}}\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{1}{3}\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2})\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{6}\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{144}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+(\frac{\tilde{M}}{6}+\frac{4\tilde{M}^{2}}{\beta})\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2},$
		$\displaystyle-\frac{1}{3}\langle{\mathcal{N}LE}^{n,(3)},\boldsymbol{e}^{n+1}\rangle\leq\frac{1}{3}\\|{\mathcal{N}LE}^{n,(3)}\\|_{\frac{8}{5}}\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{1}{3}\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2})\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{6}\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{36}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+(\frac{\tilde{M}}{6}+\frac{\tilde{M}^{2}}{\beta})\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2},$
		$\displaystyle\frac{2}{3}\langle{\mathcal{N}LE}^{n,(4)},\boldsymbol{e}^{n+1}\rangle\leq\frac{2}{3}\\|{\mathcal{N}LE}^{n,(4)}\\|_{\frac{8}{5}}\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{2}{3}\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2})\cdot\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{3}\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\frac{\beta}{48}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+(\frac{\tilde{M}}{3}+\frac{16\tilde{M}^{2}}{3\beta})\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2}.$

In turn, a substitution of (5.109), (5.110) into (5.108) leads to

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}-\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}-\frac{\beta}{48}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$		(5.111)
		$\displaystyle-\frac{\beta}{144}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}-\frac{\beta}{36}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}-\frac{1}{2}k(\\|\tau_{0}^{n}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}\\|_{2}^{2})$
		$\displaystyle\leq\frac{\tilde{M}k}{6}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+2\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2})+(\frac{\tilde{2M}}{3}+\frac{31\tilde{M}^{2}}{3\beta})\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2},$

and its combination with (5.102) yields

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}-\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}$		(5.112)
		$\displaystyle+\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\frac{\beta}{72}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{18}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{24}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$
		$\displaystyle\leq\mathcal{C}k(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2})+\hat{M}k\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2}+\frac{1}{2}k(\\|\tau_{0}^{n}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}).$

with $\hat{M}=\frac{\tilde{2M}}{3}+\frac{31\tilde{M}^{2}}{3\beta}$ . To control the terms $\hat{M}k\|\boldsymbol{e}^{n+1}\|_{\frac{8}{3}}^{2}$ , we see that an application of triangle inequality implies that

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}\leq\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}+\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}},\quad\mbox{so that}$		(5.113)
		$\displaystyle\hat{M}\\|\boldsymbol{e}^{n+1}\\|_{\frac{8}{3}}^{2}\leq 2\hat{M}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2}+\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2}).$		(5.113)

Meanwhile, an application of inverse inequality (3.24) (by taking $q=\frac{8}{3}$ , in Lemma 2) results in

		$\displaystyle\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}\leq\gamma{h}^{-\frac{3}{8}}\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2},\quad\mbox{so that}$		(5.114)
		$\displaystyle 2\hat{M}k\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2}\leq 2\hat{M}\gamma^{2}k{h}^{-\frac{3}{4}}\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}\leq\frac{1}{2}\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2},$		(5.114)

provided that $2\hat{M}\gamma^{2}k{h}^{-\frac{3}{4}}\leq\frac{1}{2}$ , which is always valid under the linear refinement requirement, $C_{1}h\leq k\leq C_{2}h$ , and the assumption that $k$ and $h$ are sufficiently small. In addition, we recall the preliminary estimate (5.100) for $\|\tilde{\boldsymbol{e}}^{n,(4)}\|_{\frac{8}{3}}$ :

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}\leq C\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{4}\leq\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{1}{4}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{3}{4}}),\quad\mbox{so that}$		(5.115)
		$\displaystyle 2\hat{M}\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{\frac{8}{3}}^{2}\leq\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{1}{2}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{\frac{3}{2}})$
		$\displaystyle\qquad\qquad\qquad\leq\mathcal{C}\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\frac{\beta}{48}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}.$

Therefore, a substitution of (5.113)-(5.115) into (5.112) gives

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}-\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}$		(5.116)
		$\displaystyle+\frac{1}{2}\\|\boldsymbol{e}^{n+1}-\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}+\frac{\beta}{72}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{18}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{48}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$
		$\displaystyle\leq\mathcal{C}k(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2})+\frac{1}{2}k(\\|\tau_{0}^{n}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}).$

Moreover, by making use of the triangle inequalities

		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq\\|\boldsymbol{e}^{n}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2},$		(5.117)
		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq\\|\boldsymbol{e}^{n}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2},$
		$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}\leq\\|\boldsymbol{e}^{n}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(2)}-\boldsymbol{e}^{n}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}-\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(4)}-\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2},$

we arrive at

		$\displaystyle\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}-\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\frac{\beta}{72}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{18}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{48}k\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(4)}\\|_{2}^{2}$		(5.118)
		$\displaystyle\leq\mathcal{C}k\\|\boldsymbol{e}^{n}\\|_{2}^{2}+\frac{1}{2}k(\\|\tau_{0}^{n}\\|_{2}^{2}+\\|\boldsymbol{e}^{n+1}\\|_{2}^{2}),$		(5.118)

provided that $k$ is sufficiently small. In turn, an application of discrete Gronwall inequality [52] yields the desired convergence estimate at the next time step

\|\boldsymbol{e}^{n+1}\|_{2}\leq\mathcal{C}(k^{2}+h^{2}),

(5.119)

based on the fact that $\|\tau_{0}^{n}\|^{2}\leq{\mathcal{C}}(k^{2}+h^{2})$ . As a result, we see that the a-priori assumption (5.79) has also been validated at the next time step $t^{n+2}$ , provided that $k$ and $h$ are sufficiently small. By mathematical induction, this completes the proof of Theorem 4.

Remark 5.4.

For the multi-step IMEX numerical schemes to various nonlinear PDEs, there have been some existing works of convergence analysis [9, 18], etc. Meanwhile, for the IMEX-RK numerical schemes, the theoretical works have been very limited, due to the theoretical challenges associated with the multi-stage nature, lack of sufficient numerical diffusion, etc. A linearized stability analysis is provided for the IMEX-RK1 and IMEX-RK2 schemes for the diffusion problem [30], as well as the error estimate for the constant-coefficient diffusion equation. The convergence analysis has been reported for the IMEX-RK numerical methods to various nonlinear PDEs, such as the convection-diffusion equation [29], the porous media equation [28, 31], etc. On the other hand, the degree of nonlinearity of the LL equation (even the simplified version (5.45), with only the damping term) is higher than these reported models. As a result, the associated theoretical analysis presented in this article contains more techniques than the existing works.

Remark 5.5.

At the intermediate Runge-Kutta stages, the $\ell^{2}$ error bound estimate (5.81), (5.92), (5.103), the associated $\ell^{\infty}$ and $W_{h}^{1,8}$ error bounds (5.82)-(5.83), (5.93)-(5.94), (5.104)-(5.105), stand for a rough error estimate. These error bound estimates do not preserve a full accuracy order; instead, the motivation of these analyses is to obtain a rough bound of the numerical error function, so that a preliminary bound becomes available for the numerical solution at the associated Runge-Kutta stages, as derived in (5.84)-(5.85), (5.95)-(5.96), (5.106)-(5.107). As a result of these preliminary bounds, the nonlinear error terms could be analyzed with the help of Lemma 5, so that the nonlinear analysis would go through.

In comparison, with these nonlinear analyses established, the numerical error inequalities (5.80), (5.91), (5.102) stand for a refined error estimate. Finally, the derived estimate (5.118) becomes an inequality in which we are able to derive the desired accuracy order for the numerical error.

A combination of rough error estimate and refined error estimate has been reported for various nonlinear PDEs, such as ternary Flory-Huggins-Cahn-Hilliard system [12], Poisson-Nernst-Planck system [22], the porous medium equation by an energetic variational approach [14, 13], the reaction-diffusion system with detailed balance [21], etc. This article reports the application of such a technique to the LL equation for the first time.

Remark 5.6.

For simplicity, we only consider the damping term in the convergence analysis, as the PDE system formulated as (5.45). For the full LL equation (2.4), the convergence estimate may still go through, under a large damping parameter assumption, combined with a technical requirement of $k={\mathcal{O}}(h^{2})$ ; the details are left to interested readers, for the sake of brevity. Of course, such a requirement only stands for a theoretical difficulty, and this restrictive time step constraint is not needed in the practical computation, as demonstrated in extensive numerical experiments.

Remark 5.7.

A theoretical analysis of the IMEX-RK3 numerical algorithm, including the linearized stability estimate and optimal rate convergence analysis, is expected to be even more challenging, due to the complicated diffusion coefficient stencil in the RK stages. This theoretical work will be studied in a forthcoming paper.

6 Conclusions

In this paper, we propose implicit-explicit Runge-Kutta numerical methods for solving the Landau-Lifshitz equation. By introducing an artificial damping term, IMEX-RK methods only solve a few linear systems with constant coefficients and SPD structures, regardless of the damping parameter in a magnetic material. Accuracy, efficiency, and the insensitive dependence on the artificial damping parameter of IMEX-RK2 and IMEX-RK3 have been verified in both the 1-D and 3-D computations. Micromagnetics simulations using the full Landau-Lifshitz equation are conducted, including three stable structures and the first benchmark problem from NIST. Reasonable results are generated by the IMEX-RK methods, in qualitative and quantitative agreements with available results. In addition, the linearized stability estimate and optimal rate convergence analysis are presented for an alternate IMEX-RK2 numerical scheme, the SSP-IMEX-RK2 algorithm. This analysis has provided a theoretical evidence of the robust performance of the IMEX-RK schemes.

Acknowledgments

This work is supported in part by the grants NSFC 11971021 (J. Chen) and NSF DMS-2012269 (C. Wang).

References

[1] L. Landau and E. Lifshitz, “On the theory of the dispersion of magnetic permeability in ferromagnetic bodies,” Physikalische Zeitschrift Der Sowjetunion, vol. 63, no. 9, pp. 153–169, 1935.
[2] T. Gilbert, “A Lagrangian formulation of the gyromagnetic equation of the magnetization field,” Physical Review, vol. 100, p. 1243, 1955.
[3] R. An, “Optimal error estimates of linearized Crank–Nicolson Galerkin method for Landau–Lifshitz equation,” Journal of Scientific Computing, vol. 69, no. 1, pp. 1–27, 2016.
[4] R. An and W. Sun, “Analysis of backward Euler projection FEM for the Landau-Lifshitz equation,” IMA Journal of Numerical Analysis, vol. 42, no. 3, pp. 2336–2360, 2022.
[5] S. Bartels and P. Andreas, “Convergence of an implicit finite element method for the Landau–Lifshitz–Gilbert equation,” SIAM journal on numerical analysis, vol. 44, no. 4, pp. 1405–1419, 2006.
[6] S. Bartels, J. Ko, and A. Prohl, “Numerical analysis of an explicit approximation scheme for the Landau–Lifshitz–Gilbert equation,” Mathematics of Computation, vol. 77, no. 262, pp. 773–788, 2008.
[7] W. Chen, Y. Liu, C. Wang, and S. Wise, “An optimal-rate convergence analysis of a fully discrete finite difference scheme for Cahn-Hilliard-Hele-Shaw equation,” Math. Comp., vol. 85, pp. 2231–2257, 2016.
[8] W. Chen, C. Wang, S. Wang, X. Wang, and S. Wise, “Energy stable numerical schemes for ternary Cahn-Hilliard system,” J. Sci. Comput., vol. 84, p. 27, 2020.
[9] K. Cheng and C. Wang, “Long time stability of high order multi-step numerical schemes for two-dimensional incompressible Navier-Stokes equations,” SIAM Journal on Numerical Analysis, vol. 54, pp. 3123–3144, 2016.
[10] P. Ciarlet, The Finite Element Method for Elliptic Problems. Elsevier Science, 1978.
[11] S. Conde, S. Gottlieb, Z. Grant, and J. Shadid, “Implicit and implicit-explicit strong stability preserving Runge-Kutta methods with high linear order,” Journal of Scientific Computing, vol. 83, no. 2, pp. 667–690, 2017.
[12] L. Dong, C. Wang, S. Wise, and Z. Zhang, “Optimal rate convergence analysis of a numerical scheme for the ternary Cahn-Hilliard system with a Flory-Huggins-deGennes energy potential,” Journal of Computational and Applied Mathematics, vol. 406, p. 114474, 2022.
[13] C. Duan, C. Liu, C. Wang, and X. Yue, “Convergence analysis of a numerical scheme for the porous medium equation by an energetic variational approach,” Numerical Mathematics: Theory, Methods and Applications, vol. 13, pp. 1–18, 2020.
[14] C. Duan, W. Chen, C. Liu, C. Wang, and X. Yue, “A second order accurate, energy stable numerical scheme for one-dimensional porous medium equation by an energetic variational approach,” Communications in Mathematical Sciences, vol. 20, no. 4, pp. 987–1024, 2022.
[15] W. E and X. Wang, “Numerical methods for the Landau–Lifshitz equation,” SIAM journal on numerical analysis, pp. 1647–1665, 2001.
[16] A. Fuwa, T. Ishiwata, and M. Tsutsumi, “Finite difference scheme for the Landau–Lifshitz equation,” Japan journal of industrial and applied mathematics, vol. 29, no. 1, pp. 83–110, 2012.
[17] S. Gottlieb, C. Shu, and E. Tadmor, “Strong stability-preserving high-order time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001.
[18] S. Gottlieb and C. Wang, “Stability and convergence analysis of fully discrete Fourier collocation spectral method for 3-D viscous Burgers’ equation,” Journal of Scientific Computing, vol. 53, pp. 102–128, 2012.
[19] Z. Guan, C. Wang, and S. Wise, “A convergent convex splitting scheme for the periodic nonlocal Cahn-Hilliard equation,” Numer. Math., vol. 128, pp. 377–406, 2014.
[20] Z. Guan, J. Lowengrub, and C. Wang, “Convergence analysis for second order accurate schemes for the periodic nonlocal Allen-Cahn and Cahn-Hilliard equations,” Math. Methods Appl. Sci., vol. 40, no. 18, pp. 6836–6863, 2017.
[21] C. Liu, C. Wang, Y. Wang, and S. Wise, “Convergence analysis of the variational operator splitting scheme for a reaction-diffusion system with detailed balance,” SIAM Journal on Numerical Analysis, vol. 60, no. 2, pp. 781–803, 2022.
[22] C. Liu, C. Wang, S. Wise, X. Yue, and S. Zhou, “A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system,” Mathematics of Computation, vol. 90, pp. 2071–2106, 2021.
[23] C. Wang and J.-G. Liu, “Convergence of gauge method for incompressible flow,” Mathematics of Computation, vol. 69, pp. 1385–1407, 2000.
[24] R. Samelson, R. Temam, C. Wang, and S. Wang, “Surface pressure Poisson equation formulation of the primitive equations: Numerical schemes,” SIAM Journal on Numerical Analysis, vol. 41, pp. 1163–1194, 2003.
[25] C. Wang, J.-G. Liu, and H. Johnston, “Analysis of a fourth order finite difference method for incompressible Boussinesq equation,” Numerische Mathematik, vol. 97, pp. 555–594, 2004.
[26] X.-P. Wang, C. Garcıa-Cervera, and W. E, “A Gauss–Seidel projection method for micromagnetics simulations,” Journal of Computational Physics, vol. 171, no. 1, pp. 357–372, 2001.
[27] C. Xie, C. Garcıa-Cervera, C. Wang, Z. Zhou, and J. Chen, “Second-order semi-implicit projection methods for micromagnetics simulations,” Journal of Computational Physics, vol. 404, p. 109104, 2020.
[28] C. Nan and H. Song, “Error estimates of local discontinuous Galerkin method with implicit-explicit Runge Kutta for two-phase miscible flow in porous media,” Applied Numerical Mathematics, vol. 169, pp. 334–350, 2021.
[29] H. Wang, S. Wang, Q. Zhang, and C. Shu, “Local discontinuous Galerkin methods with explicit-implicit time marching for multi-dimensional convection-diffusion problems,” ESAIM: Mathematical Modeling and Numerical Analysis, vol. 50, pp. 1083–1105, 2016.
[30] H. Wang, Q. Zhang, S. Wang, and C. Shu, “Local discontinuous Galerkin methods with explicit–implicit–null time discretizations for solving nonlinear diffusion problems,” Science China Mathematics, vol. 63, no. 1, pp. 183–204, 2020.
[31] H. Wang, J. Zheng, F. You, H. Guo, and Q. Zhang, “Local discontinuous Galerkin method with explicit-implicit time marching for incompressible miscible displacement problem in porous media,” Journal of Scientific Computing, vol. 78, pp. 1–28, 2019.
[32] U. Ascher, S. Ruuth, and R. Spiteri, “Implicit-explicit Runge–Kutta methods for time–dependent partial differential equations,” Applied Numerical Mathematics, vol. 25, no. 2-3, pp. 151–167, 1997.
[33] Y. Zheng and J. Zhu, “Switching field variation in patterned submicron magnetic film elements,” Journal of Applied Physics, vol. 81, no. 8, pp. 5471–5473, 1997.
[34] T. Schrefl, D. Suess, W. Scholz, H. Forster, V. Tsiantos, and J. Fidler, “Finite element micromagnetics,” Computational Electromagnetics, pp. 165–181, 2003.
[35] F. Alouges and P. Jaisson, “Convergence of a finite element discretization for the Landau–Lifshitz equations in micromagnetism,” Mathematical Models and Methods in Applied Sciences, vol. 16, no. 02, pp. 299–316, 2006.
[36] C. García-Cervera and W. E, “Improved Gauss–Seidel projection method for micromagnetics simulations,” IEEE transactions on magnetics, vol. 39, no. 3, pp. 1766–1770, 2003.
[37] I. Cimrák, “Error estimates for a semi–implicit numerical scheme solving the Landau–Lifshitz equation with an exchange field,” IMA journal of numerical analysis, vol. 25, no. 3, pp. 611–634, 2005.
[38] H. Gao, “Optimal error estimates of a linearized backward euler fem for the Landau–Lifshitz equation,” SIAM Journal on Numerical Analysis, vol. 52, no. 5, pp. 2574–2593, 2014.
[39] S. Boscarino, F. Filbet, and G. Russo, “High order semi–implicit schemes for time dependent partial differential equations,” Journal of Scientific Computing, vol. 68, no. 3, pp. 975–1001, 2016.
[40] J. Chen, C. Wang, and C. Xie, “Convergence analysis of a second–order semi–implicit projection method for Landau–Lifshitz equation,” Applied Numerical Mathematics, vol. 168, pp. 55–74, 2021.
[41] Y. Cai, J. Chen, C. Wang, and C. Xie, “A second-order numerical method for Landau–Lifshitz–Gilbert equation with large damping parameters,” Journal of Computational Physics, vol. 451, p. 110831, 2022.
[42] Micromagnetic Modeling Activity Group and Technology, “2000,” https://www.ctcms.nist.gov/~rdm/mumag.org.html.
[43] R. An, H. Gao, and W. Sun, “Optimal error analysis of Euler and Crank–Nicolson projection finite difference schemes for Landau–Lifshitz equation,” SIAM J. Numer. Anal., vol. 59, no. 3, pp. 1639–1662, 2021.
[44] M. Kruzík and A. Prohl, “Recent developments in the modeling, analysis, and numerics of ferromagnetism,” SIAM Rev., vol. 48, no. 3, pp. 439–483, 2006.
[45] I. Cimrák, “A survey on the numerics and computations for the Landau-Lifshitz equation of micromagnetism,” Arch. Comput. Methods Eng., vol. 15, no. 3, pp. 277–309, 2008.
[46] C. J. García-Cervera, “Numerical micromagnetics: a review,” Bol. Soc. Esp. Mat. Apl., vol. 39, pp. 103–135, 2007.
[47] P. Li, C. Xie, R. Du, J. Chen, and X. Wang, “Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation,” J. Comput. Phys., vol. 401, p. 109046, 2020.
[48] G. Akrivis, M. Feischl, B. Kovács, and C. Lubich, “Higher-order linearly implicit full discretization of the Landau-Lifshitz-Gilbert equation,” Math. Comp., vol. 90, pp. 995–1038, 2021.
[49] W. Brown, Micromagnetics. Interscience Tracts on Physics and Astronomy. Interscience Publishers (John Wiley and Sons), New York-London, 1963.
[50] Y. Cai, J. Chen, C. Wang, and C. Xie, “Error analysis of a linear numerical scheme for the Landau-Lifshitz equation with large damping parameters,” Math. Meth. Appl. Sci., vol. 46, pp. 18 952–18 974, 2023.
[51] C. Gear, “Hybrid methods for initial value problems in ordinary differential equations,” Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, vol. 2, no. 1, pp. 69–86, 1965.
[52] V. Girault and P. Raviart, Finite Element Methods for Navier-Stokes Equations, Theorems and Algorithms. Springer-Verlag, 1986.

	$\displaystyle\Big{\\|}\beta\|{\mathcal{A}}_{h}\nabla_{h}\tilde{\Phi}^{n,(j)}\|^{2}\tilde{\boldsymbol{e}}^{n,(j)}\Big{\\|}_{\frac{8}{5}}\leq\beta\\|\nabla_{h}\tilde{\Phi}^{n,(j)}\\|_{\infty}^{2}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\leq\beta({\mathcal{C}}^{})^{2}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\beta({\mathcal{C}}^{})^{2}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.75)
	$\displaystyle\Big{\\|}\beta\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\tilde{\Phi}^{n,(j)}+\tilde{\boldsymbol{m}}_{j})\cdot{\mathcal{A}}_{h}\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\Big{)}\tilde{\boldsymbol{m}}_{j}\Big{\\|}_{\frac{8}{5}}$
$\displaystyle\leq$	$\displaystyle\beta(\\|\nabla_{h}\tilde{\Phi}^{n,(j)}\\|_{8}+\\|\nabla_{h}\tilde{\boldsymbol{m}}_{j}\\|_{8})\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2}\cdot\\|\tilde{\boldsymbol{m}}_{j}\\|_{\infty}$
$\displaystyle\leq$	$\displaystyle C\beta({\mathcal{C}}^{*}+\tilde{\mathcal{C}})\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.76)
	$\displaystyle\Big{\\|}\alpha\tilde{\boldsymbol{m}}_{j}\times(\tilde{\boldsymbol{e}}^{n,(j)}\times\boldsymbol{f})\Big{\\|}_{\frac{8}{5}}\leq\alpha\\|\tilde{\boldsymbol{m}}_{j}\\|_{\infty}\cdot\\|\boldsymbol{f}\\|_{\infty}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\leq\frac{5\alpha}{4}C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\alpha C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.77)
	$\displaystyle\Big{\\|}\alpha\tilde{\boldsymbol{e}}^{n,(j)}\times(\tilde{\Phi}^{n,(j)}\times\boldsymbol{f})\Big{\\|}_{\frac{8}{5}}\leq\alpha\\|\tilde{\Phi}^{n,(j)}\\|_{\infty}\cdot\\|\boldsymbol{f}\\|_{\infty}\cdot\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}$
	$\displaystyle\qquad\qquad\qquad\leq\frac{9\alpha}{8}C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{\frac{8}{5}}\leq C\alpha C_{0}\\|\tilde{\boldsymbol{e}}^{n,(j)}\\|_{2},$	(5.78)

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq$	$\displaystyle\\|\boldsymbol{e}^{n}\\|_{2}\leq k^{\frac{15}{8}}+h^{\frac{15}{8}},$	(5.81)
$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}\leq$	$\displaystyle\sqrt{2}\beta^{-\frac{1}{2}}k^{-\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq\sqrt{2}\beta^{-\frac{1}{2}}\Big{(}k^{\frac{11}{8}}+\frac{h^{\frac{15}{8}}}{k^{\frac{1}{2}}}\Big{)}$
$\displaystyle\leq$	$\displaystyle C(k^{\frac{11}{8}}+h^{\frac{11}{8}})\leq k^{\frac{5}{4}}+h^{\frac{5}{4}},$

		$\displaystyle\\|{\mathcal{N}LE}^{n,(2)}\\|_{\frac{8}{5}}\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}),\quad\mbox{so that}$		(5.88)
		$\displaystyle\langle{\mathcal{N}LE}^{n,(2)},\tilde{\boldsymbol{e}}^{n,(3)}\rangle\leq\\|{\mathcal{N}LE}^{n,(2)}\\|_{\frac{8}{5}}\cdot\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\tilde{M}(\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}+\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2})\cdot\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}$
		$\displaystyle\leq\frac{\tilde{M}}{2}\\|\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+\frac{\beta}{144}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(2)}\\|_{2}^{2}+(\frac{\tilde{M}}{2}+\frac{36\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}^{2},$

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}\leq C\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{4}\leq$	$\displaystyle\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{1}{4}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{3}{4}}),\quad\mbox{so that}$
$\displaystyle(\frac{\tilde{M}}{2}+\frac{36\tilde{M}^{2}}{\beta})\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{\frac{8}{3}}^{2}\leq$	$\displaystyle\mathcal{C}(\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{1}{2}}\cdot\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{\frac{3}{2}})$
$\displaystyle\leq$	$\displaystyle\mathcal{C}\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2}+\frac{\beta}{24}\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}^{2},$	(5.89)

$\displaystyle\\|\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq$	$\displaystyle\Big{(}\frac{1+\frac{\tilde{M}k}{2}}{1-\mathcal{C}k}\Big{)}^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq 2(k^{\frac{15}{8}}+h^{\frac{15}{8}}),$	(5.92)
$\displaystyle\\|\nabla_{h}\tilde{\boldsymbol{e}}^{n,(3)}\\|_{2}\leq$	$\displaystyle\sqrt{4.8}\beta^{-\frac{1}{2}}k^{-\frac{1}{2}}(1+\frac{\tilde{M}k}{2})^{\frac{1}{2}}\\|\boldsymbol{e}^{n}\\|_{2}\leq\sqrt{5}\beta^{-\frac{1}{2}}\Big{(}k^{\frac{11}{8}}+\frac{h^{\frac{15}{8}}}{k^{\frac{1}{2}}}\Big{)}$
$\displaystyle\leq$	$\displaystyle C(k^{\frac{11}{8}}+h^{\frac{11}{8}})\leq k^{\frac{5}{4}}+h^{\frac{5}{4}},$

IMEX-RK methods for Landau-Lifshitz equation with arbitrary damping††thanks: Received date, and accepted date (The correct dates will be entered by the editor).

Abstract

keywords:

1 Introduction

2 The model and the proposed methods

2.1 Landau-Lifshitz equation

2.2 Second-order backward differentiation formula method

2.3 Implicit-explicit Runge-Kutta methods

Remark 2.1.

Remark 2.2.

3 Numerical results

Definition 3.1 (ℓ2\ell^{2} inner product, ∥⋅∥2{{\left\|\cdot\right\|}_{2}} norm).

Definition 3.2 ( ∥⋅∥∞{{\left\|\cdot\right\|}_{\infty}} and ∥⋅∥p{{\left\|\cdot\right\|}_{p}} norms in the discrete sense).

Lemma 1 (Summation by parts).

Lemma 2 (Inverse inequality).

Lemma 3 (Discrete Sobolev inequality).

3.1 Accuracy and efficiency test of IMEX-RK2

3.2 Accuracy and efficiency test of IMEX-RK3

3.3 Accuracy test of SSP-IMEX-RK2

3.4 Dependence on the damping parameter

Remark 3.3.

4 Micromagnetics simulations

4.1 Equilibrium states

4.2 Benchmark peoblem from NIST

5 Numerical stability and convergence analysis for the proposed SSP-IMEX-RK2 scheme

5.1 Linear stability estimate for the SSP-IMEX-RK2 scheme (2.21)

Remark 5.1.

5.2 Convergence analysis of the SSP-IMEX-RK2 scheme (2.21) for a simplified nonlinear LL equation

Theorem 4.

Proof 5.2.

Lemma 5.

Proof 5.3.

Remark 5.4.

Remark 5.5.

Remark 5.6.

Remark 5.7.

6 Conclusions

Acknowledgments

References

IMEX-RK methods for Landau-Lifshitz equation with arbitrary damping^†^†thanks: Received date, and accepted date (The correct dates will be entered by the editor).

Definition 3.1 ( $\ell^{2}$ inner product, ${{\left\|\cdot\right\|}_{2}}$ norm).

Definition 3.2 ( ${{\left\|\cdot\right\|}_{\infty}}$ and ${{\left\|\cdot\right\|}_{p}}$ norms in the discrete sense).