A high-order discrete energy decay and maximum-principle preserving scheme for time fractional Allen-Cahn equation

TL;DR

This paper introduces a novel finite difference scheme using the shifted fractional trapezoidal rule for the time-fractional Allen-Cahn equation, ensuring energy decay and maximum-principle preservation, with demonstrated accuracy and efficiency.

Contribution

The paper develops a new SFTR-based scheme that preserves key physical properties and effectively handles initial singularities in the time-fractional Allen-Cahn equation.

Findings

The scheme guarantees discrete energy decay.

It preserves the maximum principle.

Numerical results outperform classical methods.

Abstract

The shifted fractional trapezoidal rule (SFTR) with a special shift is adopted to construct a finite difference scheme for the time-fractional Allen-Cahn (tFAC) equation. Some essential key properties of the weights of SFTR are explored for the first time. Based on these properties, we rigorously demonstrate the discrete energy decay property and maximum-principle preservation for the scheme. Numerical investigations show that the scheme can resolve the intrinsic initial singularity of such nonlinear fractional equations as tFAC equation on uniform meshes without any correction. Comparison with the classic fractional BDF2 and L2-1$_\sigma$ method further validates the superiority of SFTR in solving the tFAC equation. Experiments concerning both discrete energy decay and discrete maximum-principle also verify the correctness of the theoretical results.

Tables

Table 1. Table 1: Results for smooth solutions in Example 2 .

$\alpha$	$\tau$	SFTR-$\frac{1}{2}$			F-BDF2
		Errer	Rates		Errer	Rates
0.3	1/20	5.9056E-04			8.6484E-05
	1/40	1.6353E-04	1.85		2.3802E-05	1.86
	1/80	4.3146E-05	1.92		6.0214E-06	1.98
	1/160	1.1200E-05	1.95		1.2785E-06	2.24
0.6	1/20	1.5359E-04			1.5720E-04
	1/40	4.0232E-05	1.93		4.2323E-05	1.89
	1/80	1.0320E-05	1.96		1.0950E-05	1.95
	1/160	2.6647E-06	1.95		2.6695E-06	2.04
0.9	1/20	4.0839E-05			2.2212E-04
	1/40	1.0424E-05	1.97		5.9128E-05	1.91
	1/80	2.6237E-06	1.99		1.5253E-05	1.95
	1/160	6.7011E-07	1.97		3.8653E-06	1.98

Table 2. Table 2: Comparisons of our scheme ( 21 ) using SFTR- 1 2 1 2 \frac{1}{2} with F-BDF2 for Example 3 .

$\alpha$	$\tau$	SFTR-$\frac{1}{2}$			F-BDF2
		Errer	Rates		Errer	Rates
0.3	1/20	1.7143E-04			4.0081E-02
	1/40	5.7584E-05	1.57		2.7325E-02	0.55
	1/80	1.8625E-05	1.63		1.7007E-02	0.68
	1/160	5.1691E-06	1.85		8.6550E-03	0.97
0.6	1/20	7.7902E-05			1.5362E-02
	1/40	1.8377E-05	2.08		9.1160E-03	0.75
	1/80	4.2841E-06	2.10		4.9807E-03	0.87
	1/160	8.8900E-07	2.27		2.2460E-03	1.15
0.9	1/20	2.7149E-05			4.6119E-03
	1/40	6.6695E-06	2.03		2.3217E-03	0.99
	1/80	1.6462E-06	2.02		1.0908E-03	1.09
	1/160	3.5626E-07	2.21		4.2951E-04	1.34

Table 3. Table 3: Comparisons of SFTR- 1 2 1 2 \frac{1}{2} with L2-1 σ for Example 4 .

	$\alpha =0.2$				$\alpha =0.4$
	$N=32$	$N=64$	$N=128$		$N=32$	$N=64$	$N=128$
SFTR-$\frac{1}{2}$	7.00E-05	2.73E-05	9.01E-06		3.14E-05	9.11E-06	2.48E-06
L2-1σ(uniform)	5.31E-04	2.40E-04	9.68E-05		7.83E-04	3.55E-04	1.43E-04
L2-1σ($\min$)	2.55E-05	5.67E-06	1.22E-06		6.79E-05	1.58E-05	3.45E-06
	$\alpha =0.6$				$\alpha =0.8$
	$N=32$	$N=64$	$N=128$		$N=32$	$N=64$	$N=128$
SFTR-$\frac{1}{2}$	1.37E-05	3.42E-06	7.70E-07		8.66E-06	2.10E-06	5.28E-07
L2-1σ(uniform)	7.31E-04	3.31E-04	1.33E-04		4.40E-04	1.99E-04	8.05E-05
L2-1σ($\min$)	8.78E-05	2.09E-05	4.34E-06		7.16E-05	1.46E-05	2.56E-06

Equations

$$\begin{array}{r}{\partial }_t^{\alpha }u(\mathbf{x},t)=-\frac{\delta E}{\delta u},\mathbf{x}\in \Omega ,t\in (0,T],\end{array}$$

$$\begin{split}E[u](t)=\int_{\Omega}\bigg{(}\frac{\varepsilon^{2}}{2}|\nabla u|^{2}+F(u)\bigg{)}\mathrm{d}\boldsymbol{x},\end{split}$$

$$\begin{array}{r}{\partial }_t^{\alpha }u(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{u^{\prime }(s)}{(t-s{)}^{\alpha }}\mathrm{d}s,\Gamma (\alpha )\text{ is the Gamma function}.\end{array}$$

$$\begin{array}{rl}{\partial }_t^{\alpha }u(\mathbf{x},t)& ={\epsilon }^2\Delta u-f(u),\mathbf{x}\in \Omega ,t\in (0,T],u(\mathbf{x},0)=u_0(\mathbf{x}),\mathbf{x}\in \overline{\Omega },\end{array}$$

$$\begin{split}E^{\prime}[u](t)+\bigg{\|}\frac{\delta E}{\delta u}\bigg{\|}^{2}=0,\quad\text{or}\quad E[u](t_{1})\leq E[u](t_{2}),\quad\text{if }t_{1}\geq t_{2},\end{split}$$

$$\begin{array}{r}{E}_w[u](t)={\int }_0^{1}w(s)E[u](st)\mathrm{d}s,\end{array}$$

$$\begin{split}E_{c}[u]=E[u]+\frac{1}{2}\mathcal{I}^{\alpha}_{t}\bigg{\|}\frac{\delta E}{\delta u}\bigg{\|}^{2},\end{split}$$

$$\begin{split}E^{\prime}_{c}[u](t)+\frac{t^{\alpha-1}}{2\Gamma(\alpha)}\bigg{\|}\frac{\delta E}{\delta u}\bigg{\|}^{2}\leq 0\quad\text{for }t>0,\end{split}$$

$$\begin{array}{r}|u(\mathbf{x},t)|\le 1,\text{if }|u_0(\mathbf{x})|\le 1.\end{array}$$

$$\begin{array}{r}{\partial }_{\tau }^{\alpha ,n-\frac{1}{2}}\varphi ={\tau }^{-\alpha }\sum _{m=0}^n{\omega }_m({\varphi }^{n-m}-{\varphi }^0),\end{array}$$

$$\omega(\xi)=\sum_{m=0}^{\infty}\omega_{m}\xi^{m}=\bigg{[}\frac{1-\xi}{\frac{1}{2}(1+\xi)+\frac{1}{2\alpha}(1-\xi)}\bigg{]}^{\alpha},$$

$$\begin{split}\omega_{0}&=\bigg{(}\frac{2\alpha}{\alpha+1}\bigg{)}^{\alpha},\quad\omega_{1}=-\alpha\bigg{(}\frac{2\alpha}{\alpha+1}\bigg{)}^{\alpha+1},\\ \omega_{m}&=\frac{2\alpha}{m(\alpha+1)}\bigg{\{}\bigg{[}\frac{1}{\alpha}(m-1)-\alpha\bigg{]}\omega_{m-1}+\frac{\alpha-1}{2\alpha}(m-2)\omega_{m-2}\bigg{\}},m\geq 2.\end{split}$$

$$\begin{split}\Phi(x,y)=\frac{1}{y}(y-1-x^{2})\bigg{(}1-\frac{x+1}{2}\frac{y-1}{y}\bigg{)}-\frac{(1-x)(y-1)}{2(y-x^{2})}\bigg{(}1+\frac{1}{y}\bigg{)}.\end{split}$$

$$\begin{array}{r}\tilde{\Phi }(x,\tilde{y})=\frac{(1-x)(9-{\tilde{y}}^2)}{6(x^2\tilde{y}-3)}-\frac{1}{18}(x^2\tilde{y}+\tilde{y}-3)(3-3x+\tilde{y}+x\tilde{y}),\end{array}$$

$$\begin{split}\omega_{2}=\frac{\alpha}{\alpha+1}\bigg{(}\frac{1}{\alpha}-\alpha\bigg{)}\omega_{1}<0.\end{split}$$

$$\frac{\omega_{m-1}}{\omega_{m-2}}\geq\sigma_{m-1}>\frac{\frac{1-\alpha}{2\alpha}(m-2)}{\frac{m-1}{\alpha}-\alpha},\quad\text{where~{}}\sigma_{m}=\frac{\frac{1-\alpha}{2\alpha}(m-1)}{\frac{m}{\alpha}-\alpha}\frac{2}{1+\alpha}\big{(}1+\frac{1}{m}\big{)},$$

$$\omega_{m}=\frac{2\alpha\omega_{m-2}}{m(\alpha+1)}\bigg{[}\frac{1}{\alpha}(m-1)-\alpha\bigg{]}\bigg{[}\frac{\omega_{m-1}}{\omega_{m-2}}-\frac{\frac{1-\alpha}{2\alpha}(m-2)}{\frac{1}{\alpha}(m-1)-\alpha}\bigg{]}<0$$

$$\frac{{\omega }_2}{{\omega }_1}=1-\alpha \ge {\sigma }_2=\frac{3}{2}\frac{1-\alpha }{(2-{\alpha }^2)(1+\alpha )}.$$

$$\begin{split}\frac{\omega_{m}}{\omega_{m-1}}&\geq\frac{2\alpha}{m(\alpha+1)}\bigg{(}\frac{m-1}{\alpha}-\alpha-\frac{1-\alpha}{2\alpha}\frac{m-2}{\sigma_{m-1}}\bigg{)}\\ &=\frac{2\alpha}{m(\alpha+1)}\bigg{(}\frac{m-1}{\alpha}-\alpha\bigg{)}\bigg{(}1-\frac{1+\alpha}{2}\frac{m-1}{m}\bigg{)},\end{split}$$

$$\begin{split}\frac{\omega_{m}}{\omega_{m-1}}&=\frac{2\alpha}{m(\alpha+1)}\bigg{[}\frac{1}{\alpha}(m-1)-\alpha-\frac{1-\alpha}{2\alpha}(m-2)\frac{\omega_{m-2}}{\omega_{m-1}}\bigg{]}\\ &\leq\frac{2\alpha}{m(\alpha+1)}\bigg{[}\frac{1}{\alpha}(m-1)-\alpha-\frac{1-\alpha}{2\alpha}(m-2)\bigg{]}\\ &=1-\frac{2\alpha}{m}<1.\end{split}$$

$$\Vert v{\Vert }_{\infty }=\underset{1\le j,k\le M}{\max }|v_{j,k}|,\Vert v{\Vert }^2=h^2\sum _{1\le j,k\le M}|v_{j,k}{|}^2.$$

$$(u,v)=h^2\sum _{1\le j,k\le M}{u}_{j,k}{v}_{j,k}.$$

$$\begin{array}{rl}{\delta }_x^2{\psi }_{j,k}=\frac{{\psi }_{j+1,k}-2{\psi }_{j,k}+{\psi }_{j-1,k}}{h^2},& {\delta }_y^2{\psi }_{j,k}=\frac{{\psi }_{j,k+1}-2{\psi }_{j,k}+{\psi }_{j,k-1}}{h^2}.\end{array}$$

$$\begin{array}{rl}{\delta }_x^2{\psi }_{j,k}& ={\partial }_x^2{\psi }_{j,k}+O(h^2),{\delta }_y^2{\psi }_{j,k}={\partial }_y^2{\psi }_{j,k}+O(h^2).\end{array}$$

$$\begin{array}{r}{\tau }^{-\alpha }\sum _{m=0}^n{\omega }_m(U_{j,k}^{n-m}-U_{j,k}^0)={\epsilon }^2({\delta }_x^2+{\delta }_y^2)U_{j,k}^{n-\frac{1}{2}}-f_{j,k}^{n-1,n},\end{array}$$

$$\begin{array}{r}{f}_{j,k}^{n-1,n}=\frac{1}{3}(U_{j,k}^n{)}^3+\frac{1}{2}(U_{j,k}^{n-1}{)}^2(U_{j,k}^n)+\frac{1}{6}(U_{j,k}^{n-1}{)}^3-\frac{1}{2}(U_{j,k}^n+U_{j,k}^{n-1}),\end{array}$$

$$\begin{array}{r}(U_{j,k}^n-U_{j,k}^{n-1})f_{j,k}^{n-1,n}\ge F(U_{j,k}^n)-F(U_{j,k}^{n-1}).\end{array}$$

$$\mathbf{v}=\mathrm{vec}(v):=[v_{1,1},v_{2,1},\cdots ,v_{M,1},v_{1,2},\cdots ,v_{M,2},\cdots ,v_{1,M},\cdots ,v_{M,M}{]}^T,$$

$$\begin{array}{r}{\tau }^{-\alpha }\sum _{m=0}^n{\omega }_m({\mathbf{U}}^{n-m}-{\mathbf{U}}^0)={\epsilon }^2{\Delta }_h{\mathbf{U}}^{n-\frac{1}{2}}-{\mathbf{f}}^{n-1,n},n\ge 1,\end{array}$$

$$D={\left(\begin{array}{ccccc}-2& 1& 0& \cdots & 1\\ 1& -2& 1& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 1& -2& 1\\ 1& \cdots & 0& 1& -2\end{array}\right)}_{M\times M}.$$

Full text

1 Introduction

For systems driven by dissipation of free energy, the gradient flow is a useful tool in diverse research areas (see, e.g., AndersonMcFaddenWheeler ; AllenCahn ; ShaoRappelLevine ; XuPrinceSnakes ; LiuChengWang ). In this work, we are interested in systems where dissipation of the free energy has memory effect and thus can be simulated by the $\alpha$th order ($0<\alpha<1$) fractional gradient flow (see HouZhuXu ; DuYangZhou ; TangYuZhou ),

[TABLE]

with the energy

[TABLE]

where $T>0$, $\Omega=(a,b)\times(a,b)$, $F(u)=\frac{1}{4}(1-u^{2})^{2}$ is a double-well potential, and the constant $\varepsilon(0<\varepsilon\ll 1)$ represents the thickness of the transition boundary between materials. The symbol $\partial_{t}^{\alpha}$ denotes the Caputo fractional operator known as

[TABLE]

Imposed the periodic boundary condition, equation (1) is then essentially the nonlocal phase field model, or specifically, the time-fractional Allen-Cahn (tFAC) equation

[TABLE]

where $f(u)=F^{\prime}(u)=u^{3}-u$ represents the nonlinear bulk force. The nonlocality of the fractional derivative and nonlinearity of equation (3) generally render solving this problem rather difficult by analytic methods. We thus consider numerical methods for this problem. In this context, an important issue is to investigate whether numerical methods can inherit some intrinsic properties of the solution, including the energy decay and maximum-principle property.

It is well known that for classical Allen-Cahn equation (with $\alpha=1$), the energy $E[u](t)$ admits a monotonic property:

[TABLE]

while for general cases (with $\alpha\in(0,1)$), the monotonicity (4) still remains an open problem, despite the fact that all numerical tests published in the literature up to now suggest its correctness. To circumvent this difficulty, on the one hand a “weak version” of $(\ref{Int.1.1})$ called energy dissipation, i.e. $E[u](t)\leq E[u](0)$ for $t\geq 0$, was studied in both continuous and discrete levels by Tang et al.TangYuZhou for time-fractional phase-field models (e.g., the time-fractional Allen-Cahn/Cahn-Hilliard/molecular-beam-epitaxy equation). See also ChenZhangZhaoCaoWangZhang . On the other hand, some modified types of energy have been proposed for the model (3), in which the variable $\alpha$ is involved. Specifically, the authors in QuanTangYang defined a weighted energy

[TABLE]

where $w(s)$ denotes some weight function satisfying $\int_{0}^{1}w(s)\mathrm{d}s=1$, and proved $E^{\prime}_{w}[u](t)\leq 0$ provided $w(s)s^{1-\alpha}(1-s)^{\alpha}$ is nonincreasing w.r.t. $s$. It is notable that the energy $E_{w}[u](t)$ involves $\alpha$ since $w(s)$ in general depends on $\alpha$. Considering the tFAC equation degrades to the classical case, it is reasonable to require the modified energy recover some aspects of (4) when $\alpha\to 1$, leading to the compatible energy LiaoTangZhou1 ; LiaoZhuWang ,

[TABLE]

where $\mathcal{I}^{\alpha}_{t}$ denotes the Riemann-Liouville fractional integral operator of order $\alpha$ defined by $(\mathcal{I}^{\alpha}_{t}u)(t)=(\kappa*u)(t)$ with the kernel function $\kappa(s)=s^{\alpha-1}/\Gamma(\alpha)$ and $*$ standing for the convolution. With this modified energy $E_{c}[u]$, one can get

[TABLE]

which degrades into the form E^{\prime}[u](t)+\big{\|}\frac{\delta E}{\delta u}\big{\|}^{2}\leq 0 if $\alpha\to 1$, c.f. (4) for the classical Allen-Cahn equation. Other types of modified energy can be found, for example, in hou2022second ; hou2021highly . In addition to energy decay, another important property of (3) is that its solution satisfies the maximum-principle TangYuZhou , i.e.,

[TABLE]

Note that the classicial model ($\alpha=1$) also holds such property as (8) ShenTangYang .

In the discrete case, it is hence necessary to develop numerical schemes that preserve the energy decay and maximum-principle for model (3). Historically, there has been enormous numerical exploration for tFAC equations or phase field models or more generally, the nonlinear subdiffustion models, where some energy stable schemes were developed (see, e.g., JiaZhangXuJiang ; ji2020adaptive ; huang2020optimal ; li2021conforming ; li2019nonconforming ; jin2018numerical ; bu2017finite ; li2018finite ; liu2009numerical ; zhang2020efficient ; mohebbi2013high ; li2018unconditionally ; gao2011compact ). However, the studies of preservation on discrete energy dissipation, or energy decay and maximum-principle are still limited. Liao et al.LiaoTangZhou presented a second-order nonuniform time-stepping scheme for (3) and demonstrated the scheme can preserve the discrete maximum-principle. In MaskariKaraa , the authors presented a finite element method for the time fractional Cahn-Hilliard equation using the fractional Euler formula in temporal direction and proved the scheme can preserve the energy dissipation. Very recently, the authors in LiaoTangZhou1 constructed a $(1+\alpha)$th order Crank-Nicolson scheme for the tFAC equation by applying an L1-type formula of the Riemann-Liouville derivative, which can preserve both the discrete version energy decay (7) and maximum-principle (8). Hou and Xu hou2022second explored the time fractional $L^{2}$ gradient flows where the fractional derivative was approximated by the well known L2-1σ method alikhanov2015new coupled with the so-called sum-of-exponentials (SOE) technique jiang2017fast and proved a modified type of energy dissipation on uniform meshes. They obtained the optimal second-order accuracy numerically by using graded meshes to compensate the low accuracy at initial time. However, to the best of our knowledge, high-order methods on uniform meshes that can preserve simultaneously the two intrinsic properties, i.e., the energy decay and maximum-principle, of tFAC equation are still missing in the literature. In this work, we shall fill this gap by a closer examination of a second-order accuracy formula called shifted fractional trapezoidal rule (SFTR) Yin01 ; YinLiuLiZhang222 , revealing some essential characteristics of this formula that are not shared with for most other second-order ones.

It has been a consensus that the solutions of time fractional models in general have intrinsic initial singularities, which may lead to order reduction for numerical methods of high order designed for smooth solutions Stynes ; zhao2014collocation . To resolve this issue among finite difference methods, two popular techniques are often considered, i.e., adding correction terms if the mesh is uniformly divided (see e.g. JinLazarovZhou ; Lubich1 ), or adopting nonuniform meshes that are much denser near the initial point (see e.g. LiaoTangZhou ; stynes2017error ; kopteva2019error , to name just a few). Some other techniques, such as constructing log type basis functions chen2020spectrally , or employing a nonlinear transformation for the time variable before discretization li2021novel , are also developed recently. Although our interest in this work mainly focuses on the study of the preservation of energy decay and maximum-principle, we demonstrate in numerical tests that our scheme, without any correction term, can result optimal accuracy at any positive time, even on uniform meshes. We argue that this is the first report of such phenomenon, i.e., optimal accuracy can be automatically recovered, for nonlinear time-fractional models. For a similar phenomenon examined for linear models, readers can refer YinLiuLiZhang222 .

We sum up the main contributions of this work:

• •

Prove theoretically some essential properties of the weights of SFTR-$\frac{1}{2}$ (i.e., SFTR with the shift parameter as $\frac{1}{2}$) and further demonstrate the discrete energy decay and maximum-principle preservation of the scheme.

• •

Demonstrate numerically the superiority of SFTR-$\frac{1}{2}$ that the intrinsic intial singularity can be resolved automatically.

The outline of the rest of the paper is as follows. In section 2, we prove the key properties of the weights of SFTR-$\frac{1}{2}$ which are crucial for the discrete maximum-principle preservation. In section 3, we first present the Crank-Nicolson finite difference scheme for the tFAC equation and then introduce a sequence closely related the initial weights of SFTR-$\frac{1}{2}$ with rigorous justification of its properties being essential to the discrete energy decay. At last, we prove the discrete maximum-princple. In section 4, we first examine our theoretical results, the discrete energy decay and maximum-principle preservation, by resorting to simulating the model (3) with a randomly distributed initial condition, and then verify the high accuracy of SFTR-$\frac{1}{2}$ formula, by comparisons with the classical fractional BDF2 and L2-1σ method, that our method can automatically resolve the initial singularity. Some concluding remarks are made in section 5.

2 SFTR and kernel properties

Divide the time interval by $0=t_{0}<t_{1}<t_{2}<\cdots<t_{N}=T$ with $t_{n}=n\tau$, $\tau=T/N$. Set $\phi^{n}=\phi(t_{n})$. At the node $t_{n-\frac{1}{2}}$, the SFTR-$\frac{1}{2}$ formula for the fractional derivative $\partial_{t}^{\alpha}\phi$ takes the form (see Yin01 ; YinLiuLiZhang222 )

[TABLE]

where the weights $\omega_{m}$ can either be obtained by the generating function

[TABLE]

or calculated directly by the recursive relation

[TABLE]

It has been shown in Yin01 the formula (9) approximates $\partial_{t}^{\alpha}\phi(t_{n-\frac{1}{2}})$ with second-order accuracy. In fact, it is a special case of general SFTR introduced in Yin01 ; YinLiuLiZhang222 . In the following, we give some peculiar characteristics of this formula.

Lemma 1

The following defined function $\Phi(x,y)$ is nonnegative for $(x,y)\in[0,1]\times[3,+\infty)$:**

[TABLE]

Proof

Replace $y$ with $3/\widetilde{y}$ for $\widetilde{y}\in(0,1]$ and set $\widetilde{\Phi}(x,\widetilde{y})=\Phi(x,y)$ to obtain

[TABLE]

for $(x,y)\in[0,1]\times(0,1]$. For fixed $\widetilde{y}$, it is easy to check $\frac{\partial\widetilde{\Phi}}{\partial x}\geq 0$, then $\widetilde{\Phi}(x,\widetilde{y})\geq\widetilde{\Phi}(0,\widetilde{y})=0$, which completes the proof of the lemma.

Lemma 2

For weights $\omega_{m}$ defined in (9), there hold (i) $\omega_{0}>0$, $\omega_{m}<0,m\geq 1$, (ii) $\omega_{m}>\omega_{m-1}$ $m\geq 2$, and (iii) $\displaystyle\sum_{m=0}^{n}\omega_{m}>0$ for any $n\geq 0$.

Proof

For (i), the correctness for $m=0,1$ can be checked directly. For $m=2$, (10) leads to

[TABLE]

For $m\geq 3$, we first show the following fact

[TABLE]

then by (10) we can get

[TABLE]

provided that $\omega_{m-2}<0$. In fact, the inequality (14) can be derived by induction arguments. Note first that $\sigma_{m-1}>\frac{\frac{1-\alpha}{2\alpha}(m-2)}{\frac{m-1}{\alpha}-\alpha}$ always holds for $m\geq 3$. For $m=3$, (14) is straightforward, as

[TABLE]

Then, for the $m+1$ case, we derive by the recursive formula (10) that

[TABLE]

where the assumption has been used within the first inequality. By Lemma 1 with $x,y$ replaced with $\alpha$ and $m$, respectively, we have $\frac{\omega_{m}}{\omega_{m-1}}\geq\sigma_{m}$ and complete the induction.

For (ii), $\omega_{2}>\omega_{1}$ can be directly derived by (13) and the fact $\omega_{m}<0$ for all $m\geq 1$. We assume $\frac{\omega_{m-1}}{\omega_{m-2}}<1$ for $m\geq 3$ and again by the recursive formula (10) we have

[TABLE]

For (iii), obviously we have $\displaystyle\sum_{m=0}^{\infty}\omega_{m}=0$, then by (i) there holds $\displaystyle\sum_{m=0}^{n}\omega_{m}>\sum_{m=0}^{\infty}\omega_{m}$ which completes the proof of the lemma.

3 Energy decay and maximum-principle

We first formulate the fully discrete scheme for model (3) by adopting finite difference methods in spatial direction and the SFTR-$\frac{1}{2}$ formula in temporal direction for the fractional derivative. The space $\Omega$ is divided by grid points $a=x_{0}<x_{1}<x_{2}<\dots<x_{M}=b$, $a=y_{0}<y_{1}<y_{2}<\dots<y_{M}=b$, where $x_{k}=y_{k}=a+kh$, $h=\frac{b-a}{M}$. Set $\psi_{j,k}=\psi(x_{j},y_{k})$. Take the grid function space $V_{h}=\{v=\{v_{j,k}\}|v_{j,k}\in\mathbb{R},1\leq j,k\leq M\}$ with norms:

[TABLE]

For any $u,v\in V_{h}$, define the discrete $L^{2}$ inner product as

[TABLE]

Introduce the following centered difference formula, for $1\leq j,k\leq M$:

[TABLE]

It is a well known result that for sufficiently smooth function $\psi$, there hold

[TABLE]

Define $\phi^{n-\frac{1}{2}}=\frac{\phi^{n}+\phi^{n-1}}{2}$ and let $u^{n}_{j,k}=u(x_{j},y_{k},t_{n})$. Take $U^{n}_{j,k}$ as the approximation to $u^{n}_{j,k}$, then for problem (3) at grid point $(x_{j},y_{k},t_{n-\frac{1}{2}})$, the second-order difference scheme for $1\leq j,k\leq M$, $1\leq n\leq N$ reads

[TABLE]

where $U^{0}_{j,k}=u_{0}(x_{j},y_{k})$ and the nonlinear term $f^{n-1,n}_{j,k}$ takes the form

[TABLE]

which immediately leads to the key property (see Appendix A in LiaoTangZhou1 ),

[TABLE]

To formulate the related matrix form, for any $v\in V_{h}$, define $\boldsymbol{v}\in\mathbb{R}^{M^{2}}$ by

[TABLE]

with $\|\boldsymbol{v}\|_{\infty}=\|v\|_{\infty}$, $\|\boldsymbol{v}\|=\|v\|$, and $(\boldsymbol{u},\boldsymbol{v})=(u,v)$ for any $u,v\in V_{h}$. The matrix form of (21) then reads

[TABLE]

where $\Delta_{h}=\frac{1}{h^{2}}[I\otimes D+D\otimes I]$ with $\otimes$ standing for the operation of Kronecker tensor product of matrices, and $D$ is a symmetric matrix:

[TABLE]

To show the unique solvability of the nonlinear scheme (24), we reformulate it as

[TABLE]

where the matrix \Upsilon_{n}={\rm Diag}\big{(}\tau^{-\alpha}\omega_{0}-\frac{1}{2}+\frac{1}{2}(\boldsymbol{U}^{n-1})^{\circ 2}\big{)}-\frac{\varepsilon^{2}}{2}\Delta_{h}, and $\mathcal{G}$ stands for the vector

[TABLE]

Here for the powers with a symbol $\circ$ we mean the calculation is performed elementwisely. It is a well known result that (25) admits a unique solution so long as the matrix $\Upsilon_{n}$ is positive definite LiaoTangZhou1 ; LiaoTangZhou ; HouTang01 , which is satisfied if we limit $\tau$ such that $\tau^{-\alpha}\omega_{0}-\frac{1}{2}>0$ as concluded in the following lemma.

Lemma 3

The nonlinear difference scheme (21) or (24) is uniquely solvable provided $0<\tau<2^{\frac{1}{\alpha}}\frac{2\alpha}{\alpha+1}$.

3.1 Energy decay

The discrete energy $E_{h}^{n}$ which is an analogue to $E[u](t_{n})$ in (1) can be defined by

[TABLE]

where $\nabla_{h}$ is an approximation to $\nabla$ satisfying $(-\Delta_{h}\boldsymbol{U},\boldsymbol{V})=(\nabla_{h}\boldsymbol{U},\nabla_{h}\boldsymbol{V})$.

The following sequence is important in deriving the energy decay property. Define $\{\theta_{m}\}_{m=0}^{\infty}$ satisfying

[TABLE]

Lemma 4

The sequence $\{\theta_{m}\}_{m=0}^{\infty}$ defined in (27) can be calculated by the recurrence formula:

[TABLE]

for $m\geq 2$, with the starting values \theta_{0}=\big{(}\frac{\alpha+1}{2\alpha}\big{)}^{\alpha}, $\theta_{1}=\frac{(\alpha-1)(2\alpha+1)}{\alpha+1}\theta_{0}$.

Proof

By the definition (27), $\theta(\xi)$ has the explicit form

[TABLE]

Take p_{1}(\xi)=\frac{1}{2}+\frac{1}{2\alpha}+\big{(}\frac{1}{2}-\frac{1}{2\alpha}\big{)}\xi and $p_{2}(\xi)=1-\xi$, then the formula (28) follows from Theorem A.1 of YinLiuLiZhang333 readily.

The weights $\theta_{m}$ actually hold a similar property as that stated in Lemma 2 for $\omega_{k}$, which is determined essentially by the nonnegativity of a special designed function

[TABLE]

Lemma 5

The function $\Psi(x,y)$ defined in (30) is nonnegative for $(x,y)\in[0,1]\times[3,\infty)$.

Proof

Let $y=3/\widetilde{y}$ for $\widetilde{y}\in(0,1]$, and $\widetilde{\Psi}(x,\widetilde{y})=\Psi(x,3/\widetilde{y})$, then one can get

[TABLE]

Some tedious but simple calculations show that $\frac{\partial\widetilde{\Psi}}{\partial\widetilde{y}}\geq 0$, and thus $\Psi(x,y)=\widetilde{\Psi}(x,\widetilde{y})\geq\widetilde{\Psi}(x,0)=0$.

Lemma 6

For weights $\theta_{m}$ defined in (27), there hold (i) $\theta_{0}>0$, $\theta_{m}<0,m\geq 1$, and (ii) $\displaystyle\sum_{m=0}^{n}\theta_{m}>0$ for any $n\geq 0$.

Proof

We first prove (i). The cases for $m=0,1$ can be verified directly. For $m=2$, by (28) we have $\theta_{2}=-\frac{2\alpha^{3}(1-\alpha)}{(1+\alpha)^{2}}\theta_{0}<0$. For $m=3$, similarly, one gets \theta_{3}=\frac{2\alpha}{3(1+\alpha)}\big{(}\alpha+\frac{3-\alpha}{2\alpha}\big{)}\theta_{2}<0. For $m\geq 4$, we first demonstrate that

[TABLE]

where

[TABLE]

with which the relation (28) yields

[TABLE]

provided that $\theta_{m-2}<0$. Actually, (31) can be proved by induction. For $m=4$, there holds

[TABLE]

We assume that (31) is valid, then from (28) one gets,

[TABLE]

which, by the assumption, leads to (for $m\geq 4$),

[TABLE]

where the last inequality is due to Lemma 5 with $x,y$ replaced by $\alpha,m$, respectively.

For (ii), it is obvious that $\displaystyle\sum_{m=0}^{\infty}\theta_{m}=0$, and by (i) we know that any finite sum with index from [math] to $n$ is positive. The proof is completed.

Remark 1

The generating function $\theta(\xi)$ in (29) in fact produces a difference formula for the Riemann-Liouville fractional differential operator ${}^{R}\mathcal{\partial}^{1-\alpha}_{t}=\partial_{t}\mathcal{I}^{\alpha}$ at time $t_{n}$ with second-order accuracy. To see that, on the one hand $\theta_{k}$ satisfies $\theta_{n}=O(n^{-(1-\alpha)-1})$ since \big{[}\frac{1}{2}+\frac{1}{2\alpha}+\big{(}\frac{1}{2}-\frac{1}{2\alpha}\big{)}\xi\big{]}^{\alpha} is analytic on the closed unite disc $|\xi|\leq 1$ for any fixed $\alpha\in(0,1)$, and on the other hand,

[TABLE]

Then, the convolution quadrature theory Lubich1 guarantees $\tau^{\alpha-1}\displaystyle\sum_{m=0}^{n}\theta_{m}\phi^{n-m}$ approximates ${}^{R}\mathcal{\partial}^{1-\alpha}_{t}\phi(t_{n})$ with second-order accuracy under mild conditions.

Remark 2

We emphasize that the properties revealed in Lemma 2 and Lemma 6 are essential to the following analysis of discrete energy decay and maximum-principle preservation. The difficulties in the proof for these two lemmas lie in the subtle design of (14) and (31), i.e., the choice of $\sigma_{m}$ and $\delta_{m}$, both of which lead to nonnegative functions described in Lemma 1 and Lemma 5, respectively.

To discretize the integral part within (6) based on the weights of SFTR-$\frac{1}{2}$, we introduce another sequence $\{\vartheta_{m}\}_{m=0}^{\infty}$ defined as

[TABLE]

Then, by Cauchy product of series $\displaystyle\sum_{m=0}^{\infty}\theta_{m}\xi^{m}$ and $\displaystyle\sum_{m=0}^{\infty}\xi^{m}$, i.e.,

[TABLE]

one gets $\vartheta_{m}=\displaystyle\sum_{s=0}^{m}\theta_{s}$, followed by the fact

[TABLE]

thanks to Lemma 6. By similar arguments stated in Remark 1, we have

[TABLE]

Then, the shifted convolution quadrature theory YinLiuLiZhang444 indicates that, under mild conditions, $\tau^{\alpha}\displaystyle\sum_{m=0}^{n}\vartheta_{m}\phi^{n-m}$ approximates $\mathcal{I}_{t}^{\alpha}\phi(t_{n+\frac{1}{2}})$ with second-order accuracy. Now, with the help of the sequence $\{\vartheta_{m}\}_{m=0}^{\infty}$, we define the discrete counterpart of $E_{c}[u]$ in (6) by

[TABLE]

where $\boldsymbol{V}^{n-\frac{1}{2}}:=\varepsilon^{2}\Delta_{h}\boldsymbol{U}^{n-\frac{1}{2}}-\boldsymbol{f}^{n-1,n}$ is an approximation to $\frac{\delta E}{\delta u}$ at time $t_{n-\frac{1}{2}}$.

Based on Lemma 6 and the monotonicity of $\vartheta_{k}$ in (37), the following energy decay analysis can be carried out in a standard manner LiaoTangZhou1 . We present the details of the analysis within our setting for integrity.

Theorem 3.1

The discrete energy $\mathcal{E}^{n}_{h}$ defined in (39) satisfies the decay property $\mathcal{E}^{n}_{h}\leq\mathcal{E}^{n-1}_{h}$, $n\geq 1$.

Proof

Multiply both hand sides of (24) by $\tau^{\alpha-1}\theta_{n-s}$ after replacing $n$ with $s$ in (24), and sum the index $s$ from $1$ to $n$ to obtain

[TABLE]

Rearranging the terms in the left hand side of the above equation, one gets

[TABLE]

where $c_{m}=\sum_{s=0}^{m}\omega_{s}\theta_{m-s}$ is the $m$-th coefficient of the series of the product $\omega(\xi)\theta(\xi)$, which by (27) yields $c_{0}=1,c_{1}=-1$ and $c_{k}=0$ for $k>1$. We then have

[TABLE]

Taking the $L^{2}$ inner product of (40) with respect to $\boldsymbol{V}^{n-\frac{1}{2}}$, we get

[TABLE]

Assume $\boldsymbol{V}^{s}={\rm vec}(v^{s})$ for some $v^{s}=\{v_{j,k}^{s}\}\in V_{h}$. Since $\theta_{s}=\vartheta_{s}-\vartheta_{s-1}<0,s\geq 1$ and $\theta_{0}=\vartheta_{0}>0$ by Lemma 6, there holds

[TABLE]

Then, for the left hand side of (41) we have

[TABLE]

Moreover, since \big{(}\boldsymbol{U}^{n}-\boldsymbol{U}^{n-1},\Delta_{h}\boldsymbol{U}^{n-\frac{1}{2}}\big{)}=-\big{(}\|\nabla_{h}\boldsymbol{U}^{n}\|^{2}-\|\nabla_{h}\boldsymbol{U}^{n-1}\|^{2}\big{)}, combined with (23), the right hand side of (41) can be estimated by

[TABLE]

which combined with (43) leads to

[TABLE]

This completes the proof of the theorem.

Remark 3

It is notable that when $\alpha\to 1$, $\theta(\xi)=1$ by (29) and thus $\vartheta_{m}=1$ for $m\geq 0$. In this case, $\mathcal{E}^{n}_{h}$ defined in (39) becomes to

[TABLE]

and (45) degrades into

[TABLE]

which is compatible with the classical case (4).

3.2 Maximum-principle

The properties of $\omega_{m}$ proved in Lemma 2 and the form of nonlinear term approximation (22) followed by some standard analysis essentially lead to the discrete maximum-principle (see (LiaoTangZhou1, , Theorem 3.2)). We shall present in the next theorem the outline of the proof to find out the restriction on time-step size $\tau$ for our scheme.

Theorem 3.2

The numerical scheme (21) or (24) preserves the discrete maximum-principle, i.e., for $\tau$ satisfying 0<\tau<\min\big{\{}2^{\frac{1}{\alpha}}\frac{2\alpha}{\alpha+1},(\frac{\alpha h^{2}}{2\varepsilon^{2}})^{\frac{1}{\alpha}}(\frac{2\alpha}{\alpha+1})^{\frac{\alpha+1}{\alpha}}\big{\}}, there holds $\|\boldsymbol{U}^{n}\|_{\infty}\leq 1$ for all $n\geq 1$ provided $u_{0}(\boldsymbol{x})\in C(\overline{\Omega})$ and $\max_{\boldsymbol{x}\in\overline{\Omega}}|u_{0}(\boldsymbol{x})|\leq 1$.

Proof

With $\Upsilon_{n}$ defined in (25), we rewrite (24) into

[TABLE]

where

[TABLE]

Assume $\|\boldsymbol{U}^{m}\|_{\infty}\leq 1$ for $m=0,1,\cdots,n-1$, and let G_{h}=(g_{jk})=\big{(}\frac{\varepsilon^{2}}{2}\Delta_{h}-\tau^{-\alpha}\omega_{1}I\big{)}. Then, there hold $g_{jk}>0$ for $j\neq k$ and $g_{jj}=-\frac{2\varepsilon^{2}}{h^{2}}-\tau^{-\alpha}\omega_{1}$, and $\sum_{k}g_{jk}=-\tau^{-\alpha}\omega_{1}$ for any $j$. If \tau\leq\big{(}\frac{\alpha h^{2}}{2\varepsilon^{2}}\big{)}^{\frac{1}{\alpha}}\big{(}\frac{2\alpha}{\alpha+1}\big{)}^{\frac{\alpha+1}{\alpha}}, i.e., $g_{jj}\geq 0$, then

[TABLE]

Thanks to Lemma 2, the $\mathcal{G}_{2}$ term can be bounded by $\|\mathcal{G}_{2}\|_{\infty}\leq\tau^{-\alpha}(\omega_{0}+\omega_{1})$. Since $|3z-z^{3}|\leq 2$ for any $z\in[-1,1]$, then one gets \|\frac{1}{6}\big{[}3\boldsymbol{U}^{n-1}-(\boldsymbol{U}^{n-1})^{\circ 3}\big{]}\|_{\infty}\leq\frac{1}{3}. Moreover, the unique solvability condition on $\tau$ in Lemma 3, i.e., $\tau<2^{\frac{1}{\alpha}}\frac{2\alpha}{\alpha+1}$ can lead to the following estimates (see (LiaoTangZhou1, , Lemma 3.5)),

[TABLE]

The above term by term estimates of (LABEL:ene.20) combined with a subtle contradiction argument finally give the desired conclusion. More details can be found in LiaoTangZhou1 .

Remark 4

In this paper, our focus is on the structure-preserving properties of the scheme. As for its convergence, here we give some remarks. If the exact solution is sufficiently smooth, it is easy to prove the second-order convergence of the scheme by using the standard energy method. Here we omit the proof because it is more reasonable to think that the solution has a weak singularity at $t=0$. In the next section, we will give a numerical investigation, and the results show that the method has second order convergence at fixed positive time when solving the tFAC equation (3). Despite the rigorous arguments of sharp error estimate is still open to us, we would like to give some insights into why SFTR-$\frac{1}{2}$ automatically preserve the optimal accuracy at fixed time in numerical tests. Actually, for the linear counterpart of (3), i.e., the linear subdiffusion problem with $f(u)=0$, rigorous arguments have been carried out in YinLiuLiZhang222 for SFTR-$\theta$ with $\theta\in(0,\frac{1}{2})$, where the correction term

[TABLE]

is added to the initial time step to compensate for the weak regularity of the solution and to restore the second-order accuracy. $\tilde{\Delta}_{h}$ denotes some discrete approximation to the Laplacian. The correction term can be arbitrarily small so long as the given $\theta$ approaches $\frac{1}{2}$ adequately while the second-order accuracy is still guaranteed, and it is then quite natural to ask what may happen if $\theta$ takes $\frac{1}{2}$, since in this case the correction term vanishes formally. However, the rigorous theoretical analysis for the limit case $\theta=\frac{1}{2}$ is still an open problem for that the arguments developed for $\theta\in(0,\frac{1}{2})$ can not be simply carried out by setting $\theta=\frac{1}{2}$ directly. As shown in YinLiuLiZhang222 , the arguments for $\theta\in(0,\frac{1}{2})$ rely essentially on the solution representation theory where a contour integral is involved. To be specific, let $U_{h}^{n}$ be the approximation to $u(t_{n})$ within some space discretization framework such as the finite element method, then $U_{h}^{n}$ can be represented by

[TABLE]

where $K(z)=-z^{-1}(z^{\alpha}-\tilde{\Delta}_{h})^{-1}$, \beta_{\tau}(\zeta)=\frac{1-\zeta}{\tau(1+\mu_{1}\zeta)}\big{(}\frac{\mu_{0}}{1-\theta+\theta\zeta}\big{)}^{\frac{1}{\alpha}} and $\mu(\zeta)=\frac{\mu_{0}^{\frac{1}{\alpha}}(3/2-\theta)\zeta(1-\mu_{2}\zeta)}{(1+\mu_{1}\zeta)(1-\theta+\theta\zeta)^{1+\frac{1}{\alpha}}}$ with the coefficients $\mu_{0},\mu_{1}$ and $\mu_{2}$ defined by \mu_{0}=\big{(}\frac{2\alpha}{\alpha+2\theta}\big{)}^{\alpha}, $\mu_{1}=\frac{\alpha-2\theta}{\alpha+2\theta}$ and $\mu_{2}=\frac{1-2\theta}{3-2\theta}$.

To obtain (52), it is required (Theorem 4.3 in YinLiuLiZhang222 ) that both $\beta_{\tau}(\zeta)$ and $\mu(\zeta)$ are analytic at $\zeta\in\{\zeta\in\mathbb{C}:|\zeta|\leq 1,\zeta\neq 1\}$. Unfortunately, when $\theta=\frac{1}{2}$, both $\beta_{\tau}(\zeta)$ and $\mu(\zeta)$ are singular at point $\zeta=-1$, which implies the underlying arguments for $\theta=\frac{1}{2}$ is non-standard and that replacing $\theta$ with $\frac{1}{2}$ directly is not acceptable. For the nonlinear subdiffusion problem such as the tFAC equation, more technique should be involved. Despite the difficulties and unknowns in the error estimates, the following numerical simulations are still valuable which, for the one hand suggests SFTR-$\frac{1}{2}$ is superior to the averaged fractional BDF2 on uniform meshes and is comparable to L2-1σ method on graded meshes, and for the other hand convince us by the idea that difference methods proposed at shifted non-integer grids can resolve the initial singularity of fractional problems to some extent.

4 Numerical experiments

We implement several examples in this section to verify the energy decay and maximum-principle preservation, and demonstrate the high efficiency of SFTR-$\frac{1}{2}$ under low solution regularity. Since the fully discrete scheme (21) or (24) is nonlinear, we adopt a simple iteration method with the stopping criterion that $\|U^{(k+1)}-U^{(k)}\|_{\infty}\leq 10^{-6}$.

4.1 Verification of energy decay and maximum-principle

Example 1. Let $\Omega=(0,1)\times(0,1)$ and $\varepsilon=0.01$ and take the initial condition $u_{0}(\boldsymbol{x})=0.5\times{\rm rand}(\boldsymbol{x})-0.25$ where the function ${\rm rand}(\boldsymbol{x})$ returns uniformly distributed random numbers in $(0,1)$ and thus $|u_{0}(\boldsymbol{x})|\leq 1$. The time space mesh is chosen with $\tau=0.05$ and $h=\frac{1}{200}$. We simulate the coarsening dynamics under the above setting and illustrate some snapshots in Fig.1 up to $T=20$ for different fractional differential order $\alpha=0.3,0.6$ and $0.9$, respectively. The maximum of $|U^{n}|$ and discrete compatible energy $\mathcal{E}_{h}^{n}$ for each time level are depicted in Fig.2 (a) and (b), for $\alpha=0.3,0.6,0.9$ respectively. Obviously, under this random initial condition, the discrete maximum-principle and energy decay property are preserved.

4.2 High accuracy under low regularity

We demonstrate in the following two examples the serious impact of low regularity of solution on temporal convergence accuracy by first simulating a model with smooth enough solution (in which case the source term is nonzero), and then exploring the original model and comparing the results between our method and the well known fractional BDF2 (F-BDF2) Lubich1 . For both examples, let $\Omega=(0,1)\times(0,1)$, $T=1$ and $\varepsilon=0.005$.

Example 2. To approximate $\partial_{t}^{\alpha}u$ at time $t_{n-\frac{1}{2}}$ by the F-BDF2, we take the following process,

[TABLE]

where $\partial_{t}^{\alpha}u(t_{n})\approx\tau^{-\alpha}\sum_{m=0}^{n}\widetilde{\omega}_{n-m}(u(t_{m})-u(t_{0}))$ with $\widetilde{\omega}_{m}$ generated by the function

[TABLE]

Assume the exact solution is

[TABLE]

The source term can be derived easily. Note that such solution is smooth in time from which nice convergence results can be expected. The temporal rates at time level $t_{n}$ is obtained by the formula:

[TABLE]

In Table 1, we report for different $\alpha=0.3,0.6,0.9$ the $L^{2}$ error at time $t=\frac{T}{2}$ and convergence rates for both SFTR-$\frac{1}{2}$ and F-BDF2, by taking different time-step size $\tau=\frac{1}{20},\frac{1}{40},\frac{1}{80}$ and $\frac{1}{160}$, respectively. For such smooth solution, the optimal second-order convergence accuracy is obtained which is in line with our expectation.

Example 3. This example is devoted to the model (3), i.e., with zero source term. Take the initial condition as $u_{0}(\boldsymbol{x})=u_{0}(x_{1},x_{2})=\sin(2\pi x_{1})\sin(2\pi x_{2})$, with no explicit exact solution known in advance. We therefore adopt the numerical solution (denoted by $U^{n}_{ref}$) obtained at the fine mesh with $\overline{h}=\frac{1}{200},\tau=\frac{1}{400}$ as the reference solution, and the following formula to calculate the temporal rates at time $t_{n}$:

[TABLE]

Similar as in Example 2, the errors and convergence rates are reported for $\alpha=0.3,0.6$ and $0.9$, respectively. However, unlike the results in Example 2, the SFTR-$\frac{1}{2}$ can still yield almost optimal convergence accuracy even through the solution is weak regular, which is much better than the classical F-BDF2 suffering a severe reduction in precision.

4.3 Comparison with the L2-1σ formula

Example 4. In this example, we compare our scheme with that in LiaoTangZhou formulated by the well-known L2-1σ formula on nonuniform meshes. For better demonstration, we calculate the $L^{2}$ norm errors at the finial time $t=T$, under the condition that the number of grid points among $[0,T]$, which is denoted by $N$, is exactly the same for both schemes.

Let $\Omega=(0,1)\times(0,1)$, $T=1$ and $\varepsilon=0.01$. Choose the same initial condition as that in Example 3. By graded meshes we mean that for some given $\gamma\geq 1$, the interval $[0,T]$ is divided by points $t_{n}=T(n/N)^{\gamma}$. Similar as Example 3, the reference solutions $U_{ref}$, which depend generally on $\alpha,h,N$ and on $\gamma$ in addition if graded meshes are considered, are obtained under the fine mesh with $h=\overline{h}=\frac{1}{200},N=400$. Fig.3 illustrates the errors in $L^{2}$ norm of the scheme LiaoTangZhou for different $\gamma$ with $\gamma\in[1,\frac{2}{\alpha}]$. Note that if $\gamma=1$, the mesh is uniform and if $\gamma=\frac{2}{\alpha}$, the uniform second-order convergence rate is arrived at according to LiaoTangZhou . One easily observes from Fig.3 that to obtain the minimum error at final time $t=T$, the mesh parameter $\gamma$ should be between $1$ and $\frac{2}{\alpha}$, e.g., from Fig.3 (a) where $\alpha=0.2$, $\gamma=2.6$ results in the minimum error. We therefore in Table 3 collect some key results and denote by “L2-1σ (uniform)” the errors on uniform mesh (i.e., $\gamma=1$) and by “L2-1σ (min)” the minimum errors among $\gamma\in[1,\frac{2}{\alpha}]$. Obviously, SFTR-$\frac{1}{2}$ is much more accurate than L2-1σ on uniform meshes and both of these two methods are evenly matched (so long as $\alpha$ is not too small) if graded meshes are equipped to L2-1σ. We also remark that SFTR-$\frac{1}{2}$ performs even better for large $\alpha$ (e.g., $\alpha\geq 0.4$ in our example), though L2-1σ is more accurate when $\alpha$ is small in which case more work must be involved to avoid the round-off errors polluting the coefficients of the formula (see (3.2) and (3.3) in LiaoTangZhou ).

5 Conclusion

In this work, we develop a finite difference scheme for time-fractional Allen-Cahn equation by SFTR-$\frac{1}{2}$ for time discretization and the centered difference formula for space variables. We first prove some key properties of the weights of SFTR-$\frac{1}{2}$ and further show that the scheme proposed is uniquely solvable. Moreover, the discrete energy decay property of the proposed scheme is exhibited by using some properties of two sets of coefficients $\theta_{k}$ and $\vartheta_{k}$, which are related to $\omega_{k}$. The discrete maximum-principle preserved of the scheme is also proved rigorously, by resorting to the negativity of the weights $\omega_{k}$ for $k\geq 1$. Numerical results confirm the correctness of our theoretical results concerning discrete energy decay law and maximum-principle, and reveal that SFTR-$\frac{1}{2}$ can resolve automatically the intrinsic singularity of tFAC equation, indicating SFTR-$\frac{1}{2}$ is superior to other high-order approximation formulas on uniform meshes (e.g., the fractional BDF2 or L2-1σ). Comparisons with L2-1σ method under graded meshes also confirm the high accuracy of our method.

There are some aspects deserving further study, e.g., proposing the theoretical analysis for the optimal convergence results of SFTR-$\frac{1}{2}$ when applied to nonlinear fractional equations, which will be our future work.

Acknowledgments

The authors are grateful to two anonymous referees and editors for their valuable suggestions which improve the presentation of this work greatly. This work is supported by the Natural Science Foundation of Inner Mongolia Autonomous Region of China (No. 2021BS01003 to G.Z.), NSF of China (Nos. 12171177 and 12011530058 to C.H., 12201322 to B.Y.), RFBR (No. 20-51-53007 to A. A.) and North-Caucasus Center for Mathematical Research under agreement No. 075-02-2021-1749 with the Ministry of Science and Higher Education of the Russian Federation.

Conflict of interest

The authors declare that they have no conflict of interest.