Designing a Matrix Collocation Method for Fractional Delay Integro-Differential Equations with Weakly Singular Kernels Based on Vieta–Fibonacci Polynomials

In the present work, the numerical solution of fractional delay integro-differential equations (FDIDEs) with weakly singular kernels is addressed by designing a Vieta–Fibonacci collocation method. These equations play immense roles in scientific fields, such as astrophysics, economy, control, biology, and electro-dynamics. The emerged fractional derivative is in the Caputo sense. By resultant operational matrices related to the Vieta–Fibonacci polynomials (VFPs) for the first time accompanied by the collocation method, the problem taken into consideration is converted into a system of algebraic equations, the solving of which leads to an approximate solution to the main problem. The existence and uniqueness of the solution of this category of fractional delay singular integro-differential equations (FDSIDEs) are investigated and proved using Krasnoselskii’s fixed-point theorem. A new formula for extracting the VFPs and their derivatives is given, and the orthogonality of the derivatives of VFPs is easily proved via it. An error bound of the residual function is estimated in a Vieta–Fibonacci-weighted Sobolev space, which shows that by properly choosing the number of terms of the series solution, the approximation error tends to zero. Ultimately, the designed algorithm is examined on four FDIDEs, whose results display the simple implementation and accuracy of the proposed scheme, compared to ones obtained from previous methods. Furthermore, the orthogonality of the VFPs leads to having sparse operational matrices, which makes the execution of the presented method easy.


Introduction
Numerous works have been devoted to numerically solving fractional integro-differential equations with weakly singular kernels. These equations emerged in diverse fields of science, such as the heat conduction problem, radiative equilibrium, elasticity, and fracture mechanics [1][2][3][4]. Determining the analytic solutions of fractional integro-differential equations with weakly singular kernels is often complicated and even infeasible. Therefore, finding new numerical methods or developing existing methods to solve this class of equations is unavoidable. Hence, different computing methods have been presented, for where I 0 = [0, 1], 0 < q, r ≤ 1, g, p, f : I 0 → R and θ i : I 0 × I 0 → R, i = 1, 2, 3 are known continuous functions, and C 0 D η t is the Caputo derivative operator of the order η ∈ (0, 1]. The initial condition is as follows: where z 0 is a real constant. Equation (1) is revealed in astrophysics, economy, control, biology, and electro-dynamics [10,11]. To solve problem (1)- (2), the operational matrices of the integration of the integer and fractional orders and the product operational matrix are derived for the VFPs. A matrix relationship between the main basis and its delay form is presented. The integral parts with weakly singular kernels are approximated in a matrix form. By substituting resultant approximations and operational matrices into the considered equation, the main problem is transformed into an algebraic equation; collocating it at the collocation points leads to a system of algebraic equations. It is worth noting that the roots of the VFPs are considered collocation nodes. By solving the resultant system, an approximate solution is obtained for the main problem. The existence and uniqueness of the solution of problem (1)-(2) is investigated, using Krasnoselskii's fixed point theorem. An error bound is calculated for the residual function in a Vieta-Fibonacci-weighted Sobolev space. For this, it is necessary to prove the orthogonality of the derivatives of the VFPs. As seen, Equation (1) has integral parts with singular kernels and some delay terms, which make it difficult to solve the problem (1)- (2). Hence, there exist few numerical and analytical approximation methods to solve Equation (1). Rezabeyk et al. [14] constructed the fractional-order Euler polynomials to obtain the numerical solution of FDIDEs in the form of Equation (1). Ezz-Eldien and Doha proposed a Chebyshev collocation method to solve pantograph Volterra integro-differential equations [10]. Zhao et al. suggested the Sinc method for pantograph Volterra delay integro-differential equations [11]. The discretization process of these three methods increased the computational cost. The VFPs are not widely taken into consideration as basis functions in the spectral methods. These reasons motivate the authors to employ the VFPs to proffer a numerical method with a smaller computational size to solve the FDIDEs with weakly singular kernels. In general, the objectives of this work are highlighted as follows: Eventually, the obtained approximate solutions demonstrate the efficiency of the VFPs as an influential tool to solve different functional equations.
The rest of the paper is organized as follows: some vital definitions and concepts related to the fractional calculus are recalled in Section 2. The conditions of the existence and uniqueness of the solutions of problems (1)-(2) are studied in Section 3. The shifted VFPs are introduced, and their operational matrices are derived in Section 4. The method of extracting the VFPs and their derivatives from Rodigues' formula is expressed in this section. Section 5 is devoted to stating the methodology. An error bound for the residual function is estimated in a Vieta-Fibonacci-weighted Sobolev space in Section 6. The proposed scheme is implemented on four illustrated examples in Section 7. The paper concludes with a conclusion in Section 8.

Fractional Operators
Definitions of the Caputo derivative and the Riemann-Liouville integral and some of their properties are given in this section. Definition 1. The fractional Caputo derivative of order µ > 0 of a n times differentiable function h is given as follows [24]: The operator C 0 D µ t possesses the following properties: where γ, λ 1 , and λ 2 are constants.

Definition 2.
The fractional Riemann-Liouville integral of order µ > 0 of a continuous function h is defined as shown below [24]: The operator RL 0 I µ t has the following characteristics:

Existence and Uniqueness of the Solution
In this section, conditions of the existence and uniqueness of the solution for the problem (1)-(2) are investigated and proved using Krasnoselskii's fixed point theorem. First, this theorem is stated.

Theorem 1.
Let Ω be a closed convex and non-empty subset of a Banach space (X, ||.||) (Krasnoselskii's fixed-point theorem [33]). Suppose that A, B maps Ω into X such that the following hypotheses are fulfilled: Ax + By ∈ Ω for all x, y ∈ Ω; (ii) A is a contraction mapping; (iii) B is compact and continuous.
Then there exists z ∈ Ω such that z = Az + Bz.
By applying the Riemann-Liouville integral operator of the order η ∈ (0, 1] on Equation (1), the following equivalent equation is achieved: Now suppose that C(I 0 , X) is a Banach space of real-valued continuous functions from I 0 = [0, 1] into X ⊆ R equipped with the norm X C = sup t∈I 0 |X (t)|, ∀X ∈ C(I 0 , X) and the following assumptions hold for any t ∈ I 0 and (t, ξ) ∈ J = I 0 × I 0 : Theorem 2. Suppose that the assumptions in (4) and the following inequalities hold: Then, problem (1)-(2) has a unique solution on C(I 0 , X).
Proof. Let B r = {u ∈ C(I 0 , X) | u C ≤ r} subject to the following: Then B r is a closed, bounded, and convex subset of C(I 0 , X). The operators K and L are defined as shown below: It is shown that K + L has a fixed point in B r . The proof is divided into four steps: Step 1 It is shown that KZ + LU ∈ B r for every Z, U ∈ B r . Using (4) and (6), one obtains the following: Therefore, KZ + LU C ≤ r, which implies that KZ + LU ∈ B r for any Z, U ∈ B r .
Step 2 The operator K is a contraction mapping on B r . For each Z, U ∈ B r and each t ∈ I 0 , one gets the following: From (5), this shows that K is a contraction mapping on B r .
Step 3 The operator L is compact and continuous. We know that the following holds: This shows that L is uniformly bounded on B r . It remains to prove the compactness of the operator L. For t 1 , t 2 ∈ I 0 such that t 1 < t 2 and Z ∈ B r , one has the following: By taking norm, one has the following: As t 2 → t 1 , the right-hand side of the above inequality tends to zero. Thus, L is equicontinuous. As a consequence of Steps 1-3 together with the Arzela-Ascoli theorem, one deduces that the operator L is compact and continuous. Finally, by Krasnoselskii's fixed point theorem, problem (1) has at least one solution on I 0 .
Step 4 For any Z, U ∈ C(I 0 , X) and t ∈ I 0 one has the following: Imposing the hypothesis in (5) implies that W is a contraction mapping. It follows that W has a unique fixed-point, which is a solution of the problem (1)-(2).

Shifted VFPs and Their Operational Matrices
The definition of the shifted VFPs is given in this section, the operational matrices of the integration of the integer and fractional orders and the product operational matrix are derived. A matrix relationship between the main basis and its delay form is presented, and the integral parts with the singular kernels are approximated, using an operational matrix.
These polynomials are orthogonal concerning the weight function On the other hand, one has the following: Every square-integrable function U ∈ L 2 (I 0 ) can be expanded in the VFPs as follows: where the coefficients u i are obtained as follows: A finite form of the series in (10) is considered an approximation to U (t), that is, the following: where U and Ψ(t) are the following vectors: Due to the need of the orthogonality of derivatives of the shifted VFPs in the subject of the approximation error, we show that the mth derivative of these polynomials is orthogonal with respect to the weight function w m (t) = t m+ 1 2 (1 − t) m+ 1 2 , m = 1, 2, . . .. First, we obtain Rodrigues' formula for extracting the VFPs and their derivatives.
The mth derivative of V F n (t) can be stated as follows: By differentiating (8) and comparing the leading terms in , the value of A m n is determined as follows: Then, from (13) one has the following: Using the integration by parts for the second integral in (14) results in the following: By performing the integration by parts (k − m − 1) times for I, one has the following: If k < n, using the integration by parts again, one obtains the following result: are orthogonal with respect to the weight function w m (t). If k = n, one has from (13) and (15) the following:

Operational Matrix of the Integration of the Integer Order
Integrating the components of vector Ψ(t) leads to the following: Now, t i+1 is approximated in terms of shifted VFPs as follows: Using the last equality in (18), the integral in (17) can be written as follows: Therefore, one obtains the following: where Π is the integral operational matrix of the integer order and its entries are as follows: , k, j = 1, 2, . . . , N + 1.

Operational Matrix of the Integration of the Fractional Order
The Riemann-Liouville integral of the kth component of vector Ψ(t) is computed as follows: The term t i+µ can be approximated as follows: where the coefficients ρ (i) j , j = 1, 2, ..., N + 1 can be computed as follows: .
Thus, relation (20) is written as the following: If the last relation is written as a matrix form, one has the following: where Π (µ) is the operational matrix of the fractional order µ > 0 and its entries are calculated as follows: , k, j = 1, 2, . . . , N + 1.

Approximating the Integral Part with the Singular Kernel
A matrix relation is derived to approximate the integral parts in Equation (1) with the weakly singular kernels as shown below: Since t k−σ+1 = t −σ t k+1 , the term t k+1 is approximated in terms of the Vieta-Fibonacci basis. So, relation (22) is written as follows: By rewriting the last approximation as a matrix form, one has the following: where S (σ) is the operational matrix for approximating the integral part with the weakly singular kernels and its entries are as follows: , k, j = 1, 2, . . . , N + 1.

Delay Operational Matrix
To obtain a matrix relation between Ψ(qt) and Ψ(t), the kth component of the vector Ψ(qt) is approximated in terms of the shifted VFPs. According to representation (8), one has the following: Now, t i is approximated as the following: Therefore, equality (24) is as follows: The matrix form of the last relation is as follows: where Q (q) is the delay operational matrix. Using approximation (25), one can approximate the following integral: If the last approximation is represented as a matrix form, one has the following: where R (r) is a (N + 1) × (N + 1) matrix with the following entries: , k, j = 1, 2, . . . , N + 1.

Operational Matrix of the Product
After substituting appropriate approximations into integral parts, the expression Ψ(t) Ψ T (t) H may appear in which H is a vector. Here, an approximation is given for this multiplication expression as follows: where H is the product operational matrix. The matrix form of (27) is as follows: So, the ith row of (28) is as follows: Taking inner product and using (9) result in the following: Now, the product ψ k (t) ψ m (t) is written as follows: where coefficients τ (k,m) l are computed from Algorithms 1 and 2. By (30), the integral on the left-hand side of (29) is computed as follows: Thus, from the right-hand side of (29) and the last equality, the entries of H are as follows: , i, m = 1, 2, . . . , N + 1.

Methodology
To better understand the implementation of the proposed method, consider Equation (1) again. To present an approximation for the function Z, first an approximation is given for Applying the Riemann-Liouville integral operator to approximation (31) and using the operational matrix in (21) and the initial condition in (2), an approximation is achieved for Z (t): where U T = C T Π (η) + F T . Using approximation (32) and the operational matrix in (25), an approximation is obtained for Z (qt) as follows: where V T = C T Π (η) Q (q) + F T Q (q) . The functions θ i , i = 1, 2, 3 are approximated as follows: where Θ i , i = 1, 2, 3 are the known (N + 1) × (N + 1) matrices for which their entries are computed by means of the inner product. Utilizing approximations (32)-(34), the operational matrices in (23) and (26), and the definition of the product operational matrix in (27), the integral parts are approximated as follows: Substituting approximations (31)-(37) into Equation (1) leads to the following residual function: Collocating Equation (38) at roots of the V F N+2 (t) (V F N+2 (t) is of degree N + 1) results in a system involving (N + 1) algebraic equations. By solving the resultant system and using (32), an approximate solution is obtained.

Error Bound
In this section, an error bound for the residual function is computed in a weighted Sobolev space. Let Ω = {X | X C ≤ 1} and w(t) = √ t − t 2 be the weight function. The Vieta-Fibonacci-weighted Sobolev space is introduced as follows: equipped with the following norm and semi-norm: where . w k stands for the known L 2 -norm.
Proof. Because U N+1 (t) is the Vieta-Fibonacci approximation to U (t), one has the following: Then, for m ∈ N ∪ {0}, one has the following: From (16), one obtains the following: Similarly, one has the following: Using the Stirling formula [34], one deduces the following: From (41)-(43), one obtains the following: Then, one has the following: Using Definition 3 implies the desired result for λ = [λ] + δ, 0 < δ < 1 as follows: To find an error bound for U (t) − U N+1 (t), the same procedure of proving in Theorem 3 is pursued. So, for m ≤ s one has the following: Therefore, the following bound is obtained for 0 < λ < s: where C † is a positive constant.
Proof. Using the second property in Definition 2, one has the following: where the star denotes the convolution. Using Young's convolution inequality f * g q ≤ C 0 f 1 g q and Corollary 1, one has the following: Now, suppose that Z (t) is the exact solution of Equation (1), and Z N+1 (t) is the approximate solution obtained from the proposed scheme. Therefore, Z N+1 (t) satisfies the following equation: where R(t) is the residual function or perturbation term. By Subtracting Equation (46) from Equation (1), the following error equation is obtained: Assume that the functions g and p are square-integrable on I 0 , the functions θ i , i = 1, 2, 3 are square-integrable on J = I 0 × I 0 , s ∈ N ∪ {0}, and 0 < λ < s. Using the obtained bounds in (40) and (45), an error bound is achieved for R(t).
The right-hand side of inequality (48) shows that by appropriately choosing the number of basis functions, the approximation error will decrease.

Illustrated Examples
In this section, the Vieta-Fibonacci collocation method is utilized to solve four FDIDEs as special cases of Equation (1). Maximum absolute errors (MAEs) of the obtained results are computed by the following formula: Example 1 ([14]). As the first example, consider the following equation: with the exact solution Z (t) = t 2 and the initial condition Z (0) = 0. For examining the effect of N on the numerical solution, this example is solved for different values of N. According to the results in Figure 1, the MAEs are decreased by increasing N. A graphical comparison of the exact and approximate solutions and the absolute error function are seen in Figure 2 for N = 25. Values of absolute errors are listed in Table 1 at the equally spaced points t j = 0.2j, j = 0, 1, . . . , 5 and compared with ones reported in [14], which have used the fractional-order Euler polynomials to find an approximate solution. The results show a good agreement of the approximate solution with the exact one.    [11,14]: with the exact solution Z (t) = exp(t) − 1 as η = 1 and the initial condition Z (0) = 0. The above equation is solved by the presented method, and the MAE is 5.3398 × 10 −6 for η = 1 and N = 6. This equation was solved in [11,14], and the reported MAEs are 1.9487 × 10 −6 for N = 6 and 4.7215 × 10 −6 for N = 20, respectively. From these values of errors, it can be seen that the suggested algorithm has given results with the same accuracy with a smaller number of basis functions. The figures of the exact and approximate solutions are plotted in Figure 3 for η = 1, N = 6. As seen, the approximate solution has good agreement with the exact one. In Figure 4, approximate solutions are observed for N = 6 and η = 0.7, 0.8, 0.9, 1. It can be seen that the numerical solutions get close to the exact one as η → 1.  Example 3. The equation in Example 2 is considered again with the exact solution Z (t) = t 2 if η = 1. In here, the function f is as below [14]: The equation is solved by the suggested scheme, and the MAE is 3.0170 × 10 −6 for η = 1, N = 6. This equation was solved in [14] using a modified algorithm based on constructing fractional-order Euler polynomials. The reported MAE is 4.5311 × 10 −7 for m = 6, α = 1 2 . As seen, the obtained results are close together. The figures of the exact and approximate solutions are plotted in Figure 5 for η = 1, N = 6. As seen, the approximate solution has good agreement with the exact one. In Figure 6, approximate solutions are observed for N = 6 and η = 0.7, 0.8, 0.9, 1. It can be seen that the numerical solutions get close to the exact one as η → 1. The values of absolute errors of the approximate solution at equally spaced points t j = 0.2j, j = 0, 1, . . . , 5 are seen in Table 2 for N = 6 and η = 0.8, 0.9, 1. It can be found that when η → 1, the values of the errors decrease.   Example 4. Consider the following FDIDE with the singular kernels [14]: The exact solution is Z (t) = t 2 if η = 1 and the initial condition is Z (0) = 0. This equation is solved by the presented method and the value of the MAE is 4.7670 × 10 −20 for η = 1 and N = 3, while the MAE reported by [14] is 5.511 × 10 −16 for ν = α = 1 and m = 3. So, the proposed method has given the more accurate result with the same number of basis functions. The figures of the exact and approximate solutions are plotted in Figure 7 for η = 1, N = 5. As seen, the approximate solution has good agreement with the exact one. In Figure 8, approximate solutions are observed for N = 5 and η = 0.6, 0.7, 0.8, 0.9, 1. It can be seen that the numerical solutions converge the exact one as η → 1.

Conclusions
In this paper, the VFPs were applied to set a new matrix spectral method for solving FDIDEs. First, operational matrices of the integration of fractional and integer orders, a pseudo-operational matrix for the integral parts with singular kernels, a delay operational matrix, and a product operational matrix were constructed, and then the collocation method along with the obtained operational matrices were used for reducing the FDIDEs to an algebraic system. The existence and uniqueness of the solution of the given problem were investigated and proved, using Krasnoselskii's fixed point theorem. An error bound of the approximation error was computed in a Vieta-Fibonacci-weighted Sobolev space, which showed that by increasing the number of basis functions, the approximation error is decreased. A Rodrigues' formula was presented to extract the VFPs and their derivatives,