Next Article in Journal
Chaotic Characteristic Analysis of Dynamic Gravity Model with Fractal Structures via an Improved Conical Volume-Delay Function
Next Article in Special Issue
Existence of Sobolev-Type Hilfer Fractional Neutral Stochastic Evolution Hemivariational Inequalities and Optimal Controls
Previous Article in Journal
Photothermal Response for the Thermoelastic Bending Effect Considering Dissipating Effects by Means of Fractional Dual-Phase-Lag Theory
Previous Article in Special Issue
Solitary Wave Solutions to a Fractional Model Using the Improved Modified Extended Tanh-Function Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Efficient Spectral Collocation Method for Tempered Fractional Differential Equations

School of Mathematics and Statistics, Shandong University of Technology, Zibo 255049, China
Fractal Fract. 2023, 7(3), 277; https://doi.org/10.3390/fractalfract7030277
Submission received: 3 February 2023 / Revised: 17 March 2023 / Accepted: 18 March 2023 / Published: 22 March 2023

Abstract

:
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 efficient 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.

1. 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 second-order 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.

2. Preliminaries

2.1. 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 I t μ , κ u ( t ) : = 1 Γ ( μ ) a t e κ ( t s ) u ( s ) ( t s ) 1 μ d s , t I b μ , κ u ( t ) : = 1 Γ ( μ ) t b e κ ( s t ) u ( s ) ( s t ) 1 μ d s ,
where Γ ( · ) is the Euler gamma function.
Definition 2.
For κ 0 , the left and right tempered Riemman–Liouville fractional derivatives of function u ( t ) on ( a , b ) of order μ > 0 are defined, respectively, by [5,8,29]
a R D t μ , κ u ( t ) : = e κ t Γ ( n μ ) d n d t n a t e κ s u ( s ) ( t s ) μ n + 1 d s ,
and
t R D b μ , κ u ( t ) : = e κ t Γ ( n μ ) d n d t n t b e κ s u ( s ) ( s t ) μ n + 1 d s ,
where n 1 < μ n .
By direct calculation, the following relations are obtained:
a R D t μ , κ u ( t ) : = d d t + κ n a I t n μ , κ u ( t ) ,
and
t R D b μ , κ u ( t ) : = ( 1 ) n d d t κ n t I b n μ , κ u ( t ) ,
where d d t ± κ n = k = 0 n k ! ( n k ) ! n ! ( ± κ ) n k d k d t k .
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]
a C D t μ , κ u ( t ) : = a I t n μ , κ d d t + κ n u ( t ) = e κ t Γ ( n μ ) a t e κ s u ( s ) ( n ) ( t s ) μ n + 1 d s ,
and
t C D b μ , κ u ( t ) : = t I b n μ , κ ( 1 ) n d d t κ n u ( t ) = ( 1 ) n e κ t Γ ( n μ ) t b e κ s u ( s ) ( n ) ( s t ) μ n + 1 d s ,
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 I t μ , t I b μ , a R D t μ , t R D b μ , a C D t μ and t C D b μ , respectively. Then, we have
a I t μ , κ u ( t ) = e κ t a I t μ ( e κ t u ( t ) ) , t I b μ , κ u ( t ) = e κ t t I b μ ( e κ t u ( t ) ) , a R D t μ , κ u ( t ) = e κ t a R D t μ ( e κ t u ( t ) ) , t R D b μ , κ u ( t ) = e κ t t R D b μ ( e κ t u ( t ) ) , a C D t μ , κ u ( t ) = e κ t a C D t μ ( e κ t u ( t ) ) , t C D b μ , κ u ( t ) = e κ t t C D b μ ( e κ t u ( t ) ) .
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
a R D t μ , κ u ( t ) = a C D t μ , κ u ( t ) + 1 Γ ( 1 μ ) e κ t ( t a ) μ u ( a ) ,
and
t R D b μ , κ u ( t ) = t C D b μ , κ u ( t ) + 1 Γ ( 1 μ ) e κ t ( b t ) μ u ( b ) .
For practical purposes, we consider the linear transform from [ a , b ] to [ 1 , 1 ] as
s = 2 ( t a ) b a 1 , t [ a , b ] ( o r t = b a 2 ( s + 1 ) + a , s [ 1 , 1 ] ) .
Then, it holds that
a I t μ , κ u ( t ) = b a 2 μ 1 I s μ , κ ˜ u ˜ ( s ) , t I b μ , κ u ( t ) = b a 2 μ s I 1 μ , κ ˜ u ˜ ( s ) ,
and
a C D t μ , κ u ( t ) = 2 b a μ 1 C D s μ , κ ˜ u ˜ ( s ) , t C D b μ , κ u ( t ) = 2 b a μ s C D 1 μ , κ ˜ u ˜ ( s ) ,
where κ ˜ = b a 2 κ , u ˜ ( s ) = u ( t ) = u ( b a 2 ( s + 1 ) + a ) .

2.2. 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 1 P n α , β ( s ) P m α , β ( s ) ω α , β ( s ) d s = γ n α , β δ m n ,
where ω α , β ( s ) = ( 1 s ) α ( 1 + s ) β ,
γ n α , β = 2 α + β + 1 Γ ( n + α + 1 ) Γ ( n + β + 1 ) ( 2 n + α + β + 1 ) n ! Γ ( n + α + β + 1 ) ,
and δ m n is the Kronecker symbol, i.e.,
δ m n = 1 , if m = n , 0 , otherwise .
The derivative relation of the Jacobi polynomials is of importance:
d k d s k P n α , β ( s ) = d n , k α , β P n k α + k , β + k ( s ) , d n , k α , β = Γ ( n + k + α + β + 1 ) 2 k Γ ( n + α + β + 1 ) , n k .
The following relation is useful:
P n α , β ( s ) = d d s A ^ n α , β P n 1 α , β ( s ) + B ^ n α , β P n α , β ( s ) + C ^ n α , β P n + 1 α , β ( s ) ,
where
A ^ n α , β = 2 ( n + α ) ( n + β ) ( n + α + β ) ( 2 n + α + β ) ( 2 n + α + β + 1 ) , B ^ n α , β = 2 ( α β ) ( 2 n + α + β ) ( 2 n + α + β + 2 ) , C ^ n α , β = 2 ( n + α + β + 1 ) ( 2 n + α + β + 1 ) ( 2 n + α + β + 2 ) .
Some additional properties of Jacobi polynomials can be referred to [30,31].
Definition 4.
Define the TFJFs by
J n , l α , β , δ , κ ( s ) = : e κ s ( 1 + s ) δ P n α , β ( s ) , J n , r α , β , δ , κ ( s ) = : e κ s ( 1 s ) δ P n α , β ( s ) ,
for s I , n = 0 , 1 , .
From the properties of Jacobi polynomials and Definition 4, we can easily derive the following lemma:
Lemma 1.
The TFJFs have the following properties.
P1. 
Three-term recurrence relation: For n 1 ,
J n + 1 , l α , β , δ , κ ( s ) = A n α , β s B n α , β J n , l α , β , δ , κ ( s ) C n α , β J n 1 , l α , β , δ , κ ( s ) , J n + 1 , r α , β , δ , κ ( s ) = A n α , β s B n α , β J n , r α , β , δ , κ ( s ) C n α , β J n 1 , r α , β , δ , κ ( s ) ,
with starting terms
J 0 , l α , β , δ , κ ( s ) = e κ s ( 1 + s ) δ , J 1 , l α , β , δ , κ ( s ) = e κ s ( 1 + s ) δ α + β + 2 2 s + α β 2 , J 0 , r α , β , δ , κ ( s ) = e κ s ( 1 s ) δ , J 1 , r α , β , δ , κ ( s ) = e κ s ( 1 s ) δ α + β + 2 2 s + α β 2 ,
where
A n α , β = ( 2 n + α + β + 1 ) ( 2 n + α + β + 2 ) 2 ( n + 1 ) ( n + α + β + 1 ) , B n α , β = ( 2 n + α + β + 1 ) ( 2 n + α + β + 2 ) 2 ( n + 1 ) ( n + α + β + 1 ) , C n α , β = ( n + α ) ( n + β ) ( 2 n + α + β + 2 ) ( n + 1 ) ( n + α + β + 1 ) ( 2 n + α + β ) .
P2. 
Orthogonality: For α , β > 1 ,
I J n , l α , β , δ , κ ( s ) J m , l α , β , δ , κ ( x ) ϖ l α , β , δ , κ ( s ) d s = γ n α , β δ n m ,
I J n , r α , β , δ , κ ( s ) J m , r α , β , δ , κ ( x ) ϖ r α , β , δ , κ ( s ) d s = γ n α , β δ n m ,
I J n , l α , β , δ , κ ( s ) J m , r α , β , δ , κ ( x ) ϖ h α , β , δ ( s ) d s = γ n α , β δ n m ,
where
ϖ l α , β , δ , κ ( s ) = e 2 κ s ω α , β 2 δ ( s ) = e 2 κ s ( 1 s ) α ( 1 + s ) β 2 δ , ϖ r α , β , δ , κ ( s ) = e 2 κ s ω α 2 δ , β ( s ) = e 2 κ s ( 1 s ) α 2 δ ( 1 + s ) β , ϖ h α , β , δ , κ ( s ) = ω α δ , β δ ( s ) = ( 1 s ) α δ ( 1 + s ) β δ .
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
I J n , l α , β , δ , κ ( s ) J m , l α , β , δ , κ ( s ) ϖ l α , β , δ , κ ( s ) d s = I P n α , β ( s ) P m α , β ( s ) ω α , β ( s ) d s .
Then, we have the orthogonality (9). We can derive the other equalities (10) and (11) similarly. This ends the proof. □
Remark 1.
The orthogonalities (9) and (10) are in the inner product of space L ϖ l α , β , δ , κ ( I ) and L ϖ r α , β , δ , κ ( I ) , respectively. The relationship (11) means that elements between the two above spaces are under measure ϖ h α , β , δ , κ .
The tempered fractional derivatives of the TFJFs can also be represented by the same class of functions.
Theorem 1.
For n 0 , there holds
1 R D s μ , κ J n , l α , β , β , κ ( s ) = λ n β , μ J n , l α + μ , β μ , β μ , κ ( s ) , α R , β > μ 1 . s R D 1 μ , κ J n , r α , β , α , κ ( s ) = λ n α , μ J n , r α μ , β + μ , α μ , κ ( s ) , α > μ 1 , β R .
where
λ n β , μ = Γ ( n + β + 1 ) Γ ( n + β μ + 1 ) .
Proof. 
From Definitions 4, we have
1 R D s μ , κ J n , l α , β , β , κ ( s ) = e κ s 1 R D s μ ( 1 + s ) β P n α , β ( s ) .
Then, by Theorem 3.1 in [28], where
1 R D s μ ( 1 + s ) β P n α , β ( s ) = λ n β , μ ( 1 + s ) β μ P n α + μ , β μ ( s ) ,
the desired result can be obtained. The second equality is derived similarly. □

3. Approximation by the TFJFs

3.1. Error Estimate of Projection Operators

We introduce the ( N + 1 ) -dimensional spaces of the TFJFs as
F N , l δ , κ : = e κ s ( 1 + s ) δ P N = s p a n J 0 , l α , β , δ , κ , J 1 , l α , β , δ , κ , , J N , l α , β , δ , κ ,
and
F N , r δ , κ : = e κ s ( 1 s ) δ P N = s p a n J 0 , r α , β , δ , κ , J 1 , r α , β , δ , κ , , J N , r α , β , δ , κ .
Let α , β > 1 . Denote P N , l δ , κ : L ϖ l α , β , δ , κ 2 ( I ) F N , l δ , κ as the orthogonal projection such that
( P N , l δ , κ u u , v ) ϖ l α , β , δ , κ = 0 v F N , l δ , κ ,
and P N , r δ , κ : L ϖ r α , β , δ , κ 2 ( I ) F N , r δ , κ as the orthogonal projection such that
( P N , r δ , κ u u , v ) ϖ r α , β , δ , κ = 0 v F N , r δ , κ .
Consider P N , l δ , κ . One can express the orthogonal projection as
P N , l δ , κ u ( s ) = k = 0 N u ^ k J k , l α , β , δ , κ ( s ) ,
with
u ^ k = 1 1 u ( s ) J k , l α , β , δ , κ ( s ) ϖ l α , β , δ , κ ( s ) d s J k , l α , β , δ , κ ϖ l α , β , δ , κ 2 .
For any u L ϖ l α , β , δ , κ 2 ( I ) , we have e κ s ( 1 + s ) δ u L ω α , β 2 ( I ) . Since the Jacobi polynomials { P n α , β ( x ) } n 0 are mutually orthogonal and complete in L ω α , β 2 ( I ) , it uniquely holds that
e κ s ( 1 + s ) δ u ( s ) = n = 0 u ^ n P n α , β ( s ) ,
where
u ^ n = 1 1 e κ s ( 1 + s ) δ u ( s ) P n α , β ( s ) ω α , β ( s ) d s P n α , β ω α , β 2 = 1 1 u ( s ) J n , l α , β , β , κ ( s ) ϖ l α , β , δ , κ ( s ) d s J n , l α , β , δ , κ ϖ l α , β , δ , κ 2 .
That is, the unique representation of u is
u ( s ) = n = 0 u ^ n e κ s ( 1 + s ) δ P n α , β ( s ) = n = 0 u ^ n J n , l α , β , δ , κ ( s ) , s I a . e .
Hence, the set { J n , l α , β , δ , κ } n 0 is complete in L ϖ l α , β , δ , κ 2 ( I ) . By the definition of the projection operator and Parseval equality, one has
P N , l δ , κ u u ϖ l α , β , δ , κ 2 = n = N + 1 | u ^ n | 2 γ n α , β 0 ( N ) .
It is clear that a parallel statement for P N , r δ , κ holds. We do not repeat that.
To describe the projection error, we define
B ˇ α , β , δ , l ν ( I ) : = u L ϖ l α , β , δ , κ 2 ( I ) : 1 R D s μ , κ u L ϖ l α + μ , β μ , δ μ , κ 2 ( I ) , 0 μ ν ,
and
B ˇ α , β , δ , r ν ( I ) : = u L ϖ r α , β , δ , κ 2 ( I ) : s R D 1 μ , κ u L ϖ r α μ , β + μ , δ μ , κ 2 ( I ) , 0 μ ν .
We have the following error estimate results about the projection operators defined in (12) and (13).
Theorem 2.
Let α > 1 , β > ν 1 , 0 μ ν . For any u B ˇ α , β , β , l ν ( I ) ,
1 R D s μ , κ ( P N , l β , κ u u ) ϖ l α + μ , β μ , β μ , κ c N μ ν 1 R D s ν , κ u ϖ l α + ν , β ν , β ν , κ .
and for α > ν 1 , β > 1 , u B ˇ α , β , α , r ν ( I ) ,
s R D 1 μ , κ ( P N , r α , κ u u ) ϖ r α μ , β + μ , α μ , κ c N μ ν s R D 1 ν , κ u ϖ r α ν , β + ν , α ν , κ .
Proof. 
We first consider P N , l β , κ . From Equation (14) and the orthogonality of Lemma 1, one has
P N , l β , κ u u = k = N + 1 u ^ k J n , l α , β , β , κ .
Now from Theorem 1, one has
1 R D s μ , κ ( P N β , κ u u ) ϖ l α + μ , μ β , μ β , κ 2 = k = N + 1 | u ^ k | 2 λ k β , μ 2 γ k α + μ , β μ , = k = N + 1 | u ^ k | 2 λ k β , ν 2 γ k α + ν , β ν Γ ( k + β ν + 1 ) Γ ( k + α + μ + 1 ) Γ ( k + β μ + 1 ) Γ ( k + α + ν + 1 ) Γ ( N + β ν + 2 ) Γ ( N + α + μ + 2 ) Γ ( N + β μ + 2 ) Γ ( N + α + ν + 2 ) k = N + 1 | u ^ k | 2 λ k β , ν 2 γ k α + ν , β ν c N 2 ( μ ν ) 1 R D s ν , κ u ϖ l α + ν , ν β , ν β , κ 2 .
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.
Theorem 3.
For α , β > 1 . If u L ϖ l α , β , 0 , κ 2 ( I ) such that v ( s ) = e κ s u ( s ) B α , β m ( I ) , we have
P N , l 0 , κ u u ϖ l α , β , 0 , κ c N m s m v B α , β m .
If u L ϖ r α , β , 0 , κ 2 ( I ) such that v ( s ) = e κ s u ( s ) B α , β m ( I ) , we have
P N , r 0 , κ u u ϖ r α , β , 0 , κ c N m s m v B α , β m ,
where the non-uniformly Jacobi-weighted Sobolev space
B α , β m ( I ) : = v : s k v L ω α + k , β + k 2 ( I ) , 0 k m ,
equipped with the inner product and norm
u , v B α , β m = k = 0 m s k u , s k v ω α + k , β + k , u B α , β m = u , u B α , β m .
Proof. 
For u L ϖ l α , β , 0 , κ 2 ( I ) , it is easy to know that P N , l 0 , κ u = P N 0 , 0 v ( P N , r 0 , κ u = P N 0 , 0 v for u L ϖ r α , β , 0 , κ 2 ( I ) ) and P N , l δ , κ u u ϖ l α , β , 0 , κ = P N 0 , 0 v v ω α , β (also P N , r δ , κ u u ϖ r α , β , 0 , κ = P N 0 , 0 v v ω α , β ). Then, by utilizing the approximation result of P N 0 , 0 (Theorem 3.35 in page 118 of [30]), we can achieve the desired estimate result. □

3.2. Error Estimate of Interpolation Operators

In the sequence, we denote { ξ i } i = 0 N as the Jacobi–Gauss–Lobatto points(JGL) in I ¯ , which are also the roots of ( 1 s 2 ) P N 1 α , β ( s ) at the same time.
For a function u ( s ) C ( I ) , which satisfies e κ s ( 1 + s ) δ u ( s ) being continuous on I ¯ for some δ > 1 , we define the interpolation operator Π N , l δ , κ : C ( I ) F N , l δ , κ at given nodes { ξ i } i = 0 N as
Π N , l δ , κ u ( s ) = u ( s ) , s = ξ i , i = 0 , 1 , , N .
Similarly, for a function u ( s ) C ( I ) , which satisfies e κ s ( 1 s ) δ u ( s ) being continuous on I ¯ for some δ > 1 , we define the interpolation operator Π N , r δ , κ : C ( I ) F N , r δ , κ at the same nodes { ξ i } i = 0 N as
Π N , r δ , κ u ( s ) = u ( s ) , s = ξ i , i = 0 , 1 , , N .
If e κ s ( 1 + s ) δ u ( s ) C ( I ¯ ) , we can define the Jacobi–Gauss–Lobatto interpolation operator Π N 0 , 0 for e κ s ( 1 + s ) δ u ( s ) . Then,
Π N , l δ , κ u ( s ) = e κ s ( 1 + s ) δ Π N 0 , 0 e κ s ( 1 + s ) δ u ( s ) .
In the same way, if e κ s ( 1 s ) δ u ( s ) C ( I ¯ ) , we can define the Jacobi–Gauss–Lobatto interpolation operator Π N 0 , 0 for e κ s ( 1 s ) δ u ( s ) . Then,
Π N , r δ , κ u ( s ) = e κ s ( 1 s ) δ Π N 0 , 0 e κ s ( 1 s ) δ u ( s ) .
Theorem 4.
For α , β > 1 . If v = e κ s ( 1 + s ) δ u B α , β m ( I ) , then
Π N , l δ , κ u u ϖ l α , β , δ , κ c N m s m v B α , β m .
and if v = e κ s ( 1 s ) δ u B α , β m ( I ) , then
Π N , r δ , κ u u ϖ r α , β , δ , κ c N m s m v B α , β m .
Proof. 
Since Π N , l δ , κ u u ϖ l α , β , δ , κ = Π N 0 , 0 v v ω α , β , by utilizing the approximation result of Π N 0 , 0 (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
Π N , l δ , κ u ϖ l α , β , δ , κ Π N , l δ , κ u u ϖ l α , β , δ , κ + u ϖ l α , β , δ , κ C u ϖ l α , β , δ , κ .
and
Π N , r δ , κ u ϖ r α , β , δ , κ C u ϖ r α , β , δ , κ .
This is the stability of the interpolations Π N , l δ , κ and Π N , r δ , κ .
Now, we give an inverse inequality.
Theorem 5.
For 1 < α , β < 1 , it holds
1 R D s μ , κ ϕ ϖ l α + μ , β μ , β μ , κ C N μ ϕ ϖ l α , β , β , κ , ϕ F N , l β , κ , s R D 1 μ , κ ϕ ϖ r α μ , β + μ , α μ , κ C N μ ϕ ϖ r α , β , α , κ , ϕ F N , r α , κ .
Proof. 
For ϕ ( s ) F N , l β , κ , let ϕ ( s ) = k = 0 N ϕ ^ k J k , l α , β , β , κ ( s ) and one has
ϕ ϖ l α , β , β , κ 2 = k = 0 N | ϕ ^ k | 2 γ k α , β .
On the other hand, from Theorem 1, one has
1 R D s μ , κ ϕ = k = 0 N ϕ ^ k 1 R D s μ , κ J k , l α , β , β , κ ( s ) = k = 0 N ϕ ^ k λ k β , μ J k , l α + μ , β μ , β μ , κ ( s ) .
Then,
1 R D s μ , κ ϕ ϖ l α + μ , β μ , β μ , κ 2 = k = 0 N | ϕ ^ k | 2 λ k β , μ 2 γ k α + μ , β μ = k = 0 N λ k β , μ 2 γ k α + μ , β μ γ k α , β | ϕ ^ k | 2 γ k α , β Γ ( N + β + 1 ) Γ ( N + α + μ + 1 ) Γ ( N + β μ + 1 ) Γ ( N + α + 1 ) k = 0 N | ϕ ^ k | 2 γ k α , β C N 2 μ ϕ ϖ l α , β , β , κ 2 .
This is the first inequality. For ϕ ( s ) F N , r α , κ , let ϕ ( s ) = k = 0 N ϕ ^ k J k , r α , β , α , κ ( s ) and one has
ϕ ϖ r α , β , α , κ 2 = k = 0 N | ϕ ^ k | 2 γ k α , β .
On the other hand, from Theorem 1, one has
s R D 1 μ , κ ϕ = k = 0 N ϕ ^ k s R D 1 μ , κ J k , r α , β , α , κ ( s ) = k = 0 N ϕ ^ k λ k α , μ J k , r α μ , β + μ , α μ , κ ( s ) .
Then
s R D 1 μ , κ ϕ ϖ r α μ , β + μ , α μ , κ 2 = k = 0 N | ϕ ^ k | 2 λ k α , μ 2 γ k α μ , β + μ = k = 0 N λ k α , μ 2 γ k α μ , β + μ γ k α , β | ϕ ^ k | 2 γ k α , β Γ ( N + α + 1 ) Γ ( N + β + μ + 1 ) Γ ( N + α μ + 1 ) Γ ( N + β + 1 ) k = 0 N | ϕ ^ k | 2 γ k α , β C N 2 μ ϕ ϖ r α , β , α , κ 2 .
This is the second inequality. □
Then, we have the following estimate.
Theorem 6.
For α , β > ν 1 and 0 μ ν . If u B ˇ α , β , β , l ν ( I ) , then
1 R D s μ , κ ( Π N , l β , κ u u ) ϖ l α + μ , β μ , β μ , κ c N μ ν 1 R D s ν , κ u ϖ l α + ν , β ν , β ν , κ ,
and if u B ˇ α , β , α , r ν ( I ) , then
s R D 1 μ , κ ( Π N , r α , κ u u ) ϖ r α μ , β + μ , α μ , κ c N μ ν s R D 1 ν , κ u ϖ r α ν , β + ν , α ν , κ .
Proof. 
Since Π N , l β , κ ( P N , l β , κ u ) = P N , l β , κ u , then, from Theorems 2, 5 and the stability (17), one has
1 R D s μ , κ ( Π N , l β , κ u u ) ϖ l α + μ , β μ , β μ , κ 1 R D s μ , κ ( Π N , l β , κ ( P N , l β , κ u u ) ) ϖ l α + μ , β μ , β μ , κ + 1 R D s μ , κ ( P N , l β , κ u u ) ϖ l α + μ , β μ , β μ , κ C 1 N μ Π N , l β , κ ( P N , l β , κ u u ) ϖ l α , β , β , κ + C 2 N μ ν 1 R D s ν , κ u ϖ l α + ν , β ν , β ν , κ C 3 N μ P N , l β , κ u u ϖ l α , β , β , κ + C 2 N μ ν 1 R D s ν , κ u ϖ l α + ν , β ν , β ν , κ c N μ ν 1 R D s ν , κ u ϖ l α + ν , β ν , β ν , κ .
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
s R D 1 μ , κ ( Π N , r α , κ u u ) ϖ r α μ , β + μ , α μ , κ s R D 1 μ , κ ( Π N , r α , κ ( P N , r α , κ u u ) ) ϖ r α μ , β + μ , α μ , κ + s R D 1 μ , κ ( P N , r α , κ u u ) ϖ r α μ , β + μ , α μ , κ C 1 N μ Π N , r α , κ ( P N , r α , κ u u ) ϖ r α , β , α , κ + C 2 N μ ν s R D 1 ν , κ u ϖ r α ν , β + ν , α ν , κ C 3 N μ P N , r α , κ u u ϖ r α , β , β , κ + C 2 N μ ν s R D 1 ν , κ u ϖ r α ν , β + ν , α ν , κ c N μ ν s R D 1 ν , κ u ϖ r α ν , β + ν , α ν , κ .
This ends the proof. □

4. Differentiation Matrix of the Collocation Method

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

4.1. Differentiation Matrix for the Left Tempered Derivative

For the interpolation operator Π N , l δ , κ , the tempered fractional Lagrange interpolation basis functions are defined by
F i , l δ , κ ( s ) : = e κ ( s ξ i ) 1 + s 1 + ξ i δ L i ( s ) , i = 0 , 1 , 2 , , N ,
which satisfies F i , l δ , κ ( ξ k ) = δ i k . 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
u ( s ) u N ( s ) = i = 0 N u ( ξ i ) F i , l δ , κ ( s ) .
Naturally, the tempered Caputo fractional derivative of u ( s ) can be approximated by its interpolation
1 C D s μ , κ u ( s ) i = 0 N u ( ξ i ) 1 C D s μ , κ F i , l δ , κ ( s ) .
Now, we compute the differentiation matrix of the left tempered Caputo derivative (DMLTCD), which is denoted as
D s , l μ , κ ( N + 1 ) × ( N + 1 ) : = ( 1 C D s μ , κ F i , l δ , κ ( ξ k ) ) k , i = 0 N .
To evaluate 1 C D s μ , κ F i , l δ , κ ( 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])
L i ( s ) = j = 0 N l i j P j α , β ( s ) ,
where
l i j = P j α , β ( ξ i ) ω i γ j α , β , j = 0 , 1 , , N 1 ; l i N = P N α , β ( ξ i ) ω i 2 + α + β + 1 N γ N α , β ,
and { ω i } i = 0 N are the weights corresponding with the nodes { ξ i } i = 0 N in the Jacobi–Gauss–Lobatto quadrature. Then, we obtain that
1 C D s μ , κ F i , l δ , κ ( s ) = e κ ξ i ( 1 + ξ i ) δ j = 0 N l i j 1 C D s μ , κ J i , l α , β , δ , κ ( s ) .
Denote
S n , α , β , l δ , κ , μ ( s ) : = 1 I s μ , κ J n , l α , β , δ , κ ( s ) = 1 Γ ( μ ) 1 s ( s t ) μ 1 e κ ( s t ) J n , l α , β , δ , κ ( t ) d t .
In the sequence, we demonstrate how to compute S n , α , β , l δ , κ , μ ( s ) .
Theorem 7.
Let α , β , δ > 1 and μ > 0 , s I . Then S n , l : = S n , α , β , l δ , κ , μ ( s ) and S ^ n , l : = μ S n , α , β , l δ , κ , μ + 1 ( s ) satisfy
S n + 1 , l = ( A n α , β s B n α , β ) S n , l C n α , β S n 1 , l A n α , β S ^ n , l , n 1
and
S ^ n + 1 , l = A ˜ n S ^ n , l B ˜ n S ^ n 1 , l + ( 1 + s ) a ˜ n S n 1 , l + b ˜ n S n , l + c ˜ n S n + 1 , l , n 1
with the starting terms
S 0 , l = Γ ( δ + 1 ) Γ ( δ + μ + 1 ) e κ s ( 1 + s ) δ + μ , S ^ 0 , l = μ δ + μ + 1 ( 1 + s ) S 0 , l , S 1 , l = ( α + β + 2 ) ( δ + 1 ) 2 ( δ + μ + 1 ) ( 1 + s ) ( β + 1 ) S 0 , l , S ^ 1 , l = ( α + β + 2 ) ( δ + 1 ) 2 ( δ + μ + 2 ) ( 1 + s ) ( β + 1 ) S ^ 0 , l ,
respectively, where by denoting f n α , β = 1 + ( δ + μ + 1 ) A n α , β C ^ n α , β ,
A ˜ n = ( A n α , β + B n α , β ) ( δ + μ + 1 ) A n α , β B ^ n α , β f n α , β , B ˜ n = C n α , β + ( δ + μ + 1 ) A n α , β A ^ n α , β f n α , β , a ˜ n = μ A n α , β A ^ n α , β f n α , β , b ˜ n = μ A n α , β B ^ n α , β f n α , β , c ˜ n = μ A n α , β C ^ n α , β f n α , β ,
in which A n α , β , B n α , β , C n α , β and A ^ n α , β , B ^ n α , β , A ^ n α , β are given in Lemma 1 and (6), respectively.
Proof. 
By making use of the relation (7), we have
S n + 1 , l = 1 I s μ , κ ( A n α , β s B n α , β ) J n , l α , β , δ , κ ( s ) C n α , β J n 1 , l α , β , δ , κ ( s ) = A n α , β 1 I s μ , κ s J n , l α , β , δ , κ ( s ) B n α , β S n , l C n α , β S n 1 , l .
Thus,
1 I s μ , κ s J n , l α , β , δ , κ ( s ) = e κ s Γ ( μ ) 1 s ( s t ) μ 1 t ( 1 + t ) δ P n α , β ( t ) d t = e κ s Γ ( μ ) 1 s ( s t ) μ ( 1 + t ) δ P n α , β ( t ) d t + s 1 s ( s t ) μ 1 ( 1 + t ) δ P n α , β ( t ) d t = S ^ n , l + s S n , l .
Then, we obtain the first equality. With integration by parts, and from (5), we have
1 s ( s t ) μ ( 1 + t ) δ + 1 P n α , β ( t ) d t = 1 s ( s t ) μ ( 1 + t ) δ + 1 A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t = 1 s ( s t ) μ ( 1 + t ) δ + 1 A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t = ( δ + μ + 1 ) 1 s ( s t ) μ ( 1 + t ) δ A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t + μ ( 1 + s ) 1 s ( s t ) μ 1 ( 1 + t ) δ A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t .
Hence, we have
S ^ n + 1 , l = e κ s Γ ( μ ) 1 s ( s t ) μ ( 1 + t ) δ P n + 1 α , β ( t ) d t = e κ s Γ ( μ ) 1 s ( s t ) μ ( 1 + t ) δ ( A n α , β t B n α , β ) P n α , β ( t ) C n α , β P n 1 α , β ( t ) d t = e κ s Γ ( μ ) A n α , β 1 s ( s t ) μ t ( 1 + t ) δ P n α , β ( t ) d t B n α , β S ^ n , l C n α , β S ^ n 1 , l = e κ s Γ ( μ ) A n α , β 1 s ( s t ) μ ( 1 + t ) δ + 1 P n α , β ( t ) d t ( A n α , β + B n α , β ) S ^ n , l C n α , β S ^ n 1 , l = ( δ + μ + 1 ) A n α , β C ^ n α , β S ^ n + 1 , l + ( δ + μ + 1 ) A n α , β B ^ n α , β ( A n α , β + B n α , β ) S ^ n , l + ( δ + μ + 1 ) A n α , β A ^ n α , β C n α , β S ^ n 1 , l + μ ( 1 + s ) A n α , β A ^ n α , β S n 1 , l + B ^ n α , β S n , l + C ^ n α , β S n + 1 , l .
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
1 C D s μ , κ J i , l α , β , δ , κ ( s ) = 1 I s n μ , κ ( e κ s J i , l α , β , δ , κ ( s ) ) ( n ) ,
we can evaluate quickly and stably by Theorem 7.

4.2. Differentiation Matrix for the Right Tempered Derivative

For the interpolation operator Π N , r δ , κ , the tempered fractional Lagrange interpolation basis functions are defined by
F i , r δ , κ ( s ) : = e κ ( s ξ i ) 1 s 1 ξ i δ L i ( s ) , i = 0 , 1 , 2 , , N ,
which satisfies F i , r δ , κ ( ξ k ) = δ i k . 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
u ( s ) u N ( s ) = i = 0 N u ( ξ i ) F i , r δ , κ ( s ) .
Naturally, the right tempered Caputo fractional derivative of u ( s ) can be approximated by its interpolation
s C D 1 μ , κ u ( s ) i = 0 N u ( ξ i ) s C D 1 μ , κ F i , r δ , κ ( s ) .
Now, we compute the differentiation matrix of the right tempered Caputo derivative(DMRTCD), which is denoted as
D s , r μ , κ ( N + 1 ) × ( N + 1 ) : = ( s C D 1 μ , κ F i , r δ , κ ( ξ k ) ) k , i = 0 N .
Similar to the above subsection, we use (20) to obtain
s C D 1 μ , κ F i , r δ , κ ( s ) = e κ ξ i ( 1 ξ i ) δ j = 0 N l i j s C D 1 μ , κ J i , r α , β , δ , κ ( s ) .
Now, denote
S n , α , β , r δ , κ , μ ( s ) : = s I 1 μ , κ J n , r α , β , δ , κ ( s ) = 1 Γ ( μ ) s 1 ( t s ) μ 1 e κ ( t s ) J n , r α , β , δ , κ ( t ) d t .
Theorem 8.
Let α , β , δ > 1 and μ > 0 , s I . Then S n , r : = S n , α , β , r δ , κ , μ ( s ) and S ^ n , r : = μ S n , α , β , r δ , κ , μ + 1 ( s ) satisfy
S n + 1 , r = ( A n α , β s B n α , β ) S n , r C n α , β S n 1 , r + A n α , β S ^ n , r , n 1
and
S ^ n + 1 , r = C ˜ n S ^ n , r B ˜ n S ^ n 1 , r + ( 1 s ) a ˜ n S n 1 , r + b ˜ n S n , r + c ˜ n S n + 1 , r , n 1
with the starting terms
S 0 , r = Γ ( δ + 1 ) Γ ( δ + μ + 1 ) e κ s ( 1 s ) δ + μ , S ^ 0 , r = μ δ + μ + 1 ( 1 s ) S 0 , r , S 1 , r = ( α + β + 2 ) ( δ + 1 ) 2 ( δ + μ + 1 ) ( 1 s ) + ( α + 1 ) S 0 , r , S ^ 1 , r = ( α + β + 2 ) ( δ + 1 ) 2 ( δ + μ + 2 ) ( 1 s ) + ( α + 1 ) S ^ 0 , r ,
respectively, where
C ˜ n = A n α , β B n α , β ( δ + μ + 1 ) A n α , β B ^ n f n α , β ,
and f n α , β , B ˜ n and a ˜ n , b ˜ n , c ˜ n are given in Theorem 7.
Proof. 
By making use of the relation (7), we have
S n + 1 , r = s I 1 μ , κ ( A n α , β s B n α , β ) J n , r α , β , δ , κ ( s ) C n α , β J n 1 , r α , β , δ , κ ( s ) = A n α , β s I 1 μ , κ s J n , r α , β , δ , κ ( s ) B n α , β S n , r C n α , β S n 1 , r .
Thus,
s I 1 μ , κ s J n , r α , β , δ , κ ( s ) = e κ s Γ ( μ ) s 1 ( t s ) μ 1 t ( 1 t ) δ P n α , β ( t ) d t = e κ s Γ ( μ ) s 1 ( t s ) μ ( 1 t ) δ P n α , β ( t ) d t + s s 1 ( t s ) μ 1 ( 1 t ) δ P n α , β ( t ) d t = S ^ n , r + s S n , r .
Then, we obtain the first equality. With integrating by parts and from (5), we have
s 1 ( t s ) μ ( 1 t ) δ + 1 P n α , β ( t ) d t = s 1 ( t s ) μ ( 1 t ) δ + 1 A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t = s 1 ( t s ) μ ( 1 t ) δ + 1 A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t = ( δ + μ + 1 ) s 1 ( t s ) μ ( 1 t ) δ A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t μ ( 1 s ) s 1 ( t s ) μ 1 ( 1 t ) δ A ^ n α , β P n 1 α , β ( t ) + B ^ n α , β P n α , β ( t ) + C ^ n α , β P n + 1 α , β ( t ) d t .
Hence, we have
S ^ n + 1 , r = e κ s Γ ( μ ) s 1 ( t s ) μ ( 1 t ) δ P n + 1 α , β ( t ) d t = e κ s Γ ( μ ) s 1 ( t s ) μ ( 1 t ) δ ( A n α , β t B n α , β ) P n α , β ( t ) C n α , β P n 1 α , β ( t ) d t = e κ s Γ ( μ ) A n α , β s 1 ( t s ) μ t ( 1 t ) δ P n α , β ( t ) d t B n α , β S ^ n , r C n α , β S ^ n 1 , r = e κ s Γ ( μ ) A n α , β s 1 ( t s ) μ ( 1 t ) δ + 1 P n α , β ( t ) d t + ( A n α , β B n α , β ) S ^ n , r C n α , β S ^ n 1 , r = ( δ + μ + 1 ) A n α , β C ^ n α , β S ^ n + 1 , r + ( δ + μ + 1 ) A n α , β B ^ n α , β + A n α , β B n α , β S ^ n , r + ( δ + μ + 1 ) A n α , β A ^ n α , β C n α , β S ^ n 1 , r + μ ( 1 s ) A n α , β A ^ n α , β S n 1 , r + B ^ n α , β S n , r + C ^ n α , β S n + 1 , r .
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 D s , l μ , κ can be obtained from D s , r μ , κ by D s , l μ , κ ( i , j ) = D s , r μ , κ ( n + 1 i , n + 1 j ) 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]):
Lemma 2.
For a continuous function u ( x ) on [ 1 , 1 ] ,
lim s 1 + 1 I s μ ( 1 + s ) δ u ( s ) = 0 , i f μ > δ , u ( 1 ) Γ ( δ + 1 ) , i f μ = δ , , i f μ < δ .
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.

5. 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:
0 C D x μ , κ u ( x ) = f ( x , u ) , 0 < x T , u ( 0 ) = u 0 .
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 u N F N , l δ , κ such that
0 C D x μ , κ ˜ u N ( x ) = f ( x , u N ( x ) ) , x = x j α , β = T ( 1 + ξ j ) 2 , j = 1 , 2 , , N ,
and
u N ( 0 ) = u 0 .
The above two equations lead to the following linear system
2 T μ D ¯ s , l μ , κ ˜ u = f ( u ) ,
where D ¯ s , l μ , κ ˜ is the submatrix of D s , l μ , κ ˜ , which deletes the first row and column of D s , l μ , κ ˜ , u = u N ( x ) , f ( u ) = f ( x , u ) + u a d 0 , x = [ x 1 α , β , x 2 α , β , , x N α , β ] T , d 0 the first column of D s , l μ , κ ˜ .
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
E 0 = max i = 1 , 2 , , N { | u ( x i ) u N ( x i ) | } ,
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
g ( x ) = u 2 + 0 C D x μ , κ u ( x ) .
The following two cases of the exact solutions are considered [33]:
C1. 
u ( x ) = e κ x ( x 2 x ) with T = 5 .
C2. 
u ( x ) = e κ x x 8 3 x 4 + μ / 2 + 9 4 x μ with T = 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:
u n + 1 = u n + D + 2 d i a g ( u n ) 1 ( D u n + u n . u n g ) , n = 0 , 1 , ,
where g = [ g ( x 1 α , β ) , g ( x 2 α , β ) , , g ( x N α , β ) ] T and D = ( 2 T ) μ D ¯ s , l μ , κ ˜ .
Since u ( x ) F 2 0 , κ , 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.

6. Applications to the Fractional Helmholtz and Burgers Equations

6.1. Fractional Helmholtz Equation

Let 1 < μ < 2 . In this subsection, we apply the spectral collocation method to the following fractional Helmholtz equation:
λ 2 u ( x ) a C D x μ , κ u ( x ) = f ( x ) , a < x < b , u ( a ) = u a , u ( b ) = u b .
The spectral collocation method based on the TFJFs for (32) is to find u N F N δ , κ such that
λ 2 u N ( x ) a C D x μ , κ u N ( x ) = f ( x ) , x = x j α , β , j = 1 , 2 , , N 1 ,
and
u N ( a ) = u a , u N ( b ) = u b .
The above two equations lead to the following linear system:
λ 2 I 2 b a μ D ˜ s , l μ , κ ˜ u = f ,
where I is the unitary matrix, D ˜ s , l μ , κ ˜ 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 μ , κ ˜ .
To measure the accuracy of the method, we define the errors by
E 0 ( N ) = max i = 1 , , N 1 { | u ( x i ) u N ( x i ) | } ,
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
O r d ( N 1 , N 2 ) = log ( E 0 ( N 1 ) ) log ( E 0 ( N 2 ) ) log ( N 1 ) log ( N 2 ) .
Example 2.
Consider (32) with [ a , b ] = [ 0 , 2 ] . The source term f ( x ) is chosen such that the problem satisfies one of the following two cases:
C1. 
Smooth solution: u ( x ) = e κ x sin π x .
C2. 
Solution with low regularity: u ( x ) = e κ x 2 σ + 1 x σ x 2 σ + 1 .
The homogeneous boundary conditions are implemented for two cases. The term a C D x μ , κ u ( x ) in C1 is approximated by
0 C D x μ , κ e κ x sin π x e κ x k = 1 L ( 1 ) k π 2 k + 1 Γ ( 2 k + 2 μ ) x 2 k + 1 μ ,
with L = 50 to compute the right-hand function (RHF) f ( x ) .
For the smooth solution C1, we solve Equation (32) by δ = 0 . We list the errors E 0 in Table 2, Table 3 and Table 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 O r d in Table 5, Table 6, Table 7 and Table 8. A possible guess of the O r d for the case C2 with δ = 0 is O ( N 2 ( σ μ + 1 ) ) from the results in Table 5 and Table 7. When δ = σ 1 , the super convergence can be observed; see the results in Table 6 and Table 8.

6.2. Fractional Burgers Equation

In the current subsection, we employ the spectral collocation method to solve the fractional Burgers equation (FBE)
t u ( x , t ) + u ( x , t ) x u ( x , t ) = ϵ a C D x μ , κ u ( x , t ) ,
subject to homogeneous Dirichlet boundary conditions and initial condition u ( x , 0 ) = u 0 ( x ) , where ϵ > 0 and 1 < μ < 2 , ( x , t ) ( a , b ) × ( 0 , 1 ] .
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
( I ϵ τ D ˜ ) u n + 1 = ( I + ϵ τ D ˜ ) u n 1 2 τ ( d i a g ( u n ) D ) u n , n 1 , u 1 = ( I + ϵ τ D ˜ ) u 0 τ ( d i a g ( u 0 ) D ) u 0 , u 0 = u 0 ( x ) ,
where D is the first-order differentiation matrix and D ˜ = ( 2 b a ) μ D ˜ s , l μ , κ .
Example 3.
Consider (36) with the initial profile as one of the following:
C1. 
u 0 ( x ) = sin ( π x ) , x [ 1 , 1 ] .
C2. 
u 0 ( x ) = e 4 x 2 , x [ 6 , 6 ] .
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.

7. 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.

Funding

This research was funded by the National Natural Science Foundation of China under Grant No. 11661048.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The author acknowledge ’Overleaf’ for easy typesetting.

Conflicts of Interest

The author declares no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
TFJFstempered fractional Jacobi functions
FCfractional calculus
CTRWcontinuous time random walk
TFCtempered fractional calculus
TFDtempered fractional differential
DMLTCDdifferentiation matrix of left tempered Caputo derivative
DMRTCDdifferentiation matrix of right tempered Caputo derivative

References

  1. Baeumera, B.; Meerschaert, M.M. Tempered stable Levy motion and transient super-diffusion. J. Comput. Appl. Math. 2010, 233, 2438–2448. [Google Scholar] [CrossRef] [Green Version]
  2. Sabzikar, F.; Meerschaerta, M.M.; Chen, J.H. Tempered fractional calculus. J. Comput. Phys. 2015, 293, 14–28. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Sibatov, R.T.; Sun, H.G. Tempered fractional equations for quantum transport in mesoscopic one-dimensional systems with fractal disorder. Fractal Fract. 2019, 3, 47. [Google Scholar] [CrossRef] [Green Version]
  4. Wang, H.H.; Li, C. Fast difference scheme for a tempered fractional Burgers equation in porous media. Appl. Math. Lett. 2022, 132, 108143. [Google Scholar] [CrossRef]
  5. Meerschaert, M.M.; Zhang, Y.; Baeumer, B. Tempered anomalous diffusion in heterogeneous systems. Geophys. Res. Lett. 2008, 35, L17403. [Google Scholar] [CrossRef]
  6. Meerschaert, M.M.; Sabzikar, F.; Phanikumar, M.; Zeleke, A. Tempered fractional time series model for turbulence in geophysical flows. J. Stat. Mech. Theory Exp. 2014, 14, 1742–5468. [Google Scholar] [CrossRef] [Green Version]
  7. Li, C.; Deng, W.H. High order schemes for the tempered fractional diffusion equations. Adv. Comput. Math. 2016, 42, 543–572. [Google Scholar] [CrossRef] [Green Version]
  8. Deng, W.H.; Zhang, Z.J. Variational formulation and efficient implementation for solving the tempered fractional problems. Numer. Methods Part. Diff. Eqns. 2018, 34, 1224–1257. [Google Scholar] [CrossRef] [Green Version]
  9. Chen, M.H.; Deng, W.H. A second-order accurate numerical method for the space-time tempered fractional diffusion-wave equation. Appl. Math. Lett. 2017, 68, 87–93. [Google Scholar] [CrossRef]
  10. Ding, H.F.; Li, C.P. A high-order algorithm for time-Caputo-tempered partial differential equation with Riesz derivatives in two spatial dimensions. J. Sci. Comput. 2019, 80, 81–109. [Google Scholar] [CrossRef]
  11. Guo, L.; Zeng, F.H.; Turner, I.; Burrage, K.; Karniadakis, G.E. Efficient multistep methods for tempered fractional calculus: Algorithms and simulations. SIAM J. Sci. Comput. 2019, 41, A2510–A2535. [Google Scholar] [CrossRef] [Green Version]
  12. Cao, J.L.; Xiao, A.G.; Bu, W.P. Finite difference/finite element method for tempered time fractional advection-dispersion equation with fast evaluation of Caputo derivative. J. Sci. Comput. 2020, 83, 48. [Google Scholar] [CrossRef]
  13. Sun, X.R.; Li, C.; Zhao, F.Q. Local discontinuous Galerkin methods for the time tempered fractional diffusion equation. Appl. Math. Comput. 2020, 365, 124725. [Google Scholar] [CrossRef] [Green Version]
  14. Safari, Z.; Loghmani, G.B.; Ahmadinia, M. Convergence analysis of a LDG method for time-space tempered fractional diffusion equations with weakly singular solutions. J. Sci. Comput. 2022, 91, 68. [Google Scholar] [CrossRef]
  15. Bu, L.L.; Oosterlee, C.W. On a multigrid method for tempered fractional diffusion equations. Fractal Fract. 2021, 5, 145. [Google Scholar] [CrossRef]
  16. Çelik, C.; Duman, M. Finite element method for a symmetric tempered fractional diffusion equation. Appl. Numer. Math. 2017, 120, 270–286. [Google Scholar] [CrossRef]
  17. Zayernouri, M.; Ainsworth, M.; Karniadakis, G.E. Tempered fractional Sturm-Liouville eigenproblems. SIAM J. Sci. Comput. 2015, 37, A1777–A1800. [Google Scholar] [CrossRef] [Green Version]
  18. Hanert, E.; Piret, C. A Chebyshev pseudospectral method to solve the space-time tempered fractional diffusion equation. SIAM J. Sci. Comput. 2015, 36, A1797–A1812. [Google Scholar] [CrossRef]
  19. Chen, S.; Shen, J.; Wang, L.L. Laguerre functions and their applications to tempered fractional differential equations on infinite intervals. J. Sci. Comput. 2018, 74, 1286–1313. [Google Scholar] [CrossRef]
  20. Luo, W.H.; Gu, X.M.; Yang, L.; Meng, J. A Lagrange-quadratic spline optimal collocation method for the time tempered fractional diffusion equation. Math. Comput. Simulat. 2021, 182, 1–24. [Google Scholar] [CrossRef]
  21. Moghaddam, B.P.; Machado, J.A.T.; Babaei, A. A computationally efficient method for tempered fractional differential equations with application. Comp. Appl. Math. 2018, 37, 3657–3671. [Google Scholar] [CrossRef]
  22. Liemert, A.; Kienle, A. Computational solutions of the tempered fractional wave-diffusion equation. Fract. Calc. Appl. Aanl. 2017, 20, 139–158. [Google Scholar] [CrossRef]
  23. Li, C.; Deng, W.H.; Zhao, L.J. Well-posedness and numerical algorithm for the tempered fractional ordinary differential equations. Disc. Cont. Dyn. Sys. B 2019, 24, 1989–2015. [Google Scholar]
  24. Deng, W.H.; Li, B.Y.; Tian, W.Y.; Zhang, P.W. Boundary problems for the fractional and tempered fractional operators. Multiscale Model. Simul. 2018, 16, 125–149. [Google Scholar] [CrossRef] [Green Version]
  25. Fahad, H.M.; Fernandez, A.; Rehman, M.U.; Siddiqi, M. Tempered and Hadamard-Type fractional calculus with respect to functions. Mediterr. J. Math. 2021, 18, 143. [Google Scholar] [CrossRef]
  26. Chen, M.H.; Deng, W.H. High order algorithms for the fractional substantial diffusion equation with truncated Lévy flights. SIAM J. Sci. Comput. 2015, 37, A890–A917. [Google Scholar] [CrossRef] [Green Version]
  27. Huang, C.; Zhang, Z.M.; Song, Q.S. Spectral method for substantial fractional differential equations. J. Sci. Comput. 2018, 74, 1554–1574. [Google Scholar] [CrossRef]
  28. Chen, S.; Shen, J.; Wang, L.L. Generalized Jacobi functions and their applications to fractional differential equations. Math. Comput. 2016, 85, 1603–1638. [Google Scholar] [CrossRef] [Green Version]
  29. Wu, X.C.; Deng, W.H.; Barkai, E. Tempered fractional Feynman-Kac equation: Theory and examples. Phys. Rev. E 2016, 93, 032151. [Google Scholar] [CrossRef] [Green Version]
  30. Shen, J.; Tang, T.; Wang, L.L. Spectral Methods: Algorithms, Analysis and Applications; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  31. Szegö, G. Orthogonal Polynomials, 4th ed.; American Mathematical Society: Providence, RI, USA, 1975. [Google Scholar]
  32. Liu, W.J.; Wang, L.L.; Wu, B.Y. Optimal error estimates for Legendre expansions of singular functions with fractional derivatives of bounded variation. Adv. Comput. Math. 2021, 47, 79. [Google Scholar] [CrossRef]
  33. Deng, J.W.; Zhao, L.J.; Wu, Y.J. Fast predictor-corrector approach for the tempered fractional differential equations. Numer. Algor. 2017, 74, 717–754. [Google Scholar] [CrossRef]
  34. Zhao, T.G.; Mao, Z.P.; Karniadakis, G.E. Multi-domain spectral collocation method for variable-order nonlinear fractional differential equations. Comput. Methods Appl. Mech. Eng. 2019, 348, 377–395. [Google Scholar] [CrossRef] [Green Version]
  35. Zhou, Y.; Zhang, Y. Noether symmetries for fractional generalized Birkhoffian systems in terms of classical and combined Caputo derivatives. Acta Mech. 2020, 231, 3017–3029. [Google Scholar] [CrossRef]
  36. Youssri, Y.H.; Atta, A.G. Spectral collocation approach via normalized shifted Jacobi polynomials for the nonlinear Lane-Emden equation with fractal-fractional derivative. Fractal Fract. 2023, 7, 133. [Google Scholar] [CrossRef]
Figure 1. Condition number of DMLTCD with 0 < μ < 1 , α = β = 0 , δ = 0 , κ = 1 .
Figure 1. Condition number of DMLTCD with 0 < μ < 1 , α = β = 0 , δ = 0 , κ = 1 .
Fractalfract 07 00277 g001
Figure 2. Condition number of DMLTCD with 1 < μ < 2 , α = β = 0 , δ = 0 , κ = 1 .
Figure 2. Condition number of DMLTCD with 1 < μ < 2 , α = β = 0 , δ = 0 , κ = 1 .
Fractalfract 07 00277 g002
Figure 3. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with μ = 1.18 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Figure 3. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with μ = 1.18 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g003
Figure 4. Numerical solution u N ( x , t ) for Example 3 of u 0 ( x ) = sin ( π x ) with μ = 1.4 , ϵ = 1 , N = 200 , α = β = 0 , τ = 10 3 .
Figure 4. Numerical solution u N ( x , t ) for Example 3 of u 0 ( x ) = sin ( π x ) with μ = 1.4 , ϵ = 1 , N = 200 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g004
Figure 5. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with N = 600 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Figure 5. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with N = 600 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g005
Figure 6. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with ϵ = 1 , N = 600 , α = β = 0 , τ = 10 3 .
Figure 6. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = sin ( π x ) with ϵ = 1 , N = 600 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g006
Figure 7. Numerical solution u N ( x , t ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 300 , μ = 1.5 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Figure 7. Numerical solution u N ( x , t ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 300 , μ = 1.5 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g007
Figure 8. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 400 , κ = 1 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Figure 8. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 400 , κ = 1 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g008
Figure 9. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 400 , μ = 1.55 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Figure 9. Numerical solution u N ( x , 1 ) for Example 3 of u 0 ( x ) = e 4 x 2 with N = 400 , μ = 1.55 , ϵ = 1 , α = β = 0 , τ = 10 3 .
Fractalfract 07 00277 g009
Table 1. The numerical errors E 0 of Example 1 with C2 solution ( α = 0 , β = 0 , κ = 1 ).
Table 1. The numerical errors E 0 of Example 1 with C2 solution ( α = 0 , β = 0 , κ = 1 ).
N ( 1 / h ) μ = 0.2  [33] μ = 0.2 μ = 0.5  [33] μ = 0.5 μ = 0.9  [33] μ = 0.9
105.811 × 10 1 1.0898 × 10 7 2.441 × 10 2 8.4114 × 10 7 7.771 × 10 3 1.3961 × 10 6
202.001 × 10 1 5.1711 × 10 10 6.381 × 10 3 2.8558 × 10 9 1.441 × 10 3 3.3180 × 10 9
406.501 × 10 2 3.8736 × 10 13 2.821 × 10 3 7.4061 × 10 12 6.001 × 10 4 7.5808 × 10 12
808.351 × 10 3 1.2290 × 10 13 3.841 × 10 4 3.6507 × 10 13 1.521 × 10 4 1.1427 × 10 13
1604.721 × 10 4 1.8164 × 10 11 2.231 × 10 5 1.5149 × 10 12 1.541 × 10 5 1.7109 × 10 13
Table 2. Errors E 0 of Example 2 with solution C1 ( α = β = 0 , δ = 0 , κ = 1 , λ = 0 ).
Table 2. Errors E 0 of Example 2 with solution C1 ( α = β = 0 , δ = 0 , κ = 1 , λ = 0 ).
N μ = 1.1 μ = 1.3 μ = 1.5 μ = 1.7 μ = 1.9 μ = 1.99
41.781 × 10 1 1.191 × 10 1 9.844 × 10 2 1.091 × 10 1 9.028 × 10 2 4.519 × 10 2
85.110 × 10 4 3.488 × 10 4 2.818 × 10 4 2.931 × 10 4 2.610 × 10 4 4.833 × 10 5
122.173 × 10 7 1.373 × 10 7 1.016 × 10 7 1.448 × 10 7 1.699 × 10 7 3.235 × 10 8
162.618 × 10 11 1.550 × 10 11 1.160 × 10 11 1.732 × 10 11 2.501 × 10 11 5.182 × 10 12
201.200 × 10 14 1.299 × 10 14 1.241 × 10 14 1.474 × 10 14 1.302 × 10 14 1.258 × 10 14
Table 3. Errors E 0 of Example 2 with solution C1 ( μ = 1.4 , δ = 0 , κ = 5 , λ = 2 ).
Table 3. Errors E 0 of Example 2 with solution C1 ( μ = 1.4 , δ = 0 , κ = 5 , λ = 2 ).
N α = β = 0 α = β = 0.5 α = β = 0.5 ( α , β ) = ( 0.3 , 0.8 ) α = β = 1
42.282 × 10 2 2.901 × 10 2 1.322 × 10 2 2.338 × 10 2 1.605 × 10 2
89.931 × 10 5 5.828 × 10 5 2.701 × 10 5 3.397 × 10 4 1.821 × 10 4
123.315 × 10 8 2.504 × 10 8 1.210 × 10 8 1.889 × 10 7 1.075 × 10 7
163.008 × 10 12 2.379 × 10 12 1.493 × 10 12 2.353 × 10 11 1.419 × 10 11
201.174 × 10 14 2.498 × 10 15 1.443 × 10 15 1.282 × 10 14 2.776 × 10 14
Table 4. Errors E 0 of Example 2 with solution C1 ( μ = 1.6 , α = β = 0 , δ = 0 , λ = 1000 ).
Table 4. Errors E 0 of Example 2 with solution C1 ( μ = 1.6 , α = β = 0 , δ = 0 , λ = 1000 ).
N κ = 1 κ = 2 κ = 3 κ = 5 κ = 8 κ = 10
42.131 × 10 6 1.509 × 10 6 1.068 × 10 6 5.353 × 10 7 1.900 × 10 7 9.521 × 10 8
83.355 × 10 8 3.035 × 10 8 2.745 × 10 8 2.247 × 10 8 1.663 × 10 8 1.361 × 10 8
123.806 × 10 11 3.632 × 10 11 3.467 × 10 11 3.158 × 10 11 2.745 × 10 11 2.500 × 10 11
169.076 × 10 15 8.868 × 10 15 8.618 × 10 15 8.188 × 10 15 7.522 × 10 15 7.119 × 10 15
202.220 × 10 16 1.110 × 10 16 1.110 × 10 16 5.551 × 10 17 2.082 × 10 17 2.776 × 10 17
Table 5. Errors E 0 and O r d of Example 2 with solution C2 ( σ = 1.4 , α = β = 0 , δ = 0 , κ = 1 , λ = 1 ).
Table 5. Errors E 0 and O r d of Example 2 with solution C2 ( σ = 1.4 , α = β = 0 , δ = 0 , κ = 1 , λ = 1 ).
μ = 1.1 μ = 1.3 μ = 1.5 μ = 1.9
N E 0 Ord E 0 Ord E 0 Ord E 0 Ord
83.744 × 10 2 -4.300 × 10 2 -4.544 × 10 2 -2.319 × 10 2 -
166.835 × 10 3 2.458.416 × 10 3 2.351.331 × 10 2 1.771.350 × 10 2 0.78
321.072 × 10 3 2.671.877 × 10 3 2.163.968 × 10 3 1.757.138 × 10 3 0.92
641.607 × 10 4 2.744.121 × 10 4 2.191.150 × 10 3 1.793.632 × 10 3 0.97
1282.482 × 10 5 2.698.961 × 10 5 2.203.327 × 10 4 1.791.829 × 10 3 0.99
2563.946 × 10 6 2.651.950 × 10 5 2.209.585 × 10 5 1.809.166 × 10 4 1.00
5126.369 × 10 7 2.634.245 × 10 6 2.202.757 × 10 5 1.804.588 × 10 4 1.00
Table 6. Errors E 0 and O r d of Example 2 with solution C2 ( σ = 1.4 , α = β = 0 , δ = 0.4 , κ = 1 , λ = 1 ).
Table 6. Errors E 0 and O r d of Example 2 with solution C2 ( σ = 1.4 , α = β = 0 , δ = 0.4 , κ = 1 , λ = 1 ).
μ = 1.1 μ = 1.3 μ = 1.5 μ = 1.9
N E 0 Ord E 0 Ord E 0 Ord E 0 Ord
84.526 × 10 5 -1.453 × 10 5 -9.694 × 10 6 -5.743 × 10 6 -
165.361 × 10 7 6.401.066 × 10 7 7.095.753 × 10 8 7.403.732 × 10 8 7.27
321.061 × 10 8 5.667.197 × 10 10 7.213.757 × 10 10 7.262.171 × 10 10 7.43
641.072 × 10 10 6.634.430 × 10 12 7.342.126 × 10 12 7.471.183 × 10 12 7.52
1287.944 × 10 14 10.406.140 × 10 14 6.175.207 × 10 14 5.354.952 × 10 14 4.58
Table 7. Errors E 0 and O r d of Example 2 with solution C2 ( μ = 1.4 , α = β = 0 , δ = 0 , κ = 1 , λ = 1 ).
Table 7. Errors E 0 and O r d of Example 2 with solution C2 ( μ = 1.4 , α = β = 0 , δ = 0 , κ = 1 , λ = 1 ).
σ = 1.2 σ = 1.6 σ = 1.9 σ = 2.4
N E 0 Ord E 0 Ord E 0 Ord E 0 Ord
88.700 × 10 2 -2.008 × 10 2 -2.608 × 10 3 -2.659 × 10 3 -
162.745 × 10 2 1.663.490 × 10 3 2.522.875 × 10 4 3.181.369 × 10 4 4.28
329.543 × 10 3 1.526.934 × 10 4 2.333.740 × 10 5 2.948.736 × 10 6 3.97
643.160 × 10 3 1.591.317 × 10 4 2.404.672 × 10 6 3.005.404 × 10 7 4.01
1281.049 × 10 3 1.592.509 × 10 5 2.395.866 × 10 7 2.993.380 × 10 8 4.00
2563.468 × 10 4 1.604.768 × 10 6 2.407.357 × 10 8 3.002.120 × 10 9 3.99
5121.145 × 10 4 1.609.050 × 10 7 2.409.206 × 10 9 3.001.476 × 10 10 3.84
Table 8. Errors E 0 and O r d of Example 2 with solution C2 ( μ = 1.33 , α = β = 0 , δ = σ 1 , κ = 2 , λ = 1 ).
Table 8. Errors E 0 and O r d of Example 2 with solution C2 ( μ = 1.33 , α = β = 0 , δ = σ 1 , κ = 2 , λ = 1 ).
σ = 1.2 σ = 1.6 σ = 1.9 σ = 2.4
N E 0 Ord E 0 Ord E 0 Ord E 0 Ord
81.163 × 10 5 -1.331 × 10 5 -1.202 × 10 5 -1.302 × 10 5 -
161.165 × 10 7 6.649.526 × 10 8 7.136.495 × 10 8 7.537.022 × 10 8 7.53
321.180 × 10 9 6.636.321 × 10 10 7.243.020 × 10 10 7.754.644 × 10 10 7.24
641.122 × 10 11 6.723.846 × 10 12 7.361.316 × 10 12 7.841.643 × 10 11 4.82
1281.047 × 10 13 6.745.762 × 10 14 6.065.529 × 10 14 4.571.715 × 10 13 6.58
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhao, T. Efficient Spectral Collocation Method for Tempered Fractional Differential Equations. Fractal Fract. 2023, 7, 277. https://doi.org/10.3390/fractalfract7030277

AMA Style

Zhao T. Efficient Spectral Collocation Method for Tempered Fractional Differential Equations. Fractal and Fractional. 2023; 7(3):277. https://doi.org/10.3390/fractalfract7030277

Chicago/Turabian Style

Zhao, Tinggang. 2023. "Efficient Spectral Collocation Method for Tempered Fractional Differential Equations" Fractal and Fractional 7, no. 3: 277. https://doi.org/10.3390/fractalfract7030277

APA Style

Zhao, T. (2023). Efficient Spectral Collocation Method for Tempered Fractional Differential Equations. Fractal and Fractional, 7(3), 277. https://doi.org/10.3390/fractalfract7030277

Article Metrics

Back to TopTop