Continuous and discrete one dimensional autonomous fractional ODEs

In this paper, we study 1D autonomous fractional ODEs $D_c^{\gamma}u=f(u), 0<\gamma<1$, where $u: [0,\infty)\mapsto\mathbb{R}$ is the unknown function and $D_c^{\gamma}$ is the generalized Caputo derivative introduced by Li and Liu ( arXiv:1612.05103). Based on the existence and uniqueness theorem and regularity results in previous work, we show the monotonicity of solutions to the autonomous fractional ODEs and several versions of comparison principles. We also perform a detailed discussion of the asymptotic behavior for $f(u)=Au^p$. In particular, based on an Osgood type blow-up criteria, we find relatively sharp bounds of the blow-up time in the case $A>0, p>1$. These bounds indicate that as the memory effect becomes stronger ($\gamma\to 0$), if the initial value is big, the blow-up time tends to zero while if the initial value is small, the blow-up time tends to infinity. In the case $A<0, p>1$, we show that the solution decays to zero more slowly compared with the usual derivative. Lastly, we show several comparison principles and Gr\"onwall inequalities for discretized equations, and perform some numerical simulations to confirm our analysis.


Introduction
The fractional calculus in continuous time has been used widely in physics and engineering for memory effect, viscoelasticity, porous media etc [1][2][3][4][5][6][7]. There are two types of fractional derivatives that are commonly used: the Riemann-Liouville derivatives and the Caputo derivatives (See [3]).
The Riemann-Liouville derivatives are named after Bernhard Riemann and Joseph Liouville. Liouville was the first to study fractional derivative rigorously (see, for example, [8,9] for a better survey). On the other hand, the Caputo's definition of fractional derivatives was first introduced in [10] to study the memory effect of energy dissipation for some anelastic materials, and soon became a useful modeling tool in engineering and physical sciences to construct physical models for nonlocal interactions in time (see [11]).
Compared with Riemann-Liouville derivatives, Caputo derivatives remove the singularities at the origin and have many properties that are similar to the ordinary derivative so that they are suitable for initial value problems [12]. However, the classical γ-th Caputo derivative of a function requires an integer order derivative no less than γ, which seems to be artificial. In [6], Allen, Caffarelli and Vasseur have introduced an alternative form of Caputo derivatives to avoid using higher regularity of the function. In [12], another extension of Caputo derivatives was proposed, by which the higher derivative of the function is not needed either and can recover the definition in [6]. Moreover, this new definition allows us to transform fractional ODEs into Volterra type integral equations, by deconvolution through an underlying group property without the higher regularity assumption. This provides a convenient framework for us to study the fractional ODEs with Caputo derivatives.
There is a huge amount of literature discussing fractional differential equations. However, few of them discuss the behavior of the solutions to fractional ODEs systematically. For reference, some results can be found in [2,5] using the traditional Caputo derivatives.
In this paper, we use the new definition of Caputo derivative in [12] (also see Definition 2.2 below) to make a detailed investigation of the nonlinear fractional ODE (1.1) D γ c u = f (u), u(0) = u 0 , for γ ∈ (0, 1). Here f is locally Lipschitz whose domain contains u 0 , and D γ c represents the Caputo derivative of order γ. In the rest of this paper we will assume u 0 0 without loss of generality (if u 0 < 0, we can do change of variables v = −u and study D γ . Studying the behavior of the solution to this fractional ODE is important for the analysis of fractional partial differential equations (fractional PDEs), as we usually need a priori estimates of certain energies of the solution to a fractional PDE, which have form D γ c E AE p . By the comparison principles in [12] or in Section 4, the energy norm may be controlled by the solution of the fractional ODE (1.1). Hence, we will focus on the particular cases f (u) = Au p in detail.
According to [12], the fractional ODE (1.1) is equivalent to a Volterra type integral equation without assuming high regularity of the solution, which is the important starting point for our study. It is well-known that the solutions of 1D autonomous ODEs with usual first order derivative are monotone, since the solution curves never cross zeros of f and f (u) has a definite sign. One of our main results is that if f ∈ C 1 and f is locally Lipschitz, the first order derivative of the solution to the fractional ODE (1.1) does not change sign and therefore the solution is monotone (see Theorem 3.4). This is based on Lemma 3.5, which is a slightly different version of [13,Theorem 1]. Lemma 3.5 ensures the positivity of the solutions of the integral equation that y = u or y = −u satisfies: where v is continuous. The idea is to use the resolvent for the kernel λt γ−1 to transform this integral equation into another integral equation (see (3.5)) so that all the functions involved are non-negative. The solution to the new integral equation (3.5) is nonnegative, implying that the first derivative of the solution to (1.1) does not change sign.
Another contribution of this paper is to discuss the special cases f (u) = Au p in detail and to reveal several interesting roles of memory. In particular, for the cases A > 0, p > 1, we find relatively sharp estimates of the blow-up time T b . The lower bound of T b is important for the inequality D γ c E AE p since it ensures that E is defined and controlled by the solution of (1.1) up to this lower bound. Through these bounds, we find that there exist u 02 > u 01 > 0 so that if u 0 < u 01 , the blow-up time T b → +∞ as γ → 0 (the memory becomes stronger) and if u 0 > u 02 , T b → 0, as γ → 0 (See Theorem 5.3). For the cases A < 0, p > 1, we show that under the memory, the solution decays to zero much more slowly compared with the usual ODE (see Theorem 5.10). By discretizing the differential equation (1.1) or the equivalent integral equation, we obtain two classes of numerical schemes or discrete equations. Using some discrete comparison principles, we show that if f is nonnegative, nondecreasing, then the numerical solutions to the explicit schemes for the integral equation are absolutely stable: u n u(nk) where k is the time step (Proposition 6.3). In the case f is nonnegative, nondecreasing and the solution to (1.1) is convex, we prove that the numerical solutions to the explicit schemes for the differential equation are also absolutely stable: u n ex u(nk) and the numerical solutions to the implicit schemes for the differential equation are bounded below as u n im u(nk) (Theorem 6.7). Hence, the explicit schemes may be used to prove the stability and convergence of some approximation schemes for fractional PDEs and thus the convergence and existence of solutions. The implicit schemes may be used to prove positivity of solutions and to estimate the blow-up time.
The rest of the paper is organized as follows: In Section 2, we introduce the basic definitions, notations and results that are mainly established in [12]. In Section 3, we study the basic properties of the solutions. In particular, (1) given f (u) is smooth, the solutions are smooth in (0, ∞) but only γ-Hölder continuous at t = 0; (2) the solutions are monotone on the interval of existence; (3) an Osgood type finite time blow-up criteria holds provided that f (u) is positive nondecreasing. In Section 4, we prove several comparison principles. In Section 5, we study the special cases f (u) = Au p . More precisely, for A > 0, p > 1, we provide relatively sharp bounds for the blow-up time, while for A < 0, p > 1, we show the slow decaying as t → ∞. These discussions reveal the roles of memory introduced by fractional derivatives. Lastly, in Section 6, we discuss the discrete equations. To be more specific, we show several discrete comparison principles and use them to study some explicit and implicit schemes. Some numerical simulations are then performed using these schemes to verify our analysis for the continuous cases.

Preliminaries
In this section we collect some notations and definitions we will use in this paper.
2.1. Fractional derivatives. First, let us make a brief introduction of the definition of fractional derivatives. Before we state the definition, we need the following clarification of notation: we call u 0 the right limit of u at t = 0, denoted as u(0+) := u 0 .
As in [12], we use the following distributions {g β } as the convolution kernels for β > −1: . Here θ(t) is the standard Heaviside step function, Γ(γ) is the gamma function, and D means the distributional derivative. g β can also be defined for β −1 (see [12]) so that these distributions form a convolution group {g β : β ∈ R}, and consequently we have where the convolution between distributions with special non-compact supports is defined through the partition of unit of R. Now we are able to give the definition of the fractional derivatives.
Definition 2.2. Let 0 < γ < 1. Consider u ∈ L 1 loc (0, T ) that has a right limit u(0+) at t = 0 in the sense of Definition 2.1. The γ-th order Caputo derivative of u is a distribution in D (−∞, T ) with support in [0, T ), given by 3. In the case T = ∞, the convolution g −γ * u is defined through partition of unit of R. In the case of T < ∞, g −γ * u should be understood as the restriction of the convolution onto D (−∞, T ). One can refer to [12] for the technical details.
Remark 2.4. As discussed in [12], if there is a version of u (i.e. modifying u on a Lebesgue measure zero set) that is absolutely continuous on (0, T ), which is denoted as u again, then the Caputo derivative is reduced to which is the traditional definition of Caputo derivative. Whenever u is γ + δ-Hölder continuous (∀δ > 0), we have is the definition for the Caputo derivative used in [6]. Intuitively, (2.4) is obtained by integration by parts from (2.3). Definition 2.2 is more useful than the traditional definition (Equation (2.3)) (see for instance [1][2][3][4][5]7]) theoretically, since it asks for little regularity and reveals the underlying group structure. With the assumption that u is locally integrable and has a right limit at t = 0, Definition 2.2 and the group property (2.2) allow one to convert (1.1) into the integral form Remark 2.5. Obtaining Equation (2.5) from the traditional Caputo derivative (2.3) needs us to assume in advance that the unknown function u has too much regularity (for example, for the definition to make sense, we have to assume that u (s) exists). Using the new definition of Caputo derivative in [12], the integral form (2.5) is equivalent to Equation (1.1) with the assumption that u is locally integrable and has a right limit at t = 0 in the sense of Definition 2.1 only. One can check [2,3,5,12] for more details. Equation (2.5) is called Volterra integral equation, which has been studied extensively. Analysis of Equation (1.1) (or equivalently (2.5)) can help us understand the time-delay properties of Caputo derivatives.

2.2.
Existence and uniqueness of solutions to (1.1). In this paper, we will use the following definition of solutions: Definition 2.6. u(·) ∈ L 1 (0, T ) that has a right limit at t = 0 in the sense of Definition 2.1 is called a weak solution of (1.1) if the equation is satisfied in the distribution sense and u(0+) = u 0 . A weak solution u is called a strong solution if D γ c u ∈ L 1 (0, T ) and (1.1) is satisfied almost everywhere with respect to Lebesgue measure.
By the equivalence of (1.1) and (2.5) established in [12], all weak solutions of (1.1) satisfy the integral equation (2.5) almost everywhere with respect to Lebesgue measure. By modifying the result in [12,Theorem 6], we have the following proposition: is locally Lipschitz continuous on an interval (α, β) ⊂ R, then ∀u 0 ∈ (α, β), there is a unique continuous strong solution with u(0) = u 0 . Either this solution exists globally on [0, ∞) or there exists T b > 0 such that either The claim is essentially the same as [12,Theorem 6], so we omit the proof. 3.1. Regularity and monotonicity of solutions. In this subsection, we present and prove the regularity and monotonicity results of solutions to (1.1). Lemma 3.1 is the result proved in [14] for integral equations. This lemma gives the regularity of the solutions to (1.1) and lays the foundation for our later discussion. Theorem 3.4 is the main result in this subsection, which states that the solutions of the autonomous equations are generally monotone. The proof of this theorem relies on Lemma 3.5, which ensures the positivity of the solutions to a certain class of integral equations. Lemma 3.5 is a slightly different version of [13, Theorem 1]: the author of [13] assumed y and h to be continuous at t = 0 but we cannot assume this for our purpose. However, the idea of the proof is the same. We present the regularity lemma ( [14, Theorem 1]): Moreover, y = u satisfies the integral equation For the idea of proof, one may fix T ∈ (0, T b ) and show that (3.1) has a unique continuous solution on (0, T ). Then using the equation for (u(t+h)−u(t))/h, one can verify that this finite difference converges to the solution of (3.1). One can refer to [14] for a detailed discussion. For the last claim, as long as we have u = O(t γ−1 ), we can show that the integral is then dominated by the first term as t → 0 + . Remark 3.2. Using the group property g nγ * g γ = g (n+1)γ (Equation (2.2)), we may find that the solution to (1.1) is a power series of t γ if f is real analytic. This observation tells us that t γ power is intrinsic to the Caputo derivative D γ c and the solution is only γ-Hölder continuous at t = 0, but smooth on (0, T b ).
In the following theorem, we will show the sign of f (u(t)) does not change: Proof. Without loss of generality, we assume f (u 0 ) > 0. Define First of all, we have that t * > 0 since f (u 0 ) > 0. To prove the theorem, we only need to show t * = ∞.
0 for t ∈ (0, t * ), and f (u(t)) is strictly positive when t is close to 0. In addition, f (u(t)) 0 for t ∈ [t * , t * + δ]. Therefore, the right hand side is strictly negative. Hence u(t) < u(t * ) and the claim is proved.
By the continuity of u, we conclude that there exist t 1 , t 2 , 0 For ordinary derivative, as long as we have shown that f (u(t)) has a definite sign, we have that the solution is monotone. For fractional derivatives, this is not obvious, however we can also show that this is true provided f is close to C 2 . More precisely, we have: for some interval (α, β) and f is locally Lipschitz on (α, β). Then, the solution u to (1.1) with u(0) = u 0 ∈ (α, β) is monotone on the interval of existence (0, T b ), where T b is given by Proposition 2.7. If f (u 0 ) = 0, the monotonicity is strict.
Before proving this theorem, let us prove a useful lemma that ensures the positivity of the solution to an integral equation, which is a slightly different version of [13, Theorem 1] (for more discussions on positivity of solutions to Volterra equations, see [13,15]) : Here r λ is the resolvent for kernel λt γ−1 satisfying Suppose v ∈ C 0 [0, T ], then the integral equation Proof. It can be computed explicitly that is the Mittag-Leffler function [16,17]. E γ (−λΓ(γ)t γ ) is completely monotone that goes from 1 to 0 on (0, ∞) [16]. For the concept of completely monotone, see [18]. As a result, r λ ∈ L 1 (0, T ) ∩ C 0 (0, T ] and r λ > 0. Then, by the fact that the convolution of two locally integrable functions is again locally integrable, all the convolutions are well-defined. (Actually, by an abstract argument, it has also been shown in [15, Lemma 2.1] that r λ is completely monotone and thus non-negative.) The existence and uniqueness of (3.3) are shown in [14,Lemma 1]. We now prove y > 0, a.e.. As Taking the difference between (3.3) and (3.4), . As a result, y also solves the integral equation λ > 0 and then y 0 a.e. on [0, T ) follows from [14, Lemma 1]. Now from the assumption of this lemma, we have h(t) − t 0 r λ (t − s)h(s)ds > 0, in addition we also have r λ > 0 and 1 − v λ > 0. As a result, y > 0 a.e. from (3.5).
The last claim is proved.
Remark 3.6. In the proof of Theorem 1 of [13], the author assumed h to be continuous and the solution to be continuous at t = 0. In Lemma 3.5, we do not assume y to be continuous, which is crucial in the case that h(t) = αt γ−1 .
Now, we are able to prove Theorem 3.4: Proof of Theorem 3.4. Clearly, if f (u 0 ) = 0, then u = u 0 is the solution by the uniqueness. This is trivially monotone.
. The derivative y = u satisfies the equation Since f (u(t)) is continuous on [0, T ] and f (u 0 ) > 0, applying Lemma 3.5, we find y is positive on (0, T ). Since T is arbitrary, y > 0 on (0, T b ). As a result, u is increasing.
If f (u 0 ) < 0, we simply consider the equation for y = −u . The argument is similar.
For usual ODEs, the solution curves do not intersect. For fractional ODEs, we can conclude directly from the integral equation (2.5) that Proposition 3.7. If f (u) is locally Lipschitz continuous, non-decreasing, then the solution curves of (1.1) do not intersect with each other.
Remark 3.8. In the case that f (u) is non-decreasing only on some interval, then as long as one can show that the solutions stay in this interval, then the curves with initial value in this interval does not intersect. For general f , it is unclear whether or not the solution curves intersect. The memory is playing a tricky role.

3.2.
Blow-up criterion. Now we present some results regarding the blow-up behavior. We first have the following observation Proof. First of all, let us show that the solution u is non-decreasing on (0, T b ). Note that the f in this lemma is less regular than the function in Theorem 3.4, hence we cannot use Theorem 3.4 directly. To show the monotonicity of u, let us consider the following sequence of functions {u n } ∞ n=0 : u 0 = u 0 , D γ c u n = f (u n−1 ), u n (0) = u 0 , n 1.
From the integral form of the fractional derivative (2.5), it is clear that u n is continuous on [0, ∞). Since f (u 0 ) > 0, we have By induction, u n (t) u n−1 (t) for all n 1.
Next, we claim that u(t) > u 0 for t ∈ (0, T b ). For this purpose, we define t * = sup{t ∈ (0, T b ) : f (u(t)) > 0, ∀t ∈ (0,t)}. We show that t * = T b . First of all, according to the continuity of u(t) and f (u), and the fact f (u 0 ) > 0, we have t * > 0. If t * < T b , then f (u(t * )) = 0 by the continuity of f and u. In addition, by the definition of t * : Since f is non-decreasing, we have f (u(t * )) f (u 0 ) > 0, which is a contradiction.
Using u(t) u 0 , we find for t ∈ [0, T b ). Again by induction, u(t) u 2 (t) and u(t) u 3 (t), etc. Moreover, since f is non-decreasing and f (u n ) is positive, we find that for any 0 t 1 < t 2 < ∞: The second equality here is achieved by a change of variable τ = s − (t 2 − t 1 ). Hence, u 1 is non-decreasing on [0, ∞). Similarly, u 2 is non-decreasing on [0, ∞). By induction, u n is non-decreasing. As a consequence, the sequence {u n (t)} converges to a non-decreasing functionū(t) for any t ∈ [0, T b ). By monotone convergence theorem and taking the limit both sides of u n (t) = u 0 + 1 Γ(γ) The next result, which is an Osgood type criterion is essentially from [19] for the Volterra type integral equations. Here, we reinterpret it for our fractional ODE (1.1), and using similar ideas we present an improved proof, which enables us to improve the bounds of blow-up time in Section 5: Proof. Consider the equivalent Volterra type equation (2.5). By Lemma 3.9, u is increasing and . There exists t n < T b so that u(t n ) = r nγ for n = 1, 2, . . .. By (2.5), we have As a result, there exist constants C(γ) > 0 and C 1 (γ, r) > 0 such that On the other hand, As a result, there exist two constantsC 1 (γ) > 0 andC 2 (γ, r) > 0 such that

Comparison principles
The following comparison principle ( [12, Theorem 7]) is useful when we study the behavior of (1.1) and derive certain Grönwall type inequalities: is a nonpositive distribution (see [12,Def. 6]). Let v 2 be the unique solution to the equation Using the idea of the proof for [12, Theorem 7], we are able to show some other versions of comparison principles. For example, an integral version is as follows then we have v u on (0, min(T, T b )).
Note that the integral version is not a pure repetition of Proposition 4.1 since we do not ) for all t ∈ (0, T ). Another version of comparison principle is a corollary of Proposition 4.1: Corollary 4.3. Suppose both f 1 (u) and f 2 (u) are locally Lipschitz on (α, β), satisfying f 1 (u) f 2 (u) for any u ∈ (α, β). Assume that one of them is non-decreasing. Let u 1 and u 2 be the solutions of D γ c u = f 1 (u) and D γ c u = f 2 (u) on the intervals (0, T 1 b ) and (0, T 2 b ) with initial values u 1 (0) and u 2 (0) respectively. If in addition α < u 2 (0) u 1 (0) < β, then First, assume that f 2 is non-decreasing. Then, we have . Applying Proposition 4.1 for D γ c u 2 f 1 (u 2 ) yields the claim.

blowup and long time behavior for a class of fractional ODEs
In this section, we will focus on the cases f (u) = Au p and (α, β) = (0, ∞) for simplicity.
This type of equations are general enough. For example, if p > 0 and there exist C 1 > 0, C 2 > 0 such that C 1 u p f (u) C 2 u p , then the solution is under control according to Corollary 4.3. We will discuss in different cases to show that (5.1) shares the regularity properties of normal time derivative ODE. Moreover, one can also prove some "time-delay properties" of (5.1).

Finite time blowup.
As an application of Proposition 3.10, we have the following theorem:

The bounds of blow-up time.
In this subsection, our main goal is to find suitable bounds of the blow-up time and to understand the effects of the memory introduced by the Caputo derivatives. Clearly, one possible lower bound is the radius of convergence of the power series u = ∞ n=0 a n t nγ , however the asymptotic behavior of a n is hard to find. In [20], the author provided some bounds for the blow-up time of the integral equation (2.5). In this paper, we have the following improved result: Proposition 5.2. Suppose γ ∈ (0, 1), p > 1, A > 0, and u 0 > 0. Let T b be the blow-up time of the solution to (5.1). Then, we have the following inequality Proof. Let r > 1. We now choose t n such that It is then clear that 0 = t 0 < t 1 < t 2 . . .. For convenience, we denote The following relation yields that As a result, To prove the upper bound, we fix m 1, and then find and u(tm) For n m + 1, we find Combining (5.5) and (5.6), we finally have the upper bound, In the proof of the upper bound, the estimate we did for t m essentially follows the method in [20]. By optimizing the constants we get in the Proposition 5.2, we have Consequently, with A > 0, p > 1 fixed, there exist u 02 > u 01 > 0 such that whenever u 0 < u 01 , Proof. For the lower bound, picking r = 2 1/γ > 1, we find sup r>1 (r γ − 1) 1/γ r p − r Similarly, picking r = (p/(p − 1)) 1/γ yields For the upper bound, we fix m > 1 p−1 , and let r → ∞: If instead we choose m = 1 and r = 2 1/(p−1) > 1, we have Hence, For the second inequality, note that γ − 2 γ + 1 is concave on (0, 1) and equals zero at γ = 0, 1, so γ > 2 γ − 1 for γ ∈ (0, 1). We find and the upper bound follows.
As long as we have these two bounds, it is clear that we can pick . Remark 5.4. From Theorem 5.3, one can clearly see how the memory plays the role. The memory is getting stronger as γ goes closer to 0. When u 0 is very small, the strong memory defers the blowup. If u 0 is large, the strong memory accelerates the blowup. For the critical value of u 0 , we believe it is determined by the limiting case γ = 0: , this algebraic equation has no solution and it means the blow-up time , there is a constant solution for t > 0 which means the blow-up time is infinity.
Remark 5.5. The estimates (p−1) p−1 p p and ( 1 p−1 ) 1/γ for the blow-up time can also be obtained by the results in [21] for the Volterra integral equations, but we have better constants here. One may observe the following two facts:  ∈ (1, 2), G(p) = 2 p in (5.8) while for p ∈ (2, ∞), G(p) = p p (p−1) p−1 . The latter gives asymptotic behavior for large p.
Remark 5.6. One may wonder the asymptotic behavior, or so-called growth rate of the solution near the blow-up time. There are a lot of references about this topic. One can check, for instance [20,21]. To find the correct power of the blow-up profile, one can plug 1 (T −t) α into (5.1) and use the heuristic calculation D γ c ( 1 (T −t) α ) ≈ 1 (T −t) α+γ , which means α + γ = pα, or α = γ p−1 . In fact, from (3.2) in [21], the solution to (5.1) satisfies One can find the proof in the appendix (Section 7). In addition, as in [20,21], one can expect explicit asymptotic behavior for more general f (u).

Other cases.
In this subsection, we discuss other choices of the parameters A and p in (5.1).
First of all, we investigate the cases A > 0 and p < 0.
By Proposition 2.7, D γ c v =f (v) has a unique solution v with an interval of existence [0, T b ). Clearly,f (v(t)) > 0 on [0, T b ), and consequently v(t) > u 0 for t ∈ (0, T b ). This implies that v is actually the solution to D γ c u = Au p on [0, T b ). We therefore identify v with u. The monotonicity of u follows from Theorem 3.4.
We compare u with the solution of D γ c w = Au p 0 , w(0) = u 0 using Corollary 4.3 and find that This implies that T b = ∞.
Since f (u) = Au p (A < 0, p ∈ R) is smooth on (0, ∞) and f (u 0 ) = 0, by Theorem 3.4, u is strictly monotone. Using the integral form (2.5), it is clear that u(t) < u 0 for t > 0. Hence u is decreasing. From Proposition 2.7, either T b = ∞ or T b < ∞ and lim t→T − b u(t) = 0. To finish the proof, we only need to show that if T b = ∞, lim t→T b u(t) = 0. Suppose for otherwise lim t→T b u(t) = 0. Then, u(t) is bounded below by δ > 0. Then, as t → T b = ∞, which is a contradiction.
Remark 5.9. In the case p < 1, it is possible that T b < ∞ and Au p is defined on R. The solution may be extended beyond T b . However, Au p may not be Lipschitz continuous at u = 0 and it makes the analysis complicated (of course p = 0 case is trivial and we have u(t) = u 0 + Ag γ+1 ).
In the case p < 1, Au p may not be Lipschitz at u = 0. Hence, for simplicity, we only consider A < 0, p 1 for further discussion. Actually, we are able to show: Theorem 5.10. Fix γ ∈ (0, 1), p 1, u 0 > 0 and A < 0. Let u(t) be the unique solution to (5.1) with initial value u 0 . Then, u(t) > 0, ∀t > 0. u(·) is decreasing and lim t→∞ u(t) = 0. Moreover, there exists C(u 0 , A, p) > 0 such that when t is large enough, Proof. First of all, by Proposition 5.8, u is decreasing. Pick r ∈ (0, 1). By the fact lim t→T − b u(t) = 0, we are able to pick disjoint intervals J n = (t n−1 , t n ) such that u stays between u(t n−1 ) and u(t n ) inside J n and u(t m ) = u 0 r mγ . Therefore,
Since u ∈ C 1 (0, T b ) ∩ C 0 [0, T b ), integrating by parts from (2.3) gives us the alternative expression for Caputo derivative (2.4). Now, u(t) u(s) for all s t and as a result, When t is large enough, Remark 5.11. The proof of Theorem 5.10 is quite indirect. The equation may be rewritten as and |A|u p is an m-accretive operator (see [22]) of u when p > 1, u > 0. This form is related to the equations studied in [22] and may yield some direct proof using functional analysis. In the case that the kernel is not L 1 , [22] requires that m-accretive operator to be coercive which does not apply here.
Remark 5.12. It is well known that γ = 1 yields u(t) ∼ Ct −1/(p−1) , which decays to zero faster than t −γ/p . The memory really gives a slow decaying rate. As γ → 1, Γ(1 − γ) → ∞ and the dominant term in (5.11) vanishes. This means t −1/(p−1) must appear in the next order and the slow decaying dominate term (5.11) is an effect of memory.
Remark 5.13. Regarding the asymptotic behavior of Caputo derivative, we may consider the derivative of (1 + t) p . If p > 0, p is smooth and one can use (2.3) to compute. In the decaying cases p < 0, This means no matter how fast the function decays, the Caputo derivative is always like −Ct −γ asymptotically, which can also be confirmed through the proof of Theorem 5.10.
Actually, −t −γ should be the intrinsic rate for the Caputo derivative of decaying functions. If, for example, D γ c u(t) −C(1 + t) −γ+δ for some C > 0 and δ > 0, then u(t) → −∞. Conversely, if D γ c u(t) ∼ −(1 + t) −γ−δ , then u, though is less than u 0 , will eventually go back to u 0 . Notice that though the Caputo derivative is negative, the function does not always decay. This is because the decaying property at the earlier stage lingers to later stage due to memory.

Discrete equations and numerical simulations
In this section, we study discrete equations obtained from discretizing the differential equation (1.1) or the integral equation (2.5). We will consider some typical numerical schemes which are useful in different situations (e.g. stability analysis for numerical schemes or the proof of existence of weak solutions to fractional PDEs).
In this section, k > 0 is the time step, and t n = nk. u n is the computed numerical value at t n and u(t n ) is the value of the solution to (1.1) evaluated at t n .
6.1. Schemes for the integral equation. Consider discretizing (2.5) with explicit schemes. We have To study these two schemes, we first prove the following discrete Grönwall inequalities: Lemma 6.1. Let f (u) be nonnegative, non-decreasing, locally Lipschitz on [0, ∞) and let u 0 > 0. Suppose w n (0 n N ) is a nonnegative sequence (w n 0) such that Proof. We prove by induction. n = 0 is clearly true. Now, let 1 n N and assume that w m u(mk) for all m n − 1. Then, by the non-decreasing property of f and the induction assumption, we have Since f (·) is non-negative, by Lemma 3.9, u(·) is increasing. As a result, The proof for the other inequality is similar due to the fact This lemma recovers the discrete Grönwall inequality in [23]: 23]). Let {a n } be a non-negative sequence. If {a n } satisfies where B > 0 and λ > 0 are independent of n, k, γ, then, a n u(nk) = BE γ (λ(nk) γ ), 0 n N.
Here, u(t) is the solution to D γ c u = λu with initial value B and E γ is the Mittag-Leffler function. We conclude the following stability result about the schemes, which is useful when studying numerical schemes of fractional PDEs. Proposition 6.3. Let f (u) be nonnegative, non-decreasing, locally Lipschitz on [0, ∞) and u 0 > 0. Suppose u n solves the numerical scheme (6.1) or (6.2), and u is the unique solution to (1.1) with initial value u 0 . Then, we have u n−1 u n u(nk), 1 n < T b /k.
Proof. u n u(nk) follows directly from Lemma 6.1. Now, let us show that {u n } is non-decreasing under the scheme (6.1) by induction. For n = 1, it is clear that u 1 u 0 by direct computation. Now, let n 2 and assume that u m u m−1 for all 1 m n − 1.
The proof for scheme (6.2) is similar. 6.2. Schemes for the differential equation. Now, let us discretize Equation (1.1) directly. We assume the solution is C 1 (0, T b ), and use the following first order scheme from [24,25] based on the explicit formula (2.3): We can determine that for all n 0, does not depend on n if m n. Hence, for simplicity, we write Using this basic discretization, we can formulate the explicit and implicit schemes respectively as First of all, we discuss the explicit scheme: The following result follows from the facts b m < 0 for m 1 and b n + b n+1 n+1 = b n n : Lemma 6.4. Suppose f (u) is nonnegative, non-decreasing on [0, ∞) and u 0 0, then u n given by the explicit scheme (6.5) is nondecreasing.
We have the following discrete comparison principles: Proof. We only prove ' ' case while the other case is similar. It follows directly from the following induction inequality We now move on to the implicit scheme (6.6), which is given by As a result, we have (D h w) n+1 D γ c u(t n ) = f (u(t n )) = f (w n ) and the result follows from Lemma 6.5.
For implicit scheme, The claim then follows from Lemma 6.6.
Remark 6.8. The explicit schemes can be used to prove the stability and convergence of some approximation schemes for fractional PDEs and thus the convergence and existence of solutions. The implicit schemes can be used to prove positivity of solutions and to estimate the blow-up time.
6.3. Numerical simulations. For f (u) = Au p , the numerical solutions using explicit schemes (6.1), (6.2) and (6.5) never break up (i.e. u n can be computed for any n 1). The implicit scheme is more suitable for the study of blowup. If we use the implicit scheme (6.6), we look for the root of the scheme (6.8) in [u n , M ] to find u n+1 where M = (k −γ /(pAΓ(2 − γ)) 1/(p−1) . Suppose that the sequence terminates at N * and numerically we set T b (k) = N * k. It is expected that T b (k) → T b as k → 0 + .
For p = 2, the implicit scheme (6.8) can be solved exactly and therefore this allows us to compute the numerical solutions accurately and fast enough. Below, we do the numerical simulations using the implicit scheme for f (u) = u 2 by choosing k sufficiently small.
In Figure 1, we sketch two typical solution curves. Figure 1 (a) shows the solution curve with u 0 = 0.12, γ = 0.6, while Figure 1 (b) shows the solution curve with u 0 = 1.2, γ = 0.6. Comparing the blow-up time in both cases, we find clearly that small u 0 defers the blowup while large u 0 accelerates the blowup.
To investigate this issue further, in Figure 2, we plot the blow-up time versus γ, meanwhile we also plot the estimated upper and lower bounds gained from Theorem 5.3. In the case u 0 = 0.12, the line of real blow-up time around γ = 0.2 is quite steep. In the case u 0 = 1.2, the line of real blow-up time around γ = 0 is approximately equal to 0. The numerical results agree with our analysis in Section 5.2. If u 0 is big enough such that u − u 2 = u 0 has no solution, i.e. u 0 > 0.25, the blow-up time decreases as γ decreases, which means samller γ accelerates the blowup. However, if u 0 is small, i.e. u 0 < 0.25, then the blow-up time increases as γ decreases, which means samller γ defers the blowup. These observations agree with intuition that as the smaller γ is, the stronger the memory effect is.

Appendix
In this section, we restate the result in [21] regarding the growth rate of (5.1), as we mentioned in Remark 5.6. The statement is tailored to our problem, and we also present the proof for convenience. In fact, we have the following statement.