Stability and Convergence Analysis of a Class of Continuous Piecewise Polynomial Approximations for Time-Fractional Differential Equations

We propose and study a class of numerical schemes to approximate time-fractional differential equations. The methods are based on the approximations of the Caputo fractional derivative of order α∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \in (0, 1)$$\end{document} by using continuous piecewise polynomials, which are strongly related to the backward differentiation formulae. We investigate their theoretical properties, such as the local truncation error and global error estimates with respect to sufficiently smooth solutions, and the numerical stability in terms of stability region and A(π2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(\frac{\pi }{2})$$\end{document}-stability. Numerical experiments are given to verify our theoretical investigations.


Introduction
Fractional calculus, as a generalization of ordinary calculus, has been an intriguing topic for many famous mathematicians since the end of the 17th century.During the last four decades, many scholars have been working on the development of theory for fractional derivatives and integrals, found their way in the world of fractional calculus and their applications.For more detailed information on the historical background, we refer the interested reader to the following books: [39,42,38,40,24,25,6] and [23].Differential equations possessing terms with fractional derivatives in the space-or time-or space-time direction have become very important in many application areas.Particularly, in recent years a huge amount of interesting and surprising fractional models have been proposed.Here, we mention just a few typical applications: in the theory of Hankel transforms [17], in financial models [44,46], in elasticity theory [5], in medical applications [43,26], in geology [8,30], in physics [12,7,37] and many more.
Similar to the work for ordinary differential equations, that has started more than a century earlier, research on numerical methods for time fractional differential equations (tfDEs) started its development.In this paper we consider approximations to tfDEs involving Caputo fractional derivatives of order 0 < α < 1 in the form of with prescribed initial condition u(0) = u 0 .According to [14], it holds that if function f (t, u(t)) is continuous and satisfies the Lipschitz condition with respect to the second variable, the problem (1.1) then possesses a unique solution u(t) ∈ C([0, T ]).In terms of the numerical approximation of formula (1.1), we mainly aim at the numerical discretisation to the Caputo fractional derivative, the definition of which we refer to Definition 2.1 in the next section.It is observed that the Caputo fractional derivative of a well-behaved function is an operator combined with the integer-order derivative and the fractional integral, which can be regarded as a convolution of the weakly singular kernel t −β (0 < β < 1) and a function.The research on numerical approximations to fractional integral was developed in numerically solving a type of Volterra integral equation which is an equivalent form of (1.1) when f satisfies the Lipschitz continuous condition.With respect to numerical approximation for (1.2), two general approaches are proposed respectively, i.e, the product integration method and fractional linear multistep methods.In these cases, a general discretized formula to (1.2) is written as w n, j f (t j , u j ), n ≥ k. (1. 3) The fractional linear multistep methods was originally proposed in [35] in the mid eighties of the last century.This type of methods is devoted to constructing the power series generated by the convolution quadrature weights {ω j } ∞ j=0 based on classical implicit linear multistep formulae (ρ, σ ) with the following relationship For the motivation behind this idea we refer to [33].For this type of methods, the accuracy and stability properties highly benefit from those of the corresponding multistep methods [32,34].Another more straightforward approach to generate the weights {ω j } and {w n, j } is based on product quadrature method applied to the underlying fractional integral, for examples, to replace the integrand f (ξ , u(ξ )) by the piecewise Lagrange interpolation polynomials of degree k(k ≥ 0), and approximate the corresponding fractional integral in (1.2).On the accuracy and efficiency of this class of methods to a type of nonlinear Volterra integral equation with irregular kernel, we can refer to [29,13,11,16] and recently [6,27].
Other useful approaches such as collocation methods on non-smooth solution are discussed in [10,9].In addition, a series of other approaches were developed such as in [20], exponential integrators are applied to fractional order problems.Generalized Adams methods and so-called m-steps methods are utilized by [1,2].
In contrast to previous methods, the numerical approach to solve tfDEs in this paper is arrived at through directly approximating fractional derivative combining with numerical differentiation and integration.Recently, a new type of numerical schemes was designed to approximate the Caputo fractional derivative for solving time fractional partial differential equations, such as L1 method [28], L1-2 method [19], L2-1 σ method [4].These methods are all based on piecewise linear or quadratic interpolating polynomials approximation.It is natural to generalise the approach by improving the degree of the piecewise polynomial to approximate function that possesses suitable smoothness, in which situation the higher order of accuracy can be obtained.In the next section, we will devote to deriving a series of numerical schemes for formula (1.1) based on constructing piecewise interpolation polynomials on interval [0,t] as the approximations to solution u(t), and consequently, the α order Caputo derivative of the polynomials as the approximation to C D α u(t).The local truncation errors of the numerical schemes are discussed correspondingly.The flexibility via choosing different interpolating points on subintervals to construct the piecewise polynomials will produce various schemes under the similar restriction of accuracy order.
In order to study the numerical stability of such methods applying to problem (1.1), we will examine the behaviour of the numerical method on the linear scalar equation with initial value u(0) = u 0 .It is already shown that the solution of (1.4) satisfies that u(t) → 0 as t → +∞ provided that |arg(λ )| > απ 2 (−π ≤ arg(λ ) ≤ π) for arbitrary bounded initial value [32,36], accordingly, it can be studied in seeking those λ for which the corresponding numerical solutions preserve the same property as true solution.In fact, several classical numerical stability theories have been constructed on solving problem (1.4) in the case of α = 1 [21,22].Furthermore, there are some efforts on generalising the numerical stability theory on linear multistep methods to integral equations, such as Volterra-type integral equation [31,32].It is known that, for example, in the case of all those λ satisfying |arg(λ , it is referred to as A(θ )-stable.Inspired by the successive work, we make use of the technique pioneered in [32], specialise and refine the results to the fractional case in this paper.We confirm the stability regions and strong stability of the proposed numerical methods, and provide the rigorous analysis on the A( π 2 )-stability of some methods.Actually, it can be observed from the numerical experiments that the class of methods possesses the property of A(θ )-stability uniformly for 0 < α < 1, and for some specific α ∈ (0, 1), the A-stability can be obtained.
The paper is organized as follows.Section 2 introduces the continuous piecewise polynomial approximations of the Caputo derivative.We also derive some useful properties of the weight coefficients and discuss the local truncation errors.Section 3 and Section 4 respectively treat the stability and convergence aspects of the numerical schemes when applied to time fractional differential equations.In section 5 numerical experiments confirm our theoretical considerations with respect to order of convergence and stability restrictions.

Continuous piecewise polynomial approximation to the Caputo fractional derivative
We first introduce the fractional derivative in the Caputo sense: 14]).Let α > 0, and n = ⌈α⌉, the α order Caputo derivative of function u(t) on [0, T ] is defined by In particular, the Caputo derivative of order α ∈ (0, 1) is defined by whenever u (1) In the sequel, we will derive a class of piecewise polynomial approximations to the Caputo fractional derivative of order α ∈ (0, 1).The main idea is as follows.Let I = [0, T ] be an interval, the M + 1 nodes 2) is assumed to be continuous on the interval I, when we think about an piecewise polynomial approximation to u(t), it is reasonable to find the approximate solution at least in the continuous piecewise polynomial space, which is defined by Specifically, denoting the space of continuous piecewise polynomial of degree at most k by we then construct a class of approximate solutions of u(t) in the space C k p (I).Here the Lagrange interpolation technique by prescribing interpolation conditions on distinct k + 1 nodes is made use of such that the coefficients {a j,l } for 0 ≤ l ≤ k on each I j are uniquely determined.In addition, suppose that u(t) is not a constant function, we only need to focus on the cases of k ≥ 1 in view of the continuity restriction.The choice of the interpolating points is provided in the following way.
We define a class of polynomials p k j,q (t) which are of degree k ≥ 1 and have a compact support I j .The coefficients of the polynomials are uniquely determined by the following k + 1 interpolation conditions Here the index q records the number of shifts of the k + 1 interpolating nodes {t n } j−1 n= j−1−k , and the sign of q indicates the direction of the shift.Based on (2.4), the polynomials can be represented by a Lagrange form In particular, if the partition (2.3) is equidistant, i.e., t n = n∆t and ∆t = T M as M ∈ N + , the alternative Newton expression is given by For the convenience of notation, we then rewrite (2.6) by changing the variable t = t j−1 + s∆t to obtain where the r-th order backward difference operator ∇ r is commonly defined by and s−q+r−1 r is the binomial coefficient.In the following, we construct a class of approximate solution P k i (t) ∈ C k p (I) to u(t) on the uniform grid for 1 ≤ i ≤ k ≤ 6.The general representations are proposed by Remark 1.The construction of polynomials P k i (t) is mainly based on the continuity requirement on interval I, i.e., the interpolation conditions should be satisfied.It yields that on each subinterval I j , according to (2.4), both conditions j + q − 1 ≥ j and j + q − k − 1 ≤ j − 1 should be satisfied, which indicates 1 ≤ q ≤ k.Therefore, in the case of k = 1, there is a unique continuous piecewise linear polynomial, denoted by P 1 1 (t), in the space C 1 p (I).According to (2.8), it is expressed by In the other case of k = 2, there are three options on each I j , that is, p 1 j,1 (t), p 2 j,1 (t) and p 2 j,2 (t) to constitute the interpolating polynomial that belongs to space C 2 p (I).It is known that the construction of P 2 (t) is therefore not unique.In order to preserve the convolution property as much as possible, here we propose three available continuous piecewise polynomials in the forms of As a consequence, the operator is proposed on t ∈ I as the approximations to C D α u(t).In the case of t = t n , formula (2.11) can also be rewritten as where u n := u(t n ).
In the following part, we will present the explicit representations of the weight coefficients {w where q, r ∈ N + and n ∈ Z.In addition, note that Then the weight coefficients can be expressed in terms of integrals , n ≥ 3, and by more complicated formulae of forms (2.16) In addition, it is observed that when α → 1, the difference operator D α k,i u n in (2.12) recovers to the k-step BDF method.
Remark 2. The construction process of operator (2.11) can be extended to the case of α > 1 as well.In a general case of ⌈α⌉ − 1 < α < ⌈α⌉, the interpolating polynomials P k i (t) ∈ C k p (I) could be constructed as the approximations to u(t) under the condition that k ≥ ⌈α⌉, and the α order Caputo derivative of P k i (t) are proposed in an analogous way by as the approximation to C D α u(t).Here the condition of k ≥ ⌈α⌉ is required such that the ⌈α⌉ order derivative of P k i (t) is nonzero a.e..We take the case of ⌈α⌉ = 2 as an example.Assume that β = α − 1 ∈ (0, 1), the polynomials P 2 i (t) denoted by (2.10) are in the space C 2 p (I), and it follows where the last equality holds based on the relation Moreover, it can be rewritten as a form analogous to (2.12), and the corresponding weights coefficients {w n,i } and {ω n } are therefore derived by (2.17) Here the integrals I r n,q = I r n,q (β ) are defined by (2.13) where the index α is replaced by β .

Complete monotonicity and error analysis
First, we explore the completely monotonic property of the sequence {I r n,q } ∞ n=0 .Lemma 2.1.Assume that I r n,q is defined by (2.13), then for n ≥ k with k ∈ N, it holds that in the case of r ≤ q, and (−1) k+q+1 ∇ k I r n,q ≥ 0 (2.19) in the case of r > q.
Moreover, we discuss the complete monotonicity of a general class of sequences.
Lemma 2.2.The sequence {s n } ∞ n=0 is defined by Proof.It is easy to check that s n ≥ 0 for all n ≥ 0, since for n ≥ 0, 0 ≤ s ≤ 1, it holds that (n+1−s) −α > 0 and ϕ(s) ≥ 0 by assumption.The definition of s n implies that 0 and ϕ(s) ≥ 0 for n ≥ 1 and 0 ≤ s, ξ ≤ 1, thus ∇s n ≤ 0 holds.Therefore an induction process yields that Next, we construct the numerical scheme as the approximation to problem (1.1) with prescribed starting values, and define the local truncation error of the n-th step by where u(t) is the exact solution of problem (1.1). ) holds uniformly for n ≥ k.
Proof.According to (2.7), it holds that where t = t j−1 + s∆t with 0 ≤ s ≤ 1 and Inspired by [19], making use of the integration by part technique, we arrive at for n ≥ k, which is based on the conditions of (2.9) and (2.29).From the general representation of , in the other case of k = i, the polynomials of degree k are chosen on each subinterval I j instead.Therefore, we next consider the two cases seperately.
Substituting (2.8) and (2.29) into the last equivalent formula of (2.30), if k = i, one obtains Since for any q ≤ k and q, k ∈ N + , the factor is bounded for 0 ≤ s ≤ 1, thus we can obtain that where C (k) is bounded relevant to u (k+1) and k.Moreover, for i < k, it holds The error accuracy and convergence rate of |u(t where C (k,i) is a constant depending on u (k) , u (k+1) and k, i.
Remark 3. It is shown from formula (2.27) that the order accuracy isn't uniform for all n ≥ k.In the case of t n being near the origin, the accuracy order of the local truncation error reduced to the (k − α) order, in view that the (k − 1) degree polynomials as shown in (1) are chosen on the subinterval ∪ k−i j=1 I j .However, replacing polynomials of degree k on the corresponding subinterval can avoid this drawback, which is shown in (2.

28).
Remark 4.There is need to point out that the local truncation error estimations (2.27) and (2.28) holds only in the case of the solution u(t) possessing proper smoothness on the closed interval [0, T ].In order to check the convergence rate of the global error when f (t, u(t)) is smooth with respect to t and u, we apply the methods (2.12) in the cases of )

Stability analysis
To consider the numerical stability of schemes (2.25) with initial value u(0) = u 0 , the analysis on the linear difference equation is given as follows.Using formulae (3.1), we construct the equivalent relationship with respect to the generating power series where z := λ (∆t) α , the formal power series are denoted by Inspired by stability anlyses in [31,32], we therefore present the following preliminary conclusions.

Lemma 3.1 ([32]
). Assume that the coefficient sequence of a(ξ 3,45]).For the moment problem to be soluble within the class of non-decreasing functions iff the inequalities Lemma 3.2.The coefficient sequences of series g (k,i) (ξ ) converge to zero.
Proof.According to the expression of ∇ k I r n,q in Lemma 2.1, it yields that for k, q, r ∈ N + that are independent of n and α > 0. Note that g the finite linear combination of ∇ k I r j,q for finite k, thus it deduces g Proof.As indicated in Lemma 2.1 and Lemma 3.2, the following relationship holds for p ≥ k ≥ 1.Therefore, according to the definition of sequence {ω which implies the result.
Proof.According to the expression of ω (k,i) (ξ ), the following series can be rewritten to On one hand, from Lemma 3.3, we have On the other hand, by the definition of Combining with we arrive at the conclusion.
Corollary 3.3.For any |ξ 0 | ≤ 1 and 1 ≤ i ≤ k ≤ 6, it holds that the sequence belongs to l 1 space, where series ϕ (k,i) (ξ ) is defined to satisfy the relation ω Proof.Based on the definition of ϕ (k,i) (ξ ), it yields that because of the absolute convergence of sequences ξ −ξ 0 and ω (k,i) (ξ ) shown in Lemma 3.3 and Lemma 3.4, we shall arrive at the result.
On the other hand, assume that for any z = ω (k,i) (ξ 0 ) with |ξ 0 | ≤ 1, according to (3.2) the solution satisfies that Note that method (2.12) is exact for constant function, which leads to and a corresponding formal power series satisfies that In the case of ω (k,i) (ξ 0 ) = 0, it yields that u(ξ ) = u 0 1−ξ , which means that u n = u 0 for any n ∈ N.And for the rest case, there is If assume that u n → 0 as n → ∞, since according to Lemma 3.4, the coefficient sequence of ξ −ξ 0 tends to zero, in addition, according to Lemma 3.1, it yields that the coefficient sequence of ξ −ξ 0 converges to zero, however, the divergence of the coefficient sequence of 1 ξ −ξ 0 for |ξ 0 | ≤ 1 leads to the contradiction.Thus, it holds that there exist some nonzero bounded initial values According to the definition of A(θ )-stability [21] in usual case, we define the A(θ )-stability in the following sense of 0 < α < 1.
is contained in the stability region.

Convergence analysis
In this section, we consider the global error estimation for the problem (1.1) when the numerical approximations (2.25) are employed.Assume that u(t n ) is the exact solution of (1.1) at t = t n , then it satisfies where the difference operator D α k,i is defined by (2.12) and the local truncation error τ is the solution of (2.25) for each pair of (k, i), we denote the global error by where e (k,i) 0 = 0. Thus subtracting (2.25) by (4.1) implies that where the notation δ f In addition, substituting (2.12) into (4.3)yields Multiplying ξ n−k on both sides of (4.4) and summing up for all n ≥ k, one obtains where n,m are denoted by (4.6).Then for all n ≥ 0, it holds that s (k,i) n,m is bounded, and there exist some bounded constants c (k,i) m > 0 , which are independent of n and α, such that for n ≥ 1 and m ≥ 1.
Proof.It is known from (2.13) that for any finite q, r ∈ N + , I r n,q is bounded for all n ∈ Z. Since the coefficients s n,m are denoted as the linear combinations of I r n,q , we can immediately obtain the boundedness of s n,0 can be expressed as a linear combination of I l and I r l,1 with l ≥ n and 1 ≤ r ≤ 3. Based on formulae (2.20) and (2.23), it yields for r ≥ 2 and n ≥ 1, respectively.Therefore, it implies that there is a uniform bound with respect to n and α, denoted by c n,m are the linear combinations of ∇I l , I r l,1 and ∇ p I r l,1 , for l ≥ n + 1, r ≥ 2 and 1 ≤ p ≤ 3.According to formulae (2.21) and (2.24), we know that , and hence there exist constants c > 0 such that the last inequality of (4.7) is satisfied.
According to the definition of the series ω (k,i) (ξ ), it is important to notice the decompositions of the form where series ϕ (k,i) (ξ ) is defined by (3.12) and denote that Formula (4.8) indicates a relationship between the proposed method and the fractional Euler method mentioned in [35].In the following part, we would like to discuss some relevant properties of the series ψ (k,i) (ξ ) as preliminaries.
Proof.According to the expression of ϕ (k,i) (ξ ) presented in Theorem 3.6, there holds In addition, together with (4.9), it follows therefore, it suffices to prove that the coefficients of series (1 − ξ ) 1−α I(ξ ) belong to l 1 space.
From the gamma function's definition of the form one obtains the asymptotically equal relation where the notation ∼ = means the ratio n β −1 /Γ(β ) (−1) n −β n → 1 as n → ∞.Furthermore, it is known from [18,32] that On the other hand, based on the definition of I n , it yields that Therefore, in combination with (4.14), it holds that hence one has which yields the desired result.
Therefore, according to the statements from Theorem 3.1, Lemma 4.3 and Lemma 4.4, we can immediately obtain the following result.
then there exist bounded positive constants, denoted by M Theorem 4.2.Let u(t) and {u n } N n=k be the solutions of equations (1.1) and (4.4), respectively.The function f (t, u(t)) in (1.1) is assumed to satisfy the Lipschitz continuous condition with respect to the second variable u, and chosen properly such that the solution of (1.1) is sufficiently smooth.Then ii) in the cases of for sufficiently small ∆t > 0, where N∆t = T is fixed and constant C (k,i) > 0 is independent of N and n.
Proof.Substituting formula (4.8) into (4.5) and using (4.17), one has which can also be written into a matrix-vector form with arbitrary N ∈ N. Therefore for any k ≤ n ≤ N, it holds that where the coefficients {g On one hand, based on the relations (4.7) and (4.16), there exist constant ck,i > 0, such that |s since it is known that ∑ where denote by , and in the cases of 1 ≤ k ≤ 3, one has from (4.24), (4.25) and (4.26) that δ and in the cases of 1 ≤ i < k ≤ 3, from formulae (4.24), (4.25) and (4.27), it yields Then assume that the non-negative sequence {p where the coefficient L(k,i) is chosen such that Therefore, according to the weakly singular discrete Gronwall inequality shown in [15], the monotonic increasing property of sequence { n } n≥1 is monotonic increasing with respect to n for each 1 ≤ i ≤ k ≤ 3, and correspondingly, it follows where E α (•) is denoted by the Mittag-Leffler function.In addition, according to (4.28) and (4.29), an induction process yields that |e Remark 6.The provided convergence order is uniform for all n ≥ k, especially suitable for the step t n near the origin.On the other hand, for those t n away from origin, the convergence result can be better.For example, it can be observed numerically if the computed starting values satisfy u m = u(t m )+O((∆t

Numerical experiments
In this section, we utilize formula (2.12) to approximate the following equations in Example 5.1 and Example 5.2, and we prescribe the starting values exactly.Since in practical computation, the starting values are normally obtained by numerical computation in advance.
Example 5.1.We consider the linear equation , where the Mittag-Leffler function [40] is defined by (d) α = 0.98, λ = 500i  n e inθ (0 ≤ θ ≤ 2π) in the cases of 1 ≤ i ≤ k ≤ 3 for different α ∈ (0, 1).It is known from Theorem 3.4 that the stability regions of methods (3.1) lie outside the corresponding curves.We introduce the points z n := λ (∆t n ) α for 1 ≤ n ≤ 5, where ∆t n = 1/2 n+6 correspond to different time stepsize.Table 3 3-4, we can see the influence of the stability of the numerical methods on the error accuracy.In Figure 1(a), the points z n for 1 ≤ n ≤ 5 all lie in the stability regions in the case of α = 0.5 and λ = −50, in which situation the reliable accuracy is obtained, and it is observed that |e M | = O(∆t k+1−α ) in Table 3-4.In Figure 1(b)-1(c), z n are choose on the half line with angle πα 2 and different λ , it is observed that when all {z n } 5 n=1 fall out of the instability region (cf. Figure 1(b)), correspondingly, as shown in Table 3-4, the the global error e M agrees with the expectation of the accuracy.On the other hand, due to points z 4 and z 5 outside the stability regions of k = 3 (cf.Figure 1(c)), the perturbation errors are magnified and accumulated significantly, which are shown in Table 3-4 as well.In Figure 1(d), z n are choose on the imaginary axis with imaginary number λ , according to Theorem 3.5, all the z n belongs to the stability region of the method methods (3.1) in the case of k = 1, 2, and the error accuracy and the convergence order are obtained (cf.Table 3).As a counter example, when z 3 doesn't belong to the stability region for α = 0.98 in Figure 1(d), the corresponding error e M shown in Table 4 can't ensure the desirable accuracy.In fact, it can be observed that for k = 3, methods (3.1) don't possess A( π 2 )-stability when α tends to 1, which appears to be predictable, since it is well known that BDF3 method for ODEs is not A( π 2 )-stable.

Conclusions
We have proposed a new higher-order approximation method for solving time-fractional initial value models of order 0 < α < 1.Furthermore, the local truncation error in terms of solution possessing sufficient smoothness is derived and demonstrated.Additionally, the stability and convergence analyses of the proposed method are provided in details, which also motivate the relevant research on applications to time-fractional partial differential equations.
and the coefficient sequence of series g (k,i) (ξ ) z converges to zero.In addition, assume that lim n→∞ n ∑ j=0 |l i | = L < +∞ and lim

Figure 1 :
Figure 1: The boundary of the stability region for different α and λ .
e M = u(t M ) − u M in Example 5.1, where t M = M∆t = 1 is fixed and M = 2 j for 7 ≤ j ≤ 11, u(t M ) and u M are the exact solution and the computed solution, respectively.By the comparison of Figures 1(a)-1(d) and Tables

Figure 2
Figure 2 and 3 plot the global error |e M | = |u(t M ) − u M | in Example 5.2 for different µ and α, where t M = 1 is fixed and ∆t = 1/M with M = 2 j , 2 ≤ j ≤ 11.It is observed that |e M | = O(∆t k+1−α ) in the cases of 1 ≤ i ≤ k ≤ 3, when using discretisation formula (2.25) in combination with Newton's method for the nonlinear equation behind the implicit method.

Table 1
and 2, the accuracy and the convergence order of the error |u(t M ) − u M | are shown for different timestep and order α.According to the numerical experiment, the high order convergence seems to reduce to the first order in the cases of1 ≤ i ≤ k ≤ 3.It is because that the solution of the problem (1.1)only possess continuity on interval I if function f (t, u(t)) is smooth on I. On the other hand, the integer order derivative of any smooth function on a compact domain I still preserves to be smooth on I, in contrast, the α order fractional derivative of the smooth function isn't smooth any more, which implies there exist some continuous functions f (t, u) such that the solution is smooth on I.

Table 2 :
The error accuracy and convergence rate of |u(t M ) − u M | in problem (2.31).

Table 3 :
The error accuracy and convergence rate of |u(t M ) − u M | in Example 5.1 for different α and λ .

Table 4 :
and Table 4 list the global error The error accuracy and convergence rate of |u(t M ) − u M | in Example 5.1 for different α and λ .