EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 1, 2024, 477-503 ISSN 1307-5543 – ejpam.com Published by New York Business Global A spectral collocation method for solving Caputo-Liouville fractional order Fredholm integro-differential equations Khaled M. Saad1,2, M. Q. Khirallah1,3,∗ 1 Department of Mathematics, College of Sciences and Arts, Najran University, Najran, State, Saudi Arabia 2 Department of Mathematics, Faculty of Applied Science, Taiz University, Taiz, Yemen 3 Department of Mathematics and Computer Science, Faculty of Science, Ibb University, Ibb, Yemen Abstract. In this paper, a numerical method for solving the fractional order Fredholm integro- differential equations via the Caputo-Liouville derivative is presented. The method uses the well- known shifted Chebyshev expansion and a truncated series to represent the unknown function. It also incorporates numerical integration techniques like the Trapezoidal, Simpson’s 1/3, and Simpson’s 8/3 methods. The paper also provides an approximation for the derivative of an integer. The procedure converts the provided problem into a system of algebraic equations using shifted Chebyshev coefficients and collocation points. The coefficients are found by solving this system using well-known techniques like Newton’s method. Numerical results are presented graphycally to illustrate the applicability, efficacy, and accuracy of the approach presented in this work. All calculations in this study were performed using the MATHEMATICA software program. 2020 Mathematics Subject Classifications: 74Sxx, 97Nxx Key Words and Phrases: Fractional order integro-differential equations, Caputo type fractional derivative, the shifted Chebyshev spectral collocation method, Trapezoidal, Simpson 1. Introduction A subfield of mathematics known as fractional calculus extends the idea of derivatives and integrals to non-integer orders. Fractional calculus uses fractional or real numbers for the order of differentiation or integration rather than whole numbers. Fractional deriva- tives and fractional integrals are two important ideas in fractional calculus. The rate at which a function changes in relation to a variable of order α is represented by Dα, the fractional derivative of the function. Similar to this, a generalization of integration is ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i1.5049 Email addresses: khaledmasd@hotmail.com (K.M.Saad), mqm73@yahoo.com (M.Q.Khirallah) https://www.ejpam.com 477 © 2024 EJPAM All rights reserved. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 478 represented by the fractional integral of a function with regard to a variable of order D−α ( [18], [17], [13]). Numerous areas of mathematical physics and engineering applications deal with frac- tional integral-differential equations. A great deal of attention has been focused on devel- oping efficient techniques for getting approximate or numerical solutions for both linear and nonlinear fractional integro-differential equations because of the difficulties in obtain- ing analytical solutions for these problems. Furthermore, using numerical or approximat- ing methods to solve fractional integro-differential equations containing realistic nonlinear elements is still a challenging undertaking. Integrals and derivatives of an unknown function are combined in the integro-differential equation . Different kinds of functional equations, such as integral and integro-differential equations, stochastic equations, and ordinary or partial differential equations, arise when real-world issues are mathematically modeled. In many different domains, including physics, astronomy, potential theory, fluid dynamics, biological models, and chemical kinetics, fractional integral-differential equations are used to mathematically formulate physical processes. Fractional integro-differential equations are sometimes difficult to solve analytically, requiring the construction of effective approximation solutions. The Jacobi spectral method [20], Runge Kutta method [24], Chebyshev collocation method [3], Laplace Power Series Method [1], rationalized Haar functions method [15], Galerkin methods with hybrid functions[14] and Laguerre collocation method [5] are just a few of the numerical techniques that have been used to solve such equations. Numerous applications can be also found for the well-known set of orthogonal poly- nomials defined on the interval [−1, 1], known as Chebyshev polynomials [11, 19]. Their advantageous qualities in function approximation are the reason for their extensive use. When it comes to Chebyshev polynomials, the wide range of qualities that orthogonal poly- nomials have is especially concise, which makes them stand out above other orthogonal polynomials. These polynomials are members of the unique class of orthogonal polynomi- als called Jacobi polynomials. Chebyshev polynomials offer advantages in terms of orthog- onality, error minimization, and convergence properties within specific intervals. However, their limited applicability outside these intervals and challenges in certain mathematical operations may be considered disadvantages in certain contexts. Jacobi polynomials are solutions to Sturm-Liouville equations and correspond to weight functions of the kind (1− β)α(1 + β)α [16]. For instance, the orthogonality condition of the Chebyshev polynomials is utilized to approximate the functions of the period [a, b]. In these techniques, which strongly rely on polynomials, (see ( [21])). There are several advantages to employee shifted Chebyshev polynomials: Chebyshev polynomial exhibit a multitude of intriguing and beneficial properties. Utilizing Chebyshev polynomials as fundamental functions yields highly precise solutions. The utilization of Chebyshev polynomials in research contributions is comparatively limited in comparison to other polynomial types. By selecting the modified set of shifted Chebyshev polynomials as the basis functions and retaining only a few terms of the modes, it becomes feasible to generate highly accurate approximations with reduced computational effort. Furthermore, Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 479 the associated errors are minimal. The structure of this study is as follows. The definitions of the fractional derivatives and shifting Chebyshev polynomials are briefly discussed in Section 2 as well as some preliminary remarks. We demonstrate the numerical application of the suggested method and applications in Sections 3 and 4. Section 5 provides the conclusion. 2. Preliminaries and notations 2.1. Some definitions of fractional derivatives Definition 1. The fractional derivative of order 0 < α ≤ 1 in the Caputo sense is provided for ϕ(β) ∈ H1(0, b) by: CDαϕ(β) = 1 Γ(1− α) ∫ β 0 ϕ ′ (τ) (β − τ)ν dτ, β > 0, Definition 2. where H1(0, b) is the Sobolev space and is given by H1(0, b) = { ϕ ∈ L2(0, b) : dϕ dβ ∈ L2(0, b), L2(0, b) = {ϕ(β) : ( ∫ b 0 ϕ(β)2dβ) 1 2 < ∞}, } Dαβm = { 0, m ∈ {0, 1, 2, . . . , ⌈α⌉ − 1}, Γ(m+1) Γ(m+1−α)β m−α, m ∈ N ∧m ≥ ⌈α⌉, where ⌈α⌉ the ceiling function of α and N = 1, 2, 3, · · ·. 2.2. The shifting Chebyshev polynomials and function approximations In this section , we give the definitions of the shifted Chebyshev polynomials (CPs), their notations, and their properties . The majority of our studies have concentrated on an orthogonal polynomial class. The recurrence relations and analytical equations of these polynomials can be used to construct a family of orthogonal polynomials called Chebyshev polynomials . Now, we will provide a quick review of the definitions and formulas related to the first-type Chebyshev polynomials in this section. It is well-known that the first-kind Chebyshev polynomials are defined on the interval [−1, 1] as follows (see, for details, [16, 23]; see also the recently-published survey-cum- expository review article [9] on the Chebyshev and related orthogonal polynomials): The range [−1, 1] is where the first-type Chebyshev polynomials are typically defined, as follows (see [16, 23] for more details; additionally, see the recently published survey and expository review in [9] on Chebyshev and similar orthogonal polynomials). Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 480 Ψn(γ) = cos(nθ) (n ∈ N0 := N ∪ {0} = 0, 1, 2, · · · ), (1) where γ = cos(θ). The Chebyshev polynomials {Ψn(γ)}n∈N0 can be obtained from the following recur- rence relation: Ψn+1(γ) = 2γΨn(γ)−Ψn−1(γ) (n ∈ N) ( Ψ0(γ) = 1; Ψ1(γ) = γ ) . (2) The Chebyshev polynomials {Ψn(γ)}n∈N0 are orthogonal over the interval [−1, 1] with the weight function (1− γ2)− 1 2 and we have the following orthogonality property:∫ 1 −1 (1− γ2)− 1 2 Ψi(γ)Ψj(γ) dξ =  0 (i ̸= j) π 2 (i = j ̸= 0) π (i = j = 0). (3) The following is the exact formula for the Chebyshev polynomial: Ψn(γ) = n 2 [n/2]∑ i=0 (−1)i (n− i− 1)! (i)! (n− 2i)! (2γ)n−2i. (4) We define the shifted Chebyshev polynomials on the interval [0, 1] by setting the vari- able γ = 2β − 1. The following expressions describe these polynomials: Φs(β) = Φs(2β − 1) = β2s( √ β), where a set of orthogonal Chebyshev polynomials over the range [0, 1] is generated by the polynomial collection {Φ2s(β)}s∈N0 . Calculating the specific expression of the shifted Chebyshev polynomial is a straight- forward task. T̄s(ζ) of degree s as follows (see [16]): Φs(β) = s s∑ k=0 (−1)s−k 22k (s+ k − 1)! (2k)! (s− k)! βk, (5) where Φ0(β) = 1 and Φ1(β) = 2β − 1. Using a linear combination of the first (m + 1) terms of Φs, we expand and evaluate the function Ω(β) spanning the interval [0, 1]. We find that: Ω(β) ≃ Ωm(β) = m∑ i=0 aiΦi(β). (6) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 481 The coefficients ai are determined by: ai =  1 π ∫ 1 0 Ω(η) Φi(β)√ β − β2 dβ (i = 0) 2 π ∫ 1 0 Ω(β)Φi(β)√ β−β2 dβ (i ∈ N). (7) The primary approximate expression for the derivative of ϕm(β) is provided in the theorem that follows. Theorem 1. [10, 22] In Eq. (6), the approximate solution of the main problem is given in terms of shifted Chebyshev polynomials polynomials. Following that, the fractional-order terms can be changed into the following algebraic equations: Dα ( Ωm(β) ) = m∑ i=⌈α⌉ i−⌈α⌉∑ k=0 ciχ (α) i,k β i−k−α, (8) χ (α) i,k = (−1)k 4i−k2iΓ(2i− k)Γ(i− k + 1) Γ(k + 1)Γ(2i− 2k + 1)Γ(i− k + 1− α) , (9) where Γ(.) is the gamma function. 2.3. Error Analysis This section focuses specifically on introducing the convergence analysis and assessing the upper limit of the error associated with the proposed formula. Theorem 2. [7] Suppose that the function Ω(β) is so constrained that Ω′′(β) ∈ L2[0, b] and |Ω′′ (β)| ≦ c, where c is a constant. Then the series (6) of the shifted Chebyshev expansion is uniformly convergent and: |aℓ| < c ℓ2 , (ℓ ∈ 1, 2, ...). (10) Theorem 3. [7] Suppose that Ω(β) ∈ Cm[0, 1]. Then the error in approximating the function Ω(β) by Ωm(β) by using the formula (6) can be bounded by: ∥ϕ(β)− ϕm(β)∥ ≦ ℘∆m+1 (m+ 1)! √ π 2 and ℘ = maxt∈[0,1]ϕ (m+1)(β) (11) (∆ = max{β0, β − β0}). Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 482 3. Approach to Fractional Fredholm Integro-Differential Equation Solving In this section, we present the schema for the following nonlinear fractional Fredholm integro-differential equation: Dαϕ(β) = G ( β, ϕ(β), ∫ 1 0 H(β, ϕ(β))dβ ) , 0 < β ≤ 1 , n− 1 < α ≤ n. (12) Here, we use the shifted Chebyshev polynomials collocation method and Theorem 1 to solve (12) as follows m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = G β, m∑ j=0 cj Φj(β), ∫ 1 0 H β, m∑ j=0 cj Φj(β)  dβ  . (13) Therefore, we use the following numerical methods for integration to analyze the system of equations given in equation (13): (i) Trapezoidal’s Method m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = G β, m∑ j=0 cj Φj(β), h 2 ( F (β0) + F (βL) + 2 L−1∑ k=1 F (βk) ) . (14) At these points, βs, s = 0, 1, ...,m− α, we collocate (14). m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = G βs, m∑ j=0 cj Φj(βs), h 2 ( F (β0) + F (βL) + 2 L−1∑ k=1 F (βk) ) , (15) where F (β) = H ( β, ∑m j=0 cj Φj(β) ) . (ii) Simpson’s 1/3 Method m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 ciχ (α) i,k β j−k−α = G β, m∑ j=0 cj Φj(β), h 2 ( F (β0) + F (βL) + 2 L 2 −1∑ k=1 F (β2k) + 4 L 2∑ k=1 F (β2k−1) ) . (16) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 483 At these points, βs, s = 0, 1, ...,m− α, we collocate (16). m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = G βs, m∑ j=0 cj Φj(β), h 2 ( F (β0) + F (βL) + 2 L 2 −1∑ k=1 F (β2k) + 4 L 2∑ k=1 F (β2k−1) ) , (17) where F (β) = H ( β, ∑m j=0 cj Φj(β) ) . (iii) Simpson’s 3/8 Method m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = G ( β, m∑ j=0 cj Φj(β), 3h 8 ( F (β0) + F (βL) + 3 L 3∑ k=1 (F (β3k−2) + F (β3k−1)) +2 n 3 −1∑ k=1 F (β3k) )) . (18) At these points, βs, s = 0, 1, ...,m− α, we collocate (18). m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = G ( βs, m∑ j=0 cj Φj(β), 3h 8 ( F (β0) + F (βL) + 3 L 3∑ k=1 (F (β3k−2) + F (β3k−1)) +2 L 3 −1∑ k=1 F (β3k) )) , (19) where F (β) = H ( β, ∑m j=0 cj Φj(β) ) . The roots of the shifted Chebyshev Polynomials are used to find appropriate collocation locations Φm+1−⌈α⌉. Additionally, we can get r equations by inserting (6) in the boundary conditions. Equations (15) or (17) or (19), when combined with the r equations of the boundary conditions, gives (m + 1) of an algebraic equation system that can be solved using the Newton iteration method for the unknowns cj , j = 0, 1, ..., .,m. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 484 4. Numerical Examples In this section we present three examples of fractional Fredholm integro-differential using the proposed approach. Example 1. Consider the following fractional Fredholm integro-differential equa- tion [8] Dαϕ(β) = βeβ + eβ − β + ∫ 1 0 β ϕ(β) dβ, 0 < α ≤ 1, (20) subject to the initial condition ϕ(0) = 0, (21) with exact solution ϕ(β) = βeβ. We apply the provided procedure and arrive at an approximation of the solution as, ϕm(β) = m∑ j=0 cjΦj(β). (22) Using the equations provided by the Trapezoidal method (15), Simpson’s method (17), and Simpson’s 3/8 method (19), we construct the schema as follows: (i) Trapezoidal’s Method Using (14) and (15), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) k βj−k−α = φ(β) + h 2 ( β0 m∑ j=0 cjΦj(β0) + 2 L−1∑ k=1 βk m∑ j=0 cjΦj(βk) + βL m∑ j=0 cjΦj(βL) ) , (23) where φ(β) = βeβ + eβ − β, and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + h 2 ( β0 m∑ j=0 cjΦj(β0) + 2 L−1∑ k=1 βk m∑ j=0 cjΦj(βk) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 485 + βL m∑ j=0 cjΦj(βL) ) . (24) (ii) Simpson’s 1/3 Method Using (16) and (17), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ(β) + h 3 ( β0 m∑ j=0 cjΦj(β0) + 2 L 2 −1∑ k=1 β2k m∑ j=0 cjΦj(β2k) + 4 L 2∑ k=1 β2k−1 m∑ j=0 cjΦj(β2k−1) + βL m∑ j=0 cjΦj(βL) ) , (25) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + h 3 ( β0 m∑ i=0 cjΦj(β0) + 2 L 2 −1∑ k=1 β2k m∑ j=0 cjΦj(β2k) + 4 L 2∑ k=1 β2k−1 m∑ j=0 cjΦj(β2k−1) + βL m∑ j=0 cjΦj(βL) ) . (26) (iii) Simpson’s 3/8 Method Using (18) and (19), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ(β) + 3h 8 ( β0 m∑ j=0 cjΦj(β0) + 3 L 3∑ k=1 ( β3k−2 m∑ j=0 ciΦj(β3k−2) + β3k−1 m∑ j=0 cjΦj(β3k−1) ) + 2 L 3 −1∑ k=1 β3k m∑ j=0 cjΦj(β3k) + βL m∑ j=0 cjΦj(βL) ) , (27) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + h 3 ( β0 m∑ j=0 cjΦj(β0) + 2 L 2 −1∑ k=1 β2k m∑ j=0 cjΦj(β2k) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 486 + 4 L 2∑ k=1 β2k−1 m∑ j=0 cjΦj(β2k−1) + βL m∑ j=0 cjΦj(βL) ) , (28) where βs are the roots of the shifted Chebyshev polynomial and s = 0, 1, 2, 3, ...,m. The initial condition (21) can be written as ϕm(0) = m∑ j=0 cjΦj(0) = m∑ j=0 (−1)jcj = 0. (29) To acquire the coefficients ’cj ’ in the preceding three cases, one can solve algebraic equations (24), (26) and (28), corresponding to equation (29) in each case. Ultimately, by replacing the coefficients ’cj ’ in equation (22), one can obtain an approximate numerical solution for equation (20). Now, we present various figures to illustrate the numerical results. A comparison of the exact and approximate solutions with α = 0.8, 0.9, 1 and m = 6 is shown in Figure 1(a). This comparison applies specifically to Trapezoidal’s case, while the remaining two cases exhibit identical behavior. In this graphical representation, we observe the trends of the approximate solutions for various α values. These solutions exhibit regular behavior, and their proximity increases as α approaches toward the integer value. The different value of α is highlighted in the Figure 1(a). Figure 1(b) illustrates the corresponding absolute error for the Trapezoidal, Simpson ’ 3/8 and Simpson ’ 1/3 methods. To further verify, considering the absence of an exact solution in the non-integer case, it becomes crucial to assess the error. Therefore, to confirm the validity of our approach, we compute the absolute error in a two-step process, i.e. |ϕm+1(β)− ϕm(β)|. The error in a two-step process in Figure 1(c) is plotted for the same values as in Figures 1(a) and 1(b). It is clear from these figures that the order of the error is very small. Example 2. Consider the following fractional Fredholm integro-differential equation Dαϕ(β) = 2− 7β2 3 + 2β + ∫ 1 0 β2 ϕ(β) dβ, (30) subject to the initial condition ϕ(0) = 1. (31) Using the suggested approach , we deduce the following approximation for the solution: ϕm(β) = m∑ j=0 cjΦj(β). (32) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 487 Using the Trapezoidal method (15), Simpson’s method (17), and Simpson’s 3/8 method (19), we construct the schema as follows: (i) Trapezoidal’s Method Using (14) and (15), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ(β) + h 2 ( β2 0 m∑ j=0 ciΦj(β0) + 2 L−1∑ k=1 β2 k m∑ j=0 cjΦj(βk) + β2 L m∑ i=0 ciΦj(βL) ) , (33) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + h 2 ( β2 0 m∑ j=0 cjΦj(β0) + 2 L−1∑ k=1 β2 k m∑ j=0 cjΦj(βk) + β2 L m∑ j=0 cjΦj(βL) ) . (34) (ii) Simpson’s 1/3 Method Using (16) and (17), we obtain m∑ j=⌈α⌉ m−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ(β) + h 3 ( β2 0 m∑ j=0 cjΦj(β0) + 2 L 2 −1∑ k=1 β2 2k m∑ j=0 cjΦj(β2k) + 4 L 2∑ k=1 β2 2k−1 m∑ j=0 cjΦj(β2k−1) + β2 L m∑ j=0 cjΦj(j, βL) ) , (35) and n∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + h 3 ( β2 0 m∑ j=0 cjΦj(β0) + 2 L 2 −1∑ k=1 β2 2k m∑ j=0 cjΦj(β2k) + 4 L 2∑ k=1 β2 2k−1 m∑ j=0 cjΦj(β2k−1) + β2 L m∑ j=0 cjΦj(βL) ) . (36) (iii) Simpson’s 3/8 Method Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 488 Using (18) and (19), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ(β) + 3h 8 ( β2 0 m∑ j=0 cjΦj(j, β0) + 3 L 3∑ k=1 ( β2 3k−2 m∑ j=0 cjΦj(j, β3k−2) + β3k−1 m∑ j=0 cjΦj(j, β3k−1) ) + 2 L 3 −1∑ k=1 β2 3k m∑ j=0 cjΦj(β3k) + β2 L m∑ j=0 cjΦj(βL) ) , (37) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ(β) + 3h 8 ( β2 0 m∑ j=0 cjΦj(β0) + 3 L 3∑ k=1 ( β2 3k−2 m∑ j=0 cjΦj(β3k−2) + β3k−1 m∑ j=0 cjΦj(β3k−1) ) + 2 L 3 −1∑ k=1 β2 3k m∑ j=0 cjΦj(β3k) + β2 L m∑ j=0 cjΦj(βL) ) . (38) We follow the same procedure as described in example one. Approximate solutions for a variety of α values are shown in Figure 3 (a). This comparison is specifically relevant to the Simpson’s 1/8 case, while the other two cases demonstrate similar behavior. The absolute error between the approximate solutions via Trapezoidal, Simpson ’ 1/8, and Simpson ’ 1/3 methods and the exact solution is shown in Figure 3 (b). In Figure 4, the absolute error between each subsequent step when the non-integer α values via Trapezoidal, Simpson ’ 1/8, and Simpson ’ 1/3 methods are applied is shown. Collectively, these numerical outcomes illustrate the precision of the approximations. It has been demonstrated that augmenting the number of steps by m enhances accuracy. Example 3. Consider the following fractional Fredholm integro-differential equa- tion [4] Dαϕ1(β) = 2 + 12 5 β − ∫ 1 0 β(ϕ2 1(β) + ϕ2(β) 2)), dβ, (39) Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 489 Dαϕ2(β) = −2 + 4 3 β − ∫ 1 0 β(ϕ2 1(β)− ϕ2(β) 2)), dβ, (40) 0 < α ≤ 2, subject to the initial conditions ϕ1(0) = 1, ϕ ′ 1(0) = 0, ϕ2(0) = 1, ϕ ′ 2(0) = 0, (41) with exact solutions ϕ1 = 1 + β2 and ϕ1 = 1− β2 (i) Trapezoidal’s Method Using (14), we get m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ1(β) + h 2 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 + 2 L−1∑ k=1 βk ( m∑ j=0 cjΦj(βk) )2 + 2 L−1∑ k=1 βk ( m∑ j=0 djΦj(βk) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (42) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α = φ1(β) + h 2 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(β0) )2 + 2 L−1∑ k=1 βk ( m∑ j=0 cjΦj(βk) )2 − 2 L−1∑ k=1 βk ( m∑ j=0 djΦj(βk) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ j=0 djΦj(βL) )2) , (43) where φ1(β) = 2 + 12 5 β, φ2(β) = −2 + 4 3 β. Using (15), we get m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ1(β) + h 2 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 + Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 490 2 L−1∑ k=1 βk ( m∑ j=0 cjΦj(βk) )2 + 2 L−1∑ k=1 βk ( m∑ j=0 djΦj(βk) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (44) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α s = φ1(β) + h 2 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(β0) )2 + 2 L−1∑ k=1 βk ( m∑ j=0 cjΦj(βk) )2 − 2 L−1∑ k=1 βk ( m∑ j=0 djΦj(βk) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ j=0 djΦj(βL) )2) . (45) (ii) Simpson’s 1/3 Method Using (16), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ1(β) + h 3 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 cjΦj(β2k) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 4 L 2∑ k=1 β2k−1 ( m∑ j=0 cjΦj(β2k−1) )2 + 4 L 2∑ k=1 β2k ( m∑ i=0 diΦj(j, β2k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (46) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α = φ2(β) + h 3 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(β0) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 cjΦj(β2k) )2 − 2 L 2 −1∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 491 4 L 2∑ k=1 β2k−1 ( m∑ j=0 cjΦj(β2k−1) )2 − 4 L 2∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ j=0 djΦj(βL) )2) . (47) Now, Using (17), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ1(βs) + h 3 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 cjΦj(β2k) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 4 L 2∑ k=1 β2k−1 ( m∑ j=0 cjΦj(β2k−1) )2 + 4 L 2∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (48) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α s = φ2(βs) + h 3 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(jβ0) )2 + 2 L 2 −1∑ k=1 β2k ( m∑ j=0 cjΦj(β2k) )2 − 2 L 2 −1∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 4 L 2∑ k=1 β2k−1 ( m∑ j=0 cjΦj(β2k−1) )2 − 4 L 2∑ k=1 β2k ( m∑ j=0 djΦj(β2k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ j=0 djΦj(βL) )2) . (49) (iii) Simpson’s 3/8 Method Using (18), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α = φ1(β) + 3h 8 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 492 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 cjΦj(β2k−2) )2 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 djΦj(β2k−2) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 cjΦj(β3k−1) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 djΦj(β3k−1) )2 2 L 3 −1∑ k=1 β3k ( m∑ j=0 cjΦj(β3k) )2 + 2 L 2∑ k=1 β3k ( m∑ j=0 djΦj(β3k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (50) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α = φ2(β) + 3h 8 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(β0) )2 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 cjΦj(β2k−2) )2 − 3 L 3∑ k=1 β2k−2 ( m∑ j=0 djΦj(β2k−2) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 cjΦj(β3k−1) )2 − 3 L 3∑ k=1 β3k−1 ( m∑ j=0 djΦj(β3k−1) )2 2 L 3 −1∑ k=1 β3k ( m∑ j=0 cjΦj(β3k) )2 − 2 L 2∑ k=1 β3k ( m∑ j=0 djΦj(β3k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ j=0 djΦj(βL) )2) . (51) Using (19), we obtain m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 cjχ (α) j,k β j−k−α s = φ1(βs) + 3h 8 ( β0 ( m∑ j=0 cjΦj(β0) )2 + β0 ( m∑ j=0 djΦj(β0) )2 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 cjΦj(β2k−2) )2 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 djΦj(β2k−2) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 cjΦj(β3k−1) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 djΦj(β3k−1) )2 Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 493 2 L 3 −1∑ k=1 β3k ( m∑ j=0 cjΦj(β3k) )2 + 2 L 2∑ k=1 β3k ( m∑ j=0 djΦj(β3k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 + βL ( m∑ j=0 djΦj(βL) )2) , (52) and m∑ j=⌈α⌉ j−⌈α⌉∑ k=0 djχ (α) j,k β j−k−α s = φ2(βs) + 3h 8 ( β0 ( m∑ j=0 cjΦj(β0) )2 − β0 ( m∑ j=0 djΦj(β0) )2 + 3 L 3∑ k=1 β2k−2 ( m∑ j=0 cjΦj(β2k−2) )2 − 3 L 3∑ k=1 β2k−2 ( m∑ j=0 djΦj(β2k−2) )2 + 3 L 3∑ k=1 β3k−1 ( m∑ j=0 cjΦj(β3k−1) )2 − 3 L 3∑ k=1 β3k−1 ( m∑ j=0 djΦj(β3k−1) )2 2 L 3 −1∑ k=1 β3k ( m∑ j=0 cjΦj(β3k) )2 − 2 L 2∑ k=1 β3k ( m∑ j=0 djΦj(β3k) )2 + βL ( m∑ j=0 cjΦj(βL) )2 − βL ( m∑ i=0 diΦj(βL) )2) . (53) Additionally, we follow the same stages outlined in examples 1 and 2. Figures 5 (a) and 5 (a) show the approximate solutions and how they coincide with the exact solution for α = 0.8, 0.9 and α = 1. The absolute error between the approximate solutions and the exact solution is shown in Figures 5(b) and 6(b). When the order of the derivative becomes close to the integer number, the approximate solutions approach to the exact solution. In this example, the comparison applies specifically to Simpson 3/8’s case, while the remaining two cases exhibit identical behavior. Furthermore, for the non-integer order in two sequential approximations, Figure 7 con- firms the accuracy of the approximate solutions with α = 0.8 and m = 6 and m = 7, via Trapezoidal, Simpson ’ 1/8, and Simpson ’ 1/3 methods. In the preceding three examples, we observe a consistent pattern in the behavior of the approximate solutions. They tend to converge towards the exact solution as the order of the non-integer derivative approaches the integer order. This observation contributes positively to the overall presentation of this work. Additionally, when dealing with an non-integer order, the absolute error between successive approximate solutions for various values of m yielded accurate results. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 494 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Β Φ HΒL HaL 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 Β  ΦH Β L-Φ 7 HΒL   HbL Trapezoidal Simpson ' 3 8 Simpson ' 1 3 Figure 1: (a) Combining approximate solutions with the exact solution for different values of alpha for example 1. (Red solid color: α = 0.8; Blue solid color: α = 0.9; Black solid color: α = 1; Green dashed color: Exact solution). (b) The absolute error between the approximate solutions and the analytical solution for example 1. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 495 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 Β  Φ 6 HΒL - Φ 7 HΒL   Figure 2: Plotting the difference between the two step of the approximate solutions with α = 0.8 and m = 6 and m = 7 for example 1. (Black dashed color: Trapezoidal Green dashed color: Simpson ’ 3/8 Simpson ’ 1/3: Red solid color). Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 496 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Β Φ HΒL HaL 0.0 0.2 0.4 0.6 0.8 1.0 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 Β  ΦH Β L-Φ 7 HΒL   HbL Trapezoidal Simpson ' 3 8 Simpson ' 1 3 Figure 3: (a) Combining approximate solutions with the exact solution for different values of alpha for example 2. (Red solid color: α = 0.8; Blue solid color: α = 0.9; Black solid color: α = 1; Green dashed color: Exact solution). (b) The absolute error between the approximate solutions and the analytical solution for example 2. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 497 0.0 0.2 0.4 0.6 0.8 1.0 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Β  Φ 6 HΒL - Φ 7 HΒL   Figure 4: Plotting the difference between the two step of the approximate solutions with α = 0.8 and m = 6 and m = 7 for example 2. (Black dashed color: Trapezoidal Green dashed color: Simpson ’ 3/8 Simpson ’ 1/3: Red solid color). Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 498 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.2 1.4 1.6 1.8 2.0 Β Φ 1 HΒL HaL 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 0.00002 0.00004 0.00006 0.00008 0.0001 Β  Φ 1 HΒL - Φ 1 ,7 HΒL   HbL Trapezoidal Simpson ' 1 3 Simpson ' 3 8 Figure 5: (a) Combining approximate solutions with the exact solution for different values of alpha for example 3. (Red solid color: α = 0.8; Blue solid color: α = 0.9; Black solid color: α = 1; Green dashed color: Exact solution). (b) The absolute error between the approximate solutions and the analytical solution for example 3. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 499 0.0 0.2 0.4 0.6 0.8 1.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Β Φ 2 HΒL HbL 0.0 0.2 0.4 0.6 0.8 1.0 0.00000 0.00005 0.00010 0.00015 0.00020 Β  Φ 2 HΒL - Φ 2 ,7 HΒL   HbL Trapezoidal Simpson ' 1 3 Simpson ' 3 8 Figure 6: (a) Combining approximate solutions with the exact solution for different values of alpha for example 3. (Red solid color: α = 0.8; Blue solid color: α = 0.9; Black solid color: α = 1; Green dashed color: Exact solution). (b) The absolute error between the approximate solutions and the analytical solution for example 3. Khaled M. Saad, M. Q. Khirallah / Eur. J. Pure Appl. Math, 17 (1) (2024), 477-503 500 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 Β  Φ 1 ,6 HΒL - Φ 1 ,7 HΒL   HaL 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 0.0005 0.0010 0.0015 Β  Φ 2 ,6 HΒL - Φ 2 ,7 HΒL   HbL Figure 7: (a) Plotting the difference between the two step of the approximate solutions u with α = 0.8 and m = 6 and m = 7 for example 3. (Black dashed color: Trapezoidal Green dashed color: Simpson ’ 3/8 Simpson ’ 1/3: Red solid color). (b) Plotting the difference between the two step of the approximate solutions v with α = 0.8 and m = 6 and m = 7 for example 3. (Black dashed color: Trapezoidal Green dashed color: Simpson ’ 3/8 Simpson ’ 1/3: Red solid color). REFERENCES 501 5. Conclusions This study used the Caputo fractional derivative in conjunction with the Chebyshev spectral approach to solve fractional integro-differential equations. The Trapezoidal, Simp- son’s 1/3, and Simpson’s 8/3 methods combined with the properties of Chebyshev poly- nomials to convert fractional integro-differential equations into algebraic equations. The resulting equations were then solved using well-known techniques like Newton’s. The numerical results was carried out using the MATHMETICA soft program. We suggest emphasizing the incorporation of fractional space-time derivatives in our forthcoming re- search. Additionally, we plan to transform the fractional time derivative into a discrete equation using unconventional finite-difference techniques. To streamline intricate models into a set of solvable differential equations, we may also utilize special additional polyno- mial functions [2, 6, 12]. References [1] A. K. Alomari, M. Alaroud, T. N. Tahat, and A. Almalki. Extended laplace power series method for solving nonlinear caputo fractional volterra integro-differential equa- tions. Symmetry, 15:1296–1296, 2023. [2] M. Asif, I. Khan, N. Haider, and Q. Al-Mdallal. Legendre multi-wavelets collocation method for numerical solution of linear and nonlinear integral equations. Alexandria Engineering Journal, 59:5099–5109, 2020. [3] A. G. Atta and Y. H. Youssri. Advanced shifted first-kind chebyshev collocation approach for solving the nonlinear time-fractional partial integro-differential equation with a weakly singular kernel. Computational Applied Mathematics, 41, 2022. [4] H. O. Bakodah, M. Al-Mazmumy, and S. O. Almuhalbedi. Solving system of inte- gro differential equations using discrete adomian decomposition method. Journal of Taibah University for Science, 13:805–812, 2019. [5] D. V. Bayram and A. Dascioglu. A method for fractional volterra integro-differential equations by laguerre polynomials. Advances in Difference Equations, 2018, 2018. [6] Ali. F. Jameel, N. R. Anakira, A. K. Alomari, and Noraziah H. Man. Solution and analysis of the fuzzy volterra integral equations via homotopy analysis method. Computer Modeling in Engineering Sciences, 127:875–899, 2021. [7] W. G and M. A. Snyder. Chebyshev methods in numerical approximation. Mathe- matics of Computation, 22:894–894, 1968. [8] B. D. Garba and S. L. Bichi. On solving linear fredholm integro-differential equations via finite difference-simpson’s approach. Malaya Journal of Matematik, 8:469–472, 2020. REFERENCES 502 [9] Hari Mohan H. M. Srivastava. A survey of some recent developments on higher tran- scendental functions of analytic number theory and applied mathematics. Symmetry, 13:2294, 2021. [10] M. M. Khader. On the numerical solutions for the fractional diffusion equation. Com- munications in Nonlinear Science and Numerical Simulation, 16:2535–2542, 2011. [11] M. M. Khader and K. M. Saad. On the numerical evaluation for studying the frac- tional kdv, kdv-burgers and burgers equations. The European Physical Journal Plus, 133, 2018. [12] I. Khan, M. Asif, R. Amin, Q. Al-Mdallal, and F. Jarad. On a new method for finding numerical solutions to integro-differential equations based on legendre multi-wavelets collocation. Alexandria Engineering Journal, 61:3037–3049, 2022. [13] A. A. Kilbas, H. M Srivastava, and J. Trujillo. Theory and applications of fractional differential equations. North-holland Mathematics Studies, 204:vii–x, 2006. [14] K. Maleknejad and M.Tavassoli Kajani. Solving linear integro-differential equation system by galerkin methods with hybrid functions. Applied Mathematics and Com- putation, 159:603–612, 2004. [15] K. Maleknejad, F. Mirzaee, and S. Abbasbandy. Solving linear integro-differential equations system by using rationalized haar functions method. Applied Mathematics and Computation, 155:317–328, 2004. [16] J. C. Mason and D. C. Handscomb. Chebyshev Polynomials. Chapman Hall/CRC, 2002. [17] K. S. Miller and B. Ross. An introduction to the fractional calculus and fractional differential equations. Wiley, 1993. [18] I. Podlubny. Fractional differential equations : an introduction to fractional deriva- tives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, 1998. [19] K. M. Saad and H. M. Srivastava. Numerical solutions of the multi-space fractional- order coupled korteweg–de vries equation with several different kernels. Fractal and fractional, 7:716–716, 09 2023. [20] Harendra Singh and Ramta Ram Pathak. Jacobi spectral method for the fractional reaction–diffusion equation arising in ecology. Mathematical Methods in the Applied Sciences, 2024. [21] H. M. Srivastava, K. M. Saad, and W. M. Hamanah. Certain new models of the multi-space fractal-fractional kuramoto-sivashinsky and korteweg-de vries equations. Mathematics, 10:1089, 2022. REFERENCES 503 [22] N. H. Sweilam and M. M. Khader. A chebyshev pseudo-spectral method for solving fractional order integro-differential equations. The ANZIAM Journal, 51:464–475, 2010. [23] G. Szego. Orthogonal polynomials. American Mathematical Society, 2003. [24] G. Zhang and R. Zhu. Runge–kutta convolution quadrature methods with conver- gence and stability analysis for nonlinear singular fractional integro–differential equa- tions. Communications in Nonlinear Science and Numerical Simulation, 84:105132– 105132, 2020.