Efﬁcient Spectral Collocation Method for Tempered Fractional Differential Equations

: Transient anomalous diffusion may be modeled by a tempered fractional diffusion equation. In this paper, we present a spectral collocation method with tempered fractional Jacobi functions (TFJFs) as basis functions and obtain an efﬁcient algorithm to solve tempered-type fractional differential equations. We set up the approximation error as O ( N µ − ν ) for projection and interpolation by the TFJFs, which shows “spectral accuracy” for a certain class of functions. We derive a recurrence relation to evaluate the collocation differentiation matrix for implementing the spectral collocation algorithm. We demonstrate the effectiveness of the new method for the nonlinear initial and boundary problems, i.e., the fractional Helmholtz equation, and the fractional Burgers equation.


Introduction
Anomalous diffusion transport is often observed in nature and well described with fractional calculus (FC).Two kinds of anomalous diffusions are super diffusion and sub diffusion, according to the superlinear or sublinear growth of the variance of the variable of interest.Anomalous diffusion models, as the limit of a continuous time random walk (CTRW) with a heavy-tailed power-law probability distribution function, lead to diverging moments which can be problematic from a physical point of view.Exponentially tempering the power-law distributions is a popular way to temper the power-law distribution, which has both mathematical and technique advantages [1][2][3].The FC involves an exponentially tempered factor, referred to as tempered fractional calculus (TFC), which describes the transition between normal and anomalous diffusions (or the anomalous motion in finite time or bounded physical space).TFC models have been applied in many scientific and engineering fields, such as poro-elasticity [4], ground water hydrology [5], geophysical flows [6], etc.
The application of TFC models to realistic problems requires numerical schemes that are available to solve the tempered fractional differential (TFD) equations.The difficulty of numerically solving the TFD equations is partially caused by the tempered fractional operators involving weak singular and exponential kernels.However, there is a lot of work that contributes to it.Li and Deng [7] presented the weighted and shifted Grünwald difference for the tempered fractional diffusion equations.Wang and Li [4] presented a fast difference scheme for a tempered fractional Burgers equation in porous media.Deng and Zhang [8] provided the variational formulation and efficient implementation for solving both spacial and temporal tempered fractional problems.Chen and Deng [9] proposed a secondorder accurate numerical method for space-time tempered fractional diffusion-wave equations.Ding and Li [10] proposed a high-order algorithm for time-Caputo-tempered partial differential equations with Riesz derivatives in two spatial dimensions.Guo et al. [11] presented efficient fractional linear multistep methods for tempered fractional calculus.Cao et al. [12] studied finite difference/finite element methods for tempered time fractional advection-dispersion equations with fast evaluation of the Caputo derivative.Local discontinuous Galerkin methods were studied by Sun et al. [13] for the time-tempered fractional diffusion equation and by Safari et al. [14] for a class of time-space tempered fractional diffusion equations.Bu and Oosterlee [15] proposed a multigrid method for tempered fractional diffusion equations.Çelik and Duman [16] studied the finite element method for a symmetric tempered fractional diffusion equation.Zayernouri et al. [17] studied two classes of regular and singular tempered fractional Sturm-Liouville problems, in which a Petrov-Galerkin spectral method for solving tempered fractional ODEs is developed, and the tempered Jacobi poly-fractonomials are used as basis functions.Hanert and Piret [18] proposed a Chebyshev pseudospectral method to solve the space-time tempered fractional diffusion equation.Chen et al. [19] considered a Laguerre spectral-Galerkin method for solving tempered fractional diffusion equations on unbounded domain.Luo et al. [20] presented a Lagrange-quadratic spline optimal collocation method for the time-tempered fractional diffusion equation.Moghaddam et al. [21] suggested sinc collocation for the TFD equation, and Liemert and Kienle [22] used the Fourier spectral method to solve the tempered fractional wave-diffusion equation.
Meanwhile, Li et al. [23] discussed the existence, uniqueness, and stability of the tempered fractional ordinary differential equations.Deng et al. [24] discussed the reflecting boundaries for tempered fractional diffusion operators in the context of demonstrating the well-posedness of PDEs with generalized boundary conditions.From a mathematical view, the tempered fractional calculus operators are equivalent to the fractional substantial calculus after a suitable range of parameters [25].One can also refer to [26,27] and references therein for the numerical methods of the substantial fractional differential equations.
Motivated by the generalized Jacobi functions in [28] and the method in [17], we aim to develop a high-accuracy spectral collocation method that uses the tempered fractional Jacobi functions (TFJFs) as the basis functions for solving TFD equations in this work.The main contributions of this paper are as follows: † We define the TFJFs and derive the approximation results of orthogonal projection and interpolation based on the TFJFs.

†
We derive the differentiation matrix of the tempered fractional Caputo derivative and give a fast and stable evaluation method based on the recurrence relationship.

†
We demonstrate the effectiveness of the proposed spectral collocation method for the initial or boundary value problems, i.e., the fractional Helmholtz equation, and the fractional Burgers equation.
The outline of the paper is as follows.In the next section, we review some definitions and properties of tempered fractional calculus.We also define the TFJFs and show their properties.We present the approximation properties based on the TFJFs in Section 3. In Section 4, we present the spectral collocation method based on the TFJFs and derive a recurrence relation for the fast generation of the collocation fractional differentiation matrix.Finally, we discuss the condition number of the fractional differentiation matrix.In Section 5, we apply the spectral collocation to the initial value problems of nonlinear fractional ordinary differential equation.In Section 6, we apply the spectral collocation to the fractional Helmholtz equation and the fractional Burgers equation.Several numerical tests are presented in this section.The paper ends with some conclusions in Section 7.

Definitions of Tempered Fractional Calculus
In the sequence, we start with the definitions of tempered fractional operators.Definition 1.For κ ≥ 0, the left and right tempered fractional integrals of function u(t) on (a, b) of order µ > 0 are defined, respectively, by [5,8,29] a By direct calculation, the following relations are obtained: and where n − 1 < µ ≤ n.
If κ ≡ 0, then the left and right tempered fractional calculus reduce to the left and right (non-tempered) fractional calculus, denoted by a Similar to the corresponding non-tempered case, for 0 < µ < 1 and u(t) be absolutely continuous in [a, b], the relations between the Riemann-Liouville and Caputo derivatives are and For practical purposes, we consider the linear transform from Then, it holds that and ,

Tempered Fractional Jacobi Functions (TFJFs)
In this subsection, we introduce the TFJFs and derive some properties for use.Denote P N (I) as the space of all algebraic polynomials with degree at most N defined on I =: (−1, 1).Let P α,β n (s) (α, β > −1), s ∈ I be the Jacobi orthogonal polynomials, satisfying 1 where , and δ mn is the Kronecker symbol, i.e., The derivative relation of the Jacobi polynomials is of importance: The following relation is useful: where Some additional properties of Jacobi polynomials can be referred to [30,31].
with starting terms where P2. Orthogonality: For α, β > −1, α,β,δ,κ r where Proof.The three-term recurrence relation given by Equation ( 7) is a straightforward result from the corresponding relation of Jacobi polynomials.We derive, from Definition 4, that Then, we have the orthogonality (9).We can derive the other equalities (10) and ( 11) similarly.This ends the proof.
The tempered fractional derivatives of the TFJFs can also be represented by the same class of functions. .
Proof.From Definitions 4, we have Then, by Theorem 3.1 in [28], where the desired result can be obtained.The second equality is derived similarly.

Error Estimate of Projection Operators
We introduce the (N + 1)-dimensional spaces of the TFJFs as N,l as the orthogonal projection such that and N,r as the orthogonal projection such that Consider P δ,κ N,l .One can express the orthogonal projection as with Since the Jacobi polynomials {P α,β n (x)} n≥0 are mutually orthogonal and complete in L 2 ω α,β (I), it uniquely holds that where That is, the unique representation of u is Hence, the set {J . By the definition of the projection operator and Parseval equality, one has It is clear that a parallel statement for P δ,κ N,r holds.We do not repeat that.To describe the projection error, we define We have the following error estimate results about the projection operators defined in ( 12) and (13).
Proof.We first consider P β,κ N,l .From Equation ( 14) and the orthogonality of Lemma 1, one has Then the estimate is derived.
In a similar way, we can obtain the error estimate for P α,κ N,r .
For the special case of δ = 0, we have the following error estimate.
, where the non-uniformly Jacobi-weighted Sobolev space equipped with the inner product and norm ).Then, by utilizing the approximation result of P 0,0 N (Theorem 3.35 in page 118 of [30]), we can achieve the desired estimate result.
, by utilizing the approximation result of Π 0,0 N (Theorem 3.43, page 137 of [30]), we can obtain the desired estimate.The second estimate is derived similarly.
From the above estimate in Theorem 4, one has and This is the stability of the interpolations Π δ,κ N,l and Π δ,κ N,r .Now, we give an inverse inequality.
(s) and one has On the other hand, from Theorem 1, one has .
(s) and one has On the other hand, from Theorem 1, one has .
. This is the second inequality.
Then, we have the following estimate.
N,l u, then, from Theorems 2, 5 and the stability (17), one has Then we obtain the first estimate.Since Π α,κ N,r (P α,κ N,r u) = P α,κ N,r u, then from Theorems 2, 5 and the stability (18), one has This ends the proof.

Differentiation Matrix of the Collocation Method
Let {ξ i } N i=0 be the JGL points defined as in the previous subsection and {L i (s)} N i=0 be the Lagrange interpolation basis functions with respect to the nodes {ξ i } N i=0 , that is,

Differentiation Matrix for the Left Tempered Derivative
For the interpolation operator Π δ,κ N,l , the tempered fractional Lagrange interpolation basis functions are defined by which satisfies F δ,κ i,l (ξ k ) = δ ik .Given function u(s) ∈ C(I) satisfies e κs (1 + s) −δ u(s) being continuous on I for some δ > −1.It can be interpolated by Naturally, the tempered Caputo fractional derivative of u(s) can be approximated by its interpolation .
Now, we compute the differentiation matrix of the left tempered Caputo derivative (DMLTCD), which is denoted as s), the link between the Lagrange interpolation function L i (s) and the Jacobi polynomials should be used (Theorem 3.28 of page 87 in [30]) where , and {ω i } N i=0 are the weights corresponding with the nodes {ξ i } N i=0 in the Jacobi-Gauss-Lobatto quadrature.Then, we obtain that In the sequence, we demonstrate how to compute S δ,κ,µ n,α,β,l (s).
respectively, where by denoting f n are given in Lemma 1 and (6), respectively.
Proof.By making use of the relation (7), we have Thus, Then, we obtain the first equality.With integration by parts, and from (5), we have Hence, we have By working out S n+1,l , we obtain the desired relation for n ≥ 1.The starting terms can be derived by direct calculation.

Now, since
we can evaluate quickly and stably by Theorem 7.

Differentiation Matrix for the Right Tempered Derivative
For the interpolation operator Π δ,κ N,r , the tempered fractional Lagrange interpolation basis functions are defined by which satisfies F δ,κ i,r (ξ k ) = δ ik .Given function u(s) ∈ C(I) satisfies that e −κs (1 − s) −δ u(s) is continuous on I for some δ > −1.It can be interpolated by Naturally, the right tempered Caputo fractional derivative of u(s) can be approximated by its interpolation Now, we compute the differentiation matrix of the right tempered Caputo derivative(DMRTCD), which is denoted as Similar to the above subsection, we use ( 20) to obtain Now, denote with the starting terms S 0,r = Γ(δ + 1) respectively, where and f Proof.By making use of the relation (7), we have Thus, Then, we obtain the first equality.With integrating by parts and from (5), we have Hence, we have By working out S n+1,r , we obtain the desired relation for n ≥ 1.The starting terms can be derived by direct calculation.
The evaluation of the DMRTCD is similar to the previous part.In fact, the matrix according to the relation between the left and right tempered fractional calculus.
In order to understand the behavior of the fractional calculus, we cite the following result (see Proposition 2.1 in [32]): From the above Lemma, it is known that the first row of the DMLTCD is all zeros when δ > µ, and all infinity when δ < µ, so is the last row of the DMRTCD.In practice, the first row and column of the DMCHD (the DMRTCD) are removed according to initial or boundary value conditions.The DMLTCD (the DMRTCD) is full for every µ > 0, and the condition number of the DMLTCD (the DMRTCD) is like the first-or second-order differentiation matrix in the collocation method, as usual.We compute the condition numbers for different values of µ and show those in Figure 1 (for 0 < µ < 1) and Figure 2 (for 1 < µ < 2).The behavior of the condition number of the DMLTCD (the DMRTCD) is like O(N 2µ ) from numerical tests.
Before providing more numerical results of tests, I would like to mention that all codes were compiled by using MatLab R2019a on pc with AMD Ryzen 7, 5800H with Radeon Graphics 3.20 GHz, RAM 16 G.

Applications to Nonlinear Initial Value Problems
Let 0 < µ < 1.In this section, we apply the spectral collocation method to the initial value problem of nonlinear fractional ordinary differential equation below: Hereafter, we assume that f (x, u) is continuous in the given domain and satisfies the Lipschitz condition with respect to the second variable, which guarantees the existence and uniqueness of the solution.
The spectral collocation method based on the TFJFs for ( 28) is to find and The above two equations lead to the following linear system where D µ, κ s,l is the submatrix of D µ, κ s,l , which deletes the first row and column of D We need to solve a nonlinear system (31) when f is nonlinear with respect to u.The iteration methods can be used, e.g., Newton method.To measure the accuracy of the method, we define the errors by where u(x) and u N (x) are the exact and numerical solutions, respectively.Example 1.Consider (28) and f (x, u) = g(x) − u 2 , where g(x) is determined by The following two cases of the exact solutions are considered [33]: x µ with = 1.Clearly, the homogeneous initial condition u(0) = 0 is fulfilled.In this example, the Newton method is applied to solving the nonlinear system (31), which takes the following form: Since u(x) ∈ F 0,κ 2 , with two degrees of freedom (N = 2, δ = 0), the method solves the problem C1 exactly without consideration of nonlinear approximation.It is noticed that the predictor-corrector method to solve the same problem in [33] achieves accuracy 10 −4 by h = 1/160, whilst our method achieves an accuracy 10 −16 by N = 2 (the iterative number is less than 8 in the Newton method).
For the exact solution u / ∈ F δ,κ N for any δ > −1 in C2, we test different N, δ and apply the Newton method with maximum iterative step 20 to solve the nonlinear system.We list the error E 0 by N = 10, 20, 40, 80, 160 and δ = µ − 1 + 5 * 10 −11 for the case C2 in Table 1.We also list the results by the predictor-corrector method to solve the same problem in [33].It is shown our method is more superior.Table 1.The numerical errors E 0 of Example 1 with C2 solution (α = 0, β = 0, κ = 1).In this subsection, we apply the spectral collocation method to the following fractional Helmholtz equation:

N(1/h)
The spectral collocation method based on the TFJFs for ( 32) is to find and The above two equations lead to the following linear system: where I is the unitary matrix, D To measure the accuracy of the method, we define the errors by where u(x) and u N (x) are exact and numerical solutions, respectively.To understand the convergence behavior of the method, we define the order of convergence by For the smooth solution C1, we solve Equation ( 32) by δ = 0. We list the errors E 0 in Tables 2-4 for different values of parameters µ, α, β, κ and λ, which show that spectral accuracy can be achieved for a large range of parameters.It is observed that there is no significant difference between the Chebyshev method (α = β = −0.5)and the Legendre method (α = β = 0) or other Jacobi methods by precision (see Table 3).The exact solution has low regularity when σ is a non-integer in C2, so we expect the order of convergence of our method to be finite.We list the errors E 0 and the order of convergence Ord in Tables 5-8.A possible guess of the Ord for the case C2 with δ = 0 is O(N −2(σ−µ+1) ) from the results in Tables 5 and 7. When δ = σ − 1, the super convergence can be observed; see the results in Tables 6 and 8.
For time advancing, we employ a semi-implicit time-discretization scheme, namely, the two-step second-order Crank-Nicolson/leapfrog scheme.Then, the full discretization scheme reads as where D is the first-order differentiation matrix and D = ( Example 3. Consider (36) with the initial profile as one of the following: In the example, we always take α = β = 0, τ = 10 −3 .
We first consider the initial profile with two peaks, i.e., case C1.We show the numerical solutions of the FBE at time t = 1 in Figure 3 for different values N = 50, 100, 300 by µ = 1.18, = 1 and κ = 0, 1 (κ = 0 for non-tempered and κ = 1 for tempered case).It is observed that the numerical solution is in good agreement with each other.The surface of the numerical solution for µ = 1.4,κ = 1, = 1 is plotted in Figure 4 with N = 200.The evolution of the numerical solution is observed.The numerical solutions of the FBE at time t = 1 are plotted for different values of fractional-order µ = 1.2, 1.3, 1.5, 1.7, 1.9 and κ = 1, = 1 in Figure 5 and for different values of tempered factor κ = 0, 0.5, 1, 1.5, 2 and µ = 1.55, = 1 in Figure 6.We observed that the tempered diffusion becomes slower with larger κ.
We consider the initial profile with single peak of case C2.We show the surface of the numerical solution for µ = 1.4,κ = 1, = 1 in Figure 7 and the evolution of the numerical solution.The numerical solutions of the FBE at time t = 1 are plotted in Figure 8 for different values of fractional-order µ = 1.2, 1.3, 1.5, 1.7, 1.9 where κ = 1, = 1 and in Figure 9 for different values of tempered factor κ = 0, 0.5, 1, 1.5, 2 where µ = 1.55, = 1 .It is shown that our method is a powerful tool to simulate the fractional-order diffusion of the Caputo derivative as well as of the tempered case.

Conclusions
Tempered fractional calculus is of importance in the family of fractional calculi with potential applications.In this paper, we present a spectral collocation method using the TFJFs as basis functions and obtain an efficient algorithm to solve tempered fractional differential equations.The key in implementing it is to stably evaluate the collocation differentiation matrix by utilizing a recurrence relation.The method is the direct collocation method in strong form and shares spectral accuracy when the solutions of FDEs are smooth from the approximate properties of the TFJFs.The effectiveness of the derived method is confirmed by numerical tests for the nonlinear initial problems, the fractional Helmholtz equation, and the fractional Burgers equation with the tempered Caputo derivative.
The method can be applied to other kinds of tempered fractional linear or nonlinear differential equations and even for the variable-order case.On the other hand, the method can be generalized to the multi-domain case, such as that in [34].We will develop its further applications in the future.Meanwhile, the symmetries of the FD equations are also an interesting topic [35] that can be used to construct a structure-preserving high-order algorithm.We expect to consider it in the future.
The present method is not able to resolve directly the fractional derivative of more generalized fractal fractional derivatives, such as that in [36].We will improve our method to tackle it in the future.

Definition 3 .
For κ ≥ 0, the left and right tempered Caputo fractional derivatives of function u(t) on (a, b) of order µ > 0 are defined by[8,23]

n
, B n and a n , b n , c n are given in Theorem 7.
is the inner part of D µ, κ s,l , i.e., without the last row and column of D µ, κ s,l , f = f (x) + u a d 0 + u b d N , d 0 and d N the first and last columns of D µ, κ s,l .

Table 3 .
Errors E 0 of Example 2 with solution C1

Table 5 .
Errors E 0 and Ord of Example 2 with solution C2

Table 7 .
Errors E 0 and Ord of Example 2 with solution C2