New Numerical Algorithm to Solve Variable-Order Fractional Integrodifferential Equations in the Sense of Hilfer-

In this article, a numerical technique based on the Chebyshev cardinal functions (CCFs) and the Lagrange multiplier technique for the numerical approximation of the variable-order fractional integrodifferential equations are shown. The variable-order fractional derivative is considered in the sense of regularized Hilfer-Prabhakar and Hilfer-Prabhakar fractional derivatives. To solve the problem, first, we obtain the operational matrix of the regularized Hilfer-Prabhakar and Hilfer-Prabhakar fractional derivatives of CCFs. Then, this matrix and collocation method are used to reduce the solution of the nonlinear coupled variable-order fractional integrodifferential equations to a system of algebraic equations which is technically simpler for handling. Convergence and error analysis are examined. Finally, some examples are given to test the proposed numerical method to illustrate its accuracy and efficiency.


Introduction
Fractional calculus as a branch of mathematical analysis due to important applications in computational science, engineering [1][2][3][4], chemistry, and many physical equations such as Burger's equation [5] and Korteweg-de Vries equation [6], and viscous fluid [7] has been broadly significant. Recently, obtaining solution of the linear and nonlinear differential equations has important roles in fractional calculus [8,9], because the behavior of the equations is determined by their analytical solutions. In this paper, a differential equation including a variable-order fractional operator is considered in which the equation is an extension of the fractional differential equations of fractional order. Getting exact solutions of the equations in mathematical sciences and physics is not easy, so approximate solutions to these equations using numerical methods are obtained. Different numerical methods for obtaining solutions of the linear and nonlinear equations in the papers are presented. For instance, in [10], a solution of the nonlinear variable-order fractional Riccati differential equations was obtained by using the spectral technique. A finite difference method for getting approximate solution of the variable-order fractional Volterra integrodifferential equations (VO-FVIDEs) has been considered in [11]. The numerical solution of the variable-order fractional integrodifferential equation consisting of cubic spline interpolation can be obtained in [12]. In paper [13], a collocation technique to obtain solution of the variable-order nonlinear cable equation was found. Obtaining the numerical solution of a class of nonlinear variable-order fractional differential equations (FDEs) by using the Legendre wavelets in [14] is studied. Also, Bhrawy and Zaky [15] used the shifted Jacobi polynomials to obtain the solution for variable-order fractional Schrödinger equations, Bhrawy and Zaky [16] studied the Jacobi-Gauss-Lobatto collocation method to obtain the solution for variable-order fractional Schrödinger equations, Mahmoud et al. [17] proposed the Jacobi wavelet collocation method to obtain the solution for variable-order fractional equations, Zaky et al. [18] proposed the shifted Chebyshev polynomials to obtain the solution for variable-order fractional equations, and in [19], a proper discrete form of fractional Gronwalltype inequality is introduced and other methods such as the shifted Jacobi collocation method [13] and shifted Jacobi method [20]. Other methods are proposed to obtain numerical solutions of the variable-order fractional integrodifferential equations in the literatures [21][22][23][24][25][26][27]. In this paper, we obtain the solution for the following fractional differential equation of order 0 < μðtÞ, ν ≤ 1 by using of the Chebyshev cardinal functions (CCFs) and the Lagrange multiplier method in which f ðx, tÞ is known and uðx, tÞ is unknown: The initial condition for this equation is defined by Here, C D γ,μðtÞ,ν ρ,ω,0 + and D γ,μðtÞ,ν ρ,ω,0 + are regularized Hilfer-Prabhakar fractional derivative [28] and Hilfer-Prabhakar fractional derivative [28] of order 0 < μðtÞ ≤ 1, 0 < ν ≤ 1, respectively. According to the definition of the Hilfer fractional derivative, this Hilfer-Prabhakar fractional derivative is defined but instead of Riemann integral operator in the Hilfer fractional derivative, Prabhakar integral operator [28] of variable-order fractional is replaced, in which this Prabhakar integral operator for 0 < μðtÞ ≤ 1 is given by where E γ ρ,μðtÞ ðωt ρ Þ is the Prabhakar function and it is a generalization of the two-parametric Mittag-Leffler function E α,β ðtÞ, α > 0, β ∈ ℂ, for γ = 1 and the classical Mittag-Leffler function E α ðtÞ, α > 0, for γ = μðtÞ = 1 is considered which is defined by [28] where ðγÞ n = ðΓðγ + nÞÞ/ΓðγÞ and Γð:Þ are the gamma function; moreover, this Prabhakar function for RðμðtÞÞ > 0 is an entire function [29]. Our interest in this type of fractional derivative is related to its application in physical phenomena, for example, linear viscoelastic con-stitutive equations [30], anomalous dielectric relaxation of Havriliak-Negami function [31,32], fractional viscoelasticity [33], spherical stellar systems, stochastic processes [34], telegraph equations [34], and anomalous relaxation in dielectrics [35]. Since this mathematical system which is given in Equation (1), due to the variableorder fractional operators and nonlinearity, is very complex, we need to reduce it by using a highly accurate and efficient expansion scheme. Thus, for this aim, we first extract an operational matrix of variable-order fractional derivatives for the CCFs; then, we use them for extending the unknown solution. In other words, the problem shown in Equation (1) is converted into an algebraic system of equations by exploiting the operational matrix of variable-order derivative. In this paper, we consider a class of variable-order fractional integrodifferential equations. For solving the given equations, operational matrices based on the CCFs are applied. First, we approximate the unknown function and its derivatives in terms of the CCFs. Then, by substituting these approximations into the equation and applying the properties of the CCFs together with the collocation points, the main problem is reduced to a set of nonlinear algebraic equations. By solving this system, the approximate solution is calculated. This article is divided into four sections as follows. In Section 2, some important and application definitions, lemma, and theorem which are applied in the next section are expressed. In Section 3, properties of the Chebyshev cardinal functions (CCFs) and the Lagrange multipliers are introduced. In Section 4, we use the properties described in Section 3 to solve the variable-order fractional differential equations. In Section 5, to illustrate the efficiency and accuracy of the proposed method in this article, four examples are displayed.

Definition and Fundamental Properties of Fractional Calculus
This section refers to the definitions and lemmas which are used for the next sections, and these definitions and lemmas are taken from articles [36,37].

Approximation Function and Properties of the Lagrange and the CCF Polynomials
Any function uðx, tÞ ∈ Cð½x min , x max × ½0, t max Þ can be approximated in terms of Chebyshev cardinal functions of order m and Lagrange interpolating polynomials of order n as follows: where φ i ðxÞ is the Chebyshev cardinal functions of order m on the interval x ∈ ½x min , x max which is defined by [38] φ and ψ j ðtÞ is the Lagrange interpolating polynomials of order n on the interval ½0, t max which is given by [38] where Also, in relation (10), t j is considered Function uðx, tÞ in Equation (10) can be rewritten as follows: where Y m ðxÞ, Δ, Λ n ðtÞ is defined by The purpose of this paper is to express the function ψ j ðtÞ, j = 0, 1, 2, ⋯, n + 1, as an operational matrix as follows [39]: where γ j , b j,k for k = 0, 1, 2, ⋯, n are given by 3 Abstract and Applied Analysis Theorem 4. Suppose ψ j ðtÞ is relation (17). Then, the Hilfer-Prabhakar derivative of a regularized type and the Hilfer-Prabhakar derivative on the function ψ j ðtÞ are obtained, respectively, as follows:     Abstract and Applied Analysis Proof. By taking the Hilfer-Prabhakar derivative of a regularized type on both sides of Equation (17), we get applying Equation (9) on Equation (21) and the result is achieved. Prove relation (20) like the same process is considered for the Hilfer-Prabhakar derivative of a regularized type. Now, using Theorem 4, the function Hilfer-Prabhakar derivative of a regularized-type ψ j ðtÞ can be approximated as a finite series in the form of the Chebyshev cardinal functions as follows: Also, using Theorem 4, the function Hilfer-Prabhakar derivative ψ j ðtÞ can be approximated as a finite series in the form of the Chebyshev cardinal functions as follows: j,r ψ r t ð Þ, j = 1, 2, ⋯, n + 1: ð23Þ So according to relations (22) and (23), we can get the following results:

A Description of the Proposed Method for Solving the System
In this section, using relations described in Section 3, we solve the proposed system stated in Section 1 which is introduced with (1). Using relation (15) and applying relations (24) and (25) on Equation (1), in this case, we obtain

Abstract and Applied Analysis
Using relation (15) for the function f ðx, tÞ, in this case, the following relation is obtained: where M = ½ f ðx i , t j , uðx i , t j ÞÞ ðm+1Þ×ðn+1Þ and W = ðΔðE μðtÞ + E μðtÞ Þ − MÞ and the initial condition for this equation is defined by Here, dðxÞ = ∑ m+1 i=1 dðx i Þφ i ðxÞ ≜ d T Y m ðxÞ, d = ½dðx 1 Þ, dðx 2 Þ, ⋯, dðx m+1 Þ T , and Ξ = ðΔΛ n ð0Þ − dÞ are considered. By solving the following equations, the unknown function uðx, tÞ is determined: Here, E 1 and E 2 show the first and the second values of the maximum absolute error (MAE) provided by the proposed method, respectively. Moreover, ðm i + 1Þðm i + 1Þ for i = 1, 2 are the number of the basis functions used in the first and the second implementations, respectively. Here, the absolute error which is the difference between a numerical and exact solution is defined by Example 6. Let μðtÞ = ð6/10Þ − ð3/10Þ sin ðπtÞ. Then, we consider the following linear fractional system of orders μðtÞ: where f ðx, t, uðx, tÞÞ is given by The exact solution is uðx, tÞ = e −ðx+tÞ . The upper part in Figure 1 illustrates the numerical solution of Example 6 with m = 4 and n = 10, and the lower part shows the absolute error of this example. Table 1 shows the absolute error values for different values of m, n. Also, The values of the C-order of the presented method are listed in Table 2. From the obtained (3.11, 0.45) 3:3110e − 9 1 :9981e − 9 1 :8032e − 9 4 :9952e − 10 Table 6: C-order of the presented method for Example 9.  Table 1, we see that the proposed method is an effective and good tool for solving this problem.
Example 7. Let μðtÞ = ðsin ð10πtÞÞ/6 + ð2/3Þ. Then, we state the following nonlinear fractional system: where f ðx, t, uðx, tÞÞ is given by The analytical solution for this question is uðx, tÞ = x 3 sin ðtÞ. The upper part in Figure 2 illustrates the numerical solution of Example 7 with m = 4 and n = 10, and the lower part shows the absolute error of this example. Table 3 shows the absolute error values for different values of m, n. Example 8. We consider a fractional system of orders μðtÞ = t 2 − t + 0:8 as follows: where f ðx, tÞ is given as where for this problem, analytical solution is uðx, tÞ = t 2 ð2x − x 2 Þ.The upper part in Figure 3 illustrates the numerical solution of Example 8 with m = 4, n = 10, and the lower part shows the absolute error of this example. Table 4 shows the absolute error values for different values of m, n.
Example 9. We consider a fractional system as follows: where f ðx, tÞ is given as where for this problem, analytical solution is uðx, tÞ = 2 l t l sin ðx + 1Þ. The upper part in Figure 4 illustrates the numerical solution of Example 9 with m = 4 and n = 10, and the lower part shows the absolute error of this example. Table 5 shows the absolute error values for different values of m, n. Also, The values of the C-order of the presented method are listed in Table 6.

Data Availability
No data were used to support the study.

Conflicts of Interest
The author declares that there are no competing interests.