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) [13]. 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 [1214, 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


These polynomials are orthogonal with respect to the weight function , that is,


By the change of variable , the shifted Chebyshev polynomials are orthogonal regarding the weight function on the interval ,


The series form of the shifted Chebyshev polynomials of the sixth kind is as follows:


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.

Step 1. Derive operational matrices , , and from (21), (29), and (36).
Step 2. Construct the operational matrices , , , , , from (46)-(49).
Step 3. Consider the approximation in (51).
Step 4. Find approximations to , , , and from (52)-(55).
Step 5. Find an approximation to from (57).
Step 6. Determine the residual function from (61).
Step 7. Obtain roots of and () using command in Maple.
Step 8. Collocate the residual function at tensor points .
Step 9. Solve the resultant non-linear system in Step 8 by Newton’s iteration method and obtain the unknown vector .
Step 10. Find from (54).

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