A discrete Grönwall inequality with application to numerical schemes for subdiffusion problems
Hong-lin Liao William McLean Jiwei Zhang

Jittering jets in cooling flow clusters
    Noam Soker1
    December 30, 2023
=======================================


We consider a class of numerical approximations to the Caputo fractional
derivative. Our assumptions permit the use of nonuniform time steps, such as is
appropriate for accurately resolving the behavior of a solution whose
temporal derivatives are singular at t=0.  The main result is a type of fractional
Grönwall inequality and we illustrate its use by outlining some stability and
convergence estimates of schemes for fractional reaction-subdiffusion problems.
This approach extends earlier work that used the familiar L1 approximation to
the Caputo fractional derivative, and will facilitate the analysis of higher
order and linearized fast schemes.



fractional subdiffusion equations,  nonuniform time mesh,
discrete Caputo derivative, discrete Grönwall inequality.



65M06, 35B65




§ INTRODUCTION

This paper builds on earlier results <cit.> for
the nonuniform L1 method applied to the time discretization of a fractional
reaction-subdiffusion problem <cit.> in a spatial domain Ω,

    α u+ u   =f(x, t, u)      for x∈Ω and 0<t≤ T,
    
    u   =u_0(x)      for x∈Ω when t=0,
    
    u   =0      for x∈∂Ω and 0<t<T.

Here, α=^C_0α denotes the Caputo fractional
derivative of order α with respect to time t, with 0<α<1,
and  is a linear, second-order, strongly-elliptic partial differential
operator in the spatial variable(s) x. We establish a discrete
Grönwall inequality intended for the error analysis of higher-order
time discretizations <cit.> and linearized fast algorithms
<cit.> for solving (<ref>)
that employ nonuniform step sizes.

In any numerical methods for solving the
reaction-subdiffusion problem (<ref>), a key consideration
is that the solution u(x,t) is typically less regular than would be the case
for a classical parabolic PDE (which arises as the limiting case
when α→1).  For example, in the simplest case f(x,t,u)≡0
when (<ref>) is linear and homogeneous, let
φ_ be a Dirichlet eigenfunction of  on Ω, with
eigenvalue λ_>0, so that
φ_=λ_φ_. Let E_α
denote the Mittag–Leffler function,

    E_α(z)∑_k=0^∞z^k/Γ(1+kα),

and choose as the initial data u_0(x)=φ_(x).  Term-by-term
differentiation shows that the solution is
u(x,t)=E_α(-λ_ t^α)φ_(x), and so
∂ u/∂ t=O(t^α-1) as t→0, whereas the solution of
the classical parabolic equation,
u(x,t)=e^-λ_tφ_(x), is a smooth function of t.
Sakamoto and Yamamoto <cit.> study the (lack of) regularity
of u for more general initial data u_0 and a linear source term f=f(x,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 nature of polynomial interpolation means that the
convergence rate of the L1 or similar approximations to α u is
limited by the smoothness of the solution u.  In the presence of a fixed
singularity at t=0 of the type described above, an established technique to
restore an optimal convergence rate is to employ a graded mesh

    t_n(n/N)^γ T  for 0≤ n≤ N,

where the parameter γ≥1 must be adapted to the strength of the
singularity.  Choosing γ=1 results in a uniform mesh, and the larger the
value of γ the more strongly the grid points are concentrated near t=0.
For example, 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 <cit.> is now well
established.

Early papers on L1 schemes <cit.> assumed a uniform
step size τ, and showed that if u is smooth then the time discretization
error is O(τ^2-α). Recently, Jin, Lazarov and
Zhou <cit.> presented a new analysis, based on generating
functions, that permitted nonsmooth initial data u_0. They showed that if
f≡0 and u_0∈ L_2(Ω), then the error in the norm of L_2(Ω)
due to the time discretization is O(τ t_n^-1).  Thus, for t_n bounded
away from zero, the method achieves first-order accuracy in time.
Yan, Khan and Ford <cit.> proposed a modified L1 scheme
and obtained error estimates for smooth and nonsmooth initial data.
It was shown that the modified L1 scheme on a uniform mesh has a convergence
rate of O(τ^2-α).
Aliknanov <cit.> introduced the L2-1_σ formula, a
modification of the L1 method that uses
piecewise-quadratic instead of piecewise-linear interpolation, and approximates
α u at an offset grid point t_j+σ=(j+σ)τ.  He
showed that if u is sufficiently smooth then the time discretization error is
O(τ^2) for the special choice σ=1-α/2.

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.  Stynes, O'Riordan
and Gracia <cit.> considered the L1 method on a graded
mesh of the form (<ref>) applied
to (<ref>) for the case u=-u_xx and a linear
reaction term f(x,t,u)=-c(x)u+g(x,t). They showed that, given the typical
singular behavior of u, the maximum error in the fully-discrete solution is
of order N^-min{2-α,γα}.  (Here we ignore the additional
error due to the spatial discretization.)  Thus, for a uniform mesh the
error is O(N^-α), but if γ=(2-α)/α
then the error is O(N^α-2). Their stability
analysis requires c(x)≥0, which prevents extending the approach to deal with
a reaction term that is nonlinear but uniformly Lipschitz in u.  This
limitation was overcome recently in the precursor <cit.> to the
present work by exploiting a novel discrete fractional Grönwall inequality for the L1
method.


Nonetheless, practical applications of the discrete
Grönwall inequality in its basic form <cit.> are still
limited because it does not apply to other numerical approximation schemes for
the Caputo derivative and excludes certain adaptive time meshes required
to resolve complex behaviors (physical oscillations, blowup and so on)
in nonlinear fractional differential equations.  Also, the proof relies on
specific properties of the L1 kernels a^(n)_n-k and their complementary
discrete kernels P^(n)_n-k, with a key
step <cit.> employing rough estimates of the
truncation error that, to a large extent, rely on the simple form
of the a^(n)_n-k.  In summary, the main novel contributions of the
present work are threefold:


  (i) to generalize the discrete Gronwall inequality, permitting its use
with a variety of discretizations of the Caputo derivative, not just the L1
scheme;

  (ii) to provide a concise proof based on two simple assumptions on the
discrete kernels, independent of their precise form;

  (iii) to permit a more general class of nonuniform meshes or adaptive
time grids, not just the graded meshes for resolving the initial singularity.

In more detail, <ref> defines a discrete fractional
derivative (<ref>) having the form of the classical L1
approximation but with general discrete kernels.  We formulate three
assumptions required for our theory.  The first two impose a monotonicity
property (1) and a lower bound (2) on the discrete kernels, and the
third (3) places a mild restriction on the local step-size ratio.  We
give some examples of schemes satisfying these assumptions, and define a family
of complementary discrete kernels, generalizing those introduced in the earlier
paper <cit.>.  <Ref> establishes a key
estimate involving the discrete kernels and the Mittag–Leffler
function (<ref>).  In <ref>,
we prove our main result, a discrete fractional Grönwall inequality
stated as <ref>, and provide, in <ref>, a
strategy to treat cases where the monotonicity assumption breaks down.
<Ref> illustrates the use of the Gronwall inequality
in conjunction with an abstract Galerkin method for the spatial discretization.
Finally, a short appendix proves two technical inequalities needed for the
stability analysis of <ref>.

The generalized results proved below will allow us to show, in two
companion papers <cit.>, that
Alikhanov's L2-1_σ formula can achieve second-order accuracy on certain
nonuniform time grids and that a linearized fast algorithm is unconditionally
convergent for nonlinear subdiffusion equations.








§ DISCRETE FRACTIONAL DERIVATIVE

Recall that the Riemann–Liouville fractional integral operator of
order β>0 is defined by <cit.>

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

and, in turn, the Caputo fractional derivative is defined by

    (αv)(t)(ℐ^1-αv')(t)
    	=∫_0^tω_1-α(t-s)v'(s)s  for t>0.

For (possibly nonuniform) time levels
0=t_0<t_1<t_2<⋯<t_N=T, we denote the nth step size
by τ_n t_n-t_n-1, fix an offset parameter θ∈[0,1) and
define

    t_n-θθ t_n-1+(1-θ)t_n
      and  
    v^n-θθ v^n-1+(1-θ)v^n,

where v^k may be any sequence. Letting v^k≈ v(t_k) and
v^k v^k-v^k-1,
we consider a discrete Caputo derivative (not necessarily a direct approximation of (<ref>),
see Remark <ref>)
given by a convolution-like sum, as follows,

    (αv)^n-θ∑_k=1^n
    	A^(n)_n-k v^k  for 1≤ n≤ N.

Here, the corresponding discrete convolution kernels are written as A^(n)_n-k instead
of A_nk to reflect the convolution structure of the fractional derivative.

Our theory requires the following three
assumptions:


  A1. The discrete kernels are positive and monotone, that is,

    A^(n)_0≥ A^(n)_1≥ A^(n)_2≥⋯ A^(n)_n-1>0
      for 1≤ n≤ N.


  A2. There is a constant π_A>0 such that the discrete kernels satisfy
the lower bound

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


  A3. There is a constant >0 such that the step size
ratios ρ_kτ_k/τ_k+1 satisfy

    ρ_k≤  for 1≤ k≤ N-1.




The boundedness and monotonicity assumptions 1 and 2
on the discrete convolution kernels A^(n)_n-k are valid for several
frequently-used discrete Caputo derivatives, at least if assumption 3 is
satisfied for appropriate ρ.  Included are the well-known L1 formula
<cit.>, the
fast L1 formula <cit.>,
the Alikhanov approximation <cit.>,
and their applications for multi-term and distributed-order Caputo derivatives
(see Remark <ref>). Here we list three examples on nonuniform grids.
Note that, the local mesh parameter ρ from 3 will always
appear in our discrete fractional Grönwall inequality and our stability estimates.


[nonuniform L1 formula]
The widespread L1 formula <cit.>
uses θ=0 and v'(s)≈ v^k/τ_k (linear interpolation) to obtain

    (αv)^n∑_k=1^na^(n)_n-k v^k  with  
    a^(n)_n-k1/τ_k∫_t_k-1^t_kω_1-α(t_n-s)s.

This sum has the desired form (<ref>) where (using the
integral mean value theorem),

    A^(n)_n-k:=a^(n)_n-k=ω_1-α(t_n-s_nk)  for some s_nk∈[t_k-1,t_k].

It follows that assumption 1 is satisfied, and
2 holds with π_A=1.

[fast L1 formula]
In the two-level fast L1 approximation <cit.>,
the sum-of-exponentials technique is applied to approximate the weakly singular
kernel ω_1-α(t-s). That is, for a user-given absolute tolerance
error ϵ≪1 and a cut-off time Δt>0, one determines a
positive integer N_q, positive quadrature nodes θ^ℓ and positive
weights ϖ^ℓ (1≤ℓ≤ N_q) such that

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

Then we use θ=0 and v'(s)≈ v^k/τ_k (linear interpolation) to obtain

    α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)
        +1/τ_k∫_t_k-1^t_k
        e^-θ^ℓ(t_k-s)∇_τu^k s,  
        k≥1,  1≤ℓ≤ N_q .

This approximation also has the form (<ref>) with θ=0,

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

If the tolerance error ϵ is small enough such that
ϵ≤min{1/3ω_1-α(T),α ω_2-α(1)},
then
<cit.> ensures that 1-2 hold true with π_A=3/2.

[nonuniform Alikhanov formula]
Let Π_1,kv be 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.
Taking a special choice θ=α/2, and applying the linear and
quadratic polynomial interpolations, we have the nonuniform Alikhanov formula
 <cit.>

    (αv)^n-θ    ∑_k=1^n-1∫_t_k-1^t_kω_1-α(t_n-θ-s)⟨Π_2,kv|'(s)s
        +∫_t_n-1^t_n-θω_1-α(t_n-θ-s)
    		⟨Π_1,nv|'(s)s  for n≥1.

This formula can be written as the form (<ref>) with
A_0^(1)â_0^(1) for n=1 and, for n≥2,

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

where the discrete coefficients
â_n-k^(n) and b̂_n-k^(n) are defined by

    â^(n)_01/τ_n∫_t_n-1^t_n-θω_1-α(t_n-θ-s)s,
       â^(n)_n-k1/τ_k∫_t_k-1^t_kω_1-α(t_n-θ-s)s  for 1≤ k≤ n-1,
       b̂^(n)_n-k2/τ_k(τ_k+τ_k+1)∫_t_k-1^t_k
    	(s-t_k-1/2)ω_1-α(t_n-θ-s)s  for 1≤ k≤ n-1.

The theoretical properties in <cit.> assure
1–2 with π_A=11/4 provided the local mesh assumption 3
holds with the maximum step size ratio ρ=7/4.





We now continue to introduce an important tool: the complementary discrete convolution kernels.
The semigroup property of the fractional integral,
ℐ^αℐ^β=ℐ^α+β, holds because
the integral kernels satisfy ω_α*ω_β=ω_α+β.  It
follows that

    ∫_s^tω_α(t-μ)ω_1-α(μ-s)μ=ω_1(t-s)=1
    	  for all 0<s<t<∞,

and it motives us to seek a family
of complementary discrete convolution kernels P^(n)_n-j having the identical property

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

In fact, by taking m=k and m=k+1,

    P^(n)_n-kA^(k)_0+∑_j=k+1^nP^(n)_n-jA^(j)_j-k
    	=1=∑_j=k+1^nP^(n)_n-jA^(j)_j-(k+1) ,   1≤ k≤ n-1,

we see that

    P^(n)_n-k=1/A^(k)_0∑_j=k+1^n P^(n)_n-j
    	A^(j)_j-k-1-A^(j)_j-k,   1≤ k≤ n-1,

and the complementary discrete kernels may be defined via the
recursion <cit.>

    P^(n)_01/A^(n)_0, 
    P^(n)_j1/A^(n-j)_0∑_k=0^j-1
    	A^(n-k)_j-k-1-A^(n-k)_j-kP^(n)_k
    	  for 1≤ j≤ n-1.


[Pictures of A_j^(n) and P_j^(n) of L1 formula]
Consider the widespread L1 approximation in Example <ref>.
Figure <ref> plots the L1 discrete kernels A^(n)_j and the complementary discrete kernels P^(n)_j
 when T=1 and n=30 for three graded meshes of the form (<ref>).
 



As a consequence of the identity (<ref>), we find that

    ∫_0^tω_α(t-s)(α v)(s)s=∫_0^t v'(s)s,

which provides the inspiration for the second part of the next lemma.


Let the assumptions 1 and 2 hold.


  * The discrete kernels P^(n)_j
in (<ref>) having the property (<ref>) satisfy

    0≤ P^(n)_n-j≤π_AΓ(2-α)τ_n^α  for 1≤ j≤ n≤ N,

and

    ∑_j=1^nP_n-j^(n) ω_1-α(t_j)≤π_A  for 1≤ n≤ N.


  * If v:[0,T]→ is any continuous, piecewise-C^1 function such that
v' is non-negative and monotone decreasing, then

    ∑_j=1^n P^(n)_n-j(α v)(t_j)≤π_A∫_0^t_nv'(s)s  for 1≤ n≤ N.





It follows at once from  the monotonicity assumption 1 that
A^(n)_0>0 and A^(n-k)_j-k-1-A^(n-k)_j-k≥0 for 0≤ k≤ j-1.
The lower bound P^(n)_j≥0 is then clear from the
recursion (<ref>).  Since all the discrete kernels are
non-negative, we have

    P^(n)_n-kA^(k)_0≤∑_j=k^nP^(n)_n-jA^(j)_j-k=1

and taking n=k in the assumption 2 gives

    A^(k)_0≥1/π_Aτ_k∫_t_k-1^t_kω_1-α(t_k-s)s
    	=ω_2-α(τ_k)/π_Aτ_k
    	=1/Γ(2-α)π_Aτ_k^α,

so the complementary discrete convolution kernels P^(n)_n-k are well-defined and satisfy
the upper bound
P^(n)_n-k≤1/A^(k)_0≤Γ(2-α)π_Aτ_k^α.
Furthermore, the assumption A2 and the identity (<ref>) imply that
ω_1-α (t_j) ≤π_A A_j-1^(j) and

    ∑_j=1^nP_n-j^(n) ω_1-α(t_j)≤π_A ∑_j=1^nP_n-j^(n) A_j-1^(j) = π_A  for n≥1,

which completes the proof of part 1.

Recall Chebyshev's sorting
inequality <cit.>: if f is
monotone increasing and g is monotone decreasing on the interval [a,b], and
if both functions are integrable, then

    (b-a)∫_a^bf(s)g(s)s≤∫_a^b f(t) t∫_a^b g(s)s.

Taking [a,b]=[t_k-1,t_k], f(s)=ω_1-α(t_j-s) and g(s)=v'(s)≥0,
and using 2, we see that

    (α v)(t_j)   =∑_k=1^j∫_t_k-1^t_kω_1-α(t_j-s)v'(s)s
       ≤∑_k=1^j1/τ_k∫_t_k-1^t_kω_1-α(t_j-t) t
    		∫_t_k-1^t_kv'(s)s≤π_A∑_k=1^jA^(j)_j-k∫_t_k-1^t_kv'(s)s.

Thus, from the identical property (<ref>) of the discrete kernels P^(n)_n-j, we conclude that

    ∑_j=1^nP^(n)_n-j(α v)(t_j)
    	   ≤∑_j=1^nP^(n)_n-jπ_A∑_k=1^jA^(j)_j-k∫_t_k-1^t_kv'(s)s
       =π_A∑_k=1^n∫_t_k-1^t_kv'(s)s∑_j=k^nP^(n)_n-jA^(j)_j-k
    	=π_A∑_k=1^n∫_t_k-1^t_kv'(s)s,

and part 2 follows.  









When 3 also holds, we have a variant of the second part of
<ref>.


Let the assumptions 1–3 hold.  If v:[0,T]→ is any continuous,
piecewise-C^1 function
such that v' is non-negative and monotone, then

    ∑_j=1^n-1P^(n)_n-j(α v)(t_j)
    	≤max(1,)π_A∫_0^t_nv'(s)s  for 1≤ n≤ N.



If v' is non-negative and monotone decreasing, then α v(t_j)≥0 and the results of
<ref> imply that

    ∑_j=1^n-1P^(n)_n-j(α v)(t_j)
    	≤∑_j=1^nP^(n)_n-j(αv)(t_j)
    	≤π_A∫_0^t_nv'(s)s .

 Otherwise, if v' is
monotonely increasing, then

    ∑_j=1^n-1P^(n)_n-j(α v)(t_j)
    	   =∑_j=1^n-1P^(n)_n-j∑_k=1^j∫_t_k-1^t_kω_1-α(t_j-s)v'(s)s
       ≤∑_j=1^n-1P^(n)_n-j∑_k=1^j v'(t_k)
    	∫_t_k-1^t_kω_1-α(t_j-s)s
       ≤π_A∑_j=1^n-1P^(n)_n-j∑_k=1^j v'(t_k)τ_kA^(j)_j-k
       =π_A∑_k=1^n-1 v'(t_k)τ_k
    		∑_j=k^n-1P^(n)_n-jA^(j)_j-k≤π_A∑_k=1^n-1v'(t_k)τ_k
       ≤π_A∑_k=1^n-1v'(t_k)ρ_kτ_k+1≤π_A∑_k=1^n-1∫_t_k^t_k+1v'(s)s,

and the desired estimate again holds.


We can use <ref> to prove the following property of the
Mittag–Leffler function (<ref>).


Let the assumptions 1–3 hold. For any real μ>0,

    ∑_j=1^n-1P^(n)_n-jE_α(μ t_j^α)
    	≤π_Amax(1,)E_α(μ t_n^α)-1/μ  for 1≤ n≤ N.



The series definition (<ref>) shows that

    E_α(μ t^α)
    	=1+∑_k=1^∞μ^kt^kα/Γ(1+kα)
            =1+∑_k=1^∞μ^kv_k(t),

where v_k(t)=ω_1+kα(t) and we have v_k'(t)=ω_kα(t)>0 for all k≥1.
If 1≤ k≤1/α, then -1≤ kα-1≤0 and
v_k”(t)=ω_kα-1(t)≤0 for all t>0.  Otherwise, if k>1/α,
then kα-1>0 and v_k”(t)>0 for all t>0.  Thus, v_k' is always
non-negative and monotone, so we may apply <ref> and deduce
that

    ∑_j=1^n-1P^(n)_n-j(α v_k)(t_j)
        ≤max(1,)π_A∫_0^t_nv_k'(s)s
    	=max(1,)π_Av_k(t_n)  for k≥1.

Multiplying both sides of this inequality by μ^k, summing over the index k, and using
the fact that

    α v_k(t)=∫_0^tω_1-α(t-s)ω_kα(s)s
    	=ω_1+(k-1)α(t)=v_k-1(t)  for all k≥1,

we have

    ∑_k=1^mμ^k∑_j=1^n-1P^(n)_n-jv_k-1(t)
    	≤max(1,)π_A∑_k=1^mμ^kv_k(t_n).

Because the series ∑_k=1^∞μ^kv_k(t) is absolutely convergent and ω_1(t)=1,
the desired inequality follows after interchanging the
sums on the left-hand side and then sending m→∞. The proof is completed.





§ DISCRETE FRACTIONAL GRÖNWALL INEQUALITY

Our main result is stated in the next theorem.
The proof is similar to that of <cit.>,
but we include it here to incorporate the nonuniform mesh parameter ρ in
3, which does not appear in discrete Grönwall inequalities
for classical parabolic equations.

Let the assumptions 1–3 hold, let 0≤θ<1, and
let (g^n)_n=1^N and (λ_l)_l=0^N-1 be given non-negative
sequences.  Assume further that there exists a constant Λ
(independent of the step sizes) such that
Λ≥∑_l=0^N-1λ_l, and that the maximum step size satisfies

    max_1≤ n≤ Nτ_n
        ≤1/√(2π_AΓ(2-α)Λ) .

Then, for any non-negative sequence (v^k)_k=0^N such that

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

it holds that

    v^n≤2E_α2max(1,)π_A Λ t_n^α(v^0+max_1≤ k≤ n∑_j=1^k P^(k)_k-jg^j
    	)  for 1≤ n≤ N.



We replace the index n with j in (<ref>), then multiply
by P^(n)_n-j and sum over j to obtain

    ∑_j=1^nP^(n)_n-j∑_k=1^jA^(j)_j-kv^k^2
    	≤∑_j=1^nP^(n)_n-j∑_k=1^jλ_j-kv^k-θ^2+∑_j=1^nP^(n)_n-jv^j-θg^j.

On the left-hand side, we exchange the order of summation and use the
identity (<ref>) to get

    ∑_j=1^nP^(n)_n-j∑_k=1^jA^(j)_j-kv^k^2
    	   =∑_k=1^nv^k^2∑_j=k^nP^(n)_n-jA^(j)_j-k
       =∑_k=1^nv^k^2=(v^n)^2-(v^0)^2.

Thus, it follows from (<ref>) that

    v^n^2≤v^0^2+∑_j=1^nP^(n)_n-j∑_k=1^jλ_j-kv^k-θ^2
    +∑_j=1^nP^(n)_n-jv^j-θg^j	,


For brevity, let us write the claimed estimate (<ref>) as v^n≤ F_nG_n where

    F_n 2E_α(2max(1,)π_AΛ t_n^α)
      and  
    G_n v^0+max_1≤ k≤ n∑_j=1^k P^(k)_k-jg^j.

We will use complete induction, noting that the Mittag–Leffler function (<ref>)
satisfies E_α(0)=1 and E_α'(z)>0 for all real z>0, so F_n≥ F_n-1≥2 for n≥2.

If v^1≤ v^0, then v^1≤ G_1≤ F_1G_1, as required.
Otherwise, if v^1>v^0, then v^1-θ≤ v^1. One deduces from (<ref>) that

    v^1^2   ≤v^0^2+P^(1)_0 v^1-θg^1
        +P^(1)_0λ_0v^1-θ^2
       ≤ v^1v^0+P^(1)_0g^1
    		+P^(1)_0λ_0v^1^2=v^1G_1+P^(1)_0λ_0v^1^2.

Part 1 of <ref> and the given restriction on the
maximum time-step imply that

    P^(1)_0λ_0≤π_AΓ(2-α)τ_1^αΛ≤1/2.

Thus, v^1^2≤2v^1G_1 and so v^1≤2G_1≤ F_1G_1,
which implies that the desired estimate holds for n=1.

For the inductive step, let 2≤ n≤ N and assume that

    v^k≤ F_kG_k
    	  for 1≤ k≤ n-1.

Choose some k(n) such that v^k(n)=max_0≤ j≤ n-1v^j.  If
v^n≤ v^k(n) then, since F_k and G_k are monotone increasing in k,

    v^n≤ v^k(n)≤ F_k(n)G_k(n)≤ F_nG_n,

as required.  Otherwise, if v^n>v^k(n), then
v^j-θ≤max(v^j-1,v^j)≤ v^n for 1≤ j≤ n. We deduce
from (<ref>) that

    v^n^2≤ v^nv^0+v^n∑_j=1^n P^(n)_n-jg^j		
    +v^n∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kv^k-θ
    +v^n^2P^(n)_0∑_k=1^nλ_n-k.

Using part 1 of <ref>,

    P^(n)_0∑_k=1^nλ_n-k≤π_AΓ(2-α)Λτ_n^α,

so the limitation on the maximum step size implies that

    ⟨v^n|^2≤ v^nG_n
    +∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kv^k-θ
    	+1/2⟨v^n|^2.

Thus, applying the induction hypothesis (<ref>), we deduce
from (<ref>) that

    v^n   ≤2G_n+2∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kθ v^k-1+(1-θ)v^k
       ≤2G_n+2∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kθ F_k-1G_k-1+(1-θ)F_kG_k
       ≤2G_n+2∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kF_kG_k
    	≤2G_n+2∑_j=1^n-1P^(n)_n-jF_jG_j∑_k=1^jλ_j-k
       ≤2G_n+4Λ G_n-1∑_j=1^n-1P^(n)_n-j
    		E_α(2max(1,)π_AΛ t_j^α).

Finally, by <ref> with μ=2max(1,)π_AΛ,

    v^n≤ 2G_n+2max(1,)π_AΛ G_n E_α2max(1,)π_AΛ t_n^α-1/max(1,)π_AΛ
    	= F_nG_n,

which completes the inductive step and the proof.



One may use the inequality (<ref>) in part 1 of
<ref> to bound the convolutional summation
∑_j=1^k P^(k)_k-jg^j, that is,

    ∑_j=1^k P^(k)_k-jg^j≤∑_j=1^k P^(k)_k-jω_1-α(t_j)max_1≤ j≤ kg^j/ω_1-α(t_j)≤π_Amax_1≤ j≤ kg^j/ω_1-α(t_j)

So the discrete solution of (<ref>) can also be bounded by

    v^n≤2E_α2max(1,)π_A Λ t_n^αv^0+π_AΓ(1-α)max_1≤ j≤ n{t_j^αg^j}  for 1≤ n≤ N.

On the other hand, if the given sequence (λ_l)_l=0^N-1 is non-positive and the constant Λ≤0,
a similar argument will show that the discrete inequality (<ref>) holds in a simpler form,
requiring only the assumptions 1-2 but no restrictions on time steps,

    v^n≤ v^0+max_1≤ k≤ n∑_j=1^k P^(k)_k-jg^j≤ v^0+π_AΓ(1-α)max_1≤ j≤ n{t_j^αg^j}  for 1≤ n≤ N.




By including the non-negative sequence (λ_l)_l=0^N-1 in
(<ref>), we are able to treat various numerical approaches
to solving linear and nonlinear subdiffusion problems.
Typically, the sequence takes only a few nonzero values.  Recent examples
include λ_l=0 for l≥1 in the time-weighted method from
<ref>,
and λ_l=0 for l≥2 in the one-step linearized scheme
<cit.> for a semilinear subdiffusion equation. Thus, the
constant Λ is always not very large and the maximum time-step
restriction max_1≤ n≤ Nτ_n
≤1/√(2π_AΓ(2-α)Λ)
is also not stringent in practical applications.


The Mittag–Leffler function E_α also arises naturally in other discrete
and continuous Grönwall inequalities for fractional diffusion and wave
equations <cit.>, and for weakly singular Volterra
equations <cit.>.
The presence of the nonuniform mesh parameter ρ in the argument
of E_α indicates that sudden, drastic reductions of the time-step should
be avoided.  Nevertheless, 
our discrete Grönwall inequality does not restrict the heterogeneous degree of time mesh, this is, fits for general nonuniform mesh.


We also have an alternative version of the above theorem.


<ref> remains valid if the
condition (<ref>) is replaced by

    ∑_k=1^nA^(n)_n-k v^k≤∑_k=1^nλ_n-kv^k-θ+g^n
          for 1≤ n≤ N.

Moreover, if the given sequence (λ_l)_l=0^N-1 is non-positive and the constant Λ≤0,

    v^n≤ v^0+∑_j=1^n P^(n)_n-jg^j≤ v^0+π_AΓ(1-α)max_1≤ j≤ n{t_j^αg^j}  for 1≤ n≤ N.



The structure of proof is as before.  However, instead of
(<ref>) and (<ref>), we now have

    ∑_j=1^nP^(n)_n-j∑_k=1^jA^(j)_j-k v^k
        ≤∑_j=1^nP^(n)_n-j∑_k=1^jλ_j-kv^k-θ
        +∑_j=1^nP^(n)_n-jg^j

and

    ∑_j=1^nP^(n)_n-j∑_k=1^jA^(j)_j-k v^k
        =∑_k=1^n v^k∑_j=k^nP^(n)_n-jA^(j)_j-k
        =∑_k=1^n v^k=v^n-v^0,

respectively, so that instead of (<ref>) we obtain

    v^n≤ v^0
        +∑_j=1^nP^(n)_n-j∑_k=1^jλ_j-kv^k-θ
        +∑_j=1^nP^(n)_n-jg^j.

As before, if v^1≤ v^0 then v^1≤ G_1.  For the alternative
case v^1>v^0, we again have v^1-θ≤ v^1 which now yields

    v^1≤ v^0+P^(1)_0g^1+P^(1)_0λ_0v^1-θ
        =G_1+P^(1)_0λ_0v^1-θ≤ G_1+12 v^1,

where the final step again relies on the step size assumption to
ensure (<ref>).  Thus, once again, v^1≤2G_1.  In the inductive
step, (<ref>) is replaced by

    v^n≤ v^0+∑_j=1^n P^(n)_n-jg^j
        +∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kv^k-θ
        +v^nP^(n)_0∑_k=1^nλ_n-k,

and by again using (<ref>) together with the limitation on the
maximum step size, we see that

    v^n≤(
    G_n+∑_j=1^n-1P^(n)_n-j∑_k=1^jλ_j-kv^k-θ)
            +v^n/2,

which is equivalent to (<ref>) so the remainder of the proof is
unchanged.



The discrete fractional Grönwall inequalities in
<ref> are valid on very general nonuniform time meshes
and differ substantially from the discrete fractional Grönwall inequality
of Jin et al. <cit.>, which is built on the uniform mesh for both the L1 scheme
and the convolution quadratures generated by backward difference formulas.



Note that our theory starts only from the discrete convolution form
(<ref>) and the three assumptions 1–3,
but not the continuous counterpart (<ref>).
Correspondingly, the complementary discrete kernels P^(n)_n-j
defined in (<ref>) are also independent of (<ref>).
In other words, the fractional order α of Caputo's
derivative α v in <ref>,
and the fractional exponent α in the Mittag–Leffler function E_α
in <ref>, are
determined only by the integrand function ω_1-α(t_n-s) of the
lower bound in 2, but are independent of
the continuous counterpart of (<ref>).

To explain this point more clearly, suppose that the discrete convolution
form (<ref>) arises from some numerical formula
for a multi-term Caputo derivative
∑_i=1^mw_iα_i v with 0<α_i<1 and the weights w_i>0,
see <cit.>. Then all of the fractional exponents α_i or
the maximum order max_1≤ i≤ mα_i can
determine a single fractional exponent α for 2 and
the Mittag–Leffler function E_α in <ref>.
Hence, the presented results would be also useful for studying numerical
approximations of multi-term Caputo derivatives and distributed-order Caputo
derivatives, since the latter can be approximated by certain multi-term
derivatives via a proper quadrature rule <cit.>.





There are other practically important formulas, such as the Caputo BDF2-like
approach <cit.>.
To start the time-stepping process, one computes the first-level solution
by the L1 approach in <ref>, (αv)^1:=a^(1)_0
v^1,
or the Alikhanov formula in <ref>,
(αv)^1:=â^(1)_0 v^1.
For any time-level t_n with n≥2, taking θ=0 and
applying the quadratic polynomial interpolation Π_2,kv,
we have a Caputo BDF2-like formula <cit.>

    (αv)^n    ∑_k=1^n-1∫_t_k-1^t_kω_1-α(t_n-s)⟨Π_2,kv|'(s)s
        +∫_t_n-1^t_nω_1-α(t_n-s)⟨Π_2,n-1v|'(s)s  for n≥2.

One can obtain the compact form (<ref>)
with the discrete kernels A_n-k^(n),

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

where the coefficients a_n-k^(n) are defined in <ref>,
and the b_n-k^(n) are defined by

    b^(n)_02/τ_n-1(τ_n-1+τ_n)∫_t_n-1^t_n
    	(s-t_n-1/2)ω_1-α(t_n-s)s,
       b^(n)_n-k2/τ_k(τ_k+τ_k+1)∫_t_k-1^t_k
    	(s-t_k-1/2)ω_1-α(t_n-s)s  for 1≤ k≤ n-1.


Notice that if the fractional order α→1, then
ω_3-α(t)→ t, ω_2-α(t)→1
and ω_1-α(t)→0, uniformly for t>0.  Thus, we have
a^(n)_0=ω_2-α(τ_n)/τ_n→1/τ_n and

    b^(n)_0=2/τ_n-1(τ_n-1+τ_n)ω_3-α(τ_n)-τ_n/2ω_2-α(τ_n)→τ_n/τ_n-1(τ_n-1+τ_n),

whereas a^(n)_n-k→0 and b^(n)_n-k→0 for k≥ 1.
So, when the fractional order α→1,

    (αv)^n→ D_2v^n:=1/τ_n+1/τ_n-1+τ_n v^n
    -τ_n/τ_n-1(τ_n-1+τ_n) v^n-1

which is the second-order BDF2 scheme for the classical diffusion equations.
We see that the second kernel A_1^(n) can be negative, at least,
when α is close to 1 (whereas the Caupto BDF2 scheme is shown in <cit.>
to preserve the discrete maximum principle and nonnegativity property when
α is close to 0).



The Caputo BDF2 formula may not meet
our a priori assumptions 1–2, which results in that our Grönwall inequality
would be not applicable directly. It is not surprising because, for a classical
parabolic equation, the standard discrete Grönwall inequality can also not
be applied directly to the second-order BDF2 scheme.
However, a weighted recombination technique works well;
see the detailed analysis by Thomée <cit.>
for a uniform time mesh, and a similar technique for nonuniform meshes
<cit.>. For the Caputo BDF2 formula, <ref> would be also useful for the stability and convergence
analysis if it can be rearranged to meet the positive and monotone assumptions
1–2. On the uniform mesh with τ_n=τ, Lv and Xu
<cit.> developed a new technique of variable-weights recombination and
achieved a new form of (αv)^n with a new variable
v̅^k:=v^k-η v^k-1 and v̅^0:=v^0; in our notations,

    (αv)^n=∑_k=1^n
    	A̅^(n)_n-kv̅^k+v^0∑_j=1^nA^(n)_n-jη^j

where the combination parameter η:=1/21-A^(n)_1/A^(n)_0.
From the substitution formulas

    v^k
    =∑_ℓ=0^kη^k-ℓv̅^ℓ  and   v^k=∑_ℓ=1^kη^k-ℓv̅^ℓ+η^kv^0,

one has a new series of discrete convolution weights

    A̅^(n)_n-k:=∑_j=k^nA^(n)_n-jη^j-k  for 1≤ k≤ n.

The results of <cit.> imply that 0<η<2/3
and the new convolution kernels A̅^(n)_n-k are positive and monotone,

    A̅^(n)_0>A̅^(n)_1>⋯>A̅^(n)_n-1>0  for 1≤ k≤ n.

Thus, our discrete Grönwall inequalities
 (and the complementary discrete convolution kernels P^(n)_n-j as well)
could be applied  for this new form (<ref>)
directly once a proper constant π_A in 2 is determined by a more careful examination.


Nonetheless, we do not know whether the variable-weights recombination
technique <cit.> works on nonuniform time grids. More precisely, it
has yet to be determined what constraints must be imposed on a nonuniform
mesh so that the new discrete form (<ref>) satisfies
the a priori assumptions 1-3 required by
<ref>. This problem could be very challenging, at
least technically, and remains open to us.









§ STABILITY AND CONSISTENCY

We will now outline how the results of Section <ref>
can be applied to study a numerical solution of problem (<ref>).
For simplicity, we restrict our attention to the case of a linear reaction term
f(x,t,u)κ u+ψ(x,t) with a constant κ≥0.
By applying the first Green identity, the
fractional PDE (<ref>) is written in a weak form as

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

where u,v denotes the inner product in L_2(Ω), and
(u,v)= u,v is the bilinear form induced by the elliptic
operator . Since the latter is strongly elliptic, by increasing
κ if necessary, we may assume that the bilinear form is coercive: there is a constant c>0 such that

    (v,v)≥ cv_H^1_0(Ω)^2  for all v∈ H^1_0(Ω).


Let X_h be a finite dimensional subspace of H^1_0(Ω); for example,
a (conforming) finite element space based on a triangulation of Ω
with the mesh size h.  Galerkin's method yields a spatially-discrete approximate
solution u_h:[0,T]→ X_h satisfying

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

with u_h(0)=u_h0≈ u_0 for a suitable u_h0∈ X_h.
To compute a fully-discrete numerical solution u^n_h∈ X_h, where
u(t_n)≈ u^n_h for 1≤ n≤ N, we apply the approximation (<ref>) to the fractional derivative term
in (<ref>) so that

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

for all χ∈ X_h and for 1≤ n≤ N.

The next lemma is a
discrete analogue of the inequality <cit.>

    αv^2(t)≤2(α v)(t),v(t)  for 0≤ t≤ T and 0<α<1,

and helps set the stage for applying our discrete fractional Grönwall
inequality.


Let the assumption 1 hold and  fix the parameter θ∈[0,1).  Then
every sequence (v^n)_n=0^N in L_2(Ω) satisfies

    ∑_k=1^nA^(n)_n-kv^k^2≤ 2(αv)^n-θ,v^n-θ-
    	d_nθ^(n)-θ(αv)^n-θ^2,

for 1≤ n≤ N, where 0<d_n<1/A^(n)_0 and 0<θ^(n)<1/2 are given
by

    d_n2A^(n)_0-A^(n)_1/A^(n)_0(A^(n)_0-A^(n)_1)>0
      and  θ^(n)A^(n)_0-A^(n)_1/2A^(n)_0-A^(n)_1<1/2.



By <ref> (see appendix <ref>),

    2(αv)^n-θ,v^n-θ
        =2θαv)^n-θ,v^n-1
        +2(1-θ)αv)^n-θ,v^n
    ≥∑_k=1^nA^(n)_n-k⟨v^k^2-v^k-1^2|+
    	(1-θ/A^(n)_0-θ/A^(n)_0-A^(n)_1)
    	(αv)^n-θ^2,

and the second term on the right side equals
d_n(θ^(n)-θ)(αv)^n-θ^2.



Let the assumption 1 hold and 0≤θ≤θ^(n)
for 1≤ n≤ N.  Then the fully-discrete solution u^n_h∈ X_h, defined
by (<ref>), satisfies

    ∑_k=1^nA^(n)_n-ku^k_h^2≤2κu^n-θ_h^2
    	+2u^n-θ_hψ(t_n-θ)  for 1≤ n≤ N.



Put χ=2u_h^n-θ in the Galerkin discrete equation (<ref>),
apply <ref> with v^n=u^n_h, and use
positive-definiteness (<ref>) of the bilinear form.


Applying the discrete fractional Grönwall inequality from
<ref> with

    v^nu^n_h,  
    g^n2ψ(t_n-θ),  λ_0 2κ  and  λ_j0 for 1≤ j≤ N-1,

we see from <ref> that the
scheme (<ref>) is stable in L_2(Ω),

    u^n_h≤2E_α(4max(1,)π_Aκ t_n^α)
    (u_0h
    	+2max_1≤ k≤ n∑_j=1^kP^(k)_k-jψ(t_j-θ)),

provided θ≤θ^(n) for 1≤ n≤ N. The inequality from
Remark <ref> yields a weaker but simpler stability estimate,

    u^n_h≤2E_α(4max(1,)π_Aκ t_n^α)
    (u_0h+2π_AΓ(1-α)max_1≤ k≤ nt_k^αψ(t_k-θ)).


To bound the error in u^n_h, we introduce the Ritz
projector R_h:H^1_0(Ω)→ X_h, which is well-defined by

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

because the bilinear form satisfies (<ref>).
Put 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 error in the Ritz projection R_hu^n is estimated in the usual way from a
study of the elliptic problem, so it suffices to deal with e^n_h.
Using the weak form (<ref>) at t=t_n-θ,
with v=χ, we see that

    (α u)(t_n-θ),χ+(u(t_n-θ),χ)
        =κu(t_n-θ),χ+ψ(t_n-θ),χ.








It follows from (<ref>) that

    (αe_h)^n-θ,χ+(e_h^n-θ,χ)
           =κu_h^n-θ,χ+ψ(t_n-θ),χ
           -(αR_hu)^n-θ,χ
        -(R_hu^n-θ,χ).

Therefore, since (<ref>) and (<ref>) imply

    (R_hu^n-θ,χ)=    (u^n-θ-u(t_n-θ),χ)
    	+(u(t_n-θ),χ)
    
        =    -(u^n-θ-u(t_n-θ)),χ
    	+κu(t_n-θ),χ
        +ψ(t_n-θ),χ
            -(α u)(t_n-θ),χ,

we have

    (α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 χ=2e^n-θ_h and arguing as before,
but now with v^n:=e^n_h and
g^n:=2ℛ^n, we see that (for appropriate θ)

    e^n_h≤2E_α(4max(1,)π_Aκ t_n^α)
    	(u_0h-u_0
        +2max_1≤ k≤ n∑_j=1^kP^(k)_k-jℛ^j)

for 1≤ n≤ N.  A complete error analysis would typically proceed by
applying the triangle inequality to obtain

    ℛ^j   ≤(α u)(t_j-θ)-(αu)^j-θ
    	+(α(u-R_hu))^j-θ
           +κ(u-R_hu)^j-θ
    	+(κ+)(u^j-θ-u(t_j-θ)),

and estimating separately the resulting convolutional sums over j,
refer to a new technique of global consistency error analysis
developed in recent works <cit.>.
The details would depend on the choice of the discrete kernels A^(n)_n-j and
of the space X_h, and would rely on some a priori estimates for the partial derivatives of u.

A similar approach
works if finite differences are used for the space
discretization <cit.>, by introducing an appropriate discrete
ℓ_2 inner product in place of the inner product u,v.




§ ACKNOWLEDGEMENTS

Hong-lin Liao and Jiwei Zhang would like to thank Prof. Ying Zhao,  Prof. Weiwei Sun,
Prof. Martin Stynes and Dr. Yonggui Yan for their valuable discussions and fruitful suggestions.
Hong-lin Liao thanks for the hospitality of Beijing
Computational Science Research Center during the period of his visit.







§ TWO TECHNICAL INEQUALITIES

The proof of <ref> relies on the following result,
essentially due to Alikhanov <cit.>.


If the assumption 1 holds, then every sequence (v^n)_n=0^N
in L_2(Ω) satisfies

    2(αv)^n-θ,v^n≥∑_k=1^nA^(n)_n-kv^k^2-v^k-1^2
    	+(αv)^n-θ^2/A^(n)_0

and

    2(αv)^n-θ,v^n-1[t]≥∑_k=1^nA^(n)_n-kv^k^2-v^k-1^2
    	-(αv)^n-θ^2/A^(n)_0-A^(n)_1

for 1≤ n≤ N, provided we set A^(1)_1=0 in the case n=1.


Fix n and consider the difference

    J_n:=2(αv)^n-θ,v^n
    	-∑_k=1^nA^(n)_n-kv^k^2-v^k-1^2.

We have

    J_n   =∑_k=1^nA^(n)_n-k(2v^k-v^k-1,v^n
    		-v^k-v^k-1,v^k+v^k-1)
       =∑_k=1^nA^(n)_n-kv^k-v^k-1,2v^n-(v^k+v^k-1)

and, using the identity
2v^n-(v^k+v^k-1)=v^k-v^k-1+2∑_j=k+1^n(v^j-v^j-1),

    J_n   =∑_k=1^nA^(n)_n-kv^k-v^k-1^2
    	+2∑_k=1^nA^(n)_n-k∑_j=k+1^n
    		v^k-v^k-1,v^j-v^j-1
       =∑_k=1^nA^(n)_n-kv^k-v^k-1^2
    	+2∑_j=2^n∑_k=1^j-1A^(n)_n-kv^k-v^k-1,v^j-v^j-1.

To continue the proof, it is convenient to introduce

    w^j:=∑_k=1^jA^(n)_n-k(v^k-v^k-1)
    	  and   Q_j:=1/A^(n)_n-j  for 1≤ j≤ n.

Notice that v^j-v^j-1=Q_j(w^j-w^j-1) for 2≤ j≤ n, and that
the assumption 1 implies Q_1≥ Q_2≥⋯≥ Q_n.
Thus, one deduces that

    J_n   =Q_1w^1^2+∑_j=2^nQ_jw^j-w^j-1^2
    	+2∑_j=2^nQ_jw^j-1,w^j-w^j-1
       =Q_1w^1^2+∑_j=2^nQ_jw^j^2
            -w^j-1^2
       =Q_nw^n^2+∑_j=1^n-1(Q_j-Q_j+1)
            w^j^2≥ Q_nw^n^2.

The first inequality now follows by noting that w^n=(αv)^n-θ and
Q_n=1/A^(n)_0. Furthermore, by using the identity
v^n-1=v^n-(v^n-v^n-1)=v^n-Q_n(w^n-w^n-1), we have

    2⟨(αv)^n-θ,   v^n-1⟩
    	-∑_k=1^nA^(n)_n-kv^k^2-v^k-1^2
    	=J_n-2Q_nw^n,w^n-w^n-1
       ≥ Q_nw^n^2
            +(Q_n-1-Q_n)w^n-1^2-2Q_nw^n,w^n-w^n-1
       =-Q_nw^n^2+2Q_nw^n,w^n-1
            +(Q_n-1-Q_n)w^n-1^2
       =1/Q_n-1-Q_nQ_nw^n+(Q_n-1-Q_n) w^n-1^2
    		-Q_nQ_n-1w^n^2
       ≥-Q_nQ_n-1/Q_n-1-Q_nw^n^2
        =-w^n^2/A^(n)_0-A^(n)_1.

Therefore the claimed second inequality follows and the proof is complete.




10

Alikhanov2010
A. A. Alikhanov, A priori estimates for solutions of boundary value
  problems for fractional-order equations, Differential Equations, 46 (2010),
  pp. 660–666. 

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

Becker:1998
J. Becker, A second order backward difference method with
variable steps for a parabolic problem,
BIT, 38(4) (1998), pp. 644–662.


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

DixonMcKee1986
J. Dixon and S. McKee, Weakly singular Gronwall inequalities,
  ZAMM Z. Angew. Math. Mech., 66 (1986).

Emmrich:2005
E. Emmrich, Convergence of the variable two-step BDF time
discretisation of nonlinear evolution problems governed by a monotone potential operator,
BIT Numer. Math., 49 (2009), pp. 297–323.


GaoSunZhang:2014
G. Gao, Z. Sun and H. Zhang, A new fractional numerical 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.
  

HardyLittlewoodPolya1934
G. H. Hardy, J. E. Littlewood, and G. Pólya, Inequalities,
  Cambridge University Press, 1934.

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


JinLiZhou2018
B. Jin, B. Li, and Z. Zhou,
Numerical analysis of nonlinear subdiffusion equations,
SIAM J. Numer. Anal., 56 (2018),  pp. 1–23.

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


LiaoLyuVongZhao:2017
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.


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


LiaoMcLeanZhang:2018
H.-L. Liao, W. McLean, J. Zhang, A second-order scheme
  with nonuniform time steps for a linear reaction-subdiffusion equation,
  arXiv:1803.09873v2, 2018. Submitted to this journal.

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

LinXu2007
Y. Lin and C. Xu, Finite difference/spectral approximations for the
  time-fractional diffusion equation, J. Comput. Phys., 225
  (2007), pp. 1533–1552.



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),
A2699–A2724.

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

OldhamSpanier1974
K. B. Oldham and J. Spanier, The Fractional Calculus, Academic
  Press, 1974.

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

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.
  

SunWu2006
Z. Sun and X. Wu, A fully discrete difference scheme for a
  diffusion-wave system, Appl. Numer. Math., 56 (2006), pp. 193–209.
  

Thomee:2006
V. Thomée, Galerkin Finite Element Methods for Parabolic Problems,
Second Edition, Springer-Verlag, Berlin Heidelberg,  2006.

YanKhanFord:2018
Y.  Yan, M. Khan and N. J. Ford, An analysis of the modified L1 scheme for time fractional
partial differential equations with nonsmooth data, SIAM J. Numer.
Anal., 56 (2018), pp. 210–227.




