# Chebyshev Polynomials of Sixth Kind for Solving Nonlinear Fractional PDEs with Proportional Delay and Its Convergence Analysis

Authors: Khadijeh Sadri; Hossein Aminikhah
Journal: Journal of Function Spaces
 (2022)

---

## 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.


---

## Body

## 1. Introduction

Mathematical modeling of some physical and biological phenomena leads to delay fractional differential equations (DFDEs) [1–3]. The independent variables t and x 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]:
(1)D0Ctσux,t+Lux,t+Nuq1x,q2t=Qx,t,x,t∈Ω,with the conditions(2)ux,0=hx,u0,t=g1t,uxx,t=g2t,whereDtσ0C, 0<σ≤1 is the Caputo operator, Ω=0,1×0,1, 0<qi≤1,i=1,2, and L and N 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 σ=1 [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: Section2 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 functionft,t>0 belongs to the space Cq,q∈ℝ if a real number p>q exists such that ft=tpf1t where f1t∈0,∞, and it belongs to the space Cqn,n∈ℕ if and only if fnt∈Cq [28].Definition 2.
Suppose thatft∈Cq,t>0 and q>−1. The Riemann-Liouville fractional integral of the order σ>0 is defined as [28]
(3)I0RLtσft=1Γσ∫0tt−sσ−1fsds,σ>0,ft,σ=0.Definition 3.
The Caputo fractional derivative of the orderσ>0 of the funcion ft∈Cq,t>0 and q>−1 is given by the following expression [28]:
(4)D0Ctσft=1Γm−σ∫0tt−sm−σ−1dmfsdsmds,m−1<σ≤m,ft,σ=0.
The following properties of these operators hold:(5)Dtσ0Ctm=Γm+1tm−σΓm−σ+1,m≥σ,0,m<σ,0RLItσtm=Γm+1tm+σΓm+σ+1,0CDtσ10CDtσ2ft=0CDtσ1+σ2ft,0CDtσ0RLItσft=ft,0RLItσ0CDtσft=ft−f0,0<σ<1.


## 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 0,1; 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(6)Y¯jz=zY¯j−1z−θj+1Y¯j−2z,j≥2,z∈−1,1,Y¯0z=1,Y¯1z=z,where(7)θj=jj+1+−1j2j+1+14jj+1.These polynomials are orthogonal with respect to the weight functionw¯z=z21−z2, that is,
(8)∫−11Y¯izY¯jzz21−z2dz=ℏ¯iδij,where(9)ℏ¯i=π22i+3,ieven,πi+322i+3i+1,iodd.By the change of variablez=2t−1, the shifted Chebyshev polynomials Yit=Y¯i2t−1 are orthogonal regarding the weight function wt=2t−12t−t2 on the interval 0,1,
(10)∫01YitYjt2t−12t−t2dt=ℏiδij,and(11)ℏ¯i=π22i+5,ieven,πi+322i+5i+1,iodd.The series form of the shifted Chebyshev polynomials of the sixth kind is as follows:(12)Yjt=∑r=0jςr,jtr,where(13)ςr,j=22r−j2r+1!∑l=r+1/2j/2−1j/2+l+r2l+r+1!2l−r+1!,jeven,2j+1∑l=r/2j−1/2−1j+1/2+l+rl+12l+r+2!2l−r+1!,jodd.Every square-integrable functionvt∈Lw2I, I=0,1 can be expanded in the shifted sixth-kind Chebyshev polynomials as
(14)vt=∑j=0∞VjYjt,t∈I,where the coefficientsVj are computed as
(15)Vj=1ℏj∫01vtYjtwtdt.The first few coefficients in (14) practically keep information of function vt. In other words, a finite series can present an approximation to vt as
(16)vt≈vNt=∑j=0NVjYjt=YTtV=VTYt,whereV and Yt are the N+1×1 vectors as follows:
(17)V=V0V1⋯VNT,Yt=Y0tY1t⋯YNtT.


### 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
(18)∫01tνYktwtdt=∑m=0kςm,kΓ324Γν+m+7/2Γν+m+5−4Γν+m+5/2Γν+m+4+Γν+m+3/2Γν+m+3.Proof.
By the series form of the shifted sixth-kind Chebyshev polynomials in (12) and the weight function wt, one has
(19)∫01tνYktwtdt=∫01tν∑m=0kςm,ktm2t−12t−t2dt=∑m=0kςm,k∫014tν+m+5/2−4tν+m+3/2+xν+m+1/21−t1/2dt=∑m=0kςm,k4Bν+m+72,32−4Bν+m+52,32+Bν+m+32,32,
whereBr,s is the well-known beta function, so, the desired result is achieved.Theorem 5.
IfYt is the basis vector in (17), the integral of Yt can be computed as
(20)∫0tYsds≈PYt,t∈I,
whereP is the N+1×N+1 integral operational matrix of the integer-order in the following form:
(21)P=π0,0π0,1⋯π0,Nπ1,0π1,1⋯π1,N⋮⋮⋱⋮πj,0πj,1⋯πj,N⋮⋮⋱⋮πN,0πN,1⋯πN,N,
where the entriesπj,k are computed as
(22)πj,k=∑r=0jςr,jΓ3/2r+1ℏk∑m=0kςm,k4Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4,j=0,1,⋯,N,k=0,1,⋯,N.Proof.
Integrating the elements of the vectorYt yields
(23)∫0tYjsds=∑r=0jςr,j∫0tsrds=∑r=0jςr,jtr+1r+1,j=0,1,⋯,N.
Now,tr+1 is approximated in terms of the shifted sixth-kind Chebyshev polynomials
(24)tr+1≈∑k=0Nρk,r+1Ykt,
where(25)ρk,r+1=1ℏk∫01tr+1Yktwtdt.
Using Lemma4, the integral part of (25) is computed as follows:
(26)∫01tr+1Yktwtdt=∑m=0kςm,kΓ324Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4.
Therefore, (23) is written as
(27)∫0tYjsds≈∑k=0N∑r=0jςr,jΓ3/2r+1ℏk∑m=0kςm,k4Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4Ykt.
By rewriting the last series as a matrix form, the desired result is achieved.Theorem 6.
Assume thatYt is the basis vector in (17) and I0RLtσ is the Riemann-Liouville integral operator of the order σ, 0<σ<1. Then, one has
(28)I0RLtσYt≈PσYt,
wherePσ is the N+1×N+1 fractional operational matrix of the order σ as follows:
(29)Pσ=πσ0,0πσ0,1⋯πσ0,Nπσ1,0πσ1,1⋯πσ1,N⋮⋮⋱⋮πσi,0πσi,1⋯πσi,N⋮⋮⋱⋮πσN,0πσN,1⋯πσN,N,
where the entriesπσi,k are computed as
(30)πσi,k=∑r=0iςr,iΓr+σ+1Γ3/2Γr+σ+2ℏk∑m=0kςm,k4Γr+m+7/2Γr+m+5−4Γr+m+5/2Γr+m+4+Γr+m+3/2Γr+m+3,i=0,1,⋯,N,k=0,1,⋯,N.Proof.
The proof process is similar to Theorem5. Noting the definition of 0RLItσ and its properties in (5), the fractional integral of Yit is computed as
(31)I0RLtσYit=∑r=0iςr,iΓr+σ+1tr+σΓr+σ+2,i=0,1,⋯,N.
Now,tr+σ is approximated by the Chebyshev polynomials of the sixth kind as
(32)tr+σ≈∑k=0Nρk,rYkt,s.t.ρk,r=1ℏk∫01tr+σYktwtdt.
By Lemma4 and pursuing the proof process in Theorem 5, Equation (31) is written as
(33)I0RLtσYit≈∑k=0N∑r=0iςr,iΓr+σ+1Γ3/2Γr+σ+2ℏk×∑m=0kςm,k4Γr+σ+m+7/2Γr+σ+m+5−4Γr+σ+m+5/2Γr+σ+m+4+Γr+σ+m+3/2Γr+σ+m+3Ykt,i=0,1,⋯,N.Theorem 7.
Assume thatYt is the basis vector and Yqt, 0<q<1 is its delay form. Yqt can be approximated in Yt as
(34)Yqt≈LYt,
whereL is a N+1×N+1 matrix as follows:
(35)L=100⋯0l01q0⋯0l02l12q⋯0⋮⋮⋮⋱⋮l0Nl1Nl2N⋯qN,
andlkj is computed from the following recurrence formulas:
(36)l0j+1=q−1l0j+qθ2l1j−θjl0j−1,j=1,2,⋯,N−1,lkj+1=qlk−1j+q−1lkj+qθk+1lk+1j−θjlkj−1,k=1,2,⋯,j−1,ljj+1=qlj−1j+q−1qj,j=1,2,⋯,N−1,
with the starting valuesljj=qj,j=0,1,⋯,N, and l01=q−1.Proof.
Consider the following recurrence formula obtained from Equation (6)
(37)Yj+1qt=2qt−1Yjqt−θj+1Yj−1qt,j≥1.
Also, the following auxiliary relation is obtained from formula (6):
(38)2tYjt=Yj+1t+Yjt+θj+1Yj−1t,j≥1.
Now,Yjqt can be expanded in terms of the Chebyshev polynomials of the sixth kind as follows:
(39)Yjqt=∑k=0jlkjYkt,j=0,1,⋯,N.
From the first few polynomials in Equation (39), it is easily obtained l01=q−1, ljj=qj,j=0,1,⋯,N. Using auxiliary relation (38) and substituting Equation (39) into Equation (37), one gets
(40)∑k=0j+1lkj+1Ykt=q∑k=0jlkjYk+1t+Ykt+θk+1Yk−1t−θj+1∑k=0j−1lkj−1Ykt−∑k=0jlkjYkt=∑k=0jlkjqYk+1t+q−1Ykt+qθk+1Yk−1t−θj+1∑k=0j−1lkj−1Ykt.
Equating coefficients ofYkt 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Ω=0,1×0,1 as
(41)Wijx,t=YixYjt,i,j=0,1,⋯,x,t∈Ω.These polynomials are orthogonal regarding the weight functionωx,t=wxwt on Ω,
(42)∫01∫01Wijx,tWklx,tωx,tdxdt=ℏiℏjδikδjl,whereℏi and ℏj are calculated by (11). The function Vx,t∈Lω2Ω is expanded as
(43)Vx,t=∑i=0∞∑j=0∞VijWijx,t,x,t∈Ω,and a truncated series of (43) is considered as an approximation to the function Vx,t,
(44)Vx,t≈VNx,t=∑i=0N∑j=0NVijWijx,t=∑i=0N+12−1Vi∗Wi∗x,t=VTWx,t=WTx,tV,whereVi∗=Vrs,Wi∗x,t=Wrsx,t such that r=i/N+1,s=i−rN+1, and V,Wx,t are the N+12×1 vectors as
(45)V=V00V01⋯V0NV10V11⋯V1N⋯VN0VN1⋯VNNT,Wx,t=W00x,tW01x,t⋯W0Nx,tW10W11x,t⋯W1Nx,t⋯WN0x,tWN1x,t⋯WNNx,tT.


### 3.4. Operational Matrices of Two-Variable Basis

Consider the two-variable basisWx,t in (45). The integral operational matrices of Wx,t with respect to variables x and t are obtained, respectively, as
(46)∫0xWs,tds≈ℙx1Wx,t=P⊗IWx,t,∫0tWx,τdτ≈ℙt1Wx,t=I⊗PWx,t,whereℙx1 and ℙt1 are the N+12×N+12 integral operational matrices related to x and t, respectively, P is the operational matrix in Theorem 5, and I is the N+1×N+1 identity matrix. Similarly, the fractional integral of Wx,t of the order σ with respect to t can be computed as
(47)I0RLtσWx,t≈ℙtσWx,t=I⊗PσWx,t,whereℙtσ is the N+12×N+12 fractional integral operational matrix of the order σ related to t, Pσ is the operational matrix in Theorem 6, and I is the N+1×N+1 identity matrix. Now, the relationship of the vectors Wx,qt, Wqx,t, and Wqx,qt to the basis vector Wx,t is specified.By settingt=qt in Wx,t, the vector Wx,qt is rewritten as follows:
(48)Wx,qt=Y0xY0qt,Y0xY1qt,⋯,Y0xYNqt,Y1xY0qt,Y1xY1qt,⋯,Y1xYNqt,⋯,YNxY0qt,YNxY1qt,⋯,YNxYNqtT=Y0xY0qt,Y1qt,⋯,YNqt,Y1xY0qt,Y1qt,⋯,YNqt,⋯,YNxY0qt,Y1qt,⋯,YNqtT=Y0xYqt,Y1xYqt,⋯,YNxYqtT≈Y0xLYt,Y1xLYt,⋯,YNxLYtT=LON+1⋯ON+1ON+1L⋯ON+1⋮⋮⋱⋮ON+1ON+1⋯LY0xY0t⋮YNxY0t⋮YNxYNt=I⊗LWx,t=L∗∗Wx,t,whereL is the matrix in Theorem 7, ON+1 is the N+1×N+1 zero matrix, and L∗∗ is the N+12×N+12 delay matrix. Similarly, it is found that
(49)Wqx,t≈L⊗IWx,t=L∗Wx,t,Wqx,qt≈L∗L∗∗Wx,t=LWx,t,L=L∗L∗∗,whereL∗ and L are N+12×N+12 matrices.


### 3.5. Solution Method

To describe the methodology, three forms of Equations (1) and (2) are considered [26, 27]:Form I:(50)D0Ctσux,t−∂2ux,t∂x2−∂ux,t/2∂xux2,t2−12ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x,u0,t=0,ux0,t=expt.According to the highest orders of derivatives regardingx and t, the following approximation is considered:
(51)∂3ux,t∂x2∂t≈WTx,tC.Integrating (51) concerning t and x, respectively, leads to the following approximations:
(52)∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC,(53)∂ux,t∂x≈WTx,tℙx1Tℙt1TC+∂u0,t∂x≈WTx,tℙx1Tℙt1TC+WTx,tF,(54)ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1TF+u0,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1TF.Now, twice integrating approximation (51) with respect to x leads to an approximation for ∂ux,t/∂t:
(55)∂2ux,t∂x∂t≈WTx,tℙx1TC+WTx,tF,∂ux,t∂t≈WTx,tℙx1T2C+WTx,tℙx1TF.To obtain an approximation to0CDtσux,t, approximation (55) is rewritten as follows:
(56)∂ux,t∂t=∂1−σ∂σux,t∂t1−σ∂tσ≈WTx,tℙx1T2C+WTx,tℙx1TF.By applying the Riemann-Liouville operator of the order1−σ to both sides of (56), one gets
(57)D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C+WTx,tℙt1−σTℙx1TF+0CDtσux,0.Fractionally differentiating approximation (54) and setting t=0 lead to
(58)D0Ctσux,0≈0CDtσWTx,tℙx1T2ℙt1TC+0CDtσWTx,tℙx1TFt=0=0.Approximation (58) equals zero because after fractionally differentiating Wx,t related to t, all components of the basis vector involve terms as tμ,1−σ<μ<N−σ or are zero.The terms with delays can be approximated as(59)∂ux,t/2∂x≈WTx,t2ℙx1Tℙt1TC+WTx,t2F≈WTx,tL∗∗Tℙx1Tℙt1TC+WTx,tL∗∗TF,(60)ux2,t2≈WTx2,t2ℙx1T2ℙt1TC+WTx2,t2ℙx1TF≈WTx,tLTℙx1T2ℙt1TC+WTx,tLTℙx1TF.Substituting approximations (52)–(60) into Equation (50) results in the following residual function:
(61)RNx,t=WTx,tℙt1−σTℙx1T2C+WTx,tℙt1−σTℙx1TF−WTx,tℙt1TC−WTx,tL∗∗Tℙx1Tℙt1TC+WTx,tL∗∗TFWTx,tLTℙx1T2ℙt1TC+WTx,tLTℙx1TF−12WTx,tℙx1T2C+WTx,tℙx1TF≈0.Form II:(62)D0Ctσux,t−∂2ux,t/2∂x2ux,t2+ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x2,u0,t=0,ux0,t=0.The functions in Equation (62) are approximated based on what was done for Equation (50):
(63)∂3ux,t∂x2∂t≈WTx,tC,∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC+WTx,tF,ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F,D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C,∂2ux,t/2∂x2≈WTx,tL∗∗Tℙt1TC+WTx,tL∗∗TF,ux,t2≈WTx,tL∗∗Tℙx1T2ℙt1TC+WTx,tL∗∗Tℙx1TF.Substituting approximations (63) into Equation (62) leads to the following residual function:
(64)RNx,t=WTx,tℙt1−σTℙx1T2C−WTx,tL∗∗Tℙt1TC+WTx,tL∗∗TF×WTx,tL∗∗Tℙx1T2ℙt1TC+WTx,tL∗∗Tℙx1TF+WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F≈0.Form III:(65)D0Ctσux,t−∂2ux/2,t/2∂x2∂ux/2,t/2∂x+∂ux,t∂x+ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x2,u0,t=0,ux0,t=0.The following approximations can be obtained for the functions in Equation (65):
(66)∂3ux,t∂x2∂t≈WTx,tC,∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC+WTx,tF,∂ux,t∂x≈WTx,tℙx1Tℙt1TC+WTx,tℙx1TF,ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F,D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C,∂2ux/2,t/2∂x2≈WTx,tLTℙt1TC+WTx,tLTF,∂ux/2,t/2∂x≈WTx,tLTℙx1Tℙt1TC+WTx,tLTℙx1TF.By substituting approximations (66) into Equation (65), one gets the following residual function:
(67)RNx,t=WTx,tℙt1−σTℙx1T2C−WTx,tLTℙt1TC+WTx,tLTF×WTx,tLTℙx1Tℙt1TC+WTx,tLTℙx1TF+WTx,tℙx1Tℙt1TC+WTx,tℙx1TF+WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F≈0.The residual functions in Equations (61), (64), and (65) are collocated at collocation nodes xi,tji,j=0∞, which xi and tj are roots of YN+1x and YN+1t, respectively. Hence, a nonlinear system involving N+12 algebraic equations is achieved that can be solved by the Newton’s iteration method. Therefore, the coefficient vector C is determined approximately; then, an approximate solution, ux,t, is achieved. To better describe the solution method, pursuing Algorithm 1 for Equation (50) is suggested.Algorithm 1: Solution method
Input:σ,NStep 1. Derive operational matrices P, Pσ, and L from (21), (29), and (36).Step 2. Construct the operational matrices ℙx1, ℙt1, ℙtσ, L, L∗, L∗∗ from (46)-(49).Step 3. Consider the approximation ∂3ux,t/∂x2∂t≈WTx,tC in (51).Step 4. Find approximations to ∂2ux,t/∂x2, ∂ux,t/∂x, ux,t, and ∂ux,t/∂t from (52)-(55).Step 5. Find an approximation to 0CDtσux,t from (57).Step 6. Determine the residual function from (61).Step 7. Obtain roots of YN+1x and YN+1t (xi,tj,i,j=0,1,⋯,N) using fsolve command in Maple.Step 8. Collocate the residual function at tensor points xi,tj,i,j=0,1,⋯,N.Step 9. Solve the resultant non-linear system in Step 8 by Newton’s iteration method and obtain the unknown vector C.Step 10. Find uNx,t from (54).Output:uNx,t


## 3.1. One-Variable Chebyshev Polynomials of the Sixth Kind

The following recurrence relation holds for the sixth-kind Chebyshev polynomials(6)Y¯jz=zY¯j−1z−θj+1Y¯j−2z,j≥2,z∈−1,1,Y¯0z=1,Y¯1z=z,where(7)θj=jj+1+−1j2j+1+14jj+1.These polynomials are orthogonal with respect to the weight functionw¯z=z21−z2, that is,
(8)∫−11Y¯izY¯jzz21−z2dz=ℏ¯iδij,where(9)ℏ¯i=π22i+3,ieven,πi+322i+3i+1,iodd.By the change of variablez=2t−1, the shifted Chebyshev polynomials Yit=Y¯i2t−1 are orthogonal regarding the weight function wt=2t−12t−t2 on the interval 0,1,
(10)∫01YitYjt2t−12t−t2dt=ℏiδij,and(11)ℏ¯i=π22i+5,ieven,πi+322i+5i+1,iodd.The series form of the shifted Chebyshev polynomials of the sixth kind is as follows:(12)Yjt=∑r=0jςr,jtr,where(13)ςr,j=22r−j2r+1!∑l=r+1/2j/2−1j/2+l+r2l+r+1!2l−r+1!,jeven,2j+1∑l=r/2j−1/2−1j+1/2+l+rl+12l+r+2!2l−r+1!,jodd.Every square-integrable functionvt∈Lw2I, I=0,1 can be expanded in the shifted sixth-kind Chebyshev polynomials as
(14)vt=∑j=0∞VjYjt,t∈I,where the coefficientsVj are computed as
(15)Vj=1ℏj∫01vtYjtwtdt.The first few coefficients in (14) practically keep information of function vt. In other words, a finite series can present an approximation to vt as
(16)vt≈vNt=∑j=0NVjYjt=YTtV=VTYt,whereV and Yt are the N+1×1 vectors as follows:
(17)V=V0V1⋯VNT,Yt=Y0tY1t⋯YNtT.


## 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
(18)∫01tνYktwtdt=∑m=0kςm,kΓ324Γν+m+7/2Γν+m+5−4Γν+m+5/2Γν+m+4+Γν+m+3/2Γν+m+3.Proof.
By the series form of the shifted sixth-kind Chebyshev polynomials in (12) and the weight function wt, one has
(19)∫01tνYktwtdt=∫01tν∑m=0kςm,ktm2t−12t−t2dt=∑m=0kςm,k∫014tν+m+5/2−4tν+m+3/2+xν+m+1/21−t1/2dt=∑m=0kςm,k4Bν+m+72,32−4Bν+m+52,32+Bν+m+32,32,
whereBr,s is the well-known beta function, so, the desired result is achieved.Theorem 5.
IfYt is the basis vector in (17), the integral of Yt can be computed as
(20)∫0tYsds≈PYt,t∈I,
whereP is the N+1×N+1 integral operational matrix of the integer-order in the following form:
(21)P=π0,0π0,1⋯π0,Nπ1,0π1,1⋯π1,N⋮⋮⋱⋮πj,0πj,1⋯πj,N⋮⋮⋱⋮πN,0πN,1⋯πN,N,
where the entriesπj,k are computed as
(22)πj,k=∑r=0jςr,jΓ3/2r+1ℏk∑m=0kςm,k4Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4,j=0,1,⋯,N,k=0,1,⋯,N.Proof.
Integrating the elements of the vectorYt yields
(23)∫0tYjsds=∑r=0jςr,j∫0tsrds=∑r=0jςr,jtr+1r+1,j=0,1,⋯,N.
Now,tr+1 is approximated in terms of the shifted sixth-kind Chebyshev polynomials
(24)tr+1≈∑k=0Nρk,r+1Ykt,
where(25)ρk,r+1=1ℏk∫01tr+1Yktwtdt.
Using Lemma4, the integral part of (25) is computed as follows:
(26)∫01tr+1Yktwtdt=∑m=0kςm,kΓ324Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4.
Therefore, (23) is written as
(27)∫0tYjsds≈∑k=0N∑r=0jςr,jΓ3/2r+1ℏk∑m=0kςm,k4Γr+m+9/2Γr+m+6−4Γr+m+7/2Γr+m+5+Γr+m+5/2Γr+m+4Ykt.
By rewriting the last series as a matrix form, the desired result is achieved.Theorem 6.
Assume thatYt is the basis vector in (17) and I0RLtσ is the Riemann-Liouville integral operator of the order σ, 0<σ<1. Then, one has
(28)I0RLtσYt≈PσYt,
wherePσ is the N+1×N+1 fractional operational matrix of the order σ as follows:
(29)Pσ=πσ0,0πσ0,1⋯πσ0,Nπσ1,0πσ1,1⋯πσ1,N⋮⋮⋱⋮πσi,0πσi,1⋯πσi,N⋮⋮⋱⋮πσN,0πσN,1⋯πσN,N,
where the entriesπσi,k are computed as
(30)πσi,k=∑r=0iςr,iΓr+σ+1Γ3/2Γr+σ+2ℏk∑m=0kςm,k4Γr+m+7/2Γr+m+5−4Γr+m+5/2Γr+m+4+Γr+m+3/2Γr+m+3,i=0,1,⋯,N,k=0,1,⋯,N.Proof.
The proof process is similar to Theorem5. Noting the definition of 0RLItσ and its properties in (5), the fractional integral of Yit is computed as
(31)I0RLtσYit=∑r=0iςr,iΓr+σ+1tr+σΓr+σ+2,i=0,1,⋯,N.
Now,tr+σ is approximated by the Chebyshev polynomials of the sixth kind as
(32)tr+σ≈∑k=0Nρk,rYkt,s.t.ρk,r=1ℏk∫01tr+σYktwtdt.
By Lemma4 and pursuing the proof process in Theorem 5, Equation (31) is written as
(33)I0RLtσYit≈∑k=0N∑r=0iςr,iΓr+σ+1Γ3/2Γr+σ+2ℏk×∑m=0kςm,k4Γr+σ+m+7/2Γr+σ+m+5−4Γr+σ+m+5/2Γr+σ+m+4+Γr+σ+m+3/2Γr+σ+m+3Ykt,i=0,1,⋯,N.Theorem 7.
Assume thatYt is the basis vector and Yqt, 0<q<1 is its delay form. Yqt can be approximated in Yt as
(34)Yqt≈LYt,
whereL is a N+1×N+1 matrix as follows:
(35)L=100⋯0l01q0⋯0l02l12q⋯0⋮⋮⋮⋱⋮l0Nl1Nl2N⋯qN,
andlkj is computed from the following recurrence formulas:
(36)l0j+1=q−1l0j+qθ2l1j−θjl0j−1,j=1,2,⋯,N−1,lkj+1=qlk−1j+q−1lkj+qθk+1lk+1j−θjlkj−1,k=1,2,⋯,j−1,ljj+1=qlj−1j+q−1qj,j=1,2,⋯,N−1,
with the starting valuesljj=qj,j=0,1,⋯,N, and l01=q−1.Proof.
Consider the following recurrence formula obtained from Equation (6)
(37)Yj+1qt=2qt−1Yjqt−θj+1Yj−1qt,j≥1.
Also, the following auxiliary relation is obtained from formula (6):
(38)2tYjt=Yj+1t+Yjt+θj+1Yj−1t,j≥1.
Now,Yjqt can be expanded in terms of the Chebyshev polynomials of the sixth kind as follows:
(39)Yjqt=∑k=0jlkjYkt,j=0,1,⋯,N.
From the first few polynomials in Equation (39), it is easily obtained l01=q−1, ljj=qj,j=0,1,⋯,N. Using auxiliary relation (38) and substituting Equation (39) into Equation (37), one gets
(40)∑k=0j+1lkj+1Ykt=q∑k=0jlkjYk+1t+Ykt+θk+1Yk−1t−θj+1∑k=0j−1lkj−1Ykt−∑k=0jlkjYkt=∑k=0jlkjqYk+1t+q−1Ykt+qθk+1Yk−1t−θj+1∑k=0j−1lkj−1Ykt.
Equating coefficients ofYkt 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Ω=0,1×0,1 as
(41)Wijx,t=YixYjt,i,j=0,1,⋯,x,t∈Ω.These polynomials are orthogonal regarding the weight functionωx,t=wxwt on Ω,
(42)∫01∫01Wijx,tWklx,tωx,tdxdt=ℏiℏjδikδjl,whereℏi and ℏj are calculated by (11). The function Vx,t∈Lω2Ω is expanded as
(43)Vx,t=∑i=0∞∑j=0∞VijWijx,t,x,t∈Ω,and a truncated series of (43) is considered as an approximation to the function Vx,t,
(44)Vx,t≈VNx,t=∑i=0N∑j=0NVijWijx,t=∑i=0N+12−1Vi∗Wi∗x,t=VTWx,t=WTx,tV,whereVi∗=Vrs,Wi∗x,t=Wrsx,t such that r=i/N+1,s=i−rN+1, and V,Wx,t are the N+12×1 vectors as
(45)V=V00V01⋯V0NV10V11⋯V1N⋯VN0VN1⋯VNNT,Wx,t=W00x,tW01x,t⋯W0Nx,tW10W11x,t⋯W1Nx,t⋯WN0x,tWN1x,t⋯WNNx,tT.


## 3.4. Operational Matrices of Two-Variable Basis

Consider the two-variable basisWx,t in (45). The integral operational matrices of Wx,t with respect to variables x and t are obtained, respectively, as
(46)∫0xWs,tds≈ℙx1Wx,t=P⊗IWx,t,∫0tWx,τdτ≈ℙt1Wx,t=I⊗PWx,t,whereℙx1 and ℙt1 are the N+12×N+12 integral operational matrices related to x and t, respectively, P is the operational matrix in Theorem 5, and I is the N+1×N+1 identity matrix. Similarly, the fractional integral of Wx,t of the order σ with respect to t can be computed as
(47)I0RLtσWx,t≈ℙtσWx,t=I⊗PσWx,t,whereℙtσ is the N+12×N+12 fractional integral operational matrix of the order σ related to t, Pσ is the operational matrix in Theorem 6, and I is the N+1×N+1 identity matrix. Now, the relationship of the vectors Wx,qt, Wqx,t, and Wqx,qt to the basis vector Wx,t is specified.By settingt=qt in Wx,t, the vector Wx,qt is rewritten as follows:
(48)Wx,qt=Y0xY0qt,Y0xY1qt,⋯,Y0xYNqt,Y1xY0qt,Y1xY1qt,⋯,Y1xYNqt,⋯,YNxY0qt,YNxY1qt,⋯,YNxYNqtT=Y0xY0qt,Y1qt,⋯,YNqt,Y1xY0qt,Y1qt,⋯,YNqt,⋯,YNxY0qt,Y1qt,⋯,YNqtT=Y0xYqt,Y1xYqt,⋯,YNxYqtT≈Y0xLYt,Y1xLYt,⋯,YNxLYtT=LON+1⋯ON+1ON+1L⋯ON+1⋮⋮⋱⋮ON+1ON+1⋯LY0xY0t⋮YNxY0t⋮YNxYNt=I⊗LWx,t=L∗∗Wx,t,whereL is the matrix in Theorem 7, ON+1 is the N+1×N+1 zero matrix, and L∗∗ is the N+12×N+12 delay matrix. Similarly, it is found that
(49)Wqx,t≈L⊗IWx,t=L∗Wx,t,Wqx,qt≈L∗L∗∗Wx,t=LWx,t,L=L∗L∗∗,whereL∗ and L are N+12×N+12 matrices.


## 3.5. Solution Method

To describe the methodology, three forms of Equations (1) and (2) are considered [26, 27]:Form I:(50)D0Ctσux,t−∂2ux,t∂x2−∂ux,t/2∂xux2,t2−12ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x,u0,t=0,ux0,t=expt.According to the highest orders of derivatives regardingx and t, the following approximation is considered:
(51)∂3ux,t∂x2∂t≈WTx,tC.Integrating (51) concerning t and x, respectively, leads to the following approximations:
(52)∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC,(53)∂ux,t∂x≈WTx,tℙx1Tℙt1TC+∂u0,t∂x≈WTx,tℙx1Tℙt1TC+WTx,tF,(54)ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1TF+u0,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1TF.Now, twice integrating approximation (51) with respect to x leads to an approximation for ∂ux,t/∂t:
(55)∂2ux,t∂x∂t≈WTx,tℙx1TC+WTx,tF,∂ux,t∂t≈WTx,tℙx1T2C+WTx,tℙx1TF.To obtain an approximation to0CDtσux,t, approximation (55) is rewritten as follows:
(56)∂ux,t∂t=∂1−σ∂σux,t∂t1−σ∂tσ≈WTx,tℙx1T2C+WTx,tℙx1TF.By applying the Riemann-Liouville operator of the order1−σ to both sides of (56), one gets
(57)D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C+WTx,tℙt1−σTℙx1TF+0CDtσux,0.Fractionally differentiating approximation (54) and setting t=0 lead to
(58)D0Ctσux,0≈0CDtσWTx,tℙx1T2ℙt1TC+0CDtσWTx,tℙx1TFt=0=0.Approximation (58) equals zero because after fractionally differentiating Wx,t related to t, all components of the basis vector involve terms as tμ,1−σ<μ<N−σ or are zero.The terms with delays can be approximated as(59)∂ux,t/2∂x≈WTx,t2ℙx1Tℙt1TC+WTx,t2F≈WTx,tL∗∗Tℙx1Tℙt1TC+WTx,tL∗∗TF,(60)ux2,t2≈WTx2,t2ℙx1T2ℙt1TC+WTx2,t2ℙx1TF≈WTx,tLTℙx1T2ℙt1TC+WTx,tLTℙx1TF.Substituting approximations (52)–(60) into Equation (50) results in the following residual function:
(61)RNx,t=WTx,tℙt1−σTℙx1T2C+WTx,tℙt1−σTℙx1TF−WTx,tℙt1TC−WTx,tL∗∗Tℙx1Tℙt1TC+WTx,tL∗∗TFWTx,tLTℙx1T2ℙt1TC+WTx,tLTℙx1TF−12WTx,tℙx1T2C+WTx,tℙx1TF≈0.Form II:(62)D0Ctσux,t−∂2ux,t/2∂x2ux,t2+ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x2,u0,t=0,ux0,t=0.The functions in Equation (62) are approximated based on what was done for Equation (50):
(63)∂3ux,t∂x2∂t≈WTx,tC,∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC+WTx,tF,ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F,D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C,∂2ux,t/2∂x2≈WTx,tL∗∗Tℙt1TC+WTx,tL∗∗TF,ux,t2≈WTx,tL∗∗Tℙx1T2ℙt1TC+WTx,tL∗∗Tℙx1TF.Substituting approximations (63) into Equation (62) leads to the following residual function:
(64)RNx,t=WTx,tℙt1−σTℙx1T2C−WTx,tL∗∗Tℙt1TC+WTx,tL∗∗TF×WTx,tL∗∗Tℙx1T2ℙt1TC+WTx,tL∗∗Tℙx1TF+WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F≈0.Form III:(65)D0Ctσux,t−∂2ux/2,t/2∂x2∂ux/2,t/2∂x+∂ux,t∂x+ux,t=0,0<σ⩽1,x,t∈Ω,ux,0=x2,u0,t=0,ux0,t=0.The following approximations can be obtained for the functions in Equation (65):
(66)∂3ux,t∂x2∂t≈WTx,tC,∂2ux,t∂x2≈WTx,tℙt1TC+∂2ux,0∂x2≈WTx,tℙt1TC+WTx,tF,∂ux,t∂x≈WTx,tℙx1Tℙt1TC+WTx,tℙx1TF,ux,t≈WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F,D0Ctσux,t≈WTx,tℙt1−σTℙx1T2C,∂2ux/2,t/2∂x2≈WTx,tLTℙt1TC+WTx,tLTF,∂ux/2,t/2∂x≈WTx,tLTℙx1Tℙt1TC+WTx,tLTℙx1TF.By substituting approximations (66) into Equation (65), one gets the following residual function:
(67)RNx,t=WTx,tℙt1−σTℙx1T2C−WTx,tLTℙt1TC+WTx,tLTF×WTx,tLTℙx1Tℙt1TC+WTx,tLTℙx1TF+WTx,tℙx1Tℙt1TC+WTx,tℙx1TF+WTx,tℙx1T2ℙt1TC+WTx,tℙx1T2F≈0.The residual functions in Equations (61), (64), and (65) are collocated at collocation nodes xi,tji,j=0∞, which xi and tj are roots of YN+1x and YN+1t, respectively. Hence, a nonlinear system involving N+12 algebraic equations is achieved that can be solved by the Newton’s iteration method. Therefore, the coefficient vector C is determined approximately; then, an approximate solution, ux,t, is achieved. To better describe the solution method, pursuing Algorithm 1 for Equation (50) is suggested.Algorithm 1: Solution method
Input:σ,NStep 1. Derive operational matrices P, Pσ, and L from (21), (29), and (36).Step 2. Construct the operational matrices ℙx1, ℙt1, ℙtσ, L, L∗, L∗∗ from (46)-(49).Step 3. Consider the approximation ∂3ux,t/∂x2∂t≈WTx,tC in (51).Step 4. Find approximations to ∂2ux,t/∂x2, ∂ux,t/∂x, ux,t, and ∂ux,t/∂t from (52)-(55).Step 5. Find an approximation to 0CDtσux,t from (57).Step 6. Determine the residual function from (61).Step 7. Obtain roots of YN+1x and YN+1t (xi,tj,i,j=0,1,⋯,N) using fsolve command in Maple.Step 8. Collocate the residual function at tensor points xi,tj,i,j=0,1,⋯,N.Step 9. Solve the resultant non-linear system in Step 8 by Newton’s iteration method and obtain the unknown vector C.Step 10. Find uNx,t from (54).Output:uNx,t


## 4. Error Analysis

It is found that the value of the error functionENx,t decreases when values of N increase. First, some error bounds are obtained for the unknown function ux,t and its derivatives.Theorem 8.
Assume that∂i+jux,t/∂xi∂tj∈CΩ,i,j=0,1,⋯,N+1,ΘN=spanWijx,t,i,j=0,1,⋯,N, and uNx,t is the approximate solution obtained from the method belonging to ΘN, and
(68)τN=supx,t∈Ω∂2N+1ux,t∂xN+1∂tN+1.
Then, the error bound of theuNx,t can be obtained as
(69)ux,t−uNx,tLω2Ω≤A0τNπ25/2N3/2ΓN+22.Proof.
Define the bivariate Taylor expansion ofuNx,t as
(70)TNx,t=∑i=0N∑j=0N−i∂i+jux,t∂xi∂tj0,0xitjΓi+1Γj+1,
and(71)ux,t−TNx,t=xN+1tN+1ΓN+22∂2N+1uξ,η∂xN+1∂tN+1,ξ,η∈Ω.
SinceuNx,t is the best approximate solution of ux,t from ΘN and according to Equation (71), one gets
(72)ux,t−uNx,tLω2Ω2≤ux,t−TNx,tLω2Ω2≤∫01∫01τN2x2N+1t2N+1ΓN+24ωx,tdxdt=τN2ΓN+24∫014x2N+9/21−x1/2−4x2N+7/21−x1/2+4x2N+5/21−x1/2dx×∫014t2N+9/21−t1/2−4t2N+7/21−t1/2+4t2N+5/21−t1/2dt=τN2π4ΓN+244Γ2N+11/2Γ2N+7−4Γ2N+9/2Γ2N+6+Γ2N+7/2Γ2N+52.
By the Stirling formula [29], some bounds are computed for the last equality:
(73)ux,t−uNx,tLω2Ω≤τNπ2ΓN+224α12N−3/2−4α22N−3/2+α32N−3/2≤A0τNπ25/2N3/2ΓN+22,
whereA0=4α1−4α2+α3. Therefore, the desired result is achieved.Theorem 9.
Assume thatux,t, uNx,t, and ∂i+jux,t/∂xi∂tj, i,j=0,1,⋯,N satisfy the condition of Theorem 8, and set
(74)θN,k=supx,t∈Ω∂2N−k+2ux,t∂x2N−k+1∂tN+1,k=1,2.
Then,(75)∂kux,t∂xk−∂kuNx,t∂xkLω2Ω≤A1,kθN,kπ25/2NN−k3/4ΓN−k+2ΓN+2,k=1,2.Proof.
The bivariate Taylor expansion of∂kux,t/∂xk leads to
(76)∂kux,t∂xk−∂kuNx,t∂xk=xN−k+1tN+1ΓN−k+2ΓN+2∂2N−k+2uξ′,η′∂xN−k+1∂tN+1,ξ′,η′∈Ω,k=1,2.
Therefore, by takingL2-norm and using Equation (74) and Stirling formula, one has
(77)∂kux,t∂xk−∂kuNx,t∂xkLω2Ω2≤∫01∫01θN,k2x2N−k+1t2N+1ΓN−k+22ΓN+22ωx,tdxdt=θN,k2π4ΓN−k+22ΓN+22×4Γ2N−k+11/2Γ2N−k+7−4Γ2N−k+9/2Γ2N−k+6+Γ2N−k+7/2Γ2N−k+5×4Γ2N+11/2Γ2N+7−4Γ2N+9/2Γ2N+6+Γ2N+7/2Γ2N+5≤θN,k2π4ΓN−k+2ΓN+24β1N−k−3/2−4β2N−k−3/2+β3N−k−3/2×4γ1N−3/2−4γ2N−3/2+γ3N−3/2≤A1,k2θN,k2π4ΓN−k+22ΓN+222N−k−3/22N−3/2,
whereA1,k2=4β1−4β2+β34γ1−4γ2+γ3.Theorem 10.
Assume thatD0Ctσux,t∈CΩ and the conditions of Theorem 8 hold. Then,
(78)D0Ctσux,t–D0CtσuNx,tLω2Ω≤A2τNπ25/2ΓN+2ΓN−σ+2NN−σ3/4.Proof.
According to Equation (71) and properties of the Caputo operator in (5), one can write
(79)D0Ctσux,t−uNx,t≤τNxN+1tN−σ+1ΓN−σ+1ΓN+2.
TakingL2-norm yields
(80)D0Ctσux,t−D0CtσuNx,tLω2Ω≤∫01∫01τN2x2N+1t2N−σ+1ΓN−σ+22ΓN+22ωx,tdxdt=τN2π4ΓN−σ+22ΓN+22×4Γ2N+11/2Γ2N+7−4Γ2N+9/2Γ2N+6+Γ2N+7/2Γ2N+5×4Γ2N−σ+11/2Γ2N−σ+7−4Γ2N−σ+9/2Γ2N−σ+6+Γ2N−σ+7/2Γ2N−σ+5≤A2τN2π4ΓN−σ+22ΓN+222N−3/22N−σ−3/2.Corollary 11.
Some error bounds for functions with proportional delays can be obtained using the resultant bounds in Theorems8, 9, and 10.
(81)ux,t2−uNx,t2Lω2Ω≤A0τNπ2N+7/2N3/2ΓN+22,ux2,t2−uNx2,t2Lω2Ω≤A0τNπ22N+9/2N3/2ΓN+22,∂ux,t/2∂x−∂uNx,t/2∂xLω2Ω≤A1,1θN,1π2N+7/2NN−13/4ΓN+1ΓN+2,∂2ux,t/2∂x2−∂2uNx,t/2∂x2Lω2Ω≤A1,2θN,2π2N+7/2NN−23/4ΓNΓN+2,∂ux/2,t/2∂x−∂uNx/2,t/2∂xLω2Ω≤A1,1θN,1π22N+9/2NN−13/4ΓN+1ΓN+2,∂2ux/2,t/2∂x2−∂2uNx/2,t/2∂x2Lω2Ω≤A1,2θN,2π22N+9/2NN−23/4ΓNΓN+2.Now, Theorems8, 9, and 10 and Corollary 11 are applied to show that the error of the method becomes sufficiently small when N is sufficiently large. For this, three equations (50), (62), and (65) in Section 4 are called again. Moreover, consider the following bounds:
(82)uNx2,t2Lω2Ω≤U1,∂uNx,t/2∂xLω2Ω≤U2,∂2uNx,t/2∂x2Lω2Ω≤U3,uNx,t2Lω2Ω≤U4,∂uNx/2,t/2∂xLω2Ω≤U5,∂2uNx/2,t/2∂x2Lω2Ω≤U6.Form I. SupposeuNx,t is the best approximate solution from ΘN which is obtained from the proposed scheme. Therefore, it satisfies the following equation:
(83)D0CtσuNx,t−∂2uNx,t∂x2−∂uNx,t/2∂xuNx2,t2−12uNx,t=ENx,t,whereENx,t is the error term. Subtracting Equation (83) from Equation (50) leads to the following equation:
(84)ENx,t=D0Ctσux,t−D0CtσuNx,t−∂2ux,t∂x2−∂2uNx,t∂x2−∂ux,t/2∂xux2,t2−∂uNx,t/2∂xuNx2,t2−12ux,t−uNx,t=D0Ctσux,t−D0CtσuNx,t−∂2ux,t∂x2−∂2uNx,t∂x2−∂ux,t/2∂x−∂uNx,t/2∂xux2,t2−uNx2,t2+uNx2,t2−∂uNx,t/2∂xux2,t2−uNx2,t2−12ux,t−uNx,t.So, one has(85)ENx,tLω2ΩA2τNπ25/2ΓN+2ΓN−σ+2NN−σ3/4+A1,2θN,2π25/2ΓNΓN+2NN−23/4+A1,1θN,1π2N+7/2ΓN+1ΓN+2NN−13/4A0τNπ22N+9/2ΓN+12N3/2+U1+U2A0τNπ22N+9/2ΓN+12N3/2+12A0τNπ25/2ΓN+12N3/2.Form II. IfuNx,t is the solution obtained from the proposed algorithm for Equation (62), then it satisfies the following equation:
(86)ENx,t=D0CtσuNx,t−∂2uNx,t/2∂x2uNx,t2+uNx,t.Subtracting Equation (86) from Equation (62) leads to
(87)ENx,t=D0Ctσux,t−D0CtσuNx,t−∂2ux,t/2∂x2ux,t2−∂2uNx,t/2∂x2uNx,t2+ux,t−uNx,t=D0Ctσux,t−D0CtσuNx,t−∂2ux,t/2∂x2−∂2uNx,t/2∂x2ux,t2−uNx,t2+uNx,t2−∂2uNx,t/2∂x2ux,t2−uNx,t2+ux,t−uNx,t.Then, one has(88)ENx,tLω2Ω≤A2τNπ25/2ΓN+2ΓN−σ+2NN−σ3/4+A1,2θN,2π2N+7/2ΓNΓN+2NN−23/4A0τNπ2N+7/2ΓN+22N3/2+U4+U3A0τNπ2N+7/2ΓN+22N3/2+A0τNπ25/2ΓN+22N3/2.Form III. IfuNx,t is an approximate solution obtained from the suggested algorithm for Equation (65), then one has
(89)ENx,t=D0CtσuNx,t−∂2uNx/2,t/2∂x2∂uNx/2,t/2∂x+∂uNx,t∂x+uNx,t.Subtracting (89) from Equation (65) leads to
(90)ENx,t=D0CtσuNx,t−D0CtσuNx,t−∂2ux/2,t/2∂x2∂ux/2,t/2∂x−∂2uNx/2,t/2∂x2∂uNx/2,t/2∂x+∂ux,t∂x−∂uNx,t∂x+ux,t−uNx,t=D0CtσuNx,t−D0CtσuNx,t−∂2ux/2,t/2∂x2−∂2uNx/2,t/2∂x2∂ux/2,t/2∂x−∂uNx/2,t/2∂x+∂uNx/2,t/2∂x−∂2uNx/2,t/2∂x2∂ux/2,t/2∂x−∂uNx/2,t/2∂x+∂ux,t∂x−∂uNx,t∂x+ux,t−uNx,t.By takingL2-norm, one gets
(91)ENx,tLω2ΩA2τNπ25/2ΓN+2ΓN−σ+2NN−σ3/4+A1,2θN,2π22N+9/2ΓNΓN+2NN−23/4U5+A1,1θN,1π22N+9/2ΓN+1ΓN+2NN−13/4+U6A1,1θN,1π22N+9/2ΓN+1ΓN+2NN−13/4+A1,1θN,1π25/2ΓN+1ΓN+2NN−13/4+A0τNπ25/2ΓN+22N3/2.As seen from the right-hand sides of inequalities (85), (88), and (91), ENx,t will decrease by choosing appropriate values of N, i.e., the error bound will be sufficiently small for the sufficiently large values of N.


## 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 ux,t=xexpt for σ=1 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 N=6 and σ=0.7,0.8,0.9,1. Also, the figures of approximate solutions at t=1 are seen in Figure 2(b) for N=6 and σ=0.7,0.8,0.9,1. 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 N=6 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 N and σ=1, and results are observed in Table 3. Increasing N 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 ux,t=x2expt if σ=1. The 3D figures of the approximate solutions are depicted in Figure 3(a) for N=6 and σ=0.7,0.8,0.9,1. The plot of the absolute error function is seen in Figure 3(b). The approximate solutions are plotted at t=1 and compared to those presented by [26] in Figure 4 for N=6 and σ=0.7,0.8,0.9,1. 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 N=6 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 σ=1 and N=4,5,6,7,8.Example 3.
Consider Form III in (65) and residual function (67). The exact solution is ux,t=x2exp−t if σ=1. The three-dimensional figures of the approximate solutions are depicted in Figures 5(a)–5(d) for N=8 and σ=0.7,0.8,0.9,1. The plot of the absolute error function is seen in Figure 5(e) for N=8 and σ=1. The obtained solutions are plotted in Figure 6 at t=1 for N=8 and σ=0.7,0.8,0.9,1 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 N=8 and various values of σ. As σ⟶1, the numerical solutions converge to the exact one. The obtained results are compared to those reported in [26, 27] in Table 8 for N=8 and σ=1. As seen, the proposed method presents better accuracy.Figure 1
Surface behavior of solutions of Example1 for N=6 and various values of σ.
(a)(b)(c)(d)Figure 2
(a) Absolute error function for Example1 for N=6, σ=1. (b) Plots of solutions for different values of σ, N=6,t=1.
(a)(b)Table 1
Numerical results of Example1 for N=6 and various values of σ.
xitjσ=0.7σ=0.8σ=0.9σ=1ExactError0.250.250.31998660.32031020.32065180.32100640.32100641.3966×10−90.500.41140500.41168340.41194570.41218030.41218038.0055×10−100.750.52856560.52881070.52903950.52925000.52925009.2143×10−101.000.67887370.67906520.67926340.67957040.67957052.5998×10−80.500.250.63375800.63626310.63897490.64201270.64201272.5203×10−90.500.81875310.82097530.82292750.82436060.82436061.4904×10−90.751.05286061.05473611.05658111.05850001.05850008.8651×10−101.001.35128391.35218271.35400141.35914091.35914095.3534×10−80.750.250.93404110.94174540.95077900.96301910.96301911.9619×10−90.501.22159251.22934261.23507741.23654101.23654107.2547×10−100.751.56726521.57296401.57937731.58775001.58775005.3519×10−101.001.99364841.98886271.99379732.03871132.03871149.5450×10−81.00.251.20853471.22277371.24372341.28402541.28402544.6333×10−90.501.62984711.65060661.66128031.64872131.64872136.2784×10−90.752.06089362.07148572.08783372.11700002.11700002.3318×10−91.002.51938192.45854622.45180892.71828162.71828182.1185×10−7Table 2
Numerical results of Example1 for N=6, σ=1.
xitjMethod in [26]Method in [27]Chebyshev methodExact0.250.250.32100420.32100420.32100640.32100640.500.41210940.41210940.41218030.41218030.750.52868650.52868650.52925000.52925001.000.67708330.67708330.67957040.67957050.500.250.64200850.64200850.64201270.64201270.500.82421880.82421870.82436060.82436060.751.05737301.05737301.05850001.05850001.001.35416671.35416671.35914091.35914090.750.250.96301270.96301270.96301910.96301910.501.2363281.23632811.23654101.23654100.751.5860601.58605961.58775001.58775001.002.03125002.03124902.03871132.0387114Table 3
Maximum absolute errors of Example1 for σ=1 and various values of N.
N45678Error7.5579×10−54.5135×10−62.1185×10−77.7009×10−93.6459×10−10Figure 3
(a) Surface solutions forN=6 and σ=1 (blue), σ=0.9 (green), σ=0.8 (gray), and σ=0.7 (red). (b) Absolute error function of Example 2 for N=6, σ=1.
(a)(b)Figure 4
Left: approximate solutions depicted in [26]. Right: obtained solutions from the proposed method for different values of σ, N=6,t=1 for Example 2.
(a)(b)Table 4
Numerical results of Example2 for various values of σ, N=6.
xitjσ=0.7σ=0.8σ=0.9σ=1ExactError0.250.250.10132350.09144260.08490360.08025160.08025168.0280×10−100.500.14931510.12660510.11257840.10304510.10304511.5944×10−90.750.21649300.17340920.14850010.13231250.13231251.9099×10−101.000.31370110.23724360.19569910.16989260.16989266.6257×10−90.500.250.40529410.36577050.33961420.32100640.32100643.2112×10−90.500.59726050.50642040.45031380.41218030.41218036.3777×10−90.750.86597200.69363680.59400050.52925000.52925007.6358×10−101.001.25480430.94897430.78279650.67957040.67957052.6504×10−80.750.250.91191170.82298370.76413210.72226430.72226437.2249×10−90.501.34383621.13944601.01320590.92740570.92740571.4348×10−80.751.94843701.56068291.33650101.19081251.19081251.7060×10−91.002.82330962.13519211.76129221.52903351.52903355.9679×10−81.00.250.33715091.46308211.35845701.28402541.28402541.2842×10−80.502.38904212.02568181.80125501.64872141.64872132.5483×10−80.753.46388802.77454742.37600192.11700002.1170002.9197×10−91.005.01921693.79589713.13118622.71828172.71828181.0649×10−7Table 5
Numerical results of Example2 for N=6, σ=1.
xitjMethod in [26]Method in [27]Chebyshev MethodExact0.250.250.08025110.08025160.08025160.08025160.500.10302730.10304510.10304510.10304510.750.13217160.13231230.13231250.13231251.000.16927080.16989090.16989260.16989260.500.250.32100420.32100640.32100640.32100640.500.41210940.41218030.41218030.41218030.750.52868650.52924930.52925000.52925001.000.67708330.67956350.67957040.67957050.750.250.7222595−0.72226430.72226430.500.9272461−0.92740570.92740570.751.1895447−1.19081251.19081251.001.5234375−1.52903351.52903351.000.25−1.28402541.28402541.28402540.50−1.64872121.64872131.64872130.75−2.11699732.11700002.11700001.00−2.71825402.71828172.7182818Table 6
Maximum absolute errors of Example2 for σ=1 and various values of N.
N45678Error2.3903×10−52.3083×10−61.0149×10−71.6453×10−82.8933×10−9Figure 5
(a)–(d) Surface solutions forN=8 and σ=1,0.9,0.8,0.7. (e) Absolute error function of Example 3 for N=8, σ=1.
(a)(b)(c)(d)(e)Figure 6
Left: approximate solutions depicted in [26]. Right: obtained solutions from the proposed method for different values of σ, N=8,x=1 for Example 3.
(a)(b)Table 7
Numerical results of Example3 for various values of σ, N=8.
xitjσ=0.7σ=0.8σ=0.9σ=1ExactError0.250.250.03577800.04200240.04565420.04867500.04867507.6753×10−120.500.02051740.02865990.03330340.03790820.03790821.6649×10−100.750.01127790.01899740.02384300.02952290.02952293.8695×10−101.000.00697930.01272590.01699060.02299250.02299255.6253×10−100.500.250.16604900.17488870.18544950.19470020.19470021.1877×10−100.500.12794830.13531760.14413040.15163270.15163277.2917×10−100.750.09463320.10521570.11242820.11809160.11809161.1956×10−91.000.06965680.08203490.08723060.09196990.09196991.2032×10−90.750.250.36608030.38890220.41490600.43807540.43807542.6901×10−100.500.29807050.29875020.31937890.34117350.34117354.0079×10−100.750.25592710.24239330.25466680.26570620.26570621.0422×10−91.000.22399210.20489910.20815810.20693220.20693223.1164×10−91.00.250.64014320.69328950.73982430.77880080.77880086.2071×10−110.500.48740580.51559390.56094730.60653070.60653071.5635×10−90.750.41458840.40157280.43318950.47236650.47236663.2396×10−91.000.39137230.33641940.34920840.36787940.36787941.8734×10−9Table 8
Numerical results of Example3 for N=8, σ=1.
xitjMethod in [26]Method in [27]Chebyshev methodExact0.250.250.04867550.04867500.04867500.04867500.500.03792320.03790810.03790820.03790820.750.02963260.02952280.02952290.02952291.000.02343750.02299110.02299250.02299250.500.250.19470210.19470020.19470020.19470020.500.15169270.15163260.15163270.15163270.750.11853030.11809110.11809160.11809161.000.09375000.09196430.09196990.09196990.750.250.4380798−0.43807540.43807540.500.3413086−0.34117350.34117350.750.2666931−0.26570620.26570621.000.2109375−0.20693220.20693221.000.25−0.77880080.77880080.77880080.50−0.60653060.60653070.60653070.75−0.47236430.47236650.47236651.00−0.36785710.36787940.3678794


## 6. Conclusion

This paper deals with numerically solving a class of fractional partial differential equations with proportional delays on the domainΩ=0,1×0,1. 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 N 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.



---

*Source: 9512048-2022-01-07.xml*