g L 1 Scheme for Solving a Class of Generalized Time-Fractional Diffusion Equations

: In this paper, a numerical scheme based on a general temporal mesh is constructed for a generalized time-fractional diffusion problem of order α . The main idea involves the generalized linear interpolation and so we term the numerical scheme the gL 1 scheme . The stability and convergence of the numerical scheme are analyzed using the energy method. It is proven that the temporal convergence order is ( 2 − α ) for a general temporal mesh. Simulation is carried out to verify the efﬁciency of the proposed numerical scheme.


Introduction
The fractional model has been shown to be a powerful tool in modeling various memory process in applications [1,2], such as diffusion process in a porous media. The essence of this tool lies in fractional derivative (or integral) which is an extension of integer derivative (or integral). Many well known mathematicians such as Euler, Lagrange, Liouville, Riemann, Grüwald, Letnikov and Caputo have devoted their efforts to fractional calculus and have made great contributions to this topic. The extension from integer derivative (or integral) to fractional derivative (or integral) may not be unique due to the different techniques applied. To unify the different approaches, Agrawal [3] has proposed some generalized fractional operators which may unify some well known fractional operators such as Caputo, Riemann-Louville and Hadamard. The generalized fractional operators incorporate a scale function z(t) and a weight function w(t). With these two functions, many equations can be written in a general form and thus can be solved in an elegant way as shown in [3]. Moreover, by choosing different z(t) and w(t), one can readily obtain the well known fractional operators.
The generalized fractional operators naturally lead to generalized fractional problems. As in the fractional situation, the analytical solutions may not be easily derived and hence the corresponding numerical solutions are both necessary and useful in applications. In fact, some pioneer works have been done on the numerical treatment of certain generalized fractional problems [4][5][6][7][8][9][10]. In the earlier papers [8][9][10], the convergence of the numerical scheme is established using the Lax-Richtmyer theorem but the order of convergence is not explicitly given. Inspired by [11,12], some generalized weighted shifted Grünwald-Letnikov (gWSGL) type approximations [4][5][6] and generalized Alikhanov's approximation [7] were recently proposed, which improve the accuracy of previous work. In fact, the convergence order of these methods are shown to be O(τ 2 z ) (or higher) based on a particular choice of the temporal mesh that closely depends on the scale function z(t). It is noteworthy that the higher convergence order comes at the expense of more computations and restrictions in reality.
In this paper, we shall consider a generalized time-fractional diffusion problem that is more general than that in [9]. As one of the pioneer numerical treatments for generalized time-fractional diffusion problems, the authors in [9] first present analytical solution of a generalized time-fractional problem involving second order spatial derivative. Then, based on a uniform temporal mesh, they construct a numerical scheme by finite difference method. The stability of the numerical scheme is investigated via an estimate of the inverse of the coefficient matrix, and the convergence then follows from the Lax-Richtmyer theorem without giving explicit convergence order. Different from the work of [9], we aim to derive a numerical scheme based on a more general temporal mesh than [9] and give the convergence order of the proposed scheme explicitly using energy method. Our major contributions in this paper are as follows: • We consider a problem involving an operator L and propose a numerical scheme based on a general temporal mesh using generalized linear interpolation. • We establish the stability and convergence of the proposed scheme, which is based on a general temporal mesh, via energy method. The analysis in the context of general temporal mesh is not trivial.
We consider the following generalized time-fractional diffusion equation with weight function w(t) ≡ 1 where L is a linear operator defined by and φ 1 (t), φ 2 (t), ψ(x), f (x, t) are given functions that are sufficiently smooth. We remark that (i) a generalized fractional problem of type (1) with any weight function w(t) (not necessarily 1) can be converted to (1) by a simple formula u(x, t) = w(t)v(x, t) (see [4][5][6][7]); and (ii) the generalized fractional equation considered in [9] is a special case of (1) when p(x, t) ≡ 1 and q(x, t) ≡ 0. The plan of the paper is as follows. In Section 2, we shall develop a numerical scheme for the problem (1) based on a more general temporal mesh than [9]. Then, the stability as well as the convergence of the proposed numerical scheme will be established rigorously using energy method in Section 3. In Section 4, we carry out experiments to verify as well as to demonstrate the efficiency of the proposed numerical scheme. Finally, a brief conclusion is given in Section 5.

Numerical Scheme
In this section, we shall derive a numerical scheme for the problem (1) using the key idea of generalized linear interpolation. To begin, let be any mesh in the temporal dimension and a uniform mesh in the spatial dimension with step size h = 1 M , respectively. Throughout, assume that the scale function z(t) is a strictly increasing function. Let Z : t → z(t) be a map from [0, 1] to [z(0), z(1)]. Denote z(t n ) = z n . Moreover, denote by U n j the exact solution of (1) at (x j , t n ), and by u n j an approximation of U n j .
Note that L j,k−1 (z(t)) is a generalized linear polynomial of z(t) and it will be used to approximate u(x j , t) over the interval [t k−1 , t k ] in (2). Indeed, we derive the following approximation scheme at (x j , t n ) where µ = 1 Γ(2−α) and the coefficients We shall call (4) the gL1 approximation of the generalized Caputo fractional derivative.

Remark 1.
(a) Note that (4) is derived for any temporal mesh ∆ and this is an extension of [9] which considers uniform temporal mesh. Moreover, the technique used to obtain (4) involves generalized linear interpolation whereas in [9] finite difference is used to approximate the derivatives in the integrand of the fractional derivative. (b) The Formula (4) is also different from the nonuniform L1 formula in [13] since a scale function z(t) such that z j = z i , i = j is considered.
To proceed further, we shall investigate the accuracy of the gL1 approximation (4) and the properties of the coefficients ω k z which are both vital in subsequent analysis. For the former one, we introduce the following definition. Definition 1 ([14]). Given the mesh ∆ : where τ z,max = max 1≤k≤N |z k − z k−1 |, τ z,min = min 1≤k≤N |z k − z k−1 | and ρ > 0 is a constant.

Remark 2.
There is some relation between the commonly used uniform mesh of t and the quasi uniform mesh of z(t).
(a) For finite intervals of t, say t ∈ [0, 1], in practice the partition of [0, 1] always results in finite number of subintervals, so we are able to find a constant ρ such that τ z,max τ z,min ≤ ρ. Hence, it is clear that any general mesh of t yields a quasi uniform mesh of z(t). In particular, we can say that for finite intervals, the commonly used uniform mesh of t is a special case of the quasi uniform mesh of z(t).
(b) For infinite intervals of t, in general the uniform mesh of t may not yield the quasi uniform mesh of z(t). However, if 0 < c < |z (t)| < C for all t in the infinite interval, then the uniform mesh of t will give the quasi uniform mesh of z(t), and so we can conclude here that the uniform mesh of t is a special case of the quasi uniform mesh of z(t).
Our next result gives the properties of the coefficients ω k z in (5) that is vital in subsequent analysis.
We are now ready to construct a numerical scheme for the generalized time-fractional diffusion Equation (1). Discretizing (1) at (x j , t n ) and using the approximation (4) together with finite difference method in the spatial dimension yields where µ = 1 Γ(2−α) and Λ is an operator given by (12) the gL1 scheme of the generalized time-fractional diffusion Equation (1).

Remark 3.
It is easy to verify that the coefficient matrix of the system (12) is strictly diagonally dominated. Therefore, it is uniquely solvable.

Theoretical Results
In this section, we shall analyze the stability as well as the convergence of the numerical scheme (12). To begin, let U h = {u = (u 0 , u 1 , · · · , u M )|u 0 = u M = 0} and for any u, v ∈ U h , define Obviously, · = (·, ·) is a norm defined over the space U h . Next, we shall present three lemmas that will be used later to establish the stability and convergence results. The first one gives a relation between |u| 1 and u ∞ .
The next lemma is from [16] which reveals a relation between |u| 1 and u .
Lemma 3 (Discrete Poincare inequality [16]). Suppose that u ∈ U h , then Remark 4. If h < 1 which is the case in our method, from the fact sin π 2 h ≥ h, the inequality (14) yields |u| 1 ≥ 2 u that will be used in subsequent analysis.
The last lemma gives an explicit expression of −(Λu, u).

Lemma 4. For any u
Proof. Using the definition in (13), it is found that where we have used u 0 = u M = 0 in the last equality. The result is immediate from the above equation.

Remark 5.
Since p(x, t) ≥ p 0 > 0 and q(x, t) ≥ 0, it is clear from Lemma 4 that −(Λu, u) ≥ p 0 |u| 2 1 for any u ∈ U h . Now, let us present the stability and convergence of the numerical scheme (12). To this aim, we shall first consider (12) with zero boundary conditions. Theorem 2. Let {u n j , 1 ≤ j ≤ M − 1, 1 ≤ n ≤ N} be the solution of the system (12) with zero boundary conditions. Then, we have Proof. Multiplying both sides of the first equation of (12) by −u n j and summing j from 1 to which is rearranged to Since we consider zero boundary conditions here, it is obvious that u n ∈ U h . Therefore, using Remark 5, we get a lower bound for the left side of (16) below Noting (10) in Lemma 1 and using xy ≤ 1 2 (x 2 + y 2 ), an upper bound for the first two terms on the right side of (16) is found as follows For the third term on the right side of (16), the Young's inequality gives the following upper bound Upon substituting the above upper and lower bounds into (16), we immediately get Next, using Lemma 3 and noting Remark 4, we have |u n | 1 ≥ 2 u n . Then, with = 2p 0 , the left side of (17) gives the following lower bound Substituting the above into (17) leads to Further, from (9) we see that ω n−1 Using the above inequality and µ = 1 Γ(2−α) in (18), we find where E is defined in the theorem. Noting (10), the above inequality readily leads to Now, we shall show by mathematical induction that In fact, let n = 1 in (19) and we get which implies that (20) holds for n = 1. Suppose that (20) is true up to (n − 1). Then, from (19), we have which immediately gives (20), or equivalently (15). This completes the proof.
Remark 6 (Stability). Using a similar argument as in [4][5][6][7], it can readily be deduced from Theorem 2 that the numerical scheme (12) is robust (or stable) with respect to the initial data ψ(x) and the non-homogeneous data f (x, t).
We are now ready to establish the convergence of the proposed scheme (12).

Theorem 3 (Convergence).
Assume that u(x, Z −1 (z)) =ū(x, z) ∈ C 4,2 ([0, 1] × [z(0), z(1)]). Suppose that the mesh ∆ z : z(0) = z 0 < z 1 < · · · < z N−1 < z N = z(1) is quasi uniform. Let {U n j = u(x j , t n )} be the exact solution of the problem (1), {u n j } be the numerical solution obtained from the scheme (12) and e n j = U n j − u n j be the error at (x j , t n ). Then, we have e n 2 + c(α, z)|e n | 2 Proof. Since {U n j } is the exact solution of (1), it is clear that where T n j is the local truncation error of the j-th equation. Noting (6) in Theorem 1, and using Taylor expansion at x = x j in (22), we find that Next, from (12) and (22) it is obvious that {e n j } is the solution of the system Hence, e n ∈ U h . Finally, applying Theorem 2 to (24) and noting e 0 j = 0, we obtain for 1 ≤ n ≤ N, which completes the proof.

Remark 7. From Theorem 3, it is easily seen that
Hence, the temporal convergence order in the norm · is (2 − α), which is optimal. On the other hand, the convergence order in · ∞ is not optimal. In fact, from Lemma 2, (21) and the properties of quasi uniform mesh, we get Remark 8. Theorem 3 is an extension of the work [9] in the sense that • the problem (1) involves an operator L and is more general than the problem considered in [9]; • we consider a general temporal mesh which is quasi uniform in terms of the scale function z(t), in view of Remark 2 this is more general than the uniform temporal mesh considered in [9]; • (21) and Remark 7 give the explicit convergence order of the proposed scheme (12), while convergence is proven without giving the explicit convergence order in [9].

Numerical Simulation
In this section, we shall use two examples to demonstrate the efficacy of the proposed numerical scheme (12) in the experiment to get optimal accuracy, refer to [17,18] for details); -Uniform z : uniform mesh with respect to z(t).
In view of Remark 2, we note that all the above types of temporal meshes are particular cases of quasi uniform mesh of z(t).
Clearly, the exact solution is required to compute e(N, h). When the exact solution is not available (which is commonly encountered in applications), we shall use 'approximate' exact solution, which is obtained by the numerical scheme (12) with sufficiently small mesh sizes (e.g., M = N = 2000 in our experiments), as 'exact' solution to compute errors. This is reasonable as the numerical scheme (12) is convergent.
Example 1 ([6,7,9]). Consider the generalized time-fractional diffusion equation where 0 < α < 1, z(t) is a strictly increasing scale function and Note that when z(t) = t and α = 0.85, the exact solution of Equation (26) is In this example, p(x, t) ≡ 1 and q(x, t) ≡ 0. Let us start with α = 0.85 and z(t) = t, t 0.5 , t 2 . Applying the numerical scheme (12) with fixed h = 1 512 and varied N, we compute the errors and temporal convergence orders for three types of temporal meshes-Uniform, Graded, Uniform z . The results are displayed in Table 1. From the table, it is easily seen that the experimental temporal convergence orders (≈1.15) are consistent with the theoretical ones (=2 − α) for various scale functions z(t) and different types of temporal meshes. It is a pleasant surprise that the numerical performance in terms of maximum norm is better than the theoretical result in Remark 7.
Next, we investigate the performance of the proposed numerical scheme for different values of α and scale functions z(t). Here, we use the uniform mesh of t (Uniform) and compute e N .

•
In Table 2, we fix h = 1 512 and let N vary. It is observed that the temporal convergence order is (2 − α) which is consistent with the theoretical one. • In Table 3, we fix N = 2000 and let h vary to investigate the spatial convergence. Obviously, the scheme (12) achieves second order spatial convergence order which is the same as that stated in Remark 7.  Our next example is modified from an example in [19] and it involves a general operator L.
First, consider (27) when α = 0.5. With fixed h = 1 2000 , we shall apply the numerical scheme (12) and compute the errors e N , e N ∞ and temporal convergence orders. The results are presented in Table 4, and it is clear that the numerical scheme (12) works well for different types of temporal meshes as well as for a wide range of problems (i.e., different z(t)). The experimental temporal convergence orders in · (≈1.5) agree with the theoretical ones (=2 − α), while once again it is a pleasant surprise that the numerical performance in terms of maximum norm is better than the theoretical result in Remark 7.
Next, we investigate the temporal convergence and spatial convergence of the numerical scheme (12) for different values of α and scale functions z(t). Here, we use the uniform mesh of t and compute e N . The results are presented in Tables 5 and 6. We observe that the experimental temporal/spatial convergence orders are consistent with the theoretical result.

Conclusions
In this paper, we derive a numerical scheme based on a general temporal mesh for the generalized time-fractional diffusion problem. The main idea involves the generalized linear interpolation. The stability and convergence of the proposed numerical scheme is established rigorously using energy method. More importantly, it is shown that the global convergence order is O(τ 2−α z,max + h 2 ) that extends the previous work [9]. For future work, we plan to investigate (i) the validity of the proposed scheme for nonsmooth data; (ii) high order methods based on general temporal mesh and spatial mesh for generalized fractional problems with smooth as well as nonsmooth data. We believe this will make the numerical scheme more applicable in reality.
Author Contributions: Formal analysis, X.L. and P.J.Y.W.; Investigation, X.L. and P.J.Y.W.; writingoriginal draft, X.L. and P.J.Y.W.; writing-review and editing, X.L. and P.J.Y.W. All authors have read and agreed to the published version of the manuscript.