NUMERICAL SOLUTION OF VARIABLE-ORDER TIME FRACTIONAL WEAKLY SINGULAR PARTIAL INTEGRO-DIFFERENTIAL EQUATIONS WITH ERROR ESTIMATION

In this paper, we apply Legendre-Laguerre functions (LLFs) and collocation method to obtain the approximate solution of 
 variable-order time-fractional partial integro-diﬀerential equations (VO-TF-PIDEs) with the weakly singular kernel. For this purpose, we derive 
 the pseudo-operational matrices with the use of the transformation matrix. The collocation method and pseudo-operational matrices transfer the 
 problem to a system of algebraic equations. Also, the error analysis of the proposed method is given. We consider several examples to illustrate 
 the proposed method is accurate.

1 Introduction radiation of heat from semi-infinite solids [16]. Various numerical and analytical methods have been used for solving integro-differential equations. For example, Baleanu et al. [3] utilized the optimum q-homotopic analysis method for solving weakly singular kernel fractional fourth-order partial integro-differential equations, Tang [34] applied a finite difference scheme for partial integrodifferential equations with a weakly singular kernel, McLean et al. [19] used Laplace transformation of an integro-differential equation of parabolic types, Nemati et al. [25] discussed the numerical solution of two-dimensional nonlinear Volterra integral equations by the Legendre polynomials, Fairweather [12] introduced spline collocation methods for a class of hyperbolic partial integro-differential equations, Dehestani et al. [9] proposed an efficient computational approach based on the Genocchi hybrid functions for solving a class of fractional Fredholm-Volterra functional integro-differential equations, Patel et al. [26] applied two-dimensional wavelets collocation scheme for linear and nonlinear Volterra weakly singular partial integro-differential equations and many other methods have been used in this type of equation, which can refer to [4,10,25,30,31,32,35].
In 1993, Samko and Ross [6], introduced an interesting generalization of fractional operators. They presented the study of fractional integration and differentiation when the order is not a constant but a function. Afterward, the variable order fractional calculus (VOFC) is presented as a useful tool with various applications in science and engineering, specifically, VOFC describes the mechanics of an oscillating mass subjected to a variable viscoelasticity damper and a linear spring [6] to characterize the dynamics of Van der Pol equation [11], to develop motion for spherical particle sedimentation in a quiescent viscous liquid [27], to analyze elastoplastic indentation problems [14], to interpolate the behavior of systems with multiple fractional terms [33] and to obtain variableorder fractional noise [29].
Several papers have been devoted to the study of variable-order fractional problems, such as variable-order fractional differential equations [17,22,24], variable-order fractional partial integro-differential equations [2,8,20,21,23] and so on. Limited work has been done in the study on VO-TF-PIDEs with the weakly singular kernel. Therefore, we consider the VO-TF-PIDEs with the weakly singular kernel as with initial conditions and boundary conditions u(0, t) = ϕ 0 (t), u(1, t) = ϕ 1 (t), (1.2) where α, β ∈ [0, 1), which are not both zero and u(x, t) is an unknown function. The known functions h(t), f 0 (x), f 1 (x), ϕ 0 (t) and ϕ 1 (t) are defined on interval Ω = [0, 1] × [0, ∞), µ ij , λ are real constants, q is a positive integer and G is an identity or differential operator. Also, D γ(x,t) t denotes the variable order Caputo fractional derivative operator that is defined as follows [6,8]: for q − 1 < γ(x, t) ≤ q, t > 0 and q ∈ Z + . To guarantee the existence and uniqueness of the solution of the proposed equation, we supposed that the functions u(x, t) and g(x, t) be sufficiently smooth functions. To solve this problem, we focus on providing a numerical scheme to solve VO-TF-PIDEs with the weakly singular kernel by LLFs and collocation method. Then, several numerical examples are presented to illustrate the effectiveness of the proposed method.
In modeling of physical phenomena, the variable order fractional derivatives are very powerful, because they have the memory with respect to time and spatial location. Therefore, the fractional derivatives statements in the proposed equation are replaced by the variable-order fractional derivative. Solving the problems by the variable-order fractional derivative with the help of an analytical method is not convenient. Hence, this fact motivates us to present the numerical algorithm base on Legendre-Laguerre functions. The Laguerre polynomials are orthogonal in the semi-infinite intervals which are appropriate for approximating problems of natural phenomena in semi-infinite intervals.
This paper is structured as follows. In Section 2, we introduce the formulation of LLFs and their properties. Section 3 is devoted to the pseudooperational matrices of integration and variable-order fractional derivative for LLFs by use of the transformation matrix of Legendre and Laguerre polynomials to Taylor polynomials. In Section 5, we will describe a successful numerical approach, which is used for making up the solution. The convergence of the approximate solution is given in Section 6. In Section 7, we report our numerical results to demonstrate the accuracy of the proposed scheme. Also, conclusions are given in Section 8.

Legendre-Laguerre functions
In this section, we introduce two-dimensional LLFs, which obtain these functions with Legendre and Laguerre polynomials. Consider the shifted Legendre polynomials P m (x) on the interval [0, 1] as [7] The shifted Legendre polynomials are orthogonal with respect to the weight function w(x) = 1 in the interval [0, 1] with the orthogonal property where δ mn is the Kronecker function. The Laguerre polynomials L n (t) on the interval [0, ∞) are defined as follows [7]: The Laguerre polynomials are orthogonal with respect to the weight function Accordingly, we construct two-dimensional LLFs as follows: ∈ Ω, m = 0, 1, . . . , M, n = 0, 1, . . . , N.
Two-dimensional LLFs ψ mn (x, t), are the orthogonal basis with respect to the weight function w(x, t) = e −t . Using the orthogonality property of Legendre and Laguerre polynomials, we obtain A function f (x, t), which is integrable in Ω, can be expanded as Then, we have truncated series for f as where F = f mn (M+1)×(N +1) , m = 0, 1, . . . , M, n = 0, 1, . . . , N,

Pseudo-operational matrices of integration for LLFs
In this section, we calculate the integral pseudo-operational matrix without coefficients and with the weakly singular coefficients for Legendre and Laguerre polynomials.

Pseudo-operational matrices of integration
To calculate the integral pseudo-operational matrices of Legendre polynomials, we use the Taylor polynomials T i (x) = x i , i = 0, 1, · · · , M [13]. The following relation holds among these polynomials and Legendre polynomials: otherwise, i, j = 0, 1, . . . , M, D 1 denotes the transformation matrix of the Legendre polynomials to the Taylor polynomials. Now, by integrating from P (x), we obtain is an integral pseudo-operational matrix of Legendre polynomials and Moreover, to deal with the problem, we utilize the following integral Moreover, we can write L(t) in the matrix form as follows: where we introduce Q 2 as an integral pseudo-operational matrix of Laguerre polynomials. In addition, we need

Pseudo-operational matrix with the weakly singular coefficients
To deal with an integral part of Equation (1.1), we present the pseudo-operational matrix with the weakly singular coefficients. Hence, we obtain integral of Legendre polynomials with the weakly singular kernel is the pseudo-operational matrix with the weakly singular coefficients for Legendre polynomials and Also, we obtain integration of Laguerre polynomials with the weakly singular kernel 1/(t − η) β on interval [0, h(t)] as: , the pseudo-operational matrix with the weakly singular coefficients for Laguerre polynomials and 4 Pseudo-operational matrix of the variable-order fractional derivative In this section, we derive the pseudo-operational matrix of the variable-order fractional derivative of Laguerre polynomials by using some properties of Caputo variable-order fractional derivative.
Theorem 1. The pseudo-operational matrix of the variable-order fractional derivative of order q − 1 < γ(x, t) ≤ q for the Laguerre polynomials is given by is called the pseudo-operational matrix of variable-order fractional derivative for the Laguerre polynomials.
Proof. By using some properties of Caputo variable-order fractional derivative and Equation (2.1), we obtain . By employing the above equations, we get The matrix form of Equation (4.1) can be written as follows:

⊓ ⊔
Also, according to the above process, we obtain is calculated.

Description of the method
In this section, we present the numerical scheme for solving the problem given in Equation (1.1) with conditions in Equation (1.2). According to the proposed equation, we expand the following function by the LLFs as: To approximate the other functions, we use the pseudo-operational matrices that presented in previous sections. By integrating the above equation with respect to t and substituting initial conditions into it, we obtain Now by integrating Equation (5.3) of order 2 with respect to x, we get is unknown, in order to obtain this function, integrating Equation (5.4) with respect to x from 0 to 1 Then, Also, we need to calculate the following expression. By integrating Equation (5.1) of order 2 with respect to x and using the boundary conditions, we have In Equations (5.6) and (5.7), ∂ 3 u(0,t) ∂x∂t 2 is an unknown function. By integrating Equation (5.6) from 0 to 1 with respect to x, we get By integrating Equation (5.7) with respect to t, we have . In addition, we calculate the variable-order fractional derivatives by applying the pseudo-operational matrix presented in Section 4. Therefore, we achieve according to the process in Theorem 1 is obtained. Consequently, by using the pseudo-operational matrix and following integration formula we approximate an integral part of Equation (1.1). In specific example dη. Q 1 andQ 2 similar toQ 1 andQ 2 are calculated, respectively. So that, We obtain an algebraic equation by substituting the above approximate functions in Equation (1.1). Then, we collocated this equation in nodal points of Newton-Cotes [7]. Ultimately, we get an unknown matrix U by solving a system of algebraic equations and using Newton's iterative method. Finally, by substituting U in Equation (5.5), we achieve the approximate solution to the problem.

Convergence and error estimation
In this section, we present the upper bound of errors by the following theorem in Sobolev space [5]. For this goal, the Sobolev norm of integer order µ ≥ 0 on the interval (a, b) in R is defined as The seminorm is defined as [5] |f | H µ;M (a,b) = µ j=min(µ,M+1) where C depends on µ.
Also, Laguerre approximations in weighted Sobolev spaces on the half-line R + of integer order µ ≥ 0 is defined as follows: Lemma 2. Let g ∈ H µ (R + ), such that P N g = N n=0 g n L n is the best approximation of g then for any µ ≥ 0 and 0 ≤ r ≤ µ, we have [5] g − P N g H k (R+) ≤ CN r− µ 2 g H µ w,µ (R+) , where C depends on µ.
Theorem 2. Suppose that u M and u N are approximated with truncated Legendre-Laguerre series u Mn and u mN , respectively. So that, Then, the following estimates hold Proof. According to the hypothesizes and Lemmas 1 and 2, we have 1] . And also, . Hence, we obtain the upper bound of the errors. ⊓ ⊔ Corollary 1. Suppose u N , u M ∈ H µ (Ω) and µ ≥ 1 then by setting N, M ≥ µ−1, and considering the previous theorem, we deduce In the first and second inequality, the values of M and N are constant, respectively. According to this fact, in the first inequality, the rate of convergence of u Mn to u M is faster than 1 N to the power of N 2 + 1 2 − r. And also, in the second inequality, the rate of convergence of u mN to u N is faster than 1 M to the power of M − 2k + 3 2 .

Illustrative examples
In order to test the validity of the present method, we consider some examples and apply the technique in Section 5 for solving them. Then, the numerical results are compared with their exact solution and other methods. The computations were performed on a personal computer and the codes were written in MATLAB 2016.

Example 1. Consider VO-TF-PIDEs with weakly singular kernel [30]
The exact solution of this example is u(x, t) = xt 2 . To evaluate the precision of the proposed method, we present the norm of residual error . In Table 1, the absolute errors are obtained between the approximate solutions and the exact solution using Legendre-Laguerre functions for various functions of h(t) and γ(x, t) with M = N = 1 (here γ 3 (x, t) = 0.8 + 0.005 cos(xt)). Also, maximum absolute errors with various functions of h(t) and γ(x, t) for x = 1 and t ∈ [0, 20] with M = N = 1 and CPU time (in seconds) are shown in Table 2. Due to the errors in figure [30], the proposed method is more accurate in comparison to the method in [30] with h(t) = t, N = 2 and γ(x, t) = 1.   1 < γ(x, t) ≤ 2, with initial conditions u(x, 0) = −x, ∂u(x,0) ∂t = 0, 0 ≤ x ≤ 1, and boundary conditions u(0, t) = 0, u(1, t) = t 3 − 1, t > 0, with g(x, t) =

Conclusions
This work introduces a new numerical approach for solving VO-TF-PIDEs with the weakly singular kernel. The method is based on expanding the derivation of the solution in terms of the LLFs and using the appropriate collocation points. The main objective is to use an integral pseudo-operational matrix of LLFs, which provides good conditions to receive the approximate solution with high accuracy. The advantage of the proposed method is that, despite defined time in an infinite interval, we achieved a good approximate solution. Also, satisfactory examples illustrate the validity and efficiency of the method.