Next Article in Journal
Fractional Modeling Applied to the Dynamics of the Action Potential in Cardiac Tissue
Next Article in Special Issue
Topological Structure of the Solution Sets for Impulsive Fractional Neutral Differential Inclusions with Delay and Generated by a Non-Compact Demi Group
Previous Article in Journal
Consensus of Fractional-Order Double-Integral Multi-Agent System in a Bounded Fluctuating Potential
Previous Article in Special Issue
Existence Results for Coupled Nonlinear Sequential Fractional Differential Equations with Coupled Riemann–Stieltjes Integro-Multipoint Boundary Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Mixed Boundary Value Problems and Chebyshev Collocation Method for Caputo-Type Fractional Ordinary Differential Equations

1
School of Sciences, Shanghai Institute of Technology, Shanghai 201418, China
2
Ocean College, Zhejiang University, Zhejiang 310012, China
3
Village 1, East China Normal University, Shanghai 200062, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(3), 148; https://doi.org/10.3390/fractalfract6030148
Submission received: 11 February 2022 / Revised: 6 March 2022 / Accepted: 7 March 2022 / Published: 9 March 2022

Abstract

:
The boundary value problem (BVP) for the varying coefficient linear Caputo-type fractional differential equation subject to the mixed boundary conditions on the interval 0 x 1 was considered. First, the BVP was converted into an equivalent differential–integral equation merging the boundary conditions. Then, the shifted Chebyshev polynomials and the collocation method were used to solve the differential–integral equation. Varying coefficients were also decomposed into the truncated shifted Chebyshev series such that calculations of integrals were only for polynomials and can be carried out exactly. Finally, numerical examples were examined and effectiveness of the proposed method was verified.

1. Introduction

In recent decades, the theory of fractional calculus has been attracting much attention partly due to its ability for describing memory and hereditary properties of various materials and processes [1,2,3,4,5,6,7]. Fractional calculus has been applied to different fields such as viscoelastic constitutive equations and related mechanical models [6,7,8,9,10,11], anomalous diffusion phenomena [4,6,12,13], hydrology [14], control and optimization theory [3,15], etc. It is worthwhile to mention that fractional calculus can be used to describe not only viscoelasticity, but also viscoinertia by different values of order [16,17]. The applications of fractional calculus lead to fractional differential equations (FDEs) in theory [2,3,4,5,18].
Let us recall some related definitions of fractional calculus used in this article. Let f ( x ) be piecewise continuous on ( 0 , + ) and integrable on any finite subinterval of ( 0 , + ) . Then, for x > 0 , the Riemann–Liouville fractional integral of f ( x ) is defined as
I x β f ( x ) = 0 x ( x τ ) β 1 Γ ( β ) f ( τ ) d τ ,
for β > 0 , and I x 0 f ( x ) = f ( x ) for β = 0 , where Γ ( · ) is the gamma function. The fractional integral satisfies the following equalities:
I x β I x ν f ( x ) = I x β + ν f ( x ) , β 0 , ν 0 ,
I x ν x μ = Γ ( μ + 1 ) Γ ( μ + ν + 1 ) x μ + ν , ν 0 , μ > 1 .
Let α be a positive real number, m 1 < α m , m N + , and f ( m ) ( x ) be piecewise continuous on ( 0 , + ) and integrable on any finite subinterval of ( 0 , + ) . Then, the Caputo fractional derivative of f ( x ) of order α is defined as
D x α f ( x ) = I x m α f ( m ) ( x ) , m 1 < α m .
For the power function x μ , μ > 0 , if 0 m 1 < α m < μ + 1 , then we have
D x α x μ = Γ ( μ + 1 ) Γ ( μ α + 1 ) x μ α , x > 0 .
The α -order integral of the α -order Caputo fractional derivative requires the knowledge of the initial values of the function and its integer-ordered derivatives,
I x α D x α f ( x ) = f ( x ) k = 0 m 1 f ( k ) ( 0 ) x k k ! , m 1 < α m .
This property enables the Caputo fractional derivative to be conveniently applied and analyzed.
In the earlier monograph [1], the Grünwald definition and the Riemann–Liouville definition of fractional calculus were introduced, where numerical differentiation and integration were considered and semi-integration was introduced by a designed electrical circuit model and semi-differentiation was applied to diffusion problems. The Weyl fractional calculus was introduced in [2] beside the Grünwald definition and the Riemann–Liouville definition. In [3], FDEs and fractional-order system and controllers were considered, where the Caputo fractional derivative was introduced. The existence, uniqueness and analytical methods of solutions for FDEs were investigated in [4]. In [18], the Caputo-type fractional derivative and FDEs were emphasized. In [6], fractional viscoelastic models and fractional wave models in viscoelastic media were introduced. In [5], numerical methods and fractional variational principle were reviewed.
Damping, deformation, vibration and dissipation arising from viscoelastic material can be modeled by FDEs [3,4,6,7]. The method of variable separation for fractional partial differential equation describing anomalous diffusion [4,6,12,14] can lead to a boundary value problem (BVP) for a fractional ordinary differential equation (ODE) [19]. The theorem of existence and uniqueness of solutions for fractional ODEs was presented in [3,4,18,20]. Some analytical and numerical methods were proposed to solve FDEs, e.g., see [3,4,5,21,22,23,24,25]. BVPs for fractional ODEs were considered in [19,26,27,28,29] by using the Adomian decomposition method, wavelet method, the method of upper and lower solutions, orthogonal polynomial method, etc. However, a fractional BVP with varying coefficients and mixed boundary conditions has hardly been considered.
In this work, we consider the BVP for the varying coefficient linear Caputo fractional ODE
D x λ u ( x ) + c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) = g ( x ) , 0 < x < 1 , 1 < λ 2 ,
Subject to the mixed boundary conditions
p 0 u ( 0 ) q 0 u ( 0 ) = b 0 ,
p 1 u ( 1 ) + q 1 u ( 1 ) = b 1 ,
where the coefficients c 1 ( x ) , c 0 ( x ) , g ( x ) are specified continuous functions, the boundary parameters satisfy p 0 , q 0 , p 1 , q 1 0 and p 0 p 1 + p 0 q 1 + q 0 p 1 0 . In the next Section 2, some preliminaries about the shifted Chebyshev polynomials are presented. In Section 3, we first convert the BVP, (7)–(9), into an equivalent fractional differential–integral equation merging the boundary conditions, then introduce the collocation method using the shifted Chebyshev polynomials of the first kind to solve the fractional differential–integral equation. Next, three numerical examples are solved by using the proposed method. Section 4 summarizes our conclusions.

2. The Shifted Chebyshev Polynomials of the First Kind

The Chebyshev polynomials of the first kind are defined by the formulae [30]
T n ( x ) = cos ( n arccos x ) , 1 x 1 , n = 0 , 1 , .
They take on the explicit expressions as
T 0 ( x ) = 1 , T n ( x ) = n 2 k = 0 [ n / 2 ] ( 1 ) k ( n k 1 ) ! k ! ( n 2 k ) ! ( 2 x ) n 2 k , n 1 .
It is well-known that the Chebyshev polynomials of the first kind are orthogonal on the interval [ 1 , 1 ] with the weight function ρ ( x ) = 1 1 x 2 , and T n ( x ) has exactly n zeros within the interval ( 1 , 1 ) : ξ i = cos 2 i + 1 2 n π , i = 0 , 1 , , n 1 . The Chebyshev polynomials of the first kind satisfy the recurrence relation
T 0 ( x ) = 1 , T 1 ( x ) = x , T n ( x ) = 2 x T n 1 ( x ) T n 2 ( x ) , n = 2 , 3 , .
It is well-known that if f ( x ) is L 2 integrable on [ 1 , 1 ] with the weight function ρ ( x ) , then its Chebyshev series expansion is L 2 convergent with respect to its weight function ρ ( x ) . If f ( x ) has better smoothness, then stronger convergence can be attained for its Chebyshev series. If the function f ( x ) has n + 1 continuous derivatives on [ 1 , 1 ] , then | f ( x ) S m f ( x ) | = O ( m n ) for all x [ 1 , 1 ] , where S m f ( x ) is the ( m + 1 ) -term truncation of the Chebyshev series expansion of f ( x ) . For more details for convergence, see [30].
In order to deal with the BVP on the interval [ 0 , 1 ] , we consider the shifted Chebyshev polynomials
T n * ( x ) = T n ( 2 x 1 ) , x [ 0 , 1 ] , n = 0 , 1 , .
They are orthogonal on the interval [ 0 , 1 ] with the weight function ρ * ( x ) = 1 x x 2 , and the zeros of T n * ( x ) are x i = 1 2 + 1 2 cos 2 i + 1 2 n π ,   i = 0 , 1 , , n 1 . As a complement to Equation (13), the shifted Chebyshev polynomials satisfy the relationship T n * ( x ) = T 2 n ( x ) . So, the explicit expressions of the shifted Chebyshev polynomials are conveniently obtained:
T 0 * ( x ) = 1 , T n * ( x ) = n k = 0 n ( 1 ) k ( 2 n k 1 ) ! k ! ( 2 n 2 k ) ! ( 4 x ) n k , n 1 .
Finally, we mention the shifted Chebyshev polynomials of the second kind, which will also be used in the next section for the representation of solutions, U n * ( x ) = U n ( 2 x 1 ) , 0 x 1 , n = 0 , 1 , , where U n ( x ) is the Chebyshev polynomials of the second kind.

3. The Equivalent Fractional Differential-Integral Equation and Chebyshev Collocation Method

First, we derive an equivalent differential–integral equation to the BVP (7)–(9). Applying the integral operator I x λ ( · ) to both sides of Equation (7) and using Equation (6) yields
u ( x ) u ( 0 ) u ( 0 ) x + I x λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) = I x λ g ( x ) .
Our aim is to solve for u ( 0 ) and u ( 0 ) from the boundary conditions (8) and (9), and then obtain an equation about the solution u ( x ) without any undetermined constants. Substituting x = 1 in Equation (15) yields
u ( 1 ) = u ( 0 ) + u ( 0 ) I x , 1 λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) + I x , 1 λ g ( x ) ,
where the value of the fractional integral is defined for a general β th order integral of a function v ( x ) at x = ξ as
I x , ξ β v ( x ) = 0 ξ ( ξ τ ) β 1 Γ ( β ) v ( τ ) d τ .
Calculating the first order derivative on the both sides of Equation (15) leads to
u ( x ) u ( 0 ) + I x λ 1 ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) = I x λ 1 g ( x ) .
Substituting x = 1 yields
u ( 1 ) = u ( 0 ) I x , 1 λ 1 ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) + I x , 1 λ 1 g ( x ) .
Substituting Equations (16) and (19) into Equation (9) yields
p 1 u ( 0 ) + ( p 1 + q 1 ) u ( 0 ) = b 1 * ,
where
b 1 * = b 1 + p 1 I x , 1 λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) p 1 I x , 1 λ g ( x ) + q 1 I x , 1 λ 1 ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) q 1 I x , 1 λ 1 g ( x ) .
Equations (8) and (20) constitute a system of algebraic equations about u ( 0 ) and u ( 0 ) . The coefficient determinant is
P = p 0 p 1 + p 0 q 1 + q 0 p 1 ,
which is positive by our assumptions. Thus, we can solve the system of algebraic Equations (8) and (20) and obtain
u ( 0 ) = ( p 1 + q 1 ) b 0 P + q 0 b 1 * P ,
u ( 0 ) = p 1 b 0 P + p 0 b 1 * P .
Substituting Equations (23) and (24) into Equation (15), we obtain
u ( x ) ( p 1 + q 1 ) b 0 P + p 1 b 0 P x p 0 x + q 0 P b 1 * + I x λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) = I x λ g ( x ) .
Replacing b 1 * by using Equation (21) and reorganizing the equation yield
u ( x ) p 1 ( p 0 x + q 0 ) P I x , 1 λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) q 1 ( p 0 x + q 0 ) P I x , 1 λ 1 ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) + I x λ ( c 1 ( x ) u ( x ) + c 0 ( x ) u ( x ) ) = h ( x ) ,
where
h ( x ) = ( p 1 + q 1 ) b 0 p 1 b 0 x P + p 0 x + q 0 P b 1 p 1 I x , 1 λ g ( x ) q 1 I x , 1 λ 1 g ( x ) + I x λ g ( x ) ,
Only involves the known boundary parameters and the known input function g ( x ) . Equation (26) is the equivalent differential–integral equation to the BVP (7)–((9). In the sequel, we seek for the solution to the differential-integral Equation (26).
We approximate the solution by an ( m + 1 ) -term truncation of the shifted Chebyshev series,
φ m ( x ) = n = 0 m a n T n * ( x ) ,
where a n , n = 0 , 1 , , m , are undetermined coefficients. Inserting φ m ( x ) into Equation (26), we obtain the linear equation about a n , n = 0 , 1 , , m ,
n = 0 m a n T n * ( x ) p 1 ( p 0 x + q 0 ) P I x , 1 λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) q 1 ( p 0 x + q 0 ) P I x , 1 λ 1 c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) + I x λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) = h ( x ) .
We note that in Equation (29), I x , 1 λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) and I x , 1 λ 1 c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) are constants, represent the values of fractional integrals.
The collocation method may be applied to determine the coefficients a n . The collocation points are taken as the zeroes of the m + 1 degree shifted Chebyshev polynomial T m + 1 * ( x ) ,
x i = 1 2 + 1 2 cos 2 i + 1 2 m + 2 π , i = 0 , 1 , , m .
Thus, the collocation equation system is
n = 0 m a n T n * ( x i ) p 1 ( p 0 x i + q 0 ) P I x , 1 λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) q 1 ( p 0 x i + q 0 ) P I x , 1 λ 1 c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) + I x , x i λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) = h ( x i ) ,
where
h ( x i ) = ( p 1 + q 1 ) b 0 p 1 b 0 x i P + p 0 x i + q 0 P b 1 p 1 I x , 1 λ g ( x ) q 1 I x , 1 λ 1 g ( x ) + I x , x i λ g ( x ) , i = 0 , 1 , , m .
The matrix form of the collocation equation system (31) is
W a = h ,
where
a = ( a 0 , a 1 , , a m ) T , h = ( h ( x 0 ) , h ( x 1 ) , , h ( x m ) ) T ,
and the entries of the matrix W are
w i j = T j * ( x i ) p 1 ( p 0 x i + q 0 ) P I x , 1 λ c 1 ( x ) T j * ( x ) + c 0 ( x ) T j * ( x ) q 1 ( p 0 x i + q 0 ) P I x , 1 λ 1 c 1 ( x ) T j * ( x ) + c 0 ( x ) T j * ( x ) + I x , x i λ c 1 ( x ) T j * ( x ) + c 0 ( x ) T j * ( x ) , i , j = 0 , 1 , , m .
The solution of the linear algebraic equation system (31) or (33) gives the coefficients a n in Equation (28).
For the Dirichlet boundary conditions u ( 0 ) = b 0 , u ( 1 ) = b 1 , the boundary parameters are simplified as p 0 = p 1 = 1 and q 0 = q 1 = 0 , and thus Equation (31) degenerates to
n = 0 m a n T n * ( x i ) x i I x , 1 λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) + I x , x i λ c 1 ( x ) T n * ( x ) + c 0 ( x ) T n * ( x ) = h ( x i ) ,
where h ( x i ) = b 0 + ( b 1 b 0 ) x i x i I x , 1 λ g ( x ) + I x , x i λ g ( x ) , i = 0 , 1 , , m .
We remark that by the relationship of the first-kind and second-kind Chebyshev polynomials T n ( x ) = n U n 1 ( x ) , we have the relationship of the shifted Chebyshev polynomials of the two kinds
T n * ( x ) = 2 n U n 1 * ( x ) .
So, the derivative T n * ( x ) in Equations (31), (35) and (36) may be replaced by the second-kind Chebyshev polynomials.
The operators I x , 1 λ ( · ) , I x , 1 λ 1 ( · ) and I x , x i λ ( · ) in Equations (31), (32) and (35) represent the values of fractional integrals of the known functions. Since the appearance of the varying coefficients g ( x ) and c k ( x ) , manual computations for these integrals are laborious in general. Here we approximate the varying coefficients again using their truncated shifted Chebyshev series as
g ( x ) = n = 0 M g n T n * ( x ) , c k ( x ) = n = 0 M c k , n T n * ( x ) , k = 0 , 1 , 0 x 1 ,
where
g n = 2 π 0 1 1 x x 2 g ( x ) T n * ( x ) d x , n = 0 , 1 , , M ,
c k , n = 2 π 0 1 1 x x 2 c k ( x ) T n * ( x ) d x , k = 0 , 1 , n = 0 , 1 , , M ,
and the superscript of ∑ denotes that the first term in the sum is halved. We note that there is no need of connections between the values of m and M in Equations (28) and (38). Utilizing the Gauss–Chebyshev quadrature formula we derive the numerical formulae for g n and c k , n as
g n = 2 M + 1 i = 0 M g ( x i ) T n * ( x i ) , n = 0 , 1 , , M ,
c k , n = 2 M + 1 i = 0 M c k ( x i ) T n * ( x i ) , k = 0 , 1 , n = 0 , 1 , , M ,
where x i are the zeroes of the M + 1 degree shifted Chebyshev polynomial T M + 1 * ( x ) ,
x i = 1 2 + 1 2 cos 2 i + 1 2 M + 2 π , i = 0 , 1 , , M .
Thus, making use of the decompositions in (38), the calculation of the integrals I x , 1 λ ( · ) , I x , 1 λ 1 ( · ) and I x , x i λ ( · ) in Equations (31), (32) and (35) only involves integrals of polynomials, so can be carried out exactly.
In the following three examples, we take M = 5 in Equation (38) to truncate the decompositions of the coefficients g ( x ) and c k ( x ) and to calculate the involved integrals I x , 1 λ ( · ) , I x , 1 λ 1 ( · ) and I x , x i λ ( · ) . Collocation equation systems are solved by using Mathematica command “LinearSolve”. Figures of approximate analytical solutions and errors are generated by using Mathematica.
Example 1.
Consider the BVP for the linear FDE
D x 1.5 u ( x ) x sin ( x ) 3 u ( x ) + sin ( x ) u ( x ) = g ( x ) , 0 < x < 1 ,
u ( 0 ) = 1 , u ( 1 ) = 1 ,
where g ( x ) = 3 π 4 + 8 x 1.5 π 2 3 x sin ( x ) + 1 2 x 1.5 sin ( x ) .
The BVP has the exact solution u * ( x ) = x + x 1.5 + x 3 . The boundary parameters are p 0 = 0 , q 0 = 1 , b 0 = 1 , p 1 = 1 , q 1 = 0 , b 1 = 1 . The collocation equation system in (31) is
n = 0 m a n T n * ( x i ) I x , 1 1.5 x sin ( x ) 3 T n * ( x ) + sin ( x ) T n * ( x ) + I x , x i 1.5 x sin ( x ) 3 T n * ( x ) + sin ( x ) T n * ( x ) = h ( x i ) ,
where h ( x i ) = 2 x i I x , 1 1.5 g ( x ) + I x , x i 1.5 g ( x ) , i = 0 , 1 , , m .
Take m = 2 , 3 , 4 and 5, respectively, the solution approximations φ m ( x ) are calculated as
φ 2 ( x ) = 0.0175774 1.10039 x + 2.05832 x 2 , φ 3 ( x ) = 0.00586106 0.689478 x + 0.925577 x 2 + 0.768798 x 3 , φ 4 ( x ) = 0.00288415 0.76024 x + 1.24242 x 2 + 0.293418 x 3 + 0.22758 x 4 , φ 5 ( x ) = 0.00167028 0.8038 x + 1.54387 x 2 0.47843 x 3 + 1.05648 x 4 0.316564 x 5 .
The error function and maximum error of the approximate solution φ m ( x ) are defined as
E R m ( x ) = φ m ( x ) u * ( x ) and M E m = max 0 x 1 E R m ( x ) .
In Figure 1, the error functions E R m ( x ) for m = 2 , 3 , 4 , 5 are depicted, where at the m + 1 collocation points of φ m ( x ) , errors are zero. The maximum errors of the four approximate solutions are 0.028696, 0.005861, 0.002884, and 0.001670, respectively.
Example 2.
Consider the BVP for the linear FDE
D x λ u ( x ) u ( x ) = 4 x e x , 0 < x < 1 , 1 < λ 2 ,
u ( 0 ) u ( 0 ) = 1 , u ( 1 ) + u ( 1 ) = e .
If λ = 2 , the BVP has the exact solution u * ( x ) = x ( 1 x ) e x .
For this example, the coefficients and parameters are c 1 ( x ) = 0 , c 0 ( x ) = 1 , g ( x ) = 4 x e x ,   p 0 = q 0 = p 1 = q 1 = 1 , b 0 = 1 and b 1 = e . The collocation equation system in Equation (31) becomes
n = 0 m a n T n * ( x i ) x i + 1 3 I x , 1 λ T n * ( x ) x i + 1 3 I x , 1 λ 1 T n * ( x ) + I x , x i λ T n * ( x ) = h ( x i ) ,
where h ( x i ) = x i 2 e x i e 3 + x i + 1 3 I x , 1 λ g ( x ) I x , 1 λ 1 g ( x ) + I x , x i λ g ( x ) , i = 0 , 1 , , m .
For the case of λ = 2 , the error functions E R m ( x ) = φ m ( x ) u * ( x ) are depicted in Figure 2 for m = 2–5. The maximum errors of the approximate solutions are 0.069103, 0.007877, 0.000620, and 0.000038, respectively. For the case of λ = 1.5 , the solution approximations φ m ( x ) , m = 2–5, are calculated as
φ 2 ( x ) = 0.578503 + 3.00467 x 2.45143 x 2 , φ 3 ( x ) = 0.659109 + 1.56065 x + 1.43569 x 2 2.61289 x 3 , φ 4 ( x ) = 0.651626 + 1.81823 x + 0.108085 x 2 0.448696 x 3 1.09621 x 4 , φ 5 ( x ) = 0.653112 + 1.75369 x + 0.599869 x 2 1.78975 x 3 + 0.411041 x 4 0.596018 x 5 .
The condition numbers of the coefficient matrices W in the derivations of the four solution approximations are 2.85, 3.29, 3.67 and 4.02, respectively. These values show that the coefficient matrices W are well conditioned. We note that the condition number is based on the l 2 -matrix norm. The four solution approximations are plotted in Figure 3, where a fast convergence is shown.
Example 3.
Consider the BVP for the linear FDE
D x λ u ( x ) x 1 + x u ( x ) 1 1 + x u ( x ) = 0 , 0 < x < 1 , 1 < λ 2 ,
u ( 0 ) 2 u ( 0 ) = 1 , u ( 1 ) + 2 u ( 1 ) = 3 e .
If λ = 2 , the BVP has the exact solution u * ( x ) = e x .
The coefficients and parameters are c 1 ( x ) = x 1 + x , c 0 ( x ) = 1 1 + x , g ( x ) = 0 , p 0 = p 1 = 1 , q 0 = q 1 = 2 , b 0 = 1 , b 1 = 3 e . The collocation equation system in Equation (31) becomes
n = 0 m a n T n * ( x i ) x i + 2 5 I x , 1 λ x 1 + x T n * ( x ) + 1 1 + x T n * ( x ) 2 ( x i + 2 ) 5 I x , 1 λ 1 x 1 + x T n * ( x ) + 1 1 + x T n * ( x ) + I x , x i λ x 1 + x T n * ( x ) + 1 1 + x T n * ( x ) = h ( x i ) ,
where h ( x i ) = 1 5 ( 6 e 3 + 3 e x i + x i ) , i = 0 , 1 , , m .
For the case of λ = 2 , the error functions E R m ( x ) = φ m ( x ) u * ( x ) for m = 2 , 3 , 4 , 5 , are depicted in Figure 4. The maximum errors of the approximate solutions are 0.011605, 0.000742, 0.000037, and 0.000002, respectively. For the case of λ = 1.5 , the solution approximations φ m ( x ) for m = 2 , 3 , 4 , 5 are calculated as
φ 2 ( x ) = 0.627196 + 0.801952 x + 0.961626 x 2 , φ 3 ( x ) = 0.614413 + 0.94835 x + 0.55162 x 2 + 0.269951 x 3 , φ 4 ( x ) = 0.615168 + 0.908026 x + 0.734701 x 2 0.0100847 x 3 + 0.135331 x 4 , φ 5 ( x ) = 0.615706 + 0.894049 x + 0.828159 x 2 0.242103 x 3 + 0.378859 x 4 0.0913166 x 5 .
The condition numbers of the coefficient matrices W in the derivations of the four solution approximations are 4.87, 7.83, 11.90 and 17.05, respectively. So the coefficient matrices W are well conditioned. The four solution approximations are plotted in Figure 5.
In the three examples, fast convergent rates are shown only using the minor term number with M = 5 in Equation (38) for the integral computation of the known functions, and the minor term number with m = 2 , 3 , 4 and 5 in Equation (28) for the truncated Chebyshev series of the unknown function.

4. Conclusions

We considered the BVP for the varying coefficient linear Caputo-type fractional ODE subject to the mixed boundary conditions on the interval 0 x 1 . The BVP was conveniently converted into an equivalent differential–integral equation merging the boundary conditions. Then, the solution was decomposed into a truncated shifted Chebyshev series. The collocation method was used to determine the solution. In order to deal with the involved integrations, the varying coefficients were again decomposed into the truncated shifted Chebyshev series. Thus, the calculations of the integrals are only for polynomials and can be carried out exactly. Three numerical examples were solved by using the proposed method, where fast convergent rates are shown only using the minor term number with M = 5 in Equation (38) for the integral computation of the known functions, and the minor term number with m = 2 , 3 , 4 and 5 in Equation (28) for the truncated Chebyshev series of the unknown function.
In the presented method, there is no need to divide the interval commonly used in numerical methods. The collocation points or the zeros of the Chebyshev polynomials have exact explicit expressions. Approximate analytical solutions in the polynomial forms are obtained, which are different from a discrete numerical solution. The obtained approximate analytical solutions in the polynomial forms can be directly checked by substitution. The convergence and effectiveness of solutions can be examined by remainder errors. Convergence order of the approximate solutions could be further consideration in this field.

Author Contributions

Conceptualization, J.-S.D. and M.L.; data curation, L.-X.J. and M.L.; formal analysis, J.-S.D., L.-X.J. and M.L.; funding acquisition, J.-S.D. and M.L.; investigation, J.-S.D. and L.-X.J.; methodology, J.-S.D. and M.L.; software, J.-S.D. and L.-X.J.; supervision, J.-S.D.; validation, M.L.; visualization, L.-X.J.; writing–original draft, J.-S.D. and L.-X.J.; writing–review and editing, J.-S.D. and M.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Nos. 11772203; 61672238).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Jun-Sheng Duan acknowledges the supports in part by the National Natural Science Foundation of China under the project grant number 11772203. Ming Li acknowledges the supports in part by the National Natural Science Foundation of China under the project grant number 61672238. The authors show their appreciation for the valuable comments from the reviewers on the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Oldham, K.B.; Spanier, J. The Fractional Calculus; Academic: New York, NY, USA, 1974. [Google Scholar]
  2. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; Wiley: New York, NY, USA, 1993. [Google Scholar]
  3. Podlubny, I. Fractional Differential Equations; Academic: San Diego, CA, USA, 1999. [Google Scholar]
  4. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  5. Băleanu, D.; Diethelm, K.; Scalas, E.; Trujillo, J.J. Fractional Calculus Models and Numerical Methods. Series on Complexity, Nonlinearity and Chaos; World Scientific: Boston, MA, USA, 2012. [Google Scholar]
  6. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity; Imperial College: London, UK, 2010. [Google Scholar]
  7. Li, M. Theory of Fractional Engineering Vibrations; De Gruyter: Berlin, Germany, 2021. [Google Scholar]
  8. Bagley, R.L.; Torvik, P.J. On the fractional calculus model of viscoelastic behavior. J. Rheol. 1986, 30, 133–155. [Google Scholar] [CrossRef]
  9. Li, M. Three classes of fractional oscillators. Symmetry 2018, 10, 40. [Google Scholar] [CrossRef] [Green Version]
  10. Ferras, L.L.; Ford, N.J.; Morgado, M.L.; Rebelo, M.; McKinley, G.H.; Nobrega, J.M. Theoretical and numerical analysis of unsteady fractional viscoelastic flows in simple geometries. Comput. Fluids 2018, 174, 14–33. [Google Scholar] [CrossRef] [Green Version]
  11. Duan, J.S.; Hu, D.C.; Chen, Y.Q. Simultaneous characterization of relaxation, creep, dissipation, and hysteresis by fractional-order constitutive models. Fractal Fract. 2021, 5, 36. [Google Scholar] [CrossRef]
  12. Jiang, X.; Xu, M.; Qi, H. The fractional diffusion model with an absorption term and modified Fick’s law for non-local transport processes. Nonlinear Anal. Real World Appl. 2010, 11, 262–269. [Google Scholar] [CrossRef]
  13. Morgado, M.L.; Rebelo, M.; Ferras, L.L.; Ford, N.J. Numerical solution for diffusion equations with distributed order in time using a Chebyshev collocation method. Appl. Numer. Math. 2017, 114, 108–123. [Google Scholar] [CrossRef]
  14. Atangana, A. Fractional Operators with Constant and Variable Order with Application to Geo-Hydrology; Elsevier: London, UK, 2018. [Google Scholar]
  15. Jiao, Z.; Chen, Y.; Podlubny, I. Distributed-Order Dynamic Systems—Stability, Simulation, Applications and Perspectives; Springer: London, UK, 2012. [Google Scholar]
  16. Li, Y.; Duan, J.S. The periodic response of a fractional oscillator with a spring-pot and an inerter-pot. J. Mech. 2021, 37, 108–117. [Google Scholar] [CrossRef]
  17. Duan, J.S.; Hu, D.C. Vibration systems with fractional-order and distributed-order derivatives characterizing viscoinertia. Fractal Fract. 2021, 5, 67. [Google Scholar] [CrossRef]
  18. Diethelm, K. The Analysis of Fractional Differential Equations; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
  19. Duan, J.S.; Wang, Z.; Liu, Y.L.; Qiu, X. Eigenvalue problems for fractional ordinary differential equations. Chaos Solitons Fractals 2013, 46, 46–53. [Google Scholar] [CrossRef]
  20. Băleanu, D.; Mustafa, O.G.; Agarwal, R.P. An existence result for a superlinear fractional differential equation. Appl. Math. Lett. 2010, 23, 1129–1132. [Google Scholar] [CrossRef] [Green Version]
  21. Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52. [Google Scholar] [CrossRef] [Green Version]
  22. Liu, F.; Zhuang, P.; Anh, V.; Turner, I.; Burrage, K. Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation. Appl. Math. Comput. 2007, 191, 12–20. [Google Scholar] [CrossRef]
  23. Wu, G.C. A fractional characteristic method for solving fractional partial differential equations. Appl. Math. Lett. 2011, 24, 1046–1050. [Google Scholar] [CrossRef] [Green Version]
  24. Atabakzadeh, M.H.; Akrami, M.H.; Erjaee, G.H. Chebyshev operational matrix method for solving multi-order fractional ordinary differential equations. Appl. Math. Model. 2013, 37, 8903–8911. [Google Scholar] [CrossRef]
  25. Duan, J.S.; Hu, D.C.; Li, M. Comparison of two different analytical forms of response for fractional oscillation equation. Fractal Fract. 2021, 5, 188. [Google Scholar] [CrossRef]
  26. Jafari, H.; Daftardar-Gejji, V. Positive solutions of nonlinear fractional boundary value problems using Adomian decomposition method. Appl. Math. Comput. 2006, 180, 700–706. [Google Scholar] [CrossRef]
  27. Rehman, M.U.; Khan, R.A. A numerical method for solving boundary value problems for fractional differential equations. Appl. Math. Model. 2012, 36, 894–907. [Google Scholar] [CrossRef]
  28. Al-Refai, M.; Hajji, M.A. Monotone iterative sequences for nonlinear boundary value problems of fractional order. Nonlinear Anal. 2011, 74, 3531–3539. [Google Scholar] [CrossRef]
  29. Doha, E.H.; Bhrawy, A.H.; Ezz-Eldien, S.S. A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order. Comput. Math. Appl 2011, 62, 2364–2373. [Google Scholar] [CrossRef] [Green Version]
  30. Mason, J.C.; Handscomb, D.C. Chebyshev Polynomials; Chapman & Hall/CRC: London, UK, 2003. [Google Scholar]
Figure 1. The error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Figure 1. The error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Fractalfract 06 00148 g001
Figure 2. For λ = 2 , the error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Figure 2. For λ = 2 , the error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Fractalfract 06 00148 g002
Figure 3. For λ = 1.5 , the solution approximations φ m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Figure 3. For λ = 1.5 , the solution approximations φ m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Fractalfract 06 00148 g003
Figure 4. For λ = 2 , the error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Figure 4. For λ = 2 , the error functions E R m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Fractalfract 06 00148 g004
Figure 5. For λ = 1.5 , the solution approximations φ m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Figure 5. For λ = 1.5 , the solution approximations φ m ( x ) for m = 2 (solid line), m = 3 (dot line), m = 4 (dash line) and m = 5 (dot-dash line).
Fractalfract 06 00148 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Duan, J.-S.; Jing, L.-X.; Li, M. The Mixed Boundary Value Problems and Chebyshev Collocation Method for Caputo-Type Fractional Ordinary Differential Equations. Fractal Fract. 2022, 6, 148. https://doi.org/10.3390/fractalfract6030148

AMA Style

Duan J-S, Jing L-X, Li M. The Mixed Boundary Value Problems and Chebyshev Collocation Method for Caputo-Type Fractional Ordinary Differential Equations. Fractal and Fractional. 2022; 6(3):148. https://doi.org/10.3390/fractalfract6030148

Chicago/Turabian Style

Duan, Jun-Sheng, Li-Xia Jing, and Ming Li. 2022. "The Mixed Boundary Value Problems and Chebyshev Collocation Method for Caputo-Type Fractional Ordinary Differential Equations" Fractal and Fractional 6, no. 3: 148. https://doi.org/10.3390/fractalfract6030148

APA Style

Duan, J. -S., Jing, L. -X., & Li, M. (2022). The Mixed Boundary Value Problems and Chebyshev Collocation Method for Caputo-Type Fractional Ordinary Differential Equations. Fractal and Fractional, 6(3), 148. https://doi.org/10.3390/fractalfract6030148

Article Metrics

Back to TopTop