Devops Kubebuilder Controllers

Devops Kubebuilder Controllers

Devops Automation Of Software Deployment Pipelin · ML Project

Python · Data Preprocessing · Model Development · Evaluation

Project focus: anomaly detection / classification using IoT sensor streams and timestamped device measurements.

A DISCRETE GRÖNWALL INEQUALITY WITH APPLICATION TO

NUMERICAL SCHEMES FOR SUBDIFFUSION PROBLEMS∗

HONG-LIN LIAO† , WILLIAM MCLEAN‡ , AND JIWEI ZHANG§

Abstract. 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 L approximation to the Caputo fractional derivative, and will facilitate the analysis of higher order and linearized fast schemes.

Key words. fractional subdiffusion equations, nonuniform time mesh, discrete Caputo deriva-

tive, discrete Grönwall inequality.

AMS subject classifications. 65M06, 35B

1. Introduction. This paper builds on earlier results for the nonuniform

L method applied to the time discretization of a fractional reaction-subdiffusion

problem in a spatial domain Ω, Dtα u + Lu = f (x, t, u) for x ∈ Ω and 0 < t ≤ T , (1.1) u = u (x) for x ∈ Ω when t = 0, u=0 for x ∈ ∂Ω and 0 < t < T .

Here, Dtα = C α

0 Dt denotes the Caputo fractional derivative of order α with respect

to time t, with 0 < α < 1, and L 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 and linearized fast algorithms for solving (1.1) that employ nonuniform step sizes.

In any numerical methods for solving the reaction-subdiffusion problem (1.1), 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 (1.1) is linear and homogeneous, let ϕL be a Dirichlet eigenfunction of L on Ω, with eigenvalue λL > 0, so that LϕL = λL ϕL . Let Eα denote the Mittag–Leffler function, ∞

X zk

(1.2) Eα (z) := ,

Γ(1 + kα)

and choose as the initial data u (x) = ϕL (x). Term-by-term differentiation shows that the solution is u(x, t) = Eα (−λL tα )ϕL (x), and so ∂u/∂t = O(tα−1 ) as t → 0, ∗ Submitted to the editors DATE.

Funding: This work was funded by a grant 1008-56SYAH180 from NUAA Scientific Research

Starting Fund of Introduced Talent and a grant DRA20155 from 3 High-level Personal Train-

ing Project of Jiangsu Province; Australian Research Council grant DP140101193; NSFC grants

11771035, 91430216, U1530401.

† Department of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing,

211106, P. R. China. ([email protected]). ‡ School of Mathematics and Statistics, University of New South Wales, Sydney 2052, Australia.

([email protected]). § Beijing Computational Science Research Center, Beijing, 100094, P. R. China. ([email protected]).

whereas the solution of the classical parabolic equation, u(x, t) = e−λL t ϕL (x), is a smooth function of t. Sakamoto and Yamamoto study the (lack of) regularity of u for more general initial data u 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 .

The nature of polynomial interpolation means that the convergence rate of the L

or similar approximations to Dtα 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

(1.3) tn := (n/N )γ T for 0 ≤ n ≤ N ,

where the parameter γ ≥ 1 must be adapted to the strength of the singularity. Choos- ing γ = 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 and Volterra integral equations, and their use for time-fractional PDEs is now well established. Early papers on L schemes [17, 25] 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 presented a new analysis, based on generating functions, that

permitted nonsmooth initial data u . They showed that if f ≡ 0 and u ∈ L (Ω), then the error in the norm of L (Ω) due to the time discretization is O(τ t−1 n ). Thus, for tn bounded away from zero, the method achieves first-order accuracy in time. Yan, Khan and Ford proposed a modified L scheme and obtained error estimates for smooth and nonsmooth initial data. It was shown that the modified L scheme on a uniform mesh has a convergence rate of O(τ 2−α ). Aliknanov introduced the L2-1σ formula,

a modification of the L method that uses piecewise-quadratic instead of piecewise- linear interpolation, and approximates Dtα u at an offset grid point tj+σ = (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 considered the L method on a graded mesh of the form (1.3) applied to (1.1) for the case Lu = −uxx 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 to the present work by exploiting a novel discrete fractional Grönwall inequality for the L method.

Nonetheless, practical applications of the discrete Grönwall inequality in its basic

form 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 (n) (n) L kernels an−k and their complementary discrete kernels Pn−k , with a key step [14, Lemma 2.1] employing rough estimates of the truncation error that, to a large extent,

(n) rely on the simple form of the an−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 L 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, section 2 defines a discrete fractional derivative (2.2) having the form of the classical L approximation but with general discrete kernels. We formulate three assumptions required for our theory. The first two impose a monotonicity property (A1) and a lower bound (A2) on the discrete kernels, and the third (A3) 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 . Lemma 2.3 establishes a key estimate involving the discrete kernels and the Mittag–Leffler function (1.2). In section 3, we prove our main result, a discrete fractional Grönwall inequality stated as Theorem 3.1, and provide, in Remark 6, a strategy to treat cases where the monotonicity assumption breaks down. Section 4 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 section 4.

The generalized results proved below will allow us to show, in two companion

papers [15, 16], 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.

2. Discrete fractional derivative. Recall that the Riemann–Liouville frac-

tional integral operator of order β > 0 is defined by [20, 21]

Z t

tβ−1 (I β v)(t) := ωβ (t − s)v(s) ds for t > 0, where ωβ (t) := , 0 Γ(β) and, in turn, the Caputo fractional derivative is defined by

Z t

(2.1) (Dtα v)(t) := (I 1−α v 0 )(t) = ω1−α (t − s)v 0 (s) ds for t > 0.

For (possibly nonuniform) time levels 0 = t < t < t < · · · < tN = T , we denote the nth step size by τn := tn − tn−1 , fix an offset parameter θ ∈ [0, 1) and define

tn−θ := θtn−1 + (1 − θ)tn and v n−θ := θv n−1 + (1 − θ)v n ,

where v k may be any sequence. Letting v k ≈ v(tk ) and Oτ v k := v k − v k−1 , we consider a discrete Caputo derivative (not necessarily a direct approximation of (2.1), see Remark 5) given by a convolution-like sum, as follows, n (n) X (2.2) (Dτα v)n−θ := An−k Oτ v k for 1 ≤ n ≤ N . k=1

Here, the corresponding discrete convolution kernels are written as An−k instead

of Ank 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, (n) (n) (n) (n) A ≥ A ≥ A ≥ · · · An−1 > 0 for 1 ≤ n ≤ N .

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

Z tk

(n) 1 An−k ≥ ω1−α (tn − s) ds for 1 ≤ k ≤ n ≤ N . πA τk tk−1

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

The boundedness and monotonicity assumptions A and A on the discrete con-

(n) volution kernels An−k are valid for several frequently-used discrete Caputo derivatives, at least if assumption A is satisfied for appropriate ρ. Included are the well-known L formula [14, 17, 20, 24, 25], the fast L formula , the Alikhanov approxima- tion [2, 12, 15], and their applications for multi-term and distributed-order Caputo derivatives (see Remark 5). Here we list three examples on nonuniform grids. Note that, the local mesh parameter ρ from A will always appear in our discrete fractional

Grönwall inequality and our stability estimates.

Example 1 (nonuniform L formula). The widespread L formula [20, p. 140]

uses θ = 0 and v 0 (s) ≈ Oτ v k /τk (linear interpolation) to obtain n Z tk

X (n) (n) 1

(Dτα v)n := an−k Oτ v k with an−k := ω1−α (tn − s) ds. τk tk−1 k=1

This sum has the desired form (2.2) where (using the integral mean value theorem), (n) (n) (2.3) An−k := an−k = ω1−α (tn − snk ) for some snk ∈ [tk−1 , tk ].

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

Example 2 (fast L formula). In the two-level fast L approximation , 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 Nq , positive quadrature nodes θ` and positive weights $` (1 ≤ ` ≤ Nq ) such that Nq

X `

ω1−α (tk − s) − $` e−θ (tk −s) ≤  ∀ tk ∈ [s + ∆t, T ]. `=1

Then we use θ = 0 and v 0 (s) ≈ Oτ v k /τk (linear interpolation) to obtain Nq (n)

X `

(Dfα u)n := a ∇τ un + $` e−θ τn H ` (tn−1 ), n ≥ 1, `=1

where H ` (tk ) satisfies H ` (t ) = 0 and the recurrence relationship

Z tk

` 1 ` H ` (tk ) = e−θ τk H ` (tk−1 ) + e−θ (tk −s) ∇τ uk ds, k ≥ 1, 1 ≤ ` ≤ Nq . τk tk−1

This approximation also has the form (2.2) with θ = 0, Nq

Z tk X

(n) (n) (n) 1 ` A := a and An−k := $` e−θ (tn −s) ds for 1 ≤ k ≤ n − 1. τk tk−1 `=1

If the tolerance error  is small enough such that  ≤ min 1 ω1−α (T ), α ω2−α (1) , 

then [16, Lemma 2.5] ensures that A1-A hold true with πA = 3/2.

Example 3 (nonuniform Alikhanov formula). Let Π1,k v be the linear interpolant

of a function v with respect to the nodes tk−1 and tk , and let Π2,k v denote the quadratic interpolant with respect to tk−1 , tk and tk+1 . Taking a special choice θ = α/2, and applying the linear and quadratic polynomial interpolations, we have the nonuniform

Alikhanov formula [12, 15]

X Z tk 0

(Dτα v)n−θ := ω1−α (tn−θ − s) (Π2,k v) (s) ds k=1 tk−1

Z tn−θ

+ ω1−α (tn−θ − s) (Π1,n v) (s) ds for n ≥ 1. tn−1

(1) (1) This formula can be written as the form (2.2) with A := â0 for n = 1 and, for n ≥ 2,  (n) (n) â0 + ρn−1 b̂1 ,  for k = n, (n)

An−k := â(n)

n−k + ρ b̂ (n) k−1 n−k+1 − b̂ (n) n−k , for 2 ≤ k ≤ n − 1,  (n)  (n) ân−1 − b̂n−1 , for k = 1, (n) (n) where the discrete coefficients ân−k and b̂n−k are defined by

Z tn−θ

(n) 1 â0 := ω1−α (tn−θ − s) ds, τn tn−1

Z tk

(n) 1 ân−k := ω1−α (tn−θ − s) ds for 1 ≤ k ≤ n − 1, τk tk−1

Z tk

(n) 2 b̂n−k := (s − tk− 1 )ω1−α (tn−θ − s) ds for 1 ≤ k ≤ n − 1. τk (τk + τk+1 ) tk−1

The theoretical properties in [15, Theorem 2.2] assure A1–A with πA = 11/4 pro-

vided the local mesh assumption A 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, I α I β = I α+β , holds because the integral kernels satisfy ωα ∗ ωβ = ωα+β . It follows that

Z t

(2.4) ωα (t − µ)ω1−α (µ − s) dµ = ω (t − s) = 1 for all 0 < s < t < ∞, s

(n) and it motives us to seek a family of complementary discrete convolution kernels Pn−j having the identical property n (n) (j) X (2.5) Pn−j Aj−m ≡ 1 for 1 ≤ m ≤ n ≤ N . j=m

In fact, by taking m = k and m = k + 1, n n (n) (k) (n) (j) (n) (j)

X X

Pn−k A + Pn−j Aj−k = 1 = Pn−j Aj−(k+1) , 1 ≤ k ≤ n − 1, j=k+1 j=k+1

we see that n (n) 1 X (n)  (j) (j) 

Pn−k = (k)

Pn−j Aj−k−1 − Aj−k , 1 ≤ k ≤ n − 1,

A j=k+1

and the complementary discrete kernels may be defined via the recursion j−1  (n) 1 (n) 1 X (n−k) (n−k)  (n) (2.6) P := (n) , Pj := (n−j) Aj−k−1 − Aj−k Pk for 1 ≤ j ≤ n − 1.

A A k=0

Example 4 (Pictures of Aj and Pj of L formula). Consider the widespread

(n) L approximation in Example 1. Figure 2.1 plots the L discrete kernels Aj and the (n) complementary discrete kernels Pj when T = 1 and n = 3 for three graded meshes of the form (1.3).

As a consequence of the identity (2.4), we find that

Z t Z t

ωα (t − s)(Dtα v)(s) ds = v 0 (s) ds, 0 0

which provides the inspiration for the second part of the next lemma. Lemma 2.1. Let the assumptions A and A hold. (n)

1. The discrete kernels Pj in (2.6) having the property (2.5) satisfy

(n) 0 ≤ Pn−j ≤ πA Γ(2 − α)τnα for 1 ≤ j ≤ n ≤ N ,

and n (n) X (2.7) Pn−j ω1−α (tj ) ≤ πA for 1 ≤ n ≤ N . j=1

2. If v : [0, T ] → R is any continuous, piecewise-C 1 function such that v 0 is

non-negative and monotone decreasing, then n Z tn (n) X (2.8) Pn−j (Dtα v)(tj ) ≤ πA v 0 (s) ds for 1 ≤ n ≤ N . j=1 0

(n) Proof. It follows at once from the monotonicity assumption A that A > 0 and (n−k) (n−k) (n) Aj−k−1 − Aj−k ≥ 0 for 0 ≤ k ≤ j − 1. The lower bound Pj ≥ 0 is then clear from the recursion (2.6). Since all the discrete kernels are non-negative, we have n (n) (k) (n) (j) X

Pn−k A ≤ Pn−j Aj−k = 1

Fig. 2.1. Top: the L discrete kernel (2.3) for three different meshes of the form (1.3) in the (n) case T = 1 and n = 30. Bottom: the complementary discrete kernels Pj .

3.0 γ = 1.0 γ = 2.0 2.5 γ = 3.0

Aj n

0.0 0 5 1 1 2 2 3

0.3 γ = 1.0 γ = 2.0 0.2 γ = 3.0

0.0 0 5 1 1 2 2 3

and taking n = k in the assumption A gives

Z tk

(k) 1 ω2−α (τk ) 1 A ≥ ω1−α (tk − s) ds = = , πA τk tk−1 πA τk Γ(2 − α)πA τkα

(n) so the complementary discrete convolution kernels Pn−k are well-defined and satisfy (n) (k) the upper bound Pn−k ≤ 1/A ≤ Γ(2 − α)πA τkα . Furthermore, the assumption A (j) and the identity (2.5) imply that ω1−α (tj ) ≤ πA Aj−1 and

n n (n) (n) (j)

X X

Pn−j ω1−α (tj ) ≤ πA Pn−j Aj−1 = πA for n ≥ 1, j=1 j=1

which completes the proof of part 1.

Recall Chebyshev’s sorting inequality [9, p. 168, item 236.]: if f is monotone

increasing and g is monotone decreasing on the interval [a, b], and if both functions

Z b Z b Z b

(b − a) f (s)g(s) ds ≤ f (t) dt g(s) ds. a a a

Taking [a, b] = [tk−1 , tk ], f (s) = ω1−α (tj − s) and g(s) = v 0 (s) ≥ 0, and using A2, we see that j Z tk X (Dtα v)(tj ) = ω1−α (tj − s)v 0 (s) ds k=1 tk−1 j Z tk Z tk j Z tk

X 1 X (j)

≤ ω1−α (tj − t) dt v 0 (s) ds ≤ πA Aj−k v 0 (s) ds. τk tk−1 tk−1 tk−1 k=1 k=1 (n) Thus, from the identical property (2.5) of the discrete kernels Pn−j , we conclude that n n j Z tk

(n) (n) (j)

X X X

Pn−j (Dtα v)(tj ) ≤ Pn−j πA Aj−k v 0 (s) ds

j=1 j=1 k=1 tk−1 n Z tk n n Z tk (n) (j)

X X X

= πA v (s) ds Pn−j Aj−k = πA v 0 (s) ds, k=1 tk−1 j=k k=1 tk−1

and part 2 follows. When A also holds, we have a variant of the second part of Lemma 2.1. Lemma 2.2. Let the assumptions A1–A hold. If v : [0, T ] → R is any continu- ous, piecewise-C 1 function such that v 0 is non-negative and monotone, then n−1

X (n) Z tn

Pn−j (Dtα v)(tj ) ≤ max(1, ρ)πA v 0 (s) ds for 1 ≤ n ≤ N . j=1 0

Proof. If v is non-negative and monotone decreasing, then Dtα v(tj ) ≥ 0 and the

results of Lemma 2.1 imply that

X (n) n Z tn

(n) X Pn−j (Dtα v)(tj ) ≤ Pn−j (Dtα v)(tj ) ≤ πA v 0 (s) ds . j=1 j=1 0

Otherwise, if v 0 is monotonely increasing, then

X (n) n−1

X (n) X j Z tk

Pn−j (Dtα v)(tj ) = Pn−j ω1−α (tj − s)v 0 (s) ds

j=1 j=1 k=1 tk−1

n−1 j Z tk (n)

X X

≤ Pn−j v 0 (tk ) ω1−α (tj − s) ds j=1 k=1 tk−1

n−1 j (n) (j)

X X

≤ πA Pn−j v 0 (tk )τk Aj−k j=1 k=1 n−1 n−1 n−1 (n) (j)

X X X

= πA v 0 (tk )τk Pn−j Aj−k ≤ πA v 0 (tk )τk k=1 j=k k=1 n−1

X n−1

X Z tk+1

≤ πA v 0 (tk )ρk τk+1 ≤ ρπA v 0 (s) ds, k=1 k=1 tk

and the desired estimate again holds.

We can use Lemma 2.2 to prove the following property of the Mittag–Leffler

function (1.2). Lemma 2.3. Let the assumptions A1–A hold. For any real µ > 0, n−1

X (n) Eα (µtα

Pn−j Eα (µtα

j ) ≤ πA max(1, ρ) for 1 ≤ n ≤ N . j=1 µ

Proof. The series definition (1.2) shows that

∞ ∞

X µk tkα X

Eα (µtα ) = 1 + =1+ µk vk (t),

Γ(1 + kα)

k=1 k=1

where vk (t) = ω1+kα (t) and we have vk (t) = ωkα (t) > 0 for all k ≥ 1. If 1 ≤ k ≤ 1/α, then −1 ≤ kα − 1 ≤ 0 and vk (t) = ωkα−1 (t) ≤ 0 for all t > 0. Otherwise, if k > 1/α, then kα − 1 > 0 and vk (t) > 0 for all t > 0. Thus, vk is always non-negative and monotone, so we may apply Lemma 2.2 and deduce that n−1 Z tn (n) X Pn−j (Dtα vk )(tj ) ≤ max(1, ρ)πA vk (s) ds = max(1, ρ)πA vk (tn ) for k ≥ 1. j=1 0

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

Z t

Dtα vk (t) = ω1−α (t − s)ωkα (s) ds = ω1+(k−1)α (t) = vk−1 (t) for all k ≥ 1,

we have m n−1 m (n)

X X X

k µ Pn−j vk−1 (t) ≤ max(1, ρ)πA µk vk (tn ). k=1 j=1 k=1

P∞ k

Because the series k=1 µ vk (t) is absolutely convergent and ω (t) = 1, the desired inequality follows after interchanging the sums on the left-hand side and then sending m → ∞. The proof is completed.

3. Discrete fractional Grönwall inequality. Our main result is stated in the

next theorem. The proof is similar to that of [14, Lemma 2.2], but we include it here to incorporate the nonuniform mesh parameter ρ in A3, which does not appear in discrete Grönwall inequalities for classical parabolic equations. Theorem 3.1. Let the assumptions A1–A hold, let 0 ≤ θ < 1, and let (g n )N n=1 −1 and (λl )N l=0 be given non-negative sequences. Assume further that there exists a

PN −1

constant Λ (independent of the step sizes) such that Λ ≥ l=0 λl , and that the maximum step size satisfies max τn ≤ p . 1≤n≤N α 2πA Γ(2 − α)Λ

Then, for any non-negative sequence (v k )N

n n k 2

X (n)

X
2

λn−k v k−θ + v n−θ g n
(3.1) An−k Oτ v ≤ for 1 ≤ n ≤ N , k=1 k=1

it holds that  k  (k) j X n ≤ 2Eα 2 max(1, ρ)πA Λtα 0
(3.2) v n v + max Pk−j g for 1 ≤ n ≤ N . 1≤k≤n j=1

(n) Proof. We replace the index n with j in (3.1), then multiply by Pn−j and sum over j to obtain n j n j n

X (n)

X (j)
2 X (n)

X
2 X (n)

(3.3) Pn−j Aj−k Oτ v k ≤ Pn−j λj−k v k−θ + Pn−j v j−θ g j . j=1 k=1 j=1 k=1 j=1

On the left-hand side, we exchange the order of summation and use the identity (2.5) to get n j n n

X (n)

X (j)
2 X
2 X (n) (j)

Pn−j Aj−k Oτ v k = Oτ v k Pn−j Aj−k

j=1 k=1 k=1 j=k (3.4) n

X
2

= Oτ v k = (v n )2 − (v 0 )2 . k=1

Thus, it follows from (3.3) that

n j n n 2 0 2 (n) k−θ 2 (n)

X X X

Pn−j v j−θ g j ,


(3.5) v ≤ v + Pn−j λj−k v + j=1 k=1 j=1

For brevity, let us write the claimed estimate (3.2) as v n ≤ Fn Gn where k (k) X Fn := 2Eα 2 max(1, ρ)πA Λtα and Gn := v 0 + max Pk−j g j .
n 1≤k≤n j=1

We will use complete induction, noting that the Mittag–Leffler function (1.2) satisfies

Eα (0) = 1 and Eα (z) > 0 for all real z > 0, so Fn ≥ Fn−1 ≥ 2 for n ≥ 2. If v 1 ≤ v 0 , then v 1 ≤ G ≤ F G , as required. Otherwise, if v 1 > v 0 , then 1−θ v ≤ v 1 . One deduces from (3.5) that
2
2 (1) (1)
2 v 1 ≤ v 0 + P v 1−θ g 1 + P λ v 1−θ (1)
(1)
2 (1)
2 ≤ v 1 v 0 + P g 1 + P λ v 1 = v 1 G + P λ v 1 . Part 1 of Lemma 2.1 and the given restriction on the maximum time-step imply that (1) (3.6) P λ ≤ πA Γ(2 − α)τ1α Λ ≤ 1/2.

Thus, (v 1 )2 ≤ 2v 1 G and so v 1 ≤ 2G ≤ F G , which implies that the desired estimate holds for n = 1.

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

(3.7) v k ≤ Fk Gk for 1 ≤ k ≤ n − 1. Choose some k(n) such that v k(n) = max0≤j≤n−1 v j . If v n ≤ v k(n) then, since Fk and Gk are monotone increasing in k, v n ≤ v k(n) ≤ Fk(n) Gk(n) ≤ Fn Gn ,

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 (3.5) that

n n−1 j n
2 X (n)

X (n)

X
2 (n) X

(3.8) v n ≤ v n v 0 +v n Pn−j g j +v n Pn−j λj−k v k−θ + v n P λn−k . j=1 j=1 k=1 k=1

Using part 1 of Lemma 2.1, n (n) X (3.9) P λn−k ≤ πA Γ(2 − α)Λτnα , k=1

so the limitation on the maximum step size implies that  n−1

X (n) X j  1

2 2 (3.10) (v n ) ≤ v n Gn + Pn−j λj−k v k−θ + (v n ) . j=1 k=1

Thus, applying the induction hypothesis (3.7), we deduce from (3.10) that

X X

n λj−k θv k−1 + (1 − θ)v k   v ≤ 2Gn + 2 Pn−j j=1 k=1 n−1 j (n)

X X  

≤ 2Gn + 2 Pn−j λj−k θFk−1 Gk−1 + (1 − θ)Fk Gk j=1 k=1 n−1 j n−1 j (n) (n)

X X X X

≤ 2Gn + 2 Pn−j λj−k Fk Gk ≤ 2Gn + 2 Pn−j Fj Gj λj−k j=1 k=1 j=1 k=1 n−1 (n) X

Pn−j Eα 2 max(1, ρ)πA Λtα

≤ 2Gn + 4ΛGn−1 j . j=1

Finally, by Lemma 2.3 with µ = 2 max(1, ρ)πA Λ,
n Eα 2 max(1, ρ)πA Λtα n −1 v ≤ 2Gn + 2 max(1, ρ)πA ΛGn = Fn Gn , max(1, ρ)πA Λ which completes the inductive step and the proof. Remark 1. One may use the inequality (2.7) in part 1 of Lemma 2.1 to bound

Pk (k)

the convolutional summation j=1 Pk−j g j , that is, k k

X (k)

X (k) gj gj

Pk−j g j ≤ Pk−j ω1−α (tj ) max ≤ πA max

1≤j≤k ω1−α (tj ) 1≤j≤k ω1−α (tj ) j=1 j=1

So the discrete solution of (3.1) can also be bounded by

 0  v n ≤ 2Eα 2 max(1, ρ)πA Λtα n v + π A Γ(1 − α) max {t α j j g } for 1 ≤ n ≤ N . 1≤j≤n

N −1

On the other hand, if the given sequence (λl )l=0 is non-positive and the constant Λ ≤ 0, a similar argument will show that the discrete inequality (3.2) holds in a simpler form, requiring only the assumptions A1-A but no restrictions on time steps, k (k) X (3.11) v n ≤ v 0 + max Pk−j g j ≤ v 0 + πA Γ(1 − α) max {tα j jg } for 1 ≤ n ≤ N . 1≤k≤n 1≤j≤n j=1

−1 Remark 2. By including the non-negative sequence (λl )N l=0 in (3.1), 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 section 4, and λl = 0 for l ≥ 2 in the one-step linearized scheme for a semilinear subdiffusion equation. Thus, the constant p Λ is always not very large and the maximum time-step restriction

max1≤n≤N τn ≤ 1/ α 2πA Γ(2 − α)Λ is also not stringent in practical applications.

Remark 3. The Mittag–Leffler function Eα also arises naturally in other discrete

and continuous Grönwall inequalities for fractional diffusion and wave equations [1, Lemma 2], and for weakly singular Volterra equations [5, Theorems 1.3 and 1.6].

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.

Theorem 3.2. Theorem 3.1 remains valid if the condition (3.1) is replaced by

n n (n)

X X

(3.12) An−k Oτ v k ≤ λn−k v k−θ + g n for 1 ≤ n ≤ N . k=1 k=1

Moreover, if the given sequence (λl )N

l=0 is non-positive and the constant Λ ≤ 0, n (n) X (3.13) vn ≤ v + Pn−j g j ≤ v 0 + πA Γ(1 − α) max {tα j jg } for 1 ≤ n ≤ N . 1≤j≤n j=1

Proof. The structure of proof is as before. However, instead of (3.3) and (3.4), we now have n j n j n (n) (j) (n) (n)

X X X X X

Pn−j Aj−k Oτ v k ≤ Pn−j λj−k v k−θ + Pn−j g j

j=1 k=1 j=1 k=1 j=1

and n j n n n (n) (j) (n) (j)

X X X X X

Pn−j Aj−k Oτ v k = Oτ v k Pn−j Aj−k = Oτ v k = v n − v 0 , j=1 k=1 k=1 j=k k=1

respectively, so that instead of (3.5) we obtain n j n (n) (n)

X X X

vn ≤ v + Pn−j λj−k v k−θ + Pn−j g j . j=1 k=1 j=1

As before, if v 1 ≤ v 0 then v 1 ≤ G . For the alternative case v 1 > v 0 , we again have v 1−θ ≤ v 1 which now yields (1) (1) (1) v 1 ≤ v 0 + P g 1 + P λ v 1−θ = G + P λ v 1−θ ≤ G + 1 v 1 , where the final step again relies on the step size assumption to ensure (3.6). Thus, once again, v 1 ≤ 2G . In the inductive step, (3.8) is replaced by n n−1 j n (n) (n) (n)

X X X X

vn ≤ v + Pn−j g j + Pn−j λj−k v k−θ + v n P λn−k , j=1 j=1 k=1 k=1

and by again using (3.9) together with the limitation on the maximum step size, we see that n−1 j vn  X (n) X  n k−θ v ≤ Gn + Pn−j λj−k v + , j=1 k=1

which is equivalent to (3.10) so the remainder of the proof is unchanged.

Remark 4. The discrete fractional Grönwall inequalities in Theorems 3.1 and 3.2

are valid on very general nonuniform time meshes and differ substantially from the discrete fractional Grönwall inequality of Jin et al. [11, Theorem 2.8], which is built on the uniform mesh for both the L scheme and the convolution quadratures generated by backward difference formulas.

Remark 5 (Multi-term and distributed-order Caputo derivatives). Note that

our theory starts only from the discrete convolution form (2.2) and the three assump- tions A1–A3, but not the continuous counterpart (2.1). Correspondingly, the comple- (n) mentary discrete kernels Pn−j defined in (2.5) are also independent of (2.1). In other words, the fractional order α of Caputo’s derivative Dtα v in Lemmas 2.1 and 2.2, and the fractional exponent α in the Mittag–Leffler function Eα in Lemma 2.3 and Theo- rems 3.1 and 3.2, are determined only by the integrand function ω1−α (tn − s) of the

lower bound in A2, but are independent of the continuous counterpart of (2.2).

To explain this point more clearly, suppose that the discrete convolution

Pmform (2.2)

arises from some numerical formula for a multi-term Caputo derivative i=1 wi Dtαi v with 0 < αi < 1 and the weights wi > 0, see . Then all of the fractional exponents αi or the maximum order max1≤i≤m αi can determine a single fractional exponent α for A and the Mittag–Leffler function Eα in Theorems 3.1 and 3.2. 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 .

Remark 6 (Caputo BDF2-like formula and an open problem). There are other

practically important formulas, such as the Caputo BDF2-like approach [7, 13, 18].

To start the time-stepping process, one computes the first-level solution by the L

(1) approach in Example 1, (Dτα v)1 := a Oτ v 1 , or the Alikhanov formula in Example 3, (1) (Dτα v)1 := â0 Oτ v 1 . For any time-level tn with n ≥ 2, taking θ = 0 and applying the quadratic polynomial interpolation Π2,k v, we have a Caputo BDF2-like formula

X Z tk 0

(Dτα v)n := ω1−α (tn − s) (Π2,k v) (s) ds k=1 tk−1

Z tn

+ ω1−α (tn − s) (Π2,n−1 v) (s) ds for n ≥ 2. tn−1

(n) One can obtain the compact form (2.2) with the discrete kernels An−k ,  (n) (n) (n)
  a + ρn−1 b + b , for k = n,  a(n) + ρn−2 b(n) − b(n) + b(n)
, for k = n − 1,   (n) 1 2 1 0

An−k := (n) (n) (n)

   an−k + ρ b k−1 n−k+1 − b n−k , for 2 ≤ k ≤ n − 2,   (n)  (n) an−1 − bn−1 , for k = 1.

(n) (n) where the coefficients an−k are defined in Example 1, and the bn−k are defined by

Z tn

(n) 2 b := (s − tn− 1 )ω1−α (tn − s) ds, τn−1 (τn−1 + τn ) tn−1

Z tk

(n) 2 bn−k := (s − tk− 1 )ω1−α (tn − s) ds for 1 ≤ k ≤ n − 1. τk (τk + τk+1 ) tk−1

Notice that if the fractional order α → 1, then ω3−α (t) → t, ω2−α (t) → 1 and (n) ω1−α (t) → 0, uniformly for t > 0. Thus, we have a = ω2−α (τn )/τn → 1/τn and

(n) 2 h τn i τn b = ω3−α (τn ) − ω2−α (τn ) → , τn−1 (τn−1 + τn ) 2 τn−1 (τn−1 + τn )

(n) (n) whereas an−k → 0 and bn−k → 0 for k ≥ 1. So, when the fractional order α → 1, 1 1  τn (Dτα v)n → D v n := + Oτ v n − Oτ v n−1 τn τn−1 + τn τn−1 (τn−1 + τn )

which is the second-order BDF scheme for the classical diffusion equations. We see (n) that the second kernel A can be negative, at least, when α is close to 1 (whereas the

Caupto BDF scheme is shown in to preserve the discrete maximum principle

and nonnegativity property when α is close to 0).

The Caputo BDF formula may not meet our a priori assumptions A1–A2, 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 BDF scheme. How- ever, a weighted recombination technique works well; see the detailed analysis by

Thomée [26, Theorem 1.7] for a uniform time mesh, and a similar technique for

nonuniform meshes [3, 6]. For the Caputo BDF formula, Theorems 3.1 and 3.2 would be also useful for the stability and convergence analysis if it can be rearranged to meet the positive and monotone assumptions A1–A2. On the uniform mesh with τn = τ , Lv and Xu developed a new technique of variable-weights recombination and achieved a new form of (Dτα v)n with a new variable v̄ k := v k −ηv k−1 and v̄ 0 := v 0 ; in our notations, n n (n) (n)

X X

(3.14) (Dτα v)n = Ān−k Oτ v̄ k + v 0 An−j η j k=1 j=1

(n) (n)
where the combination parameter η := 2 1 − A /A . From the substitution for- mulas k

X k

X vk = η k−` v̄ ` and Oτ v k = η k−` Oτ v̄ ` + η k v 0 , `=0 `=1

one has a new series of discrete convolution weights n (n) (n) X Ān−k := An−j η j−k for 1 ≤ k ≤ n. j=k

The results of [18, Lemma 3.2] imply that 0 < η < 2/3 and the new convolution (n) kernels Ān−k are positive and monotone, (n) (n) (n) Ā0 > Ā1 > · · · > Ān−1 > 0 for 1 ≤ k ≤ n.

Thus, our discrete Grönwall inequalities (and the complementary discrete convolution

(n) kernels Pn−j as well) could be applied for this new form (3.14) directly once a proper constant πA in A is determined by a more careful examination.

Nonetheless, we do not know whether the variable-weights recombination tech-

nique 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 (3.14) satisfies the a priori assumptions A1-A required by Theorems 3.1 and 3.2. This problem could be very challenging, at least technically, and remains open to us.

4. Stability and consistency. We will now outline how the results of Section 3

can be applied to study a numerical solution of problem (1.1). 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 (1.1) is written in a weak form as

(4.1) hDtα u, vi + B(u, v) = κ hu, vi + hψ(t), vi for all v ∈ H (Ω) and for 0 < t ≤ T ,

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

(4.2) B(v, v) ≥ ckvk2H 1 (Ω) for all v ∈ H (Ω).

Let Xh be a finite dimensional subspace of H (Ω); 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 uh : [0, T ] → Xh satisfying

(4.3) Dtα uh , χ + B(uh , χ) = κ uh , χ + ψ(t), χ for all χ ∈ Xh and 0 < t ≤ T ,

with uh (0) = uh ≈ u for a suitable uh ∈ Xh . To compute a fully-discrete numerical solution unh ∈ Xh , where u(tn ) ≈ unh for 1 ≤ n ≤ N , we apply the approximation (2.2) to the fractional derivative term in (4.3) so that

(Dτα uh )n−θ , χ + B un−θ , χ = κ uhn−θ , χ + ψ(tn−θ ), χ
(4.4) h

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

The next lemma is a discrete analogue of the inequality [1, Lemma 1]

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

and helps set the stage for applying our discrete fractional Grönwall inequality. Lemma 4.1. Let the assumption A hold and fix the parameter θ ∈ [0, 1). Then every sequence (v n )N n=0 in L (Ω) satisfies

X (n) 2

An−k Oτ kv k k ≤ 2 (Dτα v)n−θ , v n−θ − dn θ(n) − θ (Dτα v)n−θ ,

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

2A − A A − A 1

dn := (n) (n) (n) >0 and θ(n) := (n) (n) < .

A (A − A ) 2A − A 2

Proof. By Lemma A.1 (see appendix A),

2 (Dτα v)n−θ , v n−θ = 2θ Dτα v)n−θ , v n−1 + 2(1 − θ) Dτα v)n−θ , v n n  1 − θ 

X (n)

 2 2 θ 2 ≥ An−k v k − v k−1 + (n) − (n) (n) (Dτα v)n−θ , k=1 A 0 A 0 − A 1

and the second term on the right side equals dn (θ(n) − θ)k(Dτα v)n−θ k . Theorem 4.2. Let the assumption A hold and 0 ≤ θ ≤ θ(n) for 1 ≤ n ≤ N .

Then the fully-discrete solution unh ∈ Xh , defined by (4.4), satisfies

X (n) 2
2

An−k Oτ ukh ≤ 2κ un−θ

h + 2 un−θ h ψ(tn−θ ) for 1 ≤ n ≤ N . k=1

Proof. Put χ = 2un−θ h in the Galerkin discrete equation (4.4), apply Lemma 4.1 with v n = unh , and use positive-definiteness (4.2) of the bilinear form.

Applying the discrete fractional Grönwall inequality from Theorem 3.1 with

v n := unh , g n := 2 ψ(tn−θ ) , λ := 2κ and λj := 0 for 1 ≤ j ≤ N − 1,

we see from Theorem 4.2 that the scheme (4.4) is stable in L (Ω),  k  (k) X unh ≤ 2Eα 4 max(1, ρ)πA κtα
n u0h + 2 max Pk−j ψ(tj−θ ) , 1≤k≤n j=1

provided θ ≤ θ(n) for 1 ≤ n ≤ N . The inequality from Remark 1 yields a weaker but simpler stability estimate,   n α α
uh ≤ 2Eα 4 max(1, ρ)πA κtn u0h + 2πA Γ(1 − α) max tk ψ(tk−θ ) . 1≤k≤n

To bound the error in unh , we introduce the Ritz projector Rh : H (Ω) → Xh , which is well-defined by

(4.5) B(Rh v, χ) = B(v, χ) for all v ∈ H (Ω) and χ ∈ Xh ,

because the bilinear form satisfies (4.2). Put enh = unh − Rh un ∈ Xh where un = u(tn ), so that unh − un ≤ un − Rh un + enh . The error in the Ritz projection Rh un is estimated in the usual way from a study of the elliptic problem, so it suffices to deal with enh . Using the weak form (4.1) at t = tn−θ , with v = χ, we see that

(4.6) (Dtα u)(tn−θ ), χ + B(u(tn−θ ), χ) = κ u(tn−θ ), χ + ψ(tn−θ ), χ .

It follows from (4.4) that

(Dτα eh )n−θ , χ + B(en−θ h , χ) = κ uhn−θ , χ + ψ(tn−θ ), χ − (Dτα Rh u)n−θ , χ − B(Rh un−θ , χ).

Therefore, since (4.5) and (4.6) imply

B(Rh un−θ , χ) = B(un−θ − u(tn−θ ), χ) + B(u(tn−θ ), χ)

= − 4(un−θ − u(tn−θ )), χ + κ u(tn−θ ), χ + ψ(tn−θ ), χ − (Dtα u)(tn−θ ), χ , we have (Dτα eh )n−θ , χ + B(en−θ h , χ) = κ en−θ h , χ + hRn , χi for all χ ∈ Xh , where Rn = (Dtα u)(tn−θ ) − (Dτα Rh u)n−θ − κ u(tn−θ ) − Rh un−θ + 4 un−θ − u(tn−θ ) .

Choosing χ = 2en−θ

h and arguing as before, but now with v n := enh and g n := n

2 R , we see that (for appropriate θ)

 k  (k) X enh ≤ 2Eα 4 max(1, ρ)πA κtα j
n ku0h − u k + 2 max P k−j kR k 1≤k≤n j=1

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

Rj ≤ (Dtα u)(tj−θ ) − (Dτα u)j−θ + (Dτα (u − Rh u))j−θ

+ κ (u − Rh u)j−θ + (κ + 4) uj−θ − u(tj−θ ) ,

and estimating separately the resulting convolutional sums over j, refer to a new technique of global consistency error analysis developed in recent works [14–16]. The (n) details would depend on the choice of the discrete kernels An−j and of the space Xh , 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 discretiza-

tion , by introducing an appropriate discrete `2 inner product in place of the inner product hu, vi.

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.

Appendix A. Two technical inequalities. The proof of Lemma 4.1 relies

on the following result, essentially due to Alikhanov [2, Lemma 1]. Lemma A.1. If the assumption A holds, then every sequence (v n )N n=0 in L (Ω) satisfies 2 n

X (n)

 2 2  (Dτα v)n−θ 2 (Dτα v)n−θ , v n ≥ An−k v k − v k−1 + (n) k=1 A and n

X (n)   (Dτα v)n−θ

α n−θ n−1 k 2 k−1 2 2 (Dτ v) ,v ≥ An−k v − v − (n) (n) k=1 A − A (1) for 1 ≤ n ≤ N , provided we set A = 0 in the case n = 1.

Proof. Fix n and consider the difference

n  

X (n) 2 2

Jn := 2 (Dτα v)n−θ , v n − An−k vk − v k−1 . k=1

We have

n (n) X An−k 2 v k − v k−1 , v n − v k − v k−1 , v k + v k−1

Jn =

k=1 n (n) X = An−k v k − v k−1 , 2v n − (v k + v k−1 ) k=1 Pn and, using the identity 2v n − (v k + v k−1 ) = v k − v k−1 + 2 j=k+1 (v j − v j−1 ), n n n

X (n) 2 X (n)

X Jn = An−k v k − v k−1 +2 An−k v k − v k−1 , v j − v j−1 k=1 k=1 j=k+1 n j−1 n X

X (n) 2 X (n)

= An−k v k − v k−1 +2 An−k v k − v k−1 , v j − v j−1 . k=1 j=2 k=1

To continue the proof, it is convenient to introduce

X (n) 1

wj := An−k (v k − v k−1 ) and Qj := (n) for 1 ≤ j ≤ n. k=1