Unconditional convergence of a fast two-level linearized algorithm for semilinear subdiffusion equations
Hong-lin Liao Yonggui Yan Jiwei Zhang









 Unconditional convergence of a fast two-level linearized algorithm
for semilinear subdiffusion equations
    Hong-lin Liao
Department of Mathematics, Nanjing University of Aeronautics and Astronautics,
Nanjing, 211106, P. R. China.
E-mail: liaohl@csrc.ac.cn. He is partially supported by NSFC grant 11372354
and a grant DRA2015518 from 333 High-level Personal Training Project of Jiangsu Province.





  Yonggui Yan

Beijing Computational Science Research Center (CSRC), Beijing, 100193, P. R. China.
 E-mail: yan_yonggui@csrc.ac.cn. 
  Jiwei Zhang

Beijing Computational Science Research Center (CSRC),
Beijing, 100193, P. R. China.  E-mail: jwzhang@csrc.ac.cn.
He is partially supported by NSFC under grants 11771035, 91430216 and NSAF U1530401.

    Received date / Accepted date
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
















A fast two-level linearized scheme with unequal time-steps is constructed and analyzed
for an initial-boundary-value problem of semilinear subdiffusion equations.
The two-level fast L1 formula of the Caputo derivative is derived based on
the sum-of-exponentials technique. The resulting fast algorithm
is computationally efficient in long-time simulations because it
significantly reduces the computational cost O(MN^2) and storage O(MN)
for the standard L1 formula to O(MNlog N) and O(Mlog N),
respectively, for M grid points in space and N levels in time.
The nonuniform time mesh would be graded to handle the typical singularity
of the solution near the time t=0,
and Newton linearization is used to approximate the nonlinearity term.
Our analysis relies on three tools:
a new discrete fractional Grönwall inequality, a global
consistency analysis and a discrete H^2 energy method.
A sharp error estimate reflecting the regularity of solution is established without
any restriction on the relative diameters of the temporal and spatial mesh sizes.
Numerical examples are provided to demonstrate the effectiveness of our approach
and the sharpness of error analysis.



Keywords:  semilinear subdiffusion equation; two-level L1 formula;
discrete fractional Grönwall inequality;
discrete H^2 energy method; unconditional convergence








§ INTRODUCTION




A two-level linearized method is considered to numerically solve the following semilinear subdiffusion
equation on a bounded domain


    αu    = Δ u+f(u)  for ∈Ω and 0<t≤ T,
    
    u   =u^0()   for ∈Ω and t=0,
    
    u   = 0  for ∈∂Ω and 0<t≤ T,


where ∂Ω is the boundary of Ω:=(x_l,x_r)× (y_l,y_r),
and the nonlinear function f(u) is smooth.
In (<ref>) α=_ 0^C𝒟_t^α denotes
the  Caputo  fractional derivative of order α:   

    αv(t)
    ∫_0^tω_1-α (t-s) v'(s)s,  
    0<α<1,

where the weakly singular kernel ω_1-α (t-s) is defined by
ω_μ(t):=t^μ-1/Γ(μ). It is easy to verify
ω_μ'(t)=ω_μ-1(t) and ∫_0^tω_μ(s) s=ω_μ+1(t) for t>0.

In any numerical methods for solving nonlinear fractional diffusion equations (<ref>),
a key consideration is the singularity of the solution near the time t=0, see <cit.>.
For example, under the assumption that the nonlinear function f is Lipschitz continuous
and the initial data u^0∈ H^2(Ω)∩ H_0^1(Ω),

Jin et al. <cit.> prove that problem (<ref>) has a unique solution u for which
u∈ C⟨[0,T];H^2(Ω)∩ H_0^1(Ω)|, αu∈ C⟨[0,T];L^2(Ω)| and ∂_t u ∈ L^2(Ω)  with
∂_t u(t)_L^2(Ω)≤ C_ut^α-1 for 0<t≤ T,



where C_u>0 is a constant independent of t but may depend on T.
Their analysis of numerical methods for solving (<ref>) is applicable to both the L1 scheme
and backward Euler convolution quadrature on a uniform time grid of diameter τ;
 a lagging linearized technique is used to handle the nonlinearity f(u), and <cit.>
 shows that the discrete solution  is O(τ^α) convergent in L^∞(L^2(Ω)).

This work may be considered as a continuation of <cit.>, in which a sharp error estimate
for the L1 formula on nonuniform meshes was obtained for linear subdiffusion-reaction equations based on
a discrete fractional Grönwall inequality and a global consistency analysis. In this paper,
we combine the L1 formula and the sum-of-exponentials (SOEs) technique to develop a
one-step fast difference algorithm for the nonlinear subdiffusion problem (<ref>) by
using the Newton's linearization to approximate nonlinear term,
and present the corresponding sharp error estimate of the proposed scheme
without any restriction on the relative diameters of temporal and spatial mesh sizes.

It is known that the Caputo fractional derivative involves a convolution kernel.
The total number of operations required to evaluate the sum of L1 formula is proportional
to O(N^2), and the active memory to O(N) with N representing the total time steps,
which is prohibitively expensive for the practically large-scale and long-time simulations.
Recently, a simple fast algorithm based on SOEs approximation is proposed to significantly
reduce the computational complexity to O(N log N) and O(log N) when the final time T≫ 1,
see <cit.>.  Another fast algorithm for the evaluation of
the fractional derivative has been proposed in <cit.>, where the compression is carried
out in the Laplace domain by solving the equivalent ODE with some one-step A-stable scheme.
In this paper,  we develop a fast two-level L1 formula
by combining a nonuniform mesh suited to the initial singularity
with a fast time-stepping algorithm for the historical memory in (<ref>).
This scheme would be also useful to develop efficient parallel-in-time  algorithms
for time-fractional differential equations  <cit.>.

On the other hand, the nonlinearity of the problem also results in the difficulty for the numerical analysis.
To establish an error estimate of the two-level linearized scheme at time t_n, it requires to prove
the boundedness of the numerical solution at the previous time levels via
u^n-1_∞≤ C_u.
Traditionally it is done using mathematical induction and some inverse estimate, namely,

    u^n-1_∞≤    
    U^n-1_∞+h^-1U^n-1-u^n-1≤U^n-1_∞+C_uh^-1τ^β+h^2.

This leads to that a time-space grid
restriction τ=O(h^1/β) is required in the theoretical analysis even though
it is nonphysical and may be unnecessary in numerical simulations. In this paper, we will
extend the discrete H^2 method developed in <cit.>
to prove unconditional convergence of our fully discrete solution without the restriction conditions of
between mesh sizes τ and h comparing with the traditional method.  The main idea of
discrete H^2 energy method is to separately treat the temporal and spatial truncation errors.
This simple implementation avoids some nonphysical time-space grid restrictions in the error analysis.
A related approach in a finite element setting are discussed in <cit.>.

The convergence rate of L1 formula for the Caputo derivative is limited by the smoothness of the solution.
The analysis here is based on the following assumptions on the solution

    u_H^4(Ω)≤ C_u, ∂_tu_H^4(Ω)≤ C_u(1+t^σ-1)  and ∂_ttu_H^2(Ω)≤ C_u(1+t^σ-2)

for 0< t ≤ T, where σ∈(0,1)∪(1,2) is a regularity parameter.
To resolve the singularity at t=0, it is reasonable to use a nonuniform mesh that concentrates  grid points near t=0,
see <cit.>.
We make the following assumption on the time mesh:


  G. Let γ≥1 be a user-chosen parameter. There is a constant C_γ>0, independent of k,  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.

Since τ_1=t_1, G implies that τ_1=O(τ^γ), while for those t_k bounded
away from t=0 one has τ_k=O(τ). The parameter γ controls the extent to which the
grid points are concentrated near t=0: increasing γ will decrease the time-step sizes
near t=0 and so move mesh points closer to t=0.
A simple example of a family of meshes satisfying G is the graded grid t_k=T(k/N)^γ,
which is discussed in <cit.>.  Although nonuniform meshes are flexible and reasonably convenient for practical
implementation, they can significantly complicate the numerical analysis of
schemes, both with respect to stability and consistency. In this paper, our analysis will rely on a generalized fractional Grönwall inequality <cit.>,
which would be applicable for any discrete fractional derivatives having the discrete convolution form.

Throughout the paper, any subscripted C, such as C_u, C_γ, C_Ω,
 C_v, C_0 and C_F, denotes a generic positive constant, not necessarily
 the same at different occurrences, which is always dependent on the given data
 and the solution but independent of the time-space grid steps.
 The  paper is organized as follows. Section <ref> presents the
two-level fast L1 formula and the corresponding linearized fast scheme.
The global consistency analysis of fast L1 formula and the Newton's linearization
is presented in Section <ref>. A sharp error estimate
for the linearized fast scheme is proved in Section <ref>.
Two numerical examples in Section <ref> are given to demonstrate the sharpness of our analysis.




§ A TWO-LEVEL FAST METHOD



We approximate the Caputo fractional  derivative (<ref>) on a (possibly nonuniform) time mesh
0=t_0<⋯<t_k-1<t_k<⋯<t_N=T, with the time-step sizes
τ_k t_k-t_k-1 for 1≤ k≤ N, the maximum time-step =max_1≤ k≤ Nτ_k
and the step size ratios ρ_kτ_k/τ_k+1 for 1≤ k≤ N-1.
In space we use a standard finite difference method on a tensor product grid.
 Let M_1 and M_2 be two positive integers. Set h_1=(x_r-x_l)/M_1,  h_2=(y_r-y_l)/ M_2
 and the maximum spatial length h=max{h_1,h_2}. Then the fully discrete spatial grid
 Ω̅_h :={ _h=(x_l + ih_1, y_l + jh_2) |  0≤ i≤ M_1, 0≤ j≤ M_2}. Set Ω_h=Ω̅_h∩Ω and the boundary ∂Ω_h=Ω̅_h∩∂Ω. Given a grid function v = {v_ij}, define

    v_i-1/2,j=⟨v_i,j+v_i-1,j|/2,  δ_xv_i-1/2,j=⟨v_i,j-v_i-1,j|/h_1,  δ_x^2v_ij=δ_xv_i+1/2,j-δ_xv_i-1/2,j/h_1.

Difference operators v_i,j-1/2,  δ_yv_i,j-1/2,  δ_xδ_yv_i-1/2,j-1/2
and δ_y^2v_ij can be defined analogously.
The second-order approximation of Δ v( _h) for _h∈Ω_h
is Δ_hv_h := (δ_x^2+δ_y^2)v_h.
Let 𝒱_h be the space of grid functions,
𝒱_h={v=⟨v_h|_ _h∈Ω̅_h | 
v_h=0  _h∈∂Ω_h}.
 For v, w ∈𝒱_h, define the discrete inner product
 v,w=h_1h_2∑__h∈Ω_hv_hw_h,
 the L^2 norm v=√(v,v),
 the H^1 seminorm ∇_h v=√(δ_xv^2+δ_yv^2) and the maximum norm
v_∞=max__h∈Ω_hv_h. For any v∈𝒱_h,
by <cit.> there exists a constant C_Ω>0 such that

    v≤ C_Ω∇_hv,
      ∇_hv≤ C_ΩΔ_hv,
      v_∞≤ C_ΩΔ_hv.




 §.§ A fast variant of the L1 formula


On our nonuniform mesh, the standard L1 approximation of the Caputo derivative is

    (αv)^n∑_k=1^n1/τ_k∫_t_k-1^t_kω_1-α(t_n-s)∇_τ v^k s
    = ∑_k=1^na_n-k^(n)∇_τv^k ,

where ∇_τv^k:=v^k-v^k-1 and the convolution kernel a_n-k^(n) is defined by

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






For fixed integer n≥2, the convolution kernel a^(n)_n-k of (<ref>) satisfies

  
  (i) a_n-k-1^(n)>ω_1-α(t_n-t_k)>a_n-k^(n),   1≤ k≤ n-1;
  

  
  (ii)a^(n)_n-k-1-a^(n)_n-k>12ω_1-α(t_n-t_k)
  -ω_1-α(t_n-t_k-1),   1≤ k≤ n-1.


 The integral mean-value theorem yields (i) directly;
see <cit.>.
For any function q∈ C^2[t_k-1,t_k], let Π_1,kq be the linear interpolant
of q(t) at  t_k-1 and t_k. Let
1,kq := q - Π_1,kq be the error in this interpolant.
For q(s)=ω_1-α(t_n-s) one has q”(s)=ω_-α-1(t_n-s)>0 for 0<s<t_n, so
the Peano representation of the interpolation error <cit.>  shows that
∫_t_k-1^t_k1,kq(s) s<0.
Thus the definition (<ref>) of a^(n)_n-k yields

    a^(n)_n-k-1/2ω_1-α(t_n-t_k)
        -1/2ω_1-α(t_n-t_k-1)
        =1/τ_k∫_t_k-1^t_k1,kq(s) s<0,   1≤ k≤ n-1.

Subtract this inequality from (i) to obtain (ii) immediately.


As the L1 formula (<ref>) involves the solution at all previous time-levels,
it is computationally inefficient to directly evaluate it  when solving
the fractional diffusion problem (<ref>) using time-stepping.
We therefore use the SOEs approach
of <cit.> to develop
a fast L1 formula. A basic result of the SOE approximation
(see  <cit.>
or <cit.>) is the following:







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

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

where the number N_q  of quadrature nodes satisfies

    N_q=O⟨log1/ϵloglog1/ϵ+logT/Δt
    +log1/Δtloglog1/ϵ+log1/Δt|.



After that, we divide the fractional Caputo  derivative αv(t_n)
of (<ref>) into a sum of a local part (an integral over [t_n-1, t_n]) and
a history part (an integral over [0, t_n-1]), then approximate v' by linear interpolation in the local part (similar to the standard L1 method) and use
the SOE technique of Lemma <ref> to approximate the kernel
ω_1-α (t-s) in the history part.  It yields

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

where ℋ^ℓ(t_k) := ∫_0^t_ke^-θ^ℓ(t_k-s)u'(s) s
with ℋ^ℓ(t_0)=0 for 1≤ℓ≤ N_q.











To compute ℋ^ℓ(t_k) efficiently we apply linear interpolation
in each cell [t_k-1, t_k], obtaining

    ℋ^ℓ(t_k)=    e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +∫_t_k-1^t_ke^-θ^ℓ(t_k-s)u'(s) s
    ≈ e^-θ^ℓτ_kℋ^ℓ(t_k-1)
    +b^(k,ℓ)∇_τu^k,

where the positive coefficient is given by

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


In summary, we now have the two-level fast L1 formula


    αu^n
      a_0^(n)∇_τu^n+∑_ℓ=1^N_qϖ^ℓe^-θ^ℓτ_nH^ℓ(t_n-1),    n≥1,
    where H^ℓ(t_k) satisfies H^ℓ(t_0)=0 and the recurrence relationship
    H^ℓ(t_k)=e^-θ^ℓτ_kH^ℓ(t_k-1)
    +b^(k,ℓ)∇_τu^k,   k≥1,  1≤ℓ≤ N_q .







 §.§ The two-level linearized scheme




Write U_h^n=u(_h,t_n) for _h∈Ω̅_h, 0≤ n≤ N.
Let u_h^n be the discrete approximation of U_h^n.
Using the fast L1 formula (<ref>) and Newton linearization, we obtain a
linearized scheme for the problem (<ref>): find  {u_h^N}∈𝒱_h such that


    αu_h^n=    Δ_hu_h^n+f(u_h^n-1)
    +f'(u_h^n-1)∇_τu_h^n ,  _h∈Ω_h,   1≤ n≤ N;
    
    u_h^0=    u^0(_h),  _h∈Ω̅_h .



Note that, the Newton linearization of a general nonlinear function f=f(,t,u) at t=t_n takes
the form f(_h,t_n,u_h^n)≈ f(_h,t_n,u_h^n-1) +f_u'(_h,t_n,u_h^n-1)∇_τu_h^n .
The scheme (<ref>) is a two-level procedure for computing {u_h^n}, since
 (<ref>) can be reformulated as

    a_0^(n)-Δ_h-f'(u_h^n-1)∇_τu_h^n
    =    Δ_hu_h^n-1+f(u_h^n-1)
    -∑_ℓ=1^N_qϖ^ℓe^-θ^ℓτ_nH_h^ℓ(t_n-1),
    
    H_h^ℓ(t_n)=    e^-θ^ℓτ_nH_h^ℓ(t_n-1)
    +b^(n,ℓ)∇_τu_h^n ,    1≤ℓ≤ N_q.

Thus, once the solution {u_h^n-1, H_h^ℓ(t_n-1)} at the previous time-level t_n-1 is available,
the current solution {u_h^n} can be found by (<ref>) with a fast matrix solver
and the historic term {H_h^ℓ(t_n)} will be updated explicitly by the recurrence formula (<ref>).


At each time level the scheme (<ref>) requires O(MN_q) storage and O(MN_q) operations,
where M=M_1M_2 is the total number of spatial grid points. Given a tolerance error ϵ=ϵ_0,
by virtue of Lemma <ref>, the number of quadrature nodes N_q=O(log N)
if the final time T≫1. Hence our new method is computationally efficient since it computes
the final solution using in total O(Mlog N) storage and O(MNlog N) operations.




 §.§ Discrete fractional Grönwall inequality


Our analysis relies on a generalized discrete fractional Grönwall inequality <cit.>,
which is applicable for any discrete fractional derivative having the discrete convolution form 

    αv^n
    ≈∑_k=1^n A^(n)_n-k(v^k-v^k-1),  
    1≤ n≤ N,

provided that A^(n)_n-k and the time-steps τ_n satisfy the following three assumptions:


  1. The discrete kernel is monotone, that is, A^(n)_k-2≥ A^(n)_k-1>0 for 2≤ k≤ n≤ N.

  2. There is a constant π_A>0 such that A^(n)_n-k≥1/π_A∫_t_k-1^t_kω_1-α(t_n-s)/τ_ks for 1≤ k≤ n≤ N.

  3. There is a constant >0 such that the time-step ratios ρ_k≤
for 1≤ k≤ N-1.


The  complementary discrete kernel P^(n)_n-k was introduced by
Liao et al. <cit.>;
it satisfies the following identity

    ∑_j=k^nP^(n)_n-jA^(j)_j-k≡1  for 1≤ k≤ n≤ N.

Rearranging this identity yields a recursive formula that defines P^(n)_n-k :

    P_0^(n)1/A_0^(n),  
    P_n-j^(n)1/A_0^(j)∑_k=j+1^nA_k-j-1^(k)-A_k-j^(k)P_n-k^(n) ,   1≤ j≤ n-1.

From <cit.> we see that P^(n)_n-k is well-defined
and non-negative if the assumption 1 holds true. Furthermore, if 2 holds true, then

    ∑_j=1^nP^(n)_n-j≤π_A ω_1+α(t_n)  for 1≤ n≤ N.



Recall that the Mittag–Leffler function
E_α(z) = ∑_k=0^∞z^k/Γ(1+kα).
We state the following (slightly simplified) version of <cit.>.
This result differs substantially from
the fractional Grönwall inequality of Jin et al. <cit.>
 since it is valid on very general nonuniform time meshes.

Let 1–3 hold true. Suppose that the sequences (ξ_1^n)_n=1^N, (ξ_2^n)_n=1^N are nonnegative.
Assume that λ_0 and λ_1 are non-negative constants and the maximum step size
τ≤1/√(2max{1,ρ}π_AΓ(2-α)(λ_0+λ_1)).




If the nonnegative sequence (v^k)_k=0^N satisfies

    ∑_k=1^nA^(n)_n-k v^k≤λ_0v^n+λ_1v^n-1+ξ_1^n+ξ_2^n  for 1≤ n≤ N,

then it holds that for 1≤ n≤ N,

    v^n≤2 E_α2max{1,ρ}π_A(λ_0+λ_1)t_n^α(v^0+max_1≤ k≤ n∑_j=1^k P^(k)_k-jξ_1^j
    +π_Aω_1+α(t_n)max_1≤ j≤ nξ_2^j).





To facilitate our analysis, we now eliminate the historic term H^ℓ(t_n) from
the fast L1  formula (<ref>) for αu^n.
From the recurrence relationship (<ref>), it is easy to see that

    H^ℓ(t_k)=∑_j=1^ke^-θ^ℓ(t_k-t_j)b^(j,ℓ)∇_τu^j,   k≥1,   1≤ℓ≤ N_q.

Inserting this in (<ref>) and using
the definition (<ref>), one obtains the alternative formula

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

where the discrete convolution kernel A_n-k^(n) is henceforth 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.


The formula (<ref>) takes the form of
(<ref>), and we now verify that our A_n-k^(n)
defined by (<ref>) satisfy  1 and 2,
allowing us to apply Theorem <ref> and establish the convergence of our computed solution.
Part (I) of the next lemma ensures that 1 is valid,
while part (II) implies that 2 holds true with π_A=3/2.



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,  1≤ k≤ n-1;  (II) A_0^(n)=a_0^(n) and A_n-k^(n)≥2/3a_n-k^(n),   1≤ k≤ n-1.



The definition (<ref>) and Lemma <ref> (i) yield

    a_0^(n)-a_1^(n)>a_0^(n)-ω_1-α(τ_n)
    =ατ_nω_2-α(τ_n)≥α ω_2-α(1)≥ϵ ,

where the step size τ_n≤ 1 and our hypothesis on ϵ are used.
The definition (<ref>) and Lemma <ref> imply that A_0^(n)=a_0^(n)>a_1^(n)+ϵ>A_1^(n).
Lemma <ref> also shows that θ^ℓ, ϖ^ℓ>0 for ℓ =1,…, N_q;
the mean-value theorem now yields property (I). By Lemma <ref> (i)
and our hypothesis on ϵ we have
ϵ≤1/3ω_1-α(t_n - t_k-1)<1/3a_n-k^(n) for 1≤ k≤ n-1.

Hence Lemma <ref> gives A_n-k^(n)≥ a_n-k^(n)-ϵ≥2/3a_n-k^(n) for
1≤ k≤ n-1. The proof is complete.
















































§ GLOBAL CONSISTENCY ERROR ANALYSIS







We now proceed with the consistency error analysis of our fast linearized method, and begin with the
consistency error of the standard L1 formula  (αu)^n of (<ref>).


For v∈ C^2(0,T] with ∫_0^T t  |v”(t)| s < ∞, one has

    (αv)(t_n)-(αv)^n≤
    a_0^(n)G^n+∑_k=1^n-1a_n-k-1^(n)-a_n-k^(n)G^k,   n≥1,

where the L1 kernel a_n-k^(n) is defined by (<ref>) and G^k :=2∫_t_k-1^t_k⟨t-t_k-1|v”(t) t.



From Taylor's formula with integral remainder, the truncation error of
the standard L1 formula at time t=t_n is  (see <cit.>)

    (αv)(t_n)-(αv)^n
    =    ∑_k=1^n∫_t_k-1^t_kω_1-α(t_n-s)⟨v'(s)-∇_τ v^k/τ_k| s
    
    =    ∑_k=1^n∫_t_k-1^t_kv”(t) Π_1,kQ(t)t ,   n≥1,

where Q(t)=ω_2-α(t_n-t) and we use the notation of the proof of Lemma <ref>.










By the error formula for linear interpolation <cit.>, we have

    1,kQ(t)=∫_t_k-1^t_kχ_k(t,y)Q”(y)y,   t_k-1<t<t_k, 1≤ k≤ n,

where the Peano kernel χ_k(t,y)=max{t-y,0}-t-t_k-1/τ_k(t_k-y) satisfies

    -t-t_k-1τ_k(t_k-t)≤χ_k(t,y)<0    for any t,y∈(t_k-1,t_k).

Observing that for each fixed n≥1 the function Q is decreasing and Q”(t)=ω_-α(t_n-t)<0,
we arrive at
the interpolation error 1,kQ(t)≥0 for 1≤ k≤ n, with

    1,nQ(t)
       ≤ Q(t_n-1)-Π_1,nQ(t)=(t-t_n-1)a_0^(n),
    1,kQ(t)    ≤
    (t_k-1-t)∫_t_k-1^t_kQ”(t)t≤(t-t_k-1)ω_1-α(t_n-t_k)-ω_1-α(t_n-t_k-1)
       ≤ 2(t-t_k-1)a_n-k-1^(n)-a_n-k^(n),     t∈(t_k-1,t_k),   1≤ k≤ n-1,

where Lemma <ref> (ii) is used in the last inequality.
Thus, (<ref>) yields

    (αv)(t_n)-(αv)^n≤∫_t_n-1^t_nv”(t)Π_1,nQ(t)t
    +∑_k=1^n-1∫_t_k-1^t_kv”(t)Π_1,kQ(t)t
        ≤ a_0^(n)∫_t_n-1^t_n(t-t_n-1)v”(t)t
    +2∑_k=1^n-1a_n-k-1^(n)-a_n-k^(n)∫_t_k-1^t_k(t-t_k-1)v”(t)t,

and the desired result follows from the definition of G^k.



Compared with the previous estimate in <cit.>,
Lemma <ref> removes the time-step ratios restriction
ρ_k≤1, which is an undesirable limitation on the mesh for the problems
that allow the rapid growth of  the solution at the time far away from t=0.


We now focus on the fast L1 method by taking the initial singularity into account.
Here and hereafter, we denote T̂=max{1,T} and t̂_n=max{1,t_n} for 1≤ n≤ N.



Assume that v∈ C^2((0,T]) and that there exists a constant C_v>0 such that

    v'(t)≤ C_v (1+t^σ-1),  v”(t)≤ C_v (1+t^σ-2),    0<t≤ T,

where σ∈(0,1)∪(1,2) is a parameter. Let Υ^j :=αv(t_j)- αv^j
denote the local consistency error of the fast L1 formula (<ref>).
Assume that the SOE tolerance error ϵ satisfies
ϵ≤1/3min{ω_1-α(T),3α ω_2-α(1)}.
Then the global consistency error

    ∑_j=1^nP_n-j^(n)Υ^j
    
    ≤    C_vτ_1^σ/σ+
    1/1-αmax_2≤ k≤ n(t_k-t_1)^αt_k-1^σ-2τ_k^2-α
    +ϵ/σt_n^αt̂_n-1^ 2

for 1≤ n≤ N. Moreover, if the mesh satisfies G, then

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





The main difference between the fast L1 formula (<ref>)
and the standard L1 formula (<ref>) is that the convolution kernel
is approximated by SOEs with an absolute  tolerance error ϵ.
Thus, comparing the standard L1 formula (<ref>)
with the corresponding fast L1 formula  (<ref>),
by  Lemma <ref> and the regularity assumption (<ref>) one has

    αv^j
    -αv^j   ≤∑_k=1^j-1∇_τv^k/τ_k∫_t_k-1^t_k∑_ℓ=1^N_qϖ^ℓe^-θ^ℓ(t_j-s)
    -ω_1-α(t_j-s) s,
       ≤ϵ∑_k=1^j-1∫_t_k-1^t_kv'(s) s
    ≤ C_v t_j-1 + t_j-1^σ/σϵ≤C_v/σt̂_j-1^ 2ϵ,    j≥1.

Lemma <ref> implies that A_n-k^(n)-a_n-k^(n)≤ϵ for 1≤ k≤ n-1.
Recalling that A_0^(n)=a_0^(n), one has
a_j-k-1^(j)-a_j-k^(j)≤  A_j-k-1^(j)-A_j-k^(j)+2ϵ for 1≤ k≤ j-1.



Then Lemma <ref> and the regularity assumption (<ref>) lead to

    (αv)(t_j)-αv^j
    ≤    
    A_0^(j)G^j+∑_k=1^j-1A_j-k-1^(j)-A_j-k^(j)G^k+2ϵ∑_k=1^j-1G^k
    ≤    
    A_0^(j)G^j+∑_k=1^j-1A_j-k-1^(j)-A_j-k^(j)G^k
    +4ϵ∑_k=1^j-1∫_t_k-1^t_ktv”(t)t
    ≤    
    A_0^(j)G^j+∑_k=1^j-1A_j-k-1^(j)-A_j-k^(j)G^k
    +C_v/σt̂_j-1^ 2ϵ,   j≥1.

Now a triangle inequality gives

    Υ^j
    
    ≤    A_0^(j)G^j+∑_k=1^j-1A_j-k-1^(j)-A_j-k^(j)G^k
    + C_v/σt̂_j-1^ 2ϵ,   j≥1.

Multiplying the above inequality (<ref>) by
P_n-j^(n) and summing the index j from 1 to n, one can
exchange the order of summation and apply
the definition (<ref>) of P_n-j^(n) to obtain

    ∑_j=1^nP_n-j^(n)    Υ^j≤∑_j=1^nP_n-j^(n)A_0^(j)G^j+
    ∑_j=2^nP_n-j^(n)∑_k=1^j-1A_j-k-1^(j)-A_j-k^(j)G^k
    +C_vϵ/σ∑_j=2^nP_n-j^(n)t̂_j-1^ 2
    
    =    ∑_j=1^nG^jP_n-j^(n)A_0^(j)
    +∑_k=1^n-1G^k∑_j=k+1^nP_n-j^(n)A_j-k-1^(j)-A_j-k^(j)
    +C_vt̂_n-1^ 2ϵ/σ∑_j=2^nP_n-j^(n)
    ≤    ∑_k=1^nP_n-k^(n)A_0^(k)G^k+∑_k=1^n-1P_n-k^(n)A_0^(k)G^k
    +C_v/σt_n^αt̂_n-1^ 2ϵ,

where the property (<ref>) with π_A=3/2 is used in the last inequality.


If the SOE approximation error ϵ≤1/3min{ω_1-α(T),3α ω_2-α(1)} ,
Lemma <ref> (II) and
Lemma <ref> (i) imply that
A_0^(k)=a_0^(k)=ω_2-α(τ_k)/τ_k,
A_k-2^(k)≥23a_k-2^(k)≥23ω_1-α(t_k-t_1), and then

    A_0^(k)/A_k-2^(k)≤32(1-α)(t_k-t_1)^ατ_k^-α,   2≤ k≤ n≤ N.

Furthermore, the identical property (<ref>) for the complementary kernel P^(n)_n-j gives

    P_n-1^(n)A_0^(1)≤1    ∑_k=2^n-1P_n-k^(n)A_k-2^(k)≤∑_k=2^nP_n-k^(n)A_k-2^(k)=1.

The regularity assumption (<ref>) gives
G^1≤ C_vτ_1^σ/σ and G^k≤ C_vt_k-1^σ-2τ_k^2 (2≤ k≤ n).
Thus it follows from (<ref>)  that

    ∑_j=1^nP_n-j^(n)Υ^j≤    2G^1+2∑_k=2^nP_n-k^(n)A_0^(k)G^k
    +C_v/σt_n^αt̂_n-1^ 2ϵ
    ≤    
    C_vτ_1^σ/σ+C_v/1-α∑_k=2^nP_n-k^(n)A_k-2^(k)
    (t_k-t_1)^αt_k-1^σ-2τ_k^2-α
    +C_v/σt_n^αt̂_n-1^ 2ϵ
    ≤    
    C_vτ_1^σ/σ
    +1/1-αmax_2≤ k≤ n(t_k-t_1)^αt_k-1^σ-2τ_k^2-α
    +1/σt_n^αt̂_n-1^ 2ϵ,   1≤ n≤ N.

The claimed estimate (<ref>) is verified.
In particular, if G holds, one has

    t_k^αt_k-1^σ-2τ_k^2-α≤     C_γt_k^σ-2+ατ_k^2-α-βτ^βmin{1,t_k^β-β/γ}
    ≤     C_γt_k^σ-β/γτ_k/t_k^2-α-βτ^β≤ C_γt_k^max{0,σ-(2-α)γ}τ^β,   2≤ k≤ N,

where β=min{2-α,γσ}. The final estimate follows since
τ_1^σ≤ C_γτ^γσ≤ C_γτ^β.









Next lemma describes the global consistency error of Newton's linearized approach, which is
smaller than that generated by the above L1 approximation. In addition,
there is no error in the linearized approximation if f=f(u) is a linear function.

Assume that v∈ C([0,T])∩ C^2((0,T]) satisfies the regularity condition (<ref>),
and the nonlinear function f=f(u)∈ C^2(ℝ).
Denote v^n=v(t_n) and the local truncation error
ℛ_f^n=f(v^n)-f(v^n-1)-f'(v^n-1)∇_τv^n such that the global consistency error

    ∑_j=1^nP_n-j^(n)ℛ_f^j≤ C_vτ_1^α⟨τ_1^2+τ_1^2σ/σ^2|
    +C_vt_n^αmax_2≤ j≤ nτ_j^2+t_j-1^2σ-2τ_j^2,   1≤ n≤ N.

Moreover, if the assumption G holds, one has

    ∑_j=1^nP_n-j^(n)ℛ_f^j≤    C_vτ^min{2,2γσ}max{1,τ^γα/σ^2},   1≤ n≤ N.



Applying the formula of
Taylor expansion with integral remainder, one has

    ℛ_f^j=(∇_τv^j)^2∫_0^1f”v^j-1+s∇_τv^j(1-s)s,   j≥1.

Under the regularity conditions, one has
ℛ_f^1≤ C_v∫_t_0^t_1v'(t)t^2
≤ C_v⟨τ_1^2+τ_1^2σ/σ^2|,

    ℛ_f^j≤ C_v∫_t_j-1^t_jv'(t)t^2
    ≤ C_vτ_j^2+t_j-1^2σ-2τ_j^2,   2≤ j≤ N.

Note that, Lemma <ref> (II) and
the definition (<ref>) give
A_0^(k)=a_0^(k)=ω_2-α(τ_k)/τ_k,
so the identical property (<ref>) shows
P_n-1^(n)≤1/A_0^(1)≤Γ(2-α)τ_1^α.
Moreover, the bounded estimate (<ref>) with π_A=3/2
gives ∑_j=2^nP_n-j^(n)≤3/2ω_1+α(t_n).
Thus, it follows that

    ∑_j=1^nP_n-j^(n)ℛ_f^j≤    P_n-1^(n)ℛ_f^1
    +∑_j=2^nP_n-j^(n)ℛ_f^j≤
    C_vτ_1^αℛ_f^1+C_vt_n^αmax_2≤ j≤ nℛ_f^j
    ≤    
    C_vτ_1^α⟨τ_1^2+τ_1^2σ/σ^2|
    +C_vt_n^αmax_2≤ j≤ nτ_j^2+t_j-1^2σ-2τ_j^2,   1≤ n≤ N.

If G holds, one has
τ_j^2≤ C_γτ^2min{1,t_j^2-2/γ}≤ C_γτ^βmin{1,t_j^2-2/γ}, and

    t_j-1^2σ-2τ_j^2≤     C_γt_j^2σ-2τ_j^2-βτ^βmin{1,t_j^β-β/γ}
    ≤     C_γt_j^2σ-min{2,2γσ}/γτ_k/t_k^2-βτ^β≤ C_γt_k^max{0,2σ-2/γ}τ^β,   2≤ j≤ N,

where β=min{2,2γσ}.
The second estimate follows since
τ_1^2σ≤ C_γτ^2γσ≤ C_γτ^β.





§ UNCONDITIONAL CONVERGENCE




Assume that the time mesh fulfills 3 and G in the error analysis.
We improve the discrete H^2 energy method in
<cit.>
to prove the unconditional convergence of discrete
solution to the two-level linearized scheme (<ref>).
In this section, K_0, τ_0, τ_1, τ_0^*, h_0, ϵ_0 and
any numeric subscripted c, such as c_0, c_1, c_2 and so on,
are fixed values, which are always dependent on the given data and the solution,
but independent of the time-space grid steps and the inductive index k in the mathematical induction as well.
To make our ideas more clearly, four steps to obtain unconditional error estimate are listed in four subsections.



 §.§ STEP 1: construction of coupled discrete system

We introduce a function
wαu-f(u) with the initial-boundary values
w( ,0)Δ u^0( ) for ∈Ω and w( ,t)-f(0) for ∈∂Ω.
The problem (<ref>) can be formulated into

    w=    αu-f(u),
      ∈Ω̅,  0<t≤ T;
    
    w=    Δ u,  ∈Ω,  0≤ t≤ T.

Let w_h^n be the numerical approximation of function W_h^n=w(_h,t_n)
for _h∈Ω̅_h. As done in subsection 2.2,
one has an auxiliary discrete system: to seek {u_h^n, w_h^n} such that

    w_h^n=    αu_h^n-f(u_h^n-1)-f'(u_h^n-1)∇_τu_h^n ,
      _h∈Ω̅_h,   1≤ n≤ N;
    
    w_h^n=    Δ_hu_h^n ,
      _h∈Ω_h,   0≤ n≤ N;
    
    u_h^0=    u^0(_h),  _h∈Ω̅_h ;  
    u_h^n=0,  _h∈∂Ω_h ,1≤ n≤ N.

Obviously, by eliminating the auxiliary function w_h^n in above discrete system,
one directly arrives at the computational scheme (<ref>).
Alternately, the solution properties of two-level linearized method
(<ref>) can be studied via the auxiliary discrete system
(<ref>)-(<ref>).



 §.§ STEP 2: reduction of coupled error system




Let ũ_h^n=U_h^n-u_h^n, w̃_h^n=W_h^n-w_h^n be the solution errors for _h∈Ω̅_h.
We now have an error system with respect to the error function {w̃_h^n} as

    w̃_h^n=    αũ_h^n-𝒩_h^n+ξ_h^n ,
      _h∈Ω̅_h,   1≤ n≤ N;
    w̃_h^n=    Δ_hũ_h^n+η_h^n ,
      _h∈Ω_h,   0≤ n≤ N;
    ũ_h^0=    0,  _h∈Ω̅_h ;  ũ_h^n=0,  _h∈∂Ω_h ,1≤ n≤ N,

where ξ_h^n and η_h^n denote  temporal and spatial truncation errors, respectively, and

    𝒩_h^n    f'(u_h^n-1)∇_τũ_h^n+f(U_h^n-1)-f(u_h^n-1)
    +⟨f'(U_h^n-1)-f'(u_h^n-1)|∇_τU_h^n
    
    =    f'(u_h^n-1)∇_τũ_h^n
    +ũ_h^n-1∫_0^1f's U_h^n-1+(1-s)u_h^n-1s
        +ũ_h^n-1∇_τU_h^n∫_0^1f”s U_h^n-1+(1-s)u_h^n-1s .

Acting the difference operators Δ_h and α on the equations
(<ref>)-(<ref>), respectively, gives
 
    Δ_hw̃_h^n=    αΔ_hũ_h^n-Δ_h𝒩_h^n+Δ_hξ_h^n ,
      _h∈Ω_h,   1≤ n≤ N;
    αw̃_h^n=    αΔ_hũ_h^n+αη_h^n ,
      _h∈Ω_h,   1≤ n≤ N.

By eliminating the term αΔ_hũ_h^n in the above two equations, one gets
 
    αw̃_h^n=    Δ_hw̃_h^n+Δ_h𝒩_h^n+αη_h^n-Δ_hξ_h^n
      _h∈Ω_h,   1≤ n≤ N;
    w̃_h^0=    η_h^0,  _h∈Ω̅_h ;  w̃_h^n=0,  _h∈∂Ω_h ,1≤ n≤ N;

where the initial and boundary conditions are derived from the error system
(<ref>)-(<ref>).




 §.§ STEP 3: continuous analysis of truncation error

According to the first regularity condition in (<ref>), one has

    η^n≤
    c_1h^2,   0≤ n≤ N.

Since the spatial error η_h^n is defined uniformly
at the time t=t_n (there is no temporal error in the equation (<ref>)),
we can define a continuous function η_h(t) for _h=(x_i,y_j)∈Ω_h,

    η_h(t)=    h_1^2/6∫_0^1∂_x^(4)u(x_i-s h_1,y_j,t)
    +∂_x^(4)u(x_i+s h_1,y_j,t)(1-s)^3s
        +h_2^2/6∫_0^1∂_y^(4)u(x_i,y_j-s h_2,t)
    +∂_y^(4)u(x_i,y_j+s h_2,t)(1-s)^3s ,

such that η_h^n=η_h(t_n).
The second condition in (<ref>) implies η'(t)≤ C_uh^2(1+t^σ-1).
Hence, applying the fast L1 formula (<ref>) and
the equality (<ref>), one has

    ∑_j=1^nP_n-j^(n)αη^j≤    ∑_j=1^nP_n-j^(n)∑_k=1^jA_j-k^(j)∇_τη^k
    =∑_k=1^n∇_τη^k≤c_2/σt̂_n^ 2h^2.



Since the time truncation error ξ_h^n in (<ref>) is defined
uniformly with respect to grid point _h∈Ω̅_h,

we can define a continuous function ξ^n( )=ξ_1^n( )+ξ_2^n( ), where
ξ_1^n, ξ_2^n denotes the truncation errors of
fast L1 formula and Newton's linearized approach respectively,







    ξ_1^n=αu(t_n)-αu^n,  ξ_2^n=∇_τu(t_n)^2∫_0^1f”u(t_n-1)
    +s∇_τu(t_n)(1-s)s,

such that ξ_h^n=ξ^n(x_i,y_j) for _h∈Ω̅_h.
By the Taylor expansion formula, one has

    Δ_hξ_1^n_ij=    ∫_0^1∂_xxξ_1^n(x_i-s h_1,y_j)
    +∂_xxξ_1^n(x_i+s h_1,y_j)(1-s)s
        +∫_0^1∂_yyξ_1^n(x_i,y_j-s h_2)
    +∂_yyξ_1^n(x_i,y_j+s h_2)(1-s)s ,   1≤ n≤ N.

Applying Lemma <ref> with the second
and third regularity conditions in (<ref>), we have

    ∑_j=1^nP_n-j^(n)Δ_hξ_1^j≤    C_u/σ(1-α)τ^min{2-α,γσ}
    +C_u/σt_n^αt̂_n-1^ 2ϵ,   1≤ n≤ N.

Similarly, one can write out an integral expression of Δ_hξ_2^n_ij
by using the Taylor expansion.

Assuming f∈ C^4(ℝ) and
taking τ≤τ_1=√(σ) such that
τ^γα≤τ_1^γα=σ,
we apply Lemma <ref> with the second regularity
condition in (<ref>) to find,

    ∑_j=1^nP_n-j^(n)Δ_hξ_2^j≤ C_uτ^min{2,2γσ}max{1,τ^γα/σ^2}≤C_u/στ^min{2,2γσ},   1≤ n≤ N.

Thus, the triangle inequality leads to

    ∑_j=1^nP_n-j^(n)Δ_hξ^j≤    c_3/σ(1-α)τ^min{2-α,γσ}
    +c_4/σt_n^αt̂_n-1^ 2ϵ ,
       1≤ n≤ N.








 §.§ STEP 4: error estimate by mathematical induction




For a positive constant C_0, let ℬ(0,C_0) be a ball in the space
of grid functions on Ω̅_h such that
max{ψ_∞,∇_hψ,Δ_hψ}≤ C_0
for any grid function {ψ_h}∈ℬ(0,C_0).
Always, we need the following result to treat the nonlinear terms but leave the proof to Appendix A.

 Let F∈ C^2(ℝ) and a grid function {ψ_h}∈ℬ(0,C_0).
 Thus there is a constant C_F>0 dependent on C_0 and C_Ω such that,
 Δ_hF(ψ)v≤ C_FΔ_hv for any {v_h}∈𝒱_h.






Under the regularity assumption (<ref>) with U_h^k=u(_h,t_k), we define a constant

    K_0=1/3max_0≤ k≤ N{U^k_∞,∇_hU^k,Δ_hU^k}.

For a smooth function F∈ C^2(ℝ) and any grid function {v_h}∈𝒱_h,
we denote the maximum value of C_F in Lemma <ref> as c_0 such that

    Δ_hF(w)v≤ c_0Δ_hv  for any grid function {w_h}∈ℬ(0,K_0+1).



Let c_5 be the maximum value of C_Ω to verify the embedding inequalities
in (<ref>), and

    c_6=max{1,c_5}E_α3max{1,ρ}(2K_0+3)c_0T^α,  
    c_7=3c_1+2c_2/σT̂^2+3(2K_0+3)c_0c_1T^α.

Also let τ_0^*=1/√(3max{1,ρ}Γ(2-α)(2K_0+3)c_0) , and

    τ_0=√(σ(1-α)/6c_3c_6) ,  
    h_0=1/√(3c_6c_7) ,  ϵ_0=min{σ/6c_4c_6T̂^2T^α,1/3ω_1-α(T),α ω_2-α(1)}.

For the simplicity of presentation, denote 

    E_k E_α3max{1,ρ}(2K_0+3)c_0t_k^α,
       𝒯^k2c_3/σ(1-α)τ^min{2-α,γσ}
    +2c_1+2c_2/σt̂_k^ 2
    +3(2K_0+3)c_0c_1t_k^αh^2+2c_4/σt_k^αt̂_k-1^ 2ϵ ,

where 1≤ k≤ N. We now apply the mathematical induction to prove that

    Δ_hũ^k≤ E_k𝒯^k+c_1h^2
      for 1≤ k≤ N,

if the time-space grids and the SOE approximation satisfies

    τ≤min{τ_0,τ_1,τ_0^*},    h≤ h_0,  ϵ≤ϵ_0.

Note that, the restrictions in (<ref>) ensures
the error function {ũ_h^k}∈ℬ(0,1) for 1≤ k≤ N.




Consider k=1 firstly. Since ũ_h^0=0, {u_h^0}∈ℬ(0,K_0)⊂ℬ(0,K_0+1) and
the nonlinear term (<ref>) gives
𝒩_h^1=f'(u_h^0)ũ_h^1.
For the function f∈ C^3(ℝ), the inequality (<ref>) implies

    Δ_h𝒩^1=Δ_h⟨f'(u^0)ũ^1|≤ c_0Δ_hũ^1≤ c_0w̃^1+c_0c_1h^2 ,

where the equation (<ref>) and the estimate (<ref>) are used.
Taking the inner product of the equation (<ref>) (for n=1) by w̃_h^1, one
gets

    A_0^(1)∇_τw̃^1,w̃^1≤Δ_h𝒩^1,w̃^1
    +αη^1-Δ_hξ^1,w̃^1,

because the zero-valued boundary condition in (<ref>)
leads to Δ_hw̃^1,w̃^1≤0.
With the view of Cauchy-Schwarz inequality and (<ref>), one has
∇_τw̃^1,w̃^1≥w̃^1∇_τw̃^1
and then

    A_0^(1)∇_τw̃^1≤   Δ_h𝒩^1
    +αη^1-Δ_hξ^1≤ c_0w̃^1
    +αη^1-Δ_hξ^1+c_0c_1h^2 .

Setting τ_1≤τ_0^*≤ 1/√(3max{1,}Γ(2-α)c_0), we apply Theorem <ref>
(discrete fractional Grönwall inequality) with ξ_1^1=αη^1-Δ_hξ^1 and ξ_2^1=c_0c_1h^2 to get

    w̃^1≤     E_α3max{1,}c_0t_1^α2η^0
    +2P_0^(1)αη^1-Δ_hξ^1+3c_0c_1ω_1+α(t_1)h^2
    ≤   E_12c_3/σ(1-α)τ^min{2-α,γσ}
    +2c_1h^2+2c_2/σt̂_1^ 2h^2+3c_0c_1ω_1+α(t_1)h^2≤
    E_1𝒯^1,

where the initial condition (<ref>) and the error estimates (<ref>)-(<ref>) are used.
Thus, the equation (<ref>) and the inequality (<ref>) yield
the estimate (<ref>) for k=1,

    Δ_hũ^1≤w̃^1+η^1≤ E_1𝒯^1+c_1h^2 .



Assume that the error estimate (<ref>) holds for 1≤ k≤ n-1 (n≥2).
Thus we apply the embedding inequalities in (<ref>) to get

    max{ũ^k_∞,∇_hũ^k,Δ_hũ^k}≤max{1,c_5}E_k𝒯^k+c_1h^2,  1≤ k≤ n-1.

Under the priori settings in (<ref>), we have the error function
{ũ_h^k}∈ℬ(0,1), the discrete solution {u_h^k}∈ℬ(0,K_0+1)
for 1≤ k≤ n-1, and the continuous solution {U_h^k}∈ℬ(0,K_0)⊂ℬ(0,K_0+1).
Then, for the function f∈ C^4(ℝ),
one applies the inequality (<ref>) to find that

    Δ_hf'(u^n-1)∇_τũ^n≤ c_0Δ_h∇_τũ^n≤ c_0Δ_hũ^n+c_0Δ_hũ^n-1 ,
       Δ_hũ^n-1f's U^n-1+(1-s)u^n-1≤ c_0Δ_hũ^n-1 ,
       Δ_hũ^n-1∇_τU^nf”s U^n-1+(1-s)u^n-1≤ c_0 Δ_hũ^n-1∇_τU^n≤ 2c_0K_0Δ_hũ^n-1,

where 0≤ s≤1. From the expression (<ref>)
of 𝒩^n and the triangle inequality, one has

    Δ_h𝒩^n≤     c_0Δ_hũ^n
    +2(K_0+1)c_0Δ_hũ^n-1
    ≤     c_0w̃^n+2(K_0+1)c_0w̃^n-1
    +(2K_0+3)c_0c_1h^2 ,

where the equation (<ref>) and the estimate (<ref>)
are used.


Now, taking the inner product of (<ref>) by w̃_h^n,
one gets

    αw̃^n,w̃^n≤Δ_h𝒩^n,w̃^n
    +αη^n-Δ_hξ^n,w̃^n,

because the zero-valued boundary condition in (<ref>) leads to
Δ_hw̃^n,w̃^n≤0.
Lemma <ref> (I) says that the kernels A^(n)_n-k are decreasing,
so the Cauchy-Schwarz inequality gives

    αw̃^n,w̃^n≥    A_0^(n)w̃^n^2
    -∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)w̃^kw̃^n
    -A_n-1^(n)w̃^0w̃^n
    
    =    w̃^nA_0^(n)w̃^n
    -∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)w̃^k
    -A_n-1^(n)w̃^0
    
    =    w̃^n∑_k=1^nA_n-k^(n)w̃^k .

Thus with the help of Cauchy-Schwarz inequality and (<ref>),
it follows from (<ref>) that

    ∑_k=1^nA_n-k^(n)    w̃^k≤Δ_h𝒩^n
    +αη^n-Δ_hξ^n
    ≤    
    c_0w̃^n+2(K_0+1)c_0w̃^n-1
    +αη^n-Δ_hξ^n+(2K_0+3)c_0c_1h^2 .

Setting the maximum time-step τ≤τ_0^*=1/√(3max{1,}Γ(2-α)(2K_0+3)c_0),
we apply Theorem <ref> with
ξ_1^n=αη^n-Δ_hξ^n and ξ_2^n=(2K_0+3)c_0c_1h^2  to get

    w̃^n≤     E_n
    ⟨2η^0
    +2max_1≤ j≤ n∑_k=1^jP_j-k^(j)αη^k-Δ_hξ^k+3(2K_0+3)c_0c_1ω_1+α(t_n)h^2|
    ≤    E_n2c_3/σ(1-α)τ^min{2-α,γσ}
    +2c_4/σt_n^αt̂_n-1^ 2ϵ
        +E_n2c_1+2c_2/σt̂_n^ 2+3(2K_0+3)c_0c_1ω_1+α(t_n)h^2≤ E_n𝒯^n,

where the initial data (<ref>) and
the three estimates (<ref>)-(<ref>) are used.
Then the error equation (<ref>) with (<ref>) imply
that the claimed error estimate (<ref>) holds for k=n,

    Δ_hũ^n≤ E_n𝒯^n+c_1h^2 .

The principle of induction and the third inequality in (<ref>)
give the following result.



Assume that the solution of nonlinear subdiffusion problem (<ref>)
with the nonlinear function f∈ C^4(ℝ)
fulfills the regularity assumption (<ref>)
with σ∈(0,1)∪(1,2).
If the SOE approximation error ϵ≤ϵ_0
and the maximum step size τ≤min{τ_0,τ_1,τ_0^*},
the discrete solution of two-level linearized fast scheme
(<ref>),
on the nonuniform time mesh
satisfying 3 and G, is unconditionally convergent,

    U^k-u^k_∞≤C_u/σ(1-α)max{1,ρ}⟨τ^min{2-α,γσ}+h^2+ϵ|,    1≤ k≤ N.

It achieves an optimal time accuracy of order O(τ^2-α)
if γ≥max{1,(2-α)/σ}.







§ NUMERICAL EXPERIMENTS


Two numerical examples are reported here to support our theoretical analysis.
The two-level linearized scheme (<ref>)
runs for solving the fractional Fisher equation

    αu=Δ u+u(1-u)+g( ,t),   (,t)∈(0,π)^2×(0,T],

subject to zero-valued boundary data, with two different initial data and exterior forces:

  
  * (Example 1) u^0( )=sin x sin y and g( ,t)=0
  such that no exact solution is available;
  
  * (Example 2) g( ,t) is specified such that
  u( ,t)=ω_σ(t)sin x sin y, 0<σ<2.

Note that, Example 2 with the regularity parameter σ is set to examine
the sharpness of predicted time accuracy
on nonuniform meshes. Actually, our present theory also fits for the semilinear problem
with nonzero force g( ,t)∈ C(Ω̅×[0,T]).















In our simulations, the spatial domain Ω is divided uniformly into M parts in each direction (M_1 = M_2=M) and
the time interval [0,T] is divided into two parts [0, T_0] and [T_0, T] with total N_T subintervals.
According to the suggestion in <cit.>, the graded mesh t_k=T_0⟨k/N|^γ
is applied in the cell [0, T_0] and the uniform mesh with time step size τ≥τ_N
is used over the remainder interval.
Given certain final time T and a proper number N_T,
here we would take T_0=min{1/γ,T}, N=⌈N_T/T+1-γ^-1⌉ such that
τ=T-T_0/N_T-N≥T+1-γ^-1/N_T≥ N^-1≥τ_N.

Always, the absolute tolerance error of SOE approximation is set to ϵ=10^-12
such that the two-level L1 formula (<ref>) is comparable
with the L1 formula (<ref>) in time accuracy.


In Example 1, we investigate the asymptotic behavior of solution near t=0 and the computational efficiency
of the linearized method (<ref>).
Setting M = 100, T=1/γ and N_T= 100,
Figures <ref>-<ref>
depict, in log-log plot, the numerical behaviors
of first-order difference quotient ∇_τu_h^n/τ_n at
three spatial points near the initial time for different
fractional orders and grading parameters. Observations suggest that
logu_t( ,t)≈ C_u( )+(α-1)log t as t→0, and
the solution is weakly singular near the initial time. Compared with the uniform grid,
the graded mesh always concentrates much more points in the initial
time layer and provides better resolution for the initial singularity.





To see the effectiveness of our linearized method (<ref>),
we also consider another linearized method by replacing the two-level fast L1 formula αu_h^n with
the nonuniform L1 formula  αu_h^n defined in (<ref>).
Setting α = 0.5, γ=2, and M = 50, the two schemes are run for Example 1 to the final time T = 50
with different total numbers N_T. Figure <ref> shows
the CPU time in seconds for both linearized procedures versus the total number N_T of subintervals.
We observe that the proposed method has almost linear complexity in N_T and is
much faster than the direct scheme using traditional L1 formula.

Since the spatial error O(h^2) is standard, the time accuracy due to the numerical
approximations of Caputo derivative and nonlinear reaction is examined in Example 2 with T=1.
The maximum norm error e(N,M)=max_1≤ l≤ NU(t_l)-u^l_∞.
To test the sharpness of our error estimate,
we consider three different scenarios, respectively, in Tables 5.1-5.3:

  
  Table 5.1: σ=2-α and γ=1 with fractional orders α=0.4, 0.6 and 0.8.
  
  Table 5.2: α=0.4 and σ=0.4 with grid parameters γ=1,
  3/4γ_, γ_ and 5/4γ_.
  
  Table 5.3: α=0.4 and σ=0.8 with grid parameters γ=1,
  3/4γ_, γ_ and 5/4γ_.





Tables 5.1 lists the solution errors, for σ=2-α, on the gradually
refined grids with the coarsest grid of N=50. Numerical data indicates that
the optimal time order is of about O(τ^2-α),
which dominates the spatial error O(h^2).
Always, we take M=N in Tables 5.1-5.3 such that e(N,M)≈ e(N).
The experimental rate (listed as Order in tables) of
convergence is estimated by observing that e(N)≈ C_uτ^β  and then
β≈log_2e(N)/e(2N).



       




Numerical results in Tables 5.2-5.3 (with α=0.4 and σ<2-α)
support the predicted time accuracy
in Theorem <ref> on
the smoothly graded mesh t_k=T(k/N)^γ.
In the case of a uniform mesh (γ=1), the solution
is accurate of order O(τ^σ), and the nonuniform meshes
improve the numerical precision and convergence rate of solution evidently.
The optimal time accuracy O(τ^2-α) is observed when the grid parameter γ≥(2-α)/σ.





§ ACKNOWLEDGEMENTS

The authors gratefully thank Professor Martin Stynes for his valuable discussions
and fruitful suggestions during the preparation of this paper.
Hong-lin Liao would also thanks for the hospitality of Beijing CSRC during the period of his visit.



§ PROOF OF LEMMA <REF>

A.equation


Consider F(ψ)=ψ firstly. It is easy to check that, at point _h=(x_i,y_j)∈Ω_h,

    δ_x^2(ψ_ijv_ij)=ψ_ijδ_x^2v_ij+δ_xψ_i-1/2,jδ_xv_i-1/2,j
    +δ_xψ_i+1/2,jδ_xv_i+1/2,j+v_ijδ_x^2ψ_ij ,

so that δ_x^2(ψ v)≤ C_0⟨v+δ_xv+δ_x^2v|.
Similarly, δ_y^2(ψ v)≤ C_0⟨v+δ_yv+δ_y^2v|.



Moreover, one has
δ_yδ_x(ψ v)≤
C_0⟨v+δ_xv+δ_yv+δ_yδ_xv|,
due to the fact

    δ_yδ_x(ψ_i-1/2,j-1/2v_i-1/2,j-1/2)
    =    ψ_i-1/2,j-1/2δ_yδ_xv_i-1/2,j-1/2
    +δ_yψ_i-1/2,j-1/2δ_xv_i-1/2,j-1/2
        +δ_xψ_i-1/2,j-1/2δ_yv_i-1/2,j-1/2
    +δ_yδ_xψ_i-1/2,j-1/2v_i-1/2,j-1/2 .

Noticing that
Δ_hv^2=δ_x^2v^2+2δ_xδ_yv^2+δ_y^2v^2,
we apply the embedding inequalities in (<ref>) to obtain,
also see <cit.>,


    Δ_h(ψ v)≤ C_u⟨v+Δ_hv|≤ C_FΔ_hv,

where the constant C_F is dependent on C_0 and C_Ω.
For the general case F∈ C^2(ℝ), one has

    δ_x^2F(ψ_ij)v_ij=    F(ψ_ij)δ_x^2v_ij
    +δ_xF(ψ_i-1/2,j)δ_xv_i-1/2,j
        +δ_xF(ψ_i+1/2,j)δ_xv_i+1/2,j+v_ijδ_x^2F(ψ_ij) .

The formula of Taylor expansion with integral remainder gives

    δ_xF(ψ_i-1/2,j)=    F(ψ_ij)-F(ψ_i-1,j)/h_1
    =δ_xψ_i-1/2,j∫_0^1F's ψ_ij+(1-s)ψ_i-1,js ,
    δ_x^2F(ψ_ij)=    δ_x^2ψ_ijF'(ψ_ij)
    +δ_xψ_i-1/2,j^2∫_0^1F”s ψ_ij+(1-s)ψ_i-1,j(1-s)s
        +δ_xψ_i+1/2,j^2∫_0^1F”s ψ_ij+(1-s)ψ_i+1,j(1-s)s ,

such that δ_xF(ψ)≤ C_F and δ_x^2F(ψ)≤ C_F.
Therefore, simple calculations arrive at

    δ_x^2F(ψ) v)≤ C_F⟨v+δ_xv+δ_x^2v| .

By presenting similar arguments as those in the above simple case,
it is straightforward to get claimed estimate and complete the proof.






99




BH17 D. Baffet and J.S. Hesthaven,  A kernel compression scheme for fractional differential equations, SIAM J. Numer. Anal., 55(2) (2017), 496-520.

Brunner:book
H. Brunner,
Collocation methods for Volterra integral and related functional differential equations,
15 (2004), Cambridge University Press, Cambridge.

BrunnerLingYamamoto:2010
H. Brunner, L. Ling and M. Yamamoto,
Numerical simulations of
2D fractional subdiffusion problems, J. Comput. Phys., 229 (2010), 6613-6622.

JiangZhangZhangZhang:2017
S. Jiang, J. Zhang, Q. Zhang and Z. Zhang,
Fast evaluation of the Caputo
fractional derivative and its applications to fractional diffusion equations,
 Commun. Comput. Phys., 21(3) (2017), 650-678.

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

Hilfer:2000
R. Hilfer, ed., Applications of fractional calculus in physics,
World Scientific, Singapore, 2000.

LiSun:2013
B. Li and W. Sun,
Error analysis of linearized semi-implicit Galerkin finite elementmethods for nonlinear
parabolic equations, Int. J. Numer. Anal. Model., 10 (2013), 622-633.

LiSun:2013sinum
B. Li and W. Sun,
Unconditional convergence and optimal error estimates of a Galerkin-mixed FEM for
incompressible miscible flow in porous media, SIAM J. Numer. Anal., 51 (2013), 1959-1977.

LiWangSun:2014
B. Li, J. Wang and W. Sun,
The stability and convergence of fully discrete Galerkin FEMs for
porous medium flows, Commun. Comput. Phys., 15 (2014), 1141-1158.

LiYiChen:2016
C. Li, Q. Yi and A. Chen,
Finite difference methods with non-uniform meshes for nonlinear fractional differential
equations, J. Comput. Phys., 316 (2016), 614-631.




Li:2010
J. Li,
A fast time stepping method for evaluating fractional integrals,
SIAM J. Sci. Comput., 31 (2010), 4696-4714.












LiaoSunShi:2010
H.-L. Liao, Z. Z. Sun and H. S. Shi,
Error estimate of fourth-order compact scheme for solving
linear Schrödinger equations,  SIAM. J. Numer. Anal.,
47(6) (2010), 4381-4401.

LiaoSunShi:2010B
H.-L. Liao, Z. Z. Sun and H. S. Shi,
Maximum norm error analysis of explicit schemes for
two-dimensional nonlinear Schrödinger equations (in Chinese),
Science China Mathematics, 40(9) (2010), 827-842.

LiaoSun:2010
H.-L. Liao and Z. Z. Sun,
Maximum norm error bounds of ADI and compact ADI methods for
solving parabolic equations, Numer. Methods PDEs, 26 (2010), 37-60.



















LiaoLiZhangZhao:2017
H.-L. Liao, D. Li and J. Zhang,
Sharp error estimate of nonuniform L1 formula for linear
reaction-subdiffusion equations,  To appear in SIAM J. Numer. Anal., 2018.

LiaoMcLeanZhangZhao:2018
H.-L. Liao, W. McLean and J. Zhang,
A discrete Grönwall inequality
with application to numerical schemes for reaction-subdiffusion problems, 2018,
Preprint, DOI:10.13140/RG.2.2.34292.86408.

MustaphaMustapha:2010
K. Mustapha and H. Mustapha, A second-order accurate numerical method for a semilinear
integro-differential equation with a weakly singular kernel,
IMA J. Numer. Anal., 30(2) (2010), 555-578.

Podlubny:1999
I. Podlubny, Fractional differential equations, Academic Press, New
York, 1999.

StynesRiordanGracia:2017
M. Stynes, E. O'Riordan and J. L. Gracia,
Error analysis of a finite difference method
on graded meshes for a time-fractional diffusion equation,
SIAM J. Numer. Anal., 55(2) (2017), 1057-1079.







XuHesthavenChen:2015
Q. Xu, J.S. Hesthaven and F. Chen,
A parareal method for time-fractional differential equations,
J. Comp. Phys., 293(C) (2015), 173-183.

YanSunZhang:2017
Y. Yan, Z. Z. Sun and J. Zhang,
Fast evaluation of the Caputo fractional derivative and its applications to
fractional diffusion equations: a second-order scheme, Commun. Comput. Phys.,
22 (2017), 1028-1048.






ZhangSunLiao:2014
Y.N. Zhang, Z.Z. Sun and H.-L. Liao, Finite difference methods
for the time fractional diffusion equation on nonuniform
meshes, J. Comput. Phys., 265 (2014), 195-210.








