EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6970 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Scheme for Fractional Volterra Integro-Differential Equations: Handling Modified Atangana-Baleanu Operators Kamran1,∗, Syed Musanif Shah1, Fady Hasan2,∗, Nabil Mlaiki2,∗ 1 Department of Mathematics, Islamia College Peshawar, Peshawar 25120, Khyber Pakhtunkhwa, Pakistan 2 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, Riyadh 11586, Saudi Arabia Abstract. This study presents a comparative numerical investigation of Volterra integro-differential equations (VIDEs) incorporating the Modified Atangana-Baleanu fractional derivative in the Ca- puto sense. Fractional-order VIDEs are vital for modeling phenomena in biology, physics, and engineering. However, their non-local nature and complexity often preclude analytical solutions, necessitating efficient numerical approaches. We develop a numerical technique based on the Laplace transform (LT). The governing equation is first transformed into an algebraic equation in the Laplace domain. This transformed equation is solved algebraically, and the solution is then numerically inverted back to the time domain. Two highly effective numerical methods are employed for inverting the Laplace transform: the Gauss-Hermite quadrature and the improved Talbot’s method. The efficacy of these inversion methods is validated through numerical experi- ments, and the results are compared with exact solutions. Our findings confirm the convergence of both methods. The numerical results demonstrate the adaptability and accuracy of the suggested approaches, establishing them as promising tools for solving fractional VIDEs. 2020 Mathematics Subject Classifications: 44A10, 65R20 Key Words and Phrases: Modified Atangana-Baleanu-Caputo derivative, Volterra integro- differential equations, Gauss-Hermite quadrature, Talbot’s method, existence theory 1. Introduction Fractional calculus (FC), the branch of mathematics focused on fractional-order opera- tors, has garnered considerable interest recently due to its immense potential for modeling complicated processes across a variety of scientific fields. This field extends standard ∗Corresponding author. ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6970 Email addresses: kamran.maths@icp.edu.pk (Kamran), fhasan@psu.edu.sa (F. Hasan), nmlaiki@psu.edu.sa, nmlaiki2012@gmail.com (N. Mlaiki) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 2 of 22 calculus by including non-integer-order derivatives and integrals, which allows for more realistic representations of systems with memory and hereditary properties. For instance, FC enhances market prediction models in economics by capturing long-term dependencies and memory effects [1]. In electromagnetic theory, it offers a deeper understanding of signal processing and wave propagation [2]. Within control theory, it enables the design of efficient and adaptive control systems [3], and for viscoelastic materials, it accurately describes frequency-dependent damping and material behavior [4]. Furthermore, FC sup- ports advanced control strategies in robotics and addresses complex engineering challenges [5–8]. Equations involving fractional-order operators allow for the accurate analysis of non- local and memory-dependent processes, overcoming the limitations of classical calculus. Researchers have developed a variety of analytical and numerical methods for solving such equations. Notable techniques include the sinc-collocation method [9], the Legendre collocation method [10], Laguerre polynomials [11], the Adomian decomposition method [12, 13], the Variational Iteration Method [14], and the natural transform method [15]. These approaches offer powerful tools for solving complex fractional-order systems. Fractional-order Volterra integro-differential equations (FOVIDEs) have attracted con- siderable interest due to their capacity to model memory-dependent processes, such as population growth and decay, viscoelasticity, and heat conduction. Unlike integer-order derivatives, fractional derivatives incorporate historical information and introduce memory effects, which enhances modeling accuracy but increases the complexity of finding analyt- ical solutions. Hence, developing efficient numerical methods is necessary for solving such problems. Numerous numerical methods have been devised to obtain solutions for FOVIDEs. For example, the authors of [16] used a Bernoulli polynomials approximation technique to solve nonlinear FOVIDEs. In [17], a numerical technique based on Legendre wavelets was proposed, while [18] utilized a rationalized Haar wavelets method for a system of fractional-order Fredholm-Volterra IDEs. The reproducing kernel method was applied to FOVIDEs in [19] and for nonlinear cases in [20]. A hybrid numerical scheme was developed for FOVIDEs in [21]. Furthermore, the authors of [22] developed the least squares and shifted Legendre methods, and the Bernoulli wavelets method was utilized in [23]. Finally, a Chebyshev technique was developed to solve a system of FOVIDEs in [24]. For decades, researchers have crafted different versions of fractional-order derivatives to capture the behavior of systems with memory and hereditary properties more accurately. The Riemann-Liouville (RL) and Caputo (LC) definitions of fractional-order derivatives emerged as early frontrunners and remain among the most common today. Each derivative has proven extremely helpful in practical applications, but each also comes with its own disadvantages. For example, the RL derivative requires initial conditions to be expressed in terms of fractional integrals, which can be challenging to interpret physically and cre- ates hurdles for real-world application. The LC derivative overcomes this limitation by using integer-order initial conditions, which are more intuitive. However, the kernel of the LC derivative contains a singularity of the form (t − τ)−α, which can cause serious computational difficulties as t approaches τ . Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 3 of 22 To overcome these issues, Caputo and Fabrizio introduced a new derivative with a non- singular exponential kernel, called the Caputo-Fabrizio (CF) derivative [25]. However, this new derivative had its own limitation: the rapid decay of its exponential kernel restricts the ability to capture long-term memory effects. To circumvent this issue, Atangana and Baleanu introduced a more generalized derivative based on the Mittag-Leffler (ML) kernel [26]. Unlike the exponential kernel, the ML kernel decays slowly, enabling more precise modeling of memory effects. This derivative is known as the Atangana-Baleanu (AB) derivative and is defined in both the RL and Caputo senses. Recently, the authors of [27] proposed an extension of the ABC derivative, called the modified ABC (MABC) derivative. Their work demonstrated that the MABC derivative can solve a broader class of differential equations that were previously intractable under the ABC formulation. Building upon this advancement, our work focuses on solving FOVIDEs that incorporate the MABC derivative. For this purpose, we employ the Laplace transform method coupled with numerical inversion techniques. Specifically, we use two established approaches for the numerical inversion of the Laplace transform: the Gauss- Hermite quadrature method and the improved Talbot’s method, ensuring both accuracy and computational efficiency. 1.1. Preliminaries Definition 1. The MABC derivative of fractional order β, where 0 < β < 1, is defined as follows [28]:( MABC 0 Dβ τ y ) τ = N (β) 1− β [ y(τ)− Eβ ( − µβτ β ) y(0)− µβ ∫ τ 0 (τ − t)β−1Eβ,β ( − µβ(τ − t)β ) y(t)dt ] . Definition 2. The Laplace transform of the MABC derivative of fractional order β, where 0 < β < 1, is defined as follows [28]: L { MABC 0 Dβ τ y(τ) } = N (β) 1− β zβL{y(τ)} − y(0)zβ−1 zβ + ν , ∣∣∣∣ νzβ ∣∣∣∣ < 1. Definition 3. The one-parameter Mittag-Leffler (ML) function, as given in [28], is defined as: Eβ(τ) = ∞∑ k=0 τk Γ(kβ + 1) , where τ ∈ C and ρ is an arbitrary positive constant. Definition 4. The Laplace transform of the one-parameter Mittag-Leffler function is de- fined as follows [28]: L{Eβ(−ντβ)} = zβ−1 zβ + ν , where Re(z) > |ν|1/β. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 4 of 22 Definition 5. The two-parameter Mittag-Leffler function is defined as follows [28]: Eβ,ρ(τ) = ∞∑ k=0 τk Γ(kβ + ρ) , where τ ∈ C, and β, ρ are arbitrary positive constants. Definition 6. The Laplace transform of the two-parameter Mittag-Leffler function is de- fined as [28]: L{τρ−1Eβ,ρ(−ντβ)} = zβ−ρ zβ + ν , where Re(z) > |ν|1/β. 2. Existence and Uniqueness of Solution The study of FOVIDEs lies at the intersection of integral equation theory and FC. For modeling systems with memory and hereditary properties, such as those observed in biophysics, population dynamics, and viscoelasticity, problems of the form (1) are impor- tant. However, the combination of non-local fractional differential and integral operators introduces significant complexity. Establishing the existence of a solution confirms the well-posedness of these models, while proving the uniqueness confirms that the problem predicts a single solution, which is a fundamental requirement for any physical application. Let χ = [0, 1]. Consider the Volterra integro-differential equation:( MABC0Dβ τ y ) (τ) = f(τ) + ∫ τ 0 K(τ, t)y(t)dt, y(0) = y0, (1) where 0 < β < 1, K ∈ C(χ × χ,R), f ∈ C(χ,R), and the unknown function y ∈ C(χ,R) with the supremum norm ∥y∥ = supτ∈χ |y(τ)|. Lemma 1. The solution of (1) satisfies the integral equation: y(τ) = y0 + 1− β Q(β) f(τ) + 1− β Q(β) ∫ τ 0 K(τ, t)y(t)dt + β Q(β)Γ(β) [∫ τ 0 (τ − t)β−1f(t)dt+ ∫ τ 0 ∫ s 0 (s− t)β−1K(s, t)y(t)dtds ] . (2) Proof. Apply the MABC fractional integral operator MABC0Iβ τ to both sides of (1). Using the inversion property MABC0Iβ τ ( MABC0Dβ τ y ) (τ) = y(τ)−y(0) and linearity yields (2). Define the operator T : C(χ) → C(χ) by: (Ty)(τ) = y0 + 1− β Q(β) f(τ) + 1− β Q(β) ∫ τ 0 K(τ, t)y(t)dt+ β Q(β)Γ(β) ∫ τ 0 (τ − t)β−1f(t)dt + β Q(β)Γ(β) ∫ τ 0 ∫ s 0 (s− t)β−1K(s, t)y(t)dtds. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 5 of 22 A fixed point of T is a solution to (2) and hence to (1). Define the constants: ς1 = ∣∣∣∣1− β Q(β) ∣∣∣∣ , ς2 = sup (τ,t)∈χ×χ |K(τ, t)|, ς3 = ∣∣∣∣ β Q(β)Γ(β) ∣∣∣∣ , ς4 = sup τ∈χ |f(τ)|. Theorem 1. The problem defined in (1) has at least one solution provided that ς1ς2 + ς3ς2 Γ(β + 2) < 1. Proof. The proof proceeds in four steps using Schaefer’s fixed point theorem. Let X = C([0, 1],R) be the Banach space endowed with the supremum norm ∥y∥ = supτ∈[0,1] |y(τ)|. Step 1: Continuity of the Operator T The operator T : X → X is continuous. Indeed, let {yn} be a sequence in X such that yn → y uniformly on [0, 1]. Then,∣∣(Tyn)(τ)− (Ty)(τ) ∣∣ ≤ ς1 ∫ τ 0 |K(τ, t)||yn(t)− y(t)|dt + ς3 ∫ τ 0 ∫ s 0 (s− t)β−1|K(s, t)||yn(t)− y(t)|dtds ≤ ς1ς2∥yn − y∥ ∫ τ 0 dt+ ς3ς2∥yn − y∥ ∫ τ 0 ∫ s 0 (s− t)β−1dtds = ς1ς2∥yn − y∥τ + ς3ς2∥yn − y∥ τβ+1 Γ(β + 2) . Taking the supremum over τ ∈ [0, 1] yields: ∥Tyn − Ty∥ ≤ ( ς1ς2 + ς3ς2 Γ(β + 2) ) ∥yn − y∥. Since yn → y, it follows that ∥Tyn − Ty∥ → 0 as n → ∞. Therefore, T is continuous. Step 2: Boundedness of T We show that T maps bounded sets into bounded sets. Let Bq = {y ∈ X : ∥y∥ ≤ q}. For any y ∈ Bq and τ ∈ [0, 1], we have: |(Ty)(τ)| ≤ |y0|+ ς1|f(τ)|+ ς1 ∫ τ 0 |K(τ, t)||y(t)|dt+ ς3 ∫ τ 0 (τ − t)β−1|f(t)|dt + ς3 ∫ τ 0 ∫ s 0 (s− t)β−1|K(s, t)||y(t)|dtds ≤ |y0|+ ς1ς4 + ς1ς2qτ + ς3ς4 τβ β + ς3ς2q τβ+1 Γ(β + 2) . Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 6 of 22 Since τ ∈ [0, 1], the right-hand side is bounded by a constant M(q) independent of τ and y ∈ Bq. Hence, T (Bq) is bounded. Step 3: Equicontinuity and Compactness of T We show that T maps bounded sets into equicontinuous sets. Let y ∈ Bq and τ1, τ2 ∈ [0, 1] with τ1 < τ2. Then, |(Ty)(τ2)− (Ty)(τ1)| ≤ ς1|f(τ2)− f(τ1)|+ ς1 ∣∣∣∣∫ τ2 0 K(τ2, t)y(t)dt− ∫ τ1 0 K(τ1, t)y(t)dt ∣∣∣∣ + ς3 ∣∣∣∣∫ τ2 0 (τ2 − t)β−1f(t)dt− ∫ τ1 0 (τ1 − t)β−1f(t)dt ∣∣∣∣ + ς3 ∣∣∣∣∫ τ2 0 ∫ s 0 (s− t)β−1K(s, t)y(t)dtds− ∫ τ1 0 ∫ s 0 (s− t)β−1K(s, t)y(t)dtds ∣∣∣∣ . The first term tends to zero as |τ2 − τ1| → 0 by the uniform continuity of f on [0, 1]. The remaining terms involve integrals over intervals whose lengths tend to zero, and the integrands are bounded. Standard estimates show that each term can be made arbitrarily small independently of y ∈ Bq. Therefore, T (Bq) is equicontinuous. By the Arzelà-Ascoli theorem, T is compact. Step 4: A Priori Boundedness and Schaefer’s Theorem Let S = {y ∈ X : y = λTy for some λ ∈ (0, 1)}. For any y ∈ S, we have y = λTy, so |y(τ)| ≤ |(Ty)(τ)|. Using the estimate from Step 2 with q = ∥y∥, we get: ∥y∥ ≤ |y0|+ ς1ς4 + ς1ς2∥y∥+ ς3ς4 1 β + ς3ς2∥y∥ 1 Γ(β + 2) . Rearranging terms yields: ∥y∥ ( 1− ς1ς2 − ς3ς2 Γ(β + 2) ) ≤ |y0|+ ς1ς4 + ς3ς4 β . By the given condition, the coefficient on the left is positive. Hence, ∥y∥ is bounded by a constant independent of λ, so S is bounded. Since T is continuous and compact, and S is bounded, Schaefer’s fixed point theorem implies that T has a fixed point, which is a solution to (1). Theorem 2. The solution to (1) is unique. Proof. Let y1 and y2 be two solutions of (2). Let z(τ) = y1(τ)− y2(τ). Then, |z(τ)| = |(Ty1)(τ)− (Ty2)(τ)| ≤ ς1 ∫ τ 0 |K(τ, t)||z(t)|dt+ ς3 ∫ τ 0 ∫ s 0 (s− t)β−1|K(s, t)||z(t)|dtds ≤ ς1ς2 ∫ τ 0 |z(t)|dt+ ς3ς2 ∫ τ 0 ∫ s 0 (s− t)β−1|z(t)|dtds. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 7 of 22 To apply Grönwall’s inequality, we note that∫ τ 0 ∫ s 0 (s− t)β−1|z(t)|dtds = ∫ τ 0 |z(t)| (∫ τ t (s− t)β−1ds ) dt = 1 β ∫ τ 0 (τ − t)β|z(t)|dt. Thus, |z(τ)| ≤ ς1ς2 ∫ τ 0 |z(t)|dt+ ς3ς2 β ∫ τ 0 (τ − t)β|z(t)|dt. Since (τ − t)β ≤ 1 for τ, t ∈ [0, 1], we have |z(τ)| ≤ ( ς1ς2 + ς3ς2 β )∫ τ 0 |z(t)|dt. By the integral form of Grönwall’s inequality (with constant kernel), we conclude that |z(τ)| = 0 for all τ ∈ [0, 1]. Hence, y1 = y2, proving uniqueness. 3. Methodology In this section, we outline the proposed numerical strategy for solving FOVIDEs in- volving the MABC derivative. The approach unfolds in the following sequence: • We begin by applying the Laplace transform to turn the given FOVIDE into a more manageable algebraic equation in the Laplace domain, • This transformed equation is then solved within that domain, • To retrieve the solution in its original form, we apply the inverse Laplace transform. However, since analytical inversion is often intractable, we employ numerical techniques. Specifically, we utilize the Gauss-Hermite quadrature and the improved Talbot’s method to efficiently approximate the inverse Laplace transform and obtain the final solution. Applying the Laplace transform to (1), we obtain L { MABC 0 Dβ τ y(τ) } = L { f(τ) + ∫ τ 0 K(τ ,t)y(τ)dt } , Using the Laplace transform of the MABC derivative, we obtain ⇒ N (β) 1− β zβL{y(t)} − y(0)zβ−1 zβ + ν = L{f(τ) } + L {∫ τ 0 K(τ ,t)y(τ)dt } , which implies N (β) 1− β zβY (z)− y(0)zβ−1 zβ + ν = F (z) +K(z,t)Y (z), Y (z) { N (β) 1− β zβ zβ + ν −K(z,t) } = F (z) + N (β) 1− β y(0)zβ−1 zβ + ν , Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 8 of 22 Y (z) = F (z) + N (β) 1−β y(0)zβ−1 zβ+ν N (β) 1−β zβ zβ+ν −K(z,t) , (3) by applying the inverse Laplace transform, y(τ) = L−1 F (z) + N (β) 1−β y(0)zβ−1 zβ+ν N (β) 1−β zβ zβ+ν −K(z,t)  . (4) For many functions, it is not possible to evaluate the inverse Laplace transform in (4) analytically. In such cases, we employ numerical inversion methods. The inverse Laplace transform in (4) is given by y(τ) = 1 2πι ∫ θ+i∞ θ−i∞ ezτY (z)dz = 1 2πι ∫ Γ ezτY (z)dz. (5) Since directly finding the inverse Laplace transform in Equation (4) is often intractable, particularly due to fractional-order terms and the complex kernel, we employ numerical methods for this step. To approximate the inverse transform accurately, we utilize two numerical techniques: the Gauss-Hermite quadrature method and the modified Talbot’s method. These methods enable the efficient evaluation of the Bromwich integral given in Equation (5), providing both stability and high accuracy. In the subsequent sections, we detail the application of each method and justify their suitability for this class of problems. 3.1. Gauss-Hermite-quadrature method The Gauss-Hermite quadrature (GM) method provides an efficient numerical approx- imation for integrals of the form [29]:∫ ∞ −∞ e−t2h(t)dt ≈ ng∑ ξ=1 ηξh(tξ), (6) where tξ are the roots of the Hermite polynomial Hng(t) of degree ng, and ηξ are the corresponding weights. These weights are uniquely determined so that the approximation is exact when h(t) is a polynomial of degree up to 2ng − 1. The GM rule exhibits rapid convergence for smooth integrands, as characterized by the following error term: Theorem 3. The error term Rng(h) for the GM approximation in (6) is given by: Rng(h) = 1 2πi ∫ Γ ϕng(µ)h(µ)dµ, (7) where ϕng(µ) = Qng(µ) Hng(µ) , Qng(µ) = ∫ ∞ −∞ e−t2 Hng(t) µ− t dt. (8) Here, Γ is a contour in the complex plane enclosing the zeros of Hng(µ), and h(µ) is assumed to be analytic within Γ. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 9 of 22 The function Qng(µ), known as the Hermite function of the second kind, decays rapidly as |µ| → ∞, whereas Hng(µ) grows rapidly. Consequently, ϕng(µ) becomes negligible when Γ is chosen such that it is constrained by the region of analyticity of h(µ). To enhance the accuracy of the quadrature for functions that may not be optimally scaled, a linear transformation of the variable can be applied. Let t = ςτ , where ς ∈ R is a scaling parameter. This yields:∫ ∞ −∞ e−t2h(t)dt = ς ∫ ∞ −∞ e−ς2τ2h(ςτ)dτ = ς ∫ ∞ −∞ e−τ2 [ eτ 2(1−ς2)h(ςτ) ] dτ. (9) By defining a new function y(τ) = eτ 2(1−ς2)h(ςτ), the singularities of its analytic con- tinuation ŷ(µ) are shifted further from the real axis, which can significantly improve the convergence rate of the quadrature rule. For the numerical evaluation of the error contour integral (7), the contour Γ is param- eterized as: Γ : µ = ϑ(1 + iζ)2, −∞ < ζ < ∞, ϑ > 0. (10) This contour intersects the real axis at µ = ϑ and the imaginary axis at µ = ±2ϑi. Applying this parameterization and a subsequent scaling ζ = ςσ, the integral along the contour is transformed into a form amenable to the GM rule: y(τ) = ς 2πi ∫ ∞ −∞ e−σ2 eσ 2+µ(ςσ)τ ŷ(µ(ςσ))µ′(ςσ)dσ = ∫ ∞ −∞ e−σ2 h̃(σ)dσ, where the new integrand is defined as: h̃(σ) = ς 2πi eσ 2+µ(ςσ)τ ŷ(µ(ςσ))µ′(ςσ). The GM rule is then applied to approximate this integral. Due to the symmetry of the Gauss-Hermite nodes and weights for the real line, the approximation can be efficiently computed as: yApp(τ) ≈ 2ℜ  p∑ ξ=1 ηξh̃(σξ)  , where σξ are the positive roots of the Hermite polynomial Hng(t), and p = ng/2 (for even ng) or p = (ng + 1)/2 (for odd ng). This method leverages the rapid decay of the error term to provide a highly accurate numerical evaluation of the original integral. 3.1.1. Convergence of Gauss-Hermite quadrature method Theorem 4. Let y(τ) be bounded and analytic for τ ∈ (−∞,∞), and let exp(−τ2)y(τ) extend to a bounded analytic function in an infinite strip −c < Im(t) < c. For a fixed L > 0 and for each n ≥ 1, let In be the estimate of the integral I defined in (6) by applying the quadrature rule on the truncated interval [−Ln 1/3 g , Ln 1/3 g ]. Then, for some C > 0,∣∣I − Ing ∣∣ = O ( exp ( −Cn2/3 g )) , ng → ∞. (44) Proof. For proof, see [30]. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 10 of 22 3.2. Improved Talbot’s method The improved Talbot’s method (ITM) addresses the challenge of inverting Laplace transforms that involve highly oscillatory or slowly decaying terms by deforming the Bromwich contour. The key idea is to replace the standard vertical line with a con- tour that begins and ends in the left half-plane, ensuring it encircles all singularities of Y (z). This deformation forces the exponential term ezτ to decay rapidly along the chosen path, enabling an efficient numerical approximation of the integral via the midpoint or trapezoidal rule. In this work, we consider a Hankel contour parameterized as: Γ : z = z(ϑ), ϑ ∈ [π,−π], (11) where Rez(±π) = −∞. Substituting this into (5) yields y(τ) = 1 2πi ∫ Γ ezτY (z)dz = 1 2πι ∫ π −π ez(ϑ)τY (z(ϑ))z′(ϑ)dϑ. (12) Here, the contour z(ϑ) is defined as: z(ϑ) = nT t ζ(ϑ), ζ(ϑ) = −σ + θϑ cot(µϑ) + γιϑ, (13) where σ, θ, ϑ, and γ are user-defined parameters. The integral in (13) is approximated using the midpoint rule with nT equally spaced panels of step size h = 2π nT : yap(τ) ≈ 1 nT i nT∑ k=1 ez(ϑk)tY (z(ϑk))z ′(ϑk), ϑk = −π + ( k − 1 2 ) h. (14) 3.2.1. Convergence of Talbot’s method The ITM achieves exponential convergence when the integration contour is optimally parameterized. The convergence rate depends on the choice of parameters σ, ϑ, γ, and µ in the Hankel contour z(ϑ). While Talbot’s original contour (µ = 1) suggested in [31] yields satisfactory results, the authors of [32] modified µ to the range 0 < µ < 1, which significantly accelerates the convergence. The step size h = 2π WT also affects accuracy; a larger WT results in a smaller h, reducing the error exponentially, as evidenced by the error estimate O(e−1.358WT ). For optimal performance, the contour must balance being too close to the singularities and too far away. The parameter values found through experimentation in [32] are σ = 0.6122, θ = 0.5017, γ = 0.2645, and µ = 0.6407, which maximize the convergence efficiency and make the approach robust for inverting the Laplace transform. The error estimate of the method is given by: Errest = |yApp(τ)− y(τ)| = O(e−1.358nT ) . Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 11 of 22 4. Numerical examples This section demonstrates the performance of the proposed method through three rep- resentative examples. Detailed numerical experiments have been carried out to establish a strong foundation for comparative evaluation. To validate the accuracy of the inversion methods, we employ two error metrics: the maximum absolute error, denoted by L∞, and the root mean square error, denoted by Lrms, defined as follows: Lin = max 1≤j≤ng |y(τj)− yn(τj))|, and Lrms = √∑ng j=1(y(τj)− yn(τj))2 n , where ng denotes the number of nodes used in the GHQ method, y(τ) represents the exact solution, and yn(τ) represents the numerical solution of the problem considered. In all our experiments, the source term f(τ) and the initial conditions are derived from the exact solution. Problem 1 We consider the FOVIDE (1) with the exact solution y(τ) = τ2. Table 1: Comparison of the Improved GM and ITM methods for Problem 1, showing Lin and Lrms. GM ITM ng Lin Lrms CPU (s) nT Lin Lrms CPU(s) 15 2.55× 10−08 6.59× 10−09 0.081 12 3.22× 10−04 9.31× 10−05 0.012 16 5.61× 10−09 5.61× 10−09 0.109 14 2.91× 10−05 7.77× 10−06 0.015 17 1.30× 10−09 3.15× 10−10 0.078 16 2.91× 10−06 6.29× 10−07 0.012 18 2.83× 10−10 6.68× 10−11 0.113 18 2.11× 10−07 4.98× 10−08 0.012 19 6.48× 10−11 1.49× 10−11 0.083 20 1.73× 10−08 3.87× 10−09 0.013 20 1.40× 10−11 3.13× 10−12 0.113 22 1.39× 10−09 2.97× 10−10 0.012 21 3.09× 10−12 6.75× 10−13 0.091 24 1.10× 10−10 2.25× 10−11 0.013 22 6.46× 10−13 1.38× 10−13 0.108 26 8.60× 10−12 1.69× 10−12 0.012 23 1.40× 10−13 2.93× 10−14 0.105 28 6.66× 10−13 1.26× 10−13 0.016 24 2.66× 10−13 5.44× 10−15 0.102 30 4.81× 10−14 8.78× 10−15 0.011 32 5.78× 10−15 1.02× 10−15 0.012 Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 12 of 22 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) 0 0.2 0.4 0.6 0.8 1 10-17 10-16 10-15 10-14 10-13 (b) Figure 1: (a) Numerical and exact solutions for Problem 1, exhibiting close agreement. (b) Error comparison for Problem 1, showing the ITM outperforming the improved GM. 10 12 14 16 18 20 22 24 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 (a) 10 15 20 25 30 35 40 45 50 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 (b) Figure 2: Convergence analysis of error norms versus quadrature nodes for Problem 1: (a) GM; (b) ITM. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 13 of 22 (a) (b) Figure 3: Comparison of error norms for Problem 1 computed via the GM across the βτ plane: (a) Lin error norm; (b) Lrms error norm. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2 4 6 8 10 12 14 16 18 10-15 (a) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.5 1 1.5 2 2.5 3 10-15 (b) Figure 4: Contour plots of numerical errors for Problem 1 using the ITM in the βτ plane: (a) Lin error norm; (b) Lrms error norm . Problem 2 consider the FOVIDE 1 with the exact solution y(τ) = eτ . Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 14 of 22 Table 2: Comparison of the Improved GM and ITM methods for Problem 2, showing Lin and Lrms. GM ITM ng Lin Lrms CPU (s) nT Lin Lrms CPU(s) 15 4.5095×10−06 1.1643×10−06 0.106230 20 1.5882×10−05 3.5512×10−06 0.011312 16 1.3332×10−06 3.3330×10−07 0.163458 22 2.2628×10−06 4.8242×10−07 0.011950 17 4.0362×10−07 9.7893×10−08 0.114148 24 3.1176×10−07 6.3638×10−08 0.013142 18 1.1701×10−07 2.7580×10−08 0.084061 26 4.1707×10−08 8.1794×10−09 0.022175 19 3.4658×10−08 7.9510×10−09 0.126576 28 5.4354×10−09 1.0272×10−09 0.015180 20 9.8581×10−09 2.2043×10−09 0.087279 30 6.9187×10−10 1.2632×10−10 0.013981 21 2.8228×10−09 6.1598×10−10 0.168586 32 8.6329×10−11 1.5261×10−11 0.013615 22 7.7622×10−10 1.6549×10−10 0.093258 34 1.0887×10−11 1.8671×10−12 0.015629 23 2.1792×10−10 4.5439×10−11 0.125601 36 9.1706×10−13 1.5284×10−13 0.021020 24 5.8898×10−11 1.2022×10−11 0.124419 38 3.5357×10−13 5.7357×10−14 0.016675 40 8.1681×10−14 1.2915×10−14 0.013102 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 Exact Approx (a) 0 0.2 0.4 0.6 0.8 1 10-17 10-16 10-15 10-14 10-13 10-12 10-11 10-10 (b) Figure 5: (a) Numerical and exact solutions for Problem 2. The visual indistinguishability of the curves demonstrates the high accuracy of the numerical method. (b) Comparison of Lin and Lrms errors for the ITM and GM, demonstrating the superior performance of the ITM. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 15 of 22 10 12 14 16 18 20 22 24 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 (a) 10 15 20 25 30 35 40 45 50 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 (b) Figure 6: Convergence analysis of error norms versus quadrature nodes for Problem 2: (a) GM; (b) ITM. (a) (b) Figure 7: Surface plots of error norms for Problem 2 using the GM method in the βτ plane: (a) Lin error norm; (b) Lrms error norm. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 16 of 22 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 10-13 (a) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 10-14 (b) Figure 8: Contour plots of error norms for Problem 2 using the ITM in the βτ plane: (a) Lin error norm; (b) Lrms error norm. Problem 3 Consider the FOVIDE (1) with the exact solution y(τ) = sin(τ). Table 3: Comparison of the Improved GM and ITM methods for Problem 3, showing Lin and Lrms. GM ITM ng Lin Lrms CPU (s) nT Lin Lrms CPU(s) 15 2.0879×10−08 5.3910×10−09 0.087910 16 5.3996×10−06 1.3499×10−06 0.016834 16 8.7414×10−09 2.1853×10−09 0.098304 18 8.8526×10−07 2.0866×10−07 0.011818 17 2.9932×10−09 7.2597×10−10 0.097600 20 9.6216×10−08 2.1514×10−08 0.013080 18 8.8385×10−10 2.0833×10−10 0.129561 22 6.8963×10−09 1.4703×10−09 0.014854 19 2.4670×10−10 5.6596×10−11 0.096193 24 1.1958×10−10 2.4408×10−11 0.015823 20 6.1596×10−11 1.3773×10−11 0.093845 26 5.9034×10−11 1.1577×10−11 0.010973 21 1.4559×10−11 3.1770×10−12 0.089256 28 1.2398×10−11 2.3431×10−12 0.017625 22 3.0064×10−12 6.4096×10−13 0.088783 30 1.6377×10−12 2.9900×10−13 0.011846 23 5.7976×10−13 1.2089×10−13 0.095613 32 1.8074×10−13 3.1951×10−14 0.013646 24 9.5035×10−14 1.9399×10−14 0.103664 34 5.7458×10−14 9.8539×10−15 0.018100 36 3.6749×10−14 6.1248×10−15 0.011105 Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 17 of 22 0 0.2 0.4 0.6 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 (a) 0 0.2 0.4 0.6 0.8 1 10-17 10-16 10-15 10-14 10-13 (b) Figure 9: (a) Numerical and exact solutions for Problem 3, showing a near-perfect overlap. (b) Comparison of error norms for Problem 3, demonstrating the outperformance of the GM over the ITM. 10 12 14 16 18 20 22 24 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 (a) 10 15 20 25 30 35 40 45 50 10 -16 10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 (b) Figure 10: Convergence analysis of error norms versus quadrature nodes for Problem 3: (a) GM; (b) ITM. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 18 of 22 (a) (b) Figure 11: Surface plots of error norms for Problem 3 using the GM method in the βτ plane: (a) Lin error norm; (b) Lrms error norm. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.5 2 2.5 3 3.5 4 4.5 5 10-14 (a) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 10-15 (b) Figure 12: Contour plots of error norms for Problem 3 using the ITM in the βτ plane: (a) Lin error norm; (b) Lrms error norm.. 5. Discussion of Numerical Results The numerical performance of the GM and the ITM was rigorously evaluated across three distinct test problems. A comprehensive error analysis was conducted using both the infinite norm (Ln) and the root mean square norm (Lrms). The results, detailed in Tables 1–3, consistently demonstrate the superior accuracy and convergence behavior of the ITM over the GM method, with the ITM achieving significantly lower error values in all cases. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 19 of 22 This initial comparison was performed using a constant kernel function, K(τ, t) = 1, for all problems. The excellent agreement between the exact and numerical solutions, visually confirmed in Figures 1a, 5a, and 9a, alongside the direct error norm comparisons in Figures 1b, 5b, and 9b, underscores the high baseline accuracy of the proposed numerical inversion framework and clearly establishes the ITM’s advantage. To further investigate the robustness and versatility of the methods, the analysis was extended to problems with variable kernels: K(τ, t) = sin(τ − t) for Problem 1, K(τ, t) = τ − t for Problem 2, and K(τ, t) = e(τ−t) for Problem 3. Figures 2a, 6a, 10a and 2b, 6b, 10b illustrate the convergence of the L∞ and Lrms error norms with respect to the number of quadrature nodes for the GM and ITM, respectively. A consistent decrease in error is observed as the number of quadrature nodes increases, confirming the expected convergence for both methods. The behavior of the error landscapes was thoroughly examined through surface and contour plots of the Lin and Lrms norms across the ατ plane. The surface plots for the GM (Figures 3a, 3b, 7a, 7b, 11a, 11b) confirm its capability to maintain high accuracy even with more complex kernel functions. However, the corre- sponding contour plots for the ITM (Figures 4a, 4b, 8a, 8b, 12a, 12b) reveal a consistently more refined and superior error distribution, characterized by tighter, more concentrated contours indicating lower and more stable error levels across the parameter space. The nu- merical scheme demonstrates remarkable stability throughout the entire α ∈ (0, 1) range, with only mild degradation as α → 0. This robustness makes it suitable for practical applications where fractional orders may vary, confirming that both inversion methods effectively handle the mathematical complexities introduced by the MABC derivative op- erator Collectively, the results from all three problems, with both constant and variable ker- nels, establish that the ITM is a more robust, efficient, and highly accurate computational technique for solving this class of equations, demonstrating its reliability for a wider range of applications. 6. Conclusion This article developed and implemented an efficient numerical technique for solving a class of FOVIDEs including the MABC derivative. The proposed methodology used the Laplace transform to convert the non-local fractional order complex problem into a simple algebraic equation. The important step of inverting this solution back to the time domain was performed using two highly accurate numerical techniques: the GM and the ITM. To validate and compare the performance of two inversion methods numerical experiments were conducted. The results clearly showed the convergence and high accuracy of both methods as indicated by Lin and Lrms error norms for each test problem. However, a comprehensive analysis demonstrated that the accuracy of ITM was better compared to GM. Overall, the proposed LT-based numerical method proves out to be a very accurate and flexible approach, particularly when the ITM and GM are used for inversion. Further, in addition to offering a reliable numerical solution for complex FOVIDEs, this suggested Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 20 of 22 approach establishes an efficient foundation for further study. The efficacy of this method allows it to be a suitable choice for systems in physical, engineering, and biological sciences. In future work, we plan to extend this methodology to solve more complex problems, such as FOVIDEs with delay terms and FOVIDEs with variable coefficients. Acknowledgements The authors F. Hasan and N. Mlaiki would like to thank Prince Sultan University for paying the APC and for the support through the TAS research lab. Competing interests There are no conflicting interests, according to the authors. Author’s contributions Each author contributed equally to the writing of this work, and they have all read and approved the finished work. Declarations Ethical Approval Not applicable. Funding This work did not receive any external funding. References [1] Richard T Baillie. Long memory processes and fractional integration in econometrics. Journal of econometrics, 73(1):5–59, 1996. [2] Vasily E Tarasov. Fractional integro-differential equations for electromagnetic waves in dielectric media. Theoretical and Mathematical Physics, 158(3):355–359, 2009. [3] Gary W Bohannan. Analog fractional order controller in temperature and motor control applications. Journal of Vibration and Control, 14(9-10):1487–1498, 2008. [4] Ronald L Bagley and Peter J Torvik. A theoretical basis for the application of fractional calculus to viscoelasticity. Journal of rheology, 27(3):201–210, 1983. [5] Lokenath Debnath. Recent applications of fractional calculus to science and engineering. International Journal of Mathematics and Mathematical Sciences, 2003(54):3413–3442, 2003. [6] Kamal Shah, Manar A Alqudah, Fahd Jarad, and Thabet Abdeljawad. Semi- analytical study of pine wilt disease model with convex rate under caputo–febrizio fractional order derivative. Chaos, Solitons & Fractals, 135:109754, 2020. Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 21 of 22 [7] Zubair Ahmad, Giuliano Bonanomi, Angelamaria Cardone, Annalisa Iuorio, Gerardo Toraldo, and Francesco Giannino. Fractal-fractional sirs model for the disease dy- namics in both prey and predator with singular and nonsingular kernels. Journal of Biological Systems, 32(04):1487–1520, 2024. [8] Zubair Ahmad, Sherif A El-Kafrawy, Thamir A Alandijany, Francesco Giannino, Ahmed A Mirza, Mai M El-Daly, Arwa A Faizo, Leena H Bajrai, Mohammad Amjad Kamal, and Esam I Azhar. A global report on the dynamics of covid-19 with quar- antine and hospitalization: A fractional order model with non-local kernel. Compu- tational Biology and Chemistry, 98:107645, 2022. [9] Sertan Alkan and Veysel Fuat Hatipoglu. Approximate solutions of volterra-fredholm integro-differential equations of fractional order. 2017. [10] Abbas Saadatmandi and Mehdi Dehghan. A legendre collocation method for frac- tional integro-differential equations. Journal of Vibration and Control, 17(13):2050– 2058, 2011. [11] Dilek Varol Bayram and Ayşegül Daşcıoğlu. A method for fractional volterra integro- differential equations by laguerre polynomials. Advances in Difference Equations, 2018(1):466, 2018. [12] Dimple Rani and Vinod Mishra. Modification of laplace adomian decomposition method for solving nonlinear volterra integral and integro-differential equations based on newton raphson formula. European Journal of Pure and Applied Mathematics, 11(1):202–214, 2018. [13] Fazal Haq, Kamal Shah, Ghaus ur Rahman, and Muhammad Shahzad. Numerical so- lution of fractional order smoking model via laplace adomian decomposition method. Alexandria Engineering Journal, 57(2):1061–1069, 2018. [14] Yasir Nawaz. Variational iteration method and homotopy perturbation method for fourth-order fractional integro-differential equations. Computers & Mathematics with Applications, 61(8):2330–2341, 2011. [15] Kamal Shah, Hammad Khalil, and Rahmat Ali Khan. Analytical solutions of frac- tional order diffusion equations by natural transform method. Iranian Journal of Science and Technology, Transactions A: Science, 42(3):1479–1490, 2018. [16] Emran Tohidi, MM Ezadkhah, and S Shateyi. Numerical solution of nonlinear frac- tional volterra integro-differential equations via bernoulli polynomials. In Abstract and Applied Analysis, volume 2014, page 162896. Wiley Online Library, 2014. [17] EA Rawashdeh. Legendre wavelets method for fractional integro-differential equa- tions. Applied Mathematical Sciences, 5(2):2467–2474, 2011. [18] Ajay Kumar, Vikas Gupta, and Saurabh Kumar. A numerical approach to solve the system of fractional fredholm–volterra integro-differential equations using rationalized haar wavelets. Mathematical Methods in the Applied Sciences, 48:11067–11088, 2025. [19] Omar Abu Arqub and Banan Maayah. Fitted fractional reproducing kernel algorithm for the numerical solutions of abc–fractional volterra integro-differential equations. Chaos, Solitons & Fractals, 126:394–402, 2019. [20] Xuefei Dai, Chaoyue Guan, and Jing Niu. An efficient numerical algorithm for solving nonlinear fractional volterra integro-differential equation. Computational and Applied Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6970 22 of 22 Mathematics, 43(1):66, 2024. [21] Xiaohua Ma and Chengming Huang. Numerical solution of fractional integro- differential equations by a hybrid collocation method. Applied Mathematics and Computation, 219(12):6750–6760, 2013. [22] Amr Mahdy, Abbas S Nagdy, and Doaa Sh Mohamed. Solution of fractional integro- differential equations using least squares and shifted legendre methods. Journal of Applied Mathematics and Computational Mechanics, 23(1):59–70, 2024. [23] N Rajagopal, S Balaji, R Seethalakshmi, and VS Balaji. A new numerical method for fractional order volterra integro-differential equations. Ain Shams Engineering Journal, 11(1):171–177, 2020. [24] Pongsakorn Sunthrayuth, Roman Ullah, Adnan Khan, Rasool Shah, Jeevan Kafle, Ibrahim Mahariq, and Fahd Jarad. Numerical analysis of the fractional-order non- linear system of volterra integro-differential equations. Journal of Function Spaces, 2021(1):1537958, 2021. [25] Michele Caputo and Mauro Fabrizio. A new definition of fractional derivative without singular kernel. Progress in fractional differentiation & applications, 1(2):73–85, 2015. [26] Abdon Atangana and Dumitru Baleanu. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Thermal Science, 20:763–769, 2016. [27] Mohammed Al-Refai and Dumitru Baleanu. On an extension of the operator with mittag-leffler kernel. Fractals, 30(05):2240129, 2022. [28] Wen-Hua Huang, Muhammad Samraiz, Ahsan Mehmood, Dumitru Baleanu, Gauhar Rahman, and Saima Naheed. Modified atangana-baleanu fractional operators involv- ing generalized mittag-leffler function. Alexandria Engineering Journal, 75:639–648, 2023. [29] JAC Weideman. Gauss–hermite quadrature for the bromwich integral. SIAM Journal on Numerical Analysis, 57(5):2200–2216, 2019. [30] Lloyd N Trefethen. Exactness of quadrature formulas. Siam Review, 64(1):132–150, 2022. [31] Alan Talbot. The accurate numerical inversion of laplace transforms. IMA Journal of Applied Mathematics, 23(1):97–120, 1979. [32] Benedict Dingfelder and JAC Weideman. An improved talbot method for numerical laplace transform inversion. Numerical Algorithms, 68(1):167–183, 2015.