WELL-POSEDNESS AND NUMERICAL ALGORITHM FOR THE TEMPERED FRACTIONAL DIFFERENTIAL EQUATIONS

. Trapped dynamics widely appears in nature, e.g., the motion of particles in viscous cytoplasm. The famous continuous time random walk (CTRW) model with power law waiting time distribution ( having diverging ﬁrst moment ) describes this phenomenon. Because of the ﬁnite lifetime of biological particles, sometimes it is necessary to temper the power law measure such that the waiting time measure has convergent ﬁrst moment. Then the time operator of the Fokker-Planck equation corresponding to the CTRW model with tempered waiting time measure is the so-called tempered fractional derivative. This paper focus on discussing the properties of the time tempered fractional derivative, and studying the well-posedness and the Jacobi-predictor- corrector algorithm for the tempered fractional ordinary diﬀerential equation. By adjusting the parameter of the proposed algorithm, high convergence order can be obtained and the computational cost linearly increases with time. The numerical results show that our algorithm converges with order N I , where N I is the number of interpolating points used in the scheme.

1. Introduction. The fractional calculus has a long history. The origin of fractional calculus can be traced back to the letter between Leibniz and L'Hôpital in 1695. In the past three centuries, the development of the theories of fractional calculus is well contributed by many mathematicians and physicists. And from the last century, the books covering fractional calculus began to emerge, such as Oldham and Spanier (1974), Samko, Kilbas and Marichev (1993), Podlubny (1999), and so on. In recent years, more theories and experiments show that a broad range of non-classical phenomena appeared in the applied sciences and engineering can be described by fractional calculus [36,29,38]. Because of its good mathematical features, nowadays fractional calculus has become a powerful tool in depicting the 1990 CAN LI, WEIHUA DENG AND LIJING ZHAO anomalous kinetics which arises in physics, chemistry, biology, finance, and other complex dynamics [29]. In practical applications, several different kinds of fractional derivatives, such as Riemann-Liouville fractional derivative, Caputo fractional derivative [36,38], Riesz fractional derivative [38], and Hilfer fractional derivative [22,46] are introduced.
One of the typical applications for fractional calculus is the description of anomalous diffusion behavior of living particles; and the tempered fractional calculus describes the transition between normal and anomalous diffusions (or the anomalous diffusion in finite time or bounded physical space). In the continuous time random walk (CTRW) model, for a Lévy flight particle, the scaling limit of the CTRW with a jump distribution function φ(x) ∼ x −(1+α) (1 < α < 2) exhibits superdiffusive dynamics. The corresponding stable Lévy distribution for particle displacement contains arbitrarily large jumps and has divergent spatial moments. However, the infinite spatial moments may not be feasible for some physical processes [10]. One way to overcome the divergence of the moments of Lévy distributions in transport models is to exponentially temper the Lévy measure. Then the space fractional operator will be replaced by the spatially tempered fractional operator in the corresponding models [9,10,39]. This paper concentrates on the time tempered fractional derivative, which arises in the Fokker-Planck equation corresponding to the CTRW model with tempered power law waiting time distribution [37,42,21,19]. Tempering the power law waiting time measure makes its first moment finite and the trapped dynamics more physical. Sometimes it is necessary/reasonable to make the first moment of the waiting time measure finite, e.g., the biological particles moving in viscous cytoplasm and displaying trapped dynamical behavior just have finite lifetime. The time tempered diffusion dynamics describes the coexistence/transition of subdiffusion and normal diffusion phenomenon (or the subdiffusion in finite time) which was empirically confirmed in a number of systems [10,30]. More applications for the tempered fractional derivatives and tempered differential equations can be found, for instance, in poroelasticity [20], finance [9], ground water hydrology [30,31], geophysical flows [32] and unsteady convection-diffusion flows [2].
Tempered fractional calculus can be recognized as the generalization of fractional calculus. To the best of our knowledge, the definitions of fractional integration with weak singular and exponential kernels were firstly reported in Buschman's earlier work [5]. For the other different definitions of the tempered fractional integration, see the books [43,38,31] and references therein. Two new spectral methods on tempered fractional Sturm-Liouville problems have been recently developed by Zayernouri et al. in [49]. This work continues previous efforts [28] to explore the time tempered fractional derivative. The well-posedness, including existence, uniqueness, and stability, of the tempered fractional ordinary differential equation (ODE) is discussed, and the properties of the time tempered fractional derivative are analyzed. Then the Jacobi-predictor-corrector algorithm for the tempered fractional ODE is provided, which has the striking benefits: any desired convergence order can be obtained by simply adjusting the parameter (the number of interpolation points); the computational cost increases linearly with the time t instead of t 2 usually taken place for nonlocal time dependent problem. And extensive numerical experiments are performed to confirm these advantages.
In Section 2, we introduce the definitions and show the properties of the tempered fractional calculus, including the generalizations of the tempered fractional derivatives in the Riemann-Liouville and Caputo sense, and the composite property. More basic properties are listed and proved in Appendix A; the expression and properties of the tempered fractional calculus in Laplace space are proposed and proved in Appendix B. In Section 3, we discuss the initial value problem of the tempered fractional ODE: first derive the Volterra integral formulation of the tempered fractional ODEs; then prove the well-posedness of the considered problem. The Jacobi-Predictor-Corrector algorithm for the tempered fractional ODE is designed and discussed in Section 4, and two numerical examples are solved by the algorithm to show its powerfulness.
And denote by C n [a, b] the space of functions u(t) which are n times continuously differentiable on [a, b].
Proof. Take v(t) = e λt u(t) in the equation for the Riemann-Liouville and Caputo fractional derivatives [36,26,38] Multiplying both sides of the above equation by e −λt , we obtain Furthermore, using the definitions of Riemann-Liouville and Caputo tempered fractional derivatives, we get that Using the linearity properties presented in Proposition 4 and the formula of power function we deduce the desired result from (11).
Proof. From the definitions of Riemann-Liouville tempered fractional integral and derivative, we have Thanks to the composition formula [36,26,38] we get Inserting the above formula into (16) leads to (12). Again from the definitions of Riemann-Liouville tempered fractional integral and derivative, there exists = e −λt a D α t e λt e −λt a I α t (e λt u(t)) = e −λt a D α t a I α t (e λt u(t)) . Furthermore, using the composite properties of Riemann-Liouville fractional integral and derivative [36,26,38] we get a D α,λ t [ a I α,λ t u(t)] = e −λt a D α t a I α t (e λt u(t)) = u(t), by taking v(t) = e λt u(t) in Eq. (17).
Obviously, a D α,λ 3. Well-posedness of the tempered fractional ordinary differential equations. In this section, we consider the ODEs with Riemann-Liouville and Caputo tempered fractional derivatives, respectively, i.e., and The Cauchy problems (20) and (21) can be converted to the equivalent Volterra integral equations of the second kind under some proper conditions.
In particular, if 0 < α < 1, then u(t) satisfies the Cauchy problem (20) if and only if u(t) satisfies the following integral equation Proof. For the linear Cauchy problems (20) and (21), the conclusion is directly reached by the Laplace transform given in Appendix B. Now we prove the more general case. N ecessity. Performing the integral operator a I α,λ t on both sides of the first equation of (20), we have where we use the composite property (1) given in Proposition 2. Then Eq. (22) is obtained.
Suf f iciency. Applying the operator a D α,λ t to both sides of Eq.
where we use the fact and the composite property (13). Now we show that the solution of (22) satisfies the initial condition given in Eq. (20). Multiplying e λt and then performing the where the formula is utilized. Taking a limit t → a in the above equation, we obtain with the second term in the right hand side being equal to zero; and for its first term, we have By the similar technique in proving Lemma 3.1, we obtain the following conclusion for the Cauchy problem (21).
In particular, if 0 < α < 1, then u(t) satisfies the Cauchy problem if and only if u(t) satisfies the following integral equation 3.1. Existence and uniqueness. Many authors have considered the existence and uniqueness of the solutions to the nonlinear ODEs with fractional derivatives; see, e.g., [35,1,17,16,15,23,24,46,50]. For the global existence and uniform asymptotic stability results of fractional functional differential equations corresponding to (23), one can see [27,3]. In the following, we discuss the existence and uniqueness of the solutions of the nonlinear tempered fractional differential equations based on the equivalent Volterra equations presented in Lemmas 3.1 and 3.2. We shall employ the Banach fixed point theorem to prove it. Let f : being an open set in R. In the following, we always suppose that f (t, u) satisfies the Lipschitz type condition with respect to the second variable where C Lip is constant. We shall use the following space Proof. The proof of this theorem is similar to the references [35,16,7,23,46]. For the reader's convenience, we provide the proof in Appendix C.
By almost the same idea, we can prove the following existence and uniqueness result for the Cauchy type problem (21).

Stability.
To prove the stability of the solutions of the Cauchy type problems (20) and (21), we need the following generalized Gronwall's Lemmas.
, and Then holds If, in addition, y(·) is a nondecreasing function defined on [a, b], we have Then there exists a constant C = C(α) such that Theorem 3.7. Under the assumptions given in Theorem 3.3, let u(t) and v(t) be the solutions of the Cauchy type problem (20) with different initial conditions. Then where Similarly, under the assumptions given in Theorem 3.4, for the problem (21), there exists where Proof. Suppose that u(t) and v(t) are any two solutions of the Cauchy type problem (20) with different initial conditions. From the equivalent integral formulation (22), we have For 0 < α < 1, using the generalized Gronwall's inequality with weak singular kernel given in Lemma 3.6, we get which implies For n − 1 < α < n, n ≥ 2, taking where C(α) = max k=0,1,··· ,n−1 . With the similar method, we can prove the stability results for the problem (21).
Proof. Similar to Theorems 3.3 and 3.4. We begin our proof from the integral equations given in Lemma 3.8. We only prove the Cauchy type problem (42). Let t 1 belongs to (a, b) such that the inequality holds.The operator corresponding to (44) takes the form where From the Lipschitz condition (45) it directly follows that Furthermore, using the composition formula (12), we have where n j is the smallest integer larger than or equal to α j . From the given initial conditions, it can be checked that a D αj −kj −1 t (e λt (u(t) − v(t))) t=a = 0. Then for any t ∈ [a, b]. Taking t = t 1 in above the formula and applying (A.69), we get .
It follows that there exists a unique solution u * (t) to Eq. (44) in L([a, t 1 ]). This solution is obtained as a limit of the convergent sequence (T j u * 0 )(t) = u j (t), and holds lim i.e., lim With the same fashion of proving Theorem 3.3, we can show that there exists an unique solution u(t) ∈ L α,λ ([a, b]) to Eq. (44). In addition, which implies that a D α,λ t u(t) ∈ L([a, b]).

4.
Numerical algorithm for the tempered fractional ODE. The well-posedness of the tempered fractional ODE has been carefully discussed in the above sections. Usually, it is hard to find the analytical solutions of the tempered fractional ODE, especially for the nonlinear case. Efficient numerical algorithm naturally becomes an urgent topic for this type of equation. Now, we extend the so-called Jacobi-predictor-corrector algorithm [48] to the the tempered fractional ODE; and its striking benefits are still kept, including having any desired convergence orders and the linearly increasing computational cost with the time t.

CAN LI, WEIHUA DENG AND LIJING ZHAO
Firstly, we transform the original equation into whereg Using (N + 1)-point Jacobi-Gauss-Lobatto quadrature to approximate the integral in (52) yields where we choose ω(z) = (1 − z) α−1 (1 + z) 0 as the weight function; {z j } N j=0 and {ω j } N j=0 are the (N +1)-degree Jacobi-Gauss-Lobatto nodes and their corresponding Jacobi weights in the reference interval [−1, 1], respectively; see, e.g., [41]. Now we turn to describe the computational scheme for Eq. (52). For this purpose, we define a grid in the interval [a, b] with M + 1 equidistant nodes t j , given by where τ = (b − a)/M is the stepsize. Suppose that we have got the numerical values of u(t) at t 0 , t 1 , · · · , t n , which are denoted as u 0 , u 1 , · · · , u n ; now we are going to compute the value of u(t) at t n+1 , i.e., u n+1 . From Eq. (53), we have wherẽ To compute the second summation term of (55), we need to evaluate the value of f at the point t n+1 due tof n+1 z N ,ũ n+1 (z N ) = f t n+1 , u(t n+1 ) . It can be numerically approximated by using the piecewise linear interpolation of the term f [16,15,12]. Here, we do it using the technique given in our previous work [48]. More concretely, we interpolate the function f by the known "neighborhood" points of t n+1 . For the other values off n+1 z j ,ũ n+1 (z j ) , 0 ≤ j ≤ N − 1, we can also obtain them based on the interpolation of f at the time nodes located in the "neighborhood" points of z j (should be (1 + z j )t n+1 /2 as to variable t). Denote N I as the number of time nodes used for the interpolation. In practical application, we use the improved predictor-corrector formulas given in [12] to get the values at t 0 , t 1 , · · · , t N I −1 as the known 'initial' values. For the criterion of choosing the "neighborhood" points, refer to [48].
Collecting the above analysis, we get the predictor-corrector formulas of (51) as and where {f P n+1,j } N j=0 in (57)  whereas {f n+1,j } N −1 j=0 in (56) are obtained by using the interpolations based on the values of {f (t j , u j )} n j=0 and f (t n+1 , u P n+1 ). From the computational scheme (56)-(57), it can be clearly seen that the computational cost linearly increase with n (or time t).
With the similar methods given in [48], we can get the following estimates for the Volterra integral system (51).
Theorem 4.1. If g(t, u(t)) is Lipschitz continuous with respect to the second variable, and has the form where w k (t), 1 ≤ N I < µ 1 < ... < µ m are sufficiently smooth and m can be +∞, δ is a constant, then there exists a constant C being independent of n, τ, N , such that This theorem shows that the scheme (56)-(57) potentially have any desired convergence order by adjusting the number of interpolation points N I .

Numerical test.
In this subsection, we consider two simple numerical examples to show the numerical errors and convergence orders of the Jacobi-predictorcorrector method. The two examples are Caputo tempered ODEs; solving the Riemann-Liouville tempered ODEs can be done in the same way, so is omitted here.
error order error order error order 1/10 6.6386e-005 9.6009e-006 8.5068e-007 1/20 9.2847e-007 6.1599 1.4297e-007 6.0694 1.9943e-008 5.4147 1/40 1.5767e-008 5.8799 2.1338e-009 6.0661 3.0437e-010 6.0339 1/80 2.3505e-010 6.0678 3.5138e-011 5.9242 3.8203e-012 6.3159 1/160 3.8498e-012 5.9320 5.3434e-013 6.0391 6.7433e-014 5.8241 Then Employing the Laplace transform involving the derivative of the Mittag-Leffler function [36] L we can check that the exact solution of this initial value problem is Here the generalized Mittag-Leffler function E α,β (·) is given by [36] In this example, the solution u(t) does not have a bounded first (second) derivative at the initial time t = 0 as 0 < α < 1 (1 < α < 2). To improve the convergence order, we employ the technique given in our previous work [48]. We separately solve the equation in subintervals [0, T 0 ] and [T 0 , T ] of the interval [0, T ]. More specifically, we modify the formula (52) as Here, we suppose that the smoothness of g is weaker on the subinterval [0, whereÑ , e −λt {ω j }Ñ j=0 and {s j }Ñ j=0 correspond to the number of, the weights of, and the values of the Gauss-Lobatto nodes with the weight ω(s) = 1 in the interval [0, T 0 ], respectively. The values of g s j , v(s j ) Ñ j=0 can be computed as in the starting procedure. Since g and x are continuous in the interval [0, T 0 ], by the theory of Gauss quadrature [34] and the analysis above, we can see that ifÑ is a big number then the accuracy of the total error can still be remained. The numerical results are reported in Tables 4 and 5. And it can be seen that the desired numerical accuracy is obtained. Example 3. In the last example, we use our algorithm to solve the nonlinear problem C 0 D α,λ t u(t) = e −λt Γ(6) Γ(6 − α) t 5−α − e λt t 10 + e λt u 2 (t) , 0 < α < 1, with the initial value chosen as u(t)| t=0 = 0. Using the formula (59), we can check that the exact solution of this initial value problem is u(t) = e −λt t 5 .  Table 6 shows the maximum errors and convergence orders of Example 3 solved by the scheme (56)-(57). It can be seen that the convergence rates of the errors in maximum norm is almost of the order of 5, confirming the convergence rates given in Theorem 4.1.
5. Concluding remarks. Currently, it is widely recognized that fractional calculus is a powerful tool in describing anomalous diffusion. Because of the bounded physical space and the finite lifetime of living particles, in the CTRW model, sometimes it is necessary to truncate (temper) the measures of jump length and waiting time with divergent second and first moments, respectively. Exponential tempering offers technique advantages of remaining the infinitely divisible Lévy process at the same time the transition densities can be computed at any scale. Then the Fokker-Planck equation of the corresponding CTRW model has the time tempered fractional derivative (tempered waiting time) and/or the space tempered fractional derivative (tempered jump length). This paper focus on discussing the properties of the time tempered fractional derivatives as well as the well-posedness and numerical algorithm for the time tempered evolution equation, i.e., the tempered fractional ODEs. The proposed so-called Jacobi-predictor-corrector algorithm shows its powerfulness/advantages in solving the tempered fractional ODEs, one of which is easily getting any desired convergence orders by simply changing the parameter of the number of the interpolating points and the other one of which is linearly increasing computational cost with time t comparing with the quadratically increasing that mostly happen when numerically solving fractional evolution equations. Appendix A: Properties of tempered fractional integral and derivatives. For the Riemann-Liouville tempered fractional integral, we have the following properties. Property 1. If u(t) ∈ AC n [a, b], σ > 0, and λ ∈ R, then for all t ≥ a, ds n e λs u(s) ds.
Proof. If u(t) has continuous derivative for t ≥ a, then using integration by parts to (1), there exists and if the function has n continuous derivatives, then integrating by parts leads to Lemma .1. For a function u(t) ∈ L([a, b]), σ > 0, λ ∈ R, we have Proof. By simple calculation, we have Similar to the fractional calculus, the tempered fractional calculus is a linear operation. (1) for σ > 0, holds . Proof. The linearity of fractional integral and derivatives follows directly from the corresponding definitions. We omit the details here. Proof. After simple argument, we have From the above analysis, we get the desired result.
Appendix B: Laplace transforms of the tempered fractional calculus.
In this subsection, we discuss the Laplace transforms of the tempered fractional calculus. Define the Laplace transform of a function u(t) and its inverse as [40] L{u(t); s} = u(s) = We start with the Laplace transform of the Riemann-Liouville tempered fractional integral of order σ.
Proposition 6. The Laplace transform of the Riemann-Liouville tempered fractional integral is given by Proof. First, we rewrite the Riemann-Liouville tempered fractional integral as the form of convolution In view of the Laplace transform of the convolution [40] L{u(t) * v(t); s} = u(s) v(s), we have Recalling the Laplace transform L{e −λt t σ−1 ; s} = Γ(σ)(λ + s) −σ , we have the Laplace transform of Riemann-Liouville tempered fractional integral L{ 0 I σ,λ t u(t); s} = (λ + s) −σ u(s). (B.74) Next, we turn to consider the Laplace transform of tempered fractional derivative.
Proposition 7. The Laplace transform of the Riemann-Liouville tempered fractional derivative is given by The Laplace transform of the Caputo tempered fractional derivative is given by Proof. Using the properties of Riemann-Liouville tempered fractional calculus, we may rewrite the Riemann-Liouville tempered fractional derivative as Combining the first translation Theorem [40] L{e −λt u(t); s} = u(λ + s), λ ∈ R, Re(s) > λ, We first prove that P is a contraction operator in the subinterval [a, t 1 ] ⊂ [a, b] (a < t 1 < b). To do this, we select t 1 ∈ ([a, b]) such that the inequality (ii) For all u 1 , u 2 ∈ L([a, t 1 ]) the following inequality holds P u 1 − P u 2 L(a,b) ≤ W 1 u 1 − u 2 L([a,b]) , W 1 = C Lip (t 1 − a) α Γ(α + 1) . (C.91) In fact, since f (t, u(t)) ∈ L([a, t 1 ]) and Lemma .1 in Appendix A, the integral in the right-hand side of (C.88) belongs to L([a, t 1 ]); obviously u 0 (t) ∈ L([a, t 1 ]), hence (P u)(t) ∈ L([a, t 1 ]). Now, we prove the estimate (C.91). From Lemma .1 in Appendix A, we obtain is the known function. With the same contraction factor W 1 , we can prove that there exists a unique solution u * (t) ∈ L(t 1 , t 2 ) to Eq. (22) on the interval [t 1 , t 2 ]. By repeating this process finite times, e.g., M times, we can cover the whole interval [a, b].
To complete the proof of the theorem we must show that such a unique solution u(t) ∈ L([a, b]) belongs to the space L α,λ ([a, b]). It is sufficient to prove that