Bivariate Chebyshev polynomials of the fifth kind for variable-order time-fractional partial integro-differential equations with weakly singular kernel

The shifted Chebyshev polynomials of the fifth kind (SCPFK) and the collocation method are employed to achieve approximate solutions of a category of the functional equations, namely variable-order time-fractional weakly singular partial integro-differential equations (VTFWSPIDEs). A pseudo-operational matrix (POM) approach is developed for the numerical solution of the problem under study. The suggested method changes solving the VTFWSPIDE into the solution of a system of linear algebraic equations. Error bounds of the approximate solutions are obtained, and the application of the proposed scheme is examined on five problems. The results confirm the applicability and high accuracy of the method for the numerical solution of fractional singular partial integro-differential equations.


Introduction
Partial integro-differential equations (PIDEs) with weakly singular kernels emerge in some physical and chemical phenomena, such as the radiation of heat from semi-infinite solids, stereology, hydrodynamics, heat conduction, and theory of elasticity [1,2]. Due to possessing singular kernels, this category of functional equations is complicated, and directly obtaining exact solutions of singular PIDEs is usually extremely hard. Therefore, the study of numerical solutions of these equations is severely significant. Numerous numerical schemes have been established to numerically solve the singular PIDEs, such as the Sinc methods [3][4][5][6], the optimum q-homotopic analysis method [7], the finite difference methods [8][9][10], wavelets collocation schemes [11], and many other methods (readers can refer to [12][13][14][15]). Fractional derivative and integral operators can better describe physical phenomena of the real world, and hence researchers handle solving diverse fractional functional equations. For example, the exact solutions of the time-fractional extended (2 + 1)-dimensional Zakharov-Kuznetsov equation and a high-dimensional time-fractional KdVtype equation were obtained by means of the symmetry analysis method in [16,17]. Sing et al. suggested q-Elzakri transform method for solving multi-dimensional diffusion equations [18]. In [19], the Sumudu transform method was employed successfully for a fractional blood alcohol model. Authors in [20] presented a method based on Chebyshev polynomials to solve the fractional Bratu's equation. Alderremy et al. used the finite-difference method for the fractional two-cell cubic autocatalysis reaction model [21]. Fractional functional equations with variable order (where orders are functions of the time or space or both of them) were introduced in 1993 [22]. Due to having the memory properties, this class of equations is a powerful tool to describe mechanics of an oscillating mass subjected to a variable viscoelasticity damper and a linear spring, the motion of spherical particle sedimentation in a quiescent viscous liquid, analysis of elastoplastic indentation problems, and interpolating the behavior of systems with multiple fractional terms [22][23][24][25]. Since limited works have been done on VTFPIDEs, in the present work a numerical approach based on the bivariate Chebyshev polynomials of the fifth kind is presented to obtain approximate solutions of the PIDEs as follows: with the initial and boundary conditions where ν ij , ρ are real constants, θ , ϑ ∈ [0, 1), m ∈ Z + , and J = [0, 1] × [0, 1]. The functions h(t), φ 0 (x), φ 1 (x), ψ 0 (t), ψ 1 (t) are known continuous ones. The functions V(x, t) and f(x, t) are assumed to be sufficiently smooth, which guarantees the existence and uniqueness of the solution of Eq. (1). F is an identity or a differential operator and C 0 D σ (x,t) t designates the variable-order time-fractional derivative operator of the Caputo type. Equations of the form (1) arise in problems dealing with heat conduction in materials with memory, population dynamics, viscoelasticity, and theory of nuclear reactors (see [26,27]).
The Chebyshev polynomials of the fifth kind have been utilized for fractional ordinary differential equations with constant orders [28]. Modeling diverse phenomena in the different fields of science is performed by Eq. (1), [26,27]. The integral part of Eq. (1) possesses a singular kernel and one of its limits is a function, which makes it hard to find the exact solution of problem (1)- (2). The Chebyshev polynomials of the first to fourth kinds have been widely used to solve diverse functional equations, but the orthogonal Chebyshev polynomials of the fifth kind were taken into account less. These reasons and the importance of the equation under study motivate the authors of the current paper to present a new approach for solving the equation in (1). By generalizing the fifth-kind Chebyshev polynomials to the two-dimensional case, their efficiency as basis functions is demonstrated. Accordingly, pseudo-operational matrices of the integration and a pseudooperational matrix for the approximation of the integral part with the singular kernel are derived. Besides, the relation between the one-variable basis vector X(t) and its shifted form, X(h(t)), is exhibited as a matrix. Using the collocation method along with the resultant matrices converts the solution of the main problem to the solution of a system of algebraic equations, and solving this algebraic system leads to an approximate solution.
Authors in [31] utilized two different basis functions (Legendre and Laguerre polynomials) to construct the two-variable basis, and all computations are done for both bases, while the proposed method in the current paper constructs the two-variable basis utilizing the fifth-kind Chebyshev polynomials. It is expected that the obtained results enjoy less computational cost and more accuracy.
The rest of the paper is structured as follows: definitions of the fractional derivative and integral operators, one-and two-variable Chebyshev polynomials of the fifth kind are introduced in Sect. 2. Pseudo-operational matrices of the integration with integer and fractional orders are derived in Sect. 3, and then with the aim of them, operational matrices for the two-dimensional basis are constructed. The operational matrix of the product is constructed, the integral part in (1) is approximated, and the shifted basis vector X(h(t)) is approximated in terms of the basis vector X(t). Section 4 is dedicated to describing the suggested method, and error bounds for approximate solutions are computed in Sect. 5. The applicability and efficiency of the proposed scheme are illustrated through several experimental examples in Sect. 6. The paper concludes with a conclusion and discussion in Sect. 7.

Fractional operators
The variable-order time-fractional derivative operator of the Caputo type is defined as follows [22]: where m -1 < σ (x, t) ≤ m, m ∈ Z + and V(x, t) ∈ C m (J).
A function Y ∈ L 2 w (J) can be expanded in terms of the SCPFKs as follows: where coefficients Y j , j ≥ 0, are computed as The first few terms in (7) are practically required for approximating the function Y, i.e., where Y and X(t) are vectors as Now, BCPFKs are constructed using the one-variable SCPFKs on the domain J = [0, 1] × [0, 1] as follows: These two-variable polynomials are orthogonal with respect to the weight function where i , j are defined in (5). The two-variable function Y ∈ L 2 W (J) can be approximated as follows: where It must be mentioned that the vector (x, t) is written as (x, t) = X(x) ⊗ X(t), where ⊗ denotes the Kronecker product.

Pseudo-operational matrices
In this section, the pseudo-operational matrices of the integration (with integer and fractional orders) and the operational matrix of the product are derived. The shifted vector X(h(t)) is approximated in the vector X(t), and then a pseudo-operational matrix is computed to approximate the integral part in (1).

The integral pseudo-operational matrices
.
Proof Using series (3) and the weight function w(t), one has where B(r, s) is the well-known beta function. (9) and ε ∈ Z + , then one gets
Proof According to the definition of X(t), the left-hand side of relation (13) is calculated as follows: Now, t r is approximated in the polynomials X m (t), m = 0, 1, . . . , N , .
Substituting the last equality in the last row of (14) leads to the desired result.
is the bivariate vector in (12), then the pseudo-operational matrices of integer order regarding the variables t and x are as follows: where P ε (x) and P ε (t) are the pseudo-operational matrices in the two-dimensional case regarding x and t, respectively. (9) and RL 0 I σ (t) t is the variable-order fractional integral operator in the Riemann-Liouville sense as follows:

Theorem 3.4 Assume that X(t) is the basis vector in
where P (σ ) ε is the (N + 1)-order pseudo-operational matrix of the integration of variable order, and its entries are computed as follows: , j, i = 0, 1, . . . , N.
Proof The proof is similar to what was done for Theorem 3.2.

Remark 3.5 If (x, t) is the bivariate basis function and RL
is the variable-order Riemann-Liouville integral operator and ε ∈ Z + , then a pseudo-operational matrix of the integration with the variable order σ (x, t) is found as follows: where P (σ ,ε) t is the integral pseudo-operational matrix of variable order for the two-variable basis regarding the time variable.

The operational matrix of product Lemma 3.6 The integration of three polynomials
Proof First, the products of X j (t) and X k (t) are written as follows: where the coefficients γ (j,k) r , j, k = 0, 1, . . . , N , r = 0, 1, . . . , j + k, can be calculated from Algorithms 1 and 2. Now, the quantity q ijk is calculated as follows: By definition of the beta function, the desired result is achieved.
Suppose that X(t) is the basis vector in (9) and U is an (N + 1)-order arbitrary vector. Then one has where U is the (N + 1)-order product matrix, and its entries are computed as follows: and q ijk is given by Lemma 3.6.

Theorem 3.8 If X(t) is the one-variable basis vector, then one has
where B (β) is the (N + 1)-order pseudo-operational matrix for the integral with the weakly singular kernel, and its entries are as follows: , l, j = 0, 1, . . . , N.
Proof According to the definition of X(t) and noting that Using Lemma 3.1, t r , 0 ≤ r ≤ N , is written as So, each component of the vector in (16) is approximated as follows: Thus, one gets the following matrix representation: (2021) 2021:348 Page 11 of 26

Theorem 3.9 If h(t) ∈ C(J), then the shifted basis vector X(h(t)) is approximated in X(t) as
where h is an (N + 1)-order matrix and is found later.
Proof Noting the definition of X j (t) in (3), each X j (h(t)) is written as where coefficients ρ j,k , j, k = 0, 1, . . . , N , are given by (4). Thus, the vector X(h(t)) is written as Now, the functions h k , k = 1, 2, . . . , N , must be approximated. First, h(t) is approximated as follows: With the aid of (18), one achieves where is the (N + 1)-order operational matrix of the product, and its entries are calculated in terms of the components of the vector . Hence, the vector H(t) is approximated as follows: where e 1 = [1, 0, . . . , 0] T . So, one attains Corollary 3. 10 The following approximations are achieved by Theorems 3.2, 3.8, and 3.9:

Methodology
To find an approximate solution for problem (1)-(2), first consider the following approximation: By twice integrating (19) regarding x, the following approximations are obtained: Now, setting x = 1 in (21) leads to determining the approximate value of ∂ 3 V(0, t)/∂x∂t 2 : Twice integrating (21) regarding t leads to the following approximations: Again twice integrating (19) with respect to t yields the following approximations: Integrating (25) regarding x gives the following approximation to ∂V(x, t)/∂x: Deriving (23) regarding x and setting x = 0 lead to an approximation for ∂V(0, t)/∂x: Integrating (20) regarding t leads to the following approximation: To , t), consider the definition of the Riemann-Liouville integral operator and the approximation in (21). Using Remark 3.5, one can write where D 2 t = ∂ 2 /∂t 2 . Suppose that F(V(x, t)) in the integral part in Eq. (1) is approximated as follows: where B (θ,ϑ) = B (θ) ⊗ (B (ϑ) h ) is an (N + 1) 2 -order matrix. Substituting approximations (19)-(29) into Eq. (1) gives a residual function, and collocating it at tensor points {(x i , t j )} N i,j=0 leads to an algebraic system including (N + 1) 2 equations. x i , t j are the roots of X N+1 (x) and X N+1 (t), respectively. By solving the resultant system, the unknown coefficient vector U is determined, and an approximate solution is achieved from (23).

Error bounds
In this section, by computing error bounds of obtained approximations, an error bound for the proposed method is presented.
Set P N = Span{Z ij (t), i, j = 0, 1, . . . , N} and suppose that u N (x, t) ∈ P N is the best approximation to u(x, t) ∈ L 2 W (J), in other words, Theorem 5.1 Suppose that T N V (x, t) and V N (x, t) ∈ P N are the Taylor expansion and the best approximation to V(x, t) ∈ L 2 W (J), respectively. An approximation error can be computed as follows: where C 0 is a positive constant and Proof The error of the Taylor expansion for function V is Taking the L 2 -norm of Eq. (30) and using the definition of the beta function lead to the following inequality: Utilizing the Stirling formula in [30] for sufficiently large N leads to the following inequalities: where c i , i = 1, 2, 3, are positive constants. Using the above inequalities for (31) results in the desired result.

Theorem 5.2 Suppose that V N (x, t) and T N V (x, t) are the best approximation and Taylor expansion of V(x, t) ∈ L 2 W (J). Error bounds for derivatives of V(x, t) can be approximated as follows:
where C 1 is a positive constant and Proof According to the Taylor expansion of the function V, one has By taking the L 2 -norm of (32), one gets By the Stirling formula for the sufficiently large N , one has where c r m , m = 1, 2, 3, are positive constants. Combining the last inequalities with (33) leads to the desired result.

Theorem 5.3 Assume that σ (x, t) is a known continuous function and V(x, t) ∈ L 2 W (J). An error bound for the variable-order fractional derivative of the approximation of V(x, t) is as follows:
where C 2 is a positive constant, 0 is the same in Theorem 5.1, and σ * = max (x,t)∈J {σ (x, t)}.
Using the Stirling formula for large values of N results in the following inequalities: where c i , d i , i = 1, 2, 3, are positive constants. Combining the resultant inequalities with (35) leads to the desired results.
is its approximation obtained from the proposed method, F : R × R → R is a continuous differential operator, h(t) ∈ C(I), h 0 = max t∈I {h(t)}, and the real number > 0 exists such that

The bound of the approximation error related to the integral part in Eq.
(1) can be approximated as follows: where C 0 is a positive constant, 0 is introduced by Theorem 5.1, θ , ϑ ∈ (0, 1) and Proof According to the hypotheses of the theorem and using Theorem 5.1, one has ) .

Theorem 5.5 Suppose that V(x, t) and V N (x, t) are the exact and approximate solutions of Eq. (1) and R(x, t) is the residual function/perturbation term. Then R(x, t) → 0 when N → ∞.
Proof V N (x, t) is the approximate solution of Eq. (1), thus it satisfies the following equation: where R(x, t) is the residual function/perturbation term. Subtracting Eq. (36) from Eq. (1) leads to the following equation: By taking the L 2 -norm of Eq. (37) and using Theorems 5.1-5.4, a bound for R(x, t) can be obtained as follows: As seen, the right-hand side of (38) approaches zero if N is sufficiently large.

Numerical examples
To illustrate the efficiency and applicability of the proposed scheme, five cases of Eq. (1) are considered. Maximum absolute errors and CPU times are computed, and numerical results are compared to those reported in [31]. Computations and simulations are done by Maple 16 software.
Example 6.1 Consider the following variable-order time-fractional singular partial integro-differential equation: where (x, t) ∈ J. The initial conditions and the exact solution are as follows: Based on what was stated in Sect. 4, the following approximations are needed to calculate the residual function R(x, t).
Substituting the above approximations into Eq. (39) leads to the following residual function: The  Table 1 show that the proposed method has a reasonable computational time. Comparing results to those reported by [31] emphasizes higher accuracy and coincidence of the proposed method. The 3D plots of the exact and approximate solutions and the absolute error function are depicted in Fig. 1 for N = 4, σ (x, t) = 1e -xt , h(t) = cos(t). Absolute errors of approximate solutions at equally spaced points are listed in Table 2 and calculated for N = 4, σ i (x, t) = 1, 0.875, 0.8 + 0.005 cos(xt), i = 1, 2, 3, and h(t) = t, e t 2 . Data of this table show that the numerical results are in agreement with the exact ones.    3.00 × 10 -20 Example 6.2 Consider the following variable-order time-fractional singular partial integro-differential equation: where 1 < σ (x, t) ≤ 2, (x, t) ∈ J. The initial and boundary conditions are as follows: The exact solution is V(x, t) = x(t 3 -1) and .9203 × 10 -9 3.0399 × 10 -9 2.9478 × 10 -9 3.8210 × 10 -9 Substituting appropriate approximations into Eq. (41) leads to the following residual function: where B = B ( 1 4 ) ⊗(P 0 h ). Maximum absolute errors of obtained approximate solutions are seen in Table 3 for N = 5 and various choices of functions σ (x, t) and h(t). The results have more accuracy for h(t) = t, 1 + t 2 . For h(t) = sin(t), the errors decrease as σ (x, t) → 2. The 3D plots of the absolute error functions are depicted in Fig. 2 for N = 5, σ (x, t) = 2-0.2e -xt , and diverse cases of h(t).
The values of the absolute errors of the approximate solutions are computed at the points x j = t j = 0.1j, j = 0, 1, . . . , 10, for N = 3, σ (x, t) = 1 -0.5e -xt , 1 -cos(x)e -t , 0.25, 0.50, 1, the results are seen in Tables 4 and 5 for h(t) = t and h(t) = cos(t), respectively. Data of these tables demonstrate the agreement of the numerical results with the exact ones. The method has high accuracy even when h(t) is a continuous function. The 3D plots of the exact and approximate solutions and the absolute error function are observed in Fig. 3 for N = 3, The computational times of approximate solutions are listed in Table 6 for N = 3, h(t) = t, and various cases of σ (x, t). When values of σ (x, t) approach 1, the computational time decreases. The exact and approximate solutions are compared in Fig. 4 for N = 3, σ (x, t) = 1 -cos(x)e -t , h(t) = cos(t) at t = 0.25, 0.50, 0.75, 1. As seen, the approximate solutions are in good agreement with the exact ones.   Example 6.4 Consider the variable-order fractional partial integro-differential equation with the weakly singular kernel where 0 < σ (x, t) ≤ 1, (x, t) ∈ J. The initial and boundary conditions and the exact solution are, respectively,    Example 6.5 Consider the variable-order time-fractional partial integro-differential equation with the weakly singular kernel The initial and boundary conditions are V(x, 0) = cos(x), V(0, t) = 1 + sin(t), the exact solution is V(x, t) = cos(x) + sin(t) if σ (x, t) = 1, and f(x, t) = cos(x) + sin(t) + cos(t) + 3 2 x 2 3 sin(t).
Maximum absolute errors and CPU time are seen in Table 9 for h(t) = t, σ (x, t) = 1, and various values of N . By increasing N , the maximum absolute error decreases and CPU times show that the proposed method has a reasonable computational cost. The plots of the approximate solutions are observed in Fig. 6 for N = 5, h(t) = t, σ (x, t) = 0.8, 0.85, 0.9, 0.95, 1, and x = 1. When the values of σ (x, t) approach 1, approximate solutions get close to the exact solution in the case σ (x, t) = 1.

Conclusion
The fifth-kind Chebyshev polynomials were proposed to deal with the numerical solution of a class of partial integro-differential equations with weakly singular kernels. To this end, bivariate Chebyshev polynomials were constructed with the help of the one-variable ones, and their pseudo-operational matrices were derived by obtaining matrices for the one-dimensional case. Resultant matrices and approximations were utilized along with the collocation method to solve the problem under study. Error bounds were computed, and using them we showed that the residual function tends to zero if the number of terms in the solution series is chosen sufficiently large. The proposed method reduces the volume of computations by presenting reliable algorithms. The numerical results demonstrated the efficiency and applicability of the method since they had a good agreement with the exact ones, and the numerical results had less error in comparison with those reported in [31]. Therefore, the Chebyshev polynomials of the fifth kind are suggested to be considered as basis functions for various types of spectral and pseudo-spectral methods. By observing the stated facts, the suggested method can be employed easily to solve delay integro-partial differential equations and some nonlinear partial differential equations.