differential equations.

Over the past century, polynomial sequences have been a central focus of mathematical research, with Lucas polynomials standing out due to their profound theoretical properties and broad applicability. These polynomials are instrumental across multiple mathematical disciplines, including algebra, where they elucidate structural patterns, and geometry, where they aid in shape analysis. Furthermore, their combinatorial characteristics provide valuable insights into enumeration problems, while their utility in approximation theory offers robust mechanisms for function estimation and numerical problem-solving. In statistical modeling and number theory, Lucas polynomials reveal intricate connections to prime distribution and integer properties. As a specific subclass related to Chebyshev polynomials, they have been rigorously analyzed, highlighting their versatility and enduring significance in both theoretical and applied mathematics [2, 3].

Concurrently, fractional calculus has emerged as a powerful generalization of classical differentiation and integration, providing a sophisticated framework for modeling phenomena with non-integer order dynamics. This mathematical paradigm is particularly effective in capturing complex behaviors in science and engineering. Fractional differential equations (FDEs) are now fundamental in diverse fields. For instance, in acoustics, they offer refined models for wave propagation in heterogeneous media. In material science, fractional derivatives accurately describe the viscoelastic behavior and damping characteristics of advanced materials [4].

Moreover, fractional calculus plays a crucial role in electroanalytical chemistry, where complex electrochemical processes are intricately described by fractional differential equations, proving indispensable for characterizing electrode kinetics and diffusion-limited reactions [5]. Piecewise polynomial collocation methods have also been demonstrated to be highly effective for linear Volterra integro-differential equations with weakly singular kernels [6]. In neuroscience, fractional-order models prove advantageous in capturing the long-range dependencies and memory effects characteristic of neuronal systems [7]. Additionally, fractional calculus provides refined models for diffusion processes, accurately depicting the non-local nature of particle movement [8].

Spectral methods have played a pivotal role in solving differential equations, offering efficient and accurate numerical solutions. Researchers have extensively utilized spectral algorithms to tackle a wide range of problems, showcasing their versatility in handling complex physical phenomena, nonlinear partial differential equations, and fractional sub-diffusion equations [9, 10, 11, 12, 13, 14, 15].

Recent advancements in spectral methods and numerical analysis have significantly expanded the toolkit for solving complex fractional, stochastic, and multi-dimensional differential equations. Notably, researchers have developed enhanced collocation and Petrov–Galerkin frameworks utilizing various orthogonal polynomials—such as shifted Jacobi, hybrid fractional Gegenbauer, and second- and third-kind Chebyshev polynomials—to tackle time-fractional models ranging from the FitzHugh–Nagumo equation to Euler–Bernoulli beam dynamics and nonlinear integro-differential equations [16, 17, 18, 19]. Furthermore, the application of spectral and collocation techniques has been successfully extended to n-dimensional hyperbolic equations with trigonometric sources and heat conduction problems in cylindrical domains, demonstrating rigorous convergence and high computational efficiency [20, 21]. These contributions collectively underscore the versatility and robustness of spectral algorithms in addressing diverse challenges in modern applied mathematics and computational physics.

Numerical integration is a cornerstone of computational mathematics. Among the various techniques, Gauss quadrature stands out for its exceptional accuracy and efficiency. By strategically selecting the integration nodes as the roots of orthogonal polynomials (such as Legendre polynomials) and optimizing the corresponding weights, Gauss quadrature achieves exact integration for polynomials of high degrees using minimal points. This optimality makes it highly superior to uniform grid methods, providing rapid convergence and reliable error estimation.

Motivated by these advancements, this paper introduces a novel algorithm for solving fractional-order integro-differential equations. The proposed method leverages the Gauss quadrature formula and the spectral collocation method, incorporating Lucas polynomials to address complex mathematical challenges. Section 2 presents mathematical preliminaries, including key definitions in fractional calculus and relevant properties of Lucas polynomials. Section 3 outlines the algorithmic implementation. Section 4 provides the convergence analysis, while Section 5 offers numerical examples to demonstrate efficiency and compares the results with existing methods, including the Fibonacci-Fejér [1] approach. The paper concludes in Section 6 with a summary of

Some definitions and properties of fractional calculus. We present fundamental notations and definitions of fractional calculus theory utilized throughout this article.

Definition 1. The Riemann-Liouville fractional integral operator Iα of order α on the Lebesgue space L1[0, 1] {︄ ∫︁ is defined as: 1 t α−1 α Γ(α) 0 (t − τ) f(τ) dτ, α > 0, I f(t) = (1) f(t), α = 0.

The operator Iα possesses the following properties:

(i) IαIβ = Iα+β, (ii) (IαIβf)(t) = (IβIαf)(t), α ν Γ(ν + 1) ν+α (iii) I (t − a) = (t − a) , Γ(ν + α + 1) where f ∈ L1[0, 1], α, β ⩾ 0, and ν > −1. Γ(·) denotes the gamma function.

Definition 2. Let f(t) ∈ Cn[0, ℓ] with ℓ > 0. The Caputo fractional differential operator D∗ α is defined by:

α 1 n−α−1 (n) (D∗ f)(t) = (t − τ) f (τ) dτ, n − 1 < α ≤ n, n ∈ N. (2) Γ(n − α) 0

The operator D∗ α satisfies the following properties for n − 1 < α ≤ n:

(D∗ αIαf)(t) = f(t), n ∑︂ −1 f (k)(0+) (IαD∗ αf)(t) = f(t) − tk, t > 0,

α β Γ(β + 1) β−α D∗ t = t . (3) Γ(β + 1 − α) For further details on fractional calculus, see [22, 23, 24].

Some properties of Lucas polynomials. Lucas polynomials, denoted as Ln(t), satisfy the recurrence relation [25, 26]:

Lr(t) = t Lr−1(t) + Lr−2(t), r = 2, 3, . . . , with initial conditions L0(t) = 2 and L1(t) = t.

These polynomials possess the following analytic form:

⌊∑︂ n/2⌋ (︃ )︃ n n − k n−2k Ln(t) = t . (4) n − k k

Theorem 1. The Caputo fractional derivative of Ln(t) is given by:

⌊(n− ∑︂ α)/2⌋ (︃ )︃ α n n − k Γ(n − 2k + 1) n−2k−α D∗ Ln(t) = t . (5) n − k k Γ(n − 2k − α + 1)

∑︁⌊n/2⌋ n (︁n−k)︁ n−2k Proof : By the analytic form of Lucas polynomials in Eq. (4), we have Ln(t) = k=0 n−k k t . Applying the Caputo fractional derivative operator D∗ α to this sum and using the linearity of the operator along with property (3), specifically D∗ αtβ = Γ( Γ( β+1 β+1) −α) tβ−α, we obtain:

⌊∑︂ n/2⌋ (︃ )︃ α n n − k Γ(n − 2k + 1) n−2k−α D∗ Ln(t) = t .

n − k k Γ(n − 2k − α + 1)

The upper limit of the sum is effectively ⌊(n − α)/2⌋ because for n − 2k − α < 0, the terms either vanish or are handled by the definition of the fractional derivative of constants (which is zero for the Caputo derivative). Thus, the result follows. □

tions

In this section, we employ the collocation method aided by the Gauss quadrature formula to solve fractional integro-differential equations (FIDEs) of the form:

(︃ ∫︂ t )︃ F (t, u(t), D∗ σ1 u(t), . . . , D∗ σℓ u(t)) = G t, u(t), H (τ, u(τ)) dτ , (6)

subject to the initial conditions:

u(i)(0) = ξi, i = 0, 1, . . . , ℓ − 1, (7) where r − 1 < σr ⩽ r for r = 1, 2, . . . , ℓ.

First, we approximate u(t) using a Lucas polynomial expansion:

∑︂ N u(t) ≈ ak Lk(t) = U(t). (8)

Substituting Eq. (8) into Eq. (6) yields:

(︃ ∫︂ t )︃ F (t, U(t), D∗ σ1 U(t), . . . , D∗ σℓ U(t)) ≈ G t, U(t), H (τ, U(τ)) dτ . (9)

We apply the Gauss quadrature transformation x ↦→ 2 t (1 + ξ) to the integral, resulting in:

∫︂ t ∑︂ n (︃ )︃ ∼ t t f(x) dx = wk f (1 + ξk) , (10) 0 2 2

where ξk are the roots of the Legendre polynomial Pn+1(x), and wk are the corresponding Gauss weights. Applying this quadrature to the right-hand side of Eq. (9) gives:

(︄ n (︃ )︃)︄ t ∑︂ t (︁ t )︁ F (t, U(t), D∗ σ1 U(t), . . . , D∗ σℓ U(t)) ≈ G t, U(t), wk H (1 + ξk), U (1 + ξk) . (11) 2 2 2 1 (︁ (︁ jπ )︁)︁ Finally, we collocate Eq. (11) at the shifted Gauss-Lobatto nodes zj = 2 1 + cos N for j = 0, 1, . . . , N − ℓ, to obtain:

F (zj, U(zj), D∗ σ1 U(zj), . . . , D∗ σℓ U(zj)) (︄ n )︄ ∑︂ (︂ (︁ )︁)︂ (12) zj zj zj ≈ G zj, U(zj), wk H (1 + ξk), U (1 + ξk) , j = 0, 1, . . . , N − ℓ.

2 2 2

Incorporating the initial conditions:

U (j)(0) = ξj, j = 0, 1, . . . , ℓ − 1, (13)

Eqs. (12)–(13) generate a system of N + 1 equations in the unknown expansion coefficients {ak : k = 0, . . . , N}. This system is solved using Newton’s iterative method to obtain the approximate solution U(t).

We establish the convergence of the Lucas polynomial approximate solutions through the following theorems, providing a rigorous foundation for the reliability of our proposed algorithm.

Theorem 2. If ϕ(z) is defined on [0, ℓ] and its derivatives satisfy |ϕ(i)(r)| ≤ Qi for i ≥ 0, where r ∈ (0, ℓ) and Q is a positive constant, and if ϕ(z) has the expansion ϕ(z) = k=0 ck Lk(z), then the coefficients satisfy:

σ Qk+1 |ck| < ,

where σ is a constant dependent on Q and ℓ.

Proof : The Lucas polynomials Lk(z) can be related to Chebyshev polynomials, and they form a complete basis for functions analytic on [0, ℓ]. Since ϕ(z) has bounded derivatives |ϕ(i)(r)| ≤ Qi, it is analytic in a neighborhood of [0, ℓ]. The expansion coefficients ck can be represented via contour integration using Cauchy’s integral formula. By estimating the contour integral and utilizing the growth properties of Lucas polynomials, we obtain the bound |ck| < σQ k k !+1 , where σ is a constant dependent on the domain of analyticity and the maximum modulus of the basis functions. □

Theorem 3. If f(x) obeys the assumptions of Theorem 2, and we define the global error as ∑︁ EN = maxx∈(0,ℓ) |eN (x)|, where eN (x) = k=N+1 ck Lk(x), then the following estimation holds:

˜ ζN EN < C ,

where C ˜ and ζ are constants dependent on Q and ℓ, demonstrating spectral convergence.

Proof : The global error is defined as eN (x) = k=N+1 ckLk(x). Taking the absolute value and applying the triangle inequality, we have |eN (x)| ≤ k=N+1 |ck||Lk(x)|. Using the bound from Theorem 2 and the maximum value of L (x) on [0, ℓ], which is bounded by Cρk √ℓ2 + 4 k where ρ = , we can bound the tail of the series. Substituting the bound for ck, we get:

∑︂ ∞ σQk+1 ∑︂ ∞ (Qρ)k |eN (x)| ≤ Cρk = σCQ .

k! k!

k=N+1 k=N+1

C ˜ ζN C ˜ Recognizing the tail of the exponential series, this sum is bounded by N! , where ζ = Qρ and is a constant. This establishes the spectral convergence of the proposed method. □

In this section, we validate the applicability and accuracy of our proposed Lucas-Gauss Collocation Method (LGCM) through several numerical test problems, explicitly comparing the outcomes with existing methods, including the Fibonacci-Fejér [1] approach.

Example 1. Consider the following linear FIDE [27]:

D∗ 2 u(t) = (cos t − sin t) u(t) + t sin τ u(τ) dτ + f(t), t ∈ [0, 1], (14)

where 8t3/2 2√t 3 2 f(t) = + 2t + + t cos t − t sin t − 3t cos t, 3√π √π subject to u(0) = 0. The exact solution is u(t) = t + t2. Applying the LGCM for n = 6 and N = 2, we obtain the coefficients a0 = −1, a1 = 1, and a2 = 1. Consequently, U(t) = −L0(t) + L1(t) + L2(t) = −2 + t + (t2 + 2) = t + t2, which perfectly matches the exact solution. As illustrated in Figure 1a, the approximate solution completely overlaps with the exact one, and the absolute error in Figure 1b confirms the machine-precision accuracy of the proposed method.

Example 2. Consider the following nonlinear FIDE [27]:

D∗ u(t) = f(t)u(t) + g(t) + √t (u(τ)) dτ, t ∈ [0, 1], (15)

where f(t) = 2√t + 2t 2 − √t(1 + t) ln(1 + t), −1 (︁ )︁ 2 sinh √t 3/2 g(t) = − 2t , √π√t + 1 Example 4. Consider the following nonlinear FIDE:

D∗ 2 u(t) + u(t) = 2 (t − 2τ)u(τ) dτ, t ∈ [0, 1], (17)

subject to u(0) = 0. Since the exact solution is unavailable, we define the absolute residual error norm as: ⃓ ∫︂ t ⃓ ⃓ 1 ⃓ RE = max ⃓D∗ 2 U(t) + U(t) − 2 (t − 2τ)U(τ) dτ⃓ . (18) t∈(0,1) ⃓ ⃓

Applying our method at N = 18, n = 7 and α = 0.5 to get Table 5, which illustrates the RE at different values of t. Since the exact solution is unavailable, the approximate solution U(t) obtained by the LGCM is shown in Figure 4a. The absolute residual error (RE) at different time points is visualized in Figure 4b, confirming the extreme precision of the method.

This paper introduced a novel spectral algorithm designed to compute numerical solutions for fractional integro-differential equations. The methodology seamlessly integrates the Gauss quadrature rule with the spectral collocation method, utilizing Lucas polynomials as the basis functions. A key advantage of this approach is its broad applicability to both linear and nonlinear FIDEs, offering a robust framework for diverse mathematical challenges.

The developed algorithm demonstrates exceptional accuracy while requiring only a modest number of retained modes in the Lucas expansion. This efficiency significantly reduces computational overhead, making it highly practical for complex systems. Comparative analyses with existing techniques, including the Fibonacci-Fejér [1] collocation method, CAS wavelets, and SCW, consistently highlight the superior precision and faster convergence of the proposed Lucas-Gauss scheme.

Future research may extend this framework to higher-dimensional problems or incorporate adaptive mesh strategies to further enhance computational efficiency. All computational codes were developed and executed in Mathematica 11 on an HP Z420 Workstation (Intel Xeon E5-1620, 3.6 GHz, 16 GB RAM).

The authors would like to express their sincere gratitude to the editors and anonymous reviewers for their invaluable comments and constructive feedback, which significantly contributed to the enhancement of this paper.

The authors declare that no external funding or support was received for the research presented in this paper, including administrative, technical, or in-kind contributions.