Simple maximum-principle preserving time-stepping methods for time-fractional Allen-Cahn equation
Bingquan Ji Hong-lin Liao Luming Zhang





Simple maximum-principle preserving time-stepping methods for time-fractional Allen-Cahn equation
    Bingquan JiDepartment of Mathematics, Nanjing University of Aeronautics and Astronautics,
211101, P. R. China. Bingquan Ji (jibingquanm@163.com).
  Hong-lin LiaoCorresponding author. Department of Mathematics,
Nanjing University of Aeronautics and Astronautics,
Nanjing 211106, P. R. China. Hong-lin Liao (liaohl@nuaa.edu.cn, liaohl@csrc.ac.cn)
is supported by a grant 1008-56SYAH18037
from NUAA Scientific Research Starting Fund of Introduced Talent.
  Luming ZhangDepartment of Mathematics, Nanjing University of Aeronautics and Astronautics,
211101, P. R. China. Luming Zhang (zhanglm@nuaa.edu.cn)
is supported by the research grants No. 11571181 from National Natural Science Foundation of China.
    
================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




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önwall inequality and global consistency analysis.
 Numerical experiments are presented to show the effectiveness of our methods and to confirm our analysis.


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




§ 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 <cit.>.
Relevant numerical methods and simulations  are also increasing substantially <cit.>. 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 <cit.>, 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
<cit.>
to accurately describe anomalous diffusion problems.
Li et al. <cit.> 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. <cit.> 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 <cit.> 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. <cit.> 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 <cit.>

    ∂_t^αu
    =ε^2Δu-f(u),  𝐱∈Ω,  0<t≤T,
       u(𝐱,0)=u_0(𝐱),  𝐱∈Ω̅,

where 𝐱=(x,y)^T and
Ω=(a,b)×(c,d) with its closure Ω̅.
The notation ∂_t^α:=_0^CD_t^α in (<ref>) denotes
the fractional Caputo derivative of order α with respect to t,

    (∂_t^αv)(t)
    :=(ℐ_t^1-αv')(t)=∫_0^tω_1-α(t-s)v'(s)s,   0<α<1,

involving the fractional  Riemann-Liouville integral ℐ_t^μ of order μ>0, that is,

    (ℐ_t^μv)(t)
    :=∫_0^tω_μ(t-s)v(s)s,  where  ω_μ(t):=t^μ-1/Γ(μ).


The nonlinear bulk force f(u)=u^3-u, and the small constant ε>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 <cit.>.
In comparison to the classical physical model, the energy dissipation law of the
time-fractional Allen-Cahn equation (<ref>) is

    E(t)≤E(0),

where

    E(t):=∫_Ω[ε^2/2|∇u|^2
    +F(u)]𝐱,  
    F(u)=1/4(1-u^2)^2.

Also, it possesses a maximum principle, namely,

    |u(𝐱,t)|≤1 for t>0  if  |u(𝐱,0)|≤1.

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 (<ref>).

We consider the nonuniform time levels 0=t_0<t_1<⋯<t_k-1<t_k<⋯<t_N=T with the time-step sizes τ_k:=t_k-t_k-1 for 1≤k≤N and the maximum time-step size τ:=max_1≤k≤Nτ_k. Also, let the local time-step ratio ρ_k:=τ_k/τ_k+1 and the maximum step ratio ρ:=max_k≥ 1ρ_k. Given a grid function {v^k}, put ▿_τv^k:=v^k-v^k-1, ∂_τv^k-1/2:=▿_τv^k/τ_k and
v^k-1/2:=(v^k+v^k-1)/2 for k≥1. Always,
let (Π_1,kv)(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

    Π_1v:=Π_1,kv  so that  (Π_1v)'(t)=∂_τv^k-1/2  for t_k-1<t≤ t_k and k≥1.









As an essential mathematical feature of linear and nonlinear subdiffusion problems
including the time-fractional Allen-Cahn problem (<ref>)-(<ref>),
the solution always lacks the smoothness near the initial time
although it would be smooth away from  t=0, see <cit.>.
Actually, assuming the nonlinear function f is Lipschitz continuous
and the initial data u^0∈ H^2(Ω)∩ H_0^1(Ω),
Jin et al. <cit.> proved the subdiffusion problem has an unique solution u for which
u∈ C⟨[0,T];H^2(Ω)∩ H_0^1(Ω)|, ∂_t^αu∈ C⟨[0,T];L^2(Ω)| and ∂_t u ∈ L^2(Ω)  with
∂_t u(t)_L^2(Ω)≤ C_ut^α-1 for 0<t≤ T.




The L1 scheme with a lagging linearized technique for handling the nonlinearity f(u)
has been analyzed, and <cit.> showed that the discrete solution is O(τ^α)
convergent in L^∞(L^2(Ω)). 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 <cit.>. More directly, we consider the L1 scheme for the time-fractional problem (<ref>)-(<ref>)
describing the coalescence of two kissing bubbles, see more details in Example <ref>.
Fig. <ref> plots the discrete time derivative ∂_τv^k-1/2
  near t=0 on the graded mesh t_k=(k/N)^γ.
They suggest that

    log|u_t(𝐱,t)|≈(α-1)logt+C(𝐱)  as t→0.

It says that the solution possesses weak singularity like u_t=O(t^α-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 (<ref>) 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 (∂_f^αu)^n described in (<ref>).
Theorems <ref> and <ref> show that both the backward Euler method (<ref>)-(<ref>)
and stabilized semi-implicit method (<ref>)-(<ref>)
preserve the maximum principle (<ref>) in the discrete level
such that they are unconditionally stable in the maximum norm.
By using the recently proposed discrete fractional Grönwall inequality <cit.>
and the global consistency analysis <cit.> of L1 formula,
 we prove that, see Theorems <ref> and <ref>, the fully implicit method
 (<ref>)-(<ref>) is convergent
with an optimal order of O(τ^2-α) and the stabilized scheme (<ref>)-(<ref>) is convergent
with an optimal order of O(τ) in time on the graded meshes with a grading parameter γ≥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 <ref>).


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.




§ FAST L1 TIME-STEPPING METHODS



The well-known L1 formula of Caputo derivative (<ref>) is given by

    (∂_τ^αv)^n
    :=∫_t_0^t_nω_1-α(t_n-s)(Π_1v)'(s)s
    =∑_k=1^na_n-k^(n)▿_τv^k,

where the corresponding discrete convolution kernels a_n-k^(n) are defined by

    a_n-k^(n)
    :=1/τ_k∫_t_k-1^t_kω_1-α(t_n-s)s  for 1≤k≤n.


Obviously, the discrete convolutional kernels a_n-k^(n) are positive and decreasing, see also <cit.>,

    a_n-k^(n)>0   and   a_n-k-1^(n)>a_n-k^(n)  for 1≤k≤n-1.

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





 §.§ Fast L1 formula

It is well known that the standard L1 formula (<ref>) 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^-α efficiently on the interval [Δt, T],
see <cit.>.

For the given α∈(0, 1), an absolute tolerance error ϵ≪1, a cut-off time Δt>0 and a finial time T, there exists a positive integer N_q, positive quadrature nodes θ^ℓ and corresponding positive weights ϖ^ℓ (1≤ℓ≤N_q) such that

    |
    ω_1-α(t)
    -∑_ℓ=1^N_qϖ^ℓe^-θ^ℓt|≤ϵ,
      ∀ t∈[Δt, T].




To be more precise, the Caputo derivative (<ref>) 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,

    (∂_t^αv)(t_n)
       ≈∫_0^t_n-1v^'(s)
    ∑_ℓ=1^N_qϖ^ℓe^-θ^ℓ(t_n-s)s
    +∫_t_n-1^t_nω_1-α(t_n-s)
    ▿_τv^n/τ_ns
       =∑_ℓ=1^N_qϖ^ℓe^-θ^ℓτ_nℋ^ℓ(t_n-1)
    +a_0^(n)▿_τv^n,   n≥1,

where ℋ^ℓ(t_0):=0 and ℋ^ℓ(t_k):=∫_0^t_ke^-θ^ℓ(t_k-s)v^'(s)s.
By utilizing the linear interpolation and a recursive formula, we can approximate ℋ^ℓ(t_k) by

    ℋ^ℓ(t_k)
       ≈∫_0^t_k-1e^-θ^ℓ(t_k-s)v^'(s)s
    +∫_t_k-1^t_ke^-θ^ℓ(t_k-s)▿_τv^k/τ_ks
       =e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +b^(k,l)▿_τv^k,

where the positive coefficients

    b^(k,l):=1/τ_k∫_t_k-1^t_ke^-θ^ℓ(t_k-s)s,   k≥1.

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

    (∂_f^αv)^n
    =
    a_0^(n)▿_τv^n
    +∑_ℓ=1^N_qϖ^ℓe^-θ^ℓτ_nℋ^ℓ(t_n-1),
       n≥1,

in which ℋ^ℓ(t_k) is computed by using the recursive relationship

    ℋ^ℓ(t_k)
    =e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +b^(k,l)▿_τv^k,   k≥1,   1≤ℓ≤N_q.


For the convenience of numerical analysis, we now eliminate the historic term ℋ^ℓ(t_k)
from the fast L1
formula (<ref>). From the recursive equation (<ref>),
a direct calculation yields

    ℋ^ℓ(t_k)
    =∑_j=1^ke^-θ^ℓ(t_k-t_j)b^(j,l)▿_τv^j,   k≥1,   1≤l≤N_q.

By substituting (<ref>) into (<ref>), we get an alternative definition

    (∂_f^αv)^n
    =a_0^(n)▿_τv^n
    +∑_k=1^n-1▿_τv^k/τ_k∫_t_k-1^t_k∑_ℓ=1^N_qϖ^ℓe^-θ^ℓ(t_n-s)s
    =∑_k=1^nA_n-k^(n)▿_τv^k,   n≥1,

where the corresponding discrete convolution coefficient A_n-k^(n) is defined by

    A_0^(n):=a_0^(n),  
    A_n-k^(n):=1/τ_k∫_t_k-1^t_k∑_ℓ=1^N_qϖ^ℓe^-θ^ℓ(t_n-s)s,  1≤k≤n-1, n≥1.

For the discrete kernels A_j^(n),
we have the following result <cit.>.

If the tolerance error ϵ of SOE satisfies ϵ≤min{1/3ω_1-α(T), αω_2-α(1)}, then the discrete convolutional kernel A_n-k^(n) of (<ref>) satisfies

  
  (i) A_k-1^(n)>A_k^(n)>0 for 1≤k≤n-1;
 
  (ii) A_0^(n)=a_0^(n) and  A_n-k^(n)≥2/3a_n-k^(n) for 1≤k≤n-1.





 §.§ 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≤ i≤ M_1, 0≤ j≤ M_2.
Also, denote Ω̅_h:={𝐱_h=(x_i,y_j) | 0≤ i≤ M_1, 0≤ j≤ M_2}
and put Ω_h:=Ω̅_h∩Ω.
For any grid function {v_h | 𝐱_h∈Ω̅_h}, denote a grid function space

    𝕍_h:={v | v=(v_j)^T  for  0≤j≤M_2-1,
      with  
    v_j=(v_i,j)^T  for  0≤i≤M_1-1},

where v^T is the transpose of the vector v.
The maximum norm v_∞:=max_𝐱_h∈Ω_h|v_h|.






Let D_h be the discrete matrix of Laplace operator Δ subject to periodic boundary conditions.
With the Kronecker tensor product ⊗, the matrix D_h=I_1⊗D_1+D_2⊗I_2,



in which I_1 and I_2 are the identity matrices of order M_2×M_2 and M_1×M_1, respectively,
 and the matrices D_1 and D_2 are of forms

    D_1
    =1/h_1^2(
    [ -2  1  0  ⋯  1;  1 -2  1  ⋯  0;  ⋮  ⋱  ⋱  ⋱  ⋮;  0  ⋯  1 -2  1;  1  ⋯  0  1 -2;    ])_M_1×M_1, 
    D_2
    =1/h_2^2(
    [ -2  1  0  ⋯  1;  1 -2  1  ⋯  0;  ⋮  ⋱  ⋱  ⋱  ⋮;  0  ⋯  1 -2  1;  1  ⋯  0  1 -2;    ])_M_2×M_2.

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.


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∈𝕍_h, v^TD_hv≤0, i.e., the matrix D_h is negative semi-definite.
  
  (c) The elements of D_h=(d_ij) fulfill d_ii=-max_i∑_j≠i|d_ij| for each i.
  





Now we have the backward Euler-type scheme on irregular meshes
for (<ref>)-(<ref>),


    ∂_f^αu^n   =ε^2D_hu^n
    -f(u^n),  n≥1,
    
    u_h^0   =u_0(𝐱_h),  𝐱_h∈Ω̅_h,

where f(u^n):=(u^n)^3-u^n with the vector (u^n)^3=⟨(u_1^n)^3,(u_2^n)^3,⋯,(u_M_2-1^n)^3|^T and

    (u_j^n)^3
    =⟨(u_1,j^n)^3,(u_2,j^n)^3,⋯,(u_M_1-1,j^n)^3|^T  for j=0,1,⋯,M_2-1.



Now we prove that the fully discrete scheme (<ref>)-(<ref>) preserves the maximum principle numerically.
Always, we need the following result <cit.>.


Let B be a real M×M matrix and A=aI-B with a>0.
If the elements of B=(b_ij) fulfill b_ii=-max_i∑_j≠i|b_ij|,
then for any c>0 and V∈ℝ^M we have

    AV_∞≥aV_∞,  AV+c(V)^3_∞≥aV_∞+cV_∞^3.





If u^0_∞≤1 and the maximum time-step size τ≤1/√(Γ(2-α)),
then the solution of backward Euler scheme (<ref>)-(<ref>) satisfies u^k_∞≤1
for 0≤k≤N.



So it preserves the maximum principle (<ref>) numerically
and is unconditionally stable.



We use the complete mathematical induction.
Obviously, the claimed inequality holds for k=0.
For 1≤ n≤ N, assume that

    u^k_∞≤1  for 0≤k≤n-1.

It remains to verify that u^n_∞≤1. From the definition (<ref>), one has

    (∂_f^αu)^n
    =A_0^(n)u^n-L^n-1  where   L^n-1:=∑_k=1^n-1(A_n-k-1^(n)-A_n-k^(n))u^k+A_n-1^(n)u^0.

Thanks to the decreasing property in Lemma <ref> (i), the induction hypothesis (<ref>)
and the triangle inequality yield

    L^n-1_∞≤∑_k=1^n-1(A_n-k-1^(n)-A_n-k^(n))u^k_∞+A_n-1^(n)u^0_∞≤ A_0^(n) .

Then, from the numerical scheme (<ref>), it is easy to obtain

    (A_0^(n)-1)u^n+(u^n)^3
    -ε^2D_hu^n_∞=L^n-1_∞≤ A_0^(n).

For the left hand side of (<ref>), we apply Lemma <ref> (c) and Lemma <ref>
to find that

    (A_0^(n)-1)u^n+(u^n)^3
    -ε^2D_hu^n_∞≥
    (A_0^(n)-1)u^n_∞
    +u^n_∞^3.

Then it follows from (<ref>) that (A_0^(n)-1)u^n_∞
+u^n_∞^3≤ A_0^(n).




If A_0^(n)≥1 or the maximum step size τ≤1/√(Γ(2-α)),
the above inequality implies u^n_∞≤1 immediately.
Otherwise, we have (A_0^(n)-1)u^n_∞
+u^n_∞^3>A_0^(n),
because the function

    g(z):=(A_0^(n)-1)z
    +z^3-A_0^(n)  for z>0

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.





 §.§ Stabilized semi-implicit scheme


The backward Euler scheme (<ref>)-(<ref>) 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 <cit.>.
The resulting stabilized semi-implicit scheme for the problem (<ref>)-(<ref>) reads

    ∂_f^αu^n   =ε^2D_hu^n
    -f(u^n-1)-S(u^n-u^n-1),  n≥1,
    
    u_h^0   =u_0(𝐱_h),  𝐱_h∈Ω̅_h.


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

If u^0_∞≤1 and the scalar stabilized parameter S≥2,
then the solution of stabilized semi-implicit scheme (<ref>)-(<ref>) satisfies

    u^k_∞≤1  for 0≤k≤N.

So it preserves the maximum principle (<ref>) numerically
and is unconditionally stable.


It only needs to verify that u^n_∞≤1 under the induction hypothesis

    u^k_∞≤1  for 0≤k≤n-1.

From the linearized scheme (<ref>), one has

    (A_0^(n)+S)u^n-ε^2D_hu^n_∞=(1+S)u^n-1-(u^n-1)^3+L^n-1_∞,

where L^n-1 is defined in (<ref>).
Thanks to the decreasing property in Lemma <ref> (i),
the induction hypothesis
and the triangle inequality yield L^n-1_∞≤ A_0^(n).
Furthermore, it is easy to check that 

    (1+S)z-z^3≤ S  if z≤1 and S≥2,

thus the right hand side of (<ref>) can be bounded by

    (1+S)u^n-1-(u^n-1)^3+L^n-1_∞≤ A_0^(n)+
    (1+S)u^n-1-(u^n-1)^3_∞≤ A_0^(n)+S.

For the left hand side of (<ref>), we apply Lemma <ref> (c) and Lemma <ref>
to find that

    (A_0^(n)+S)u^n
    -ε^2D_hu^n_∞≥
    (A_0^(n)+S)u^n_∞.

Then the desired estimate u^n_∞≤1 follows from (<ref>) directly.


Due to the presence of the stabilized term S(u^n-u^n-1), the numerical solution generated by
the semi-implicit scheme
(<ref>)-(<ref>) 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.





§ 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

    u(t)_W^4,∞(Ω)≤ C_u,  u^(ℓ)(t)_W^0,∞(Ω)≤ C_u1+t^σ-ℓ  for 0<t≤ T and ℓ=1,2,

where a regularity parameter σ∈(0,1) is introduced to make our analysis extendable.

In <cit.>, the local consistency error
Υ^j:=(∂_t^αu)(t_j)-(∂_f^αu)^j
of fast L1 formula (<ref>) 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 ∑_j=1^np_n-j^(n)Υ^j,
where p_n-j^(n) are the discrete complementary convolution kernels
defined via (<ref>).








Note that, the definition (<ref>) and Lemma <ref> (i)
show that the discrete convolutional kernels A_n-k^(n) fulfill two assumptions
Ass1-Ass2 in Appendix <ref> with π_a=3/2. In this section,
we will use the results
of Lemma <ref> without further declarations.




Under the condition of Lemma <ref>,
the global consistency error is bounded by

    ∑_j=1^np_n-j^(n)Υ^j≤∑_k=1^np_n-k^(n)A_0^(k)G^k
    +∑_k=1^n-1p_n-k^(n)A_0^(k)G^k
    +C_u/σt_n^αt̂_n-1^2ϵ  for n≥1,

where the local quantities G^k :=2∫_t_k-1^t_k⟨t-t_k-1|u_tt t for 1≤ k≤ n and t̂_n:=max{1,t_n}.

 On the basis of the upper bound of (∂_t^α)u(t_n)-(∂_τ^αu)^n given in <cit.>,
  the estimate (3.5) in the proof of <cit.> gives the desired result.



To resolve such a solution u efficiently, it is appropriate to choose the time
mesh such that the following
condition <cit.> holds.


  AssG. Let γ≥1 be a user-chosen parameter. There is a mesh-independent constant C_γ>0 such that
τ_k≤ C_γτmin{1,t_k^1-1/γ}
for 1≤ k≤ N and
t_k≤ C_γt_k-1 for 2≤k≤N.

Here, the parameter γ≥1 controls the extent to which the time levels
are concentrated near t=0. If the mesh is quasi-uniform, then
AssG holds with γ=1.  As γ 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)^γ with the maximum step ratio ρ=1.


It is to note that, the global consistency error in
Lemma <ref>
gives a superconvergence estimate of nonuniform L1 formula.
Consider the first time level n=1, the regularity setting
(<ref>) gives
Υ^1≤ C_u A_0^(1)∫_0^t_1t^σ-1 t≤
C_uτ_1^σ-α/σ,




implying that the L1 formula is always inconsistent if 0<σ≤α,
also see Table <ref> in Section 4.
However, we have the global consistency error of order O(τ_1^σ), because
p_0^(1)Υ^1≤ G^1≤ C_uτ_1^σ/σ.
In general, we have the following result from <cit.>.

Under the regularity (<ref>), the global consistency error can be bounded by

    ∑_j=1^np_n-j^(n)Υ^j≤ C_u τ_1^σ/σ
    +1/1-αmax_2≤k≤nt_k^αt_k-1^σ-2τ_k^2-α
    +ϵ/σt_n^αt̂_n-1^2  for 1≤n≤N.

Specifically, if the mesh satisfies AssG, then

    ∑_j=1^np_n-j^(n)Υ^j≤C_u/σ(1-α)τ^min{2-α,γσ}
    +C_uϵ/σt_n^αt̂_n-1^2  for 1≤n≤N.






Assume that u_0_L^∞(Ω)≤1 and
the solution of (<ref>)-(<ref>) satisfies the regular assumption (<ref>).
If the maximum step size τ≤1/√(6Γ(2-α)), then the numerical solution u_h^n
of the backward Euler scheme (<ref>)-(<ref>) is convergent in the maximum norm, that is,

    u(𝐱_h,t_n)-u_h^n_∞≤ C_uτ_1^σ/σ
    +1/1-αmax_2≤k≤nt_k^αt_k-1^σ-2τ_k^2-α
    +ϵ/σt_n^αt̂_n-1^2
    +h_1^2+h_2^2

for 1≤n≤N. Moreover, when the time mesh satisfies AssG, it holds that

    u(𝐱_h,t_n)-u_h^n_∞≤C_u/σ(1-α)⟨τ^min{2-α,γσ}+ϵ|
    +C_u⟨h_1^2+h_2^2|  for 1≤n≤N,

which achieves the optimal accuracy O(τ^2-α) if the graded parameter γ≥max{1, (2-α)/σ}.


 Let U_h^n:=u(𝐱_h,t_n) and the error function e_h^n:=U_h^n-u_h^n∈𝕍_h
for 𝐱_h∈Ω̅_h and 0≤n≤N.
It is easy to find that the exact solution U_h^n satisfies the governing equations

    ⟨∂_f^αU|^n-ε^2D_hU^n   =-f(U^n)+⟨R_t|^n+⟨R_s|^n,  1≤n≤N,
    
    U_h^0   =u_0(𝐱_h),  𝐱_h∈Ω_h,

where ⟨R_t|^n and ⟨R_s|^n denote the truncation errors in time and space, respectively.
Subtracting (<ref>)-(<ref>) from the above two equations, respectively,  one gets

    ⟨∂_f^αe|^n-ε^2D_he^n   =
    -f(U^n)+f(u^n)+⟨R_t|^n+⟨R_s|^n,  1≤n≤N,
    
    e_h^0   =0,  𝐱_h∈Ω_h.

Recalling the elementary inequality |(a^3-a)-(b^3-b)|≤2|a-b| for ∀  a,b∈[-1,1],



we apply Theorem <ref> (discrete maximum principle) to get

    f(U^n)-f(u^n)_∞≤ 2e^n_∞.

Thus the triangle inequality with the error equation (<ref>) gives

    ⟨∂_f^αe|^n-ε^2D_he^n_∞≤ 2e^n_∞+⟨R_t|^n_∞+⟨R_s|^n_∞.

Applying the decreasing property (i) of the kernels A_n-k^(n) and the triangle inequality,
we can bound the left hand side of (<ref>) by

    ⟨∂_f^αe|^n-ε^2D_he^n_∞   =(A_0^(n)-ε^2D_h)e^n
    -∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)e^k-A_0^(n)e^0_∞
       ≥(A_0^(n)-ε^2D_h)e^n_∞
    -∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)e^k_∞
    -A_n-1^(n)e^0_∞
       ≥∑_k=1^nA_n-k^(n)▿_τe^k_∞,

where Lemma <ref> (c) and Lemma <ref> have been used.
Then it follows from (<ref>) that

    ∑_k=1^nA_n-k^(n)▿_τe^k_∞≤ 2e^n_∞+⟨R_t|^n_∞+⟨R_s|^n_∞,

which takes the form of (<ref>) with the substitutions λ:=2, v^k:=e^k_∞,
ξ^n:=⟨R_t|^n_∞ and η^n:=⟨R_s|^n_∞.




Lemma <ref> (the discrete fractional Grönwall inequality) says that, if the maximum step size τ≤1/√(6Γ(2-α)), then it holds that

    e^n_∞≤
    2E_α6max(1,ρ)t_n^α(max_1≤k≤n∑_j=1^kp_k-j^(k)⟨R_t|^j_∞
    +ω_1+α(t_n)max_1≤k≤n⟨R_s|^k_∞).

Then Corollary <ref> yields the claimed estimate and completes the proof.





For the semi-implicit scheme (<ref>)-(<ref>), the global error is dominated
by the stabilized term S(u^k-u^k-1). Under the regular assumption (<ref>),
the local consistency error is about ∫_t_k-1^t_ku_tt. One can follow the proof of
<cit.> to bound
the corresponding global error as follows (also see the case of m=0 in the estimate (<ref>))

    ∑_j=1^np_n-j^(n)∫_t_j-1^t_ju_tt≤ C_uτ_1^σ/σ
    +1/1-αmax_2≤k≤nt_k^αt_k-1^σ-1τ_k.

Then, a similar proof of Theorem <ref> leads to the following result.

Assume that u_0_L^∞(Ω)≤1 and
the exact solution of (<ref>)-(<ref>) satisfies the regular assumption (<ref>).
If the stabilized parameter S≥2 and the maximum time-step size τ≤1/√(6Γ(2-α)), then
the numerical solution u_h^n
of the semi-implicit scheme (<ref>)-(<ref>) is convergent in the maximum norm, that is,

    u(𝐱_h,t_n)-u_h^n_∞≤ C_uτ_1^σ/σ
    +1/1-αmax_2≤k≤nt_k^αt_k-1^σ-1τ_k
    +ϵ/σt_n^αt̂_n-1^2
    +h_1^2+h_2^2

for 1≤n≤N. Moreover, when the time mesh satisfies AssG, it holds that

    u(𝐱_h,t_n)-u_h^n_∞≤C_u/σ(1-α)τ^min{1,γσ}+ϵ
    +C_u⟨h_1^2+h_2^2|  for 1≤n≤N,

which achieves the optimal accuracy O(τ) if the graded parameter γ≥max{1, 1/σ}.



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



    a_n-k^(n)=a_n-k=1/τ^αω_2-α(n-k+1)-ω_2-α(n-k)  for 1≤k≤n,

the semi-implicit stabilized scheme (<ref>)-(<ref>) using the L1 formula inherits a discrete energy dissipation law,
see  <cit.> for details.
As seen, the proof of discrete energy dissipation law relies on
the property of a quadratic form ∑_k=1^nw_k∑_j=1^ka_k-jw_j≥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=▿_τv^k)

    ∑_k=1^n▿_τv^k(∂_τ^αv)^k
    =∑_k=1^nw_k∑_j=1^ka_k-j^(k)w_j≥0.

More generally, it has yet to be determined what restrictions must be imposed on the discrete convolution coefficients {A_n-k^(n) |1≤ k≤ n} so that the quadratic form ∑_k=1^nw_k∑_j=1^kA_k-j^(k)w_j is positive semi-definite.
This problem could be challenging and remains open to us.
















§ NUMERICAL EXAMPLES



The nonuniform fast L1 time-stepping methods (<ref>)-(<ref>) and
(<ref>)-(<ref>) are examined
for solving the Allen-Cahn problem (<ref>)-(<ref>).
Always, we set the absolute tolerance error ϵ=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
(<ref>)-(<ref>), 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≤n≤NU^n-u^n_∞ is recorded in each run,
and the experimental convergence order in time is computed by

    Order:=log⟨e(M,N)/e(M,2N)|/log⟨τ(N)/τ(2N)|

where τ(N) denotes the maximum time-step size for total N subintervals.



To examine the temporal accuracy of our time-stepping schemes, consider the time-fractional Allen-Cahn equation
∂_t^αu=1/8π^2Δu-f(u)+g(𝐱,t) for 𝐱∈(0,1)^2 and 0<t<1
such that it has an exact solution u=ω_1+σ(t)sin(2πx)sin(2π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/γ,T}, and apply the graded grid t_k=T_0(k/N_0)^γ
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

    τ_N_0+k=(T-T_0)ϵ_k/∑_k=1^N_1ϵ_k  for 1≤ k≤ N_1

where ϵ_k∈(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 (<ref>)-(<ref>)
for two different cases σ=α and σ<α are listed in Tables <ref>-<ref>, respectively. They suggest the time accuracy is of order
  O(τ^min{γσ,2-α}) and confirm Theorem <ref> experimentally.
We also run the stabilized semi-implicit scheme (<ref>)-(<ref>) by setting a variety of regularity parameters.
Tables <ref>-<ref> report
the numerical results in the case σ=α and a worse case of σ<α.
It seen that it is accurate of order O(τ^min{γσ,1}) on the graded meshes,
confirming Theorem <ref> experimentally.


[Coalescence of two kissing bubbles]
Consider the time-fractional Allen-Cahn problem (<ref>)-(<ref>)
describing the coalescence of two kissing bubbles inside the spatial domain Ω=(-π, π)^2,
by taking ε=0.1 and the initial data

    u_0(𝐱)
    =
    0.5,    (x+1)^2+y^2≤1 or (x-1)^2+y^2≤1,
    
    -0.5,    otherwise.

























This example is used to examine the physical effect of the fractional order α in the original problem
and the physical property of our suggested methods.
Theorems <ref> and <ref> 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 <cit.>

    τ_k=min{max{τ_min, tol/1+βu^k-u^k-1_∞},τ_max}  for k≥1,

where the constant 1 is set to avoid the possible singularity as the model reaches the steady state.
The parameters tol and β are used to adjust the level of adaptively and would be chosen in experience.
A small tol or a large β will generate time steps close to τ_min, which a large tol or a small β will give time steps close to τ_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)^γ with the graded parameter γ=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 (<ref>)-(<ref>) with τ_min=τ_N_0=0.001,τ_max=0.1,tol=0.15,β=200
and stabilized semi-implicit scheme (<ref>)-(<ref>) with S=0.1, τ_min=τ_N_0=0.001,τ_max=1,tol=1.5,β=200 in the remainder interval, are quite identical. Fig. <ref> 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 α, see <cit.>.
The larger the fractional order α, the faster the bubbles coalesce.



Fig.<ref> 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 (<ref>)-(<ref>).
It is obvious that the solutions are uniformly bounded by the value 1 for different fractional orders α,
as predicted by Theorem <ref>.
Moreover, the larger the fractional order α, the faster it approaches the maximum value.
The middle of Fig. <ref> says that the discrete energy
is also decreasing as the time escapes, although we can not verify it theoretically.
The right side of Fig. <ref> 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.<ref> shows analogous plots for the stabilized scheme (<ref>)-(<ref>),
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 τ_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.






§ 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.





§ DISCRETE FRACTIONAL GRÖNWALL LEMMA


The recently developed discrete fractional Grönwall inequality
in <cit.> 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_α(z):=∑_k=0^∞z^k/Γ(1+kα),
gathers three previous (slightly simplified) results from
<cit.>.

For n=1,2,⋯,N, assume that the discrete convolution kernels {A_n-k^(n)}_k=1^n satisfy the following two assumptions:

Ass1. There is a constant π_a>0 such that A_n-k^(n)≥1/π_a∫_t_k-1^t_kω_1-α(t_n-s)/τ_k s for 1≤ k≤ n.

Ass2. The discrete kernels are monotone, i.e. A_n-k-1^(n)-A_n-k^(n)≥ 0 for 1≤ k≤ n-1.

Define also a sequence of discrete complementary convolution kernels {p_n-j^(n)}_j=1^n by

    p_0^(n):=1/A_0^(n),  
    p_n-j^(n):=
    1/p_0^(j)∑_k=j+1^nA_k-j-1^(k)-A_k-j^(k)p_n-k^(n),
           1≤ j≤ n-1.

Then the discrete  complementary kernels p^(n)_n-j≥0 are well-defined and fulfill

    ∑_j=k^np^(n)_n-jA_j-k^(j)= 1,
      for 1≤ k≤ n≤ N.
       ∑_j=1^np^(n)_n-jω_1+mα-α(t_j)≤π_aω_1+mα(t_n),
      for m=0,1 and 1≤ n≤ N.

Suppose that the offset parameter 0 ≤ν <1, λ is a non-negative constant independent of the time-steps
and the maximum step size τ≤1/√(2Γ(2-α)λπ_a).
If the non-negative sequences (v^k)_k=0^N, (ξ^k)_k=1^N and (η^k)_k=1^N satisfy

    ∑_k=1^nA_n-k^(n)▿_τ v^k≤λ v^n-ν+ξ^n+η^n  for 1≤ n≤ N,

or

    ∑_k=1^na^(n)_n-k▿_τ (v^k)^2≤λ (v^n-ν)^2+v^n-ν(ξ^n+η^n)  for 1≤ n≤ N,

then it holds that, for 1≤ n≤ N,

    v^n   ≤2E_α(2max{1,ρ}λπ_a t_n^α)
    	(v^0+max_1≤ k≤ n∑_j=1^kp^(k)_k-j(ξ^j+η^j))
       ≤2E_α(2max{1,ρ}λπ_a t_n^α)
    	(v^0+Γ(1-α)π_amax_1≤ k≤ n{t_k^αξ^k}+π_aω_1+α(t_n)max_1≤ k≤ nη^k).











 



10

Allen1979A
M. Allen and W. Cahn.
A microscopic theory for antiphase boundary motion and its
  application to antiphase domain coarsening.
Acta Metall., 27:1085–1095, 1979.

Cahn1958Free
J. Cahn and J. Hilliard.
Free energy of a nonuniform system I. interfacial free energy.
J. Chem. Phys., 28:258–267, 1958.

Clarke1987Origin
S. Clarke and D. Vvedensky.
Origin of reflection high-energy electron-diffraction intensity
  oscillations during molecular-beam epitaxy: a computational modeling
  approach.
Phys. Rev. Let., 58:2235–2238, 1987.

Zhao2016A
J. Zhao, Y. Shen, M. Haapasalo, Z. Wang, and Q. Wang.
A 3D numerical study of antimicrobial persistence in heterogeneous
  multi-species biofilms.
J. Theor. Bio., 392:83–98, 2016.

Shen2012Numerical
J. Shen and X. Yang.
Numerical approximations of Allen-Cahn and Cahn-Hilliard equations.
Dis. Contin. Dyn. Sys. Ser. A, 28:1669–1691, 2012.

Lee2016Comparison
D. Lee and J. Kim.
Comparison study of the conservative Allen-Cahn and the
  Cahn-Hilliard equations.
Math. Comput. Simu., 119:35–56, 2016.

Kim2017A
J. Kim, D. Jeong, S. Yang, and Y. Choi.
A finite difference method for a conservative Allen-Cahn equation
  on non-flat surfaces.
J. Comput. Phys., 334:170–181, 2017.

QiaoAn2011
Z. Qiao, Z. Zhang, and T. Tang.
An adaptive time-stepping strategy for the molecular beam epitaxy
  models.
SIAM J. Sci. Comput., 33:1395–1414, 2011.

Zhang2013An
Z. Zhang, Y. Ma, and Z. Qiao.
An adaptive time-stepping strategy for solving the phase field
  crystal model.
J. Comput. Phys., 249:204–215, 2013.

Li2017Computationally
Y. Li, Y. Choi, and J. Kim.
Computationally efficient adaptive time step method for the
  Cahn-Hilliard equation.
Comput. Math. Appl., 73:1855–1864, 2017.

Inc2018Time
M. Inc, A. Yusuf, A. Aliyu, and D. Baleanu.
Time-fractional Cahn-Allen and time-fractional Klein-Gordon
  equations: Lie symmetry analysis, explicit solutions and convergence
  analysis.
Physica A Stat. Mech. Appl., 493:94–106, 2018.

Akagi2016Fractional
G. Akagi, G. Schimperna, and A. Segatti.
Fractional Cahn-Hilliard, Allen-Cahn and porous medium equations.
J. Differ. Equations, 261:2935–2985, 2016.

Hou2017Numerical
T. Hou, T. Tang, and J. Yang.
Numerical analysis of fully discretized Crank-Nicolson scheme for
  fractional-in-space Allen-Cahn equations.
J. Sci. Comput., 72:1–18, 2017.

Zheng2017A
Z. Li, H. Wang, and D. Yang.
A space-time fractional phase-field model with tunable sharpness and
  decay behavior and its efficient numerical simulation.
J. Comput. Phys., 347:20–38, 2017.

Liu2018Time
H. Liu, A. Cheng, H. Wang, and J. Zhao.
Time-fractional Allen-Cahn and Cahn-Hilliard phase-field
  models and their numerical investigation.
Comp. Math. Appl., 76:1876–1892, 2018.

Zhao2019On
J. Zhao, L. Chen, and H. Wang.
On power law scaling dynamics for time-fractional phase field models
  during coarsening.
Comm. Non. Sci. Numer. Simu., 70:257–270, 2019.

Tang2018On
T. Tang, H. Yu, and T. Zhou.
On energy dissipation theory and numerical stability for
  time-fractional phase field equations.
arXiv:1808.01471v1, 2018.

Jin2016An
B. Jin, R. Lazarov, and Z. Zhou.
An analysis of the L1 scheme for the subdiffusion equation with
  nonsmooth data.
IMA J. Numer. Anal., 36:197–221, 2016.

JinLiZhou:2017
B. Jin, B. Li, and Z. Zhou.
Numerical analysis of nonlinear subdiffusion equations.
SIAM J. Numer. Anal., 56:1–23, 2018.

Liao2018Sharp
H.-L. Liao, D. Li, and J. Zhang.
Sharp error estimate of nonuniform L1 formula for time-fractional
  reaction-subdiffusion equations.
SIAM J. Numer. Anal., 56:1112–1133, 2018.

Liao2018Unconditional
H.-L. Liao, Y. Yan, and J. Zhang.
Unconditional convergence of a fast two-level linearized algorithm
  for semilinear subdiffusion equations.
J. Sci. Comput., 2019.
DOI: 10.1007/s10915-019-00927-0.

Liao2018second
H.-L. Liao, W. Mclean, and J.Zhang.
A second-order scheme with nonuniform time steps for a linear
  reaction-sudiffusion problem.
arXiv:1803.09873v2, 2018.
in review.

Liao2018discrete
H.-L. Liao, W. Mclean, and J. Zhang.
A discrete Grönwall inequality with application to numerical
  schemes for subdiffusion problems.
SIAM J. Numer. Anal., 57:218–237, 2019.

Jiang2017Fast
S. Jiang, J. Zhang, Z. Qian, and Z. Zhang.
Fast evaluation of the Caputo fractional derivative and its
  applications to fractional diffusion equations.
Comm. Comput. Phys., 21:650–678, 2017.

William2007A
W. McLean and K. Mustapha.
A second-order accurate numerical method for a fractional wave
  equation.
Numer. Math., 105:481–510, 2007.




