A second-order scheme with nonuniform time steps for a linear reaction-sudiffusion problem
Hong-lin Liao William McLean Jiwei Zhang

Diversity Regularized Spatiotemporal Attention
for Video-based Person Re-identification
    
Shuang Li^1Work done while employed at Disney Research Pitsburgh.Slawomir Bak^2[1] Peter Carr^2[1] Xiaogang Wang^1 
     ^1CUHK-SenseTime Joint Lab,
 The Chinese University of Hong Kong     ^2Argo AI {sli,xgwang}@ee.cuhk.edu.hk     {sbak,pcarr}@argo.ai


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


It is reasonable to assume that a discrete convolution structure dominates
the local truncation error of any numerical Caputo formula because the
fractional time derivative and its discrete approximation have the same
convolutional form. We suggest an error convolution structure (ECS) analysis
for a class of interpolation-type approximations to the Caputo fractional
derivative. Our assumptions permit the use of adaptive time steps, such as is
appropriate for accurately resolving the initial singularity of the solution
and also certain complex behavior away from the initial time. The ECS analysis
of numerical approximations has two advantages: (i) to localize (and
simplify) the analysis of the approximation error of a discrete convolution
formula on general nonuniform time grids; and (ii) to reveal the error
distribution information in the long-time integration via the global
consistency error. The core result in this paper is an ECS bound and a global
consistency analysis of the nonuniform Alikhanov approximation, which is
constructed at an offset point by using linear and quadratic polynomial
interpolation.  Using this result, we derive a sharp L^2-norm error estimate
of a second-order Crank-Nicolson-like scheme for linear reaction-subdiffusion
problems.  An example is presented to show the sharpness of our analysis.

















Caputo fractional derivative, nonuniform time mesh,
error convolution structure, global consistency error,
stability and convergence



65M06, 35B65




§ INTRODUCTION

The time-fractional diffusion equation provides a valuable tool for modeling
complex systems such as glassy and disordered media <cit.>. This
paper builds on our recent
results <cit.>
for the nonuniform mesh technique applied to the time discretization of the
following reaction-subdiffusion problem in a bounded domain Ω⊂^d
(d=1, 2, 3),

    αu- u   =κ u+f(,t)      for ∈Ω and 0<t<T,
    
    u   =u_0()      for ∈Ω when t=0,

subject to the homogeneous Dirichlet boundary condition u=0
on ∂Ω. Here, the reaction coefficient κ is a real constant,
and α=_ 0^C𝒟_t^α denotes the Caputo
fractional derivative of order α (0<α<1) with respect to time t,
that is,

    (αv)(t)∫_0^tω_1-α(t-s)v'(s)s  for t>0,  where  ω_β(t) t^β-1/Γ(β).




 §.§ Initial singularity and the nonuniform time meshes technique

In developing numerical methods for solving the subdiffusion
problem (<ref>), an important issue to be considered is that the
solution u is typically less regular than in the case of a classical
parabolic PDE (as the limiting case α→1). Sakamoto and
Yamamoto <cit.> showed that if the initial
data u^0∈ H^2(Ω)∩ H_0^1(Ω), then the unique solution
u∈ C⟨[0,T];H^2(Ω)∩ H_0^1(Ω)|, with
αu∈ C⟨[0,T];L^2(Ω)| and ∂_t u ∈ L^2(Ω).
However, ∂_t u(t)_L^2(Ω)≤ C_ut^α-1
for 0<t≤ T, where the constant C_u>0 is independent of t but may depend
on T. In fact, u can only be a smooth function of t if the initial data
and source term satisfy some restrictive compatibility
conditions <cit.>.

The focus of this paper is on a second-order time discretization
of (<ref>).  The spatial discretization is of less interest: we apply
the standard Galerkin finite element method based on the weak form of the
fractional PDE,

    αu,v+∇ u,∇ v
    	=κu,v+f(t),v  for all v∈ H^1_0(Ω) and for 0<t≤ T,

where u,v denotes the usual inner product in L_2(Ω).
Thus, we construct the usual space of continuous, piecewise-linear functions
with respect to a partition of Ω into subintervals (in 1D), triangles (in
2D) or tetrahedrons (in 3D) with the maximum diameter h, and let X_h denote
the subspace of functions satisfying the homogeneous Dirichlet boundary
condition. In the usual way, the (semidiscrete) Galerkin finite element
solution u_h:[0,T]→ X_h is then defined by requiring that

    αu_h,χ+∇ u_h,∇χ
    	=κu_h,χ+f(t),χ  for all χ∈ X_h and for 0<t≤ T,

with u_h(0)=u_0h≈ u_0 for a suitable u_0h∈ X_h.

Consider (generally nonuniform) time levels 0=t_0<t_1<t_2<⋯<t_N=T
and define a fractional time
level t_n-θθ t_n-1+(1-θ)t_n for an off-set parameter
θ∈[0,1/2). We denote the kth time-step size
by τ_k t_k-t_k-1 for 1≤ k≤ N, and the maximum step size by
max_1≤ k≤ Nτ_k. We also define the local step-size
ratios

    ρ_kτ_k/τ_k+1  for 1≤ k≤ N-1,  and put  ρmax_1≤ k≤ N-1ρ_k.


For any time sequence v^k_k=0^N, define the backward
difference v^k:=v^k-v^k-1 and the interpolated
value v^n-θθ v^n-1+(1-θ)v^n. We consider a numerical
Caputo formula approximating (αv)(t_n-θ) of the form

    (αv)^n-θ:=∑_k=1^nA_n-k^(n) v^k≈∑_k=1^n∫_t_k-1^min{t_k,t_n-θ}ω_1-α(t_n-θ-s)v'(s) s,

for appropriate discrete convolution kernels A_n-k^(n). Our fully-discrete
numerical solution, u^n_h()≈ u(,t_n)
for ∈Ω, is then defined by a time-stepping scheme: we require
that u^n_h∈ X_h satisfies

    (αu_h)^n-θ,χ
    	+∇ u_h^n-θ,∇χ
            =κu_h^n-θ,χ+f(t_n-θ),χ

for all χ∈ X_h and for 1≤ n≤ N, with u^0_h=u_0h.

In the literature, several high-order numerical Caputo formulas have a discrete
convolution form like (<ref>), such as the L1-2 schemes
<cit.> and the
L2-1_σ formula <cit.> that applied
the piecewise quadratic polynomial interpolation. They achieve second-order
temporal accuracy for sufficiently smooth solutions when applied to time
approximation of the pure subdiffusion equation (<ref>) with
κ=0. This article considers the L2-1_σ formula of
Alikhanov <cit.>, which employs a quadratic interpolant
in each subinterval [t_k-1,t_k] for 1≤ k≤ n-1, and a linear
interpolant in the final subinterval [t_n-1,t_n-θ]. The offset
parameter is chosen as θ=α/2 (in our notation). As described below,
in the limit as α→1, this scheme reduces to the well-known
Crank–Nicolson method (θ→1/2) for the classical diffusion equation.  We therefore refer to the time-stepping
scheme (<ref>) as a fractional Crank–Nicolson method.

In the special case of uniform time steps τ_n=τ,
the discrete kernels A^(n)_n-k=A_n-k depend only on the
difference n-k, and were shown to be positive and monotonically decreasing,
leading to a proof that the resulting fractional Crank–Nicolson scheme is
stable and convergent of order O(τ^2+h^2) in the L_2-norm assuming that
the solution u is sufficiently smooth <cit.>. However, as
remarked above, in practice the time derivative ∂_tu typically behaves
like O(t^α-1) as t→0 <cit.>, and so
this error bound breaks down.

In resolving a fixed singularity at t=0 of the type described above, a simple
but useful technique to recover an optimal convergence order is to employ a
smoothly graded mesh t_k=T(k/N)^γ, where the grading
parameter γ≥1 is adapted to the strength of the singularity. The larger
the value of γ the more strongly the mesh points are concentrated
near t=0. Actually, such meshes have long been used in the numerical solution
of Fredholm <cit.> and Volterra <cit.> integral
equations, and their use for time-fractional PDEs is now well established
<cit.>.
By using such a nonuniform mesh we will restore the second-order convergence in
time of the fractional Crank–Nicolson scheme in <cit.> when the
solution is not smooth at t=0. This idea was tested recently
in <cit.> to resolve the initial singularity for the
subdiffusion problem, corresponding to κ=0 in (<ref>).
However, this is only a part of our story.

We will establish the stability and convergence theory for the fractional
Crank–Nicolson scheme on a wider class of nonuniform time meshes, not just the
standard graded mesh described above.  In this way, the theory could be
applied in advanced studies on adaptive time grids required to resolve certain
complex behavior (such as physical oscillations, blowup and so on) in nonlinear
time-fractional PDEs. These goals are natural, at least for linear
reaction-subdiffusion equations, since the backward Euler and Crank–Nicolson
schemes for the linear parabolic equation are stable and convergent (provided
τ→0) on arbitrary nonuniform grids with ρ=O(1).

We refer the reader to other high-order time approximations in
<cit.> and the recent survey
paper <cit.>, which describes some useful approaches other than the
nonuniform grids technique to achieve second-order accuracy in time.



 §.§ Error convolution structure (ECS) analysis and a new problem

Generally, our goals are theoretically challenging because the numerical Caputo formula
always has a form of discrete convolutional summation (<ref>).

Actually, the consistency analysis over the whole time interval [t_0,t_n-θ] becomes too cumbersome to implement in practice when there has not enough grid information.
To evade this difficulty, we propose an
error convolution structure (ECS) analysis which begins by recasting
the discrete Caputo formula (<ref>) as

    (αv)^n-θ=A_0^(n)v^n
    -∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)v^k-A_n-1^(n)v^0 .

Consider the local truncation error
Υ^n-θ:=(αv)(t_n-θ)-(αv)^n-θ.
Given the construction of (αv)^n-θ via local interpolation
of v, and provided the discrete convolution kernels A_k^(n) are
decreasing, it is reasonable to conjecture that a discrete convolution structure
dominates the local truncation error:

    (ECS hypothesis)    |Υ^n-θ|
        ≤ A_0^(n)^n
    +∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)^k.

Here, ^n arises from the interpolation error on the local
subinterval [t_n-1,t_n-θ] whereas the ^k (1≤ k≤ n-1)
arise from the interpolation errors over the history [t_0,t_n-1].
Obviously, this ECS hypothesis localizes the consistency analysis of discrete Caputo formulas,
and makes it possible to analyze the numerical approximations on a general class of nonuniform time grids.

Always, there is a loss of accuracy for Υ^n-θ due to the initial
singularity of solution. Actually, Υ^1-θ=O(1) holds on any
mesh and a superconvergence analysis should be required. For example, Stynes et
al. <cit.> showed that, on a graded mesh, the truncation
error of the well-known L1 formula (θ=0) behaves like
Υ^n=On^-min{2-α,γα}. Building on the
ideas first introduced by Liao et al. <cit.>, we
will prove a sharp error estimate via a fractional discrete Gronwall inequality
(Theorem <ref>) that provides a global consistency error in
the form

    ℰ_glob^n:=∑_k=1^nP^(n)_n-k|Υ^k-θ|,
      for 1≤ k≤ n≤ N,

where the complementary discrete convolution kernels P^(n)_n-k
are chosen to enforce the identity

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

In fact, rearranging this identity yields a recursive formula (in effect, a
definition)

    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.

In our recent paper <cit.>, we showed that this approach
is not limited to the L1 rule, but applies to a general class of discrete
convolution kernels A^(n)_n-k satisfying the following two assumptions:


  1. There is a constant π_A>0 such that

    A^(n)_n-k≥1/π_Aτ_k∫_t_k-1^t_kω_1-α(t_n-s)s  for 1≤ k≤ n≤ N;


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

In this case, the complementary kernels P^(n)_n-k in
(<ref>) are well-defined and non-negative, and
satisfy <cit.>

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

From the ECS hypothesis, one can exchange the order of summation
to find

    ℰ_glob^n
           ≤∑_k=1^nP^(n)_n-kA_0^(k)^k
        +∑_k=1^nP^(n)_n-k∑_j=1^k-1A_k-j-1^(k)-A_k-j^(k)^j
       ≤∑_k=1^nP^(n)_n-kA_0^(k)^k
        +∑_j=1^n-1^j∑_k=j+1^nA_k-j-1^(k)-A_k-j^(k)P_n-k^(n),

and then, by using the definition (<ref>)
directly, arrive at

    ℰ_glob^n≤    ∑_k=1^nP^(n)_n-kA_0^(k)^k
    +∑_k=1^n-1P^(n)_n-kA_0^(k)^k .

Thus, by using the properties in (<ref>) and
(<ref>), it is possible to obtain some useful error estimates on
a variety of nonuniform grids, not limited to t_k=T(k/N)^γ.

Obviously, the first term on the right-hand side
of (<ref>) represents the total error contributions
from discretization errors over the n current cells [t_k,t_k-θ]
(1≤ k≤ n), whereas the second term represents those from discretization
errors over the 1/2n(n-1) small cells in the historic intervals
[t_0,t_k-1] (2≤ k≤ n). This observation is very interesting: the
local error in the current cell [t_n,t_n-θ] and the historic errors
in the (long-time) interval [t_0,t_n-1] make almost the same contribution,
in the sense of convolutional summation, to the global consistency error of the
discrete Caputo derivative. If some appropriate time grid is chosen to make ^k=O⟨^k|
according to the error equidistribution  principle,
the error bound (<ref>) becomes

    ℰ_glob^n≲∑_k=1^nP^(n)_n-kA_0^(k)^k .

It suggests that the global approximation error of the numerical Caputo
formula (<ref>) depends mainly on the local error ^k. In this sense,
the error of numerical Caputo
formula  is “local" despite its overtly nonlocal nature.



The ECS hypothesis plays a key role in our analysis. Actually, it has been used implicitly for the nonuniform L1 (fractional backward
Euler-type) method employing a linear interpolant
in each subinterval [t_k-1,t_k] for 1≤ k≤ n.
That analysis <cit.>
showed that the ECS hypothesis
is valid for ρ=1, or in other words provided τ_k≤τ_k+1 for
all k. In a further study <cit.> on the two-level fast L1
scheme (which includes the original L1 scheme as a special case by setting the
SOE approximation error ϵ≡0), the ECS hypothesis is
shown to be valid for any nonuniform mesh with ρ=O(1)
<cit.>, that is,

    |Υ^n|≤
    a_0^(n)G^n+∑_k=1^n-1a_n-k-1^(n)-a_n-k^(n)G^k

where the L1 kernels a_n-k^(n)1/τ_k∫_t_k-1^t_kω_1-α(t_n-s)s for 1≤ k≤ n,
and

    G^k :=2∫_t_k-1^t_k⟨t-t_k-1|v”(t) t  for 1≤ k≤ n.

This local step ratio restriction is the same as
that for the backward Euler scheme for a classical diffusion equation.
Considering Table <ref>, it is then natural to ask
an elementary problem: what restriction on ρ will suffice to ensure that the  fractional
Crank–Nicolson time-stepping scheme (<ref>) is stable and
convergent? We address this problem in the condition M1 below.



 §.§ The discrete fractional Grönwall inequality and our answer


Our answer relies also on a discrete fractional Grönwall inequality suited
to general nonuniform time meshes, proved in our recent
paper <cit.> and
stated below (in a simplified form).  This result involves the aforementioned
complementary discrete convolution kernels P^(n)_n-k, which are
well-defined thanks to our assumptions 1–2 on the discrete
convolution kernels A^(n)_n-k in the numerical Caputo
formula (<ref>).  The Mittag–Leffler function
E_α(z)∑_k=0^∞ z^k/Γ(1+kα) also appears.

Let the criteria 1–2 hold, and the offset parameter θ∈[0,1).
Suppose that λ>0 is a constant independent of the time steps
and that the maximum time-step size

    ≤1/√(2Γ(2-α)π_Aλ).

If the non-negative time sequences (ξ^k)_k=1^N and (v^k)_k=0^N satisfy

    ∑_k=1^nA^(n)_n-k⟨v^k|^2≤λ⟨v^n-θ|^2+v^n-θξ^n  for 1≤ n≤ N,

then the solution v^k_k=0^N satisfies, for 1≤ n≤ N,

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



Thus, we need to complete the following three tasks:


  Task 1. Verify the assumptions 1–2 for the nonuniform
Alikhanov kernels A^(n)_n-k (defined in <ref>)
so that we can use the complementary kernels P^(n)_n-k and apply the
fractional Grönwall inequality to establish the stability of the fully
discrete scheme (<ref>).

  Task 2. Verify the ECS hypothesis on nonuniform time meshes
and determine the corresponding expressions for ^k and ^k to
insert in the bound (<ref>) for the global
consistency error ℰ_glob^n.

  Task 3. Establish a sharp error estimate in L^2 for the fully discrete
scheme (<ref>) for the subdiffusion problem (<ref>)
taking the initial singularity into account.


In more detail, we complete Task 1 in <ref>. We describe the
fractional Crank–Nicolson scheme and the corresponding discrete Alikhanov
kernels A^(n)_n-k, and show in <ref> (The lengthy and technical proofs for these properties of the discrete
kernels A^(n)_n-k are postponed until <ref>) that the criteria
1–2 hold given the following assumption on the mesh.



  1. The parameter θ=α/2, and the maximum time-step
ratio is ρ=7/4.


The special choice of θ in 1 is needed in any case to achieve
second-order accuracy (see <ref>). At the end
of <ref>, the discrete fractional Grönwall inequality is applied to
establish stability for the time-stepping scheme (<ref>).
Actually, by showing that v^n=u^n_h satisfies (<ref>),
the a priori estimate with respect to initial and external
perturbations in the
forms (<ref>)–(<ref>),
follows.

To verify the ECS hypothesis in Task 2 we make use of a proper lower
bound for A^(n)_n-k-1-A^(n)_n-k, already proved in part (II)
of <ref> to ensure the criterion 2 directly. In the first
part of <ref>, an interpolation error formula for quadratic
polynomials is derived in <ref>. Then we
complete Task 2 in <ref> by showing that the ECS
hypothesis and the bound (<ref>) for the global
consistency error ℰ_glob^n are valid under the
condition 1.


Task 3 is completed in the second part of <ref>.
To make our analysis extendable (such as, for distributed-order
subdiffusion problems), we assume that there is
a constant C_u>0 such that the continuous solution u satisfies

    u^(l)(t)_H^2(Ω)≤ C_u1+t^σ-l    for l = 0,1,2,3, and 0<t≤ T,

where σ∈(0,1)∪(1,2) is a regularity parameter. For
example <cit.>,
the assumption (<ref>) holds with σ=α for the subdiffusion problem
(<ref>) if f(,t)=0 and u_0∈ H^1_0(Ω)∩ H^2(Ω).
To resolve such a solution u efficiently, it is appropriate to choose the time
mesh in such a way that the following
condition <cit.> holds.


  2. There is a constant C_γ>0 such that
τ_k≤ C_γτmin{1,t_k^1-1/γ}
for 1≤ k≤ N, with t_k≤ C_γt_k-1 and
τ_k/t_k ≤  C_γτ_k-1/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
2 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 2 is the graded
mesh t_k=T(k/N)^γ.

When the offset parameter θ=0 and (<ref>) is the nonuniform L1
method, our previous work <cit.> proved the
following error bound for the fully discrete scheme (<ref>),

    u(t_n)-u^n_h≤C_u/σ(1-α)τ^min{γσ,2-α}+C_uh^2,  
    1≤ n≤ N.

In particular, the error is of order O(τ^2-α+h^2) if
γ≥(2-α)/σ. When θ=α/2 and
(<ref>) is the Alikhanov formula, <ref>
establishes an error bound

    u(t_n)-u^n_h≤C_u/σ(1-α)τ^min{γσ,2}+C_uh^2,
           1≤ n≤ N,

which is of order O(τ^2+h^2) if γ≥2/σ.  Thus, in comparison
to the L1 scheme, the Alikhanov formula leads to a higher convergence rate with
respect to τ; however, both methods achieve only order
O(τ^σ+h^2) convergence on a uniform mesh. Numerical experiments in
<ref> confirm that our error bound (<ref>) is
sharp.




§ NUMERICAL CAPUTO FORMULA AND STABILITY



Let Π_1,kv denote the linear interpolant of a function v with
respect to the nodes t_k-1 and t_k, and let Π_2,kv denote the
quadratic interpolant with respect to t_k-1, t_k and t_k+1.
The corresponding interpolation errors are denoted by

    (p,kv)(t) v(t)-⟨Π_p,kv|(t)  for   p∈{1,2}.





Recalling that ρ_k=τ_k/τ_k+1, it is easy to find (for instance, by using
the Newton forms of the interpolating polynomials) that

    ⟨Π_1,kv|'(t)= v^k/τ_k  and  ⟨Π_2,kv|'(t)= v^k/τ_k
    	+2(t-t_k-1/2)/τ_k(τ_k+τ_k+1)⟨ρ_k v^k+1- v^k|.

Throughout this paper, we will always use the notation

    ϖ_n(t):=-ω_2-α(t_n-θ-t)≤0
          for 0≤ t≤ t_n-θ.

If 0≤ t<t_n-θ, then
ϖ_n'(t)=ω_1-α(t_n-θ-t)>0,
ϖ_n”(t)=-ω_-α(t_n-θ-t)>0
and ϖ_n”'(t)=ω_-α-1(t_n-θ-t)>0.


 §.§ Discrete Caputo formula

The nonuniform Alikhanov approximation to the Caputo
derivative (αv)(t_n-θ) is defined by

    (αv)^n-θ   ∫_t_n-1^t_n-θϖ_n'(s)
    		⟨Π_1,nv|'(s)s	+∑_k=1^n-1∫_t_k-1^t_kϖ_n'(s)⟨Π_2,kv|'(s)s
       =a^(n)_0 v^n+∑_k=1^n-1
    	a^(n)_n-k v^k+ρ_k b^(n)_n-k v^k+1
    		-b^(n)_n-k v^k,

where the discrete coefficients
a_n-k^(n) and b_n-k^(n) are defined by

    a^(n)_n-k1/τ_k∫_t_k-1^min{t_k,t_n-θ}ϖ_n'(s)s,   1≤ k≤ n;
       b^(n)_n-k2/τ_k(τ_k+τ_k+1)∫_t_k-1^t_k
    	(s-t_k-1/2)ϖ_n'(s)s,    1≤ k≤ n-1.


When θ=0, the coefficients a^(n)_n-k in (<ref>) are just the discrete convolution kernels
in the L1 formula <cit.>. Notice that if α→1, then ω_2-α(t)→1
whereas ω_1-α(t)→0, uniformly for t in any compact
subinterval of the open half-line (0,∞).  Thus,

    a^(n)_0=ω_2-α((1-θ)τ_n)/τ_n→1/τ_n

whereas a^(n)_n-k→0 and b^(n)_n-k→0 for 1≤ k≤ n-1.
It follows that (αv)^n-θ→v^n/τ_k
and θ=α/2→1/2
so the scheme (<ref>) tends to the
Crank–Nicolson method for a linear
parabolic equation. This is why we also call (<ref>) a fractional Crank–Nicolson time-stepping method.

Rearranging the terms in (<ref>), we obtain
the compact form (<ref>)



where the discrete convolution kernels A_n-k^(n) are defined as follows:
A_0^(1) a_0^(1) if n=1 and, for n≥2,

    A^(n)_n-k
    	a^(n)_0+ρ_n-1b^(n)_1,
    	   for k=n,
    
    	a^(n)_n-k+ρ_k-1b^(n)_n-k+1-b^(n)_n-k,
    	   for 2≤ k≤ n-1,
    
    	a^(n)_n-1-b^(n)_n-1,
    	   for k=1.


Before studying the kernels A^(n)_n-k, we present two alternative formulas
for b_n-k^(n). Recall the integral form of error term for the
trapezoidal rule, which can be derived by the Taylor expansion with the
integral remainder.






Integration by parts yields the following identities.

For any function q∈ C^2([t_k-1,t_k]),

    ∫_t_k-1^t_k(s-t_k-1/2)q'(s)s   =-∫_t_k-1^t_k1,kq(s)s
       =1/2∫_t_k-1^t_k(s-t_k-1)(t_k-s)q”(s)s.


Taking q:=ϖ_n in Lemma <ref>, the definition (<ref>) of b_n-k^(n) gives

    b_n-k^(n)=    -2
    ∫_t_k-1^t_kΠ_1,kϖ_n(s)s/τ_k(τ_k+1+τ_k)
    
    =    ∫_t_k-1^t_k(t_k-s)(s-t_k-1)/τ_k(τ_k+1+τ_k)ϖ_n”(s)s,   1≤ k≤ n-1.


The following theorem gathers some useful properties of the
discrete kernels A^(n)_n-k, but the rigorous proof is left to <ref>.
It should be noted here that this proof is quite different from the previous analysis
<cit.> for the discrete convolution kernels
in high-order numerical Caputo formulas with uniform time-steps.


Let 1 hold and consider the discrete kernels defined
in (<ref>).


  (I) The discrete kernels A^(n)_n-k are bounded,

    A_0^(n)≤24/11τ_n∫_t_n-1^t_nω_1-α(t_n-s) s

and
 
    A_n-k^(n)≥4/11τ_k∫_t_k-1^t_kω_1-α(t_n-s) s,   1≤ k≤ n;


  (II) The discrete kernels A^(n)_n-k are monotone,

    A^(n)_n-k-1-A^(n)_n-k≥
    (1+ρ_k)b^(n)_n-k+1/5τ_k∫_t_k-1^t_k⟨t_k-s|ϖ_n”(s)s,
       1≤ k≤ n-1;


  (III) And the first kernel A^(n)_0 is appropriately larger than the second one,
 
    1-2θ/1-θA^(n)_0-A^(n)_1>0   for n≥2.




The first part (I) implies that the criterion 1 holds with
π_A=11/4, the second part (II) ensures that the criterion 2 is
valid and the third part (III) is used to prove the following corollary. These
results allow us to apply <ref> and establish the stability of
the time-stepping scheme (<ref>). Also, the second part (II)
establishes a stronger estimate used in obtaining an ECS bound for the
error analysis (see <ref>).
 Under the condition 1, the discrete Caputo formula (<ref>)
with the discrete kernels (<ref>) satisfies

    ⟨αv|^n-θ,v^n-θ≥1/2∑_k=1^nA^ (n)_n-kv^k^2  for 1≤ n≤ N.



The inequality is known to hold <cit.>
provided 2 is satisfied and
θ^(n)≥θ for 1≤ n≤ N, where

    θ^(1)=1/2  and  θ^(n)=A^(n)_0-A^(n)_1/2A^(n)_0-A^(n)_1  for n≥2.

Obviously, <ref> (II) ensures that 2 holds,
and the condition 1 leads to θ^(1)≥θ.
From <ref> (III), θ^(n)≥θ holds also for n≥2.




 §.§ Unconditional stability

By taking the χ=u_h^n-θ in (<ref>), one has

    ⟨αu_h|^n-θ,u_h^n-θ≤κ_+ u^n-θ_h^2+f(t_n-θ),u^n-θ_h  for 1≤ n≤ N,

where κ_+:=max{κ,0} and the property
 ∇ u_h^n-θ,∇ u_h^n-θ≥0
was used.  Therefore, applying the above <ref> along with the Cauchy–Schwarz
and triangle inequalities, one gets

    ∑_k=1^nA^(n)_n-ku^k_h^2≤2κ_+(1-θ)u^n_h+θu^n-1_h^2
    
    +2(1-θ)u^n_h+θu^n-1_hf(t_n-θ),  1≤ n≤ N,

which has the form of (<ref>) with


    λ:=2κ_+,   v^k:=u^k  and  ξ^k:=2f^k-θ  for 1≤ k≤ N.

Note that <ref> shows the criteria 1–2 of <ref> are
satisfied with π_A=11/4, and the condition 1 gives ρ=7/4. Therefore, applying <ref>, we see that the
time-stepping method (<ref>) is stable in the following sense.


If 1 holds with the maximum time step
≤1/√(11Γ(2-α)κ_+)
(there is no limit to the maximum time step if κ≤0),
then the solution u^n_h of the time-stepping scheme (<ref>)
is stable, that is,

    u^n_h   ≤2E_α(20κ_+t_n^α)
    u_0h+
    	2max_1≤ k≤ n∑_j=1^kP^(k)_k-jf(t_j-θ)
       ≤2E_α(20κ_+t_n^α)
    u_0h+
    	6Γ(1-α)max_1≤ j≤ n{t_j^αf(t_j-θ)}  for 1≤ n≤ N.








§ GLOBAL CONSISTENCY ERROR AND CONVERGENCE



We now derive a representation for the consistency error of the discrete Caputo
derivative (<ref>) with the discrete kernels in (<ref>).
Fix a function v(t) and decompose the local consistency error into n terms
corresponding to the n subintervals, writing

    Υ^n-θ:=(αv)(t_n-θ)-(αv)^n-θ
        =∑_k=1^nΥ^n-θ_k,	    1≤ n≤ N,

where, recalling the notations ϖ_n(s), (1,kv) and (2,kv)
from <ref>,

    Υ^n-θ_k∫_t_k-1^t_kϖ_n'(s)2,kv'(s)s,  1≤ k≤ n-1≤ N-1,
       Υ^n-θ_n∫_t_n-1^t_n-θϖ_n'(s)1,nv'(s)s,   1≤ n≤ N.


Compared with the traditional technique using direct estimation of the local
error Υ^n-θ,
the stability estimate in <ref> suggests
that one can consider the global consistency error ℰ_glob^n,
defined in (<ref>),
accumulated from t=t_1-θ to t=t_n-θ with the complementary discrete
kernel P_n-j^(n). To exploit this convolution structure,
we will control Υ^n-θ
by an ECS bound in terms of the discrete
kernels A_n-k^(n) defined in (<ref>),
and the following quantities

    ^k   3/2∫_t_k-1^t_k-1/2(s-t_k-1)^2|v”'(s)|s
    	+3τ_k/2∫_t_k-1/2^t_k(t_k-s)|v”'(s)|s,
    ^k   5/2∫_t_k-1^t_k(s-t_k-1)^2|v”'(s)|s
    	+5/2∫_t_k^t_k+1(t_k+1-s)^2|v”'(s)|s,

assuming in what follows that v is such that these integrals exist and are
finite.



 §.§ Global consistency error



For any function v∈ C^3((0,T]), the local consistency
error Υ_n^n-θ in (<ref>) satisfies

    Υ_n^n-θ≤ a_0^(n)G_loc^n≤
    A_0^(n)G_loc^n  for 1≤ n≤ N.



Taylor expansion (with integral remainder) about t_n-1/2 shows that

    v'(s)=v'(t_n-1/2)+v”(t_n-1/2)(s-t_n-1/2)
    	+∫_t_n-1/2^s(s-y)v”'(y)y,

and

    1,nv'(s)=    v”(t_n-1/2)(s-t_n-1/2)
    	+∫_t_n-1/2^s(s-y)v”'(y)y
        -1/2τ_n∫_t_n-1^t_n-1/2(y-t_n-1)^2v”'(y)y
    	-1/2τ_n∫_t_n-1/2^t_n(t_n-y)^2v”'(y)y.

Inserting these four terms in (<ref>) yields the splitting
Υ^n-θ_n=∑_ℓ=1^4Υ^n-θ_n,ℓ. After
integrating by parts, we find that

    Υ^n-θ_n,1
        =(α-2θ) (1-θ)^1-α/2Γ(3-α) 
    	v”(t_n-1/2)τ_n^2-α,

which vanishes for θ=α/2.
For the term Υ_n,2^n-θ,


we split the integration interval [t_n-1,t_n-θ] into two parts:
[t_n-1,t_n-1/2] and [t_n-1/2,t_n-θ].
Since t_n-1/2<t_n-θ<t_n, the second term reads

    Υ^n-θ_n,2=   ∫_t_n-1^t_n-θϖ_n'(s)∫_t_n-1/2^s(s-y)v”'(y)ys
    
    =   ∫_t_n-1^t_n-1/2ϖ_n'(s)∫_s^t_n-1/2(y-s)v”'(y)ys
    	+∫_t_n-1/2^t_n-θϖ_n'(s)∫_t_n-1/2^s(s-y)v”'(y)ys.

Reversing the order of integration, then integrating by parts in the second
term and using ϖ_n(t_n-θ)=0, we have

    Υ^n-θ_n,2   =∫_t_n-1^t_n-1/2v”'(y)∫_t_n-1^y(y-s)ϖ_n'(s)sy
    +∫_t_n-1/2^t_n-θv”'(y)
        ∫_y^t_n-θ(s-y)ϖ_n'(s)sy
       =∫_t_n-1^t_n-1/2v”'(y)∫_t_n-1^y(y-s)ϖ_n'(s)sy
    -∫_t_n-1/2^t_n-θv”'(y)∫_y^t_n-θϖ_n(s)sy.

The inner integrals can be estimated as

    |∫_t_n-1^y(y-s)ϖ_n'(s)s|
    	≤ϖ_n'(t_n-1/2)(y-t_n-1)^2/2  for t_n-1<y<t_n-1/2,
       |∫_y^t_n-θϖ_n(s)s|≤ϖ_n(t_n-1/2)(t_n-θ-y)
    	  for t_n-1/2<y<t_n-θ.













Recalling the definition (<ref>) of a^(n)_0, we see that
ω_2-α(t_n-θ-t_n-1)=τ_na^(n)_0 and then

    ϖ_n(t_n-1/2)   =ω_2-α(t_n-θ-t_n-1/2)
    	≤ω_2-α(t_n-θ-t_n-1)=τ_na^(n)_0,
    ϖ_n'(t_n-1/2)   =ω_1-α(t_n-θ-t_n-1/2)
    	=2/τ_nω_2-α(t_n-θ-t_n-1/2)
        ≤ 2a^(n)_0,

where we used the fact that t_n-θ-t_n-1/2=(1-α)τ_n/2.
Hence, it follows that

    Υ^n-θ_n,2≤ a^(n)_0∫_t_n-1^t_n-1/2
        (y-t_n-1)^2|v”'(y)|y+a^(n)_0τ_n
        ∫_t_n-1/2^t_n-θ(t_n-θ-y)|v”'(y)|y,

and finally

    ∑_ℓ=3^4Υ^n-θ_n,ℓ≤a^(n)_0/2∫_t_n-1^t_n-1/2(y-t_n-1)^2|v”'(y)|y+
    	a^(n)_0/2∫_t_n-1/2^t_n(t_n-y)^2|v”'(y)|y.

Thus the triangle inequality yields |Υ^n-θ_n|≤ a^(n)_0^n
where ^n is defined in (<ref>).
The definition (<ref>) implies a^(n)_0≤ A^(n)_0 and completes the proof.



If we were to choose θ≠α/2, the
term (<ref>) would limit the consistency error
to an order of O(τ_n^2-α), even for smooth solutions.




To estimate the remaining terms in (<ref>), we present an
interpolation error formula for the quadratic polynomial Π_2,kv
employed in the Alikhanov formula (<ref>),
but leave the proof to Appendix <ref>.
This formula is crucial for verifying the ECS hypothesis for the local
consistency error Υ^n-θ.


If v∈ C^3([t_k-1,t_k+1]) and q∈ C^2([t_k-1,t_k]), then

    ∫_t_k-1^t_kq'(t)Π_2,kv'(t)t =    ∫_t_k^t_k+1(t_k+1-s)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k+1
        -∫_t_k-1^t_k(s-t_k-1)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k
        +∫_t_k-1^t_kv”'(s)s∫_t_k-1^sΠ_1,kq(t) t,   1≤ k≤ n-1.
























Assume that the mesh condition 1 holds and
v∈ C^3((0,T]). For the nonuniform Alikhanov formula (<ref>) with
the discrete kernels (<ref>), an ECS dominates
the local consistency error Υ^n-θ in (<ref>), that is,

    Υ^n-θ≤ A_0^(n)^n+
    ∑_k=1^n-1A_n-k-1^(n)-A_n-k^(n)^k  for 1≤ n≤ N,

and consequently the global consistency error satisfies

    ℰ_glob^n≤∑_k=1^nP_n-k^(n)A_0^(k)^k+
    ∑_k=1^n-1P_n-k^(n)A_0^(k)^k  for 1≤ n≤ N,

where ^k and ^k are defined by (<ref>) and
(<ref>), respectively.








According to the arguments in <ref>, it suffices to
verify the first inequality (the ECS bound). The definition 
(<ref>) of Υ^n-θ_k and
<ref> (taking qϖ_n) yield

    Υ^n-θ_k
    =b^(n)_n-k/2∫_t_k-1^t_k(s-t_k-1)^2v”'(s)s
    -ρ_kb^(n)_n-k/2∫_t_k^t_k+1(t_k+1-s)^2v”'(s)s
    
    	+∫_t_k-1^t_kv”'(s)∫_t_k-1^s
    		(1,kϖ_n)(t)ts,   1≤ k≤ n-1,

where the alternative definition (<ref>) of b^(n)_n-k has been used.
Recall the error formula of linear interpolation <cit.>,

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

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

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

The inner integral in the last term of (<ref>) can be bounded by

    |∫_t_k-1^s(1,kϖ_n)(t)t|
    ≤1/2(s-t_k-1)^2∫_t_k-1^t_kt_k-s/τ_k ϖ_n”(s)s,   t_k-1<s<t_k.

By the definition (<ref>) of ^n and the triangle inequality, we obtain from (<ref>) that

    Υ^n-θ_k≤1/5(1+ρ_k) b^(n)_n-k
    +∫_t_k-1^t_kt_k-s/τ_kϖ_n”(s)s^k
    ≤A_n-k-1^(n)-A_n-k^(n)^k,

where <ref> (II) was used in the second inequality.
Then the definition (<ref>) and <ref>
yield the first inequality immediately. The proof is completed.


Traditionally, the global approximation error would be estimated by using the truncation error Υ^n-θ
directly. Once an upper bound of Υ^n-θ is available,
the inequality (<ref>) with m=0 will give the global approximate error

    ℰ_glob^n≤∑_j=1^nP_n-j^(n)ω_1-α(t_j)max_1≤ l≤ nΥ^l-θ/ω_1-α(t_l)≤π_AΓ(1-α)max_1≤ l≤ nt_l^αΥ^l-θ .

Nonetheless, the local and global consistency errors described in Theorem <ref>
present a new understanding of the error contributions generated by the two
different polynomial approximations, respectively, in the local cell [t_n-1,t_n-θ]
and the historical interval [0,t_n-1]
of the fractional Caputo derivative.

Originally, our ECS bound for Υ^n-θ is constructed to
preserve the convolution structure of the Caputo fractional derivative as much
as possible.  A direct estimate of the global consistency
error (<ref>) would lead to the double
sum ∑_k=1^n∑_j=1^k|Υ^k-θ_j|, whereas the ECS
bound leads to a single
sum ∑_k=1^nP_n-k^(n)A_0^(k)^k+^k.  This
simplification assists for proving sharp error bounds even with quite general
nonuniform meshes.






Nonetheless, an explicit bound for the
complementary discrete kernel P^(n)_n-j remains an open problem
until now,  and we will make full use of the
identity (<ref>) and the
upper bound (<ref>) in the subsequent analysis.



Assume that v∈ C^3((0,T]),
and there exists a positive constant C_v such that v”'(t)≤ C_v(1+t^σ-3) for 0<t≤ T,
where σ∈(0,1)∪(1,2) is a regularity parameter.
If the mesh condition  1 holds, for 1≤ n≤ N, then the global
consistency error satisfies

    ℰ_glob^n
    ≤ C_vτ_1^σ/σ+t_1^σ-3τ_2^3
    +1/1-αmax_2≤ k≤ nt_k^αt_k-1^σ-3τ_k^3/τ_k-1^α.









The bounds on the discrete kernel A_n-k^(n) in <ref> (I) yield
the inequalities

    A_0^(k)≤24/11ω_2-α(τ_k)/τ_k,  
    A_k-2^(k)≥4/11ω_1-α(t_k-t_1),

and

    A_0^(k)/A_k-2^(k)<6 ω_2-α(τ_k)/τ_k ω_1-α(t_k-t_1)≤6/1-α(t_k-t_1)^α/τ_k^α,   2≤ k≤ n≤ N.

Furthermore, the identity (<ref>) for the complementary discrete
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.

Applying the definition (<ref>) with the regularity assumption,
it is not difficult to get

    ^1≤ C_vτ_1^σ/σ  and  ^k≤ C_v t_k-1^σ-3τ_k^3   for 2≤ k≤ N.

Similarly, by using the formula (<ref>), one gets

    ^1≤
    C_v(τ_1^σ/σ+t_1^σ-3τ_2^3)  and  ^k≤ C_v ⟨t_k-1^σ-3τ_k^3
    +t_k^σ-3τ_k+1^3|  for 2≤ k≤ N-1.

Then it follows from <ref>  that

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

The first term on the right is bounded
by C_v(τ_1^σ/σ+t_1^σ-3τ_2^3), and the remaining
terms can be bounded by

    6/1-α(∑_k=2^nP_n-k^(n)A_k-2^(k)
        t_k^ατ_k^-α^k
    +∑_k=2^n-1P_n-k^(n)A_k-2^(k)
        t_k^ατ_k^-α^k)
       ≤C_v/1-αmax_2≤ k ≤ n
        t_k^α t_k-1^σ-3τ_k^3-α
        +C_v/1-αmax_2≤ k≤ n-1(t_k^αt_k-1^σ-3τ_k^3-α
        + t_k^α+σ-3τ_k+1^3τ_k^-α) 
       ≤C_v/1-αmax_2≤ k ≤ n
        t_k^α t_k-1^σ-3τ_k^3/τ_k-1^α(1+ρ_k-1^α),

implying the claimed estimate.



The proof of Lemma <ref> and the ECS bound in Theorem <ref> give

    Υ^1-θ≤ A_0^(1)^1≤
    C_vτ_1^σ-α/σ,

implying that Υ^1-θ=O(1) when σ=α, and
if 0<σ<α then the situation becomes worse. The global consistency
analysis seems therefore to be also a superconvergence analysis.


Now we describe the contribution to the global truncation error from the time
weighted terms. The next lemma suggests that the temporal error introduced by
the time weighted approach is smaller than that generated by the Alikhanov
approximation of the Caputo derivative.


Assume that v∈ C^2((0,T]),
and there exists a positive constant C_v such that v”(t)≤ C_v(1+t^σ-2) for 0<t≤ T,
where σ∈(0,1)∪(1,2) is a regularity parameter.
Denote the local truncation error of v^n-θ by


    ℛ^n-θ=v(t_n-θ)-v^n-θ  for 1≤ n≤ N.

If the mesh condition 1 holds, then the global consistency error satisfies

    ∑_j=1^nP_n-j^(n)ℛ^j-θ≤    C_vτ_1^σ+α/σ
    +t_n^αmax_2≤ k≤ nt_k-1^σ-2τ_k^2  for 1≤ n≤ N.









The following integral representation of ℛ^j-θ can be easily
verified, for example using the Taylor formula with integral
remainder <cit.>,

    ℛ^j-θ=-θ∫_t_j-1^t_j-θ(s-t_j-1)v”(s)s
    -(1-θ)∫^t_j_t_j-θ(t_j-s)v”(s)s ,   1≤ j≤ N.

Under the regularity assumption,

one has



    ℛ^1-θ≤ C_v τ_1^σ/σ  and  ℛ^j-θ≤ C_vt_j-1^σ-2τ_j^2,
       2≤ j≤ N.


Note that <ref> (I) implies A_0^(1)≥4/11ω_2-α(τ_1)/τ_1,
and then the identity (<ref>) shows that

    P_n-1^(n)≤1/A_0^(1)≤ 3Γ(2-α)τ_1^α.

Therefore we obtain

    ∑_j=1^nP_n-j^(n)ℛ^j-θ=    P_n-1^(n)ℛ^1-θ
    +∑_j=2^nP_n-j^(n)ℛ^j-θ
    ≤    3Γ(2-α)τ_1^αℛ^1-θ
    +max_2≤ k≤ nℛ^k-θ∑_j=1^nP_n-j^(n)
    ≤    C_vτ_1^σ+α/σ
    +t_n^αmax_2≤ k≤ nt_k-1^σ-2τ_k^2,
         1≤ n≤ N,

where the estimate (<ref>) with π_A=11/4 has been used in the last inequality.





 §.§ Convergence


We now establish the convergence of the numerical solution under the regularity
conditions (<ref>) and the assumptions 1–2.
To deal with the spatial error, we introduce the Ritz
projector R_h:H^1_0(Ω)→ X_h, defined by

    ∇ R_hv,∇χ=∇ v,∇χ  for v∈ H^1_0(Ω) and χ∈ X_h.


Suppose that the solution u of (<ref>) has the regularity
property (<ref>) for the parameter σ∈(0,1)∪(1,2), and
consider the time-stepping method (<ref>) using the nonuniform
Alikhanov formula (<ref>) with the discrete
kernels (<ref>). If 1 holds with the
maximum step size ≤1/√(11Γ(2-α)κ_+), then
the discrete solution u_h^n is convergent with respect to the L_2-norm,

    u(t_n)-u^n_h≤ C_uE_α(20κ_+t_n^α) (
    τ_1^σ/σ
    	+1/1-αmax_2≤ k≤ n
    	t_k^α t_k-1^σ-3 τ_k^3/τ_k-1^α
    +t_n^αmax_2≤ k≤ nt_k-1^σ-2τ_k^2	
    
        +u_0h-R_hu_0 +(t_n+t_n^α+t_n^σ)h^2 )
      for 1≤ n≤ N.

In particular, if 2 also holds  and if we choose u_0h=R_hu_0, then

    u(t_n)-u^n_h≤C_u/σ(1-α) τ^min{γσ,2}+C_uh^2
    	  for 1≤ n≤ N,

where C_u may depend on u and T, but is uniformly bounded with respect to α and σ.


Let e^n_h=u^n_h-R_hu^n∈ X_h where u^n=u(t_n), so that

    u^n_h-u^n≤u^n-R_hu^n+e^n_h.

The usual analysis of the elliptic problem shows that, under the first
regularity assumption in (<ref>),

    u^n-R_hu^n≤ C_Ωh^2u^n_H^2(Ω)≤ C_uh^2,

so it suffices to deal with e^n_h.  We find
<cit.> that

    (αe_h)^n-θ,χ
    	+∇ e_h^n-θ,∇χ
            =κe_h^n-θ,χ+ℛ^n,χ,

for all χ∈ X_h, where

    ℛ^n=(αu)(t_n-θ)-(αR_hu)^n-θ
    	-κ(u(t_n-θ)-R_hu^n-θ)
    	+(u^n-θ-u(t_n-θ)).

Choosing χ=u^n-θ_h yields an inequality
of the form (<ref>) with u^n-θ_h and f(t_n-θ)
replaced by e^n-θ_h and ℛ^n, respectively.
Hence, the argument leading to Theorem <ref> shows that

    e^n_h≤2E_α(20κ_+t_n^α)(e^0_h
    	+2max_1≤ k≤ n∑_j=1^k P^(k)_k-jℛ^j)  for 1≤ n≤ N.

Write ℛ^j=ℛ^j_1+ℛ^j_2
+ℛ^j_3+ℛ^j_4, where

    ℛ^j_1   =
    	(αu)(t_j-θ)-(αu)^j-θ,   ℛ^j_2   =
    	(κ+)(u^j-θ-u(t_j-θ)),
    ℛ^j_3   =(α(u-R_hu))^j-θ,   ℛ^j_4   =κ(R_hu-u)^j-θ.

Applying <ref>,
combined with the regularity assumption (<ref>),
one obtains

    max_1≤ k≤ n∑_j=1^k P^(k)_k-jℛ^j_1+ℛ^j_2≤ C_u(τ_1^σ/σ
        +1/1-αmax_2≤ k≤ nt_k^α
    		t_k-1^σ-3 τ_k^3/τ_k-1^α
        +t_n^αmax_2≤ k≤ nt_k-1^σ-2τ_k^2).

Since

    ℛ^j_3=∑_ℓ=1^jA^(j)_j-ℓ∇_τ(u-R_hu)^ℓ≤∑_ℓ=1^j A^(j)_j-ℓ∫_t_ℓ-1^t_ℓ(u-R_hu)'(t) dt,

the identity (<ref>), the error bound (<ref>)
for the Ritz projection and the regularity assumption (<ref>) give

    max_1≤ k≤ n∑_j=1^k P^(k)_k-jℛ^j_3   ≤max_1≤ k≤ n∑_ℓ=1^k
    		(∑_j=ℓ^k P^(k)_k-j A^(j)_j-ℓ)
    		∫_t_ℓ-1^t_ℓ(u-R_hu)'(t) dt
       ≤ C_uh^2∫_0^t_nu'(t)_H^2(Ω) dt
    	≤ C_u(t_n+t_n^σ)h^2.

Recalling the upper bound (<ref>) and the Ritz projection error (<ref>), we see that

    max_1≤ k≤ n∑_j=1^k P^(k)_k-jℛ^j_4≤ C_Ωh^2max_1≤ k≤ n∑_j=1^k P^(k)_k-ju^j-θ_H^2(Ω)≤ C_u t_n^α h^2,

so the first estimate for u^n_h-u(t_n) follows.
If the mesh assumption 2 holds, then
τ_1≤ C_γτ^γ and, with β:=min{2,γσ},

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

In addition, we have

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

so the claimed second result follows immediately by noting that t_n≤ T.



Replacing f(t_n-θ) with f^n-θ
in (<ref>) would introduce an additional
term f^n-θ-f(t_n-θ) in the
definition (<ref>) of ℛ^n, but would not affect the
final error bound, assuming f has the regularity properties needed to apply
<ref>.
Also, instead of u_0h=R_hu_0 we could choose the
interpolant or the L_2-projection of u_0 and still maintain
second-order accuracy in space.



By an argument similar to that in (<ref>), it is not difficult
to show

    t_k^αt_k-1^σ-3τ_k^3/τ_k-1^α≤ C_γt_k^σ-(3-α)/γτ^3-α,

which means that the Alikhanov formula (αv)^n-θ approximates
(αu)(t_n-θ) to order O(τ^3-α) if
γ≥(3-α)/σ.  However, the term (<ref>)
arising from the difference u(t_n-θ)-u^n-θ
in (<ref>) would still limit the convergence rate for the
overall scheme to order O(τ^2).






§ PROOF OF <REF> (DISCRETE CONVOLUTION KERNELS)




Our aim is to prove the boundedness and monotonicity of the convolution
kernels A_n-k^(n). Since the coefficients a_n-k^(n),
b_n-k^(n) and A_n-k^(n) in (<ref>), (<ref>) and
(<ref>) are defined on nonuniform meshes, it is a technically
challenging task and some new techniques will be necessary.

Note that, some techniques using Taylor expansion or function monotonicity
have been applied in <cit.> to investigate
the discrete convolution kernels in high-order numerical Caputo formulas. These techniques
can not be directly applied here although they would be well suited for the uniform case with τ_k=τ.

We start from the alternative definition (<ref>) of b_n-k^(n).
Compared with the original definition (<ref>),
the new formula (<ref>) has a nonnegative integrand
because Γ(-α)<0 and ω_-α(t_n-θ-s)<0,
motivating us to consider the integrals

    ∫_t_k-1^t_kϖ_n”(t)t,  ∫_t_k-1^t_kt_k-t/τ_kϖ_n”(t)t  and  ∫_t_k-1^t_kt-t_k-1/τ_kϖ_n”(t)t .

Actually, they are very close to the values of b_n-k^(n) and a_n-k-1^(n)-a_n-k^(n)
if ρ=O(1). In studying theoretical properties of the discrete kernels
A^(n)_n-k, these integrals will play a bridging role in establishing some
useful links between the underlaying discrete coefficients a_n-k^(n),
b_n-k^(n) and A_n-k^(n); see
<ref>.



 §.§ Proof of <ref> (I)



The discrete coefficients a^(n)_n-k defined in (<ref>) satisfy


  (i)
a^(n)_n-k>ω_1-α(t_n-θ-t_k-1)>a^(n)_n-k+1 for 1≤ k≤ n;

  (ii)
  a^(n)_0>3/4∫_t_n-1^t_nω_1-α(t_n-s)/τ_ns and
a^(n)_n-k>∫_t_k-1^t_kω_1-α(t_n-s)/τ_ks for 1≤ k≤ n-1.


(i) If k=n, one has

    a^(n)_0=1-θ/1-αω_1-α(t_n-θ-t_n-1)>ω_1-α(t_n-θ-t_n-1).

For 1≤ k<n, the claimed inequalities follow directly from the integral mean value theorem and
the fact that ϖ_n'(s)=ω_1-α(t_n-θ-s) is a strictly increasing function.

(ii) Also, the lower bounds of a_n-k^(n) for 1≤ k< n follow from
the definition (<ref>) immediately. For the remaining coefficient
a_0^(n), since e^x>1+x for all real x and since ln(1-x/2)>-x
for 0<x<1, we find that

    (1-θ)^1-α=e^(1-α)ln(1-α/2)
    >1+(1-α)ln(1-α/2)>1-α(1-α)≥3/4,

and then
a_0^(n)=(1-θ)^1-αω_2-α(τ_n)/τ_n
>3/4τ_n∫_t_n-1^t_nω_1-α(t_n-s)s.




The discrete coefficients b^(n)_n-k defined in (<ref>) satisfy

    0<b^(n)_n-k≤ρ_k/4(1+ρ_k)∫_t_k-1^t_kϖ_n”(s)s,
          1≤ k≤ n-1.



Since 0<(s-t_k-1)(t_k-s)<τ_k^2/4 for t_k-1<s<t_k, the alternative definition
(<ref>) of b^(n)_n-k yields the result and completes the proof.



As an application of <ref>, the next lemma builds up a link between
a_n-k^(n) and b_n-k^(n).
For a uniform mesh t_n=nτ, this lemma gives


    0<b_n-k^(n)<θ/4(n-θ-k)a^(n)_n-k  for 1≤ k≤ n-1.

By comparison, Alikhanov <cit.> gives
b_n-k^(n)<θ/2(1-θ)a^(n)_n-k.




Obviously, the new bound is much sharper.


The positive coefficients a^(n)_n-k, b^(n)_n-k defined in (<ref>) and (<ref>) satisfy

    b^(n)_n-k<θτ_k/2(t_n-θ-t_k) ρ_k/1+ρ_ka^(n)_n-k,   1≤ k≤ n-1.


For fixed n and 1≤ k≤ n-1, consider an auxiliary function

    φ_k(z)∫_t_k-1^t_k-1+zϖ_n”(s)s
        -2θ/t_n-θ-t_k∫_t_k-1^t_k-1+zϖ_n'(s)s,   0<z<τ_k.

Since ϖ_n”(t)=αϖ_n'(t)/(t_n-θ-t), the first derivative

    φ_k'(z)   =ϖ_n'(t_k-1+z)(
        α/t_n-θ-t_k-1-z
            -2θ/t_n-θ-t_k)
       <ϖ_n'(t_k-1+z)α-2θ/t_n-θ-t_k=0,
              0<z<τ_k, 1≤ k≤ n-1.

Hence the definition (<ref>) of a_n-k^(n) yields

    ∫_t_k-1^t_kϖ_n”(s)s-2θτ_k/t_n-θ-t_ka^(n)_n-k
        =φ_k(τ_k)<φ_k(0)=0,  1≤ k≤ n-1.

<ref> gives the claimed inequality and completes the proof.


Now we verify <ref> (I) by using <ref>.

Under the assumption 1, one has
θ<1-θ, ρ_k≤7/4  and
t_n-θ-t_k≥(1-θ)τ_k+1 for 1≤ k≤ n-1.
Thus, by using <ref>, one has

    b^(n)_n-k
    	<θτ_k/2(1-θ)τ_k+1 ρ_k /1+ρ_ka^(n)_n-k≤7ρ_k/8(1+ρ_k)a^(n)_n-k≤7/11a^(n)_n-k,   1≤ k≤ n-1,

since the function t/(1+t) is increasing for any t>0.
By <ref> (i), a_1^(n)<a_0^(n), then the definition (<ref>) yields

    A^(n)_0=a^(n)_0+ρ_n-1b^(n)_1≤ a^(n)_0+4944a^(n)_1≤24/11a_0^(n).

So the definition (<ref>)
of a_0^(n) gives the upper bound

    A^(n)_0≤24/11τ_n(1-θ)^1-α∫_t_n-1^t_nω_1-α(t_n-s) s ≤24/11τ_n∫_t_n-1^t_nω_1-α(t_n-s) s.

The lower bounds of A^(n)_n-k for 1≤ k≤ n-1 follow from
<ref> (ii) because
the definition (<ref>) implies that

    A^(n)_n-k≥ a^(n)_n-k-b^(n)_n-k≥4/11a^(n)_n-k.

The proof of <ref> (I) is complete.





 §.§ Proof of <ref> (II)–(III)

For simplicity of presentation, in this subsection we let

    I_n-k^(n):=∫_t_k-1^t_kt_k-t/τ_kϖ_n”(t)t    
    J_n-k^(n):=∫_t_k-1^t_kt-t_k-1/τ_kϖ_n”(t) t

for 1≤ k≤ n-1.

For 1≤ k≤ n-2, the positive coefficients b^(n)_n-k in (<ref>) satisfy


  (i)
I_n-k^(n)≥1+ρ_k/ρ_kb^(n)_n-k;
  (ii)  
J_n-k^(n)≥2(1+ρ_k)/ρ_kb^(n)_n-k;
  (iii)  
J_n-k^(n)≥ I_n-k^(n).


The alternative definition
(<ref>) of b^(n)_n-k gives the result (i) directly since
0<s-t_k-1<τ_k for s∈(t_k-1,t_k).
Since ϖ_n”'(t)>0 for 0<t<t_n-θ,
we take q:=ϖ_n' in Lemma <ref> to find

    ∫_t_k-1^t_k⟨s-t_k-1/τ_k-1/2|ϖ_n”(s)s
    =1/2τ_k∫_t_k-1^t_k(s-t_k-1)(t_k-s)ϖ_n”'(s)s>0,

and then J_n-k^(n)>1/2∫_t_k-1^t_kϖ_n”(s)s for  1≤ k≤ n-1.



So the inequality (ii) follows immediately from <ref>. Moreover,
2J^n_n-k>∫_t_k-1^t_kϖ_n”(s)s
=I_n-k^(n)+J_n-k^(n) so the claimed result (iii) follows directly.


For any fixed n (3≤ n≤ N) and 1≤ k≤ n-2, it holds that


  (i)
I_n-k-1^(n)≥1/ρ_kI_n-k^(n);   (ii)   J_n-k-1^(n)≥1/ρ_kJ_n-k^(n).









For fixed n≥2, introduce an auxiliary function with respect to z∈[0,1],

    ψ_k(z)1/τ_k∫_t_k-1^t_k-1+zτ_k⟨t_k-1+zτ_k-s|ϖ_n”(s)s,   1≤ k≤ n-1,

with the first and second derivatives

    ψ_k'(z)=∫_t_k-1^t_k-1+zτ_kϖ_n”(s)s,  ψ_k”(z)=τ_kϖ_n”(t_k-1+zτ_k),   1≤ k≤ n-1.

Note that ψ_k(0)=ψ_k'(0)=0 for 1≤ k≤ n-1,
and ψ_k+1(0)=ψ_k+1'(0)=0 for 0≤ k≤ n-2.
Thanks to the Cauchy differential
mean-value theorem, there exist z_1k, z_2k∈(0,1) such that

    I^(n)_n-k-1/I^(n)_n-k   =ψ_k+1(1)/ψ_k(1)
    	=ψ_k+1(1)-ψ_k+1(0)/ψ_k(1)-ψ_k(0)
    	=ψ_k+1'(z_1k)/ψ_k'(z_1k)
    	=ψ_k+1'(z_1k)-ψ_k+1'(0)/ψ_k'(z_1k)-ψ_k'(0)
       =ψ_k+1”(z_2k)/ψ_k”(z_2k)
         =τ_k+1ϖ_n”(t_k+z_2kτ_k+1)/τ_kϖ_n”(t_k-1+z_2kτ_k)≥1/ρ_k,
            1≤ k≤ n-2,

because ϖ_n”(t)>0 is increasing and
t_k>t_k-1+z_2kτ_k. The inequality (i) follows.
We now introduce another auxiliary function for z∈[0,1],

    ϕ_k(z)1/τ_k∫_t_k-1^t_k-1+zτ_k
    	(s-t_k-1)ϖ_n”(s)s,   1≤ k≤ n-1,

with the first derivative ϕ_k'(z)=zτ_kϖ_n”(t_k-1+zτ_k) for 1≤ k≤ n-1.
Then a similar argument yields the desired result (ii) and completes the proof.




The positive coefficients a^(n)_n-k in (<ref>) satisfy
 
    a_n-k-1^(n)-a_n-k^(n)
    =I_n-k-1^(n)+J_n-k^(n),   1≤ k≤ n-2  (3≤ n≤ N),

and for k=n-1 (2≤ n≤ N),

 
    a_0^(n)-a_1^(n)
    =θ/1-2θϖ_n'(t_n-1)+J_1^(n).


 For fixed n (3≤ n≤ N), applying the definition (<ref>),
we exchange the order of integration to find

    a_n-k-1^(n)-ϖ_n'(t_k)
    =    ∫_t_k^t_k+1ϖ_n'(s)-ϖ_n'(t_k)/τ_k+1 s
    =∫_t_k^t_k+1∫_t_k^sϖ_n”(t)/τ_k+1t s=I_n-k-1^(n)

for 0≤ k≤ n-2, and similarly,

    a_n-k^(n)-ϖ_n'(t_k)
    =    ∫_t_k-1^t_kϖ_n'(s)-ϖ_n'(t_k)/τ_k s
    
    =-J_n-k^(n)

for 1≤ k≤ n-1 (2≤ n≤ N). Hence the desired first equality is obtained by a simple subtraction.
For the case of k=n-1 (2≤ n≤ N), the above equality gives

    a_1^(n)-ϖ_n'(t_n-1)=-J_1^(n).

We have a^(n)_0=1-θ/1-αϖ_n'(t_n-1) such that
 a^(n)_0-ϖ_n'(t_n-1)=θ/1-2θϖ_n'(t_n-1).
 Thus a simple subtraction yields the second equality and completes the proof.
 



If 1 holds, the positive coefficients a^(n)_n-k in (<ref>) satisfy
 
    a_n-k-1^(n)-a_n-k^(n)≥
    	 b_n-2^(n)+6/5I_n-1^(n),
    	   k=1,
    
    
    	 b_n-k-1^(n)+ρ_k-1b^(n)_n-k+1+1/5I_n-k^(n),
    	   k>1,

for 1≤ k≤ n-2 (3≤ n≤ N), and for k=n-1 (2≤ n≤ N),
 
    a_0^(n)-a_1^(n)≥
    	 I_1^(2),
    	   n=2,
    
    ρ_n-2b^(n)_2+I_1^(n),
    	   n>2.



For fixed n, applying <ref> (i) and <ref> (i), we
obtain

    I_n-k-1^(n)   =ρ_k+1I_n-k-1^(n)/1+ρ_k+1
    	+I_n-k-1^(n)/1+ρ_k+1≥ b_n-k-1^(n)+I_n-k^(n)/ρ_k(1+ρ_k+1)
       ≥ b_n-k-1^(n)+I_n-k^(n)/ρ(1+ρ)
    =b_n-k-1^(n)+16/77 I_n-k^(n),
       1≤ k≤ n-2,

where the assumption 1 was used. By using <ref> (ii) and
<ref> (ii),

    ρ_k-1^3/2(1+ρ_k-1)J_n-k^(n)≥ρ_k-1^2/2(1+ρ_k-1)J_n-k+1^(n)≥ρ_k-1b^(n)_n-k+1,   2≤ k≤ n-1.

Then, noting that 2+2x-x^3≥ 9/64 for x∈[0,7/4], we apply
<ref> (iii) and the assumption 1 to get

    J_n-k^(n)   =ρ_k-1^3/2(1+ρ_k-1)J_n-k^(n)
    	+2+2ρ_k-1-ρ_k-1^3/2(1+ρ_k-1)J_n-k^(n)≥ρ_k-1b^(n)_n-k+1+9/352I_n-k^(n),

where 2≤ k≤ n-1. Hence, with help of
(<ref>)–(<ref>), we apply
<ref> to find

    a_n-k-1^(n)-a_n-k^(n)   =I_n-k-1^(n)+J_n-k^(n)≥ b_n-k-1^(n)+ρ_k-1b_n-k+1
    	+9/352+16/77I_n-k^(n)
       >b_n-k-1^(n)+ρ_k-1b_n-k+1^(n)+1/5I_n-k^(n),
    	   2≤
    k≤ n-2.

If k=1, by applying <ref> with the bound (<ref>)
and <ref> (iii), one has

    a_n-2^(n)-a_n-1^(n)=    I_n-2^(n)+J_n-1^(n)≥ b_n-2^(n)+16/77I_n-1^(n)+I_n-1^(n)≥
    b_n-2^(n)+6/5I_n-1^(n).


To complete the proof, it remains to consider the case of k=n-1 (2≤ n≤ N).
If n=2, <ref> and <ref> (iii) yield
 
    a_0^(2)-a_1^(2)
    =θ/1-2θϖ_n'(t_1)+J_1^(2)>J_1^(2)>I_1^(2).

Now treat the last case of n≥3.
We apply <ref> (by taking k=n-2), <ref> (i) and
the given condition 1 to get

    ρ_n-2b_2^(n)≤θτ_n-2/2(t_n-θ-t_n-2)ρ_n-2^2/1+ρ_n-2a_2^(n)≤θρ_n-2/2ρ_n-2^2/1+ρ_n-2a_2^(n)≤θρ^3/2(1+ρ)a_0^(n)
    
    	=343/352 θ a^(n)_0
        <θ/τ_nω_2-α(t_n-θ-t_n-1)
    	=θ(1-θ)/1-2θϖ_n'(t_n-1)
    	≤θ/1-2θϖ_n'(t_n-1).

Therefore  <ref> and <ref> (iii) lead to
 
    a_0^(n)-a_1^(n)
    =θ/1-2θϖ_n'(t_n-1)+J_1^(n)>ρ_n-2b_2^(n)+I_1^(n).

The proof is completed.


Recalling the definition (<ref>), we proceed to
apply <ref>.


With the notation I_n-k^(n) defined in (<ref>),
we can write the desired inequality as

    A^(n)_n-k-1-A^(n)_n-k≥(1+ρ_k)b^(n)_n-k+1/5I_n-k^(n),
          1≤ k≤ n-1,

and treat four separate cases covering all possibilities.


Indeed, from the definition (<ref>) of A^(n)_n-k, it is not
difficult to verify that

  
  (1) If k=1 for n=2,
  
    A^(2)_0-A^(2)_1=(1+ρ_1)b^(2)_1+a^(2)_0-a^(2)_1;

  
  (2) If k=n-1 for n≥3,
  
    A^(n)_0-A^(n)_1=(1+ρ_n-1)b^(n)_1+a^(2)_0-a^(2)_1-ρ_n-2b^(n)_2;

  
  (3) If k=1 for n≥3,
  
    A^(n)_n-2-A^(n)_n-1=(1+ρ_1)b^(n)_n-1+a^(n)_n-2-a^(n)_n-1-b^(n)_n-2;

  
  (4) If 2≤ k≤ n-2 for n≥4,
  
    A^(n)_n-k-1-A^(n)_n-k
    	=(1+ρ_k)b^(n)_n-k+a^(n)_n-k-1-a^(n)_n-k-b^(n)_n-k-1-ρ_k-1b^(n)_n-k+1 .


The claimed inequality follows from <ref> directly
and completes the proof.



The proof of <ref> shows that

    1-2θ/1-θa^(n)_0-a^(n)_1=J^(n)_1>0
    	  for 2≤ n≤ N.

In the case n=2, the definition (<ref>) gives

    1-2θ/1-θ A^(2)_0-A^(2)_1
    	   =1-2θ/1-θ (a^(2)_0+ρ_1b^(2)_1)
    	-(a^(2)_1-b^(2)_1)
       =J^(2)_1+1-2θ/1-θ ρ_1b^(2)_1+b^(2)_1>0.

For n≥3, one has

    1-2θ/1-θ A^(n)_0-A^(n)_1
    	   =1-2θ/1-θ (a^(n)_0+ρ_n-1b^(n)_1)
    	-(a^(n)_1+ρ_n-2b^(n)_2-b^(n)_1)
       =J^(n)_1-ρ_n-2b^(n)_1
    		+1-2θ/1-θ ρ_n-1b^(n)_1+b^(n)_1>0

because J^(n)_1≥ρ_n-1b^(n)_2+9/352I^(n)_1 from the
case k=n-1 of (<ref>).






§ NUMERICAL EXPERIMENTS


A numerical example is reported here to support our theory.
The fully discrete scheme (<ref>) is used to
solve the subdiffusion problem (<ref>)
in the domain Ω=(0,π) and T=1.
We take κ=2 and set the exact solution
 u(x,t)=1+ω_1+σ(t)sin(x).
This solution satisfies a stronger estimate than (<ref>), namely,
u^(ν)(t)_H^2(Ω)≤ C_ut^σ-ν for 0<t≤ T and ν∈{1,2,3}.
As noted in <cit.>, the  graded
mesh t_n=T⟨n/N|^γ satisfying 1–2,
is optimal in resolving the initial singularity.












In our computations, a linear finite element approximation is applied on a uniform mesh for Ω with M nodes.
As done in an earlier paper <cit.>, we split the interval [0,T]
into two parts [0,T_0]∪[T_0,T]. In first part [0,T_0] we used the smoothly graded
mesh t_n=(n/N_0)^γ T_0 for 0≤ n≤ N_0, while a uniform mesh with step size 
is used in the second part [T_0,T].  For a given total
number N of time levels, we put

    T_0 2^-γ  and  
    N_0⌈γ N/2^γ-1+γ⌉  so that  τT-T_0/N-N_0≥γ T_0/N_0≥τ_N_0.

To avoid problems with roundoff, the discrete coefficients a^(n)_n-k and
b^(n)_n-k from (<ref>) and (<ref>), respectively,
were computed using adaptive Gauss–Kronrod quadrature.











Since the O(h^2) behaviour of the spatial error is standard, we
fixed M=10^4 so that the temporal error dominates when N≤2,048.  Thus,
from <ref>, we expect the L_∞(L_2)-error e(N)max_1≤ n≤ Nu^n_h-u(t_n) to
behave like O(τ^min{γσ,2}).  We tested the sharpness of this
prediction by four scenarios:


  Table <ref>: σ=1+α and γ=1 with
fractional orders α=0.4, 0.6 and 0.8.

  Table <ref>: σ=1.2 and α=0.4 with mesh
parameters γ=1, 5/3 and 2.

  Table <ref>: σ=0.8 and α=0.4 with mesh
parameters γ=2, 5/2 and 3.

  Table <ref>: σ=0.4 and α=0.4 with mesh
parameters γ=2, 5/2 and 5.





The empirical order of convergence, listed as “Order” in the tables,
was computed in the usual way by supposing that e(N)≈ C_uτ^q and
evaluating the convergence rate q≈log_2[e(N)/e(2N)].  The
optimal mesh parameter 2/σ is the smallest value
of γ for which we expect second-order convergence; for
γ>
we still expect second-order convergence but with a constant factor that grows
with γ.  The convergence behaviour is always as expected, but it
is interesting to observe that, for larger values of σ (corresponding to
a less singular solution), the order can be close to 2 on the coarser grids.
In such cases, the predicted convergence order is not observed until the total number N of time levels is
quite large.



§ CONCLUDING REMARKS


An ECS analysis, including a reasonable ECS hypothesis
and a global consistency error, is proposed for investigating a class of
numerical approximations to the Caputo fractional derivative, employing
piecewise polynomial interpolation on general nonuniform time meshes. The
global consistency error bound (<ref>)
reveals an interesting behavior: the global approximate error of numerical
Caputo formula would be “local" despite being naturally
nonlocal, by choosing the time mesh
according to the error equidistribution principle. The effectiveness of the ECS analysis is shown for the
familiar L1 formula in <cit.>, and for the
higher-order Alikhanov formula in this paper. In the latter case, the
theoretical properties of the discrete convolution kernels (<ref>)
and a new interpolation error formula with integral remainder for quadratic
polynomials (<ref>) are crucial to obtaining a
useful ECS bound of the local truncation error.

As our answer to the problem formulated in <ref>, the
fractional Crank–Nicolson time-stepping scheme (<ref>)
for the linear reaction-subdiffusion equation (<ref>) is stable
(<ref>) and convergent (<ref>) on general
nonuniform grids satisfying a mild restriction (ρ=7/4) on the local time
step-size ratio. Consequently, some adaptive step-size control
with τ_next/τ_current≥4/7
is permitted to resolve certain complex behaviors of the solution arising in
nonlinear time-fractional PDEs, but a sudden, drastic reduction of the step size
should be avoided to ensure stability.  For the linear case, a sharp
L^2-norm error estimate of the numerical solution is also presented
(<ref>), demonstrating that a graded mesh can effectively
resolve the initial singularity.

The present results lead naturally to another question: can an
ECS analysis be applied to other high-order numerical Caputo formulas
such as the BDF2-like (L1-2) approximation
<cit.>?
We plan to address this issue in further research.



§ ACKNOWLEDGEMENTS

Hong-lin Liao thanks Prof. Ying Zhao for her valuable discussions and fruitful
suggestions, and the hospitality of Beijing Computational Science
Research Center during the period of his visit.









§ PROOF OF <REF>





For fixed n and 1≤ k≤ n-1, let ℓ_k,j(t) (j=k-1,k,k+1) be the
basis functions of quadratic Lagrange interpolation Π_2,kv at the
points t_k-1, t_k and t_k+1.




Firstly, we will express the interpolation error
Π_2,kv(t)=v(t)-⟨Π_2,kv|(t) in an integral form.
To do so, recall two basic properties of basis functions, ℓ_k,j(t_l)=δ_jl and

    ∑_j=k-1^k+1ℓ_k,j(t)(t_j-t)^ν=δ_0ν,  ν∈{0,1,2} ,

where δ_jl and δ_0ν are Kronecker delta functions.
Now applying the Taylor's expansion with integral remainder, one has

    v(t_j)   =∑_m=0^2v^(m)(t)/m!(t_j-t)^m
    +1/2∫_t^t_j(t_j-s)^2v”'(s)s ,   j∈{k-1,k,k+1} .

Since ⟨Π_2,kv|(t)=∑_j=k-1^k+1ℓ_k,j(t)v(t_j),
a simple combination with the three weights (basis functions) ℓ_k,j(t) (j=k-1,k,k+1) gives
the interpolation error

    Π_2,kv(t)=1/2∑_j=k-1^k+1∫_t_j^tℓ_k,j(t)(t_j-s)^2v”'(s)s,   t_k-1≤ t≤ t_k+1,

because the property (<ref>) implies that

    ∑_j=k-1^k+1ℓ_k,j(t)∑_m=0^2v^(m)(t)/m!(t_j-t)^m
    =∑_m=0^2v^(m)(t)/m!∑_j=k-1^k+1ℓ_k,j(t)(t_j-t)^m=v(t) .

Furthermore, differentiating both sides of (<ref>), one
applies (<ref>) again to get

    Π_2,kv'(t)=    1/2v”'(t)∑_j=k-1^k+1ℓ_k,j(t)(t_j-t)^2
    +1/2∑_j=k-1^k+1∫_t_j^tℓ_k,j'(t)(t_j-s)^2v”'(s)s
    
    =    1/2∑_j=k-1^k+1∫_t_j^tℓ_k,j'(t)(t_j-s)^2v”'(s)s=∑_j=k-1^k+1L_j(v)

for t_k-1≤ t≤ t_k+1, where

    L_j(v):=1/2∫_t_j^tℓ_k,j'(t)(t_j-s)^2v”'(s)s ,  
    j∈{k-1,k,k+1} .

Secondly, we express the required integration error
∫_t_k-1^t_kq'(t)Π_2,kv'(t)t
in terms of Π_1,kq by using
(<ref>).





Since

    ℓ_k,k+1'(t)=2/(τ_k+1+τ_k)τ_k+1⟨t-t_k-1/2|,

<ref> yields

    ∫_t_k-1^t_kℓ_k,k+1'(t)q'(t) t=    
    2∫_t_k-1^t_kt-t_k-1/2q'(t) t/(τ_k+1+τ_k)τ_k+1
    
    =-2∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k+1 .

Thus applying the formula for L_k+1(v), we
exchange the order of integration to find

    ∫_t_k-1^t_kq'(t)L_k+1(v)t =
    1/2∫_t_k-1^t_kℓ_k,k+1'(t)q'(t)t∫_t_k+1^t(t_k+1-s)^2v”'(s)s
        =-1/2∫_t_k-1^t_kℓ_k,k+1'(t)q'(t)t∫_t_k^t_k+1(t_k+1-s)^2v”'(s)s
        -1/2∫_t_k-1^t_kℓ_k,k+1'(t)q'(t)t∫_t_k^t(t_k+1-s)^2v”'(s)s
        =∫_t_k^t_k+1(t_k+1-s)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k+1
        -1/2∫_t_k-1^t_k(t_k+1-s)^2∫_t_k-1^sℓ_k,k+1'(t)q'(t)tv”'(s)s .

Similarly, it is easy to check the following equality

    ℓ_k,k-1'(t)=    2t-t_k-t_k+1/(τ_k+1+τ_k)τ_k
    =2/(τ_k+1+τ_k)τ_kt-t_k-1/2-1/τ_k ,

and it follows that

    ∫_t_k-1^t_kℓ_k,k-1'(t)q'(t) t=    
    2∫_t_k-1^t_kt-t_k-1/2q'(t) t/(τ_k+1+τ_k)τ_k
    -1/τ_k∫_t_k-1^t_kq'(t) t
    
    
    =    -2∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k
    -1/τ_k∫_t_k-1^t_kq'(t) t .

So applying the formula of L_k-1(v), we
exchange the order of integration to find

    ∫_t_k-1^t_kq'(t)L_k-1(v)t =
    1/2∫_t_k-1^t_kq'(t)∫_t_k-1^tℓ_k,k-1'(t)(t_k-1-s)^2v”'(s)st
       =1/2∫_t_k-1^t_k(t_k-1-s)^2v”'(s)s∫_s^t_kℓ_k,k-1'(t)q'(t) t
       =1/2∫_t_k-1^t_k(t_k-1-s)^2v”'(s)s∫_t_k-1^t_kℓ_k,k-1'(t)q'(t) t
       -1/2∫_t_k-1^t_k(t_k-1-s)^2v”'(s)s∫_t_k-1^sℓ_k,k-1'(t)q'(t) t
       =-∫_t_k-1^t_k(t_k-1-s)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1+τ_k)τ_k
       -1/2∫_t_k-1^t_k∫_t_k-1^t_kq'(t)/τ_k t
    +∫_t_k-1^sℓ_k,k-1'(t)q'(t) t(t_k-1-s)^2v”'(s)s.

For the remaining term involving L_k(v), one has

    ∫_t_k-1^t_kq'(t)L_k(v)t
    =    1/2∫_t_k-1^t_kq'(t)∫_t_k^tℓ_k,k'(t)(t_k-s)^2v”'(s)st
    
    =    -1/2∫_t_k-1^t_k(t_k-s)^2∫_t_k-1^sℓ_k,k'(t)q'(t) tv”'(s)s.

Then collecting the three equalities
(<ref>)–
(<ref>),
one applies the formula (<ref>) to get

    ∫_t_k-1^t_kq'(s)Π_2,kv'(s) s
    =∫_t_k^t_k+1(t_k+1-s)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1
        +τ_k)τ_k+1
    
    -∫_t_k-1^t_k(s-t_k-1)^2v”'(s)s∫_t_k-1^t_kΠ_1,kq(t) t/(τ_k+1
        +τ_k)τ_k+∫_t_k-1^t_k𝒦_q(s)v”'(s)s,

where the integral kernel

    𝒦_q(s):=
    -(t_k-1-s)^2/2τ_k∫_t_k-1^t_kq'(t) t
    -1/2∑_j=k-1^k+1(t_j-s)^2∫_t_k-1^sℓ_k,j'(t)q'(t)t .

Finally, to complete the proof, it remains to verify


    𝒦_q(s)=∫_t_k-1^sΠ_1,kq(t)t ,   t_k-1≤ s≤ t_k.

Differentiating the identity
(t-s)^2=∑_j=k-1^k+1(t_j-s)^2ℓ_k,j(t), we have

    ∑_j=k-1^k+1(t_j-s)^2ℓ_k,j'(t)=2(t-s) ,   t_k-1≤ s≤ t_k+1.

Thus it follows from (<ref>) that

    𝒦_q(s)=
    -(t_k-1-s)^2/2τ_k∫_t_k-1^t_kq'(t) t
    -∫_t_k-1^s(t-s)q'(t)t .

We see that 𝒦_q(t_k-1)=0 and

    𝒦_q'(s)=
    q(s)-q(t_k-1)-q(t_k)-q(t_k-1)/τ_k⟨s-t_k-1|
    =Π_1,kq(s),   t_k-1≤ s≤ t_k,

which leads to the desired result (<ref>) immediately since
𝒦_q(s)=∫_t_k-1^s𝒦_q'(t)t.


siamplain

99

Alikhanov2015
A. A. Alikhanov, A new difference scheme for the time fractional
  diffusion equation, J. Comput. Phys., 280 (2015), pp. 424–438.

  Brunner1985
H. Brunner, The numerical solution of weakly singular Volterra integral
  equations by collocation on graded meshes, Math. Comp., 45 (1985), pp. 417–437.









GaoSunZhang:2014
G.H. Gao, Z.Z. Sun and H.W. Zhang, A new fractional differentiation
formula to approximate the Caputo fractional derivative and its
applications, J. Comput. Phys., 259 (2014), pp. 33-50.


Graham1982
I. G. Graham, Galerkin methods for second kind integral equations
  with singularities, Math. Comp., 39 (1982), pp. 519–533.


  FordYan:2017
N. Ford and Y. Yan.,
An approach to construct higher order time discretisation schemes for
time fractional partial differential equations with nonsmooth data,
Fract. Calc. Appl. Anal., 20 (2017), pp. 1076–1105.

JinLiZhou:2017
 B. Jin, B. Li and Z. Zhou,
 Correction of high-order BDF convolution quadrature for fractional evolution equations,
SIAM J. Sci. Comput., 39 (2017), pp. A3129-A3152.

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

LiaoLiZhang:2018
H.-L. Liao, D. Li and J. Zhang,
Sharp error estimate of nonuniform L1 formula for linear
reaction-subdiffusion equations,,
  SIAM J. Numer. Anal.,  56 (2018), pp. 1112–1133.


LiaoMcLeanZhang:2019
H.-L. Liao, W. McLean, and J. Zhang, A discrete
  Grönwall inequality with application to numerical schemes for
  fractional reaction-subdiffusion problems, SIAM J. Numer. Anal., 2019,
  arXiv:1803.09879v2.
  


LiaoLyuVongZhao:2016
H.-L. Liao, P. Lyu, S. W. Vong and Y. Zhao,
Stability of fully discrete schemes with interpolation-type fractional
formulas for distributed-order subdiffusion equations,
Numer. Algor., 75 (2017), pp. 845-878.

 LiaoYanZhang:2019
H.-L. Liao, Y. Yan, and J. Zhang, Unconditional
  convergence of a two-level linearized fast algorithm for nonlinear
  subdiffusion equations, J. Sci. Comput., 2019, arXiv:1803.09858v2.

LiaoZhaoTeng:2016
H.-L. Liao, Y. Zhao, and X. Teng, A weighted ADI scheme for
  subdiffusion equations, J. Sci. Comput., 69 (2016), pp. 1144–1164.



  LvXu:2016
C. Lv and C. Xu, Error analysis of a high order method
for time-fractional diffusion equations, SIAM J. Sci. Comput., 38 (2016),
pp. A2699-A2724.

McLean2010
W. McLean, Regularity of solutions to a time-fractional
diffusion equation, ANZIAM J., 52 (2010), pp. 123–138.

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

McLeanMustapha:2015
W. McLean and K. Mustapha, Time-stepping error bounds for fractional diffusion problems
with nonsmooth initial data, J. Comput. Phys., 293 (2015), pp. 201-217.


SakamotoYamamoto2011
K. Sakamoto and M. Yamamoto, Initial value/boundary value problems
  for fractional diffusion-wave equations and applications to some inverse
  problems, J. Math. Anal. Appl., 382 (2011), pp. 426–447.

Stynes2016
M. Stynes, Too much regularity may force too much uniqueness, Frac.
  Calc. Appl. Anal., 19 (2016), pp. 1554–1562.

StynesEtAl2017
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 (2017), pp. 1057–1079.

ZengEtAl2015
F. Zeng, C. Li, F. Liu, and I. Turner, Numerical algorithms for a
  time-fractional subdiffusion equation with second-order accuracy, SIAM J.
  Sci. Comput., 37 (2015), pp. A55–A78.




