Simple maximum-principle preserving time-stepping methods for time-fractional Allen-Cahn equation

TL;DR

This paper introduces fast, adaptive time-stepping methods for the time-fractional Allen-Cahn equation that preserve the maximum principle and are supported by rigorous error analysis and numerical validation.

Contribution

The paper develops and analyzes fast L1 time-stepping schemes with adaptive meshes for the fractional Allen-Cahn equation, ensuring maximum principle preservation and providing sharp error estimates.

Findings

The proposed schemes preserve the discrete maximum principle.

Error estimates reflect the solution's time regularity.

Numerical experiments confirm effectiveness and theoretical analysis.

Abstract

Two fast L1 time-stepping methods, including the backward Euler and stabilized semi-implicit schemes, are suggested for the time-fractional Allen-Cahn equation with Caputo's derivative. The time mesh is refined near the initial time to resolve the intrinsically initial singularity of solution, and unequal time-steps are always incorporated into our approaches so that an adaptive time-stepping strategy can be used in long-time simulations. It is shown that the proposed schemes using the fast L1 formula preserve the discrete maximum principle. Sharp error estimates reflecting the time regularity of solution are established by applying the discrete fractional Gr\"{o}nwall inequality and global consistency analysis. Numerical experiments are presented to show the effectiveness of our methods and to confirm our analysis.

Equations

$${\partial }_t^{\alpha }u={\epsilon }^2\Delta u-f(u),\mathbf{x}\in \Omega ,0<t⩽T,$$

$$u(\mathbf{x},0)=u_0(\mathbf{x}),\mathbf{x}\in \bar{\Omega },$$

$$({\partial }_t^{\alpha }v)(t):=({\mathcal{I}}_t^{1-\alpha }{v}^{\prime })(t)={\int }_0^t{\omega }_{1-\alpha }(t-s)v^{\prime }(s)\mathrm{d}s,0<\alpha <1,$$

$$({\mathcal{I}}_t^{\mu }v)(t):={\int }_0^t{\omega }_{\mu }(t-s)v(s)\mathrm{d}s,\text{where}{\omega }_{\mu }(t):=t^{\mu -1}/\Gamma (\mu ).$$

$$E(t)⩽E(0),$$

$$E(t):={\int }_{\Omega }\left[\frac{{\epsilon }^2}{2}|\nabla u{|}^2+F(u)\right]\mathrm{d}\mathbf{x},F(u)=\frac{1}{4}(1-u^2{)}^2.$$

$$|u(\mathbf{x},t)|⩽1\text{for t>0}\text{if}|u(\mathbf{x},0)|⩽1.$$

$${\Pi }_{1}v:={\Pi }_{1,k}v\text{so that}({\Pi }_{1}v{)}^{\prime }(t)={\partial }_{\tau }{v}^{k-\frac{1}{2}}\text{for tk−1<t⩽tk and k≥1}.$$

$$\log |u_t(\mathbf{x},t)|\approx (\alpha -1)\log t+C(\mathbf{x})\text{as t→0.}$$

$$({\partial }_{\tau }^{\alpha }v{)}^n:={\int }_{t_0}^{t_n}{\omega }_{1-\alpha }(t_n-s)({\Pi }_{1}v{)}^{\prime }(s)\mathrm{d}s=\sum _{k=1}^n{a}_{n-k}^{(n)}{▽}_{\tau }{v}^k,$$

$$a_{n-k}^{(n)}:=\frac{1}{{\tau }_k}{\int }_{t_{k-1}}^{t_k}{\omega }_{1-\alpha }(t_n-s)\mathrm{d}s\text{for 1⩽k⩽n.}$$

$$a_{n-k}^{(n)}>0\text{and}{a}_{n-k-1}^{(n)}>a_{n-k}^{(n)}\text{for 1⩽k⩽n−1.}$$

$$\displaystyle\bigg{|}\omega_{1-\alpha}(t)-\sum_{\ell=1}^{N_{q}}\varpi^{\ell}e^{-\theta^{\ell}t}\bigg{|}\leqslant\epsilon,\quad\forall\,{t}\in[\Delta{t},\,T].$$

$$({\partial }_t^{\alpha }v)(t_n)$$

$$=\sum _{\ell =1}^{N_q}{\varpi }^{\ell }{e}^{-{\theta }^{\ell }{\tau }_n}{\mathcal{H}}^{\ell }(t_{n-1})+a_0^{(n)}{▽}_{\tau }{v}^n,n\ge 1,$$

$${\mathcal{H}}^{\ell }(t_k)$$

$$=e^{-{\theta }^{\ell }{\tau }_k}{\mathcal{H}}^{\ell }(t_{k-1})+b^{(k,l)}{▽}_{\tau }{v}^k,$$

$$b^{(k,l)}:=\frac{1}{{\tau }_k}{\int }_{t_{k-1}}^{t_k}{e}^{-{\theta }^{\ell }(t_k-s)}\mathrm{d}s,k\ge 1.$$

$$({\partial }_f^{\alpha }v{)}^n=a_0^{(n)}{▽}_{\tau }{v}^n+\sum _{\ell =1}^{N_q}{\varpi }^{\ell }{e}^{-{\theta }^{\ell }{\tau }_n}{\mathcal{H}}^{\ell }(t_{n-1}),n\ge 1,$$

$${\mathcal{H}}^{\ell }(t_k)=e^{-{\theta }^{\ell }{\tau }_k}{\mathcal{H}}^{\ell }(t_{k-1})+b^{(k,l)}{▽}_{\tau }{v}^k,k\ge 1,1⩽\ell ⩽N_q.$$

$${\mathcal{H}}^{\ell }(t_k)=\sum _{j=1}^k{e}^{-{\theta }^{\ell }(t_k-t_j)}{b}^{(j,l)}{▽}_{\tau }{v}^j,k\ge 1,1⩽l⩽N_q.$$

$$({\partial }_f^{\alpha }v{)}^n=a_0^{(n)}{▽}_{\tau }{v}^n+\sum _{k=1}^{n-1}\frac{{▽}_{\tau }{v}^k}{{\tau }_k}{\int }_{t_{k-1}}^{t_k}\sum _{\ell =1}^{N_q}{\varpi }^{\ell }{e}^{-{\theta }^{\ell }(t_n-s)}\mathrm{d}s=\sum _{k=1}^n{A}_{n-k}^{(n)}{▽}_{\tau }{v}^k,n\ge 1,$$

$$A_0^{(n)}:=a_0^{(n)},A_{n-k}^{(n)}:=\frac{1}{{\tau }_k}{\int }_{t_{k-1}}^{t_k}\sum _{\ell =1}^{N_q}{\varpi }^{\ell }{e}^{-{\theta }^{\ell }(t_n-s)}\mathrm{d}s,1⩽k⩽n-1,n\ge 1.$$

$$\mathbb{V}_{h}:=\big{\{}v\,|\,v=(v_{j})^{T}\;\;\text{for}\;\;0\leqslant{j}\leqslant{M_{2}-1},\;\;\text{with}\;\;v_{j}=(v_{i,j})^{T}\;\;\text{for}\;\;0\leqslant{i}\leqslant{M_{1}-1}\big{\}},$$

$$D_{1}=\frac{1}{h_{1}^{2}}\left(\begin{array}[]{ccccc}-2&1&0&\cdots&1\\ 1&-2&1&\cdots&0\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\cdots&1&-2&1\\ 1&\cdots&0&1&-2\\ \end{array}\right)_{M_{1}\times{M}_{1},}\,D_{2}=\frac{1}{h_{2}^{2}}\left(\begin{array}[]{ccccc}-2&1&0&\cdots&1\\ 1&-2&1&\cdots&0\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\cdots&1&-2&1\\ 1&\cdots&0&1&-2\\ \end{array}\right)_{M_{2}\times{M}_{2}.}$$

$$\displaystyle\big{(}\partial_{f}^{\alpha}u\big{)}^{n}$$

$$u_h^0$$

$$(u_j^n{)}^3={\left((u_{1,j}^n{)}^3,(u_{2,j}^n{)}^3,\cdots ,(u_{M_1-1,j}^n{)}^3\right)}^T\text{for j=0,1,⋯,M2−1.}$$

$$\Vert AV{\Vert }_{\infty }\ge a\Vert V{\Vert }_{\infty },\Vert AV+c(V{)}^3{\Vert }_{\infty }\ge a\Vert V{\Vert }_{\infty }+c\Vert V{\Vert }_{\infty }^3.$$

$$\Vert u^k{\Vert }_{\infty }⩽1\text{for 0⩽k⩽n−1.}$$

$$\displaystyle(\partial_{f}^{\alpha}u)^{n}=A_{0}^{(n)}u^{n}-L^{n-1}\quad\text{where}\quad L^{n-1}:=\sum_{k=1}^{n-1}\big{(}A_{n-k-1}^{(n)}-A_{n-k}^{(n)}\big{)}u^{k}+A_{n-1}^{(n)}u^{0}.$$

$$\displaystyle\big{\|}L^{n-1}\big{\|}_{\infty}\leqslant\sum_{k=1}^{n-1}\big{(}A_{n-k-1}^{(n)}-A_{n-k}^{(n)}\big{)}\|u^{k}\|_{\infty}+A_{n-1}^{(n)}\|u^{0}\|_{\infty}\leqslant A_{0}^{(n)}.$$

Full text

Keywords:   Time-fractional Allen-Cahn equation; fast L1 formula; discrete maximum principle; sharp error estimate; adaptive time-stepping strategy

AMS subject classiffications.   35Q99, 65M06, 65M12, 74A50

1 Introduction

The phase field models have become popular to describe a host of free-boundary problems in various areas, including material, physical and biology systems [1, 2, 3, 4]. Relevant numerical methods and simulations are also increasing substantially [5, 6, 7]. It is well known that the phase field models permit multiple time scales, i.e. an initial dynamics evolves on a fast time scale and later coarsening evolves on a very slow time scale. It is therefore to consider the adaptive time-stepping strategy [8, 9, 10], namely, small time steps are utilized when the energy dissipates rapidly and large time steps are employed otherwise. These works suggest that nonuniform time meshes are preferable in the numerical simulations of phase field models.

In comparison with the bright achievement of classical phase field models, in recent years, there are many researches on building fractional phase field models, such as time, space and time-space fractional Allen-Cahn equations [11, 12, 13, 14, 15, 16] to accurately describe anomalous diffusion problems. Li et al. [14] investigated a space-time fractional Allen-Cahn phase-field model that describes the transport of the fluid mixture of two immiscible fluid phases. They concluded that the alternative model could provide more accurate description of anomalous diffusion processes and sharper interfaces than the classical model. Hou et al. [13] showed that a fractional in space Allen-Cahn equation could be viewed an $L^{2}$ gradient flow for the fractional analogue version of Ginzburg-Landau free energy function. They proved the energy decay property and the maximum principle of continuous problem. Recently, the authors of [11] considered the symmetry analysis, explicit solution and convergence analysis of the time-fractional Allen-Cahn and Klein-Gordon equations with Riemann-Liouville derivative. Zhao et al. [15, 16] studied a series of the time fractional phase field models numerically. The considerable numerical evidences indicate that the effective free energy of the time fractional phase field models obeys a similar power law as the integer ones.

The multi-scale nature of time-fractional phase field models prompts us to construct reliable time-stepping methods on general nonuniform meshes. In this paper, two nonuniform time-stepping schemes are investigated for the time-fractional Allen-Cahn equation [14, 15, 16]

[TABLE]

where $\mathbf{x}=(x,y)^{T}$ and $\Omega=(a,b)\times(c,d)$ with its closure $\bar{\Omega}$. The notation $\partial_{t}^{\alpha}:={}_{0}^{C}\!D_{t}^{\alpha}$ in (1.1) denotes the fractional Caputo derivative of order $\alpha$ with respect to $t$,

[TABLE]

involving the fractional Riemann-Liouville integral $\mathcal{I}_{t}^{\mu}$ of order $\mu>0$, that is,

[TABLE]

The nonlinear bulk force $f(u)=u^{3}-u$, and the small constant $\varepsilon>{0}$, called the interaction length, describes the thickness of the transition boundary between materials. Boundary conditions are set to be periodic so as not to complicate the analysis with unwanted details.

Very recently, the energy decay laws of time-fractional phase field models, involving time-fractional Allen-Cahn equation, time-fractional Cahn-Hilliard equation and time-fractional molecular beam epitaxy models, are reported in [17]. In comparison to the classical physical model, the energy dissipation law of the time-fractional Allen-Cahn equation (1.1) is

[TABLE]

where

[TABLE]

Also, it possesses a maximum principle, namely,

[TABLE]

To our knowledge, there are few results in the literature on the discrete energy decay law or maximum principle of numerical approaches for the time-fractional phase field models, especially on nonuniform time meshes. One of our interests in this paper is to build two nonuniform L1 schemes preserving the maximum principle of the problem (1.1).

We consider the nonuniform time levels $0=t_{0}<t_{1}<\cdots<t_{k-1}<t_{k}<\cdots<t_{N}=T$ with the time-step sizes $\tau_{k}:=t_{k}-t_{k-1}$ for $1\leqslant{k}\leqslant{N}$ and the maximum time-step size $\tau:=\max_{1\leqslant{k}\leqslant{N}}\tau_{k}$. Also, let the local time-step ratio $\rho_{k}:=\tau_{k}/\tau_{k+1}$ and the maximum step ratio $\rho:=\max_{k\geq 1}\rho_{k}$. Given a grid function $\{v^{k}\}$, put $\triangledown_{\tau}v^{k}:=v^{k}-v^{k-1}$, $\partial_{\tau}v^{k-\frac{1}{2}}:=\triangledown_{\tau}v^{k}/\tau_{k}$ and $v^{k-\frac{1}{2}}:=(v^{k}+v^{k-1})/2$ for $k\geq{1}$. Always, let $(\Pi_{1,k}v)(t)$ denote the linear interpolant of a function $v(t)$ at two nodes $t_{k-1}$ and $t_{k}$, and define a piecewise linear approximation

[TABLE]

As an essential mathematical feature of linear and nonlinear subdiffusion problems including the time-fractional Allen-Cahn problem (1.1)-(1.2), the solution always lacks the smoothness near the initial time although it would be smooth away from $t=0$, see [18, 19]. Actually, assuming the nonlinear function $f$ is Lipschitz continuous and the initial data $u^{0}\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$, Jin et al. [19, Theorem 3.1] proved the subdiffusion problem has an unique solution $u$ for which $u\in C\left([0,T];H^{2}(\Omega)\cap H_{0}^{1}(\Omega)\right)$, $\partial_{t}^{\alpha}u\in C\left([0,T];L^{2}(\Omega)\right)$ and $\partial_{t}u\in L^{2}(\Omega)$ with $\|\partial_{t}u(t)\|_{L^{2}(\Omega)}\leq C_{u}t^{\alpha-1}$ for $0<t\leqslant T$. The L1 scheme with a lagging linearized technique for handling the nonlinearity $f(u)$ has been analyzed, and [19, Theorem 4.5] showed that the discrete solution is $O(\tau^{\alpha})$ convergent in $L^{\infty}(L^{2}(\Omega))$. It formally implies that, in any numerical methods for solving time-fractional diffusion equations, a key consideration is the singularity of the solution near the time $t=0$, see also [20, 21, 22]. More directly, we consider the L1 scheme for the time-fractional problem (1.1)-(1.2) describing the coalescence of two kissing bubbles, see more details in Example 4.2. Fig. 1 plots the discrete time derivative $\partial_{\tau}v^{k-\frac{1}{2}}$ near $t=0$ on the graded mesh $t_{k}=(k/N)^{\gamma}$. They suggest that

[TABLE]

It says that the solution possesses weak singularity like $u_{t}=O(t^{\alpha-1})$ near initial time, which can be alleviated by using the graded meshes. Thus the second interest of this paper is to resolve the essentially weak singularity in the equation (1.1) by refining time mesh near $t=0$. Actually, we will show that the graded mesh can recover the optimal time accuracy of L1 formula when the solution $u$ does not have the required regularity.

In the next section, we construct the backward Euler and stabilized semi-implicit schemes by using the nonuniform fast L1 formula $(\partial_{f}^{\alpha}u)^{n}$ described in (2.6). Theorems 2.1 and 2.2 show that both the backward Euler method (2.11)-(2.12) and stabilized semi-implicit method (2.16)-(2.17) preserve the maximum principle (1.7) in the discrete level such that they are unconditionally stable in the maximum norm. By using the recently proposed discrete fractional Grönwall inequality [23] and the global consistency analysis [21] of L1 formula, we prove that, see Theorems 3.1 and 3.2, the fully implicit method (2.11)-(2.12) is convergent with an optimal order of $O(\tau^{2-\alpha})$ and the stabilized scheme (2.16)-(2.17) is convergent with an optimal order of $O(\tau)$ in time on the graded meshes with a grading parameter $\gamma\geq 1$. Unfortunately, we are not able to establish any discrete energy dissipation laws on general nonuniform meshes and leave it as an open problem (see Remark 1).

In summary, the main contributions of this paper are the following: (i) develop two fast L1 time-stepping methods with unequal time-steps preserving the discrete maximum principle, (ii) prove the unconditional convergence with the optimal accuracy in time. Extensive numerical experiments are curried out in section 4 to support our analysis. Some further remarks conclude the article.

2 Fast L1 time-stepping methods

The well-known L1 formula of Caputo derivative (1.3) is given by

[TABLE]

where the corresponding discrete convolution kernels $a_{n-k}^{(n)}$ are defined by

[TABLE]

Obviously, the discrete convolutional kernels $a_{n-k}^{(n)}$ are positive and decreasing, see also [20, 21],

[TABLE]

Note that, this property (2.3) is essential to the preservation of maximum principle for the proposed L1-type schemes described below.

2.1 Fast L1 formula

It is well known that the standard L1 formula (2.1) is prohibitively expensive for long time simulations. Therefore, to reduce the computational cost and storage requirements incurred by employing the L1 formula directly, we apply the sum-of-exponentials (SOE) technique to speed up the evaluation of the original problem. A core result is to approximate the kernel function $t^{-\alpha}$ efficiently on the interval $[\Delta{t},\,T]$, see [24, Theorem 2.5].

Lemma 2.1

For the given $\alpha\in(0,\,1)$, an absolute tolerance error $\epsilon\ll{1}$, a cut-off time $\Delta{t}>0$ and a finial time $T$, there exists a positive integer $N_{q}$, positive quadrature nodes $\theta^{\ell}$ and corresponding positive weights $\varpi^{\ell}\,(1\leqslant{\ell}\leqslant{N_{q}})$ such that

[TABLE]

To be more precise, the Caputo derivative (1.3) is split into the sum of a history part (an integral over $[0,\,t_{n-1}]$) and a local part (an integral over $[t_{n-1},\,t_{n}]$) at the time $t_{n}$. Then, the local part will be approximated by linear interpolation directly, the history part can be evaluated via the SOE technique, that is,

[TABLE]

where $\mathcal{H}^{\ell}(t_{0}):=0$ and $\mathcal{H}^{\ell}(t_{k}):=\int_{0}^{t_{k}}e^{-\theta^{\ell}(t_{k}-s)}v^{\prime}(s)\,\mathrm{d}{s}$. By utilizing the linear interpolation and a recursive formula, we can approximate $\mathcal{H}^{\ell}(t_{k})$ by

[TABLE]

where the positive coefficients

[TABLE]

Having taken this excursion through (2.1)-(2.1), we arrive at the fast algorithm of L1 formula

[TABLE]

in which $\mathcal{H}^{\ell}(t_{k})$ is computed by using the recursive relationship

[TABLE]

For the convenience of numerical analysis, we now eliminate the historic term $\mathcal{H}^{\ell}(t_{k})$ from the fast L1 formula (2.6). From the recursive equation (2.7), a direct calculation yields

[TABLE]

By substituting (2.8) into (2.6), we get an alternative definition

[TABLE]

where the corresponding discrete convolution coefficient $A_{n-k}^{(n)}$ is defined by

[TABLE]

For the discrete kernels $A_{j}^{(n)}$, we have the following result [21, Lemma 2.5].

Lemma 2.2

If the tolerance error $\epsilon$ of SOE satisfies $\epsilon\leqslant\min{\{\frac{1}{3}\omega_{1-\alpha}(T),\,\alpha\omega_{2-\alpha}(1)\}}$, then the discrete convolutional kernel $A_{n-k}^{(n)}$ of (2.10) satisfies

• (i)

$A_{k-1}^{(n)}>A_{k}^{(n)}>0$* for $1\leqslant{k}\leqslant{n-1}$;*

• (ii)

$A_{0}^{(n)}=a_{0}^{(n)}$* and  $A_{n-k}^{(n)}\geq\frac{2}{3}a_{n-k}^{(n)}$ for $1\leqslant{k}\leqslant{n-1}$.*

2.2 Backward Euler scheme

We recall briefly the difference approximation in space. For two positive integers $M_{1},\,M_{2}$, let the spatial lengths $h_{1}:=(b-a)/M_{1},\,h_{2}:=(d-c)/M_{2}$ and $x_{i}=a+ih_{1}$, $y_{j}=c+jh_{2}$ for $0\leqslant i\leqslant M_{1}$, $0\leqslant j\leqslant M_{2}$. Also, denote \bar{\Omega}_{h}:=\big{\{}\mathbf{x}_{h}=(x_{i},y_{j})\,|\,0\leqslant i\leqslant M_{1},0\leqslant j\leqslant M_{2}\} and put $\Omega_{h}:=\bar{\Omega}_{h}\cap\Omega$. For any grid function $\{v_{h}\,|\,\mathbf{x}_{h}\in\bar{\Omega}_{h}\}$, denote a grid function space

[TABLE]

where $v^{T}$ is the transpose of the vector $v$. The maximum norm $\|v\|_{\infty}:=\max_{\mathbf{x}_{h}\in\Omega_{h}}|v_{h}|$.

Let $D_{h}$ be the discrete matrix of Laplace operator $\Delta$ subject to periodic boundary conditions. With the Kronecker tensor product $\otimes$, the matrix $D_{h}=I_{1}\otimes{D_{1}}+D_{2}\otimes{I_{2}}$, in which $I_{1}$ and $I_{2}$ are the identity matrices of order $M_{2}\times{M}_{2}$ and $M_{1}\times{M}_{1}$, respectively, and the matrices $D_{1}$ and $D_{2}$ are of forms

[TABLE]

Then we have some primary properties of the discrete matrix $D_{h}$ in the next lemma, which is straightforward to check and we thus omit the proof here.

Lemma 2.3

Under the periodic boundary condition, the discrete matrix $D_{h}$ of the Laplace operator possesses the following properties

• (a)

The discrete matrix $D_{h}$ is symmetric.

• (b)

For any nonzero $v\in{\mathbb{V}_{h}}$, $v^{T}D_{h}v\leqslant{0}$, i.e., the matrix $D_{h}$ is negative semi-definite.

• (c)

*The elements of $D_{h}=(d_{ij})$ fulfill $d_{ii}=-\max_{i}\sum_{j\neq{i}}|d_{ij}|$ for each $i$. *

Now we have the backward Euler-type scheme on irregular meshes for (1.1)-(1.2),

[TABLE]

where $f(u^{n}):=(u^{n})^{3}-u^{n}$ with the vector $(u^{n})^{3}=\left((u_{1}^{n})^{3},(u_{2}^{n})^{3},\cdots,(u_{M_{2}-1}^{n})^{3}\right)^{T}$ and

[TABLE]

Now we prove that the fully discrete scheme (2.11)-(2.12) preserves the maximum principle numerically. Always, we need the following result [13, Lemma 3.2].

Lemma 2.4

Let $B$ be a real $M\times{M}$ matrix and $A=aI-B$ with $a>0$. If the elements of $B=(b_{ij})$ fulfill $b_{ii}=-\max_{i}\sum_{j\neq{i}}|b_{ij}|$, then for any $c>0$ and $V\in{\mathbb{R}^{M}}$ we have

[TABLE]

Theorem 2.1

If $\|u^{0}\|_{\infty}\leqslant{1}$ and the maximum time-step size $\tau\leqslant 1/\sqrt[\alpha]{\Gamma(2-\alpha)},$ then the solution of backward Euler scheme (2.11)-(2.12) satisfies $\|u^{k}\|_{\infty}\leqslant{1}$ for $0\leqslant{k}\leqslant{N}.$ So it preserves the maximum principle (1.7) numerically and is unconditionally stable.

**Proof ** We use the complete mathematical induction. Obviously, the claimed inequality holds for $k=0$. For $1\leqslant n\leq N$, assume that

[TABLE]

It remains to verify that $\|u^{n}\|_{\infty}\leqslant{1}$. From the definition (2.9), one has

[TABLE]

Thanks to the decreasing property in Lemma 2.2 (i), the induction hypothesis (2.13) and the triangle inequality yield

[TABLE]

Then, from the numerical scheme (2.11), it is easy to obtain

[TABLE]

For the left hand side of (2.15), we apply Lemma 2.3 (c) and Lemma 2.4 to find that

[TABLE]

Then it follows from (2.15) that (A_{0}^{(n)}-1)\big{\|}u^{n}\big{\|}_{\infty}+\big{\|}u^{n}\big{\|}_{\infty}^{3}\leqslant A_{0}^{(n)}. If $A_{0}^{(n)}\geq 1$ or the maximum step size $\tau\leqslant 1/\sqrt[\alpha]{\Gamma(2-\alpha)},$ the above inequality implies $\|u^{n}\|_{\infty}\leqslant 1$ immediately. Otherwise, we have $(A_{0}^{(n)}-1)\|u^{n}\|_{\infty}+\|u^{n}\|_{\infty}^{3}>A_{0}^{(n)},$ because the function

[TABLE]

is monotonically increasing for any $z>0$. This leads to a contradiction and then the claimed result holds for $k=n$. The principle of induction completes the proof.

2.3 Stabilized semi-implicit scheme

The backward Euler scheme (2.11)-(2.12) is a fully nonlinear implicit scheme and some inner iteration will be needed. To accelerate the time-stepping process, we build a linearized scheme here by using the well-known stabilized technique via a stabilized term $S(u^{n}-u^{n-1})$ for a properly large scalar parameter $S>0$ , see also the recent work [17]. The resulting stabilized semi-implicit scheme for the problem (1.1)-(1.2) reads

[TABLE]

We have the following result on discrete maximum principle and stability.

Theorem 2.2

If $\|u^{0}\|_{\infty}\leqslant{1}$ and the scalar stabilized parameter $S\geq 2$, then the solution of stabilized semi-implicit scheme (2.16)-(2.17) satisfies

[TABLE]

So it preserves the maximum principle (1.7) numerically and is unconditionally stable.

**Proof ** It only needs to verify that $\|u^{n}\|_{\infty}\leqslant{1}$ under the induction hypothesis

[TABLE]

From the linearized scheme (2.16), one has

[TABLE]

where $L^{n-1}$ is defined in (2.14). Thanks to the decreasing property in Lemma 2.2 (i), the induction hypothesis and the triangle inequality yield \big{\|}L^{n-1}\big{\|}_{\infty}\leqslant A_{0}^{(n)}. Furthermore, it is easy to check that

[TABLE]

thus the right hand side of (2.18) can be bounded by

[TABLE]

For the left hand side of (2.18), we apply Lemma 2.3 (c) and Lemma 2.4 to find that

[TABLE]

Then the desired estimate $\|u^{n}\|_{\infty}\leqslant{1}$ follows from (2.18) directly.

Due to the presence of the stabilized term $S(u^{n}-u^{n-1})$, the numerical solution generated by the semi-implicit scheme (2.16)-(2.17) will be limited to first-order accurate in time even if the solution is sufficiently smooth. We address the error analysis in the next section.

3 Global consistency analysis and convergence

To facilitate the error analysis of difference approximations in space, we assume that the continuous solution $u$ is sufficiently smooth in space and satisfies

[TABLE]

where a regularity parameter $\sigma\in(0,1)$ is introduced to make our analysis extendable.

In [21], the local consistency error $\Upsilon^{j}:=(\partial_{t}^{\alpha}u)(t_{j})-(\partial_{f}^{\alpha}u)^{j}$ of fast L1 formula (2.9) was bounded by a discrete convolution structure, which is valid for any time meshes. It provides us an opportunity to give the global error via the global consistency error \sum_{j=1}^{n}p_{n-j}^{(n)}\big{|}\Upsilon^{j}\big{|}, where $p_{n-j}^{(n)}$ are the discrete complementary convolution kernels defined via (A.1). Note that, the definition (2.10) and Lemma 2.2 (i) show that the discrete convolutional kernels $A_{n-k}^{(n)}$ fulfill two assumptions Ass1-Ass2 in Appendix A with $\pi_{a}=\frac{3}{2}$. In this section, we will use the results of Lemma A.1 without further declarations.

Lemma 3.1

Under the condition of Lemma 2.2, the global consistency error is bounded by

[TABLE]

where the local quantities $G^{k}:=2\int_{t_{k-1}}^{t_{k}}\left(t-t_{k-1}\right)\left|u_{tt}\right|\,\mathrm{d}t$ for $1\leqslant k\leqslant n$ and $\hat{t}_{n}:=\max\{1,t_{n}\}$.

**Proof ** On the basis of the upper bound of $(\partial_{t}^{\alpha})u(t_{n})-(\partial_{\tau}^{\alpha}u)^{n}$ given in [21, Lemma 3.1], the estimate (3.5) in the proof of [21, Lemma 3.3] gives the desired result.

To resolve such a solution $u$ efficiently, it is appropriate to choose the time mesh such that the following condition [20, 21, 22, 25] holds.

1. AssG.

Let $\gamma\geq{1}$ be a user-chosen parameter. There is a mesh-independent constant $C_{\gamma}>0$ such that $\tau_{k}\leq C_{\gamma}\tau\min\{1,t_{k}^{1-1/\gamma}\}$ for $1\leq k\leq N$ and $t_{k}\leqslant C_{\gamma}t_{k-1}$ for $2\leqslant{k}\leqslant{N}$.

Here, the parameter $\gamma\geq 1$ controls the extent to which the time levels are concentrated near $t=0$. If the mesh is quasi-uniform, then AssG holds with $\gamma=1$. As $\gamma$ increases, the initial step sizes become smaller compared to the later ones. A simple example of a family of meshes satisfying AssG is the graded mesh $t_{k}=T(k/N)^{\gamma}$ with the maximum step ratio $\rho=1$.

It is to note that, the global consistency error in Lemma 3.1 gives a superconvergence estimate of nonuniform L1 formula. Consider the first time level $n=1$, the regularity setting (3.1) gives $\left|\Upsilon^{1}\right|\leqslant C_{u}A_{0}^{(1)}\int_{0}^{t_{1}}t^{\sigma-1}\,\mathrm{d}t\leqslant C_{u}\tau_{1}^{\sigma-\alpha}/\sigma,$ implying that the L1 formula is always inconsistent if $0<\sigma\leq\alpha$, also see Table 1 in Section 4. However, we have the global consistency error of order $O(\tau_{1}^{\sigma})$, because $p_{0}^{(1)}\left|\Upsilon^{1}\right|\leq G^{1}\leq C_{u}\tau_{1}^{\sigma}/\sigma$. In general, we have the following result from [21, Lemma 3.3].

Corollary 3.1

Under the regularity (3.1), the global consistency error can be bounded by

[TABLE]

Specifically, if the mesh satisfies AssG, then

[TABLE]

Theorem 3.1

Assume that $\|u_{0}\|_{L^{\infty}(\Omega)}\leqslant{1}$ and the solution of (1.1)-(1.2) satisfies the regular assumption (3.1). If the maximum step size $\tau\leqslant{1}/\sqrt[\alpha]{6\Gamma(2-\alpha)}$, then the numerical solution $u_{h}^{n}$ of the backward Euler scheme (2.11)-(2.12) is convergent in the maximum norm, that is,

[TABLE]

for $1\leqslant{n}\leqslant{N}.$ Moreover, when the time mesh satisfies AssG, it holds that

[TABLE]

which achieves the optimal accuracy $O(\tau^{2-\alpha})$ if the graded parameter $\gamma\geq\max{\{1,\,(2-\alpha)/\sigma\}}$.

**Proof ** Let $U_{h}^{n}:=u(\mathbf{x}_{h},t_{n})$ and the error function $e_{h}^{n}:=U_{h}^{n}-u_{h}^{n}\in{\mathbb{V}_{h}}$ for $\mathbf{x}_{h}\in\bar{\Omega}_{h}$ and $0\leqslant{n}\leqslant{N}$. It is easy to find that the exact solution $U_{h}^{n}$ satisfies the governing equations

[TABLE]

where $\left(R_{t}\right)^{n}$ and $\left(R_{s}\right)^{n}$ denote the truncation errors in time and space, respectively. Subtracting (2.11)-(2.12) from the above two equations, respectively, one gets

[TABLE]

Recalling the elementary inequality $|(a^{3}-a)-(b^{3}-b)|\leqslant{2}|a-b|$ for $\forall\,a,b\in[-1,1]$, we apply Theorem 2.1 (discrete maximum principle) to get

[TABLE]

Thus the triangle inequality with the error equation (3.2) gives

[TABLE]

Applying the decreasing property $(i)$ of the kernels $A_{n-k}^{(n)}$ and the triangle inequality, we can bound the left hand side of (3.4) by

[TABLE]

where Lemma 2.3 (c) and Lemma 2.4 have been used. Then it follows from (3.4) that

[TABLE]

which takes the form of (A.4) with the substitutions $\lambda:=2$, v^{k}:=\big{\|}e^{k}\big{\|}_{\infty}, \xi^{n}:=\big{\|}\left(R_{t}\right)^{n}\big{\|}_{\infty} and \eta^{n}:=\big{\|}\left(R_{s}\right)^{n}\big{\|}_{\infty}. Lemma A.1 (the discrete fractional Grönwall inequality) says that, if the maximum step size $\tau\leqslant{1}/\sqrt[\alpha]{6\Gamma(2-\alpha)}$, then it holds that

[TABLE]

Then Corollary 3.1 yields the claimed estimate and completes the proof.

For the semi-implicit scheme (2.16)-(2.17), the global error is dominated by the stabilized term $S(u^{k}-u^{k-1})$. Under the regular assumption (3.1), the local consistency error is about $\int_{t_{k-1}}^{t_{k}}\left|u_{t}\right|\,\mathrm{d}{t}$. One can follow the proof of [21, Lemma 3.3] to bound the corresponding global error as follows (also see the case of $m=0$ in the estimate (A.3))

[TABLE]

Then, a similar proof of Theorem 3.1 leads to the following result.

Theorem 3.2

Assume that $\|u_{0}\|_{L^{\infty}(\Omega)}\leqslant{1}$ and the exact solution of (1.1)-(1.2) satisfies the regular assumption (3.1). If the stabilized parameter $S\geq 2$ and the maximum time-step size $\tau\leqslant{1}/\sqrt[\alpha]{6\Gamma(2-\alpha)}$, then the numerical solution $u_{h}^{n}$ of the semi-implicit scheme (2.16)-(2.17) is convergent in the maximum norm, that is,

[TABLE]

for $1\leqslant{n}\leqslant{N}$. Moreover, when the time mesh satisfies AssG, it holds that

[TABLE]

which achieves the optimal accuracy $O(\tau)$ if the graded parameter $\gamma\geq\max{\{1,\,1/\sigma\}}$.

Remark 1

*(An open problem) It is interesting to mention that, on the uniform mesh, the discrete L1 kernels (2.2) reads *

[TABLE]

the semi-implicit stabilized scheme (2.16)-(2.17) using the L1 formula inherits a discrete energy dissipation law, see [17, Theorem 3.1] for details. As seen, the proof of discrete energy dissipation law relies on the property of a quadratic form $\sum_{k=1}^{n}w_{k}\sum_{j=1}^{k}a_{k-j}w_{j}\geq 0$. However, it seems rather difficult to extend the positive semi-definite property to a general class of nonuniform meshes. More precisely, we are not able to verify the positive semi-definite property of the following quadratic form (by taking $w_{k}=\triangledown_{\tau}v^{k}$)

[TABLE]

More generally, it has yet to be determined what restrictions must be imposed on the discrete convolution coefficients $\{A_{n-k}^{(n)}\,|1\leqslant k\leqslant n\}$ so that the quadratic form $\sum_{k=1}^{n}w_{k}\sum_{j=1}^{k}A_{k-j}^{(k)}w_{j}$ is positive semi-definite. This problem could be challenging and remains open to us.

4 Numerical examples

The nonuniform fast L1 time-stepping methods (2.11)-(2.12) and (2.16)-(2.17) are examined for solving the Allen-Cahn problem (1.1)-(1.2). Always, we set the absolute tolerance error $\epsilon=10^{-12}$ for the SOE approximation. The second-order centered difference scheme is used to approximate the Laplace operator with the same length $h=1/M$ in each spatial direction. For the nonlinear scheme (2.11)-(2.12), a simple iteration is employed to solve the nonlinear algebra equations at each time level with the termination error $10^{-12}$. The maximum norm error $e(M,N):=\max_{1\leqslant{n}\leqslant{N}}\|U^{n}-u^{n}\|_{\infty}$ is recorded in each run, and the experimental convergence order in time is computed by

[TABLE]

where $\tau(N)$ denotes the maximum time-step size for total $N$ subintervals.

Example 4.1

To examine the temporal accuracy of our time-stepping schemes, consider the time-fractional Allen-Cahn equation $\partial_{t}^{\alpha}u=\frac{1}{8\pi^{2}}\Delta{u}-f(u)+g(\mathbf{x},t)$ for $\mathbf{x}\in(0,1)^{2}$ and $0<t<1$ such that it has an exact solution $u=\omega_{1+\sigma}(t)\sin(2\pi{x})\sin(2\pi{y})$.

The time interval $[0,T]$ is always divided into two parts $[0,T_{0}]$ and $[T_{0},T]$ with total $N$ subintervals. We will take $T_{0}=\min\{1/\gamma,T\}$, and apply the graded grid $t_{k}=T_{0}(k/N_{0})^{\gamma}$ in $[0,T_{0}]$ to resolve the initial singularity. In the remainder interval $[T_{0},T]$, we put $N_{1}:=N-N_{0}$ cells with random time-steps

[TABLE]

where $\epsilon_{k}\in(0,1)$ are the random numbers.

We take the spatial grid points $M=1024$ in each direction such that the temporal error dominates the spatial error in each run. Numerical results of the backward Euler scheme (2.11)-(2.12) for two different cases $\sigma=\alpha$ and $\sigma<\alpha$ are listed in Tables 1-2, respectively. They suggest the time accuracy is of order $O(\tau^{\min\{\gamma\sigma,2-\alpha\}})$ and confirm Theorem 3.1 experimentally. We also run the stabilized semi-implicit scheme (2.16)-(2.17) by setting a variety of regularity parameters. Tables 3-4 report the numerical results in the case $\sigma=\alpha$ and a worse case of $\sigma<\alpha$. It seen that it is accurate of order $O(\tau^{\min\{\gamma\sigma,1\}})$ on the graded meshes, confirming Theorem 3.2 experimentally.

Example 4.2** (Coalescence of two kissing bubbles)**

Consider the time-fractional Allen-Cahn problem (1.1)-(1.2) describing the coalescence of two kissing bubbles inside the spatial domain $\Omega=(-\pi,\,\pi)^{2}$, by taking $\varepsilon=0.1$ and the initial data

[TABLE]

This example is used to examine the physical effect of the fractional order $\alpha$ in the original problem and the physical property of our suggested methods. Theorems 2.1 and 2.2 suggest that variable time-steps are always allowed in our time-stepping approaches. As a matter of fact, the temporal evolution of phase models involve multiple time scales which initial data evolves on a fast time scale at the early stage of dynamics and then the coarsening evolves rather slowly until it reaches a steady state. Hence, to capture the fast dynamics and reduce the cost of computation, we adapt the variant adaptive time-stepping strategy [10]

[TABLE]

where the constant 1 is set to avoid the possible singularity as the model reaches the steady state. The parameters $tol$ and $\beta$ are used to adjust the level of adaptively and would be chosen in experience. A small $tol$ or a large $\beta$ will generate time steps close to $\tau_{\min}$, which a large $tol$ or a small $\beta$ will give time steps close to $\tau_{\max}$. The problem is simulated to the final time $T=100$ by taking $M=128,T_{0}=0.1$, $t_{k}=T_{0}(k/N_{0})^{\gamma}$ with the graded parameter $\gamma=3$ in the initial interval $[0,T_{0}]$ and adopting adaptive time steps in the remainder interval.

We find that the solution profiles, generated by the backward Euler scheme (2.11)-(2.12) with $\tau_{\min}=\tau_{N_{0}}=0.001,\tau_{\max}=0.1,tol=0.15,\beta=200$ and stabilized semi-implicit scheme (2.16)-(2.17) with $S=0.1,\tau_{\min}=\tau_{N_{0}}=0.001,\tau_{\max}=1,tol=1.5,\beta=200$ in the remainder interval, are quite identical. Fig. LABEL:SnapshotsNumAPPL gathers some snapshots at four different times. It is seen that the two bubbles coalesce into a single bubble as the time escapes, while the rate of coalescence is deeply affected by the fractional order $\alpha$, see [15, 16]. The larger the fractional order $\alpha$, the faster the bubbles coalesce.

Fig.2 depicts the solution in the maximum norm and the discrete energy ($E^{n}$ is the discrete counterpart of the energy functional defined in the model in spite of no theoretical proof is available in current work) of the backward Euler scheme (2.11)-(2.12). It is obvious that the solutions are uniformly bounded by the value 1 for different fractional orders $\alpha$, as predicted by Theorem 2.1. Moreover, the larger the fractional order $\alpha$, the faster it approaches the maximum value. The middle of Fig. 2 says that the discrete energy is also decreasing as the time escapes, although we can not verify it theoretically. The right side of Fig. 2 depicts the adopted time-steps, and we observe that the time-steps are always small at the early stage, implying the fast evolution dynamics near the initial time. Fig.3 shows analogous plots for the stabilized scheme (2.16)-(2.17), where we see the similar behaviors on the maximum norm value, the discrete energy and the adaptive time-steps. Note that the maximum time step $\tau_{\max}=0.1$ of the backward Euler scheme is to ensure the convergence of iterative method, thus we can expect the stabilized scheme to be more efficient than the nonlinear one.

5 Concluding remarks

In simulating the time-fractional phase field equations including the Allen-Cahn equation considered in this paper, the initial singularity should be treated properly because it always destroys the time accuracy of numerical algorithms especially near the initial time. We consider two fast L1 time-stepping methods on a general class of nonuniform time meshes such that they will be suitable for both the refined mesh near $t=0$ and certain adaptive time-stepping strategy to resolve the multiple time scales away from $t=0$.

We show that the nonuniform fast L1 formula can be employed to construct some time-stepping methods preserving the discrete maximum principle by virtue of the uniform monotonicity of discrete kernels. By using the discrete fractional Grönwall inequality and global consistency analysis, we established obtain sharp maximum norm error estimates of proposed schemes and validated them numerically.

It seems challenging to build time-stepping approaches maintaining the discrete energy dissipation law, especially on general nonuniform time meshes. Nonetheless, the energy stable schemes permitting adaptive time-stepping strategies are very attractive because they would be applicable for other time-fractional phase-field models and for long-time simulations approaching the steady state. These issues will be addressed in the forthcoming reports.

Acknowledgements

The authors would like to thank Prof. Jia Zhao and Prof. Yuezheng Gong for their valuable discussions and fruitful suggestions.

Appendix A Discrete fractional Grönwall lemma

The recently developed discrete fractional Grönwall inequality in [23] is applicable for any nonuniform time meshes and suitable for a variety of discrete fractional derivatives. The following lemma, involving the Mittag–Leffler function $E_{\alpha}(z):=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(1+k\alpha)}$, gathers three previous (slightly simplified) results from [23, Lemma 2.2, Theorems 3.1 and 3.2].

Lemma A.1

*For $n=1,2,\cdots,N$, assume that the discrete convolution kernels $\{A_{n-k}^{(n)}\}_{k=1}^{n}$ satisfy the following two assumptions:

Ass1. There is a constant $\pi_{a}>0$ such that $A_{n-k}^{(n)}\geq\frac{1}{\pi_{a}}\int_{t_{k-1}}^{t_{k}}\frac{\omega_{1-\alpha}(t_{n}-s)}{\tau_{k}}\,\mathrm{d}s$ for $1\leqslant k\leqslant n$.

Ass2. The discrete kernels are monotone, i.e. $A_{n-k-1}^{(n)}-A_{n-k}^{(n)}\geq 0$ for $1\leqslant k\leqslant n-1$.

Define also a sequence of discrete complementary convolution kernels $\{p_{n-j}^{(n)}\}_{j=1}^{n}$ by*

[TABLE]

Then the discrete complementary kernels $p^{(n)}_{n-j}\geq 0$ are well-defined and fulfill

[TABLE]

Suppose that the offset parameter $0\leqslant\nu<1$, $\lambda$ is a non-negative constant independent of the time-steps and the maximum step size $\tau\leq 1/\sqrt[\alpha]{2\Gamma(2-\alpha)\lambda\pi_{a}}.$ If the non-negative sequences $(v^{k})_{k=0}^{N}$, $(\xi^{k})_{k=1}^{N}$ and $(\eta^{k})_{k=1}^{N}$ satisfy

[TABLE]

or

[TABLE]

then it holds that, for $1\leq n\leq N$,

[TABLE]