A Discrete Grönwall Inequality and Energy Estimates in the Analysis of a Discrete Model for a Nonlinear Time-Fractional Heat Equation

: In the present work, we investigate the efﬁciency of a numerical scheme to solve a nonlinear time-fractional heat equation with sufﬁciently smooth solutions, which was previously reported in the literature [Fract. Calc. Appl. Anal. 16 : 892–910 (2013)]. In that article, the authors established the stability and consistency of the discrete model using arguments from Fourier analysis. As opposed to that work, in the present work, we use the method of energy inequalities to show that the scheme is stable and converges to the exact solution with order O ( τ 2 − α + h 4 ) , in the case that 0 < α < 1 satisﬁes 3 α ≥ 32 , which means that 0.369 (cid:47) α ≤ 1. The novelty of the present work lies in the derivation of suitable energy estimates, and a discrete fractional Grönwall inequality, which is consistent with the discrete approximation of the Caputo fractional derivative of order 0 < α < 1 used for that scheme at t k + 1/2 .


Introduction
In recent years, various works that discuss fractional heat flux models that lead to fractional heat equations have been published. For example, Compte and Metzler [1] proposed time-fractional generalized Cattaneo equations. In their models, various physical contexts were taken into account, including continuous-time random walks, non-local transport theory and delayed flux-force relations. On the other hand, a theory of thermal stresses was formulated by Povstenko [2] following [1]. Metzler and Klafterin [3] also presented a general physical introduction to fractional diffusion equations. Later on, that work motivated the use of the fractional Cattaneo model to simulate laser short-pulse heating of a solid surface [4], as well as the modeling of salt transfer processes in fractal structured media on the base of fractional derivative equations with Caputo derivatives [5].
Motivated by vast applications of fractional heat flux models, this work is devoted to the numerical solution of the fractional differential equation with initial and boundary conditions Here, K, L and T are positive constants, and f (x, t, u(x, t)) is a sufficiently smooth function that will satisfy the requirements for consistency and convergence. In the present work, the fractional derivative of order 0 < α < 1 is defined in the Caputo sense, that is, Problems of the form (1) have been investigated in the literature from the analytical point of view. As examples, the authors of [6] studied a multidimensional initial-boundary-value problem for a fractional diffusion equation with a Caputo time-fractional derivative where the coefficients depend on the spatial and temporal variables. Homogeneous Dirichlet boundary conditions were considered therein, and the authors proved the regularity and the unique existence of weak solutions. Meanwhile, a distributed-order form of the problem studied in the previous work was tackled in [7]. The existence, uniqueness and regularity of solutions were proved in that article using a classical variational approach. However, it is worth pointing out that those papers assume that the function f above depends only on the spatial and temporal variables. An interesting analytical question is whether the results of those papers can be extended to consider source functions that depend on the solution, as in the present case. In the scheme analyzed in the present work, the one-dimensional setting is used conveniently to include the dependency of the source on the solution.
The nonlinear time-fractional heat equation represents is a nonlinear integro-differential equation, which can hardly be solved analytically, whence the development of reliable numerical methods to solve it is justified. On those grounds, a widespread approximation for fractional derivatives (2) with order O(τ 2−α ) is the so-called L1 method [8,9]. In turn, the L1-2 formula [10] is a new difference analog of the Caputo fractional derivative with the order of approximation O(τ 3−α ). Both methods have been widely used for solving fractional partial differential equations with temporal Caputo derivatives [11][12][13][14]. Later on, the L2-1 σ scheme was introduced by Alikhanov [15]. This type of scheme is used to design numerical methods for different kinds of fractional equations [16][17][18]. It is worth pointing out that other approaches have been employed successfully to solve numerically time-fractional diffusive equations. As an example, some implicit-explicit (IMEX) difference methods have been proposed to solve time-fractional subdiffusion equations [19], reaction-diffusion equations [20], multi-term time-fractional diffusion models [21], multidimensional fractional diffusion equations [22] and nonlinear stiff fractional partial differential equations [23].
It is important to point out that higher order schemes for Caputo derivatives have been built in recent years based on the L1, L1-2 and L2-1 σ approaches [24]. In the case of L1 schemes, the improvements have relied on the use of piece-wise high-degree interpolating polynomials instead of linear interpolations. For the L1-2 approach, the new schemes were obtained by means of linear interpolating polynomials on the first subinterval, and quadratic polynomials on the remaining ones. However, in spite of the popularity and wide applicability of L1, L1-2 and L2-1 σ schemes, the analysis of those methods for time-fractional nonlinear problems is limited mainly due to the lack of a fundamental Grönwall-type inequality [25]. In [26], the authors established such a fundamental inequality for the L1 approximation of the Caputo derivative, and used that result to provide optimal error estimates of several fully discrete linear Galerkin finite element methods for nonlinear problems. In the present work, we will put emphasis on the one-dimensional scenario. It is worth pointing out that there are other works in the literature -such as reference [27]-, which tackle the one-, two-and three-dimensional cases.
Motivated by those and other efforts [27], we depart from a Crank-Nicholson scheme to solve linear time-fractional heat equations previously reported in the literature [28]. In that work, the authors used arguments based on Fourier analysis to prove that the proposed method is stable, and that the numerical approximations converge conditionally to the exact solution, with order O(τ 2−α + h 2 ). In our present work, we carry out the stability and the convergence analyses by using discrete energy inequalities on a compact difference scheme. We apply that compact approach Crank-Nicolson discretization proposed by [28] and establish thus the stability and the convergence properties of the scheme. The novelty of the present approach lies on the introduction of a novel discrete fractional Grönwall inequality appropriate to the modified L1 approximation formula. It is worth pointing out that some reports published in the literature for the integer case α = 1 have provided numerical schemes with orders of convergence of O(τ + h 2 ) (see reference [29]). In the present case, however, we report on a finite-difference scheme with order of convergence O(τ 2−α + h 2 ), which is temporally inferior to [29] in terms of convergence rates.

Compact Difference Scheme
In this section, we will present a Crank-Nicholson numerical model to solve (1a)-(1c). To that end, let M, N ∈ N, and define h = L/M and τ = T/N. Define also x i = ih, t k = kτ and t k+ 1 The space of all real functions defined on Ω hτ will be denoted by W.
Finally, the real number U k i will denote a numerical approximation to the exact value of u(x i , t k ), for each 0 ≤ i ≤ M and 0 ≤ k ≤ N. Obviously, U ∈ W.
Our first result yields a consistent formula to approximate second-order spatial derivatives with a high order of consistency. The proof can be found in Lemma 4.1 of Liao and Sun [30].
The following result provides a consistent formula to approximate Caputo fractional derivatives. We provide it without proof, but the derivation is summarized as Equation (2.1) of Karatay et al. [28].

Lemma 2.
If 0 < α < 1 then the following approximation formula is satisfied: where Substituting the approximation formula (7) into Equation (1a) at the point (x i , t k+ 1 2 ), we readily obtain The following result provides consistency properties of the method proposed in this work. Its proof is based on a straightforward use of Taylor's theorem.
where |R Proof. Using Taylor's expansions, the following identities hold: for each 1 ≤ i ≤ M − 1 and 0 ≤ k ≤ N − 1. Substituting these expressions into (10), using the continuity of the derivatives of f in its third component and rearranging terms, we obtain On the other hand, Lemma 1 implies that there is Applying now the operator A to (15) and using then the regularity assumption on u, we readily reach the conclusion of this result.
The finite-difference scheme to approximate the solutions of (1) is obtained by neglecting R k+ 1 2 i in the formula of Lemma 3, and replacing U k i by u k i in (12). In that way, we obtain the discrete equations with the discrete initial and boundary conditions

Auxiliary Lemmas
In the present section, we will establish some technical results to prove the stability and the convergence of the difference scheme, including a new discrete fractional Grönwall inequality which will be the cornerstone of this report. Some additional discrete nomenclature will be needed firstly.
be a grid function space on Ω h , and let u, v ∈ V h . It follows that u = (u 0 , u 1 , . . . , u M ) and v = (v 0 , v 1 , . . . , v M ). Using these conventions, we define the discrete inner product and the discrete norms and let · be the Euclidean norm induced by ·, · .
The inequalities in the next result were proved in Lemma 4.2 of [31]. They will be useful to provide bounds across discrete norms.

Lemma 4.
For any u ∈ V h the following inequalities are satisfied: The next lemma is a discrete analog of the formula for integration by parts. ).
Proof. The proof is a result of some easy (though tedious) algebraic calculations in which the formula for telescoping sums is used. We refer to [32] for a proof of this and other interesting results on summation by parts in more general settings.
Note that the time-fractional approximation (7) can be rewritten as where a α 0 = 1/2 1−α , and a α l = (l + 1 2 ) 1−α − (l − 1 2 ) 1−α for each l ≥ 1. As a consequence, we reach Here where g k+1 For convenience, we will agree that It is worthwhile to point out that the approximation of the time-fractional derivative provided in Lemma 2 coincides with that of Definition 2 (see [28]).
The proof of the following result is an easy modification of arguments used in [28]. We provide the prove herein for the sake of completeness. The result provides properties of the coefficients C k+1 m .
The following result is a straightforward consequence of Lemma 6 and Equation (28). It provides some monotonicity properties of the coefficients in (27).
and g k+1 The following result is Lemma 1 in [15] with new notation. We refer to that work for the proof.
The next result readily follows from Lemma 8.

Corollary 1.
If the hypotheses of the previous theorem hold then More technical properties are summarized in the following result. In particular, they provide some estimations required in the sequel. The result is a small modification of Lemma 3.2 in [26].
Then the following estimations are valid: (a) If 1 ≤ k ≤ j then 0 < p j < 1 and ∑ j r=k p j−r c k+1 The proof of the next result follows as in Lemma 3.3 of [26].
for each m ∈ N ∪ {0} and i ≥ j + 1. Here, The proof of the following result follows in the same lines as Lemma 3.1 of [26].

Numerical Properties
In this section, we establish mathematically the most important properties of the difference scheme. More precisely, we will prove that the discrete model has a unique solution, and that it is stable and convergent.
Proof. The solution at the time t 0 is prescribed by the initial conditions. So, let 0 ≤ k ≤ N − 1 and suppose that u k i is known, for each 0 ≤ i ≤ M. From (17a), we obtain a system of linear algebraic equations in terms of u l i . The coefficient matrix of this system is strictly diagonally dominant, so there exists a unique solution. Indeed, note that (17) can be rewritten as (57) The coefficient matrix takes the form A = (a ij ), where a ii = 10 12 and σ/2 1−α > 0. The matrix A is clearly strictly diagonally dominant.
Using the consistent approximation formula (27), the difference scheme (17) can be rewritten as with initial and boundary conditions In the following, we define ε k i = u(x i , t k ) − u k i for each 0 ≤ i ≤ M and 0 ≤ k ≤ N, and convey that Using these conventions, it is easy to check that the next identity holds: Definition 3. Let q 1 , q 2 ≥ 1. The difference model (17) is convergent with order h q 1 + τ q 2 if there exists a constant C independent of ε k i , h and ∆ such that |ε k i | ≤ C(h q 1 + τ q 2 ), for all 0 ≤ i ≤ M and 0 ≤ k ≤ N.

Theorem 3 (Convergence).
Let 0 < α < 1 be such that 3 α ≥ 3 2 , and let (u j ) N j=0 be the unique solution of the system (17). If f is Lipschitz continuous in its third component, then the method is convergent with convergence order h 4 + τ 2−α .
Proof. Multiplying each side of (63) by hε k+ 1 2 i and adding over all indexes i from 1 to M − 1, we obtain Using Corollary 1 on the left-hand side of (64), and that ε k 0 = ε k M = 0 for each k = 0, 1, . . . , N, we obtain Next, we use Lemmas 4 and 5, the Cauchy-Schwarz inequality and algebraic simplifications on the right-hand side of (64). In that way, we reach where L 1 is a Lipschitz constant for the function f in the third component. Substitute now (65) and (66) into (64), use again Lemmas 4 and 5 and employ the error of approximation to note that there exists a positive constant C 1 independent of h and τ, such that Applying now Theorem 1, we obtain Lemmas 4 and 5 yield now the conclusion of this theorem.
Next, we analyze the numerical stability of the compact difference scheme (17). More precisely, we will show that small perturbations on the initial conditions lead to small changes on the final solutions. To this end, we will suppose that v i k is the solution of subject to the initial and boundary conditions where ρ 0 i is a small perturbation of ψ(x i , t 0 ). (17) is stable if the discrete numerical solutions u k i satisfying (17), together with the solutions v k i satisfying (69), are such that

Definition 4. The numerical scheme
whereC is a nonnegative constant which is independent of h and τ.
Proof. Subtracting (69) form (60) we obtain A∆ α η as well as the conditions Multiply each side of (63) by hη k+ 1 2 i and add over all indexes i from 1 to N − 1. From this point on, we mimic the proof of Theorem 3 to complete the proof.
In the following example, we examine the stability of the scheme for various values of α. Example 1. Let L = T = 1, and consider the time-fractional partial differential equation with initial-boundary conditions It is well known that this problem has the exact solution U(x, t) = x(1 − x)(t 3 + 1) when α = 1 2 . Figure 1 shows the numerical solution of this problem when α = 1 2 , using the parameters h = 0.01 and τ = 0.001. The results exhibit the stability of the finite-difference scheme. We have used other values of the parameter α in (0, 1), though no analytical solution is known in exact form for those cases. Our results (not shown here to avoid redundancy) have shown that the solutions are stable when 0.369 α ≤ 1. If 0 < α 0.369, then the solutions are still stable despite the fact that the theoretical results are not valid for such case. Before closing this section, it is worth pointing out that the results presented herein are only theoretically valid for smooth functions. Indeed, this is a potentially strong limitation of our approach in view that some important classes of solutions to the problem (1) are excluded (see [33]). Also, our theoretical analysis uses the techniques of energy estimates employed in [28], and it is only valid for smooth functions. This approximation can be extended for non-smooth functions also [34], and its analysis could be a topic of work in the near future.

Conclusions
In this work, we departed from a difference scheme previously reported in the literature, which was used to solve a nonlinear time-fractional heat equation. The efficiency analysis carried out in that previous work relied in the use of Fourier analysis. Meanwhile, the present approach uses the method of energy inequalities to show that the scheme is stable and converges to the exact solution with order O(τ 2−α + h 4 ) if 3 α ≥ 3 2 and 0 < α < 1, that is, when 1 > α 0.369. The novelty of the present work lies in the derivation of a new discrete fractional Grönwall inequality, which is consistent with the discrete approximation of the Caputo fractional derivative of order 0 < α < 1. The numerical experiments presented in this work show that the restrictions on α may be entirely theoretical, and that a wider range of values of α may lead to stable results. To that end, the further development of the theoretical background is necessary. In that sense, the investigation of the case when the number 0 < α < 1 satisfies 3 α < 3 2 is an open problem of research. To this day, the authors have not been able to establish an affirmative solution to this case at t k+1/2 . However, we have been able to check that the arguments used in this work to establish the stability and convergence of the scheme unfortunately do not carry over to the case when α 0.369. We are convinced that new analytical tools (in the form of novel Grönwall-type inequalities) are required to that end. That may be the topic of a future work.
Author Contributions: All authors contributed equally. All authors have read and agreed to the published version of the manuscript.
Funding: The first author wishes to acknowledge the support of RFBR Grant 19-01-00019. Meanwhile, the second author would like to acknowledge the financial support of the National Council for Science and Technology of Mexico (CONACYT). The second author acknowledges financial support from CONACYT through grant A1-S-45928.