Next Article in Journal
Some Properties of Fractal Tsallis Entropy
Next Article in Special Issue
Total Controllability for a Class of Fractional Hybrid Neutral Evolution Equations with Non-Instantaneous Impulses
Previous Article in Journal
Fractional-Order Windkessel Boundary Conditions in a One-Dimensional Blood Flow Model for Fractional Flow Reserve (FFR) Estimation
Previous Article in Special Issue
The Rates of Convergence for Functional Limit Theorems with Stable Subordinators and for CTRW Approximations to Fractional Evolutions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Applications of Fractional Differentiation Matrices in Solving Caputo Fractional Differential Equations

Department of Mathematics, Shanghai University, Shanghai 200444, China
*
Author to whom correspondence should be addressed.
Submission received: 22 March 2023 / Revised: 16 April 2023 / Accepted: 27 April 2023 / Published: 30 April 2023

Abstract

:
This paper pursues obtaining Jacobi spectral collocation methods to solve Caputo fractional differential equations numerically. We used the shifted Jacobi–Gauss–Lobatto or Jacobi–Gauss–Radau quadrature nodes as the collocation points and derived the fractional differentiation matrices for Caputo fractional derivatives. With the fractional differentiation matrices, the fractional differential equations were transformed into linear systems, which are easier to solve. Two types of fractional differential equations were used for the numerical simulations, and the numerical results demonstrated the fast convergence and high accuracy of the proposed methods.

1. Introduction

The Caputo fractional derivative is used to model phenomena that involve interactions with the past and also problems with nonlocal properties [1,2,3,4,5]. The most-popular definition of the left-side Caputo derivative [1] is given below.
Definition 1.
Suppose that function u is defined in ( 0 , l ) , n 1 ν < n , n N , then the left-side Caputo fractional derivative of order ν is defined by
x ν u ( x , t ) = 0 C D x ν u ( x ) : = 1 Γ ( n ν ) 0 x ( x t ) n ν 1 u ( n ) ( t ) d t .
Given the above definition, it can be observed that the fractional derivative’s value at x 0 depends on the function values at x < x 0 . In this sense, one can think of the derivative as having “memory”. Thus, the fractional partial differential equations (FDEs) are more suitable than the traditional partial differential equations to study the processes of memory, hereditary properties, and heterogeneous materials.
For the fractional differential equations with an initial value, the existence and uniqueness of the solutions were studied in [6]. To achieve the analytical solutions, the Adomian decomposition methods [7] and the variational iteration methods [8,9] are the two common methods. The above literature showed that the exact solutions of FDEs are usually very hard to acquire or contain some special functions. Therefore, more and more researchers turn to numerical methods, for example finite-difference methods [10,11,12] and spectral methods [13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Spectral methods are well-known high-accuracy methods [27,28,29,30,31,32,33,34,35], normally carried out by a Galerkin, Tau, or collocation approach in practical problems. Because fractional derivative operators are nonlocal and spectral methods are global methods, it is very natural to apply spectral methods to solve FDEs.
Recently, the methods of operational matrices combined with spectral techniques have been widely used to solve FDEs due to their efficiency, stability, and usability. Suppose that u ( x ) is the solution of an FDE and { φ n ( x ) } n = 0 N are orthogonal polynomials. The solution u ( x ) can be approximated by u N ( x ) = n = 0 N c n φ n ( x ) . The operational matrix O of the fractional derivative is a matrix satisfying
x ν Φ ( x ) O Φ ( x ) ,
with Φ ( x ) = [ φ 0 ( x ) , φ 1 ( x ) , , φ N ( x ) ] T . By the operational matrices,
x ν u N ( x ) [ c 0 , c 1 , , c N ] O Φ ( x ) ,
the FDEs, hence, can be transformed into algebraic equations of coefficients c n , n = 0 , 1 , , N . The operational matrices based on all kinds of shifted orthogonal polynomials have been discussed by scholars, such as Laguerre operational matrices [14,15], Legendre operational matrices [16], Bernstein operational matrices [17], Chebyshev operational matrices [18,19], Jacobi operational matrices [20,21], and Genocchi operational matrices [22]. The analysis of the exponential convergence of the operational matrix methods for some FDEs can be found in [36,37,38].
Suppose that { x n } n = 0 N are a collection of N + 1 collocation points and u n = u N ( x n ) , n = 0 , 1 , , N . The approximation u N can be written in the Lagrange interpolating formula:
u N ( x ) = n = 0 N u n h n ( x ) ,
with h j ( x i ) = δ i j .
Let u = [ u 0 , u 1 , , u N ] T and x u = [ x ν u N ( x 0 ) , x ν u N ( x 1 ) , , x ν u N ( x N ) ] T ; it then follows that
x u = c D s ( ν ) u ,
where the matrix c D s ( ν ) = ( x ν h j ( x i ) ) is the fractional differentiation matrix (FDM) of the Caputo fractional derivative of order ν . The Lagrange cardinal polynomial can be formulated by h j ( x ) = i = 0 N t i j φ i ( x ) . The corresponding operational matrix, therefore, can be used to obtain the fractional differentiation matrix as
c D s ( ν ) [ Φ ( x 0 ) , Φ ( x 1 ) , , Φ ( x N ) ] O T ,
with T = ( t i j ) . However, to achieve the operational matrices, one often uses a polynomial expansion to approximate the non-integer power of x resulting from the fractional order derivative, which lowers the accuracy of the approaches. In this paper, we focused on deriving a new fractional differentiation matrix based on Jacobi polynomials to solve FDEs numerically.
This paper is structured as follows. Section 2 provides the preliminaries for standard Jacobi spectral collocation methods. In Section 3, the shifted Jacobi polynomials are introduced and the expressions of the fractional differentiation matrices for the fractional derivative are deduced using the shifted Jacobi–Gauss–Lobatto (JGL) and shifted Jacobi–Gauss–Radau (JGR) collocation points. The fractional differentiation matrix based on the shifted JGL points is used to numerically solve the generalized space–fractional Burgers’ equations in Section 4, and the fractional differentiation matrix based on the shifted JGR points is used to numerically solve fractional integro-differential equations in Section 5. The conclusion is made in the last section.

2. Jacobi Polynomials and Differentiation Matrix

Let P n α , β ( · ) stand for the common nth-degree Jacobi polynomials defined on [ 1 , 1 ] . Jacobi polynomials P n α , β ( · ) for n N 0 are orthogonal in regard to the weight function ω α , β ( x ) = ( 1 x ) α ( 1 + x ) β with α > 1 and β > 1 . From [34], P n α , β ( x ) has the following representation:
P n α , β ( x ) = Γ ( n + α + 1 ) n ! Γ ( n + α + β + 1 ) k = 0 n ( 1 ) k n k Γ ( n + k + α + β + 1 ) Γ ( k + α + 1 ) 1 x 2 k .
Suppose that { x j α , β } j = 0 N are the corresponding Jacobi–Gauss–Lobatto (JGL) quadrature nodes ( x 0 α , β = 1 , x N α , β = 1 , { x i α , β } i = 1 N 1 are the zeros of x P N α , β ( x ) ) or Jacobi–Gauss–Radau (JGR) quadrature nodes ( x 0 α , β = 1 and { x i α , β } i = 1 N are the zeros of P N α , β + 1 ( x ) ) and taken as the collocation points.
Let function u be defined on [ 1 , 1 ] and { h j α , β ( x ) } j = 0 N be Lagrange cardinal polynomials; the interpolation of u then is
u N ( x ) = j = 0 N u ( x j ) h j α , β ( x ) .
The Lagrange cardinal polynomials h j α , β ( x ) can be expressed by the linear combination of Jacobi polynomials as
h j α , β ( x ) = i = 0 N t i j α , β P i α , β ( x ) ,
with the coefficients:
t i j α , β = ω j α , β γ i α , β P i α , β ( x j ) ,
where
ω j α , β = 1 1 h j α , β ( x ) ω α , β ( x ) d x ,
and
γ i α , β = 1 1 [ P i α , β ( x ) ] 2 ω α , β ( x ) d x , i = 0 , , N 1 , 1 1 [ P N α , β ( x ) ] 2 ω α , β ( x ) d x , i = N , if JGR quadrature nodes are used , ( 2 + α + β + 1 N ) 1 1 [ P N α , β ( x ) ] 2 ω α , β ( x ) d x , i = N , if JGL quadrature nodes are used .
The coefficients { γ i α , β } i = 0 N and { ω j α , β } j = 0 N can be calculated as below.
When the JGL collocation points are used,
γ i α , β = 2 α + β + 1 Γ ( i + α + 1 ) Γ ( i + β + 1 ) ( 2 i + α + β + 1 ) i ! Γ ( i + α + β + 1 ) , i = 0 , , N 1 , ( 2 + α + β + 1 N ) 2 α + β + 1 Γ ( N + α + 1 ) Γ ( N + β + 1 ) ( 2 N + α + β + 1 ) N ! Γ ( N + α + β + 1 ) , i = N ,
and
ω j α , β = 2 α + β + 1 ( β + 1 ) Γ 2 ( β + 1 ) Γ ( N ) Γ ( N + α + 1 ) Γ ( N + β + 1 ) Γ ( N + α + β + 2 ) , j = 0 , G ˜ N 2 α + 1 , β + 1 ( 1 x j 2 ) 2 ( x P N 1 α + 1 , β + 1 ( x j ) ) 2 , j = 1 , , N 1 , 2 α + β + 1 ( α + 1 ) Γ 2 ( α + 1 ) Γ ( N ) Γ ( N + β + 1 ) Γ ( N + α + 1 ) Γ ( N + α + β + 2 ) , j = N ,
with G ˜ N 2 α + 1 , β + 1 = 2 α + β + 3 Γ ( N + α + 1 ) Γ ( N + β + 1 ) ( N 1 ) ! Γ ( N + α + β + 2 ) ; when the JGR collocation points are used,
γ i α , β = 2 α + β + 1 Γ ( i + α + 1 ) Γ ( i + β + 1 ) ( 2 i + α + β + 1 ) i ! Γ ( i + α + β + 1 ) , i = 0 , , N
and
ω j α , β = 2 α + β + 1 ( β + 1 ) Γ 2 ( β + 1 ) Γ ( N + 1 ) Γ ( N + α + 1 ) Γ ( N + β + 2 ) Γ ( N + α + β + 2 ) , j = 0 , G ˜ N 1 α , β + 1 ( 1 + x j ) ( 1 x j 2 ) ( x P N α , β + 1 ( x j ) ) 2 , j = 1 , , N ,
with G ˜ N 1 α , β + 1 = 2 α + β + 2 Γ ( N + α + 1 ) Γ ( N + β + 2 ) N ! Γ ( N + α + β + 2 ) .
After differentiation on Equation (3), one has
x u N ( x ) = j = 0 N u ( x j ) x h j α , β ( x ) .
The vector with its components being the values of the first-order derivative of u N at the collocation points can be retrieved by a matrix D multiplying a vector composed of function values of u N , i.e.,
u = D u ,
where
u = u ( x 0 α , β ) , u ( x 1 α , β ) , , u ( x N α , β ) T , u = x u ( x 0 α , β ) , x u ( x 1 α , β ) , , x u ( x N α , β ) T , D = d i j = x h j α , β ( x i α , β ) .
The matrix D is the first-order differentiation matrix, and the exact formula of D for the JGL or JGR collocation points can be found in [34]. The differentiation matrix of higher-integer-order derivatives of u N can be obtained by calculating the corresponding power of D . For example, the differentiation matrix of the second-order derivative is D 2 .

3. Jacobi Fractional Differentiation Matrix

The Caputo fractional differentiation matrix at the shifted Jacobi collocation points is derived in this section. The matrix will be used to solve FDEs in the next sections.
The shifted Jacobi polynomials, P 0 , l , n α , β ( x ) , n = 0 , 1 , 2 , , which are defined on [ 0 , l ] , can be obtained by a translation and dilation of general Jacobi polynomials, i.e.,
P 0 , l , n α , β ( x ) = P n α , β ( 2 x l 1 ) , n = 0 , 1 , 2 , .
The shifted JGL collocation points then become
y j α , β = l 2 ( x j α , β + 1 ) , j = 0 , 1 , , N .
According to (2), the analytic expression of shifted Jacobi polynomials is
P 0 , l , n α , β ( x ) = k = 0 n ( 1 ) n k Γ ( n + β + 1 ) Γ ( n + k + α + β + 1 ) k ! ( n k ) ! Γ ( n + α + β + 1 ) Γ ( k + β + 1 ) 1 l k x k .
Similarly, the approximation u N of function u defined on [ 0 , l ] has the following representation:
u N ( x ) = j = 0 N u ( y j α , β ) i = 0 N t i j α , β P 0 , l , i α , β ( x ) .
The differentiation matrix D s associated with the shifted Jacobi polynomials P 0 , l , n α , β ( x ) can be calculated by 2 D / l with D defined in Section 2.
Take the Caputo derivative of order ν on Equation (12); one has
x ν u N ( x ) = j = 0 N u ( y j α , β ) i = 0 N t i j α , β x ν P 0 , l , i α , β ( x ) .
According to [1,26], we have the following lemma for the Caputo derivative.
Lemma 1.
Suppose that C is a constant, then x ν C = 0 ; x ( 0 , l ) , ν ( n 1 , n ] , n N , then
x ν x m = 0 i f m N 0 , m < ν , Γ ( m + 1 ) Γ ( m + 1 ν ) x m ν , i f m N 0 , m ν , o r m N , m > ν .
It follows from the above lemma that
x ν P 0 , l , i α , β ( x ) = k = ν i ( 1 ) i k Γ ( i + β + 1 ) Γ ( i + k + α + β + 1 ) l k ( i k ) ! Γ ( k + 1 ν ) Γ ( i + α + β + 1 ) Γ ( k + β + 1 ) x k ν .
The formulas of the FDMs for Caputo fractional derivatives are, hence, presented in the following theorem.
Theorem 1.
The Jacobi–Caputo fractional differentiation matrix of order ν with the shifted JGL or JGR collocation points { y j α , β } j = 0 N is formulated as follows:
c D s ( ν ) : = ZA T ,
where T : = ( t i j α , β ) with t i j α , β defined in (4), Z : = z j m with for j = 0 , 1 , , N :
z j m = 0 , m = 0 , , ν 1 , ( y j α , β ) m ν , m = ν , , N ,
and A : = a m i with for i = 0 , 1 , , N ,
a m i = l m ( 1 ) i m Γ ( i + β + 1 ) Γ ( i + m + α + β + 1 ) ( i m ) ! Γ ( m + 1 ν ) Γ ( i + α + β + 1 ) Γ ( m + β + 1 ) m = ν , , i , 0 m = 0 , , ν 1 , 0 i < m N .
The Chebyshev polynomial of first kind is related to the Jacobi polynomials of α = β = 1 / 2 as
P n 1 2 , 1 2 ( x ) = Γ ( n + 1 2 ) Γ ( n + 1 ) Γ ( 1 2 ) T n ( x ) ,
where T n ( x ) = cos ( n arccos x ) is the nth-degree Chebyshev polynomial. The Chebyshev–Gauss–Lobatto (CGL) collocation points are
x j 1 2 , 1 2 = cos ( j π N ) , j = 0 , , N ,
and the Chebyshev–Gauss–Radau (CGR) collocation points are
x j 1 2 , 1 2 = cos ( 2 j π 2 N + 1 ) , j = 0 , , N .
The following corollary gives the formulas of the Chebyshev fractional differentiation matrix for the Caputo fractional derivative.
Corollary 1.
The Chebyshev–Caputo fractional differentiation matrix of order ν with the shifted CGL or CGR collocation points { y j 1 2 , 1 2 } j = 0 N is formulated as follows:
c D s ( ν ) : = ZA T C ,
where Z : = z j m with z j m defined in (14), A : = a m i with for i = 0 , 1 , , N :
a m i = l m ( 1 ) i m i Γ ( 1 / 2 ) Γ ( i + m ) ( i m ) ! Γ ( m + 1 ν ) Γ ( m + 1 / 2 ) m = ν , , i , 0 m = 0 , , ν 1 , 0 i < m N ,
and T C : = ( t i j ) with for i = 0 , 1 , , N :
t i j = ( 1 ) i 2 cos ( i j π N ) c i c j N , c 0 = c N = 2 , c i = 1 , i = 1 , , N 1
for the CGL collocation points, and
t i j = ( 1 ) i 4 cos ( 2 i j π 2 N + 1 ) c ˜ i c ˜ j ( 2 N + 1 ) , c ˜ 0 = 2 , c ˜ i = 1 , i = 1 , , N
for the CGR collocation points.
Proof. 
Using (4), (12) and (16), we have
u N ( x ) = j = 0 N u ( y j α , β ) i = 0 N C i 2 ω j γ i T i ( x j 1 2 , 1 2 ) T 0 , l , i ( x ) ,
where C i = Γ ( i + 1 2 ) Γ ( i + 1 ) Γ ( 1 2 ) and T 0 , l , i ( x ) is the shifted ith-degree Chebyshev polynomial. By the definition of γ i in (6), we have
γ i = C i 2 1 1 [ T i ( x ) ] 2 ( 1 x 2 ) 1 2 d x : = C i 2 γ ˜ i ,
thus for the CGL collocation points,
u N ( x ) = j = 0 N u ( y j α , β ) i = 0 N ω j γ ˜ i cos ( i j π N ) T 0 , l , i ( x ) ,
and for the CGR collocation points,
u N ( x ) = j = 0 N u ( y j α , β ) i = 0 N ω j γ ˜ i cos ( 2 i j π 2 N + 1 ) T 0 , l , i ( x ) .
According to [34], when the CGL collocation points are used,
ω j = π N , j = 1 , , N 1 , ω 0 = ω N = π 2 N , and γ ˜ i = π 2 , i = 1 , , N 1 , γ ˜ 0 = γ ˜ N = π ;
when the CGR collocation points are used,
ω 0 = π 2 N + 1 , ω j = 2 π 2 N + 1 , j = 1 , , N , and γ ˜ i = π 2 , i = 1 , , N , γ ˜ 0 = π .
Hence, the formula of the elements of T C is verified.
The Caputo fractional derivative of order ν of T 0 , l , i ( x ) is
x ν T 0 , l , i ( x ) = Γ ( 1 2 ) k = ν i l k ( 1 ) i k i Γ ( i + k ) ( i k ) ! Γ ( k + 1 ν ) Γ ( k + 1 2 ) x k ν ,
which verifies the the formula of the elements of Z , A and completes the proof. □
When α = β = 0 , the Jacobi polynomials P i 0 , 0 ( x ) become Legendre polynomials L i ( x ) . The Legendre–Gauss–Lobatto (LGL) quadrature nodes { x j 0 , 0 } j = 0 N are zeros of ( 1 x 2 ) x L N + 1 ( x ) , and the Legendre–Gauss–Radau (LGR) quadrature nodes { x j 0 , 0 } j = 0 N are zeros of L N ( x ) + L N + 1 ( x ) . The following corollary delivers the expression of Legendre–Caputo fractional differentiation matrix.
Corollary 2.
The Legendre–Caputo fractional differentiation matrix of order ν with the shifted LGL or LGR collocation points { y j 0 , 0 } j = 0 N is formulated as follows:
c D s ( ν ) : = ZA T L ,
where Z : = z j m with z j m defined in (14), A : = a m i with for i = 0 , 1 , , N :
a m i = l m ( 1 ) i m Γ ( i + m + 1 ) ( i m ) ! Γ ( m + 1 ν ) Γ ( m + 1 ) m = ν , , i , 0 m = 0 , , ν 1 , 0 i < m N ,
and T L : = ( t i j ) with for j = 0 , 1 , , N :
t i j = 2 i + 1 N ( N + 1 ) L N 2 ( x j ) L i ( x j ) , i = 0 , 1 , , N 1 , 1 ( N + 1 ) L N ( x j ) , i = N ,
for the LGL collocation points, and
t i j = ( 2 i + 1 ) ( 1 x j ) 2 ( N + 1 ) 2 L N 2 ( x j ) L i ( x j ) , i = 0 , 1 , , N ,
for the LGR collocation points.
Proof. 
The formula of the elements of Z and A can be easily verified by taking α = β = 0 in Theorem 1. According to [34], when the LGL collocation points are used,
ω j = 2 N ( N + 1 ) L N 2 ( x j ) , j = 0 , 1 , , N , and γ i = 2 2 i + 1 , i = 0 , , N 1 , γ N = 2 N ;
when the LGR collocation points are used,
ω j = 1 x j ( N + 1 ) 2 L N 2 ( x j ) , j = 0 , 1 , , N , and γ i = 2 2 i + 1 , i = 0 , , N .
Thus, the formula of the elements in T L is proven. □
Corollary 3.
The Caputo fractional differentiation matrix of order ν with the shifted JGL or JGR collocation points { y j α , 0 } j = 0 N is formulated as follows:
c D s ( ν ) : = R T ,
where T : = ( t i j α , 0 ) with t i j α , 0 defined in (4) by taking β = 0 and R : = ( r i j ) with
r i j = Γ ( i + 1 ) Γ ( i ν + 1 ) y j α , 0 ν P 0 , l , i α + ν , ν ( y j α , 0 ) k = 0 ν ( 1 ) i k i k Γ ( i + k + α + 1 ) l k Γ ( k + 1 ν ) Γ ( i + α + 1 ) y j α , 0 k ν , if i > ν , j > 0 , 0 , otherwise .
Proof. 
When x = 0 , x ν P 0 , l , i α , 0 ( 0 ) = 0 . When β = 0 , x 0 , and i > ν , we have
x ν P 0 , l , i α , 0 ( x ) = k = ν i ( 1 ) i k Γ ( i + 1 ) Γ ( i + k + α + 1 ) l k ( i k ) ! k ! Γ ( k + 1 ν ) Γ ( i + α + 1 ) x k ν = x ν Γ ( i + 1 ) Γ ( i ν + 1 ) k = ν i ( 1 ) i k Γ ( i ν + 1 ) Γ ( i + k + α + 1 ) l k ( i k ) ! k ! Γ ( k + 1 ν ) Γ ( i + α + 1 ) x k = x ν Γ ( i + 1 ) Γ ( i ν + 1 ) P 0 , l , i α + ν , ν ( x ) k = 0 ν ( 1 ) i k i k Γ ( i + k + α + 1 ) l k Γ ( k + 1 ν ) Γ ( i + α + 1 ) x k .
The proof is complete according to (12) and (18). □
When β = 0 and if 0 < ν < 1 , according to the proof of the above corollary, we have
x ν P 0 , l , i α , 0 ( x ) = x ν Γ ( i + 1 ) Γ ( i ν + 1 ) P 0 , l , i α + ν , ν ( x ) ( 1 ) i Γ ( 1 ν ) ,
then the element r i j , i > ν , j > 0 of matrix R in the above corollary has the following form:
Γ ( i + 1 ) Γ ( i ν + 1 ) y j α , 0 ν P 0 , l , i α + ν , ν ( y j α , 0 ) ( 1 ) i Γ ( 1 ν ) ( y j α , 0 ) ν .

4. Algorithm and Numerical Examples for the Generalized Space–Fractional Burgers’ Equations

The generalized space–fractional Burgers’ equation depicts the acoustic phenomenon of a sound wave with weak nonlinearity propagating through a gas-filled tube in one direction [12]. In this section, we intended to apply the Jacobi–Caputo fractional differentiation matrix derived in the previous section to solve the generalized space–fractional Burgers’ equation defined as follows:
t u ( x , t ) η x ν u ( x , t ) + ϵ u ( x , t ) x u ( x , t ) μ x 2 u ( x , t ) = g ( x ) , 0 < t < T , 0 < x < 1 , u ( 0 , t ) = u ( 1 , t ) = 0 , 0 < t < T , u ( x , 0 ) = u 0 ( x ) , 0 < x < 1 ,
where η , μ , and ϵ are arbitrary positive numbers and the fractional order of the Caputo derivative ν is a positive number less than 1.
Firstly, let y j α , β j = 0 N be a set of the shifted JGL collocation points. The shifted Jacobi spectral collocation method for (18) is to find u N ( · , t ) P N ( [ 0 , 1 ] ) , such that
t u N y j α , β , t = η x ν u N y j α , β , t ϵ u N y j α , β , t x u N y j α , β , t + μ x 2 u N y j α , β , t + g ( y j α , β , t ) , j = 1 , 2 , , N 1 , 0 < t < T , u N y 0 α , β , t = 0 , u N y N α , β , t = 0 , 0 < t < T , u N y j α , β , 0 = u 0 y j α , β , j = 0 , 1 , , N .
Let c D 1 : N 1 ( ν ) , D 1 : N 1 , and D 1 : N 1 ( 2 ) be ( N 1 ) × ( N 1 ) square matrices by removing the 1st, N + 1 st columns and rows, respectively from c D s ( ν ) , D , and D 2 . In accordance with (11), (19) can be rewritten in the following matrix–vector representation:
t u ( t ) = η c D 1 : N 1 ( ν ) u ( t ) ϵ u ( t ) · D 1 : N 1 u ( t ) + μ D 1 : N 1 ( 2 ) u ( t ) + g ( t ) , u ( 0 ) = u 0 ,
where
t u ( t ) = t u N y 1 α , β , t , t u N y 2 α , β , t , , t u N y N 1 α , β , t T , u ( t ) = u N y 1 α , β , t , u N y 2 α , β , t , , u N y N 1 α , β , t T , u 0 = u 0 y 1 α , β , u 0 y 2 α , β , , u 0 y N 1 α , β T , g ( t ) = g y 1 α , β , t , g y 2 α , β , t , , g y N 1 α , β , t T .
Two examples of the generalized space–fractional Burgers’ equation are given below.
Example 1.
In this example, the function on the right-hand side of (18) is chosen to be
g ( x , t ) = g 1 ( x , t ) + g 2 ( x , t ) + g 3 ( x , t ) + g 4 ( x , t ) , g 1 ( x , t ) = e σ x ( e σ 1 ) x 1 e t , g 2 ( x , t ) = ε e σ x ( e σ 1 ) x 1 σ e σ x e σ + 1 e 2 t , g 3 ( x , t ) = μ σ 2 e t + σ x , g 4 ( x , t ) = η σ ν x ν ν E 1 , ν + 1 ν ( σ x ) Γ ( 2 ) Γ ( 2 ν ) ( e σ 1 ) x 1 ν e t ,
where E a , b ( x ) is the two-parameter Mittag–Leffler function. The exact solution of (18) is e σ x ( e σ 1 ) x 1 e t .
Example 2.
Let g ( x , t ) in (18) be chosen as
g ( x , t ) = g 1 ( x , t ) + g 2 ( x , t ) + g 3 ( x , t ) + g 4 ( x , t ) , g 1 ( x , t ) = sin ( π x ) e t , g 2 ( x , t ) = ϵ π 2 sin ( 2 π x ) e 2 t , g 3 ( x , t ) = μ π 2 sin ( π x ) e t , g 4 ( x , t ) = η π x 1 ν ( 1 ν ) Γ ( 1 ν ) 1 F 2 ( 1 ; 1 ν 2 , 3 2 ν 2 ; 1 4 π 2 x 2 ) e t ,
where 1 F 2 is the generalized hypergeometric function. The exact solution of (18) is sin ( π x ) e t .
In the numerical simulations, we used the fourth-order Runge–Kutta method, a highly accurate algorithm, to solve the system of ODEs (20) and chose Δ t = 1.0 × 10 4 to make sure that the error resulting in the direction of time was small enough. We performed the numerical simulations with the fractional orders ν being 0.3 , 0.6 , and 0.9 and chose the CGL, LGL, and JGL( α = ν , β = 0 ) quadrature nodes as collocation points. In Table 1 and Table 2, we show the root-mean-squared (RMS) errors of the numerical solution at the collocation points at T = 1 of Examples 1 and 2, respectively. The RMS errors of these two numerical examples with ν = 0.6 using our method and the Jacobi operational matrix method proposed in [21] are plotted in Figure 1. The numerical results indicated the fast convergence of the methods and illustrated that our proposed algorithm in this section performed much better than the methods in [21].

5. Algorithm and Numerical Examples for Initial Fractional Integro-Differential Equations

In this section, we numerically solved the initial fractional integro-differential equation (FIDE) in the following form:
t ν u ( t ) = f ( t ) u ( t ) + 0 t K ( t , s ) u ( s ) d s + g ( t ) , 0 < t T , u ( 0 ) = u 0 ,
with the fractional order of the Caputo derivative ν ( 0 , 1 ) .
The algorithm developed in the section for (21) is to find u N ( t ) P N ( [ 0 , T ] ) , such that
t ν u N ( y m α , β ) = f ( y m α , β ) u N ( y m α , β ) + 0 y m α , β K ( y m α , β , s ) u N ( s ) d s + g ( y m α , β ) , m = 1 , , N , u ( 0 ) = u ( y 0 α , β ) = u 0 .
where y m α , β , m = 0 , , N , are the JGR points.
Let s ( θ , y m α , β ) = y m α , β 2 ( θ + 1 ) , 1 θ 1 . By the change of variable, one has
0 y m α , β K ( y m α , β , s ) u N ( s ) d s = y m α , β 2 1 1 K ( y m α , β , y m α , β 2 ( θ + 1 ) ) u N ( y m α , β 2 ( θ + 1 ) ) d θ .
Suppose that θ k 0 , 0 , k = 0 . , N is the LGR quadrature nodes and the integration can be approximated by a corresponding LGR quadrature rule with respect to the Jacobi weights ω k 0 , 0 , k = 0 . , N ; we hence have
t ν u N ( y m α , β ) = f ( y m α , β ) u N ( y m α , β ) + g ( y m α , β ) + y m α , β 2 k = 0 N K y m α , β , y m α , β 2 ( θ k 0 , 0 + 1 ) u N y m α , β 2 ( θ k 0 , 0 + 1 ) ω k 0 , 0 , m = 1 , , N , u ( y 0 α , β ) = u 0 .
Let u j = u ( y j α , β ) , f j = f ( y j α , β ) , g j = g ( y j α , β ) , j = 0 , , N . According to (12) and (13), we have
i = 0 N j = 1 N u j t i j k = ν i ( 1 ) i k Γ ( i + β + 1 ) Γ ( i + k + α + β + 1 ) L k ( i k ) ! Γ ( k + 1 ν ) Γ ( i + α + β + 1 ) Γ ( k + β + 1 ) y m α , β k ν = y m α , β 2 i = 0 N j = 1 N u j t i j k = 0 N K y m α , β , y m α , β ( θ k 0 , 0 + 1 ) 2 P 0 , L , i α , β y m α , β ( θ k 0 , 0 + 1 ) 2 ω k 0 , 0 + y m α , β 2 u 0 i = 0 N t i 0 k = 0 N K y m α , β , y m α , β ( θ k 0 , 0 + 1 ) 2 P 0 , L , i α , β y m α , β ( θ k 0 , 0 + 1 ) 2 ω k 0 , 0 u 0 i = 0 N t i 0 k = ν i ( 1 ) i k Γ ( i + β + 1 ) Γ ( i + k + α + β + 1 ) L k ( i k ) ! Γ ( k + 1 ν ) Γ ( i + α + β + 1 ) Γ ( k + β + 1 ) y m α , β k ν + f m u m + g m , m = 1 , , N .
In the next step, we worked on converting (24) into a vector–matrix form and started with the quadrature part first.
Let W = ( w m i ) with
w m i = y m α , β 2 k = 0 N K y m α , β , y m α , β ( θ k 0 , 0 + 1 ) 2 P 0 , L , i α , β y m α , β ( θ k 0 , 0 + 1 ) 2 ω k 0 , 0 .
The operational matrix for the quadrature part can, hence, be written as
M = W T ,
where the elements in T are defined in (4). Let f = ( f 0 , f 1 , , f N ) , u = u 0 , u 1 , , u N , and g = g 0 , g 1 , , g N . Suppose that c D 1 : N ( ν ) and M 1 : N are N × N square matrices by removing the first row and column, respectively, from c D s ( ν ) and M , while c d i ( ν ) and m i are the ( i + 1 ) th column of c D s ( ν ) and M without the first element; f 1 : N , g 1 : N , and u 1 : N are vectors formed by the last N elements of f , g , and u , respectively. The vector–matrix form of (24) is
c D 1 : N ( ν ) u 1 : N = f 1 : N · u 1 : N + M 1 : N u 1 : N + g 1 : N + u 0 m 0 u 0 c d 0 ( ν ) .
Example 3.
t ν u ( t ) = u ( t ) + 0 t ( t + s ) u ( s ) d s t 3 + t 2 Γ ( 3 ) Γ ( 3 ν ) t 2 ν + Γ ( 4 ) Γ ( 4 ν ) t 3 ν + 7 12 t 4 9 20 t 5 , u ( 0 ) = 0 ,
in association with the analytic solution u ( t ) = t 2 ( t 1 ) .
Example 4.
t ν u ( t ) = u ( t ) + 0 t ( t s ) u ( s ) d s + t 1 ν ( 1 ν ) Γ ( 1 ν ) 1 F 2 ( 1 ; 1 ν 2 , 3 2 ν 2 ; t 2 4 ) t , u ( 0 ) = 0 ,
with 1 F 2 being the generalized hypergeometric function and the exact solution u ( t ) = sin ( t ) .
We took the fractional order to be ν = 0.6 , used the LGR collocation points, and chose N = 16 in the numerical simulations. The absolute errors at the collocation points of Example 3 and 4 are plotted in the right panel in Figure 2 together with the numerical results produced by the methods proposed in [39]. The FIDE (21) was converted into an integral equation by taking the fractional-order integral and then approximated by the Gaussian quadrature rule with Jacobi collocation points in [39]. One can observe from the figure that the algorithm proposed in the section performs better compared to the methods in [39]. The numerical results matched well with the analytic solutions of the two examples, which demonstrated that the proposed scheme (26) worked accurately for the initial fractional integro-differential Equation (21).

6. Conclusions

In summary, Jacobi collocation spectral methods with the fractional differentiation matrices using shifted JGL or JGR collocation points were presented in this paper. We employed the matrices to numerically solve the generalized space–fractional Burgers’ equations and the initial fractional integro-differential equations. The numerical results of these two types of fractional differential equations fit very well with the analytic solutions, which indicated that the proposed approach can effectively solve the FDE problems.

Author Contributions

Conceptualization, X.Z. (Xiaoyan Zeng); methodology, X.Z. (Xiaoyan Zeng); software, X.Z. (Xiaoyan Zeng) and Z.W.; validation, X.Z. (Xiaoyan Zeng), Z.W., X.Z. (Xinxia Zhang) and J.W.; formal analysis, X.Z. (Xiaoyan Zeng); writing—original draft preparation, X.Z. (Xiaoyan Zeng), Z.W.; writing—review and editing, X.Z. (Xiaoyan Zeng) and Z.W.; supervision, X.Z. (Xiaoyan Zeng). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 11701358, 11774218.

Data Availability Statement

Not applicable.

Acknowledgments

The work was discussed with Heping Ma and Changpin Li.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

In this manuscript, the following abbreviations are used:
FDEsfractional partial differential equations
JGLJacobi–Gauss–Lobatto
JGRJacobi–Gauss–Radau
CGLChebyshev–Gauss–Lobatto
CGRChebyshev–Gauss–Radau
LGLLegendre–Gauss–Lobatto
LGRLegendre–Gauss–Radau
RMSroot mean square
FIDEfractional integro-differential equation

References

  1. Podulubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
  2. Chester, W.N. Resonant ocillations in closed tubes. J. Fluid Mech. 1964, 18, 44–64. [Google Scholar] [CrossRef]
  3. Keller, J.J. Propagation of simple non-linear waves in gas filled tubes with friction. Z. Ang. Math. Phy. 1982, 32, 170–181. [Google Scholar] [CrossRef]
  4. Kakutanil, T.; Matsuuchi, K. Effect of viscosity on long gravity waves. J. Phys. Soc. Jpn. 1975, 39, 237–246. [Google Scholar] [CrossRef]
  5. Miksis, M.J.; Ting, L. Effective equations for multiphase flows-Waves in a bubbly liquid. Adv. Appl. Mech. 1991, 28, 141–260. [Google Scholar]
  6. Biler, P.; Funaki, T.; Woyczynski, W.A. Fractal Burgers equations. J. Diff. Equ. 1998, 148, 9–46. [Google Scholar] [CrossRef] [Green Version]
  7. El-Shahed, M. Adomian decomposition method for solving Burgers equation with fractional derivative. J. Fract. Calc. 2003, 24, 23–28. [Google Scholar]
  8. Momani, S. Non-perturbative analytical solutions of the space- and time-fractional Burgers equations. Chaos Soliton. Fract. 2006, 28, 930–937. [Google Scholar] [CrossRef]
  9. Inc, M. The approximate and exact solutions of the space- and time-fractional Burgers equations with initial conditions by variational iteration method. J. Math. Anal. Appl. 2008, 345, 476–484. [Google Scholar] [CrossRef] [Green Version]
  10. Esen, A.; Bulut, F.; Oruç, Ö. A unified approach for the numerical solution of time fractional Burgers’ type equations. Euro. Phys. J. Plus 2016, 131, 116. [Google Scholar] [CrossRef]
  11. Li, D.; Zhang, C.; Ran, M. A linear finite difference scheme for generalized time fractional Burgers equation. Appl. Math. Model. 2016, 40, 6069–6081. [Google Scholar] [CrossRef]
  12. Sugimoto, N. Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves. J. Fluid Mech. 1991, 225, 631–653. [Google Scholar] [CrossRef]
  13. Khatera, A.H.; Temsaha, R.S.; Hassanb, M.M. A Chebyshev spectral collocation method for solving Burgers’-type equations. J. Comput. Appl. Math. 2008, 222, 333–350. [Google Scholar] [CrossRef] [Green Version]
  14. Hwang, C.; Shih, Y.-P. Laguerre operational matrices for fractional calculus and applications. Int. J. Control 1981, 34, 577–584. [Google Scholar] [CrossRef]
  15. Parand, K.; Dehghan, M.; Taghavi, A. Modified generalized Laguerre function Tau method for solving laminar viscous flow: The Blasius equation. Int. J. Numer. Method Heat 2010, 20, 728–743. [Google Scholar] [CrossRef]
  16. Saadatmandi, A.; Dehghan, M. A new operational matrix for solving fractional-order differential equations. Comput. Math. Appl. 2010, 59, 1326–1336. [Google Scholar] [CrossRef] [Green Version]
  17. Maleknejad, K.; Hashemizadeh, E.; Basirat, B. Computational method based on Bernstein operational matrices for nonlinear Volterra–Fredholm–Hammerstein integral equations. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 52–61. [Google Scholar] [CrossRef]
  18. Bhrawy, A.H.; Alofi, A.S. The operational matrix of fractional integration for shifted Chebyshev polynomials. Appl. Math. Lett. 2013, 26, 25–31. [Google Scholar] [CrossRef] [Green Version]
  19. Jafari, H.; Nemati, S.; Ganji, R.M. Operational matrices based on the shifted fifth-kind Chebyshev polynomials for solving nonlinear variable order integro-differential equations. Adv. Differ. Equ. 2021, 2021, 435. [Google Scholar] [CrossRef]
  20. Doha, E.H.; Bhrawy, A.H.; Ezz-Eldien, S.S. A new Jacobi operational matrix: An application for solving fractional differential equations. Appl. Math. Model. 2012, 36, 4931–4943. [Google Scholar] [CrossRef]
  21. Wu, Q.; Zeng, X. Jacobi Collocation Methods for Solving Generalized Space-Fractional Burgers’ Equations. Commun. Appl. Math. Comput. 2020, 36, 4931–4943. [Google Scholar] [CrossRef] [Green Version]
  22. Sadeghi, S.; Jafari, H.; Nemati, S. Operational matrix for Atangana–Baleanu derivative based on Genocchi polynomials for solving FDEs. Chaos Solit. 2020, 135, 106973. [Google Scholar] [CrossRef]
  23. Zayernouri, M.; Karniadakis, G.E. Fractional spectral collocation method. SIAM J. Sci. Comput. 2014, 36, A40–A62. [Google Scholar] [CrossRef]
  24. Bhrawy, A.H.; Taha, M.T.; Machado, J.A.T. A review of operational matrices and spectral techniques for fractional calculus. Nonlinear Dyn. 2015, 81, 1023–1052. [Google Scholar] [CrossRef]
  25. Yang, Y.B.; Ma, H.P. The Legendre Galerkin-Chebyshev Collocation method for generalized space–fractional Burgers equations. J. Numer. Meth. Comp. Appl. 2017, 38, 236–244. [Google Scholar]
  26. Afsane, S.; Mahmoud, B. Application Jacobi spectral method for solving the time-fractional differential equation. J. Comput. Appl. Math. 2018, 399, 49–68. [Google Scholar]
  27. Guo, B.Y. Spectral Methods and Their Applications; World Scientific: Singapore, 1998. [Google Scholar]
  28. Guo, B.Y.; Wang, L.L. Jacobi interpolation approximations and their applications to singular differential equations. Adv. Comput. Math. 2001, 14, 227–276. [Google Scholar]
  29. Ma, H.P.; Sun, W.W. Optimal Error Estimates of the Legendre-Petrov-Galerkin Method for the Korteweg-de Vries Equation. SIAM J. Numer. Anal. 2001, 39, 1380–1394. [Google Scholar] [CrossRef]
  30. Wu, H.; Ma, H.P.; Li, H.Y. Optimal error estimates of the Chebyshev-Legendre method for solving the generalized Burgers equation. SIAM J. Numer. Anal. 2003, 41, 659–672. [Google Scholar] [CrossRef]
  31. Guo, B.Y.; Wang, L.L. Jacobi approximations in non-uniformly Jacobi-weighted Sobolev spaces. J. Approx. Theory 2004, 128, 1–41. [Google Scholar] [CrossRef] [Green Version]
  32. Shen, J.; Tang, T. Spectral and High-Order Methods with Applications; Science Press: Beijing, China, 2006. [Google Scholar]
  33. Canuto, C.G.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods: Fundamentals in Single Domains; Springer: New York, NY, USA, 2010. [Google Scholar]
  34. Shen, J.; Tang, T.; Wang, L.L. Spectral Methods Algorithms, Analysis and Applications; Springer Series in Computational Mathematics: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  35. Zhao, T.G.; Wu, Y.J.; Ma, H.P. Error Analysis of Chebyshev-Legendre Pseudo-spectral Method for a Class of Nonclassical Parabolic Equation. SIAM J. Sci. Comput. 2012, 52, 588–602. [Google Scholar] [CrossRef]
  36. Mockary, S.; Babolian, E.; Vahidi, A.R. A fast numerical method for fractional partial differential equations. Adv. Differ. Equ. 2019, 2019, 452. [Google Scholar] [CrossRef] [Green Version]
  37. Wu, C.; Wang, Z. The spectral collocation method for solving a fractional integro-differential equation. AIMS Math. 2022, 7, 9577–9587. [Google Scholar] [CrossRef]
  38. Doha, E.H.; Abdelkawy, M.A.; Amin, A.Z.M.; Lopes, A.M. Shifted Jacobi–Gauss-collocation with convergence analysis for fractional integro-differential equations. Commun. Nonlinear Sci. Numer. Simul. 2019, 72, 342–359. [Google Scholar] [CrossRef]
  39. Wu, Q.; Wu, Z.; Zeng, X. A Jacobi Spectral Collocation Method for Solving Fractional Integro-Differential Equations. Commun. Appl. Math. Comput. 2021, 3, 509–526. [Google Scholar] [CrossRef]
Figure 1. Numerical results of Examples 1 and 2 when ν = 0.6 and α = β = 0 . I: results by using methods proposed in this paper; II: results by using methods proposed in [21].
Figure 1. Numerical results of Examples 1 and 2 when ν = 0.6 and α = β = 0 . I: results by using methods proposed in this paper; II: results by using methods proposed in [21].
Fractalfract 07 00374 g001
Figure 2. Numerical results of Examples 3 (top) and 4 (bottom) when ν = 0.6 and N = 16 . left: results by using methods proposed in this paper; right: results by using methods proposed in [39].
Figure 2. Numerical results of Examples 3 (top) and 4 (bottom) when ν = 0.6 and N = 16 . left: results by using methods proposed in this paper; right: results by using methods proposed in [39].
Fractalfract 07 00374 g002aFractalfract 07 00374 g002b
Table 1. RMS errors of numerical solutions of Example 1 at T = 1 when ϵ = μ = η = 1 , σ = 1 , and τ = 1.0 × 10 4 .
Table 1. RMS errors of numerical solutions of Example 1 at T = 1 when ϵ = μ = η = 1 , σ = 1 , and τ = 1.0 × 10 4 .
N ν = 0.3 ν = 0.6 ν = 0.9
46.0362 × 10 6 5.9249 × 10 6 6.2110 × 10 6
α = 1 2 83.4641 × 10 12 3.4508 × 10 12 3.5460 × 10 12
β = 1 2 121.0282 × 10 15 9.6164 × 10 16 9.4303 × 10 16
45.5072 × 10 7 1.3705 × 10 7 6.5986 × 10 7
α = 0 81.8781 × 10 13 2.6166 × 10 14 4.2505 × 10 13
β = 0 127.9866 × 10 16 7.7894 × 10 16 7.9360 × 10 16
42.5969 × 10 6 3.8212 × 10 6 5.6797 × 10 6
α = ν 81.7825 × 10 12 2.2046 × 10 12 3.1671 × 10 12
β = 0 128.8606 × 10 16 1.1035 × 10 15 1.5604 × 10 15
Table 2. RMS errors of numerical solutions of Example 2 at T = 1 when ϵ = μ = η = 1 and τ = 1.0 × 10 4 .
Table 2. RMS errors of numerical solutions of Example 2 at T = 1 when ϵ = μ = η = 1 and τ = 1.0 × 10 4 .
N ν = 0.3 ν = 0.6 ν = 0.9
48.8219 × 10 4 8.9389 × 10 4 9.2472 × 10 4
α = 1 2 81.8921 × 10 8 1.8993 × 10 8 1.9393 × 10 8
β = 1 2 121.7631 × 10 13 1.7707 × 10 13 1.8036 × 10 13
162.9591 × 10 15 3.6478 × 10 15 5.5501 × 10 15
43.1279 × 10 4 3.1274 × 10 4 3.2005 × 10 4
α = 0 89.6568 × 10 9 9.6444 × 10 9 9.7385 × 10 9
β = 0 129.6773 × 10 14 9.6694 × 10 14 9.7314 × 10 14
162.4675 × 10 15 2.6880 × 10 15 3.4963 × 10 15
44.0505 × 10 4 4.9560 × 10 4 6.1909 × 10 4
α = ν 81.2005 × 10 8 1.3271 × 10 8 1.4142 × 10 8
β = 0 121.2032 × 10 13 1.2645 × 10 13 1.3179 × 10 13
162.5721 × 10 15 4.3107 × 10 15 1.8257 × 10 14
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

Wu, Z.; Zhang, X.; Wang, J.; Zeng, X. Applications of Fractional Differentiation Matrices in Solving Caputo Fractional Differential Equations. Fractal Fract. 2023, 7, 374. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7050374

AMA Style

Wu Z, Zhang X, Wang J, Zeng X. Applications of Fractional Differentiation Matrices in Solving Caputo Fractional Differential Equations. Fractal and Fractional. 2023; 7(5):374. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7050374

Chicago/Turabian Style

Wu, Zhongshu, Xinxia Zhang, Jihan Wang, and Xiaoyan Zeng. 2023. "Applications of Fractional Differentiation Matrices in Solving Caputo Fractional Differential Equations" Fractal and Fractional 7, no. 5: 374. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7050374

Article Metrics

Back to TopTop