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 two-variable Chebyshev polynomials of the sixth kind, as basis functions in the proposed method, are constructed by the one-variable 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 one-dimensional fractional parabolic differential equations with delay. Hendy et al. [6] introduced a Crank–Nicolson difference approximation for solving multiterm time-fractional diffusion equations with delay. Nandal and Pandey constructed a linearized compact difference scheme for fourth-order 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 variable-order fractional differential equations. The Chebyshev polynomials of the first kind were used by Vlasic et al. [13] as basis functions to introduce a spline-like parametric model for compressive imaging. Nemati et al. [14] applied the second-kind 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 three-layer 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 liquid-loaded 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. Abd-Elhameed and Youssri presented a new numerical algorithm based on the sixth-kind Chebyshev polynomials for solving some linear and nonlinear fractional-order differential equations [18]. In [19], a new orthogonal wavelet based on the sixth-kind Chebyshev polynomials was constructed to obtain the solution of fractional optimal control problems. Abd-Elhameed used the sixth-kind Chebyshev polynomials for obtaining a numerical solution of nonlinear one-dimensional Burgers’ equations [20]. Atta et al. [21] employed shifted fifth-kind Chebyshev polynomials for the numerical solution of one-dimensional linear hyperbolic partial differential equations. A sixth-kind Chebyshev collocation method was considered in [22] for solving a class of variable-order fractional nonlinear quadratic integrodifferential equations. Bivariate Chebyshev polynomials of the fifth kind were utilized in [23] for variable-order time-fractional partial integrodifferential equations with the weakly singular kernel. Sadri and Aminikhah [24] employed fifth-kind Chebyshev polynomials for solving multiterm variable-order time-fractional diffusion-wave 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 time-fractional 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 above-mentioned 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 sixth-kind 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 two-variable 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 two-variable 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 Riemann-Liouville 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. Sixth-Kind 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 sixth-kind 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, two-variable Chebyshev polynomials of the sixth kind are constructed using them.
3.1. One-Variable Chebyshev Polynomials of the Sixth Kind
The following recurrence relation holds for the sixth-kind 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 square-integrable function , can be expanded in the shifted sixth-kind 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 One-Variable Basis
In this subsection, integral operational matrices of integral and fractional orders are obtained for the one-variable 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 sixth-kind Chebyshev polynomials in (12) and the weight function , one has where is the well-known 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 integer-order 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 sixth-kind 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 Riemann-Liouville 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. Two-Variable Chebyshev Polynomials of the Sixth Kind
Two-variable Chebyshev polynomials are constructed by one-variable 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 Two-Variable Basis
Consider the two-variable 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 Riemann-Liouville 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 , and is the approximate solution obtained from the method belonging to , and Then, the error bound of the can be obtained as
Proof. Define the bivariate Taylor expansion of as and Since is the best approximate solution of from and according to Equation (71), one gets By the Stirling formula [29], some bounds are computed for the last equality: where Therefore, the desired result is achieved.
Theorem 9. Assume that , , and , satisfy the condition of Theorem 8, and set Then,
Proof. The bivariate Taylor expansion of leads to Therefore, by taking -norm and using Equation (74) and Stirling formula, one has where
Theorem 10. Assume that and the conditions of Theorem 8 hold. Then,
Proof. According to Equation (71) and properties of the Caputo operator in (5), one can write Taking -norm yields
Corollary 11. Some error bounds for functions with proportional delays can be obtained using the resultant bounds in Theorems 8, 9, and 10.
Now, Theorems 8, 9, and 10 and Corollary 11 are applied to show that the error of the method becomes sufficiently small when is sufficiently large. For this, three equations (50), (62), and (65) in Section 4 are called again. Moreover, consider the following bounds:
Form I. Suppose is the best approximate solution from which is obtained from the proposed scheme. Therefore, it satisfies the following equation:
where is the error term. Subtracting Equation (83) from Equation (50) leads to the following equation:
So, one has
Form II. If is the solution obtained from the proposed algorithm for Equation (62), then it satisfies the following equation:
Subtracting Equation (86) from Equation (62) leads to
Then, one has
Form III. If is an approximate solution obtained from the suggested algorithm for Equation (65), then one has
Subtracting (89) from Equation (65) leads to
By taking -norm, one gets
As seen from the right-hand sides of inequalities (85), (88), and (91), will decrease by choosing appropriate values of , i.e., the error bound will be sufficiently small for the sufficiently large values of .
5. Numerical Examples
To demonstrate the accuracy and validity of the proposed method, three given examples in Refs. [26, 27] are solved in this section. These equations have been solved using the homotopy perturbation and natural transformation decomposition methods in [26, 27], respectively. The approximate solutions are compared to the exact ones and those reported by [26, 27], maximum absolute errors are calculated, and results are reported in tables and figures.
Example 1. Consider Equation (50) with the exact solution for and its corresponding residual function Equation (61) that is collocated at roots of the Chebyshev polynomials of the sixth kind. Therefore, an approximate solution is obtained. The surface plots of approximate solutions are depicted in Figure 1, and the absolute error function is seen in Figure 2(a) for and . Also, the figures of approximate solutions at are seen in Figure 2(b) for and . The values of the approximate solution at the selected points are listed in Table 1 which compared to the values of the exact one for and various values of . The results are compared to those reported in [26, 27] in Table 2. As seen, the proposed method presents better accuracy. Besides, the maximum absolute errors of approximate solutions are computed for different values of and , and results are observed in Table 3. Increasing leads to the decrease of the values of the errors.
Example 2. Consider the nonlinear fractional partial differential equation with the proportional delay in (62) and its corresponding residual function in Equation (64). The exact solution is if . The 3D figures of the approximate solutions are depicted in Figure 3(a) for and . The plot of the absolute error function is seen in Figure 3(b). The approximate solutions are plotted at and compared to those presented by [26] in Figure 4 for and . As seen, our obtained solutions converge faster to the exact one. The values of the resultant solution at the selected points are listed in Table 4 which compared to the values of the exact solution for and various values of . The obtained results are compared to those reported in [26, 27] in Table 5. As seen, the proposed method presents better accuracy. The maximum absolute errors are seen in Table 6 for and .
Example 3. Consider Form III in (65) and residual function (67). The exact solution is if . The three-dimensional figures of the approximate solutions are depicted in Figures 5(a)–5(d) for and . The plot of the absolute error function is seen in Figure 5(e) for and . The obtained solutions are plotted in Figure 6 at for and which are compared to those presented by [27]. The values of the resultant solution at the selected points are listed in Table 7 which are compared to the values of the exact one for and various values of . As , the numerical solutions converge to the exact one. The obtained results are compared to those reported in [26, 27] in Table 8 for and . As seen, the proposed method presents better accuracy.
(a)
(b)
(c)
(d)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
6. Conclusion
This paper deals with numerically solving a class of fractional partial differential equations with proportional delays on the domain . A spectral collocation approach, based on the sixth-kind Chebyshev polynomials as basis functions, has been considered to solve this class of equations. The two-variable Chebyshev polynomials of the sixth kind were introduced, and their integral operational matrices were derived. The relationship between the delay Chebyshev polynomials and the original basis was stated in a matrix form called delay operational matrix. The numerical results were reported in tables and figures and confirmed the accuracy and good agreement of the approximate solutions with exact ones. An error analysis has been presented which showed that the method error becomes small when is properly selected. The picked examples were also solved by homotopy perturbation and natural decomposition methods in [26, 27], and values of approximate solutions were reported at some selected points. It was clear from Tables 2, 5, and 8; the proposed method is more efficient. Therefore, the sixth-kind Chebyshev polynomials can be used to numerically solve other fractional functional equations.
Data Availability
All results have been obtained by conducting the numerical procedure, and the ideas can be shared for the researchers.
Conflicts of Interest
The authors declare that they have no conflicts of interest.