Abstract
This work devotes to solving a class of delay fractional partial differential equations that arises in physical, biological, medical, and climate models. For this, a numerical scheme is implemented that applies operational matrices to convert the main problem into a system of algebraic equations; then, solving the resultant system leads to an approximate solution. The twovariable Chebyshev polynomials of the sixth kind, as basis functions in the proposed method, are constructed by the onevariable ones, and their operational matrices are derived. Error bounds of approximate solutions and their fractional and classical derivatives are computed. With the aid of these bounds, a bound for the residual function is estimated. Three illustrative examples demonstrate the simplicity and efficiency of the proposed method.
1. Introduction
Mathematical modeling of some physical and biological phenomena leads to delay fractional differential equations (DFDEs) [1–3]. The independent variables and represent time and position in space or size of cells, and so on. The solutions can stand for temperature densities of cells, chemicals, etc. Hardly obtaining exact solutions to these equations necessitates mathematicians to construct some vigorous numerical and semianalytical schemes to handle solving these problems. Nevertheless, few methods exist for solving delay partial differential equations. The interest of scientists and mathematicians in DFDEs has resulted in the presentation of efficient schemes to solve this category of equations. For example, Pimenov and Hendy presented a difference scheme for a class of fractional diffusion equations with fixed time delay [4]. A compact difference scheme was constructed in [5] for the numerical solution of onedimensional fractional parabolic differential equations with delay. Hendy et al. [6] introduced a Crank–Nicolson difference approximation for solving multiterm timefractional diffusion equations with delay. Nandal and Pandey constructed a linearized compact difference scheme for fourthorder nonlinear fractional subdiffusion with time delay [7].
One of the most popular methods for solving diverse functional equations is the spectral method. The nature of the spectral methods has been joined with the orthogonal polynomials and functions. Orthogonal polynomials are utilized as basis functions in many numerical methods; hence, addressing the properties of orthogonal polynomials is important. For example, some properties of the generalized Gegenbauer polynomials were studied in [8, 9]. Bracciali et al. dealt with a class of Sobolev orthogonal polynomials and Hahn polynomials on the unit circle in [10]. Asymptotic approximations of Jacobi polynomials and their zeros were given in [11]. The shifted Chebyshev polynomials of the third kind were proposed in [12] to solve multiterm variableorder fractional differential equations. The Chebyshev polynomials of the first kind were used by Vlasic et al. [13] as basis functions to introduce a splinelike parametric model for compressive imaging. Nemati et al. [14] applied the secondkind Chebyshev polynomials for fractional integrodifferential equations with weakly singular kernels. Dahmen and Glorieux applied an extension of the Legendre polynomial method to model coupled Lamb wave parameters for defect detection in anisotropic composite threelayer with Kelvin–Voigt viscoelasticity [15]. A Legendre orthogonal polynomial method was proposed to calculate the reflection and transmission coefficients of plane wave at the liquid interface of a liquidloaded functionally gradient material plate [16]. Masjed–Jamei [17] presented two new classes of orthogonal polynomials which are called Chebyshev polynomials of fifth and sixth kinds. AbdElhameed and Youssri presented a new numerical algorithm based on the sixthkind Chebyshev polynomials for solving some linear and nonlinear fractionalorder differential equations [18]. In [19], a new orthogonal wavelet based on the sixthkind Chebyshev polynomials was constructed to obtain the solution of fractional optimal control problems. AbdElhameed used the sixthkind Chebyshev polynomials for obtaining a numerical solution of nonlinear onedimensional Burgers’ equations [20]. Atta et al. [21] employed shifted fifthkind Chebyshev polynomials for the numerical solution of onedimensional linear hyperbolic partial differential equations. A sixthkind Chebyshev collocation method was considered in [22] for solving a class of variableorder fractional nonlinear quadratic integrodifferential equations. Bivariate Chebyshev polynomials of the fifth kind were utilized in [23] for variableorder timefractional partial integrodifferential equations with the weakly singular kernel. Sadri and Aminikhah [24] employed fifthkind Chebyshev polynomials for solving multiterm variableorder timefractional diffusionwave equations.
Recently, spectral methods coupled with operational matrices have attracted the attention of many mathematicians and researchers. The advantage of applying operational matrices is to express the derivatives of the orthogonal polynomials as basis functions in terms of the linear combinations of original polynomials and rewrite these combinations as a sparse matrix form which decreases the computational costs [14, 20, 21, 24, 25]. In the current work, a class of timefractional partial differential equations with the proportional delay as the following form is considered [26, 27]:
with the conditions
where , is the Caputo operator, , , and and are linear and nonlinear differential operators, respectively. In [26, 27], the homotopy perturbation and natural decomposition methods have been applied for solving problems (1) and (2). The two abovementioned methods provided approximate solutions based on the Taylor expansions of time parts of the solutions which have only good accuracy for the classical case [26, 27]. The goal of the present paper is to construct a scheme using the sixthkind Chebyshev polynomials; hence, integral operational matrices of integer and fractional orders are derived. Moreover, an operational matrix is constructed to show the relation between the original basis and its delay form. Then, obtained matrices are utilized to obtain corresponding operational matrices for the twovariable basis. Resultant matrices accompanying the collocation method convert the main problem (1) and (2) into a system of algebraic equations, the solving of which leads to an approximate solution. It is worth noting that the obtained nonlinear algebraic system can be solved using Newton’s iteration method.
The rest of the paper is structured as follows: Section 2 recalls some basic definitions of fractional calculus and its properties. The one and twovariable Chebyshev polynomials of the sixth kind are introduced, and their operational matrices are constructed in Section 3. The idea of the proposed method is described and the error analysis is presented in Section 4. The accuracy and efficiency of the scheme are successfully demonstrated by implementing the algorithm on three examples in Section 5. Finally, a conclusion is given in Section 6.
2. Preliminaries
In this section, some definitions that are useful throughout the paper are presented.
Definition 1. A real function belongs to the space if a real number exists such that where , and it belongs to the space if and only if [28].
Definition 2. Suppose that and The RiemannLiouville fractional integral of the order is defined as [28]
Definition 3. The Caputo fractional derivative of the order of the funcion and is given by the following expression [28]: The following properties of these operators hold:
3. SixthKind Chebyshev Polynomials and Their Operational Matrices
The family of Chebyshev polynomials has found popularity in different spectral and pseudospectral methods [12–14, 18, 19, 21]. A class of the Chebyshev polynomials, called sixthkind Chebyshev polynomials, was proposed for the first time in [18] to solve fractional ordinary differential equations. In this section, first, the shifted form of these polynomials is introduced over ; then, twovariable Chebyshev polynomials of the sixth kind are constructed using them.
3.1. OneVariable Chebyshev Polynomials of the Sixth Kind
The following recurrence relation holds for the sixthkind Chebyshev polynomials
where
These polynomials are orthogonal with respect to the weight function , that is,
where
By the change of variable , the shifted Chebyshev polynomials are orthogonal regarding the weight function on the interval ,
and
The series form of the shifted Chebyshev polynomials of the sixth kind is as follows:
where
Every squareintegrable function , can be expanded in the shifted sixthkind Chebyshev polynomials as
where the coefficients are computed as
The first few coefficients in (14) practically keep information of function . In other words, a finite series can present an approximation to as
where and are the vectors as follows:
3.2. Operational Matrices of OneVariable Basis
In this subsection, integral operational matrices of integral and fractional orders are obtained for the onevariable basis. Furthermore, the relation between the main basis and its delay form is given as a matrix. For this, some useful lemma and theorems are stated and proved.
Lemma 4. If , then
Proof. By the series form of the shifted sixthkind Chebyshev polynomials in (12) and the weight function , one has where is the wellknown beta function, so, the desired result is achieved.
Theorem 5. If is the basis vector in (17), the integral of can be computed as where is the integral operational matrix of the integerorder in the following form: where the entries are computed as
Proof. Integrating the elements of the vector yields Now, is approximated in terms of the shifted sixthkind Chebyshev polynomials where Using Lemma 4, the integral part of (25) is computed as follows: Therefore, (23) is written as By rewriting the last series as a matrix form, the desired result is achieved.
Theorem 6. Assume that is the basis vector in (17) and is the RiemannLiouville integral operator of the order , . Then, one has where is the fractional operational matrix of the order as follows: where the entries are computed as
Proof. The proof process is similar to Theorem 5. Noting the definition of and its properties in (5), the fractional integral of is computed as Now, is approximated by the Chebyshev polynomials of the sixth kind as By Lemma 4 and pursuing the proof process in Theorem 5, Equation (31) is written as
Theorem 7. Assume that is the basis vector and , is its delay form. can be approximated in as where is a matrix as follows: and is computed from the following recurrence formulas: with the starting values and .
Proof. Consider the following recurrence formula obtained from Equation (6) Also, the following auxiliary relation is obtained from formula (6): Now, can be expanded in terms of the Chebyshev polynomials of the sixth kind as follows: From the first few polynomials in Equation (39), it is easily obtained , . Using auxiliary relation (38) and substituting Equation (39) into Equation (37), one gets Equating coefficients of on both sides of the last equality leads to recurrence formula (36).
3.3. TwoVariable Chebyshev Polynomials of the Sixth Kind
Twovariable Chebyshev polynomials are constructed by onevariable ones on the domain as
These polynomials are orthogonal regarding the weight function on ,
where and are calculated by (11). The function is expanded as
and a truncated series of (43) is considered as an approximation to the function ,
where such that and are the vectors as
3.4. Operational Matrices of TwoVariable Basis
Consider the twovariable basis in (45). The integral operational matrices of with respect to variables and are obtained, respectively, as
where and are the integral operational matrices related to and , respectively, is the operational matrix in Theorem 5, and is the identity matrix. Similarly, the fractional integral of of the order with respect to can be computed as
where is the fractional integral operational matrix of the order related to , is the operational matrix in Theorem 6, and is the identity matrix. Now, the relationship of the vectors , , and to the basis vector is specified.
By setting in , the vector is rewritten as follows:
where is the matrix in Theorem 7, is the zero matrix, and is the delay matrix. Similarly, it is found that
where and are matrices.
3.5. Solution Method
To describe the methodology, three forms of Equations (1) and (2) are considered [26, 27]:
Form I:
According to the highest orders of derivatives regarding and , the following approximation is considered:
Integrating (51) concerning and , respectively, leads to the following approximations:
Now, twice integrating approximation (51) with respect to leads to an approximation for :
To obtain an approximation to , approximation (55) is rewritten as follows:
By applying the RiemannLiouville operator of the order to both sides of (56), one gets
Fractionally differentiating approximation (54) and setting lead to
Approximation (58) equals zero because after fractionally differentiating related to , all components of the basis vector involve terms as or are zero.
The terms with delays can be approximated as
Substituting approximations (52)–(60) into Equation (50) results in the following residual function:
Form II:
The functions in Equation (62) are approximated based on what was done for Equation (50):
Substituting approximations (63) into Equation (62) leads to the following residual function:
Form III:
The following approximations can be obtained for the functions in Equation (65):
By substituting approximations (66) into Equation (65), one gets the following residual function:
The residual functions in Equations (61), (64), and (65) are collocated at collocation nodes , which and are roots of and , respectively. Hence, a nonlinear system involving algebraic equations is achieved that can be solved by the Newton’s iteration method. Therefore, the coefficient vector is determined approximately; then, an approximate solution, , is achieved. To better describe the solution method, pursuing Algorithm 1 for Equation (50) is suggested.

4. Error Analysis
It is found that the value of the error function decreases when values of increase. First, some error bounds are obtained for the unknown function and its derivatives.
Theorem 8. Assume that