Operational matrices based on the shifted fifth-kind Chebyshev polynomials for solving nonlinear variable order integro-differential equations

In this research, we study a general class of variable order integro-differential equations (VO-IDEs). We propose a numerical scheme based on the shifted fifth-kind Chebyshev polynomials (SFKCPs). First, in this scheme, we expand the unknown function and its derivatives in terms of the SFKCPs. To carry out the proposed scheme, we calculate the operational matrices depending on the SFKCPs to find an approximate solution of the original problem. These matrices, together with the collocation points, are used to transform the original problem to form a system of linear or nonlinear algebraic equations. We discuss the convergence of the method and then give an estimation of the error. We end by solving numerical tests, which show the high accuracy of our results.

Since the order of fractional integrals and derivatives may take any arbitrary value, a new extension of these operators has been proposed such that the order of these operators is not a constant but a function of some independent variables such as time or space. In 1993, Samko and Ross [16] were the first researchers who have suggested the study of VO operators. Then theory-based studies of VO calculus have been more deeply investigated by Lorenzo and Hartley [17]. Soon after, many definitions of VO derivative operators have been introduced by some researchers such as Riemann-Liouville (RL) [18,19], Lorenzo-Hartley [17], Coimbra [20], and Caputo [2,21] derivatives. These operators have been used to describe some models in a variety of science fields including biochemical tumorous bone remodeling models [22], characterizing the dynamics of van der Pol oscillators [23]; see also [24,25]. Since in this type of problems, we confront with a kernel of VO [26], computing analytical solutions is very difficult. Hence developing effective numerical techniques for finding approximate solution for such problems is very important and necessary. In recent years, many researchers have proposed different schemes to solve this kind of problems. To mention a few, we refer to [27][28][29][30], where the authors have applied operational matrices based on various polynomials to get approximate solutions of different problems of VO.
A significant aim of this research is to express a numerical scheme to solve the following VO-IDEs: with initial conditions where p -1 < υ(t) ≤ p, p is a positive integer number, F is a given continuous function, λ and a i , i = 0, 1, . . . , p -1, are real constants, K 1 , K 2 , φ 1 , and φ 2 are given known functions, z(t) is the unknown solution, and C D υ(t) t denotes the variable-order derivative operator in the Caputo sense.
Many researchers in various fields of science employ orthogonal basis functions to get approximate solutions for many problems [31][32][33]. The fifth-kind Chebyshev polynomials consist a special class of symmetric orthogonal polynomials, which are created with the help of the extended Sturm-Liouville theorem for symmetric functions. In this work, with the help of these polynomials, we reduce problem (1)- (2) to the solution of a system of nonlinear algebraic equations, which greatly simplifies the problem under study.
The design of this research is as follows. In Sect. 2, we introduce some essential definitions of variable fractional calculus and some basic properties of the SFKCPs. Section 3 is devoted to proposing a numerical scheme to solve problem (1)-(2). In Sect. 4, we study an error bound of the proposed scheme. Section 5 includes some examples. In the end, we give concluding remarks in Sect. 6.

Perliminaries
In this section, we present the definitions of VO RL-integral and Caputo derivative. Then, some basic properties of the SFKCPs are given which are used later.

VO fractional calculus
Two main properties of these operators are given as follows:

Definition of the SFKCPs and function approximation
The SFKCPs on the interval [0, 1] are defined by [28,35] where C m (t) is the fifth-kind Chebyshev polynomial defined on [-1, 1] as follows: Furthermore, the analytic form of the SFKCPs of degree m is given by and Also, the orthogonality condition is given for these polynomials as follows: Lemma 2.1 (See [35]) The SFKCPs satisfy the following boundedness property on [0, 1] for all s ≥ 0: Suppose that r 1 , r 2 ∈ L 2 w * (0, 1). Then the inner product and norm in L 2 w * (0, 1) are, respectively, defined by Any arbitrary function z(t) ∈ L 2 w * (0, 1) can be expanded by the SFKCPs as By considering only the first M + 1 terms in (5), we can approximate z(t) as where and in the vector Z = [z 0 , z 1 , . . . , z M ] T , the entries z i , i = 0, 1, . . . , M, are given by In a similar way, a bivariate function f (t, τ ) ∈ L 2 w * ((0, 1) × (0, 1)) can be approximated based on the SFKCPs as where F is an (M + 1) × (M + 1) matrix given by We can consider the vector ϕ(t) in a matrix form as where ς i,j are given by (4), and be its expansion using the SFKCPs. Then, for i > 3, the coefficient z i is bounded as (7). By applying the first-order derivative on this vector we get

Lemma 2.2 Consider the basis vector ϕ(t) defined by
where D is the operational matrix of derivative based on the SFKCPs given by Also, for m ≥ 2, we can write Proof It can be easily proved in a similar way as that of the corresponding theorem in [36]. (7), the dual operational matrix Q is given by

Lemma 2.3 For the vector ϕ(t) given by
where Q = AHA T with the well-known Hilbert matrix H.
Proof The proof process is similar to that given in [36]. (7) can be approximated as

Lemma 2.4 The integral of the vector ϕ(t) given by
where P is called the operational matrix of integration for the SFKCPs.

Theorem 2.2 Let ϕ(t) be the SFKCPs vector given in
where and otherwise.

Numerical scheme
The aim of this section is to propose a numerical scheme for solving problem (1)- (2). To do this, we first consider an approximate solution of equation (1) in terms of the SFKCPs as By employing C D υ(t) t to both sides of (15) and using (12), we have Now we must approximate the Fredholm and Volterra parts of equation (1). To do this, the functions K 1 , K 2 , φ 1 , and φ 2 are expanded using the SFKCPs as From (9)- (11) and (17), we obtain Substituting (15), (16), (18), and (19) into equation (1) yields Taking (8) and (15) into account, we can rewrite the initial conditions (2) as follows: On the other hand, by introducing the approximation z(t) Z T ϕ(t) into the functions φ 1 and φ 1 given by (17), we get To calculate the approximate solution, we put the collocation points r M+2 for r = 1, . . . , M + 1p into equation (20). By solving simultaneously the resulting system and system (21), we get an approximation of the solution using (15).

Convergence analysis
Here we consider the convergence of the approximate solution obtained by the proposed scheme in Sect. 3 to the analytical solution of problem (1)-(2).
is the universal error, then E n (t) can be evaluated as Proof By applying RL I υ(t) t to both sides of equation (1), we can rewrite equation (1) as follows: So z M (t) satisfies the following equation: where R M (t) is the residual function given by Then we have Using Theorem 4.1, we have On the other hand, since p -1 < υ(t) ≤ p, we have Since F, φ 1 , and φ 2 satisfy the Lipschitz conditions, we can write where k 1 = max (t,τ )∈(0,1) 2 |K 1 (t, τ )| and k 2 = max (t,τ )∈(0,1) 2 |K 2 (t, τ )|. Substituting (24)- (26) into (1) yields Therefore it is clear that R M (t) tends to zero as M → ∞.

Numerical examples
Now we apply the proposed scheme to some examples. For solving these examples, we used the Mathematica software.
Example 5.1 Consider the following VO problem:   in which (·, ·) is the incomplete gamma function. We have solved this problem by different values of M for υ(t) = sin 2 (t) + 2, υ(t) = t 2 + 2, and the analytical solution z(t) = e t . Figure 1 and Table 1 display the numerical results. As it can be seen from these results, the approximate solution obtained by the proposed scheme converges to the analytical one by increasing the number of basis functions.
Example 5.2 Consider the following VO problem [37]: where t ∈ [0, 1]. By considering υ(t) = t and carrying out the proposed scheme, the outputs obtained for this problem are depicted together with the analytical solution (z(t) = t 19 4 + t 31 5 ) in Fig. 2. From Fig. 2 it is clear that increasing the number of basis functions improves the accuracy. Furthermore, in Table 2, we have compared the outputs obtained by the proposed scheme with the method of [37] based on the Bernstein polynomials.   [38,39]: which gives the analytical solution. As it is seen, the proposed scheme gives the analytical solution with M = 2 (only three basis functions) compared to the methods introduced in [38][39][40]. Table 3 reports the maximum absolute errors (MAE) (E ∞ (M)) obtained in [38][39][40]. 1.42e-3 3.34e-9 1.06e-1 16 3.47e-4 1.08e-10 5.30e-2

Conclusion
In this research, we have generalized a collocation method including the shifted fifth-kind Chebyshev polynomials to numerically solve variable order integro-differential equations in the Caputo sense. For finding approximate solutions of the considered equations, we have used the properties of the shifted fifth-kind Chebyshev polynomials. In addition, by applying the collocation points, we have changed the primary problem to solving a system of algebraic equations to get an approximate solution. Also, we have discussed the convergence of the numerical solution obtained by the proposed scheme. Eventually, the efficiency and suitability of the proposed scheme are displayed by solving some problems of variable order.